WebSVN – f9daq – Diff – /lab/sipmscan/trunk/src/analysis.cpp

 #include "../include/sipmscan.h"
 #include "../include/workstation.h"
 #include <stdio.h>
 #include <stdlib.h>
 // Peak detection function
 int npeaks = 20;
 double FindPeaks(double *x, double *par)
+{
    double result = 0;
    for(int i = 0; i < npeaks; i++)
+   {
       double norm = par[3*i];
       double mean = par[3*i+1];
       double sigma = par[3*i+2];
       result += norm*TMath::Gaus(x[0], mean, sigma);
+   }
    return result;
+}
 // File browser for selecting the dark run histogram
 void TGAppMainFrame::SelectDarkHist()
+{
    TGFileInfo file_info;
    const char *filetypes[] = {"Histograms",histextall,0,0};
    char *cTemp;
    file_info.fFileTypes = filetypes;
 //   cTemp = new char[1024];
 //   sprintf(cTemp, "%s/results", rootdir);
 //   file_info.fIniDir = StrDup(cTemp);
    file_info.fIniDir = StrDup(currentAnalDir);
    file_info.fMultipleSelection = kFALSE;
    new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
 //   delete[] cTemp;
    if(file_info.fFilename != NULL)
+   {
       darkRun->widgetTE->SetText(file_info.fFilename);
       fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
       fileList->Layout();
+   }
    else
       darkRun->widgetTE->SetText("");
+}
 // Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
 void TGAppMainFrame::AnalysisDefaults()
+{
    printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
    if(analTab->GetCurrent() == 0)       // Integrate spectrum
+   {
       intSpect->widgetChBox[0]->SetState(kButtonUp);
       intSpect->widgetChBox[1]->SetState(kButtonUp);
       intSpect->widgetChBox[2]->SetState(kButtonDown);
       resol2d->widgetNE[0]->SetNumber(40);
       resol2d->widgetNE[1]->SetNumber(40);
+   }
    else if(analTab->GetCurrent() == 1)  // Relative PDE
+   {
       relPde->widgetChBox[0]->SetState(kButtonDown);
+      relPde->widgetChBox[1]->SetState(kButtonUp);
       midPeak->widgetChBox[0]->SetState(kButtonUp);
       zeroAngle->widgetNE[0]->SetNumber(0.00);
+   }
    else if(analTab->GetCurrent() == 2)  // Breakdown voltage
+   {
       minPeak->widgetNE[0]->SetNumber(2);
       minPeak->widgetNE[0]->SetNumber(1);
+   }
    else if(analTab->GetCurrent() == 3)  // Surface scan
+   {
       surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
       surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
       resolSurf->widgetNE[0]->SetNumber(40);
       resolSurf->widgetNE[1]->SetNumber(40);
+   }
+}
 // Analysis handle function
 void TGAppMainFrame::AnalysisHandle(int type)
+{
    TList *files;
    bool createTab = true;
    int tabid = -1;
    int analtab = analTab->GetCurrent();
    int analtype;
    if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
       analtype = 0;
    else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) )
       analtype = 1;
    else if( analtab == 1 )
       analtype = 2;
    else if( analtab == 2 )
       analtype = 3;
    else if( analtab == 3 )
       analtype = 4;
+   // Save analysis settings any time we run a new analysis
+   SaveAnalSettings();
    // Only integrate spectrum or make relative PDE
    if(type == 0)
+   {
       files = new TList();
       fileList->GetSelectedEntries(files);
       if( analtype == 0 )
          IntegSpectrum(files, 0, 0);
       if( intSpect->widgetChBox[0]->IsDown() )
          IntegSpectrum(files, 1, 0);
       if( intSpect->widgetChBox[1]->IsDown() )
          IntegSpectrum(files, 2, 0);
       if( analtype == 2 )
          PhotonMu(files, 0);
       if( analtype == 3 )
          BreakdownVolt(files, 0);
       if( analtype == 4 )
          SurfaceScan(files, 0);
+   }
    // Integrate spectrum or make relative PDE and open edit window
    else if(type == 1)
    {
       files = new TList();
       fileList->GetSelectedEntries(files);
       // Prepare a new analysis edit tab, if we want to edit plots
       for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
+      {
          if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
+         {
             createTab = false;
             tabid = i;
+         }
          if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
+      }
       if(files->GetSize() > 0)
+      {
          TempAnalysisTab(fTab, createTab, &tabid, analtype);
          // Integrate spectra
          if( analtype == 0 )
             IntegSpectrum(files, 0, 1);
          if( intSpect->widgetChBox[0]->IsDown() )
             IntegSpectrum(files, 1, 1);
          if( intSpect->widgetChBox[1]->IsDown() )
             IntegSpectrum(files, 2, 1);
          if( analtype == 2 )
             PhotonMu(files, 1);
          if( analtype == 3 )
             BreakdownVolt(files, 1);
          if( analtype == 4 )
             SurfaceScan(files, 1);
          fTab->SetTab(tabid);
+      }
+   }
    delete files;
+}
 // Analysis functions ----------------------------------
 void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
+{
    unsigned int nrfiles = files->GetSize();
    char ctemp[1024];
    char exportname[1024];
    int j, k = 0, m = 0;
    TTree *header_data, *meas_data;
    double *integralCount, *integralAcc;
    integralCount = new double[nrfiles];
    integralAcc = new double[nrfiles];
    double *surfxy, *surfz;
    surfxy = new double[nrfiles];
    surfz = new double[nrfiles];
    double minInteg, maxInteg;
    bool norm = intSpect->widgetChBox[2]->IsDown();
    double curzval;
    bool edge2d = false;
    TCanvas *gCanvas;
    float progVal = 0;
    analysisProgress->widgetPB->SetPosition(progVal);
    gVirtualX->Update(1);
    // Zero the integral count and accumulated vaules
    for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
    // Start if we select at least one file
    if(nrfiles > 0)
+   {
       for(int i = 0; i < (int)nrfiles; i++)
+      {
          if(files->At(i))
+         {
             // Read the X,Y and Z positions from header and ADC and TDC values from the measurements
             sprintf(ctemp, "%s", files->At(i)->GetTitle());
             inroot = new TFile(ctemp, "READ");
             inroot->GetObject("header_data", header_data);
             inroot->GetObject("meas_data", meas_data);
             header_data->SetBranchAddress("xpos", &evtheader.xpos);
             header_data->GetEntry(0);
             header_data->SetBranchAddress("ypos", &evtheader.ypos);
             header_data->GetEntry(0);
             header_data->SetBranchAddress("zpos", &evtheader.zpos);
             header_data->GetEntry(0);
             char rdc[256];
             j = selectCh->widgetNE[0]->GetNumber();
             double rangetdc[2];
             rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
             rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
             k = 0;
             m = 0;
             // Reading the data
             for(int e = 0; e < meas_data->GetEntries(); e++)
+            {
                sprintf(rdc, "ADC%d", j);
                meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
                meas_data->GetEntry(e);
                sprintf(rdc, "TDC%d", j);
                meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
                meas_data->GetEntry(e);
                // Use data point only if it is inside the TDC window
                if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
+               {
                   k++;
                   m += evtdata.adcdata[j];
+               }
+            }
             // X, Y and Z values from each file (table units or microns)
             if(posUnitsPlot->widgetCB->GetSelected() == 0)
+            {
                if(direction == 1)
                   surfxy[i] = (double)(evtheader.xpos);
                else if(direction == 2)
                   surfxy[i] = (double)(evtheader.ypos);
                surfz[i] = (double)(evtheader.zpos);
+            }
             else if(posUnitsPlot->widgetCB->GetSelected() == 1)
+            {
                if(direction == 1)
                   surfxy[i] = (double)(evtheader.xpos*lenconversion);
                else if(direction == 2)
                   surfxy[i] = (double)(evtheader.ypos*lenconversion);
                surfz[i] = (double)(evtheader.zpos*lenconversion);
+            }
             // 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
             if(i == 0) curzval = surfz[i];
             else
+            {
                if(surfz[i] != curzval)
                   edge2d = true;
+            }
             // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
             if(direction == 0)
+            {
                if(norm)
+               {
                   integralCount[i] += ((double)m)/((double)k);
                   printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
+               }
                else
+               {
                   integralCount[i] += m;
                   printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
+               }
+            }
             else
+            {
                if(norm)
                   integralCount[i] += ((double)m)/((double)k);
                else
                   integralCount[i] += m;
+            }
             inroot->Close();
             delete inroot;
+         }
          // Update the progress bar
          progVal = (float)(75.00/nrfiles)*i;
          analysisProgress->widgetPB->SetPosition(progVal);
          gVirtualX->Update(1);
+      }
       printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
       // If only an integral is needed, do not plot and exit here
       if( direction == 0 )
       {
          delete[] integralCount;
          delete[] surfxy;
          delete[] surfz;
          return;
+      }
       // Current z value and the accumulated counter
       curzval = surfz[0];
       j = 0;
       int acc = 0;
       int zb;
       for(int i = 0; i <= (int)nrfiles; i++)
+      {
          // Collect the accumulated integral in order to produce a PDF from a CDF
          // While we are at the same Z value, save under one set
          if( (surfz[i] == curzval) && (acc != nrfiles) )
+         {
             integralAcc[j] = integralCount[i];
             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]);
             j++;
             acc++;
+         }
          // 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
          else
+         {
             // Find minimal and maximal integral values to subtract the offset and normate PDF to 1
             NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
             if(acc != nrfiles)
+            {
                curzval = surfz[i];
                // PDF and CDF plot
                PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
                i--;
                j = 0;
+            }
             else
+            {
                // PDF and CDF plot
                PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
                i--;
                break;
+            }
+         }
          // Update the progress bar
          progVal = (float)(15.00/nrfiles)*i+75.00;
          analysisProgress->widgetPB->SetPosition(progVal);
          gVirtualX->Update(1);
+      }
       // Make the 2D edge plot
       if(edge2d)
