Subversion Repositories f9daq

Rev

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