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

       resol2d->widgetNE[1]->SetNumber(40);
+   }
    else if(analTab->GetCurrent() == 1)  // Relative PDE
+   {
       relPde->widgetChBox[0]->SetState(kButtonDown);
+      relPde->widgetChBox[1]->SetState(kButtonUp);
       midPeak->widgetChBox[0]->SetState(kButtonUp);
       zeroAngle->widgetNE[0]->SetNumber(0.00);
+   }
    else if(analTab->GetCurrent() == 2)  // Breakdown voltage
+   {
       analtype = 2;
    else if( analtab == 2 )
       analtype = 3;
    else if( analtab == 3 )
       analtype = 4;
+   // Save analysis settings any time we run a new analysis
+   SaveAnalSettings();
    // Only integrate spectrum or make relative PDE
    if(type == 0)
+   {
       files = new TList();
       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( 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);
+{
    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;
       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;
+      }
             gCanvas = tempAnalysisCanvas->GetCanvas();
          double range[4];
          TGraph2D *gScan2D;
          gScan2D = new TGraph2D();
-         gScan2D->SetName("edgescan");
          range[0] = TMath::MinElement(nrfiles, surfxy);
          range[1] = TMath::MaxElement(nrfiles, surfxy);
          range[2] = TMath::MinElement(nrfiles, surfz);
          range[3] = TMath::MaxElement(nrfiles, surfz);
    double meanparam, paramsigma;
    int sortindex[20];
    int adcpedestal[2];
    int zeromu = 0;
    int darkhist = -1;
+   int nopeaks = -1;
    double pointest[12];
    bool exclude = false;
    // Zero the parameter values
    for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
    float progVal = 0;
    analysisProgress->widgetPB->SetPosition(progVal);
    gVirtualX->Update(1);
+   // Check if the checkbox for no peaks is selected - TODO: Still need the situation when we do not have a peaked ADC spectrum
+   if(relPde->widgetChBox[1]->IsDown())
+   {
+      printf("PhotonMu(): ADC spectrum has no peak structure.\n");
+      nopeaks = 1;
+      // Error if there is no darkhist
+      if(strcmp("", darkRun->widgetTE->GetText()) == 0)
+      {
+         printf("PhotonMu(): Error! The no peak structure option needs a dark histogram.\n");
+         delete[] integralCount;
+         delete[] integralPedestal;
+         delete[] angle;
+         delete[] pdeval;
+         delete[] muval;
+         return;
+      }
+   }
    // Start if we select at least one file
    if(nrfiles > 0)
    {
+      // Find the pedestal peak for the dark histogram, and use it for all if there are no peaks
+      if(nopeaks != -1)
+      {
+         // Replot the spectrum on analysisCanvas and do not close the input file
+         DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1);
+         analysisCanvas->GetCanvas()->Modified();
+         analysisCanvas->GetCanvas()->Update();
+         // Get the spectrum
+         histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
+         npeaks = 15;
+         double par[300];
+         spec = new TSpectrum(npeaks);
+         // Find spectrum background
+         histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
+         // Clone histogram and subtract background from it if we select that option
+         h2 = (TH1F*)histtemp->Clone("h2");
+         if(fitChecks->widgetChBox[0]->IsDown())
+            h2->Add(histback, -1);
+         // Search for the peaks
+         int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
+         printf("PhotonMu(): Found %d candidates to fit.\n",found);
+         npeaks = found;
+         // Set initial peak parameters
+         xpeaks = spec->GetPositionX();
+         for(j = 0; j < found; j++)
+         {
+            float xp = xpeaks[j];
+            int bin = h2->GetXaxis()->FindBin(xp);
+            float yp = h2->GetBinContent(bin);
+            par[3*j] = yp;
+            par[3*j+1] = xp;
+            par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
+         }
+         // Fit the histogram
+         fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
+         TVirtualFitter::Fitter(histtemp, 3*npeaks);
+         fit->SetParameters(par);
+         fit->SetNpx(300);
+         h2->Fit("fit","Q");
+         // Get the fitted parameters
+         fittingfunc = h2->GetFunction("fit");
+         fparam = fittingfunc->GetParameters();
+         fparamerr = fittingfunc->GetParErrors();
+         // Gather the parameters (mean peak value for now)
+         int j = 1;
+         int nrfit = 0;
+         while(1)
+         {
+            if( (fparam[j] < 1.E-30) || (nrfit > 8) )
+               break;
+            else
+            {
+               // Check if pedestal is above the lower limit and sigma is smaller than the mean
+               if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+               {
+                  // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
+                  meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
+                  sigmasel[nrfit] = fparam[j+1];
+                  nrfit++;
+               }
+            }
+            j+=3;
+         }
+         TMath::Sort(nrfit, meansel, sortindex, kFALSE);
+         meanparam = meansel[sortindex[0]];
+         paramsigma = sigmasel[sortindex[0]];
+         for(j = 0; j < nrfit; j++)
+            if(DBGSIG)
+               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+         j = 0;
+         adcpedestal[0] = 0;
+         adcpedestal[1] = -1;
+         while(1)
+         {
+            int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
+            int yp = histtemp->GetBinContent(bin);
+            // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak)
+            if(adcpedestal[1] == -1)
+            {
+               adcpedestal[0] = j+meanparam+paramsigma;
+               adcpedestal[1] = yp;
+            }
+            else
+            {
+               if( (npeaks > 1) && (adcpedestal[1] >= yp) )
+               {
+                  adcpedestal[0] = j+meanparam+paramsigma;
+                  adcpedestal[1] = yp;
+               }
+               else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
+               {
+                  adcpedestal[0] = j+meanparam+paramsigma;
+                  adcpedestal[1] = yp;
+               }
+               else
+                  break;
+            }
+            j++;
+            if(j > 50) break;
+         }
+         if(midPeak->widgetChBox[0]->IsDown())
+         {
+            if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
+               m = TMath::Floor(meanparam);
+            else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
+               m = TMath::Ceil(meanparam);
+            int bin = histtemp->GetXaxis()->FindBin(m);
+            adcpedestal[0] = m;
+            printf("midpeak x = %d, ", adcpedestal[0]);
+            adcpedestal[1] = histtemp->GetBinContent(bin);
+         }
+         // Option to show the fit
+         fittingfunc->Draw("L SAME");
+         analysisCanvas->GetCanvas()->Modified();
+         analysisCanvas->GetCanvas()->Update();
+         printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
+         // Delete the opened histogram and spectrum
+         delete spec;
+         delete inroot;
+//       return;
+      }
+      printf("PhotonMu(): Continuing with the rest of the spectra.\n");
+      // Check all histograms for pedestal peak values
       for(int i = 0; i < (int)nrfiles; i++)
+      {
          if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
          {
             printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n");
             // Replot the spectrum on analysisCanvas and do not close the input file
             DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
             analysisCanvas->GetCanvas()->Modified();
                j+=3;
+            }
             TMath::Sort(nrfit, meansel, sortindex, kFALSE);
             fittingfunc->Draw("SAME");
-            analysisCanvas->GetCanvas()->Modified();
-            analysisCanvas->GetCanvas()->Update();
-            meanparam = meansel[sortindex[0]];
-            paramsigma = sigmasel[sortindex[0]];
-            for(j = 0; j < nrfit; j++)
-               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
-            return;
-         }
-         if(files->At(i))
-         {
-            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
-            {
-               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
-               darkhist = i;
-            }
-            // Replot the spectrum on analysisCanvas and do not close the input file
-            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
-            // Get the spectrum
-            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
-            npeaks = 15;
-            double par[300];
-            spec = new TSpectrum(npeaks);
-            // Find spectrum background
-            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
-            // Clone histogram and subtract background from it if we select that option
-            h2 = (TH1F*)histtemp->Clone("h2");
-            if(fitChecks->widgetChBox[0]->IsDown())
-               h2->Add(histback, -1);
-            // Search for the peaks
-            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
-            printf("PhotonMu(): Found %d candidates to fit.\n",found);
-            npeaks = found;
-            // Set initial peak parameters
-            xpeaks = spec->GetPositionX();
-            for(j = 0; j < found; j++)
-            {
-               float xp = xpeaks[j];
-               int bin = h2->GetXaxis()->FindBin(xp);
-               float yp = h2->GetBinContent(bin);
-               par[3*j] = yp;
-               par[3*j+1] = xp;
-               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
-            }
-            // Fit the histogram
-            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
-            TVirtualFitter::Fitter(histtemp, 3*npeaks);
-            fit->SetParameters(par);
-            fit->SetNpx(300);
-            h2->Fit("fit","Q");
-            // Get the fitted parameters
-            fittingfunc = h2->GetFunction("fit");
-            fparam = fittingfunc->GetParameters();
-            fparamerr = fittingfunc->GetParErrors();
-            // Gather the parameters (mean peak value for now)
-            int j = 1;
-            int nrfit = 0;
-            while(1)
-            {
-               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
-                  break;
-               else
-               {
-                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
-                  if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
-                  {
-                     // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
-                     meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
-                     sigmasel[nrfit] = fparam[j+1];
-                     nrfit++;
-                  }
-               }
-               j+=3;
-            }
-            TMath::Sort(nrfit, meansel, sortindex, kFALSE);
             meanparam = meansel[sortindex[0]];
             paramsigma = sigmasel[sortindex[0]];
             for(j = 0; j < nrfit; j++)
-               if(DBGSIG)
-                  printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+            return;
+         }
-            j = 0;
+         if(files->At(i))
+         {
-            adcpedestal[0] = 0;
+            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
+            {
+               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
-            adcpedestal[1] = -1;
+               darkhist = i;
+            }
-            while(1)
+            if(nopeaks == -1)
             {
+               // Replot the spectrum on analysisCanvas and do not close the input file
+               DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
-               int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
+               analysisCanvas->GetCanvas()->Modified();
-               int yp = histtemp->GetBinContent(bin);
+               analysisCanvas->GetCanvas()->Update();
+               // Get the spectrum
-               // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak)
+               histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
-               if(adcpedestal[1] == -1)
+               npeaks = 15;
-               {
+               double par[300];
-                  adcpedestal[0] = j+meanparam+paramsigma;
+               spec = new TSpectrum(npeaks);
+               // Find spectrum background
+               histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
+               // Clone histogram and subtract background from it if we select that option
+               h2 = (TH1F*)histtemp->Clone("h2");
+               if(fitChecks->widgetChBox[0]->IsDown())
-                  adcpedestal[1] = yp;
+                  h2->Add(histback, -1);
-               }
+               // Search for the peaks
+               int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
+               printf("PhotonMu(): Found %d candidates to fit.\n",found);
-               else
+               npeaks = found;
+               // Set initial peak parameters
+               xpeaks = spec->GetPositionX();
+               for(j = 0; j < found; j++)
+               {
+                  float xp = xpeaks[j];
+                  int bin = h2->GetXaxis()->FindBin(xp);
+                  float yp = h2->GetBinContent(bin);
+                  par[3*j] = yp;
+                  par[3*j+1] = xp;
+                  par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
+               }
+               // Fit the histogram
+               fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
+               TVirtualFitter::Fitter(histtemp, 3*npeaks);
+               fit->SetParameters(par);
+               fit->SetNpx(300);
+               h2->Fit("fit","Q");
+               // Get the fitted parameters
+               fittingfunc = h2->GetFunction("fit");
+               fparam = fittingfunc->GetParameters();
+               fparamerr = fittingfunc->GetParErrors();
+               // Gather the parameters (mean peak value for now)
+               int j = 1;
+               int nrfit = 0;
+               while(1)
+               {
-                  if( (npeaks > 1) && (adcpedestal[1] >= yp) )
+                  if( (fparam[j] < 1.E-30) || (nrfit > 8) )
+                     break;
+                  else
+                  {
+                     // Check if pedestal is above the lower limit and sigma is smaller than the mean
+                     if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+                     {
+                        // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
+                        meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
-                     adcpedestal[0] = j+meanparam+paramsigma;
+                        sigmasel[nrfit] = fparam[j+1];
-                     adcpedestal[1] = yp;
+                        nrfit++;
+                     }
+                  }
+                  j+=3;
+               }
+               TMath::Sort(nrfit, meansel, sortindex, kFALSE);
+               meanparam = meansel[sortindex[0]];
+               paramsigma = sigmasel[sortindex[0]];
+               for(j = 0; j < nrfit; j++)
+                  if(DBGSIG)
+                     printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
+               j = 0;
+               adcpedestal[0] = 0;
+               adcpedestal[1] = -1;
+               while(1)
+               {
+                  int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
+                  int yp = histtemp->GetBinContent(bin);
-                  else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
+                  // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak)
+                  if(adcpedestal[1] == -1)
+                  {
                      adcpedestal[0] = j+meanparam+paramsigma;
                      adcpedestal[1] = yp;
+                  }
                   else
+                  {
+                     if( (npeaks > 1) && (adcpedestal[1] >= yp) )
+                     {
+                        adcpedestal[0] = j+meanparam+paramsigma;
+                        adcpedestal[1] = yp;
+                     }
+                     else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
+                     {
+                        adcpedestal[0] = j+meanparam+paramsigma;
+                        adcpedestal[1] = yp;
+                     }
+                     else
-                     break;
+                        break;
+                  }
+                  j++;
+                  if(j > 50) break;
+               }
-               j++;
-               if(j > 50) break;
-            }
-            if( (npeaks > 1) && (nrfit > 1) )
-            {
-               int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
-               adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
-               printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
-               adcpedestal[1] = histtemp->GetBinContent(bin);
-            }
-            if(midPeak->widgetChBox[0]->IsDown())
-            {
-               if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
-                  m = TMath::Floor(meanparam);
-               else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
-                  m = TMath::Ceil(meanparam);
-               int bin = histtemp->GetXaxis()->FindBin(m);
-               adcpedestal[0] = m;
-               printf("midpeak x = %d, ", adcpedestal[0]);
-               adcpedestal[1] = histtemp->GetBinContent(bin);
-            }
+               if(npeaks > 1)
+               {
+                  int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
+                  adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
+                  printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
+                  adcpedestal[1] = histtemp->GetBinContent(bin);
+               }
+               if(midPeak->widgetChBox[0]->IsDown())
+               {
+                  if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
+                     m = TMath::Floor(meanparam);
+                  else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
+                     m = TMath::Ceil(meanparam);
+                  int bin = histtemp->GetXaxis()->FindBin(m);
+                  adcpedestal[0] = m;
+                  printf("midpeak x = %d, ", adcpedestal[0]);
+                  adcpedestal[1] = histtemp->GetBinContent(bin);
+               }
-/*          // Option to show the fit
+/*             // Option to show the fit
-            fittingfunc->Draw("L SAME");
+               fittingfunc->Draw("L SAME");
-            analysisCanvas->GetCanvas()->Modified();
+               analysisCanvas->GetCanvas()->Modified();
-            analysisCanvas->GetCanvas()->Update();*/
+               analysisCanvas->GetCanvas()->Update();*/
-            printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
+               printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
-            // Delete the opened histogram and spectrum
+               // Delete the opened histogram and spectrum
-            delete spec;
+               delete spec;
-            delete inroot;
+               delete inroot;
+            }
             // Open the input file and read header, ADC and TDC values
             sprintf(ctemp, "%s", files->At(i)->GetTitle());
             inroot = new TFile(ctemp, "READ");
             integralPedestal[i] += (double)m2;
             printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]);
             if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) )
                zeromu = i;