+      {
          if(edit == 0)
             gCanvas = new TCanvas("canv","canv",900,900);
          else
             gCanvas = tempAnalysisCanvas->GetCanvas();
          double range[4];
          TGraph2D *gScan2D;
          gScan2D = new TGraph2D();
-         gScan2D->SetName("edgescan");
          range[0] = TMath::MinElement(nrfiles, surfxy);
          range[1] = TMath::MaxElement(nrfiles, surfxy);
          range[2] = TMath::MinElement(nrfiles, surfz);
          range[3] = TMath::MaxElement(nrfiles, surfz);
          for(int i = 0; i < nrfiles; i++)
+         {
             if(DBGSIG)
                printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
             gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
             // Update the progress bar
             progVal = (float)(9.00/nrfiles)*i+90.00;
             analysisProgress->widgetPB->SetPosition(progVal);
             gVirtualX->Update(1);
+         }
          gCanvas->cd();
          gStyle->SetPalette(1);
          gCanvas->SetLeftMargin(0.15);
          gCanvas->SetRightMargin(0.126);
          gCanvas->SetTopMargin(0.077);
          gScan2D->Draw("COLZ");
          gCanvas->Modified();
          gCanvas->Update();
          gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
          gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
          gCanvas->Modified();
          gCanvas->Update();
          if(direction == 1)
             gScan2D->GetXaxis()->SetTitle("X [#mum]");
          else if(direction == 2)
             gScan2D->GetXaxis()->SetTitle("Y [#mum]");
          gScan2D->GetXaxis()->SetTitleOffset(1.3);
          gScan2D->GetXaxis()->CenterTitle(kTRUE);
          gScan2D->GetXaxis()->SetLabelSize(0.027);
          gScan2D->GetXaxis()->SetLabelOffset(0.02);
          gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
          gScan2D->GetXaxis()->SetNoExponent(kTRUE);
          gScan2D->GetYaxis()->SetTitleOffset(1.9);
          gScan2D->GetYaxis()->CenterTitle(kTRUE);
          gScan2D->GetYaxis()->SetLabelSize(0.027);
          gScan2D->GetYaxis()->SetLabelOffset(0.02);
          gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
          gScan2D->GetYaxis()->SetNoExponent(kTRUE);
 /*         TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
          yax->SetMaxDigits(4);*/
          if(!cleanPlots)
+         {
             if(direction == 1)
+            {
                if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
                else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
+            }
             else if(direction == 2)
+            {
                if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
                else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
+            }
+         }
          else
+         {
             if(direction == 1)
+            {
                if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle(";X [table units];Z [table units]");
                else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle(";X [#mum];Z [#mum]");
+            }
             else if(direction == 2)
+            {
                if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle(";Y [table units];Z [table units]");
                else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle(";Y [#mum];Z [#mum]");
+            }
+         }
          gCanvas->Modified();
          gCanvas->Update();
          remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
          sprintf(exportname, "%s", ctemp);
          remove_from_last(exportname, '_', ctemp);
          if(direction == 1)
             sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
          else if(direction == 2)
             sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
          gCanvas->SaveAs(exportname);
          // Update the progress bar
          analysisProgress->widgetPB->SetPosition(100.0);
          gVirtualX->Update(1);
          if(edit == 0)
+         {
             delete gScan2D;
             delete gCanvas;
+         }
+      }
       else
+      {
          // Update the progress bar
          analysisProgress->widgetPB->SetPosition(100.0);
          gVirtualX->Update(1);
+      }
+   }
+}
 void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit)
+{
    TGraph *gScan[2];
    int pdfmax = -1;
    int count = 0;
    char ctemp[1024];
    char exportname[1024];
    TCanvas *gCanvas;
    // Prepare the CDF plot
    gScan[1] = new TGraph();
    for(int i = 0; i < points; i++)
+   {
       count = filenr - points + i;
       gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max));
       if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", i, (double)xy[count], (double)integAcc[i]/(*max));
       if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) )
          pdfmax = (integAcc[i+1]-integAcc[i])/(*max);
+   }
    pdfmax = (TMath::Ceil(pdfmax*10))/10.;
    // Prepare the PDF plot
    gScan[0] = new TGraph();
    for(int i = points-1; i >= 0; i--)
+   {
       count = (filenr-1) - (points-1) + i;
       // Set any negative values of the PDF to 0
       if( (integAcc[i]-integAcc[i-1])/(*max) < 0 )
          gScan[0]->SetPoint(i, (double)xy[count], 0);
       else
          gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max));
       if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max));
+   }
    remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp);
    sprintf(exportname, "%s_edge.pdf", ctemp);
    if(edit == 0)
       gCanvas = new TCanvas("canv1d","canv1d",1200,900);
    else
       gCanvas = tempAnalysisCanvas->GetCanvas();
    // Fit the PDF with a gaussian
    gScan[0]->Fit("gaus","Q");
    gScan[0]->GetFunction("gaus")->SetNpx(400);
    gStyle->SetOptFit(1);
    gScan[1]->Draw("AL");
    gPad->Update();
    gScan[0]->Draw("LP");
    gCanvas->Modified();
    gCanvas->Update();
    TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats");
    if(!cleanPlots)
+   {
       stats->SetX1NDC(0.86); stats->SetX2NDC(1.0);
       stats->SetY1NDC(0.87); stats->SetY2NDC(1.0);
+   }
    else
