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(currentOpenDir);
   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[1]->SetState(kButtonDown);
      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;
 
   // 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();
         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;
 
   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);
 
   // 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("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]]);
       
            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(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
            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;
 
            // 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;
 
            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;
 
            // 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());
                     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;
            }
 
            // 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);
         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);
         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();
      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();
      }
   }
}
 

Rev 146	Rev 167
1	#include "../include/sipmscan.h"	1	#include "../include/sipmscan.h"
2	#include "../include/workstation.h"	2	#include "../include/workstation.h"
3		3
4	#include <stdio.h>	4	#include <stdio.h>
5	#include <stdlib.h>	5	#include <stdlib.h>
6		6
7	// Peak detection function	7	// Peak detection function
8	int npeaks = 20;	8	int npeaks = 20;
9	double FindPeaks(double x, double par)	9	double FindPeaks(double x, double par)
10	{	10	{
11	double result = 0;	11	double result = 0;
12	for(int i = 0; i < npeaks; i++)	12	for(int i = 0; i < npeaks; i++)
13	{	13	{
14	double norm = par[3*i];	14	double norm = par[3*i];
15	double mean = par[3*i+1];	15	double mean = par[3*i+1];
16	double sigma = par[3*i+2];	16	double sigma = par[3*i+2];
17	result += norm*TMath::Gaus(x[0], mean, sigma);	17	result += norm*TMath::Gaus(x[0], mean, sigma);
18	}	18	}
19	return result;	19	return result;
20	}	20	}
21		21
22	// File browser for selecting the dark run histogram	22	// File browser for selecting the dark run histogram
23	void TGAppMainFrame::SelectDarkHist()	23	void TGAppMainFrame::SelectDarkHist()
24	{	24	{
25	TGFileInfo file_info;	25	TGFileInfo file_info;
26	const char *filetypes[] = {"Histograms",histextall,0,0};	26	const char *filetypes[] = {"Histograms",histextall,0,0};
27	char *cTemp;	27	char *cTemp;
28	file_info.fFileTypes = filetypes;	28	file_info.fFileTypes = filetypes;
29	cTemp = new char[1024];	29	// cTemp = new char[1024];
30	sprintf(cTemp, "%s/results", rootdir);	30	// sprintf(cTemp, "%s/results", rootdir);
31	file_info.fIniDir = StrDup(cTemp);	31	// file_info.fIniDir = StrDup(cTemp);
-		32	file_info.fIniDir = StrDup(currentOpenDir);
32	file_info.fMultipleSelection = kFALSE;	33	file_info.fMultipleSelection = kFALSE;
33	new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);	34	new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
34	delete[] cTemp;	35	// delete[] cTemp;
35		36
36	if(file_info.fFilename != NULL)	37	if(file_info.fFilename != NULL)
37	{	38	{
38	darkRun->widgetTE->SetText(file_info.fFilename);	39	darkRun->widgetTE->SetText(file_info.fFilename);
39	fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());	40	fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
40	fileList->Layout();	41	fileList->Layout();
41	}	42	}
42	else	43	else
43	darkRun->widgetTE->SetText("");	44	darkRun->widgetTE->SetText("");
44	}	45	}
45		46
46	// Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)	47	// Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
47	void TGAppMainFrame::AnalysisDefaults()	48	void TGAppMainFrame::AnalysisDefaults()
48	{	49	{
49	printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());	50	printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
50	if(analTab->GetCurrent() == 0) // Integrate spectrum	51	if(analTab->GetCurrent() == 0) // Integrate spectrum
51	{	52	{
52	intSpect->widgetChBox[0]->SetState(kButtonUp);	53	intSpect->widgetChBox[0]->SetState(kButtonUp);
53	intSpect->widgetChBox[1]->SetState(kButtonUp);	54	intSpect->widgetChBox[1]->SetState(kButtonUp);
54	intSpect->widgetChBox[2]->SetState(kButtonDown);	55	intSpect->widgetChBox[2]->SetState(kButtonDown);
55	resol2d->widgetNE[0]->SetNumber(40);	56	resol2d->widgetNE[0]->SetNumber(40);
56	resol2d->widgetNE[1]->SetNumber(40);	57	resol2d->widgetNE[1]->SetNumber(40);
57	}	58	}
58	else if(analTab->GetCurrent() == 1) // Relative PDE	59	else if(analTab->GetCurrent() == 1) // Relative PDE
59	{	60	{
60	relPde->widgetChBox[1]->SetState(kButtonDown);	61	relPde->widgetChBox[1]->SetState(kButtonDown);
61	midPeak->widgetChBox[0]->SetState(kButtonUp);	62	midPeak->widgetChBox[0]->SetState(kButtonUp);
62	zeroAngle->widgetNE[0]->SetNumber(0.00);	63	zeroAngle->widgetNE[0]->SetNumber(0.00);
63	}	64	}
64	else if(analTab->GetCurrent() == 2) // Breakdown voltage	65	else if(analTab->GetCurrent() == 2) // Breakdown voltage
65	{	66	{
66	minPeak->widgetNE[0]->SetNumber(2);	67	minPeak->widgetNE[0]->SetNumber(2);
67	minPeak->widgetNE[0]->SetNumber(1);	68	minPeak->widgetNE[0]->SetNumber(1);
68	}	69	}
69	else if(analTab->GetCurrent() == 3) // Surface scan	70	else if(analTab->GetCurrent() == 3) // Surface scan
70	{	71	{
71	surfScanOpt->widgetChBox[0]->SetState(kButtonDown);	72	surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
72	surfScanOpt->widgetChBox[1]->SetState(kButtonUp);	73	surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
73	resolSurf->widgetNE[0]->SetNumber(40);	74	resolSurf->widgetNE[0]->SetNumber(40);
74	resolSurf->widgetNE[1]->SetNumber(40);	75	resolSurf->widgetNE[1]->SetNumber(40);
75	}	76	}
76	}	77	}
77		78
78	// Analysis handle function	79	// Analysis handle function
79	void TGAppMainFrame::AnalysisHandle(int type)	80	void TGAppMainFrame::AnalysisHandle(int type)
80	{	81	{
81	TList *files;	82	TList *files;
82	bool createTab = true;	83	bool createTab = true;
83	int tabid = -1;	84	int tabid = -1;
84	int analtab = analTab->GetCurrent();	85	int analtab = analTab->GetCurrent();
85		86
86	int analtype;	87	int analtype;
87	if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )	88	if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
88	analtype = 0;	89	analtype = 0;
89	else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() \|\| intSpect->widgetChBox[1]->IsDown()) )	90	else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() \|\| intSpect->widgetChBox[1]->IsDown()) )
90	analtype = 1;	91	analtype = 1;
91	else if( analtab == 1 )	92	else if( analtab == 1 )
92	analtype = 2;	93	analtype = 2;
93	else if( analtab == 2 )	94	else if( analtab == 2 )
94	analtype = 3;	95	analtype = 3;
95	else if( analtab == 3 )	96	else if( analtab == 3 )
96	analtype = 4;	97	analtype = 4;
97		98
98	// Only integrate spectrum or make relative PDE	99	// Only integrate spectrum or make relative PDE
99	if(type == 0)	100	if(type == 0)
100	{	101	{
101	files = new TList();	102	files = new TList();
102	fileList->GetSelectedEntries(files);	103	fileList->GetSelectedEntries(files);
103		104
104	if( analtype == 0 )	105	if( analtype == 0 )
105	IntegSpectrum(files, 0, 0);	106	IntegSpectrum(files, 0, 0);
106		107
107	if( intSpect->widgetChBox[0]->IsDown() )	108	if( intSpect->widgetChBox[0]->IsDown() )
108	IntegSpectrum(files, 1, 0);	109	IntegSpectrum(files, 1, 0);
109		110
110	if( intSpect->widgetChBox[1]->IsDown() )	111	if( intSpect->widgetChBox[1]->IsDown() )
111	IntegSpectrum(files, 2, 0);	112	IntegSpectrum(files, 2, 0);
112		113
113	if( analtype == 2 )	114	if( analtype == 2 )
114	PhotonMu(files, 0);	115	PhotonMu(files, 0);
115		116
116	if( analtype == 3 )	117	if( analtype == 3 )
117	BreakdownVolt(files, 0);	118	BreakdownVolt(files, 0);
118		119
119	if( analtype == 4 )	120	if( analtype == 4 )
120	SurfaceScan(files, 0);	121	SurfaceScan(files, 0);
121	}	122	}
122	// Integrate spectrum or make relative PDE and open edit window	123	// Integrate spectrum or make relative PDE and open edit window
123	else if(type == 1)	124	else if(type == 1)
124	{	125	{
125	files = new TList();	126	files = new TList();
126	fileList->GetSelectedEntries(files);	127	fileList->GetSelectedEntries(files);
127		128
128	// Prepare a new analysis edit tab, if we want to edit plots	129	// Prepare a new analysis edit tab, if we want to edit plots
129	for(int i = 0; i < fTab->GetNumberOfTabs(); i++)	130	for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
130	{	131	{
131	if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)	132	if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
132	{	133	{
133	createTab = false;	134	createTab = false;
134	tabid = i;	135	tabid = i;
135	}	136	}
136		137
137	if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );	138	if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
138	}	139	}
139		140
140	if(files->GetSize() > 0)	141	if(files->GetSize() > 0)
141	{	142	{
142	TempAnalysisTab(fTab, createTab, &tabid, analtype);	143	TempAnalysisTab(fTab, createTab, &tabid, analtype);
143		144
144	// Integrate spectra	145	// Integrate spectra
145	if( analtype == 0 )	146	if( analtype == 0 )
146	IntegSpectrum(files, 0, 1);	147	IntegSpectrum(files, 0, 1);
147		148
148	if( intSpect->widgetChBox[0]->IsDown() )	149	if( intSpect->widgetChBox[0]->IsDown() )
149	IntegSpectrum(files, 1, 1);	150	IntegSpectrum(files, 1, 1);
150		151
151	if( intSpect->widgetChBox[1]->IsDown() )	152	if( intSpect->widgetChBox[1]->IsDown() )
152	IntegSpectrum(files, 2, 1);	153	IntegSpectrum(files, 2, 1);
153		154
154	if( analtype == 2 )	155	if( analtype == 2 )
155	PhotonMu(files, 1);	156	PhotonMu(files, 1);
156		157
157	if( analtype == 3 )	158	if( analtype == 3 )
158	BreakdownVolt(files, 1);	159	BreakdownVolt(files, 1);
159		160
160	if( analtype == 4 )	161	if( analtype == 4 )
161	SurfaceScan(files, 1);	162	SurfaceScan(files, 1);
162		163
163	fTab->SetTab(tabid);	164	fTab->SetTab(tabid);
164	}	165	}
165	}	166	}
166		167
167	delete files;	168	delete files;
168	}	169	}
169		170
170	// Analysis functions ----------------------------------	171	// Analysis functions ----------------------------------
171	void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)	172	void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
172	{	173	{
173	unsigned int nrfiles = files->GetSize();	174	unsigned int nrfiles = files->GetSize();
174	char ctemp[1024];	175	char ctemp[1024];
175	char exportname[1024];	176	char exportname[1024];
176	int j, k = 0, m = 0;	177	int j, k = 0, m = 0;
177		178
178	TTree header_data, meas_data;	179	TTree header_data, meas_data;
179	double integralCount, integralAcc;	180	double integralCount, integralAcc;
180	integralCount = new double[nrfiles];	181	integralCount = new double[nrfiles];
181	integralAcc = new double[nrfiles];	182	integralAcc = new double[nrfiles];
182	double surfxy, surfz;	183	double surfxy, surfz;
183	surfxy = new double[nrfiles];	184	surfxy = new double[nrfiles];
184	surfz = new double[nrfiles];	185	surfz = new double[nrfiles];
185	double minInteg, maxInteg;	186	double minInteg, maxInteg;
186	bool norm = intSpect->widgetChBox[2]->IsDown();	187	bool norm = intSpect->widgetChBox[2]->IsDown();
187	double curzval;	188	double curzval;
188	bool edge2d = false;	189	bool edge2d = false;
189	TCanvas *gCanvas;	190	TCanvas *gCanvas;
190		191
191	float progVal = 0;	192	float progVal = 0;
192	analysisProgress->widgetPB->SetPosition(progVal);	193	analysisProgress->widgetPB->SetPosition(progVal);
193	gVirtualX->Update(1);	194	gVirtualX->Update(1);
194		195
195	// Zero the integral count and accumulated vaules	196	// Zero the integral count and accumulated vaules
196	for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }	197	for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
197		198
198	// Start if we select at least one file	199	// Start if we select at least one file
199	if(nrfiles > 0)	200	if(nrfiles > 0)
200	{	201	{
201	for(int i = 0; i < (int)nrfiles; i++)	202	for(int i = 0; i < (int)nrfiles; i++)
202	{	203	{
203	if(files->At(i))	204	if(files->At(i))
204	{	205	{
205	// Read the X,Y and Z positions from header and ADC and TDC values from the measurements	206	// Read the X,Y and Z positions from header and ADC and TDC values from the measurements
206	sprintf(ctemp, "%s", files->At(i)->GetTitle());	207	sprintf(ctemp, "%s", files->At(i)->GetTitle());
207	inroot = new TFile(ctemp, "READ");	208	inroot = new TFile(ctemp, "READ");
208		209
209	inroot->GetObject("header_data", header_data);	210	inroot->GetObject("header_data", header_data);
210	inroot->GetObject("meas_data", meas_data);	211	inroot->GetObject("meas_data", meas_data);
211		212
212	header_data->SetBranchAddress("xpos", &evtheader.xpos);	213	header_data->SetBranchAddress("xpos", &evtheader.xpos);
213	header_data->GetEntry(0);	214	header_data->GetEntry(0);
214	header_data->SetBranchAddress("ypos", &evtheader.ypos);	215	header_data->SetBranchAddress("ypos", &evtheader.ypos);
215	header_data->GetEntry(0);	216	header_data->GetEntry(0);
216	header_data->SetBranchAddress("zpos", &evtheader.zpos);	217	header_data->SetBranchAddress("zpos", &evtheader.zpos);
217	header_data->GetEntry(0);	218	header_data->GetEntry(0);
218		219
219	char rdc[256];	220	char rdc[256];
220	j = selectCh->widgetNE[0]->GetNumber();	221	j = selectCh->widgetNE[0]->GetNumber();
221	double rangetdc[2];	222	double rangetdc[2];
222	rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();	223	rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
223	rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();	224	rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
224		225
225	k = 0;	226	k = 0;
226	m = 0;	227	m = 0;
227		228
228	// Reading the data	229	// Reading the data
229	for(int e = 0; e < meas_data->GetEntries(); e++)	230	for(int e = 0; e < meas_data->GetEntries(); e++)
230	{	231	{
231	sprintf(rdc, "ADC%d", j);	232	sprintf(rdc, "ADC%d", j);
232	meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);	233	meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
233	meas_data->GetEntry(e);	234	meas_data->GetEntry(e);
234		235
235	sprintf(rdc, "TDC%d", j);	236	sprintf(rdc, "TDC%d", j);
236	meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);	237	meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
237	meas_data->GetEntry(e);	238	meas_data->GetEntry(e);
238		239
239	// Use data point only if it is inside the TDC window	240	// Use data point only if it is inside the TDC window
240	if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )	241	if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
241	{	242	{
242	k++;	243	k++;
243	m += evtdata.adcdata[j];	244	m += evtdata.adcdata[j];
244	}	245	}
245	}	246	}
246		247
247	// X, Y and Z values from each file (table units or microns)	248	// X, Y and Z values from each file (table units or microns)
248	if(~~posUnits~~->widgetCB->GetSelected() == 0)	249	if(posUnitsPlot->widgetCB->GetSelected() == 0)
249	{	250	{
250	if(direction == 1)	251	if(direction == 1)
251	surfxy[i] = (double)(evtheader.xpos);	252	surfxy[i] = (double)(evtheader.xpos);
252	else if(direction == 2)	253	else if(direction == 2)
253	surfxy[i] = (double)(evtheader.ypos);	254	surfxy[i] = (double)(evtheader.ypos);
254	surfz[i] = (double)(evtheader.zpos);	255	surfz[i] = (double)(evtheader.zpos);
255	}	256	}
256	else if(~~posUnits~~->widgetCB->GetSelected() == 1)	257	else if(posUnitsPlot->widgetCB->GetSelected() == 1)
257	{	258	{
258	if(direction == 1)	259	if(direction == 1)
259	surfxy[i] = (double)(evtheader.xpos*lenconversion);	260	surfxy[i] = (double)(evtheader.xpos*lenconversion);
260	else if(direction == 2)	261	else if(direction == 2)
261	surfxy[i] = (double)(evtheader.ypos*lenconversion);	262	surfxy[i] = (double)(evtheader.ypos*lenconversion);
262	surfz[i] = (double)(evtheader.zpos*lenconversion);	263	surfz[i] = (double)(evtheader.zpos*lenconversion);
263	}	264	}
264		265
265	// 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	266	// Check if we have different Z values: if no, just make the edge plots; if yes, save edge plots and make a 2d edge plot
266	if(i == 0) curzval = surfz[i];	267	if(i == 0) curzval = surfz[i];
267	else	268	else
268	{	269	{
269	if(surfz[i] != curzval)	270	if(surfz[i] != curzval)
270	edge2d = true;	271	edge2d = true;
271	}	272	}
272		273
273	// Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting	274	// Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
274	if(direction == 0)	275	if(direction == 0)
275	{	276	{
276	if(norm)	277	if(norm)
277	{	278	{
278	integralCount[i] += ((double)m)/((double)k);	279	integralCount[i] += ((double)m)/((double)k);
279	printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);	280	printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
280	}	281	}
281	else	282	else
282	{	283	{
283	integralCount[i] += m;	284	integralCount[i] += m;
284	printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);	285	printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
285	}	286	}
286	}	287	}
287	else	288	else
288	{	289	{
289	if(norm)	290	if(norm)
290	integralCount[i] += ((double)m)/((double)k);	291	integralCount[i] += ((double)m)/((double)k);
291	else	292	else
292	integralCount[i] += m;	293	integralCount[i] += m;
293	}	294	}
294		295
295	inroot->Close();	296	inroot->Close();
296	delete inroot;	297	delete inroot;
297	}	298	}
298		299
299	// Update the progress bar	300	// Update the progress bar
300	progVal = (float)(75.00/nrfiles)*i;	301	progVal = (float)(75.00/nrfiles)*i;
301	analysisProgress->widgetPB->SetPosition(progVal);	302	analysisProgress->widgetPB->SetPosition(progVal);
302	gVirtualX->Update(1);	303	gVirtualX->Update(1);
303	}	304	}
304		305
305	printf("IntegSpectrum(): %d files were selected.\n", nrfiles);	306	printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
306		307
307	// If only an integral is needed, do not plot and exit here	308	// If only an integral is needed, do not plot and exit here
308	if( direction == 0 )	309	if( direction == 0 )
309	{	310	{
310	delete[] integralCount;	311	delete[] integralCount;
311	delete[] surfxy;	312	delete[] surfxy;
312	delete[] surfz;	313	delete[] surfz;
313	return;	314	return;
314	}	315	}
315		316
316	// Current z value and the accumulated counter	317	// Current z value and the accumulated counter
317	curzval = surfz[0];	318	curzval = surfz[0];
318	j = 0;	319	j = 0;
319	int acc = 0;	320	int acc = 0;
320	int zb;	321	int zb;
321	for(int i = 0; i <= (int)nrfiles; i++)	322	for(int i = 0; i <= (int)nrfiles; i++)
322	{	323	{
323	// Collect the accumulated integral in order to produce a PDF from a CDF	324	// Collect the accumulated integral in order to produce a PDF from a CDF
324	// While we are at the same Z value, save under one set	325	// While we are at the same Z value, save under one set
325	if( (surfz[i] == curzval) && (acc != nrfiles) )	326	if( (surfz[i] == curzval) && (acc != nrfiles) )
326	{	327	{
327	integralAcc[j] = integralCount[i];	328	integralAcc[j] = integralCount[i];
328	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]);	329	if(DBGSIG) printf("IntegSpectrum(): Integral check 1 (i=%d,j=%d,z=%.2lf): %lf\t%lf\n", i, j, surfz[i], integralCount[i], integralAcc[j]);
329	j++;	330	j++;
330	acc++;	331	acc++;
331	}	332	}
332	// 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	333	// When we switch to a new set of Z values and at the end, we must save the previous ones to make 1D edge plots
333	else	334	else
334	{	335	{
335	// Find minimal and maximal integral values to subtract the offset and normate PDF to 1	336	// Find minimal and maximal integral values to subtract the offset and normate PDF to 1
336	NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);	337	NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
337		338
338	if(acc != nrfiles)	339	if(acc != nrfiles)
339	{	340	{
340	curzval = surfz[i];	341	curzval = surfz[i];
341	// PDF and CDF plot	342	// PDF and CDF plot
342	PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);	343	PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
343	i--;	344	i--;
344	j = 0;	345	j = 0;
345	}	346	}
346	else	347	else
347	{	348	{
348	// PDF and CDF plot	349	// PDF and CDF plot
349	PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);	350	PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
350	i--;	351	i--;
351	break;	352	break;
352	}	353	}
353	}	354	}
354		355
355	// Update the progress bar	356	// Update the progress bar
356	progVal = (float)(15.00/nrfiles)*i+75.00;	357	progVal = (float)(15.00/nrfiles)*i+75.00;
357	analysisProgress->widgetPB->SetPosition(progVal);	358	analysisProgress->widgetPB->SetPosition(progVal);
358	gVirtualX->Update(1);	359	gVirtualX->Update(1);
359	}	360	}
360		361
361	// Make the 2D edge plot	362	// Make the 2D edge plot
362	if(edge2d)	363	if(edge2d)
363	{	364	{
364	if(edit == 0)	365	if(edit == 0)
365	gCanvas = new TCanvas("canv","canv",900,900);	366	gCanvas = new TCanvas("canv","canv",900,900);
366	else	367	else
367	gCanvas = tempAnalysisCanvas->GetCanvas();	368	gCanvas = tempAnalysisCanvas->GetCanvas();
368		369
369	double range[4];	370	double range[4];
370	TGraph2D *gScan2D;	371	TGraph2D *gScan2D;
371	gScan2D = new TGraph2D();	372	gScan2D = new TGraph2D();
372	range[0] = TMath::MinElement(nrfiles, surfxy);	373	range[0] = TMath::MinElement(nrfiles, surfxy);
373	range[1] = TMath::MaxElement(nrfiles, surfxy);	374	range[1] = TMath::MaxElement(nrfiles, surfxy);
374	range[2] = TMath::MinElement(nrfiles, surfz);	375	range[2] = TMath::MinElement(nrfiles, surfz);
375	range[3] = TMath::MaxElement(nrfiles, surfz);	376	range[3] = TMath::MaxElement(nrfiles, surfz);
376		377
377	for(int i = 0; i < nrfiles; i++)	378	for(int i = 0; i < nrfiles; i++)
378	{	379	{
379	if(DBGSIG)	380	if(DBGSIG)
380	printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);	381	printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
381	gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);	382	gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
382		383
383	// Update the progress bar	384	// Update the progress bar
384	progVal = (float)(9.00/nrfiles)*i+90.00;	385	progVal = (float)(9.00/nrfiles)*i+90.00;
385	analysisProgress->widgetPB->SetPosition(progVal);	386	analysisProgress->widgetPB->SetPosition(progVal);
386	gVirtualX->Update(1);	387	gVirtualX->Update(1);
387	}	388	}
388		389
389	gCanvas->cd();	390	gCanvas->cd();
390	gStyle->SetPalette(1);	391	gStyle->SetPalette(1);
391	gCanvas->SetLeftMargin(0.15);	392	gCanvas->SetLeftMargin(0.15);
392	gCanvas->SetRightMargin(0.126);	393	gCanvas->SetRightMargin(0.126);
393	gCanvas->SetTopMargin(0.077);	394	gCanvas->SetTopMargin(0.077);
394	gScan2D->Draw("COLZ");	395	gScan2D->Draw("COLZ");
395		396
396	gCanvas->Modified();	397	gCanvas->Modified();
397	gCanvas->Update();	398	gCanvas->Update();
398		399
399	gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());	400	gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
400	gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());	401	gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
401		402
402	gCanvas->Modified();	403	gCanvas->Modified();
403	gCanvas->Update();	404	gCanvas->Update();
404		405
405	if(direction == 1)	406	if(direction == 1)
406	gScan2D->GetXaxis()->SetTitle("X [#mum]");	407	gScan2D->GetXaxis()->SetTitle("X [#mum]");
407	else if(direction == 2)	408	else if(direction == 2)
408	gScan2D->GetXaxis()->SetTitle("Y [#mum]");	409	gScan2D->GetXaxis()->SetTitle("Y [#mum]");
409		410
410		411
411	gScan2D->GetXaxis()->SetTitleOffset(1.3);	412	gScan2D->GetXaxis()->SetTitleOffset(1.3);
412	gScan2D->GetXaxis()->CenterTitle(kTRUE);	413	gScan2D->GetXaxis()->CenterTitle(kTRUE);
413	gScan2D->GetXaxis()->SetLabelSize(0.027);	414	gScan2D->GetXaxis()->SetLabelSize(0.027);
414	gScan2D->GetXaxis()->SetLabelOffset(0.02);	415	gScan2D->GetXaxis()->SetLabelOffset(0.02);
415	gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);	416	gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
416	gScan2D->GetXaxis()->SetNoExponent(kTRUE);	417	gScan2D->GetXaxis()->SetNoExponent(kTRUE);
417	gScan2D->GetYaxis()->SetTitleOffset(1.9);	418	gScan2D->GetYaxis()->SetTitleOffset(1.9);
418	gScan2D->GetYaxis()->CenterTitle(kTRUE);	419	gScan2D->GetYaxis()->CenterTitle(kTRUE);
419	gScan2D->GetYaxis()->SetLabelSize(0.027);	420	gScan2D->GetYaxis()->SetLabelSize(0.027);
420	gScan2D->~~GetXaxis~~()->SetLabelOffset(0.02);	421	gScan2D->GetYaxis()->SetLabelOffset(0.02);
421	gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);	422	gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
422	gScan2D->GetYaxis()->SetNoExponent(kTRUE);	423	gScan2D->GetYaxis()->SetNoExponent(kTRUE);
423		424
424	/* TGaxis yax = (TGaxis)gScan2D->GetYaxis();	425	/* TGaxis yax = (TGaxis)gScan2D->GetYaxis();
425	yax->SetMaxDigits(4);*/	426	yax->SetMaxDigits(4);*/
426		427
427	if(!cleanPlots)	428	if(!cleanPlots)
428	{	429	{
429	if(direction == 1)	430	if(direction == 1)
430	{	431	{
431	if(~~posUnits~~->widgetCB->GetSelected() == 0)	432	if(posUnitsPlot->widgetCB->GetSelected() == 0)
432	gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");	433	gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
433	else if(~~posUnits~~->widgetCB->GetSelected() == 1)	434	else if(posUnitsPlot->widgetCB->GetSelected() == 1)
434	gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");	435	gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
435	}	436	}
436	else if(direction == 2)	437	else if(direction == 2)
437	{	438	{
438	if(~~posUnits~~->widgetCB->GetSelected() == 0)	439	if(posUnitsPlot->widgetCB->GetSelected() == 0)
439	gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");	440	gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
440	else if(~~posUnits~~->widgetCB->GetSelected() == 1)	441	else if(posUnitsPlot->widgetCB->GetSelected() == 1)
441	gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");	442	gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
442	}	443	}
443	}	444	}
444	else	445	else
445	{	446	{
446	if(direction == 1)	447	if(direction == 1)
447	{	448	{
448	if(~~posUnits~~->widgetCB->GetSelected() == 0)	449	if(posUnitsPlot->widgetCB->GetSelected() == 0)
449	gScan2D->SetTitle(";X [table units];Z [table units]");	450	gScan2D->SetTitle(";X [table units];Z [table units]");
450	else if(~~posUnits~~->widgetCB->GetSelected() == 1)	451	else if(posUnitsPlot->widgetCB->GetSelected() == 1)
451	gScan2D->SetTitle(";X [#mum];Z [#mum]");	452	gScan2D->SetTitle(";X [#mum];Z [#mum]");
452	}	453	}
453	else if(direction == 2)	454	else if(direction == 2)
454	{	455	{
455	if(~~posUnits~~->widgetCB->GetSelected() == 0)	456	if(posUnitsPlot->widgetCB->GetSelected() == 0)
456	gScan2D->SetTitle(";Y [table units];Z [table units]");	457	gScan2D->SetTitle(";Y [table units];Z [table units]");
457	else if(~~posUnits~~->widgetCB->GetSelected() == 1)	458	else if(posUnitsPlot->widgetCB->GetSelected() == 1)
458	gScan2D->SetTitle(";Y [#mum];Z [#mum]");	459	gScan2D->SetTitle(";Y [#mum];Z [#mum]");
459	}	460	}
460	}	461	}
461		462
462	gCanvas->Modified();	463	gCanvas->Modified();
463	gCanvas->Update();	464	gCanvas->Update();
464		465
465	remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);	466	remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
466	sprintf(exportname, "%s", ctemp);	467	sprintf(exportname, "%s", ctemp);
467	remove_from_last(exportname, '_', ctemp);	468	remove_from_last(exportname, '_', ctemp);
468	if(direction == 1)	469	if(direction == 1)
469	sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);	470	sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
470	else if(direction == 2)	471	else if(direction == 2)
471	sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);	472	sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
472	gCanvas->SaveAs(exportname);	473	gCanvas->SaveAs(exportname);
473		474
474	// Update the progress bar	475	// Update the progress bar
475	analysisProgress->widgetPB->SetPosition(100.0);	476	analysisProgress->widgetPB->SetPosition(100.0);
476	gVirtualX->Update(1);	477	gVirtualX->Update(1);
477		478
478	if(edit == 0)	479	if(edit == 0)
479	{	480	{
480	delete gScan2D;	481	delete gScan2D;
481	delete gCanvas;	482	delete gCanvas;
482	}	483	}
483	}	484	}
484	else	485	else
485	{	486	{
486	// Update the progress bar	487	// Update the progress bar
487	analysisProgress->widgetPB->SetPosition(100.0);	488	analysisProgress->widgetPB->SetPosition(100.0);
488	gVirtualX->Update(1);	489	gVirtualX->Update(1);
489	}	490	}
490	}	491	}
491	}	492	}
492		493
493	void TGAppMainFrame::PlotEdgeDistribution(TList files, int filenr, int points, double min, double max, double xy, double *integAcc, int axis, bool edge2d, int edit)	494	void TGAppMainFrame::PlotEdgeDistribution(TList files, int filenr, int points, double min, double max, double xy, double *integAcc, int axis, bool edge2d, int edit)
494	{	495	{
495	TGraph *gScan[2];	496	TGraph *gScan[2];
496	int pdfmax = -1;	497	int pdfmax = -1;
497	int count = 0;	498	int count = 0;
498	char ctemp[1024];	499	char ctemp[1024];
499	char exportname[1024];	500	char exportname[1024];

Subversion Repositories f9daq

(root)/lab/sipmscan/trunk/src/analysis.cpp – Rev 146 → 167