Subversion Repositories f9daq

Rev

Rev 167 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. #include "../include/sipmscan.h"
  2. #include "../include/workstation.h"
  3.  
  4. #include <stdio.h>
  5. #include <stdlib.h>
  6.  
  7. // Peak detection function
  8. int npeaks = 20;
  9. double FindPeaks(double *x, double *par)
  10. {
  11.    double result = 0;
  12.    for(int i = 0; i < npeaks; i++)
  13.    {
  14.       double norm = par[3*i];
  15.       double mean = par[3*i+1];
  16.       double sigma = par[3*i+2];
  17.       result += norm*TMath::Gaus(x[0], mean, sigma);
  18.    }
  19.    return result;
  20. }
  21.  
  22. // File browser for selecting the dark run histogram
  23. void TGAppMainFrame::SelectDarkHist()
  24. {
  25.    TGFileInfo file_info;
  26.    const char *filetypes[] = {"Histograms",histextall,0,0};
  27.    char *cTemp;
  28.    file_info.fFileTypes = filetypes;
  29.    cTemp = new char[1024];
  30.    sprintf(cTemp, "%s/results", rootdir);
  31.    file_info.fIniDir = StrDup(cTemp);
  32.    file_info.fMultipleSelection = kFALSE;
  33.    new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
  34.    delete[] cTemp;
  35.  
  36.    if(file_info.fFilename != NULL)
  37.    {
  38.       darkRun->widgetTE->SetText(file_info.fFilename);
  39.       fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
  40.       fileList->Layout();
  41.    }
  42.    else
  43.       darkRun->widgetTE->SetText("");
  44. }
  45.  
  46. // Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
  47. void TGAppMainFrame::AnalysisDefaults()
  48. {
  49.    printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
  50.    if(analTab->GetCurrent() == 0)       // Integrate spectrum
  51.    {
  52.       intSpect->widgetChBox[0]->SetState(kButtonUp);
  53.       intSpect->widgetChBox[1]->SetState(kButtonUp);
  54.       intSpect->widgetChBox[2]->SetState(kButtonDown);
  55.       resol2d->widgetNE[0]->SetNumber(40);
  56.       resol2d->widgetNE[1]->SetNumber(40);
  57.    }
  58.    else if(analTab->GetCurrent() == 1)  // Relative PDE
  59.    {
  60.       relPde->widgetChBox[1]->SetState(kButtonDown);
  61.       midPeak->widgetChBox[0]->SetState(kButtonUp);
  62.       zeroAngle->widgetNE[0]->SetNumber(0.00);
  63.    }
  64.    else if(analTab->GetCurrent() == 2)  // Breakdown voltage
  65.    {
  66.       minPeak->widgetNE[0]->SetNumber(2);
  67.       minPeak->widgetNE[0]->SetNumber(1);
  68.    }
  69.    else if(analTab->GetCurrent() == 3)  // Surface scan
  70.    {
  71.       surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
  72.       surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
  73.       resolSurf->widgetNE[0]->SetNumber(40);
  74.       resolSurf->widgetNE[1]->SetNumber(40);
  75.    }
  76. }
  77.  
  78. // Analysis handle function
  79. void TGAppMainFrame::AnalysisHandle(int type)
  80. {
  81.    TList *files;
  82.    bool createTab = true;
  83.    int tabid = -1;
  84.    int analtab = analTab->GetCurrent();
  85.  
  86.    int analtype;
  87.    if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
  88.       analtype = 0;
  89.    else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) )
  90.       analtype = 1;
  91.    else if( analtab == 1 )
  92.       analtype = 2;
  93.    else if( analtab == 2 )
  94.       analtype = 3;
  95.    else if( analtab == 3 )
  96.       analtype = 4;
  97.  
  98.    // Only integrate spectrum or make relative PDE
  99.    if(type == 0)
  100.    {
  101.       files = new TList();
  102.       fileList->GetSelectedEntries(files);
  103.  
  104.       if( analtype == 0 )
  105.          IntegSpectrum(files, 0, 0);
  106.  
  107.       if( intSpect->widgetChBox[0]->IsDown() )
  108.          IntegSpectrum(files, 1, 0);
  109.  
  110.       if( intSpect->widgetChBox[1]->IsDown() )
  111.          IntegSpectrum(files, 2, 0);
  112.  
  113.       if( analtype == 2 )
  114.          PhotonMu(files, 0);
  115.  
  116.       if( analtype == 3 )
  117.          BreakdownVolt(files, 0);
  118.  
  119.       if( analtype == 4 )
  120.          SurfaceScan(files, 0);
  121.    }
  122.    // Integrate spectrum or make relative PDE and open edit window
  123.    else if(type == 1)
  124.    {
  125.       files = new TList();
  126.       fileList->GetSelectedEntries(files);
  127.  
  128.       // Prepare a new analysis edit tab, if we want to edit plots
  129.       for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
  130.       {
  131.          if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
  132.          {
  133.             createTab = false;
  134.             tabid = i;
  135.          }
  136.          
  137.          if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
  138.       }
  139.    
  140.       if(files->GetSize() > 0)
  141.       {
  142.          TempAnalysisTab(fTab, createTab, &tabid, analtype);
  143.  
  144.          // Integrate spectra
  145.          if( analtype == 0 )
  146.             IntegSpectrum(files, 0, 1);
  147.  
  148.          if( intSpect->widgetChBox[0]->IsDown() )
  149.             IntegSpectrum(files, 1, 1);
  150.  
  151.          if( intSpect->widgetChBox[1]->IsDown() )
  152.             IntegSpectrum(files, 2, 1);
  153.  
  154.          if( analtype == 2 )
  155.             PhotonMu(files, 1);
  156.  
  157.          if( analtype == 3 )
  158.             BreakdownVolt(files, 1);
  159.  
  160.          if( analtype == 4 )
  161.             SurfaceScan(files, 1);
  162.  
  163.          fTab->SetTab(tabid);
  164.       }
  165.    }
  166.  
  167.    delete files;
  168. }
  169.  
  170. // Analysis functions ----------------------------------
  171. void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
  172. {
  173.    unsigned int nrfiles = files->GetSize();
  174.    char ctemp[1024];
  175.    char exportname[1024];
  176.    int j, k = 0, m = 0;
  177.  
  178.    TTree *header_data, *meas_data;
  179.    double *integralCount, *integralAcc;
  180.    integralCount = new double[nrfiles];
  181.    integralAcc = new double[nrfiles];
  182.    double *surfxy, *surfz;
  183.    surfxy = new double[nrfiles];
  184.    surfz = new double[nrfiles];
  185.    double minInteg, maxInteg;
  186.    bool norm = intSpect->widgetChBox[2]->IsDown();
  187.    double curzval;
  188.    bool edge2d = false;
  189.    TCanvas *gCanvas;
  190.  
  191.    float progVal = 0;
  192.    analysisProgress->widgetPB->SetPosition(progVal);
  193.    gVirtualX->Update(1);
  194.  
  195.    // Zero the integral count and accumulated vaules
  196.    for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
  197.  
  198.    // Start if we select at least one file
  199.    if(nrfiles > 0)
  200.    {
  201.       for(int i = 0; i < (int)nrfiles; i++)
  202.       {
  203.          if(files->At(i))
  204.          {
  205.             // Read the X,Y and Z positions from header and ADC and TDC values from the measurements
  206.             sprintf(ctemp, "%s", files->At(i)->GetTitle());
  207.             inroot = new TFile(ctemp, "READ");
  208.          
  209.             inroot->GetObject("header_data", header_data);
  210.             inroot->GetObject("meas_data", meas_data);
  211.          
  212.             header_data->SetBranchAddress("xpos", &evtheader.xpos);
  213.             header_data->GetEntry(0);
  214.             header_data->SetBranchAddress("ypos", &evtheader.ypos);
  215.             header_data->GetEntry(0);
  216.             header_data->SetBranchAddress("zpos", &evtheader.zpos);
  217.             header_data->GetEntry(0);
  218.  
  219.             char rdc[256];
  220.             j = selectCh->widgetNE[0]->GetNumber();
  221.             double rangetdc[2];
  222.             rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
  223.             rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
  224.    
  225.             k = 0;
  226.             m = 0;
  227.          
  228.             // Reading the data
  229.             for(int e = 0; e < meas_data->GetEntries(); e++)
  230.             {
  231.                sprintf(rdc, "ADC%d", j);
  232.                meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
  233.                meas_data->GetEntry(e);
  234.          
  235.                sprintf(rdc, "TDC%d", j);
  236.                meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
  237.                meas_data->GetEntry(e);
  238.    
  239.                // Use data point only if it is inside the TDC window
  240.                if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
  241.                {
  242.                   k++;
  243.                   m += evtdata.adcdata[j];
  244.                }
  245.             }
  246.  
  247.             // X, Y and Z values from each file (table units or microns)
  248.             if(posUnits->widgetCB->GetSelected() == 0)
  249.             {
  250.                if(direction == 1)
  251.                   surfxy[i] = (double)(evtheader.xpos);
  252.                else if(direction == 2)
  253.                   surfxy[i] = (double)(evtheader.ypos);
  254.                surfz[i] = (double)(evtheader.zpos);
  255.             }
  256.             else if(posUnits->widgetCB->GetSelected() == 1)
  257.             {
  258.                if(direction == 1)
  259.                   surfxy[i] = (double)(evtheader.xpos*lenconversion);
  260.                else if(direction == 2)
  261.                   surfxy[i] = (double)(evtheader.ypos*lenconversion);
  262.                surfz[i] = (double)(evtheader.zpos*lenconversion);
  263.             }
  264.  
  265.             // Check if we have different Z values: if no, just make the edge plots; if yes, save edge plots and make a 2d edge plot
  266.             if(i == 0) curzval = surfz[i];
  267.             else
  268.             {
  269.                if(surfz[i] != curzval)
  270.                   edge2d = true;
  271.             }
  272.  
  273.             // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
  274.             if(direction == 0)
  275.             {
  276.                if(norm)
  277.                {
  278.                   integralCount[i] += ((double)m)/((double)k);
  279.                   printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
  280.                }
  281.                else
  282.                {
  283.                   integralCount[i] += m;
  284.                   printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
  285.                }
  286.             }
  287.             else
  288.             {
  289.                if(norm)
  290.                   integralCount[i] += ((double)m)/((double)k);
  291.                else
  292.                   integralCount[i] += m;
  293.             }
  294.            
  295.             inroot->Close();
  296.             delete inroot;
  297.          }
  298.  
  299.          // Update the progress bar
  300.          progVal = (float)(75.00/nrfiles)*i;
  301.          analysisProgress->widgetPB->SetPosition(progVal);
  302.          gVirtualX->Update(1);
  303.       }
  304.  
  305.       printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
  306.  
  307.       // If only an integral is needed, do not plot and exit here
  308.       if( direction == 0 )
  309.       {
  310.          delete[] integralCount;
  311.          delete[] surfxy;
  312.          delete[] surfz;
  313.          return;
  314.       }
  315.  
  316.       // Current z value and the accumulated counter
  317.       curzval = surfz[0];
  318.       j = 0;
  319.       int acc = 0;
  320.       int zb;
  321.       for(int i = 0; i <= (int)nrfiles; i++)
  322.       {
  323.          // Collect the accumulated integral in order to produce a PDF from a CDF
  324.          // While we are at the same Z value, save under one set
  325.          if( (surfz[i] == curzval) && (acc != nrfiles) )
  326.          {
  327.             integralAcc[j] = integralCount[i];
  328.             if(DBGSIG) printf("IntegSpectrum(): Integral check 1 (i=%d,j=%d,z=%.2lf): %lf\t%lf\n", i, j, surfz[i], integralCount[i], integralAcc[j]);
  329.             j++;
  330.             acc++;
  331.          }
  332.          // When we switch to a new set of Z values and at the end, we must save the previous ones to make 1D edge plots
  333.          else
  334.          {
  335.             // Find minimal and maximal integral values to subtract the offset and normate PDF to 1
  336.             NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
  337.  
  338.             if(acc != nrfiles)
  339.             {
  340.                curzval = surfz[i];
  341.                // PDF and CDF plot
  342.                PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
  343.                i--;
  344.                j = 0;
  345.             }
  346.             else
  347.             {
  348.                // PDF and CDF plot
  349.                PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
  350.                i--;
  351.                break;
  352.             }
  353.          }
  354.  
  355.          // Update the progress bar
  356.          progVal = (float)(15.00/nrfiles)*i+75.00;
  357.          analysisProgress->widgetPB->SetPosition(progVal);
  358.          gVirtualX->Update(1);
  359.       }
  360.  
  361.       // Make the 2D edge plot
  362.       if(edge2d)
  363.       {
  364.          if(edit == 0)
  365.             gCanvas = new TCanvas("canv","canv",900,900);
  366.          else
  367.             gCanvas = tempAnalysisCanvas->GetCanvas();
  368.  
  369.          double range[4];
  370.          TGraph2D *gScan2D;
  371.          gScan2D = new TGraph2D();
  372.          range[0] = TMath::MinElement(nrfiles, surfxy);
  373.          range[1] = TMath::MaxElement(nrfiles, surfxy);
  374.          range[2] = TMath::MinElement(nrfiles, surfz);
  375.          range[3] = TMath::MaxElement(nrfiles, surfz);
  376.  
  377.          for(int i = 0; i < nrfiles; i++)
  378.          {
  379.             if(DBGSIG)
  380.                printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
  381.             gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
  382.  
  383.             // Update the progress bar
  384.             progVal = (float)(9.00/nrfiles)*i+90.00;
  385.             analysisProgress->widgetPB->SetPosition(progVal);
  386.             gVirtualX->Update(1);
  387.          }
  388.  
  389.          gCanvas->cd();
  390.          gStyle->SetPalette(1);
  391.          gCanvas->SetLeftMargin(0.15);
  392.          gCanvas->SetRightMargin(0.126);
  393.          gCanvas->SetTopMargin(0.077);
  394.          gScan2D->Draw("COLZ");
  395.          
  396.          gCanvas->Modified();
  397.          gCanvas->Update();
  398.  
  399.          gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
  400.          gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
  401.          
  402.          gCanvas->Modified();
  403.          gCanvas->Update();
  404.          
  405.          if(direction == 1)
  406.             gScan2D->GetXaxis()->SetTitle("X [#mum]");
  407.          else if(direction == 2)
  408.             gScan2D->GetXaxis()->SetTitle("Y [#mum]");
  409.  
  410.  
  411.          gScan2D->GetXaxis()->SetTitleOffset(1.3);
  412.          gScan2D->GetXaxis()->CenterTitle(kTRUE);
  413.          gScan2D->GetXaxis()->SetLabelSize(0.027);
  414.          gScan2D->GetXaxis()->SetLabelOffset(0.02);
  415.          gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
  416.          gScan2D->GetXaxis()->SetNoExponent(kTRUE);
  417.          gScan2D->GetYaxis()->SetTitleOffset(1.9);
  418.          gScan2D->GetYaxis()->CenterTitle(kTRUE);
  419.          gScan2D->GetYaxis()->SetLabelSize(0.027);
  420.          gScan2D->GetXaxis()->SetLabelOffset(0.02);
  421.          gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
  422.          gScan2D->GetYaxis()->SetNoExponent(kTRUE);
  423.  
  424. /*         TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
  425.          yax->SetMaxDigits(4);*/
  426.  
  427.          if(!cleanPlots)
  428.          {
  429.             if(direction == 1)
  430.             {
  431.                if(posUnits->widgetCB->GetSelected() == 0)
  432.                   gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
  433.                else if(posUnits->widgetCB->GetSelected() == 1)
  434.                   gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
  435.             }
  436.             else if(direction == 2)
  437.             {
  438.                if(posUnits->widgetCB->GetSelected() == 0)
  439.                   gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
  440.                else if(posUnits->widgetCB->GetSelected() == 1)
  441.                   gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
  442.             }
  443.          }
  444.          else
  445.          {
  446.             if(direction == 1)
  447.             {
  448.                if(posUnits->widgetCB->GetSelected() == 0)
  449.                   gScan2D->SetTitle(";X [table units];Z [table units]");
  450.                else if(posUnits->widgetCB->GetSelected() == 1)
  451.                   gScan2D->SetTitle(";X [#mum];Z [#mum]");
  452.             }
  453.             else if(direction == 2)
  454.             {
  455.                if(posUnits->widgetCB->GetSelected() == 0)
  456.                   gScan2D->SetTitle(";Y [table units];Z [table units]");
  457.                else if(posUnits->widgetCB->GetSelected() == 1)
  458.                   gScan2D->SetTitle(";Y [#mum];Z [#mum]");
  459.             }
  460.          }
  461.  
  462.          gCanvas->Modified();
  463.          gCanvas->Update();
  464.  
  465.          remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
  466.          sprintf(exportname, "%s", ctemp);
  467.          remove_from_last(exportname, '_', ctemp);
  468.          if(direction == 1)
  469.             sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
  470.          else if(direction == 2)
  471.             sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
  472.          gCanvas->SaveAs(exportname);
  473.  
  474.          // Update the progress bar
  475.          analysisProgress->widgetPB->SetPosition(100.0);
  476.          gVirtualX->Update(1);
  477.  
  478.          if(edit == 0)
  479.          {
  480.             delete gScan2D;
  481.             delete gCanvas;
  482.          }
  483.       }
  484.       else
  485.       {
  486.          // Update the progress bar
  487.          analysisProgress->widgetPB->SetPosition(100.0);
  488.          gVirtualX->Update(1);
  489.       }
  490.    }
  491. }
  492.  
  493. void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit)
  494. {
  495.    TGraph *gScan[2];
  496.    int pdfmax = -1;
  497.    int count = 0;
  498.    char ctemp[1024];
  499.    char exportname[1024];
  500.    TCanvas *gCanvas;
  501.  
  502.    // Prepare the CDF plot
  503.    gScan[1] = new TGraph();
  504.    for(int i = 0; i < points; i++)
  505.    {
  506.       count = filenr - points + i;
  507.       gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max));
  508.       if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", i, (double)xy[count], (double)integAcc[i]/(*max));
  509.  
  510.       if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) )
  511.          pdfmax = (integAcc[i+1]-integAcc[i])/(*max);
  512.    }
  513.  
  514.    pdfmax = (TMath::Ceil(pdfmax*10))/10.;
  515.  
  516.    // Prepare the PDF plot
  517.    gScan[0] = new TGraph();
  518.    for(int i = points-1; i >= 0; i--)
  519.    {
  520.       count = (filenr-1) - (points-1) + i;
  521.       // Set any negative values of the PDF to 0
  522.       if( (integAcc[i]-integAcc[i-1])/(*max) < 0 )
  523.          gScan[0]->SetPoint(i, (double)xy[count], 0);
  524.       else
  525.          gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max));
  526.       if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max));
  527.    }
  528.  
  529.    remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp);
  530.    sprintf(exportname, "%s_edge.pdf", ctemp);
  531.  
  532.    if(edit == 0)
  533.       gCanvas = new TCanvas("canv1d","canv1d",1200,900);
  534.    else
  535.       gCanvas = tempAnalysisCanvas->GetCanvas();
  536.  
  537.    // Fit the PDF with a gaussian
  538.    gScan[0]->Fit("gaus","Q");
  539.    gScan[0]->GetFunction("gaus")->SetNpx(400);
  540.  
  541.    gStyle->SetOptFit(1);
  542.  
  543.    gScan[1]->Draw("AL");
  544.    gPad->Update();
  545.    gScan[0]->Draw("LP");
  546.  
  547.    gCanvas->Modified();
  548.    gCanvas->Update();
  549.  
  550.    TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats");
  551.    if(!cleanPlots)
  552.    {
  553.       stats->SetX1NDC(0.86); stats->SetX2NDC(1.0);
  554.       stats->SetY1NDC(0.87); stats->SetY2NDC(1.0);
  555.    }
  556.    else
  557.    {
  558.       stats->SetX1NDC(1.1); stats->SetX2NDC(1.3);
  559.       stats->SetY1NDC(1.1); stats->SetY2NDC(1.3);
  560.    }
  561.  
  562.    gCanvas->SetGridx(1);
  563.    gCanvas->SetGridy(1);
  564.    if(axis == 1)
  565.       gScan[1]->GetXaxis()->SetTitle("X [#mum]");
  566.    else if(axis == 2)
  567.       gScan[1]->GetXaxis()->SetTitle("Y [#mum]");
  568.    gScan[1]->GetXaxis()->SetTitleOffset(1.3);
  569.    gScan[1]->GetXaxis()->CenterTitle(kTRUE);
  570.    gScan[1]->GetXaxis()->SetLabelSize(0.027);
  571.    gScan[1]->GetXaxis()->SetLabelOffset(0.02);
  572.    gScan[1]->GetXaxis()->SetNoExponent(kTRUE);
  573.    gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral");
  574.  
  575.    gScan[1]->GetYaxis()->CenterTitle(kTRUE);
  576.    gScan[1]->GetYaxis()->SetLabelSize(0.027);
  577.    gScan[1]->GetYaxis()->SetLabelOffset(0.02);
  578.    gScan[1]->GetYaxis()->SetRangeUser(0,1);
  579.    gScan[1]->GetYaxis()->SetTitleOffset(1.4);
  580.    gScan[1]->GetYaxis()->SetTitleSize(0.030);
  581.  
  582.    if(!cleanPlots)
  583.    {
  584.       if(axis == 1)
  585.       {
  586.          if(posUnits->widgetCB->GetSelected() == 0)
  587.             gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
  588.          else if(posUnits->widgetCB->GetSelected() == 1)
  589.             gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
  590.       }
  591.       else if(axis == 2)
  592.       {
  593.          if(posUnits->widgetCB->GetSelected() == 0)
  594.             gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
  595.          else if(posUnits->widgetCB->GetSelected() == 1)
  596.             gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
  597.       }
  598.    }
  599.    else
  600.    {
  601.       if(axis == 1)
  602.       {
  603.          if(posUnits->widgetCB->GetSelected() == 0)
  604.             gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
  605.          else if(posUnits->widgetCB->GetSelected() == 1)
  606.             gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
  607.       }
  608.       else if(axis == 2)
  609.       {
  610.          if(posUnits->widgetCB->GetSelected() == 0)
  611.             gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
  612.          else if(posUnits->widgetCB->GetSelected() == 1)
  613.             gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral");
  614.       }
  615.    }
  616.    gScan[1]->SetLineColor(kBlue);
  617.    gScan[0]->SetLineWidth(2);
  618.    gScan[1]->SetLineWidth(2);
  619.  
  620.    gCanvas->Modified();
  621.    gCanvas->Update();
  622.  
  623.    gCanvas->SaveAs(exportname);
  624.  
  625.    // If 2D edge plot, delete the 1D edge plots as we go
  626.    if(edge2d)
  627.    {
  628.       delete gScan[0];
  629.       delete gScan[1];
  630.       if(edit == 0)
  631.          delete gCanvas;
  632.    }
  633.    else
  634.    {
  635.       if(edit == 0)
  636.       {
  637.          delete gScan[0];
  638.          delete gScan[1];
  639.          delete gCanvas;
  640.       }
  641.    }
  642. }
  643.  
  644. void TGAppMainFrame::PhotonMu(TList *files, int edit)
  645. {
  646.    unsigned int nrfiles = files->GetSize();
  647.    char ctemp[1024];
  648.    char exportname[1024];
  649.    int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0;
  650.  
  651.    TCanvas *gCanvas;
  652.    TTree *header_data, *meas_data;
  653.    double *integralCount, *integralPedestal;
  654.    integralCount = new double[nrfiles];
  655.    integralPedestal = new double[nrfiles];
  656.    double *angle, *pdeval, *muval;
  657.    angle = new double[nrfiles];
  658.    pdeval = new double[nrfiles];
  659.    muval = new double[nrfiles];
  660.  
  661.    // Zero the integral count and accumulated vaules
  662.    for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; }
  663.  
  664.    // Fitting variables
  665.    TSpectrum *spec;
  666.    TH1F *histtemp;
  667.    TH1 *histback;
  668.    TH1F *h2;
  669.    float *xpeaks;
  670.    TF1 *fit;
  671.    TF1 *fittingfunc;
  672.    double *fparam, *fparamerr;
  673.    double meansel[20];
  674.    double sigmasel[20];
  675.    double meanparam, paramsigma;
  676.    int sortindex[20];
  677.    int adcpedestal[2];
  678.    int zeromu = 0;
  679.    int darkhist = -1;
  680.  
  681.    double pointest[12];
  682.    bool exclude = false;
  683.  
  684.    // Zero the parameter values
  685.    for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
  686.  
  687.    float progVal = 0;
  688.    analysisProgress->widgetPB->SetPosition(progVal);
  689.    gVirtualX->Update(1);
  690.  
  691.    // Start if we select at least one file
  692.    if(nrfiles > 0)
  693.    {
  694.       for(int i = 0; i < (int)nrfiles; i++)
  695.       {
  696.          if(files->At(i))
  697.          {
  698.             if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
  699.             {
  700.                printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
  701.                darkhist = i;
  702.             }
  703.  
  704.             // Replot the spectrum on analysisCanvas and do not close the input file
  705.             DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
  706.             analysisCanvas->GetCanvas()->Modified();
  707.             analysisCanvas->GetCanvas()->Update();
  708.        
  709.             // Get the spectrum
  710.             histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
  711.             npeaks = 15;
  712.             double par[300];
  713.             spec = new TSpectrum(npeaks);
  714.             // Find spectrum background
  715.             histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
  716.             // Clone histogram and subtract background from it if we select that option
  717.             h2 = (TH1F*)histtemp->Clone("h2");
  718.             if(fitChecks->widgetChBox[0]->IsDown())
  719.                h2->Add(histback, -1);
  720.             // Search for the peaks
  721.             int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
  722.             printf("PhotonMu(): Found %d candidates to fit.\n",found);
  723.             npeaks = found;
  724.    
  725.             // Set initial peak parameters
  726.             xpeaks = spec->GetPositionX();
  727.             for(j = 0; j < found; j++)
  728.             {
  729.                float xp = xpeaks[j];
  730.                int bin = h2->GetXaxis()->FindBin(xp);
  731.                float yp = h2->GetBinContent(bin);
  732.                par[3*j] = yp;
  733.                par[3*j+1] = xp;
  734.                par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
  735.             }
  736.          
  737.             // Fit the histogram
  738.             fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
  739.             TVirtualFitter::Fitter(histtemp, 3*npeaks);
  740.             fit->SetParameters(par);
  741.             fit->SetNpx(300);
  742.             h2->Fit("fit","Q");
  743.             // Get the fitted parameters
  744.             fittingfunc = h2->GetFunction("fit");
  745.             fparam = fittingfunc->GetParameters();
  746.             fparamerr = fittingfunc->GetParErrors();
  747.    
  748.             // Gather the parameters (mean peak value for now)
  749.             int j = 1;
  750.             int nrfit = 0;
  751.             while(1)
  752.             {
  753.                if( (fparam[j] < 1.E-30) || (nrfit > 8) )
  754.                   break;
  755.                else
  756.                {
  757.                   // Check if pedestal is above the lower limit and sigma is smaller than the mean
  758.                   if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
  759.                   {
  760.                      // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
  761.                      meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
  762.                      sigmasel[nrfit] = fparam[j+1];
  763.                      nrfit++;
  764.                   }
  765.                }
  766.    
  767.                j+=3;
  768.             }
  769.             TMath::Sort(nrfit, meansel, sortindex, kFALSE);
  770.  
  771.             meanparam = meansel[sortindex[0]];
  772.             paramsigma = sigmasel[sortindex[0]];
  773.  
  774.             for(j = 0; j < nrfit; j++)
  775.                if(DBGSIG)
  776.                   printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
  777.        
  778.             j = 0;
  779.             adcpedestal[0] = 0;
  780.             adcpedestal[1] = -1;
  781.  
  782.             while(1)
  783.             {
  784.                int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
  785.                int yp = histtemp->GetBinContent(bin);
  786.  
  787.                // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak)
  788.                if(adcpedestal[1] == -1)
  789.                {
  790.                   adcpedestal[0] = j+meanparam+paramsigma;
  791.                   adcpedestal[1] = yp;
  792.                }
  793.                else
  794.                {
  795.                   if( (npeaks > 1) && (adcpedestal[1] >= yp) )
  796.                   {
  797.                      adcpedestal[0] = j+meanparam+paramsigma;
  798.                      adcpedestal[1] = yp;
  799.                   }
  800.                   else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
  801.                   {
  802.                      adcpedestal[0] = j+meanparam+paramsigma;
  803.                      adcpedestal[1] = yp;
  804.                   }
  805.                   else
  806.                      break;
  807.                }
  808.        
  809.                j++;
  810.                if(j > 50) break;
  811.             }
  812.  
  813.             if(npeaks > 1)
  814.             {
  815.                int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
  816.                adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
  817.                printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
  818.                adcpedestal[1] = histtemp->GetBinContent(bin);
  819.             }
  820.  
  821.             if(midPeak->widgetChBox[0]->IsDown())
  822.             {
  823.                if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
  824.                   m = TMath::Floor(meanparam);
  825.                else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
  826.                   m = TMath::Ceil(meanparam);
  827.                int bin = histtemp->GetXaxis()->FindBin(m);
  828.                adcpedestal[0] = m;
  829.                printf("midpeak x = %d, ", adcpedestal[0]);
  830.                adcpedestal[1] = histtemp->GetBinContent(bin);
  831.             }
  832.  
  833. /*          // Option to show the fit
  834.             fittingfunc->Draw("L SAME");
  835.             analysisCanvas->GetCanvas()->Modified();
  836.             analysisCanvas->GetCanvas()->Update();*/
  837.        
  838.             printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
  839.  
  840.             // Delete the opened histogram and spectrum
  841.             delete spec;
  842.             delete inroot;
  843.  
  844.             // Open the input file and read header, ADC and TDC values
  845.             sprintf(ctemp, "%s", files->At(i)->GetTitle());
  846.             inroot = new TFile(ctemp, "READ");
  847.          
  848.             inroot->GetObject("header_data", header_data);
  849.             inroot->GetObject("meas_data", meas_data);
  850.          
  851.             // Reading the header
  852.             if( header_data->FindBranch("angle") )
  853.             {
  854.                header_data->SetBranchAddress("angle", &evtheader.angle);
  855.                header_data->GetEntry(0);
  856.             }
  857.             else
  858.             {
  859.                printf("PhotonMu(): Error! Selected file has no angle header value. Please edit header to add the angle header value.\n");
  860.                break;
  861.             }
  862.    
  863.             char rdc[256];
  864.             j = selectCh->widgetNE[0]->GetNumber();
  865.             double rangetdc[2];
  866.             rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
  867.             rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
  868.    
  869.             k = 0;
  870.             k2 = 0;
  871.             m = 0;
  872.             m2 = 0;
  873.  
  874.             // Reading the data
  875.             for(int e = 0; e < meas_data->GetEntries(); e++)
  876.             {
  877.                sprintf(rdc, "ADC%d", j);
  878.                meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
  879.                meas_data->GetEntry(e);
  880.          
  881.                sprintf(rdc, "TDC%d", j);
  882.                meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
  883.                meas_data->GetEntry(e);
  884.    
  885.                // If our data point is inside the TDC window
  886.                if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
  887.                {
  888.                   // Gather only the integral of the pedestal
  889.                   if((double)evtdata.adcdata[j] < (double)adcpedestal[0]+0.5 )
  890.                   {
  891.                      k2++;
  892.                      m2 += evtdata.adcdata[j];
  893.                   }
  894.  
  895.                   // Gather the complete integral
  896.                   k++;
  897.                   m += evtdata.adcdata[j];
  898.                }
  899.             }
  900.  
  901.             // Determine the angle, mu and relative PDE
  902.             angle[i] = (double)(evtheader.angle);
  903.             integralCount[i] += (double)m;
  904.             printf("PhotonMu(): %lf: Complete integral (%d evts) = %lf\n", angle[i], k, integralCount[i]);
  905.             integralPedestal[i] += (double)m2;
  906.             printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]);
  907.             if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) )
  908.                zeromu = i;
  909.  
  910.             muval[i] = -TMath::Log((double)k2/(double)k);
  911.             printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]);
  912.  
  913.             inroot->Close();
  914.             delete inroot;
  915.          }
  916.  
  917.          // Update the progress bar
  918.          progVal = (float)(90.00/nrfiles)*i;
  919.          analysisProgress->widgetPB->SetPosition(progVal);
  920.          gVirtualX->Update(1);
  921.       }
  922.  
  923.       printf("PhotonMu(): %d files were selected.\n", nrfiles);
  924.  
  925.       printf("PhotonMu(): angle\tmu\trelative PDE\n");
  926.       m = 0;
  927.      
  928.       // Set the 0 degree muval, reuse meansel[1]
  929.       meansel[1] = muval[zeromu];
  930.       printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]);
  931.  
  932.       // TODO - point estimation still not working correctly!
  933.       for(int i = 0; i < (int)nrfiles; i++)
  934.       {
  935.          // Estimate next point and check error (5 point least square fit estimation)
  936.          if( ((i > 4) && (m == i)) || ((i > 5) && (m < i)) )
  937.          {
  938.             // Set exclude signal to false
  939.             exclude = false;
  940.  
  941.             // Get next point values (if zero value -> need to add the dark hist value again)
  942.             pointest[10] = angle[i];
  943.             pointest[11] = muval[i];
  944.  
  945.             // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0]
  946.             meansel[0] = PointEstimate(5, pointest);    // PointEstimate only works with very small step size
  947.             if(meansel[0] > accError->widgetNE[0]->GetNumber())
  948.             {
  949.                printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]);
  950.                exclude = true;
  951.             }
  952.  
  953.             // Value with 0 angle and dark histogram are always needed, so should not be excluded
  954.             if(i == darkhist)
  955.                exclude = false;
  956.  
  957.             // If nothing excluded, pass the points in pointest variable like in a FIFO
  958.             if(!exclude)
  959.             {
  960.                for(int j = 0; j < 10; j++)
  961.                {
  962.                   if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]);
  963.                   pointest[j] = pointest[j+2];
  964.                }
  965.             }
  966.             else
  967.             {
  968.                for(int j = 0; j < 10; j++)
  969.                   if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]);
  970.             }
  971.          }
  972.          else
  973.          {
  974.             // First 5 points act as estimator points for next one
  975.             pointest[2*m] = angle[i];
  976.             pointest[2*m+1] = muval[i];
  977.          }
  978.  
  979.          // Run only if we have a dark run histogram and middle pedestal peak estimation
  980.          if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() )
  981.          {
  982.             if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]);
  983.  
  984.             // Subtract the dark value from all values
  985.             angle[m] = angle[i];
  986.             muval[m] = muval[i] - muval[darkhist];
  987.  
  988.             if(DBGSIG) printf("angle = %lf, newmuval = %lf, darkmuval = %lf, ", angle[m], muval[m], muval[darkhist]);
  989.  
  990.             // Calculate relative PDE
  991. //          pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
  992.             pdeval[m] = muval[m]/((meansel[1]-muval[darkhist])*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
  993.             if(DBGSIG) printf("pdeval = %lf\n", pdeval[m]);
  994.  
  995.             // Only increase counter if error is not too big
  996.             if( (darkhist != i) && (!exclude) )
  997.                m++;
  998.          }
  999.          else
  1000.          {
  1001.             // Relative PDE calculation
  1002.             angle[m] = angle[i];
  1003.             muval[m] = muval[i];
  1004. //            pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
  1005.             pdeval[m] = muval[m]/(meansel[1]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
  1006.  
  1007.             // Only increase counter if error is not too big
  1008.             if(!exclude)
  1009.                m++;
  1010.          }
  1011.          printf("PhotonMu(): %lf\t%lf\t%lf\n", angle[i], muval[i], pdeval[i]);
  1012.       }
  1013.  
  1014.       if(DBGSIG) printf("\n");
  1015.       if(darkhist != -1)
  1016.          printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-1-m));
  1017.       else
  1018.          printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-m));
  1019.  
  1020.       // Plot mu and PDE angle dependance plots
  1021.       if(edit == 0)
  1022.          gCanvas = new TCanvas("canv","canv",1200,900);
  1023.       else if(edit == 1)
  1024.          gCanvas = tempAnalysisCanvas->GetCanvas();
  1025.       gCanvas->cd();
  1026.       gCanvas->SetGrid();
  1027.  
  1028.       TGraph *pde = new TGraph(m, angle, pdeval);
  1029.       pde->SetMarkerStyle(21);
  1030.       pde->SetMarkerSize(0.7);
  1031.       pde->SetMarkerColor(2);
  1032.       pde->SetLineWidth(1);
  1033.       pde->SetLineColor(2);
  1034.       pde->GetXaxis()->SetLabelSize(0.030);
  1035.       pde->GetXaxis()->CenterTitle();
  1036. //      pde->GetXaxis()->SetRange(angle[0],angle[nrfiles-1]);
  1037. //      pde->GetXaxis()->SetRangeUser(angle[0],angle[nrfiles-1]);
  1038.       pde->GetXaxis()->SetRange(-90,90);
  1039.       pde->GetXaxis()->SetRangeUser(-90,90);
  1040.       pde->GetXaxis()->SetLimits(-90,90);
  1041.       pde->GetYaxis()->SetTitleOffset(1.2);
  1042.       pde->GetYaxis()->SetLabelSize(0.030);
  1043.       pde->GetYaxis()->CenterTitle();
  1044.       pde->GetYaxis()->SetRangeUser(0., 1.2);
  1045.       pde->SetName("pde");
  1046.       pde->Draw("ALP");
  1047.  
  1048.       pde->SetTitle(";Incidence angle (#circ);Relative PDE(#theta) / #mu(#theta)");
  1049.  
  1050.       TGraph *mugr = new TGraph(m, angle, muval);
  1051.       mugr->SetMarkerStyle(20);
  1052.       mugr->SetMarkerSize(0.7);
  1053.       mugr->SetMarkerColor(4);
  1054.       mugr->SetLineWidth(1);
  1055.       mugr->SetLineColor(4);
  1056.       mugr->SetName("muval");
  1057.       mugr->Draw("SAME;LP");
  1058.  
  1059.       gCanvas->Modified();
  1060.       gCanvas->Update();
  1061.  
  1062.       if(edit == 0)
  1063.       {
  1064.          remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
  1065.          sprintf(exportname, "%s_pde.pdf", ctemp);
  1066.          gCanvas->SaveAs(exportname);
  1067.  
  1068.          delete mugr;
  1069.          delete pde;
  1070.          delete gCanvas;
  1071.       }
  1072.       else if(edit == 1)
  1073.       {
  1074.          gCanvas->Modified();
  1075.          gCanvas->Update();
  1076.       }
  1077.  
  1078.       // Update the progress bar
  1079.       analysisProgress->widgetPB->SetPosition(100.);
  1080.       gVirtualX->Update(1);
  1081.    }
  1082. }
  1083.  
  1084. void TGAppMainFrame::BreakdownVolt(TList *files, int edit)
  1085. {
  1086.    unsigned int nrfiles = files->GetSize();
  1087.    char ctemp[1024];
  1088.    char exportname[1024];
  1089.    char paramname[1024];
  1090.    int j, k = 0;
  1091.  
  1092.    TCanvas *gCanvas;
  1093.    TTree *header_data, *meas_data;
  1094.  
  1095.    // Fitting variables
  1096.    TSpectrum *spec;
  1097.    TH1F *histtemp;
  1098.    TH1 *histback;
  1099.    TH1F *h2;
  1100.    float *xpeaks;
  1101.    TF1 *fit;
  1102.    TF1 *fittingfunc;
  1103.    TLatex *latex;
  1104.    double *fparam, *fparamerr;
  1105.    double meansel[20];
  1106.    double meanselerr[20];
  1107.    double sigmasel[20];
  1108.    double meanparam, meanparamerr, paramsigma;
  1109.    int sortindex[20];
  1110.    bool exclude = false;
  1111.  
  1112.    int p = 0;
  1113.    double volt[nrfiles], volterr[nrfiles], sep[3][nrfiles], seperr[3][nrfiles];
  1114.    int first = 1;
  1115.  
  1116.    FILE *fp;
  1117.    remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
  1118.    sprintf(paramname, "%s_fitresult.txt", ctemp);
  1119.    fp = fopen(paramname, "w");
  1120.    fclose(fp);
  1121.  
  1122.    int peaklimit = minPeak->widgetNE[0]->GetNumber()+1; // +1 to account for the pedestal peak
  1123.  
  1124.    // Zero the parameter values
  1125.    for(int i = 0; i < nrfiles; i++)
  1126.    {
  1127.       volt[i] = 0; volterr[i] = 0;
  1128.       sep[0][i] = 0; sep[1][i] = 0; sep[2][i] = 0;
  1129.       seperr[0][i] = 0; seperr[1][i] = 0; seperr[2][i] = 0;
  1130.       if(i < 20) { meansel[i] = 0; meanselerr[i] = 0; sigmasel[i] = 0; }
  1131.    }
  1132.  
  1133.    float progVal = 0;
  1134.    analysisProgress->widgetPB->SetPosition(progVal);
  1135.    gVirtualX->Update(1);
  1136.  
  1137.    // Start if we select at least one file
  1138.    if(nrfiles > 0)
  1139.    {
  1140.       for(int i = 0; i < (int)nrfiles; i++)
  1141.       {
  1142.          if(files->At(i))
  1143.          {
  1144.             // Replot the spectrum on analysisCanvas and do not close the input file
  1145.             DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
  1146.             analysisCanvas->GetCanvas()->Modified();
  1147.             analysisCanvas->GetCanvas()->Update();
  1148.        
  1149.             // Get the spectrum
  1150.             histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
  1151.             npeaks = 15;
  1152.             double par[300];
  1153.             spec = new TSpectrum(npeaks);
  1154.             // Find spectrum background
  1155.             histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
  1156.             // Clone histogram and subtract background from it if we select that option
  1157.             h2 = (TH1F*)histtemp->Clone("h2");
  1158.             if(fitChecks->widgetChBox[0]->IsDown())
  1159.                h2->Add(histback, -1);
  1160.             // Search for the peaks
  1161.             int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
  1162.             printf("BreakdownVolt(): Found %d candidates to fit.\n",found);
  1163.             npeaks = found;
  1164.    
  1165.             // Set initial peak parameters
  1166.             xpeaks = spec->GetPositionX();
  1167.             for(j = 0; j < found; j++)
  1168.             {
  1169.                float xp = xpeaks[j];
  1170.                int bin = h2->GetXaxis()->FindBin(xp);
  1171.                float yp = h2->GetBinContent(bin);
  1172.                par[3*j] = yp;
  1173.                par[3*j+1] = xp;
  1174.                par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
  1175.             }
  1176.          
  1177.             // Fit the histogram
  1178.             fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
  1179.             TVirtualFitter::Fitter(histtemp, 3*npeaks);
  1180.             fit->SetParameters(par);
  1181.             fit->SetNpx(300);
  1182.             h2->Fit("fit","Q");
  1183.             // Get the fitted parameters
  1184.             fittingfunc = h2->GetFunction("fit");
  1185.             if(nrfiles == 1)    // TODO: Show the fit if only one file is selected
  1186.                fittingfunc->Draw("SAME");
  1187.             fparam = fittingfunc->GetParameters();
  1188.             fparamerr = fittingfunc->GetParErrors();
  1189.    
  1190.             // Gather the parameters (mean peak value for now)
  1191.             int j = 1;
  1192.             int nrfit = 0;
  1193.             while(1)
  1194.             {
  1195.                if( (fparam[j] < 1.E-30) || (fparamerr[j] < 1.E-10) )// TODO: Maybe not correct for the error
  1196.                   break;
  1197.                else
  1198.                {
  1199.                   // Check if pedestal is above the lower limit and sigma is smaller than the mean
  1200.                   if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
  1201.                   {
  1202.                      meansel[nrfit] = fparam[j];
  1203.                      meanselerr[nrfit] = fparamerr[j];
  1204.                      sigmasel[nrfit] = fparam[j+1];
  1205.                      nrfit++;
  1206.                   }
  1207.                }
  1208.    
  1209.                j+=3;
  1210.             }
  1211.             TMath::Sort(nrfit, meansel, sortindex, kFALSE);
  1212.  
  1213.             meanparam = meansel[sortindex[0]];
  1214.             meanparamerr = meanselerr[sortindex[0]];
  1215.             paramsigma = sigmasel[sortindex[0]];
  1216.  
  1217.             for(j = 0; j < nrfit; j++)
  1218.             {
  1219.                if(DBGSIG)
  1220.                   printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]);
  1221.             }
  1222.  
  1223.             // Delete the opened histogram and spectrum
  1224.             delete spec;
  1225.             delete inroot;
  1226.  
  1227.             // Open the input file and read header, ADC and TDC values
  1228.             sprintf(ctemp, "%s", files->At(i)->GetTitle());
  1229.             inroot = new TFile(ctemp, "READ");
  1230.          
  1231.             inroot->GetObject("header_data", header_data);
  1232.             inroot->GetObject("meas_data", meas_data);
  1233.          
  1234.             // Reading the header
  1235.             if( header_data->FindBranch("biasvolt") )
  1236.             {
  1237.                header_data->SetBranchAddress("biasvolt", &evtheader.biasvolt);
  1238.                header_data->GetEntry(0);
  1239.             }
  1240.             else
  1241.             {
  1242.                printf("BreakdownVolt(): Error! Selected file has no bias voltage header value. Please edit header to add the bias voltage header value.\n");
  1243.                break;
  1244.             }
  1245.  
  1246. //            h2->SetStats(0);
  1247.  
  1248.             analysisCanvas->GetCanvas()->Modified();
  1249.             analysisCanvas->GetCanvas()->Update();
  1250.    
  1251.             // Save each fitting plot
  1252.             if(fitChecks->widgetChBox[1]->IsDown())     // TODO: Check if this works
  1253.             {
  1254.                remove_ext((char*)files->At(i)->GetTitle(), ctemp);
  1255.                sprintf(exportname, "%s_fit.pdf", ctemp);
  1256.                analysisCanvas->GetCanvas()->SaveAs(exportname);
  1257.             }
  1258.  
  1259.             // Calculate the separation between peaks
  1260.             volt[p] = evtheader.biasvolt;
  1261.             volterr[p] = 1.e-5;
  1262.  
  1263.             if(nrfit == 3)
  1264.             {
  1265.                sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
  1266.                seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
  1267.  
  1268.                exclude = (seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber());
  1269.             }
  1270.             else if(nrfit == 4)
  1271.             {
  1272.                sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
  1273.                sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]];
  1274.                seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
  1275.                seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]);
  1276.  
  1277.                exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber()));
  1278.             }
  1279.             else if(nrfit > 4)
  1280.             {
  1281.                sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
  1282.                sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]];
  1283.                sep[2][p] = meansel[sortindex[4]] - meansel[sortindex[3]];
  1284.                seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
  1285.                seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]);
  1286.                seperr[2][p] = TMath::Abs(meanselerr[sortindex[4]]) + TMath::Abs(meanselerr[sortindex[3]]);
  1287.  
  1288.                exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber()) && (seperr[2][p]/sep[2][p] < accError->widgetNE[0]->GetNumber()));
  1289.             }
  1290.             else
  1291.             {
  1292.                printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n");
  1293.                exclude = false;
  1294.             }
  1295.  
  1296.             // Write out parameters to a file
  1297.             fp = fopen(paramname, "a");
  1298.             if(exclude)
  1299.             {
  1300.                if(first == 1)
  1301.                {
  1302.                   fprintf(fp, "%le\t%d\t", evtheader.biasvolt, nrfit);
  1303.                   for(j = 0; j < nrfit; j++)
  1304.                      fprintf(fp, "%le\t%le\t", meansel[sortindex[j]], meanselerr[sortindex[j]]);
  1305.                   fprintf(fp,"\n");
  1306.                   first = 0;
  1307.                }
  1308.                p++;
  1309.             }
  1310.             else
  1311.             {
  1312.                if(nrfit == 3)
  1313.                   printf("Point (at %.2lfV) rejected due to too large errors: %lf\n", volt[p], seperr[0][p]/sep[0][p]);
  1314.                else if(nrfit == 4)
  1315.                   printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p]);
  1316.                else if(nrfit > 4)
  1317.                   printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p], seperr[2][p]/sep[2][p]);
  1318.             }
  1319.             fclose(fp);
  1320.          }
  1321.  
  1322.          if(nrfiles == 1) break;
  1323.          first = 1;
  1324.  
  1325.          // Update the progress bar
  1326.          progVal = (float)(90.00/nrfiles)*i;
  1327.          analysisProgress->widgetPB->SetPosition(progVal);
  1328.          gVirtualX->Update(1);
  1329.       }
  1330.  
  1331.       printf("BreakdownVolt(): %d files were selected.\n", nrfiles);
  1332.       printf("BreakdownVolt(): Number of points to plot is %d\n", p);
  1333.          
  1334.       // Plot and fit breakdown voltage plot
  1335.       if(edit == 0)
  1336.          gCanvas = new TCanvas("canv","canv",1200,900);
  1337.       else if(edit == 1)
  1338.          gCanvas = tempAnalysisCanvas->GetCanvas();
  1339.       gCanvas->cd();
  1340.       gCanvas->SetGrid();
  1341.  
  1342.       TGraphErrors* bdplot;
  1343.       k = peakSepCalc->widgetNE[0]->GetNumber();
  1344.       if(k < 4)
  1345.          bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]);
  1346.       else
  1347.       {
  1348.          printf("BreakdownVold(): Unsupported peak separation selected (%d).\n", k);
  1349.          return;
  1350.       }
  1351.  
  1352.       bdplot->SetMarkerStyle(20);
  1353.       bdplot->SetMarkerSize(0.4);
  1354.       bdplot->SetMarkerColor(kBlue);
  1355.       bdplot->SetLineWidth(1);
  1356.       bdplot->SetLineColor(kBlue);
  1357.       bdplot->GetXaxis()->SetLabelSize(0.030);
  1358.       bdplot->GetXaxis()->CenterTitle();
  1359. //      bdplot->GetXaxis()->SetRange(-90,90);
  1360. //      bdplot->GetXaxis()->SetRangeUser(-90,90);
  1361. //      bdplot->GetXaxis()->SetLimits(-90,90);
  1362.       bdplot->GetYaxis()->SetTitleOffset(1.2);
  1363.       bdplot->GetYaxis()->SetLabelSize(0.030);
  1364.       bdplot->GetYaxis()->CenterTitle();
  1365. //      bdplot->GetYaxis()->SetRangeUser(0., 1.2);
  1366.       bdplot->SetName("bdplot");
  1367.       bdplot->Draw("AP");
  1368.       bdplot->SetTitle(";Bias voltage (V);Peak separation");
  1369.  
  1370.       // Fit the breakdown voltage plot with a line
  1371.       bdplot->Fit("pol1","Q");
  1372.  
  1373.       TF1 *fit1 = bdplot->GetFunction("pol1");
  1374.       meansel[0] = fit1->GetParameter(0);
  1375.       meanselerr[0] = fit1->GetParError(0);
  1376.       meansel[1] = fit1->GetParameter(1);
  1377.       meanselerr[1] = fit1->GetParError(1);
  1378.  
  1379.       meansel[2] = -meansel[0]/meansel[1];
  1380.       if(!cleanPlots)
  1381.       {
  1382.          sprintf(ctemp, "#splitline{#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf}", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
  1383.          latex = new TLatex();
  1384.          latex->SetTextSize(0.039);
  1385.          latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[p-1]], ctemp);
  1386.       }
  1387.       else
  1388.          printf("#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
  1389.  
  1390.       if(edit == 0)
  1391.       {
  1392.          remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
  1393.          sprintf(exportname, "%s_breakdown.pdf", ctemp);
  1394.          gCanvas->SaveAs(exportname);
  1395.  
  1396.          delete bdplot;
  1397.          delete gCanvas;
  1398.       }
  1399.       else if(edit == 1)
  1400.       {
  1401.          gCanvas->Modified();
  1402.          gCanvas->Update();
  1403.       }
  1404.  
  1405.       // Update the progress bar
  1406.       analysisProgress->widgetPB->SetPosition(100.);
  1407.       gVirtualX->Update(1);
  1408.    }
  1409. }
  1410.  
  1411. void TGAppMainFrame::SurfaceScan(TList *files, int edit)
  1412. {
  1413.    unsigned int nrfiles = files->GetSize();
  1414.    char ctemp[1024];
  1415.    char exportname[1024];
  1416.    int j, k = 0, m = 0, n = 0;
  1417.  
  1418.    TTree *header_data, *meas_data;
  1419.    double *integralCount;
  1420.    integralCount = new double[nrfiles];
  1421.    double *surfx, *surfy;
  1422.    surfx = new double[nrfiles];
  1423.    surfy = new double[nrfiles];
  1424.    double xsurfmin = 0, ysurfmin = 0;
  1425.    int nrentries;
  1426. //   double minInteg, maxInteg;
  1427.    bool norm = surfScanOpt->widgetChBox[0]->IsDown();
  1428. //   double curzval;
  1429. //   bool edge2d = false;
  1430.    
  1431.    TCanvas *gCanvas;
  1432.  
  1433.    float progVal = 0;
  1434.    analysisProgress->widgetPB->SetPosition(progVal);
  1435.    gVirtualX->Update(1);
  1436.  
  1437.    // Zero the integral count and accumulated vaules
  1438.    for(int i = 0; i < (int)nrfiles; i++) integralCount[i] = 0;
  1439.  
  1440.    // Start if we select at more than one file
  1441.    if( (nrfiles > 1) && (multiSelect->widgetChBox[0]->IsOn()) )
  1442.    {
  1443.       printf("SurfaceScan(): Creating a surface plot. Please wait...\n");
  1444.       for(int i = 0; i < (int)nrfiles; i++)
  1445.       {
  1446.          if(files->At(i))
  1447.          {
  1448.             n++;
  1449.             // Read the X and Y positions from header and ADC and TDC values from the measurements
  1450.             sprintf(ctemp, "%s", files->At(i)->GetTitle());
  1451.             inroot = new TFile(ctemp, "READ");
  1452.          
  1453.             inroot->GetObject("header_data", header_data);
  1454.             inroot->GetObject("meas_data", meas_data);
  1455.          
  1456.             header_data->SetBranchAddress("xpos", &evtheader.xpos);
  1457.             header_data->GetEntry(0);
  1458.             header_data->SetBranchAddress("ypos", &evtheader.ypos);
  1459.             header_data->GetEntry(0);
  1460.  
  1461.             char rdc[256];
  1462.             j = selectCh->widgetNE[0]->GetNumber();
  1463.             double rangetdc[2];
  1464.             rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
  1465.             rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
  1466.    
  1467.             k = 0;
  1468.             m = 0;
  1469.          
  1470.             // Reading the data
  1471.             for(int e = 0; e < meas_data->GetEntries(); e++)
  1472.             {
  1473.                sprintf(rdc, "ADC%d", j);
  1474.                meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
  1475.                meas_data->GetEntry(e);
  1476.          
  1477.                sprintf(rdc, "TDC%d", j);
  1478.                meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
  1479.                meas_data->GetEntry(e);
  1480.    
  1481.                // Use data point only if it is inside the TDC window
  1482.                if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
  1483.                {
  1484.                   k++;
  1485.                   m += evtdata.adcdata[j];
  1486.                }
  1487.             }
  1488.  
  1489.             // X, Y and Z values from each file (table units or microns)
  1490.             if(posUnits->widgetCB->GetSelected() == 0)
  1491.             {
  1492.                if(n == 1)
  1493.                {
  1494.                   xsurfmin = (double)(evtheader.xpos);
  1495.                   ysurfmin = (double)(evtheader.ypos);
  1496.                }
  1497.  
  1498.                surfx[i] = (double)(evtheader.xpos);
  1499.                surfy[i] = (double)(evtheader.ypos);
  1500.             }
  1501.             else if(posUnits->widgetCB->GetSelected() == 1)
  1502.             {
  1503.                if(n == 1)
  1504.                {
  1505.                   xsurfmin = (double)(evtheader.xpos*lenconversion);
  1506.                   ysurfmin = (double)(evtheader.ypos*lenconversion);
  1507.                }
  1508.  
  1509.                surfx[i] = (double)(evtheader.xpos*lenconversion);
  1510.                surfy[i] = (double)(evtheader.ypos*lenconversion);
  1511.             }
  1512.  
  1513.             // Save result for later plotting
  1514.             if(norm)
  1515.                integralCount[i] += ((double)m)/((double)k);
  1516.             else
  1517.                integralCount[i] += m;
  1518.            
  1519.             inroot->Close();
  1520.             delete inroot;
  1521.          }
  1522.  
  1523.          // Update the progress bar
  1524.          progVal = (float)(75.00/nrfiles)*i;
  1525.          analysisProgress->widgetPB->SetPosition(progVal);
  1526.          gVirtualX->Update(1);
  1527.       }
  1528.  
  1529.       nrentries = n;
  1530.       printf("SurfaceScan(): %d files were selected.\n", nrfiles);
  1531.  
  1532. //      // If only an integral is needed, do not plot and exit here
  1533. //      if( direction == 0 )
  1534. //      {
  1535. //         delete[] integralCount;
  1536. //         delete[] surfxy;
  1537. //         delete[] surfz;
  1538. //         return;
  1539. //      }
  1540. //
  1541. //      // Current z value and the accumulated counter
  1542. //      curzval = surfz[0];
  1543. //      j = 0;
  1544. //      int acc = 0;
  1545. //      int zb;
  1546. //      for(int i = 0; i <= (int)nrfiles; i++)
  1547. //      {
  1548. //         // Collect the accumulated integral in order to produce a PDF from a CDF
  1549. //         // While we are at the same Z value, save under one set
  1550. //         if( (surfz[i] == curzval) && (acc != nrfiles) )
  1551. //         {
  1552. //            integralAcc[j] = integralCount[i];
  1553. //            if(DBGSIG) printf("IntegSpectrum(): Integral check 1 (i=%d,j=%d,z=%.2lf): %lf\t%lf\n", i, j, surfz[i], integralCount[i], integralAcc[j]);
  1554. //            j++;
  1555. //            acc++;
  1556. //         }
  1557. //         // When we switch to a new set of Z values and at the end, we must save the previous ones to make 1D edge plots
  1558. //         else
  1559. //         {
  1560. //            // Find minimal and maximal integral values to subtract the offset and normate PDF to 1
  1561. //            NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
  1562. //
  1563. //            if(acc != nrfiles)
  1564. //            {
  1565. //               curzval = surfz[i];
  1566. //               // PDF and CDF plot
  1567. //               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
  1568. //               i--;
  1569. //               j = 0;
  1570. //            }
  1571. //            else
  1572. //            {
  1573. //               // PDF and CDF plot
  1574. //               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
  1575. //               i--;
  1576. //               break;
  1577. //            }
  1578. //         }
  1579. //
  1580. //         // Update the progress bar
  1581. //         progVal = (float)(15.00/nrfiles)*i+75.00;
  1582. //         analysisProgress->widgetPB->SetPosition(progVal);
  1583. //         gVirtualX->Update(1);
  1584. //      }
  1585.  
  1586.       // Make the 2D surface plot
  1587.       if(edit == 0)
  1588.          gCanvas = new TCanvas("canv","canv",1100,900);
  1589.       else
  1590.          gCanvas = tempAnalysisCanvas->GetCanvas();
  1591.  
  1592.       double range[4];
  1593.       TGraph2D *gScan2D;
  1594.       gScan2D = new TGraph2D();
  1595.       range[0] = TMath::MinElement(nrfiles, surfx);
  1596.       range[1] = TMath::MaxElement(nrfiles, surfx);
  1597.       range[2] = TMath::MinElement(nrfiles, surfy);
  1598.       range[3] = TMath::MaxElement(nrfiles, surfy);
  1599.  
  1600.       for(int i = 0; i < nrfiles; i++)
  1601.       {
  1602.          if(DBGSIG)
  1603.          {
  1604.             if(surfScanOpt->widgetChBox[1]->IsDown())
  1605.                printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i]-range[0], surfy[i]-range[2], integralCount[i]);
  1606.             else
  1607.                printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i], surfy[i], integralCount[i]);
  1608.          }
  1609.  
  1610.          if(surfScanOpt->widgetChBox[1]->IsDown())
  1611.             gScan2D->SetPoint(i, surfx[i]-range[0], surfy[i]-range[2], integralCount[i]);
  1612.          else
  1613.             gScan2D->SetPoint(i, surfx[i], surfy[i], integralCount[i]);
  1614.  
  1615.          // Update the progress bar
  1616.          progVal = (float)(9.00/nrfiles)*i+90.00;
  1617.          analysisProgress->widgetPB->SetPosition(progVal);
  1618.          gVirtualX->Update(1);
  1619.       }
  1620.  
  1621.       if(surfScanOpt->widgetChBox[1]->IsDown())
  1622.       {
  1623.          range[1] -= range[0];
  1624.          range[3] -= range[2];
  1625.          range[0] -= range[0];
  1626.          range[2] -= range[2];
  1627.       }
  1628.  
  1629.       gCanvas->cd();
  1630.       gStyle->SetPalette(1);
  1631.       gCanvas->SetLeftMargin(0.15);
  1632.       gCanvas->SetRightMargin(0.126);
  1633.       gCanvas->SetTopMargin(0.077);
  1634.       gScan2D->Draw("COLZ");
  1635.      
  1636.       gCanvas->Modified();
  1637.       gCanvas->Update();
  1638.      
  1639.       if(!cleanPlots)
  1640.       {
  1641.          if(posUnits->widgetCB->GetSelected() == 0)
  1642.             gScan2D->SetTitle("Surface scan;X [table units];Y [table units]");
  1643.          else if(posUnits->widgetCB->GetSelected() == 1)
  1644.             gScan2D->SetTitle("Surface scan;X [#mum];Y [#mum]");
  1645.       }
  1646.       else
  1647.       {
  1648.          if(posUnits->widgetCB->GetSelected() == 0)
  1649.             gScan2D->SetTitle(";X [table units];Y [table units]");
  1650.          else if(posUnits->widgetCB->GetSelected() == 1)
  1651.             gScan2D->SetTitle(";X [#mum];Y [#mum]");
  1652.       }
  1653. /*      TGaxis *xax = (TGaxis*)gScan2D->GetXaxis();
  1654.       xax->SetMaxDigits(4);
  1655.       TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
  1656.       yax->SetMaxDigits(4);*/
  1657.  
  1658.       gScan2D->SetNpx((int)resolSurf->widgetNE[0]->GetNumber());
  1659.       gScan2D->SetNpy((int)resolSurf->widgetNE[1]->GetNumber());
  1660.      
  1661.       gCanvas->Modified();
  1662.       gCanvas->Update();
  1663.  
  1664.       gScan2D->GetXaxis()->SetTitleOffset(1.3);
  1665.       gScan2D->GetXaxis()->CenterTitle(kTRUE);
  1666.       gScan2D->GetXaxis()->SetLabelSize(0.027);
  1667.       gScan2D->GetXaxis()->SetLabelOffset(0.02);
  1668.       gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
  1669.       gScan2D->GetXaxis()->SetNoExponent(kTRUE);
  1670.       gScan2D->GetYaxis()->SetTitleOffset(1.9);
  1671.       gScan2D->GetYaxis()->CenterTitle(kTRUE);
  1672.       gScan2D->GetYaxis()->SetLabelSize(0.027);
  1673.       gScan2D->GetXaxis()->SetLabelOffset(0.02);
  1674.       gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
  1675.       gScan2D->GetYaxis()->SetNoExponent(kTRUE);
  1676.      
  1677.       gCanvas->Modified();
  1678.       gCanvas->Update();
  1679.  
  1680.       remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
  1681.       sprintf(exportname, "%s_surfscan.pdf", ctemp);
  1682.       gCanvas->SaveAs(exportname);
  1683.  
  1684.       // Update the progress bar
  1685.       analysisProgress->widgetPB->SetPosition(100.0);
  1686.       gVirtualX->Update(1);
  1687.  
  1688.       if(edit == 0)
  1689.       {
  1690.          delete gScan2D;
  1691.          delete gCanvas;
  1692.       }
  1693.       else if(edit == 1)
  1694.       {
  1695.          gCanvas->Modified();
  1696.          gCanvas->Update();
  1697.       }
  1698.    }
  1699. }
  1700.