diff --git a/CRVReco/src/CrvCalibration_module.cc b/CRVReco/src/CrvCalibration_module.cc index e0408befb3..6101e6f641 100644 --- a/CRVReco/src/CrvCalibration_module.cc +++ b/CRVReco/src/CrvCalibration_module.cc @@ -46,11 +46,13 @@ namespace mu2e fhicl::Atom histMaxPulseHeight{Name("histMaxPulseHeight"), Comment("end range of pulseArea histogram"), 150.0}; fhicl::Atom fitRangeStart{Name("fitRangeStart"), Comment("low end of the 1PE fit range as fraction of peak"), 0.8}; fhicl::Atom fitRangeEnd{Name("fitRangeEnd"), Comment("high end of the 1PE fit range as fraction of peak"), 1.2}; + fhicl::Atom minPeakPulseArea{Name("minPeakPulseArea"), Comment("minimum accepted SPE peak for pulseArea histogram"), 250.0}; + fhicl::Atom minPeakPulseHeight{Name("minPeakPulseHeight"), Comment("minimum accepted SPE peak for pulseHeight histogram"), 10.0}; fhicl::Atom minHistEntries{Name("minHistEntries"), Comment("minimum number of entries required for a fit"), 100}; - fhicl::Atom spectrumNPeaks{Name("spectrumNPeaks"), Comment("maximum number of peaks searched by TSpectrum. numbers less then 4 result in warnings"), 6}; + fhicl::Atom spectrumNPeaks{Name("spectrumNPeaks"), Comment("maximum number of peaks searched by TSpectrum. numbers less then 4 result in warnings"), 100}; fhicl::Atom spectrumPeakSigma{Name("spectrumPeakSigma"), Comment("TSpectrum search parameter sigma"), 4.0}; fhicl::Atom spectrumPeakThreshold{Name("spectrumPeakThreshold"), Comment("TSpectrum search parameter threshold"), 0.01}; - fhicl::Atom peakRatioTolerance{Name("peakRatioTolerance"), Comment("allowed deviation of the ratio between 1PE peak and 2PE peak from 2.0"), 0.1}; + fhicl::Atom peakRatioTolerance{Name("peakRatioTolerance"), Comment("allowed deviation of the ratio between 1PE peak and 2PE peak from 2.0"), 0.2}; fhicl::Atom tmpDBfileName{Name("tmpDBfileName"), Comment("name of the tmp. DB file name for the pedestals")}; }; @@ -60,7 +62,8 @@ namespace mu2e void analyze(const art::Event& e); void beginRun(const art::Run&); void endJob(); - bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak); + template + bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, std::array &functions, double &SPEpeak, double minPeak); private: std::string _crvRecoPulsesModuleLabel; @@ -68,6 +71,7 @@ namespace mu2e double _histMaxPulseArea, _histMaxPulseHeight; double _fitRangeStart, _fitRangeEnd; int _minHistEntries; + double _minPeakPulseArea, _minPeakPulseHeight; int _spectrumNPeaks; double _spectrumPeakSigma; double _spectrumPeakThreshold; @@ -96,6 +100,8 @@ namespace mu2e _fitRangeStart(conf().fitRangeStart()), _fitRangeEnd(conf().fitRangeEnd()), _minHistEntries(conf().minHistEntries()), + _minPeakPulseArea(conf().minPeakPulseArea()), + _minPeakPulseHeight(conf().minPeakPulseHeight()), _spectrumNPeaks(conf().spectrumNPeaks()), _spectrumPeakSigma(conf().spectrumPeakSigma()), _spectrumPeakThreshold(conf().spectrumPeakThreshold()), @@ -143,7 +149,7 @@ namespace mu2e treePedestal->Branch("channel", &channel); treePedestal->Branch("pedestal", &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; @@ -160,22 +166,13 @@ namespace mu2e if(i==1) hist=_calibHistsPulseArea.at(channel); else hist=_calibHistsPulseHeight.at(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, "0QR"); - calibValue[i]=funcCalib.GetParameter(1); + calibValue[i]=SPEpeak; } pedestal=_pedestals.at(channel); @@ -241,29 +238,41 @@ namespace mu2e } } - bool CrvCalibration::FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak) + template + bool CrvCalibration::FindSPEpeak(TH1F *hist, TSpectrum &spectrum, std::array &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 > 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 - //never seen peaks at 0, but still checking to avoid division by 0. + std::vector fittedPeaks; + for(size_t iPeak=0; iPeak_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 + 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)); } - SPEpeak = peaks[peakToUse].first; + if(fittedPeaks.size()==0) return false; + + 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; iPeakmake("crvPEsMPV","crvPEsMPV", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC); - _hist2DDigiRatesROC=tfs->make("crvDigiRates","crvDigiRates", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC); - _hist2DDigiRatesROCNZS=tfs->make("crvDigiRatesNZS","crvDigiRatesNZS", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC); - _histCoincidenceClusters=tfs->make("crvCoincidencesClusters","crvCoincidenceClusters",10,0,10); + _hist2DPEsMPVROC=tfs->make("crvPEsMPV","crvPEsMPV:FEBchannel:FEB", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC); + _hist2DDigiRatesROC=tfs->make("crvDigiRates","crvDigiRates:FEBchannel:FEB", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC); + _hist2DDigiRatesROCNZS=tfs->make("crvDigiRatesNZS","crvDigiRatesNZS:FEBchannel:FEB", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC); + _histCoincidenceClusters=tfs->make("crvCoincidencesClusters","crvCoincidenceClusters:sectorType",10,0,10); _treeMetaData=tfs->make("crvMetaData","crvMetaData"); _treeMetaData->Branch("runNumberStart",&_firstRunSubrun.first); diff --git a/CRVReco/src/CrvTimingStudies_module.cc b/CRVReco/src/CrvTimingStudies_module.cc index 039dc8cc68..13b5338aa2 100644 --- a/CRVReco/src/CrvTimingStudies_module.cc +++ b/CRVReco/src/CrvTimingStudies_module.cc @@ -111,11 +111,8 @@ namespace mu2e uint16_t ROC = onlineChannel.ROC(); uint16_t feb = onlineChannel.FEB(); uint16_t febChannel = onlineChannel.FEBchannel(); -ROC--; -feb--; - uint16_t fpgaIndex = (ROC*CRVId::nFEBPerROC*CRVId::nChanPerFEB+feb*CRVId::nChanPerFEB+febChannel)/(CRVId::nChanPerFEB/CRVId::nFPGAPerFEB); -//uint16_t fpgaIndex = (ROC*CRVId::nFEBPerROC*CRVId::nChanPerFEB+feb*CRVId::nChanPerFEB+febChannel)/(8); -//uint16_t fpgaIndex = (ROC*CRVId::nFEBPerROC*CRVId::nChanPerFEB+feb*CRVId::nChanPerFEB+febChannel); + + uint16_t fpgaIndex = ((ROC-1)*CRVId::nFEBPerROC*CRVId::nChanPerFEB+(feb-1)*CRVId::nChanPerFEB+febChannel)/(CRVId::nChanPerFEB/CRVId::nFPGAPerFEB); fpgaTimes[fpgaIndex].push_back(recoPulseTime); } diff --git a/CRVReco/src/MakeCrvRecoPulses.cc b/CRVReco/src/MakeCrvRecoPulses.cc index 7f5040712f..cc12280c26 100644 --- a/CRVReco/src/MakeCrvRecoPulses.cc +++ b/CRVReco/src/MakeCrvRecoPulses.cc @@ -31,7 +31,9 @@ MakeCrvRecoPulses::MakeCrvRecoPulses(float minADCdifference, float defaultBeta, _pulseThreshold(pulseThreshold), _pulseAreaThreshold(pulseAreaThreshold), _doublePulseSeparation(doublePulseSeparation) -{} +{ + if(_pulseAreaThreshold>_minADCdifference) _pulseAreaThreshold=_minADCdifference; +} void MakeCrvRecoPulses::FillGraphAndFindPeaks(const std::vector &waveform, uint16_t startTDC, float digitizationPeriod, float pedestal,