+   {
       stats->SetX1NDC(1.1); stats->SetX2NDC(1.3);
       stats->SetY1NDC(1.1); stats->SetY2NDC(1.3);
+   }
    gCanvas->SetGridx(1);
    gCanvas->SetGridy(1);
    if(axis == 1)
       gScan[1]->GetXaxis()->SetTitle("X [#mum]");
    else if(axis == 2)
       gScan[1]->GetXaxis()->SetTitle("Y [#mum]");
    gScan[1]->GetXaxis()->SetTitleOffset(1.3);
    gScan[1]->GetXaxis()->CenterTitle(kTRUE);
    gScan[1]->GetXaxis()->SetLabelSize(0.027);
    gScan[1]->GetXaxis()->SetLabelOffset(0.02);
    gScan[1]->GetXaxis()->SetNoExponent(kTRUE);
    gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral");
    gScan[1]->GetYaxis()->CenterTitle(kTRUE);
    gScan[1]->GetYaxis()->SetLabelSize(0.027);
    gScan[1]->GetYaxis()->SetLabelOffset(0.02);
    gScan[1]->GetYaxis()->SetRangeUser(0,1);
    gScan[1]->GetYaxis()->SetTitleOffset(1.4);
    gScan[1]->GetYaxis()->SetTitleSize(0.030);
    if(!cleanPlots)
+   {
       if(axis == 1)
+      {
          if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
          else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
+      }
       else if(axis == 2)
+      {
          if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
          else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
+      }
+   }
    else
+   {
       if(axis == 1)
+      {
          if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
          else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
+      }
       else if(axis == 2)
+      {
          if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
          else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral");
+      }
+   }
    gScan[1]->SetLineColor(kBlue);
    gScan[0]->SetLineWidth(2);
    gScan[1]->SetLineWidth(2);
    gCanvas->Modified();
    gCanvas->Update();
    gCanvas->SaveAs(exportname);
    // If 2D edge plot, delete the 1D edge plots as we go
    if(edge2d)
+   {
       delete gScan[0];
       delete gScan[1];
       if(edit == 0)
          delete gCanvas;
+   }
    else
+   {
       if(edit == 0)
+      {
          delete gScan[0];
          delete gScan[1];
          delete gCanvas;
+      }
+   }
+}
 void TGAppMainFrame::PhotonMu(TList *files, int edit)
+{
    unsigned int nrfiles = files->GetSize();
    char ctemp[1024];
    char exportname[1024];
    int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0;
    TCanvas *gCanvas;
    TTree *header_data, *meas_data;
    double *integralCount, *integralPedestal;
    integralCount = new double[nrfiles];
    integralPedestal = new double[nrfiles];
    double *angle, *pdeval, *muval;
    angle = new double[nrfiles];
    pdeval = new double[nrfiles];
    muval = new double[nrfiles];
    // Zero the integral count and accumulated vaules
    for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; }
    // Fitting variables
    TSpectrum *spec;
    TH1F *histtemp;
    TH1 *histback;
    TH1F *h2;
    float *xpeaks;
    TF1 *fit;
    TF1 *fittingfunc;
    double *fparam, *fparamerr;
    double meansel[20];
    double sigmasel[20];
    double meanparam, paramsigma;
    int sortindex[20];
    int adcpedestal[2];
    int zeromu = 0;
    int darkhist = -1;
+   int nopeaks = -1;
    double pointest[12];
    bool exclude = false;
    // Zero the parameter values
    for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
    float progVal = 0;
    analysisProgress->widgetPB->SetPosition(progVal);
    gVirtualX->Update(1);
+   // Check if the checkbox for no peaks is selected - TODO: Still need the situation when we do not have a peaked ADC spectrum
+   if(relPde->widgetChBox[1]->IsDown())
+   {
+      printf("PhotonMu(): ADC spectrum has no peak structure.\n");
+      nopeaks = 1;
+      // Error if there is no darkhist
+      if(strcmp("", darkRun->widgetTE->GetText()) == 0)
+      {
+         printf("PhotonMu(): Error! The no peak structure option needs a dark histogram.\n");
+         delete[] integralCount;
+         delete[] integralPedestal;
+         delete[] angle;
+         delete[] pdeval;
+         delete[] muval;
+         return;
+      }
+   }
    // Start if we select at least one file
    if(nrfiles > 0)
    {
+      // Find the pedestal peak for the dark histogram, and use it for all if there are no peaks
+      if(nopeaks != -1)
+      {
+         // Replot the spectrum on analysisCanvas and do not close the input file
+         DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1);
+         analysisCanvas->GetCanvas()->Modified();
+         analysisCanvas->GetCanvas()->Update();
+         // Get the spectrum
+         histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
+         npeaks = 15;
+         double par[300];
+         spec = new TSpectrum(npeaks);
+         // Find spectrum background
+         histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
+         // Clone histogram and subtract background from it if we select that option
+         h2 = (TH1F*)histtemp->Clone("h2");
+         if(fitChecks->widgetChBox[0]->IsDown())
+            h2->Add(histback, -1);
+         // Search for the peaks
+         int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
+         printf("PhotonMu(): Found %d candidates to fit.\n",found);
+         npeaks = found;
+         // Set initial peak parameters
+         xpeaks = spec->GetPositionX();
+         for(j = 0; j < found; j++)
+         {
+            float xp = xpeaks[j];
+            int bin = h2->GetXaxis()->FindBin(xp);
+            float yp = h2->GetBinContent(bin);
+            par[3*j] = yp;
+            par[3*j+1] = xp;
+            par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
+         }
+         // Fit the histogram
+         fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
+         TVirtualFitter::Fitter(histtemp, 3*npeaks);
+         fit->SetParameters(par);
+         fit->SetNpx(300);
+         h2->Fit("fit","Q");
+         // Get the fitted parameters
+         fittingfunc = h2->GetFunction("fit");
+         fparam = fittingfunc->GetParameters();
+         fparamerr = fittingfunc->GetParErrors();
+         // Gather the parameters (mean peak value for now)
+         int j = 1;
+         int nrfit = 0;
+         while(1)
+         {
+            if( (fparam[j] < 1.E-30) || (nrfit > 8) )
+               break;
+            else
+            {
+               // Check if pedestal is above the lower limit and sigma is smaller than the mean
+               if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+               {
+                  // 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
+                  meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
+                  sigmasel[nrfit] = fparam[j+1];
+                  nrfit++;
+               }
+            }
+            j+=3;
+         }
+         TMath::Sort(nrfit, meansel, sortindex, kFALSE);
+         meanparam = meansel[sortindex[0]];
+         paramsigma = sigmasel[sortindex[0]];
+         for(j = 0; j < nrfit; j++)
+            if(DBGSIG)
+               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+         j = 0;
+         adcpedestal[0] = 0;
+         adcpedestal[1] = -1;
+         while(1)
+         {
+            int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
+            int yp = histtemp->GetBinContent(bin);
+            // 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)
+            if(adcpedestal[1] == -1)
+            {
+               adcpedestal[0] = j+meanparam+paramsigma;
+               adcpedestal[1] = yp;
+            }
+            else
+            {
+               if( (npeaks > 1) && (adcpedestal[1] >= yp) )
+               {
+                  adcpedestal[0] = j+meanparam+paramsigma;
+                  adcpedestal[1] = yp;
+               }
+               else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
+               {
+                  adcpedestal[0] = j+meanparam+paramsigma;
+                  adcpedestal[1] = yp;
+               }
+               else
+                  break;
+            }
+            j++;
+            if(j > 50) break;
+         }
+         if(midPeak->widgetChBox[0]->IsDown())
+         {
+            if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
+               m = TMath::Floor(meanparam);
+            else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
+               m = TMath::Ceil(meanparam);
+            int bin = histtemp->GetXaxis()->FindBin(m);
+            adcpedestal[0] = m;
+            printf("midpeak x = %d, ", adcpedestal[0]);
+            adcpedestal[1] = histtemp->GetBinContent(bin);
+         }
+         // Option to show the fit
+         fittingfunc->Draw("L SAME");
+         analysisCanvas->GetCanvas()->Modified();
+         analysisCanvas->GetCanvas()->Update();
+         printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
+         // Delete the opened histogram and spectrum
+         delete spec;
+         delete inroot;
+//       return;
+      }
+      printf("PhotonMu(): Continuing with the rest of the spectra.\n");
+      // Check all histograms for pedestal peak values
       for(int i = 0; i < (int)nrfiles; i++)
