diff --git a/CRV/CrvCalibration.C b/CRV/CrvCalibration.C index 66a0fe0..50ab44f 100644 --- a/CRV/CrvCalibration.C +++ b/CRV/CrvCalibration.C @@ -1,3 +1,13 @@ +const double fitRangeStart=0.8; +const double fitRangeEnd=1.2; +const int minHistEntries=100; +const int spectrumNPeaks=6; +const double spectrumPeakSigma=4.0; +const double spectrumPeakThreshold=0.01; +const double peakRatioTolerance=0.1; + +bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak); + void CrvCalibration(const std::string &inputFileName, const std::string &outputFileName) { TFile *inputFile = TFile::Open(inputFileName.c_str(),"update"); @@ -14,7 +24,8 @@ void CrvCalibration(const std::string &inputFileName, const std::string &outputF pedestals[channel]=pedestal; } - TF1 funcCalib("f0", "gaus"); + TF1 funcCalib("SPEpeak", "gaus"); + TSpectrum spectrum(spectrumNPeaks); //any value of 3 or less results in a "Peak buffer full" warning. std::ofstream outputFile; outputFile.open(outputFileName); @@ -33,33 +44,21 @@ void CrvCalibration(const std::string &inputFileName, const std::string &outputF if(i==1) hist=(TH1F*)gDirectory->FindObjectAny(Form("crvCalibrationHistPulseArea_%zu",channel)); else hist=(TH1F*)gDirectory->FindObjectAny(Form("crvCalibrationHistPulseHeight_%zu",channel)); - if(hist->GetEntries()<100) //not enough data - { - calibValue[i]=-1; - continue; - } - -/* - int n=hist->GetNbinsX(); - double overflow=hist->GetBinContent(0)+hist->GetBinContent(n+1); - if(overflow/((double)hist->GetEntries())>0.1) //too much underflow/overflow. something may be wrong. + double peakCalib=0; + if(!FindSPEpeak(hist, spectrum, peakCalib)) { calibValue[i]=-1; continue; } -*/ - int maxbinCalib = hist->GetMaximumBin(); - double peakCalib = hist->GetBinCenter(maxbinCalib); -//FIXME funcCalib.SetRange(peakCalib*0.8,peakCalib*1.2); - funcCalib.SetRange(peakCalib*0.7,peakCalib*1.3); - if(hist->FindBin(peakCalib*0.7)==hist->FindBin(peakCalib*1.3)) + 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, "0QR"); + hist->Fit(&funcCalib, "QR"); calibValue[i]=funcCalib.GetParameter(1); } @@ -91,3 +90,28 @@ void CrvCalibration(const std::string &inputFileName, const std::string &outputF outputFile.close(); inputFile->Close(); } + +bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak) +{ + if(hist->GetEntries() > peaks; + for(int iPeak=0; iPeak &a, const std::pair &b) {return a.first1 && peaks[0].first>0) //if more than one peak is found, the first peak could be due to baseline fluctuations + { + 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 + } + SPEpeak = peaks[peakToUse].first; + return true; +} diff --git a/CRV/CrvPass1_extracted.fcl b/CRV/CrvPass1_extracted.fcl index e56df3a..373dec4 100644 --- a/CRV/CrvPass1_extracted.fcl +++ b/CRV/CrvPass1_extracted.fcl @@ -84,6 +84,7 @@ physics.producers.CrvCoincidenceClusterFinder.sectorConfig: maxSlopeDifference : 4 coincidenceLayers : 3 minClusterPEs : 0 + initialClusterMaxDistance : 250 }, { CRVSector : "T1" @@ -97,6 +98,7 @@ physics.producers.CrvCoincidenceClusterFinder.sectorConfig: maxSlopeDifference : 4 coincidenceLayers : 3 minClusterPEs : 0 + initialClusterMaxDistance : 250 }, { CRVSector : "T2" @@ -110,6 +112,7 @@ physics.producers.CrvCoincidenceClusterFinder.sectorConfig: maxSlopeDifference : 4 coincidenceLayers : 3 minClusterPEs : 0 + initialClusterMaxDistance : 250 } ]