WebSVN - f9daq - Path Comparison - /lab/sipmscan/trunk/src/analysis.cpp Rev 172 and /lab/sipmscan/trunk/src/analysis.cpp Rev 173

Ignore whitespace Rev 172 → Rev 173

 /lab/sipmscan/trunk/src/analysis.cpp
 ,6 → 59,7
    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);
    }
 ,6 → 97,9
    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)
    {
 ,7 → 374,6
          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);
 ,6 → 682,7
    int adcpedestal[2];
    int zeromu = 0;
    int darkhist = -1;
+   int nopeaks = -1;
    double pointest[12];
    bool exclude = false;
 ,9 → 694,170
    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()) )
 ,145 → 953,148
                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();
+            if(nopeaks == -1)
+            {
+               // 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;
+               // 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();
-            }
+               // 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();
+               // 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()) )
+               // 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
                   {
-                     // 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++;
+                     // 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);
+                  j+=3;
+               }
+               TMath::Sort(nrfit, meansel, sortindex, kFALSE);
-            meanparam = meansel[sortindex[0]];
-            paramsigma = sigmasel[sortindex[0]];
+               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]]);
+               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;
+               j = 0;
+               adcpedestal[0] = 0;
+               adcpedestal[1] = -1;
-            while(1)
-            {
-               int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
-               int yp = histtemp->GetBinContent(bin);
+               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) )
+                  // 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[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
+                  else
                   {
-                     adcpedestal[0] = j+meanparam+paramsigma;
-                     adcpedestal[1] = yp;
+                     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;
                   }
-                  else
-                     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(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);
-            }
+               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();*/
+/*             // 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]);
+               printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
-            // Delete the opened histogram and spectrum
-            delete spec;
-            delete inroot;
+               // 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());
 ,14 → 1162,7
             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();
 ,7 → 1179,7
       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]);
 ,10 → 1209,6
             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)
             {
 ,6 → 1469,7
                   {
                      // 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++;
                   }
 ,6 → 1646,8
                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)
 ,7 → 1700,7
          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;
       }
 ,7 → 1737,7
          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
 ,7 → 1912,6
       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);
 ,8 → 2014,6
       {
          gCanvas->Modified();
          gCanvas->Update();
-         UpdateIntegrateSurface(-1);
       }
    }
 }

Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 172 → Rev 173