+      {
          if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
          {
             printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n");
             // Replot the spectrum on analysisCanvas and do not close the input file
             DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
             // Get the spectrum
             histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
             npeaks = 15;
             double par[300];
             spec = new TSpectrum(npeaks);
             // Find spectrum background
             histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
             // Clone histogram and subtract background from it if we select that option
             h2 = (TH1F*)histtemp->Clone("h2");
             if(fitChecks->widgetChBox[0]->IsDown())
                h2->Add(histback, -1);
             // Search for the peaks
             int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
             printf("PhotonMu(): Found %d candidates to fit.\n",found);
             npeaks = found;
             // Set initial peak parameters
             xpeaks = spec->GetPositionX();
             for(j = 0; j < found; j++)
+            {
                float xp = xpeaks[j];
                int bin = h2->GetXaxis()->FindBin(xp);
                float yp = h2->GetBinContent(bin);
                par[3*j] = yp;
                par[3*j+1] = xp;
                par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
+            }
             // Fit the histogram
             fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
             TVirtualFitter::Fitter(histtemp, 3*npeaks);
             fit->SetParameters(par);
             fit->SetNpx(300);
             h2->Fit("fit","Q");
             // Get the fitted parameters
             fittingfunc = h2->GetFunction("fit");
             fparam = fittingfunc->GetParameters();
             fparamerr = fittingfunc->GetParErrors();
             // Gather the parameters (mean peak value for now)
             int j = 1;
             int nrfit = 0;
             while(1)
+            {
                if( (fparam[j] < 1.E-30) || (nrfit > 8) )
                   break;
                else
+               {
                   // Check if pedestal is above the lower limit and sigma is smaller than the mean
                   if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+                  {
                      // 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
                      meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
                      sigmasel[nrfit] = fparam[j+1];
                      nrfit++;
+                  }
+               }
                j+=3;
+            }
             TMath::Sort(nrfit, meansel, sortindex, kFALSE);
             fittingfunc->Draw("SAME");
-            analysisCanvas->GetCanvas()->Modified();
-            analysisCanvas->GetCanvas()->Update();
-            meanparam = meansel[sortindex[0]];
-            paramsigma = sigmasel[sortindex[0]];
-            for(j = 0; j < nrfit; j++)
-               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
-            return;
-         }
-         if(files->At(i))
-         {
-            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
-            {
-               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
-               darkhist = i;
-            }
-            // Replot the spectrum on analysisCanvas and do not close the input file
-            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
-            // Get the spectrum
-            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
-            npeaks = 15;
-            double par[300];
-            spec = new TSpectrum(npeaks);
-            // Find spectrum background
-            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
-            // Clone histogram and subtract background from it if we select that option
-            h2 = (TH1F*)histtemp->Clone("h2");
-            if(fitChecks->widgetChBox[0]->IsDown())
-               h2->Add(histback, -1);
-            // Search for the peaks
-            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
-            printf("PhotonMu(): Found %d candidates to fit.\n",found);
-            npeaks = found;
-            // Set initial peak parameters
-            xpeaks = spec->GetPositionX();
-            for(j = 0; j < found; j++)
-            {
-               float xp = xpeaks[j];
-               int bin = h2->GetXaxis()->FindBin(xp);
-               float yp = h2->GetBinContent(bin);
-               par[3*j] = yp;
-               par[3*j+1] = xp;
-               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
-            }
-            // Fit the histogram
-            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
-            TVirtualFitter::Fitter(histtemp, 3*npeaks);
-            fit->SetParameters(par);
-            fit->SetNpx(300);
-            h2->Fit("fit","Q");
-            // Get the fitted parameters
-            fittingfunc = h2->GetFunction("fit");
-            fparam = fittingfunc->GetParameters();
-            fparamerr = fittingfunc->GetParErrors();
-            // Gather the parameters (mean peak value for now)
-            int j = 1;
-            int nrfit = 0;
-            while(1)
-            {
-               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
-                  break;
-               else
-               {
-                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
-                  if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
-                  {
-                     // 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
-                     meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
-                     sigmasel[nrfit] = fparam[j+1];
-                     nrfit++;
-                  }
-               }
-               j+=3;
-            }
-            TMath::Sort(nrfit, meansel, sortindex, kFALSE);
             meanparam = meansel[sortindex[0]];
             paramsigma = sigmasel[sortindex[0]];
             for(j = 0; j < nrfit; j++)
-               if(DBGSIG)
-                  printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+            return;
+         }
-            j = 0;
+         if(files->At(i))
+         {
-            adcpedestal[0] = 0;
+            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
+            {
+               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
-            adcpedestal[1] = -1;
+               darkhist = i;
+            }
-            while(1)
+            if(nopeaks == -1)
             {
+               // Replot the spectrum on analysisCanvas and do not close the input file
+               DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
-               int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
+               analysisCanvas->GetCanvas()->Modified();
-               int yp = histtemp->GetBinContent(bin);
+               analysisCanvas->GetCanvas()->Update();
+               // Get the spectrum
-               // 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)
+               histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
-               if(adcpedestal[1] == -1)
+               npeaks = 15;
-               {
+               double par[300];
-                  adcpedestal[0] = j+meanparam+paramsigma;
+               spec = new TSpectrum(npeaks);
+               // Find spectrum background
+               histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
+               // Clone histogram and subtract background from it if we select that option
+               h2 = (TH1F*)histtemp->Clone("h2");
+               if(fitChecks->widgetChBox[0]->IsDown())
-                  adcpedestal[1] = yp;
+                  h2->Add(histback, -1);
-               }
+               // Search for the peaks
+               int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
+               printf("PhotonMu(): Found %d candidates to fit.\n",found);
-               else
+               npeaks = found;
+               // Set initial peak parameters
+               xpeaks = spec->GetPositionX();
+               for(j = 0; j < found; j++)
+               {
+                  float xp = xpeaks[j];
+                  int bin = h2->GetXaxis()->FindBin(xp);
+                  float yp = h2->GetBinContent(bin);
+                  par[3*j] = yp;
+                  par[3*j+1] = xp;
+                  par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
+               }
+               // Fit the histogram
+               fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
+               TVirtualFitter::Fitter(histtemp, 3*npeaks);
+               fit->SetParameters(par);
+               fit->SetNpx(300);
+               h2->Fit("fit","Q");
+               // Get the fitted parameters
+               fittingfunc = h2->GetFunction("fit");
+               fparam = fittingfunc->GetParameters();
+               fparamerr = fittingfunc->GetParErrors();
+               // Gather the parameters (mean peak value for now)
+               int j = 1;
+               int nrfit = 0;
+               while(1)
+               {
-                  if( (npeaks > 1) && (adcpedestal[1] >= yp) )
+                  if( (fparam[j] < 1.E-30) || (nrfit > 8) )
+                     break;
+                  else
+                  {
+                     // Check if pedestal is above the lower limit and sigma is smaller than the mean
+                     if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+                     {
+                        // 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
+                        meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
-                     adcpedestal[0] = j+meanparam+paramsigma;
+                        sigmasel[nrfit] = fparam[j+1];
-                     adcpedestal[1] = yp;
+                        nrfit++;
+                     }
+                  }
+                  j+=3;
+               }
+               TMath::Sort(nrfit, meansel, sortindex, kFALSE);
+               meanparam = meansel[sortindex[0]];
+               paramsigma = sigmasel[sortindex[0]];
+               for(j = 0; j < nrfit; j++)
+                  if(DBGSIG)
+                     printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+               j = 0;
+               adcpedestal[0] = 0;
+               adcpedestal[1] = -1;
+               while(1)
+               {
+                  int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
+                  int yp = histtemp->GetBinContent(bin);
-                  else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
+                  // 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)
+                  if(adcpedestal[1] == -1)
+                  {
                      adcpedestal[0] = j+meanparam+paramsigma;
                      adcpedestal[1] = yp;
+                  }
                   else
+                  {
+                     if( (npeaks > 1) && (adcpedestal[1] >= yp) )
+                     {
+                        adcpedestal[0] = j+meanparam+paramsigma;
+                        adcpedestal[1] = yp;
+                     }
+                     else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
+                     {
+                        adcpedestal[0] = j+meanparam+paramsigma;
+                        adcpedestal[1] = yp;
+                     }
+                     else
-                     break;
+                        break;
+                  }
+                  j++;
+                  if(j > 50) break;
+               }
-               j++;
-               if(j > 50) break;
-            }
-            if( (npeaks > 1) && (nrfit > 1) )
-            {
-               int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
-               adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
-               printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
-               adcpedestal[1] = histtemp->GetBinContent(bin);
-            }
-            if(midPeak->widgetChBox[0]->IsDown())
-            {
-               if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
-                  m = TMath::Floor(meanparam);
-               else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
-                  m = TMath::Ceil(meanparam);
-               int bin = histtemp->GetXaxis()->FindBin(m);
-               adcpedestal[0] = m;
-               printf("midpeak x = %d, ", adcpedestal[0]);
-               adcpedestal[1] = histtemp->GetBinContent(bin);
-            }
+               if(npeaks > 1)
+               {
+                  int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
+                  adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
+                  printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
+                  adcpedestal[1] = histtemp->GetBinContent(bin);
+               }
+               if(midPeak->widgetChBox[0]->IsDown())
+               {
+                  if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
+                     m = TMath::Floor(meanparam);
+                  else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
+                     m = TMath::Ceil(meanparam);
+                  int bin = histtemp->GetXaxis()->FindBin(m);
+                  adcpedestal[0] = m;
+                  printf("midpeak x = %d, ", adcpedestal[0]);
+                  adcpedestal[1] = histtemp->GetBinContent(bin);
+               }
-/*          // Option to show the fit
+/*             // Option to show the fit
-            fittingfunc->Draw("L SAME");
+               fittingfunc->Draw("L SAME");
-            analysisCanvas->GetCanvas()->Modified();
+               analysisCanvas->GetCanvas()->Modified();
-            analysisCanvas->GetCanvas()->Update();*/
+               analysisCanvas->GetCanvas()->Update();*/
-            printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
+               printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
-            // Delete the opened histogram and spectrum
+               // Delete the opened histogram and spectrum
-            delete spec;
+               delete spec;
-            delete inroot;
+               delete inroot;
+            }
             // Open the input file and read header, ADC and TDC values
             sprintf(ctemp, "%s", files->At(i)->GetTitle());
             inroot = new TFile(ctemp, "READ");
             inroot->GetObject("header_data", header_data);
             inroot->GetObject("meas_data", meas_data);
             // Reading the header
             if( header_data->FindBranch("angle") )
+            {
                header_data->SetBranchAddress("angle", &evtheader.angle);
                header_data->GetEntry(0);
+            }
             else
+            {
                printf("PhotonMu(): Error! Selected file has no angle header value. Please edit header to add the angle header value.\n");
                break;
+            }
             char rdc[256];
             j = selectCh->widgetNE[0]->GetNumber();
             double rangetdc[2];
             rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
             rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
             k = 0;
             k2 = 0;
             m = 0;
             m2 = 0;
             // Reading the data
             for(int e = 0; e < meas_data->GetEntries(); e++)
+            {
                sprintf(rdc, "ADC%d", j);
                meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
                meas_data->GetEntry(e);
                sprintf(rdc, "TDC%d", j);
                meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
                meas_data->GetEntry(e);
                // If our data point is inside the TDC window
                if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
+               {
                   // Gather only the integral of the pedestal
                   if((double)evtdata.adcdata[j] < (double)adcpedestal[0]+0.5 )
+                  {
                      k2++;
                      m2 += evtdata.adcdata[j];
+                  }
                   // Gather the complete integral
                   k++;
                   m += evtdata.adcdata[j];
+               }
+            }
             // Determine the angle, mu and relative PDE
             angle[i] = (double)(evtheader.angle);
             integralCount[i] += (double)m;
             printf("PhotonMu(): %lf: Complete integral (%d evts) = %lf\n", angle[i], k, integralCount[i]);
             integralPedestal[i] += (double)m2;
             printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]);
             if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) )
                zeromu = i;