-            // Checking for errors when fitting a histogram
-            if(k2 == 0)
-            {
-               printf("PhotonMu(): No pedestal entries found. Check the fitting results.\n");
-               muval[i] = -1;
-            }
-            else
-               muval[i] = -TMath::Log((double)k2/(double)k);
+            muval[i] = -TMath::Log((double)k2/(double)k);
             printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]);
             inroot->Close();
             delete inroot;
+         }
          // Update the progress bar
          progVal = (float)(90.00/nrfiles)*i;
          analysisProgress->widgetPB->SetPosition(progVal);
          gVirtualX->Update(1);
+      }
       printf("PhotonMu(): %d files were selected.\n", nrfiles);
       printf("PhotonMu(): angle\tmu\trelative PDE\n");
       m = 0;
       // Set the 0 degree muval, reuse meansel[1]
       meansel[1] = muval[zeromu];
       printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]);
       // TODO - point estimation still not working correctly!
             exclude = false;
             // Get next point values (if zero value -> need to add the dark hist value again)
             pointest[10] = angle[i];
             pointest[11] = muval[i];
             // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0]
             meansel[0] = PointEstimate(5, pointest);    // PointEstimate only works with very small step size
             if(meansel[0] > accError->widgetNE[0]->GetNumber())
+            {
                printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]);
                exclude = true;
