Subversion Repositories f9daq

Rev

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