-            // Checking for errors when fitting a histogram
-            if(k2 == 0)
-            {
-               printf("PhotonMu(): No pedestal entries found. Check the fitting results.\n");
-               muval[i] = -1;
-            }
-            else
-               muval[i] = -TMath::Log((double)k2/(double)k);
+            muval[i] = -TMath::Log((double)k2/(double)k);
             printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]);
             inroot->Close();
             delete inroot;
+         }
          // Update the progress bar
          progVal = (float)(90.00/nrfiles)*i;
          analysisProgress->widgetPB->SetPosition(progVal);
          gVirtualX->Update(1);
+      }
       printf("PhotonMu(): %d files were selected.\n", nrfiles);
       printf("PhotonMu(): angle\tmu\trelative PDE\n");
       m = 0;
       // Set the 0 degree muval, reuse meansel[1]
       meansel[1] = muval[zeromu];
       printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]);
       // TODO - point estimation still not working correctly!
       for(int i = 0; i < (int)nrfiles; i++)
+      {
          // Estimate next point and check error (5 point least square fit estimation)
          if( ((i > 4) && (m == i)) || ((i > 5) && (m < i)) )
+         {
             // Set exclude signal to false
             exclude = false;
             // Get next point values (if zero value -> need to add the dark hist value again)
             pointest[10] = angle[i];
             pointest[11] = muval[i];
             // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0]
             meansel[0] = PointEstimate(5, pointest);    // PointEstimate only works with very small step size
             if(meansel[0] > accError->widgetNE[0]->GetNumber())
+            {
                printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]);
                exclude = true;
+            }
             // Value with 0 angle and dark histogram are always needed, so should not be excluded
             if(i == darkhist)
                exclude = false;
-            // Wrong fit
-            if(muval[i] == -1)
-               exclude = true;
             // If nothing excluded, pass the points in pointest variable like in a FIFO
             if(!exclude)
+            {
                for(int j = 0; j < 10; j++)
+               {
                   if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]);
                   pointest[j] = pointest[j+2];
+               }
+            }
             else
+            {
                for(int j = 0; j < 10; j++)
                   if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]);
+            }
+         }
          else
+         {
             // First 5 points act as estimator points for next one
             pointest[2*m] = angle[i];
             pointest[2*m+1] = muval[i];
+         }
          // Run only if we have a dark run histogram and middle pedestal peak estimation
          if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() )
+         {
             if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]);
             // Subtract the dark value from all values
             angle[m] = angle[i];
             muval[m] = muval[i] - muval[darkhist];
             if(DBGSIG) printf("angle = %lf, newmuval = %lf, darkmuval = %lf, ", angle[m], muval[m], muval[darkhist]);
             // Calculate relative PDE
 //          pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
             pdeval[m] = muval[m]/((meansel[1]-muval[darkhist])*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
             if(DBGSIG) printf("pdeval = %lf\n", pdeval[m]);
             // Only increase counter if error is not too big
             if( (darkhist != i) && (!exclude) )
                m++;
+         }
          else
+         {
             // Relative PDE calculation
             angle[m] = angle[i];
             muval[m] = muval[i];
 //            pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
             pdeval[m] = muval[m]/(meansel[1]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
             // Only increase counter if error is not too big
             if(!exclude)
                m++;
+         }
          printf("PhotonMu(): %lf\t%lf\t%lf\n", angle[i], muval[i], pdeval[i]);
+      }
       // Check for range of values to plot
       double plotMax = 0.;
       for(int i = 0; i < (int)nrfiles; i++)
