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

+{
    TGFileInfo file_info;
    const char *filetypes[] = {"Histograms",histextall,0,0};
    char *cTemp;
    file_info.fFileTypes = filetypes;
-   cTemp = new char[1024];
+//   cTemp = new char[1024];
-   sprintf(cTemp, "%s/results", rootdir);
+//   sprintf(cTemp, "%s/results", rootdir);
-   file_info.fIniDir = StrDup(cTemp);
+//   file_info.fIniDir = StrDup(cTemp);
+   file_info.fIniDir = StrDup(currentOpenDir);
    file_info.fMultipleSelection = kFALSE;
    new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
-   delete[] cTemp;
+//   delete[] cTemp;
    if(file_info.fFilename != NULL)
+   {
       darkRun->widgetTE->SetText(file_info.fFilename);
       fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
                   m += evtdata.adcdata[j];
+               }
+            }
             // X, Y and Z values from each file (table units or microns)
-            if(posUnits->widgetCB->GetSelected() == 0)
+            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(posUnits->widgetCB->GetSelected() == 1)
+            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);
          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()->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)
+               if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
-               else if(posUnits->widgetCB->GetSelected() == 1)
+               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
             }
             else if(direction == 2)
+            {
-               if(posUnits->widgetCB->GetSelected() == 0)
+               if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
-               else if(posUnits->widgetCB->GetSelected() == 1)
+               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
             }
          }
          else
+         {
             if(direction == 1)
+            {
-               if(posUnits->widgetCB->GetSelected() == 0)
+               if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle(";X [table units];Z [table units]");
-               else if(posUnits->widgetCB->GetSelected() == 1)
+               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle(";X [#mum];Z [#mum]");
+            }
             else if(direction == 2)
+            {
-               if(posUnits->widgetCB->GetSelected() == 0)
+               if(posUnitsPlot->widgetCB->GetSelected() == 0)
                   gScan2D->SetTitle(";Y [table units];Z [table units]");
-               else if(posUnits->widgetCB->GetSelected() == 1)
+               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
                   gScan2D->SetTitle(";Y [#mum];Z [#mum]");
+            }
+         }
          gCanvas->Modified();
          if(edit == 0)
+         {
             delete gScan2D;
             delete gCanvas;
+         }
       }
       else
+      {
          // Update the progress bar
          analysisProgress->widgetPB->SetPosition(100.0);
          gVirtualX->Update(1);
    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
    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);
+   }
    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);
    if(!cleanPlots)
+   {
       if(axis == 1)
+      {
-         if(posUnits->widgetCB->GetSelected() == 0)
+         if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
-         else if(posUnits->widgetCB->GetSelected() == 1)
+         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
+      }
       else if(axis == 2)
+      {
-         if(posUnits->widgetCB->GetSelected() == 0)
+         if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
-         else if(posUnits->widgetCB->GetSelected() == 1)
+         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
       }
+   }
    else
+   {
       if(axis == 1)
       {
-         if(posUnits->widgetCB->GetSelected() == 0)
+         if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
-         else if(posUnits->widgetCB->GetSelected() == 1)
+         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
+      }
       else if(axis == 2)
+      {
-         if(posUnits->widgetCB->GetSelected() == 0)
+         if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
-         else if(posUnits->widgetCB->GetSelected() == 1)
+         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);
    // Start if we select at least one file
    if(nrfiles > 0)
+   {
       for(int i = 0; i < (int)nrfiles; i++)
+      {
-         if(files->At(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();
+                  {
                      // 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]];
             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);
             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(!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->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()->SetRange(-90.0,90.0);
-      pde->GetXaxis()->SetRangeUser(-90,90);
+      pde->GetXaxis()->SetRangeUser(-90.0,90.0);
-      pde->GetXaxis()->SetLimits(-90,90);
+      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.2);
+      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)");
    // 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();
                   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;
       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);
+      {
          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", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
+         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);
    double *surfx, *surfy;
    surfx = new double[nrfiles];
    surfy = new double[nrfiles];
    double xsurfmin = 0, ysurfmin = 0;
    int nrentries;
-//   double minInteg, maxInteg;
+   double minInteg, maxInteg;
    bool norm = surfScanOpt->widgetChBox[0]->IsDown();
-//   double curzval;
+   double curyval;
 //   bool edge2d = false;
    TCanvas *gCanvas;
    float progVal = 0;
                   m += evtdata.adcdata[j];
+               }
+            }
             // X, Y and Z values from each file (table units or microns)
-            if(posUnits->widgetCB->GetSelected() == 0)
+            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(posUnits->widgetCB->GetSelected() == 1)
+            else if(posUnitsPlot->widgetCB->GetSelected() == 1)
+            {
                if(n == 1)
+               {
                   xsurfmin = (double)(evtheader.xpos*lenconversion);
                   ysurfmin = (double)(evtheader.ypos*lenconversion);
+      }
       nrentries = n;
       printf("SurfaceScan(): %d files were selected.\n", nrfiles);
-//      // If only an integral is needed, do not plot and exit here
+      // Check for Minimum and Maximum values and normalize to 1, if normalization is selected
-//      if( direction == 0 )
+      if(norm)
-//      {
+      {
-//         delete[] integralCount;
+         minInteg = TMath::MinElement(nrfiles, integralCount);
-//         delete[] surfxy;
+         for(int i = 0; i < nrfiles; i++)
+         {
-//         delete[] surfz;
+            integralCount[i] -= minInteg;
-//         return;
+            if(DBGSIG) printf("Subtraction: %lf\n", integralCount[i]);
-//      }
+         }
-//
-//      // Current z value and the accumulated counter
+         maxInteg = TMath::MaxElement(nrfiles, integralCount);
-//      curzval = surfz[0];
-//      j = 0;
-//      int acc = 0;
-//      int zb;
-//      for(int i = 0; i <= (int)nrfiles; i++)
+         for(int i = 0; i < 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];
+            integralCount[i] = integralCount[i]/maxInteg;
-//            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(DBGSIG) printf("Normalization: %lf\n", integralCount[i]);
-//
-//            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
+      {
          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);
       gCanvas->Modified();
       gCanvas->Update();
       if(!cleanPlots)
+      {
-         if(posUnits->widgetCB->GetSelected() == 0)
+         if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan2D->SetTitle("Surface scan;X [table units];Y [table units]");
-         else if(posUnits->widgetCB->GetSelected() == 1)
+         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
             gScan2D->SetTitle("Surface scan;X [#mum];Y [#mum]");
+      }
       else
+      {
-         if(posUnits->widgetCB->GetSelected() == 0)
+         if(posUnitsPlot->widgetCB->GetSelected() == 0)
             gScan2D->SetTitle(";X [table units];Y [table units]");
-         else if(posUnits->widgetCB->GetSelected() == 1)
+         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();
       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()->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);

Subversion Repositories f9daq

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