Subversion Repositories f9daq

#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.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(posUnits->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(posUnits->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->GetXaxis()->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(posUnits->widgetCB->GetSelected() == 0)
                  gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
               else if(posUnits->widgetCB->GetSelected() == 1)
                  gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
            }
            else if(direction == 2)
            {
               if(posUnits->widgetCB->GetSelected() == 0)
                  gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
               else if(posUnits->widgetCB->GetSelected() == 1)
                  gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
            }
         }
         else
         {
            if(direction == 1)
            {
               if(posUnits->widgetCB->GetSelected() == 0)
                  gScan2D->SetTitle(";X [table units];Z [table units]");
               else if(posUnits->widgetCB->GetSelected() == 1)
                  gScan2D->SetTitle(";X [#mum];Z [#mum]");
            }
            else if(direction == 2)
            {
               if(posUnits->widgetCB->GetSelected() == 0)
                  gScan2D->SetTitle(";Y [table units];Z [table units]");
               else if(posUnits->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(posUnits->widgetCB->GetSelected() == 0)
            gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
         else if(posUnits->widgetCB->GetSelected() == 1)
            gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
      }
      else if(axis == 2)
      {
         if(posUnits->widgetCB->GetSelected() == 0)
            gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
         else if(posUnits->widgetCB->GetSelected() == 1)
            gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
      }
   }
   else
   {
      if(axis == 1)
      {
         if(posUnits->widgetCB->GetSelected() == 0)
            gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
         else if(posUnits->widgetCB->GetSelected() == 1)
            gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
      }
      else if(axis == 2)
      {
         if(posUnits->widgetCB->GetSelected() == 0)
            gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
         else if(posUnits->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(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]);
      }

      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,90);
      pde->GetXaxis()->SetRangeUser(-90,90);
      pde->GetXaxis()->SetLimits(-90,90);
      pde->GetYaxis()->SetTitleOffset(1.2);
      pde->GetYaxis()->SetLabelSize(0.030);
      pde->GetYaxis()->CenterTitle();
      pde->GetYaxis()->SetRangeUser(0., 1.2);
      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(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);
      }
      else
         printf("#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])) );

      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 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;

   // 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(posUnits->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(posUnits->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);

//      // 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 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(posUnits->widgetCB->GetSelected() == 0)
            gScan2D->SetTitle("Surface scan;X [table units];Y [table units]");
         else if(posUnits->widgetCB->GetSelected() == 1)
            gScan2D->SetTitle("Surface scan;X [#mum];Y [#mum]");
      }
      else
      {
         if(posUnits->widgetCB->GetSelected() == 0)
            gScan2D->SetTitle(";X [table units];Y [table units]");
         else if(posUnits->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->GetXaxis()->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();
      }
   }
}