Subversion Repositories f9daq

Rev

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