updated CRV calibration macro

CRV/CrvCalibration.C

@@ -1,12 +1,15 @@
 const double fitRangeStart=0.8;
 const double fitRangeEnd=1.2;
 const int    minHistEntries=100;
-const int    spectrumNPeaks=6;
+const double minPeakPulseHeight=10.0;
+const double minPeakPulseArea=250.0;
+const int    spectrumNPeaks=100;
 const double spectrumPeakSigma=4.0;
-const double spectrumPeakThreshold=0.01;
-const double peakRatioTolerance=0.2;
+const double spectrumPeakThreshold=0.001;
+const double peakRatioTolerance=0.3;
 
-bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak);
+template<size_t N>
+bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, std::array<TF1,N> &functions, double &SPEpeak, double minPeak);
 
 void CrvCalibration(const std::string &inputFileName, const std::string &outputFileName)
 {
@@ -24,7 +27,7 @@ void CrvCalibration(const std::string &inputFileName, const std::string &outputF
       pedestals[channel]=pedestal;
     }
 
-    TF1 funcCalib("SPEpeak", "gaus");
+    std::array functions={TF1("calibPeak1","gaus"), TF1("calibPeak2","gaus"), TF1("calibPeak3","gaus")}; //only need to fit three peaks
     TSpectrum spectrum(spectrumNPeaks); //any value of 3 or less results in a "Peak buffer full" warning.
 
     std::ofstream outputFile;
@@ -39,27 +42,18 @@ void CrvCalibration(const std::string &inputFileName, const std::string &outputF
 
       TH1F *hist;
       double calibValue[2];
-      for(int i=0; i<2; ++i) //loop over hisograms with pulse areas and pulse heights
+      for(int i=0; i<2; ++i) //loop over histograms with pulse areas and pulse heights
       {
         if(i==1) hist=(TH1F*)gDirectory->FindObjectAny(Form("crvCalibrationHistPulseArea_%zu",channel));
         else hist=(TH1F*)gDirectory->FindObjectAny(Form("crvCalibrationHistPulseHeight_%zu",channel));
 
-        double peakCalib=0;
-        if(!FindSPEpeak(hist, spectrum, peakCalib))
+        double SPEpeak=-1;
+        if(!FindSPEpeak(hist, spectrum, functions, SPEpeak, (i==0?minPeakPulseHeight:minPeakPulseArea)))
         {
           calibValue[i]=-1;
           continue;
         }
-
-        funcCalib.SetRange(peakCalib*fitRangeStart,peakCalib*fitRangeEnd);
-        if(hist->FindBin(peakCalib*fitRangeStart)==hist->FindBin(peakCalib*fitRangeEnd)) //fit range start/end are in the same bin
-        {
-          calibValue[i]=-1;
-          continue;
-        }
-        funcCalib.SetParameter(1,peakCalib);
-        hist->Fit(&funcCalib, "QR");
-        calibValue[i]=funcCalib.GetParameter(1);
+        calibValue[i]=SPEpeak;
       }
 
       outputFile<<channel<<","<<pedestal<<","<<calibValue[0]<<","<<calibValue[1]<<std::endl;
@@ -91,27 +85,40 @@ void CrvCalibration(const std::string &inputFileName, const std::string &outputF
     inputFile->Close();
 }
 
-bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak)
+template<size_t N>
+bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, std::array<TF1,N> &functions, double &SPEpeak, double minPeak)
 {
     if(hist->GetEntries()<minHistEntries) return false; //not enough data
 
-    int nPeaks = spectrum.Search(hist,spectrumPeakSigma,"nodraw",spectrumPeakThreshold);
+    size_t nPeaks = spectrum.Search(hist,spectrumPeakSigma,"nodraw",spectrumPeakThreshold);
     if(nPeaks==0) return false;
 
-    //peaks are not returned sorted
+    //peaks are returned sorted by Y
     double *peaksX = spectrum.GetPositionX();
-    double *peaksY = spectrum.GetPositionY();
-    std::vector<std::pair<double,double> > peaks;
-    for(int iPeak=0; iPeak<nPeaks; ++iPeak) peaks.emplace_back(peaksX[iPeak],peaksY[iPeak]);
-    std::sort(peaks.begin(),peaks.end(), [](const std::pair<double,double> &a, const std::pair<double,double> &b) {return a.first<b.first;});
+    std::vector<double> fittedPeaks;
+    for(size_t iPeak=0; iPeak<nPeaks && iPeak<functions.size(); ++iPeak)
+    {
+      double x=peaksX[iPeak];
+      if(hist->FindBin(x*fitRangeStart)==hist->FindBin(x*fitRangeEnd)) continue; //fit range start/end are in the same bin
+      functions[iPeak].SetRange(x*fitRangeStart,x*fitRangeEnd);
+      functions[iPeak].SetParameter(1,x);
+      hist->Fit(&functions[iPeak], "QR+");
+      fittedPeaks.emplace_back(functions[iPeak].GetParameter(1));
+    }
+    if(fittedPeaks.size()==0) return false;
 
-    int peakToUse=0;
-    if(nPeaks>1 && peaks[0].first>0)   //if more than one peak is found, the first peak could be due to baseline fluctuations
+    int peakToUse=-1;
+    //only need to test two highest peaks (=first two entries in vector)
+    //one of the peaks could be due to baseline fluctuations
+    for(size_t iPeak=0; iPeak<fittedPeaks.size() && iPeak<2; ++iPeak)
+    for(size_t jPeak=iPeak+1; jPeak<fittedPeaks.size(); ++jPeak)
     {
-      if(fabs(peaks[1].first/peaks[0].first-2.0)>peakRatioTolerance) peakToUse=1; //2nd peak is not twice the 1st peak, so the 1st peak is not the SPE peak
-                                                                                  //assume that the 2nd peak is the SPE peak
-                                                                                  //we have never seen that the 3rd peak was the SPE peak - no need to test it
+      if(fabs(fittedPeaks.at(jPeak)/fittedPeaks.at(iPeak)-2.0)<peakRatioTolerance) peakToUse=iPeak;
     }
-    SPEpeak = peaks[peakToUse].first;
+    if(peakToUse==-1) return false;
+    if(fittedPeaks.at(peakToUse)<minPeak) return false;
+
+    SPEpeak = fittedPeaks.at(peakToUse);
+
     return true;
 }