+      {
          plotMax = TMath::Max(plotMax, muval[i]);
          plotMax = TMath::Max(plotMax, pdeval[i]);
+      }
       if(plotMax <= 0.)
          plotMax = 1.1;
       printf("PhotonMu(): Maximum value: %lf\n", plotMax);
       if(DBGSIG) printf("\n");
       if(darkhist != -1)
          printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-1-m));
       else
          printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-m));
       // Plot mu and PDE angle dependance plots
       if(edit == 0)
          gCanvas = new TCanvas("canv","canv",1200,900);
       else if(edit == 1)
          gCanvas = tempAnalysisCanvas->GetCanvas();
       gCanvas->cd();
       gCanvas->SetGrid();
       TGraph *pde = new TGraph(m, angle, pdeval);
       pde->SetMarkerStyle(21);
       pde->SetMarkerSize(0.7);
       pde->SetMarkerColor(2);
       pde->SetLineWidth(1);
       pde->SetLineColor(2);
       pde->GetXaxis()->SetLabelSize(0.030);
       pde->GetXaxis()->CenterTitle();
 //      pde->GetXaxis()->SetRange(angle[0],angle[nrfiles-1]);
 //      pde->GetXaxis()->SetRangeUser(angle[0],angle[nrfiles-1]);
       pde->GetXaxis()->SetRange(-90.0,90.0);
       pde->GetXaxis()->SetRangeUser(-90.0,90.0);
       pde->GetXaxis()->SetLimits(-90.0,90.0);
       pde->GetYaxis()->SetTitleOffset(1.2);
       pde->GetYaxis()->SetLabelSize(0.030);
       pde->GetYaxis()->CenterTitle();
       pde->GetYaxis()->SetRange(0., 1.1*plotMax);
       pde->GetYaxis()->SetRangeUser(0., 1.1*plotMax);
       pde->GetYaxis()->SetLimits(0., 1.1*plotMax);
       pde->SetName("pde");
       pde->Draw("ALP");
       pde->SetTitle(";Incidence angle (#circ);Relative PDE(#theta) / #mu(#theta)");
       TGraph *mugr = new TGraph(m, angle, muval);
       mugr->SetMarkerStyle(20);
       mugr->SetMarkerSize(0.7);
       mugr->SetMarkerColor(4);
       mugr->SetLineWidth(1);
       mugr->SetLineColor(4);
       mugr->SetName("muval");
       mugr->Draw("SAME;LP");
       gCanvas->Modified();
       gCanvas->Update();
       if(edit == 0)
+      {
          remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
          sprintf(exportname, "%s_pde.pdf", ctemp);
          gCanvas->SaveAs(exportname);
          delete mugr;
          delete pde;
          delete gCanvas;
+      }
       else if(edit == 1)
+      {
          gCanvas->Modified();
          gCanvas->Update();
+      }
       // Update the progress bar
       analysisProgress->widgetPB->SetPosition(100.);
       gVirtualX->Update(1);
+   }
+}
 void TGAppMainFrame::BreakdownVolt(TList *files, int edit)
+{
    unsigned int nrfiles = files->GetSize();
    char ctemp[1024];
    char exportname[1024];
    char paramname[1024];
    int j, k = 0;
    TCanvas *gCanvas;
    TTree *header_data, *meas_data;
    // Fitting variables
    TSpectrum *spec;
    TH1F *histtemp;
    TH1 *histback;
    TH1F *h2;
    float *xpeaks;
    TF1 *fit;
    TF1 *fittingfunc;
    TLatex *latex;
    double *fparam, *fparamerr;
    double meansel[20];
    double meanselerr[20];
    double sigmasel[20];
    double meanparam, meanparamerr, paramsigma;
    int sortindex[20];
    bool exclude = false;
    int p = 0;
    double volt[nrfiles], volterr[nrfiles], sep[3][nrfiles], seperr[3][nrfiles];
    int first = 1;
    FILE *fp;
    remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
    sprintf(paramname, "%s_fitresult.txt", ctemp);
    fp = fopen(paramname, "w");
    fclose(fp);
    int peaklimit = minPeak->widgetNE[0]->GetNumber()+1; // +1 to account for the pedestal peak
    // Zero the parameter values
    for(int i = 0; i < nrfiles; i++)
+   {
       volt[i] = 0; volterr[i] = 0;
       sep[0][i] = 0; sep[1][i] = 0; sep[2][i] = 0;
       seperr[0][i] = 0; seperr[1][i] = 0; seperr[2][i] = 0;
       if(i < 20) { meansel[i] = 0; meanselerr[i] = 0; sigmasel[i] = 0; }
+   }
    float progVal = 0;
    analysisProgress->widgetPB->SetPosition(progVal);
    gVirtualX->Update(1);
    // Start if we select at least one file
    if(nrfiles > 0)
+   {
       for(int i = 0; i < (int)nrfiles; i++)
+      {
          if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
+         {
             printf("BreakdownVolt(): Only one file selected. Not running analysis, just showing the fit.\n");
             // Replot the spectrum on analysisCanvas and do not close the input file
             DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
             // Get the spectrum
             histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
             npeaks = 15;
             double par[300];
             spec = new TSpectrum(npeaks);
             // Find spectrum background
             histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
             // Clone histogram and subtract background from it if we select that option
             h2 = (TH1F*)histtemp->Clone("h2");
             if(fitChecks->widgetChBox[0]->IsDown())
                h2->Add(histback, -1);
             // Search for the peaks
             int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
             printf("PhotonMu(): Found %d candidates to fit.\n",found);
             npeaks = found;
             // Set initial peak parameters
             xpeaks = spec->GetPositionX();
             for(j = 0; j < found; j++)
+            {
                float xp = xpeaks[j];
                int bin = h2->GetXaxis()->FindBin(xp);
                float yp = h2->GetBinContent(bin);
                par[3*j] = yp;
                par[3*j+1] = xp;
                par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
+            }
             // Fit the histogram
             fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
             TVirtualFitter::Fitter(histtemp, 3*npeaks);
             fit->SetParameters(par);
             fit->SetNpx(300);
             h2->Fit("fit","Q");
             // Get the fitted parameters
             fittingfunc = h2->GetFunction("fit");
             fparam = fittingfunc->GetParameters();
             fparamerr = fittingfunc->GetParErrors();
             // Gather the parameters (mean peak value for now)
             int j = 1;
             int nrfit = 0;
             while(1)
+            {
                if( (fparam[j] < 1.E-30) || (nrfit > 8) )
                   break;
                else
+               {
                   // Check if pedestal is above the lower limit and sigma is smaller than the mean
                   if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+                  {
                      // 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
                      meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
+                     meanselerr[nrfit] = fparamerr[j];
                      sigmasel[nrfit] = fparam[j+1];
                      nrfit++;
+                  }
+               }
                j+=3;
+            }
             TMath::Sort(nrfit, meansel, sortindex, kFALSE);
             fittingfunc->Draw("SAME");
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
             meanparam = meansel[sortindex[0]];
             meanparamerr = meanselerr[sortindex[0]];
             paramsigma = sigmasel[sortindex[0]];
             for(j = 0; j < nrfit; j++)
                printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]);
             return;
+         }
          if(files->At(i))
+         {
             // Replot the spectrum on analysisCanvas and do not close the input file
             DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
             // Get the spectrum
             histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
             npeaks = 15;
             double par[300];
             spec = new TSpectrum(npeaks);
             // Find spectrum background
             histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
             // Clone histogram and subtract background from it if we select that option
             h2 = (TH1F*)histtemp->Clone("h2");
             if(fitChecks->widgetChBox[0]->IsDown())
                h2->Add(histback, -1);
             // Search for the peaks
             int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
             printf("BreakdownVolt(): Found %d candidates to fit.\n",found);
             npeaks = found;
             // Set initial peak parameters
             xpeaks = spec->GetPositionX();
             for(j = 0; j < found; j++)
+            {
                float xp = xpeaks[j];
                int bin = h2->GetXaxis()->FindBin(xp);
                float yp = h2->GetBinContent(bin);
                par[3*j] = yp;
                par[3*j+1] = xp;
                par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
+            }
             // Fit the histogram
             fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
             TVirtualFitter::Fitter(histtemp, 3*npeaks);
             fit->SetParameters(par);
             fit->SetNpx(300);
             h2->Fit("fit","Q");
             // Get the fitted parameters
             fittingfunc = h2->GetFunction("fit");
             if(nrfiles == 1)    // TODO: Show the fit if only one file is selected
                fittingfunc->Draw("SAME");
             fparam = fittingfunc->GetParameters();
             fparamerr = fittingfunc->GetParErrors();
             // Gather the parameters (mean peak value for now)
             int j = 1;
             int nrfit = 0;
             while(1)
+            {
                if( (fparam[j] < 1.E-30) || (fparamerr[j] < 1.E-10) )// TODO: Maybe not correct for the error
                   break;
                else
+               {
                   // Check if pedestal is above the lower limit and sigma is smaller than the mean
                   if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+                  {
                      meansel[nrfit] = fparam[j];
                      meanselerr[nrfit] = fparamerr[j];
                      sigmasel[nrfit] = fparam[j+1];
                      nrfit++;
+                  }
+               }
                j+=3;
+            }
             TMath::Sort(nrfit, meansel, sortindex, kFALSE);
             meanparam = meansel[sortindex[0]];
             meanparamerr = meanselerr[sortindex[0]];
             paramsigma = sigmasel[sortindex[0]];
             for(j = 0; j < nrfit; j++)
+            {
                if(DBGSIG)
                   printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]);
+            }
             // Delete the opened histogram and spectrum
             delete spec;
             delete inroot;
             // Open the input file and read header, ADC and TDC values
             sprintf(ctemp, "%s", files->At(i)->GetTitle());
             inroot = new TFile(ctemp, "READ");
             inroot->GetObject("header_data", header_data);
             inroot->GetObject("meas_data", meas_data);
             // Reading the header
             if( header_data->FindBranch("biasvolt") )