+            }
             // Value with 0 angle and dark histogram are always needed, so should not be excluded
             if(i == darkhist)
                exclude = false;
-            // Wrong fit
-            if(muval[i] == -1)
-               exclude = true;
             // If nothing excluded, pass the points in pointest variable like in a FIFO
             if(!exclude)
+            {
                for(int j = 0; j < 10; j++)
                   // Check if pedestal is above the lower limit and sigma is smaller than the mean
                   if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
+                  {
                      // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
                      meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
+                     meanselerr[nrfit] = fparamerr[j];
                      sigmasel[nrfit] = fparam[j+1];
                      nrfit++;
+                  }
+               }
             else
+            {
                printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n");
                exclude = false;
+            }
+            printf("BreakdownVolt(): Calculated separation between peaks %d and %d: sep = %lf, seperr = %lf\n", (int)peakSepCalc->widgetNE[0]->GetNumber()-1, (int)peakSepCalc->widgetNE[0]->GetNumber(), sep[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p], seperr[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p]);
             // Write out parameters to a file
             fp = fopen(paramname, "a");
             if(exclude)
+            {
       k = peakSepCalc->widgetNE[0]->GetNumber();
       if(k < 4)
          bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]);
       else
+      {
-         printf("BreakdownVold(): Unsupported peak separation selected (%d).\n", k);
+         printf("BreakdownVolt(): Unsupported peak separation selected (%d).\n", k);
          return;
+      }
       bdplot->SetMarkerStyle(20);
       bdplot->SetMarkerSize(0.4);
       if(!cleanPlots)
+      {
          sprintf(ctemp, "#splitline{#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf}", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
          latex = new TLatex();
          latex->SetTextSize(0.039);
-         latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[p-1]], ctemp);
+         latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[1]], ctemp);
          printf("#Delta_{p}(U) = (%.6lf #pm %.8lf)#timesU + (%.6lf #pm %.8lf)}{U_{0} = %.6lf #pm %.8lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
+      }
       else
          printf("#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
          gCanvas = tempAnalysisCanvas->GetCanvas();
       double range[4];
       TGraph2D *gScan2D;
       gScan2D = new TGraph2D();
-      gScan2D->SetName("surfscan");
       range[0] = TMath::MinElement(nrfiles, surfx);
       range[1] = TMath::MaxElement(nrfiles, surfx);
       range[2] = TMath::MinElement(nrfiles, surfy);
       range[3] = TMath::MaxElement(nrfiles, surfy);
+      }
       else if(edit == 1)
+      {
          gCanvas->Modified();
          gCanvas->Update();
-         UpdateIntegrateSurface(-1);
+      }
+   }
+}

Subversion Repositories f9daq

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