Subversion Repositories f9daq

Rev

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