+            {
                header_data->SetBranchAddress("biasvolt", &evtheader.biasvolt);
                header_data->GetEntry(0);
+            }
             else
+            {
                printf("BreakdownVolt(): Error! Selected file has no bias voltage header value. Please edit header to add the bias voltage header value.\n");
                break;
+            }
 //            h2->SetStats(0);
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
             // Save each fitting plot
             if(fitChecks->widgetChBox[1]->IsDown())     // TODO: Check if this works
+            {
                remove_ext((char*)files->At(i)->GetTitle(), ctemp);
                sprintf(exportname, "%s_fit.pdf", ctemp);
                analysisCanvas->GetCanvas()->SaveAs(exportname);
+            }
             // Calculate the separation between peaks
             volt[p] = evtheader.biasvolt;
             volterr[p] = 1.e-5;
             if(nrfit == 3)
+            {
                sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
                seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
                exclude = (seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber());
+            }
             else if(nrfit == 4)
+            {
                sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
                sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]];
                seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
                seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]);
                exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber()));
+            }
             else if(nrfit > 4)
+            {
                sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
                sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]];
                sep[2][p] = meansel[sortindex[4]] - meansel[sortindex[3]];
                seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
                seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]);
                seperr[2][p] = TMath::Abs(meanselerr[sortindex[4]]) + TMath::Abs(meanselerr[sortindex[3]]);
                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()));
+            }
             else
+            {
                printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n");
                exclude = false;
+            }
+            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]);
             // Write out parameters to a file
             fp = fopen(paramname, "a");
             if(exclude)
+            {
                if(first == 1)
+               {
                   fprintf(fp, "%le\t%d\t", evtheader.biasvolt, nrfit);
                   for(j = 0; j < nrfit; j++)
                      fprintf(fp, "%le\t%le\t", meansel[sortindex[j]], meanselerr[sortindex[j]]);
                   fprintf(fp,"\n");
                   first = 0;
+               }
                p++;
+            }
             else
+            {
                if(nrfit == 3)
                   printf("Point (at %.2lfV) rejected due to too large errors: %lf\n", volt[p], seperr[0][p]/sep[0][p]);
                else if(nrfit == 4)
                   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]);
                else if(nrfit > 4)
                   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]);
+            }
             fclose(fp);
+         }
          if(nrfiles == 1) break;
          first = 1;
          // Update the progress bar
          progVal = (float)(90.00/nrfiles)*i;
          analysisProgress->widgetPB->SetPosition(progVal);
          gVirtualX->Update(1);
+      }
       printf("BreakdownVolt(): %d files were selected.\n", nrfiles);
       printf("BreakdownVolt(): Number of points to plot is %d\n", p);
       // Plot and fit breakdown voltage plot
       if(edit == 0)
          gCanvas = new TCanvas("canv","canv",1200,900);
       else if(edit == 1)
          gCanvas = tempAnalysisCanvas->GetCanvas();
       gCanvas->cd();
       gCanvas->SetGrid();
       TGraphErrors* bdplot;
       k = peakSepCalc->widgetNE[0]->GetNumber();
       if(k < 4)
          bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]);
       else
+      {
-         printf("BreakdownVold(): Unsupported peak separation selected (%d).\n", k);
+         printf("BreakdownVolt(): Unsupported peak separation selected (%d).\n", k);
          return;
+      }
       bdplot->SetMarkerStyle(20);
       bdplot->SetMarkerSize(0.4);
       bdplot->SetMarkerColor(kBlue);
       bdplot->SetLineWidth(1);
       bdplot->SetLineColor(kBlue);
       bdplot->GetXaxis()->SetLabelSize(0.030);
       bdplot->GetXaxis()->CenterTitle();
 //      bdplot->GetXaxis()->SetRange(-90,90);
 //      bdplot->GetXaxis()->SetRangeUser(-90,90);
 //      bdplot->GetXaxis()->SetLimits(-90,90);
       bdplot->GetYaxis()->SetTitleOffset(1.2);
       bdplot->GetYaxis()->SetLabelSize(0.030);
       bdplot->GetYaxis()->CenterTitle();
 //      bdplot->GetYaxis()->SetRangeUser(0., 1.2);
       bdplot->SetName("bdplot");
       bdplot->Draw("AP");
       bdplot->SetTitle(";Bias voltage (V);Peak separation");
       // Fit the breakdown voltage plot with a line
       bdplot->Fit("pol1","Q");
       TF1 *fit1 = bdplot->GetFunction("pol1");
       meansel[0] = fit1->GetParameter(0);
       meanselerr[0] = fit1->GetParError(0);
       meansel[1] = fit1->GetParameter(1);
       meanselerr[1] = fit1->GetParError(1);
       meansel[2] = -meansel[0]/meansel[1];
       if(!cleanPlots)
+      {
          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])) );
          latex = new TLatex();
          latex->SetTextSize(0.039);
-         latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[p-1]], ctemp);
+         latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[1]], ctemp);
          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])) );
+      }
       else
          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])) );
       if(edit == 0)
+      {
          remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
          sprintf(exportname, "%s_breakdown.pdf", ctemp);
          gCanvas->SaveAs(exportname);
          delete bdplot;
          delete gCanvas;
+      }
       else if(edit == 1)
+      {
          gCanvas->Modified();
          gCanvas->Update();
+      }
       // Update the progress bar
       analysisProgress->widgetPB->SetPosition(100.);
       gVirtualX->Update(1);
+   }
+}
 void TGAppMainFrame::SurfaceScan(TList *files, int edit)
+{
    unsigned int nrfiles = files->GetSize();
    char ctemp[1024];
    char exportname[1024];
    int j, k = 0, m = 0, n = 0;
    TTree *header_data, *meas_data;
    double *integralCount;
    integralCount = new double[nrfiles];
    double *surfx, *surfy;
    surfx = new double[nrfiles];
    surfy = new double[nrfiles];
    double xsurfmin = 0, ysurfmin = 0;
    int nrentries;
    double minInteg, maxInteg;
    bool norm = surfScanOpt->widgetChBox[0]->IsDown();
    double curyval;
 //   bool edge2d = false;
    TCanvas *gCanvas;
    float progVal = 0;
    analysisProgress->widgetPB->SetPosition(progVal);
    gVirtualX->Update(1);
    // Zero the integral count and accumulated vaules
    for(int i = 0; i < (int)nrfiles; i++) integralCount[i] = 0;
    // Start if we select at more than one file
    if( (nrfiles > 1) && (multiSelect->widgetChBox[0]->IsOn()) )
+   {
       printf("SurfaceScan(): Creating a surface plot. Please wait...\n");
       for(int i = 0; i < (int)nrfiles; i++)
+      {
          if(files->At(i))
+         {
             n++;
             // Read the X and Y positions from header and ADC and TDC values from the measurements
             sprintf(ctemp, "%s", files->At(i)->GetTitle());
             inroot = new TFile(ctemp, "READ");
             inroot->GetObject("header_data", header_data);
             inroot->GetObject("meas_data", meas_data);
             header_data->SetBranchAddress("xpos", &evtheader.xpos);
             header_data->GetEntry(0);
             header_data->SetBranchAddress("ypos", &evtheader.ypos);
             header_data->GetEntry(0);
             char rdc[256];
             j = selectCh->widgetNE[0]->GetNumber();
             double rangetdc[2];
             rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
             rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
             k = 0;
             m = 0;
             // Reading the data
             for(int e = 0; e < meas_data->GetEntries(); e++)
+            {
                sprintf(rdc, "ADC%d", j);
                meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
                meas_data->GetEntry(e);
                sprintf(rdc, "TDC%d", j);
                meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
                meas_data->GetEntry(e);
                // Use data point only if it is inside the TDC window
                if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
+               {
                   k++;
                   m += evtdata.adcdata[j];
+               }
+            }
             // X, Y and Z values from each file (table units or microns)
             if(posUnitsPlot->widgetCB->GetSelected() == 0)
+            {
                if(n == 1)
+               {
                   xsurfmin = (double)(evtheader.xpos);
                   ysurfmin = (double)(evtheader.ypos);
+               }
                surfx[i] = (double)(evtheader.xpos);
                surfy[i] = (double)(evtheader.ypos);
+            }
             else if(posUnitsPlot->widgetCB->GetSelected() == 1)
+            {
                if(n == 1)
+               {
                   xsurfmin = (double)(evtheader.xpos*lenconversion);
                   ysurfmin = (double)(evtheader.ypos*lenconversion);
+               }
                surfx[i] = (double)(evtheader.xpos*lenconversion);
                surfy[i] = (double)(evtheader.ypos*lenconversion);
+            }
             // Save result for later plotting
             if(norm)
                integralCount[i] += ((double)m)/((double)k);
             else
                integralCount[i] += m;
             inroot->Close();
             delete inroot;
+         }
          // Update the progress bar
          progVal = (float)(75.00/nrfiles)*i;
          analysisProgress->widgetPB->SetPosition(progVal);
          gVirtualX->Update(1);
+      }
       nrentries = n;
       printf("SurfaceScan(): %d files were selected.\n", nrfiles);
       // Check for Minimum and Maximum values and normalize to 1, if normalization is selected
       if(norm)
+      {
          minInteg = TMath::MinElement(nrfiles, integralCount);
          for(int i = 0; i < nrfiles; i++)
+         {
             integralCount[i] -= minInteg;
             if(DBGSIG) printf("Subtraction: %lf\n", integralCount[i]);
+         }
          maxInteg = TMath::MaxElement(nrfiles, integralCount);
          for(int i = 0; i < nrfiles; i++)
+         {
             integralCount[i] = integralCount[i]/maxInteg;
             if(DBGSIG) printf("Normalization: %lf\n", integralCount[i]);
+         }
+      }
       // Make the 2D surface plot
       if(edit == 0)
          gCanvas = new TCanvas("canv","canv",1100,900);
       else
          gCanvas = tempAnalysisCanvas->GetCanvas();
       double range[4];
       TGraph2D *gScan2D;
       gScan2D = new TGraph2D();
-      gScan2D->SetName("surfscan");
       range[0] = TMath::MinElement(nrfiles, surfx);
       range[1] = TMath::MaxElement(nrfiles, surfx);
       range[2] = TMath::MinElement(nrfiles, surfy);
       range[3] = TMath::MaxElement(nrfiles, surfy);
       for(int i = 0; i < nrfiles; i++)
+      {
          if(DBGSIG)
+         {
             if(surfScanOpt->widgetChBox[1]->IsDown())
                printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i]-range[0], surfy[i]-range[2], integralCount[i]);
             else
                printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i], surfy[i], integralCount[i]);
+         }
          if(surfScanOpt->widgetChBox[1]->IsDown())
             gScan2D->SetPoint(i, surfx[i]-range[0], surfy[i]-range[2], integralCount[i]);
          else
             gScan2D->SetPoint(i, surfx[i], surfy[i], integralCount[i]);
          // Update the progress bar
          progVal = (float)(9.00/nrfiles)*i+90.00;
          analysisProgress->widgetPB->SetPosition(progVal);
          gVirtualX->Update(1);
+      }
       if(surfScanOpt->widgetChBox[1]->IsDown())
+      {
          range[1] -= range[0];
          range[3] -= range[2];
          range[0] -= range[0];
          range[2] -= range[2];
+      }
       gCanvas->cd();
       gStyle->SetPalette(1);
       gCanvas->SetLeftMargin(0.15);
       gCanvas->SetRightMargin(0.126);
       gCanvas->SetTopMargin(0.077);
       gScan2D->Draw("COLZ");
       gCanvas->Modified();
       gCanvas->Update();
       if(!cleanPlots)
+      {
          if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan2D->SetTitle("Surface scan;X [table units];Y [table units]");
          else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan2D->SetTitle("Surface scan;X [#mum];Y [#mum]");
+      }
       else
+      {
          if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan2D->SetTitle(";X [table units];Y [table units]");
          else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan2D->SetTitle(";X [#mum];Y [#mum]");
+      }
 /*      TGaxis *xax = (TGaxis*)gScan2D->GetXaxis();
       xax->SetMaxDigits(4);
       TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
       yax->SetMaxDigits(4);*/
       gScan2D->SetNpx((int)resolSurf->widgetNE[0]->GetNumber());
       gScan2D->SetNpy((int)resolSurf->widgetNE[1]->GetNumber());
       gCanvas->Modified();
       gCanvas->Update();
       gScan2D->GetXaxis()->SetTitleOffset(1.3);
       gScan2D->GetXaxis()->CenterTitle(kTRUE);
       gScan2D->GetXaxis()->SetLabelSize(0.027);
       gScan2D->GetXaxis()->SetLabelOffset(0.02);
       gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
       gScan2D->GetXaxis()->SetNoExponent(kTRUE);
       gScan2D->GetYaxis()->SetTitleOffset(1.9);
       gScan2D->GetYaxis()->CenterTitle(kTRUE);
       gScan2D->GetYaxis()->SetLabelSize(0.027);
       gScan2D->GetYaxis()->SetLabelOffset(0.02);
       gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
       gScan2D->GetYaxis()->SetNoExponent(kTRUE);
       gCanvas->Modified();
       gCanvas->Update();
       remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
       sprintf(exportname, "%s_surfscan.pdf", ctemp);
       gCanvas->SaveAs(exportname);
       // Update the progress bar
       analysisProgress->widgetPB->SetPosition(100.0);
       gVirtualX->Update(1);
       if(edit == 0)
+      {
          delete gScan2D;
          delete gCanvas;
+      }
       else if(edit == 1)
+      {
          gCanvas->Modified();
          gCanvas->Update();
-         UpdateIntegrateSurface(-1);
+      }
+   }
+}

Subversion Repositories f9daq

(root)/lab/sipmscan/trunk/src/analysis.cpp – Rev 172 → 173