Crv Calibration Updates

CRVConditions/data/status_extracted_20260214.txt

@@ -0,0 +1,10 @@
+TABLE CRVBadChan
+#A channel can have more than one status bit. Add all status bits for the complete status.
+#status 1 (bit 0): not connected
+#status 2 (bit 1): ignore channel in reconstruction
+#status 4 (bit 2): no data
+#status 8 (bit 3): no pedestal
+#status 16 (bit 4): no calibration constant
+#status 32 (bit 5): noisy
+#channel,status
+1107,3

CRVConditions/inc/CRVOrdinal.hh

@@ -32,6 +32,13 @@ class CRVOrdinal : virtual public ProditionsEntity {
   CRVOrdinal(OnlineMap const& onMap, OfflineMap const& offMap) :
       ProditionsEntity(cxname), _onMap(onMap), _offMap(offMap) {}
 
+  // check, if a channel exists
+  bool onlineExists(std::uint16_t channel) const {
+    if (channel >= _onMap.size()) return false;
+    if (_onMap.at(channel).FEBchannel() >= CRVId::nChanPerFEB) return false;
+    return true;
+  }
+
   // online numbering triplet for an offline channel number
   const CRVROC& online(std::uint16_t channel) const {
     if (_onMap.at(channel).FEBchannel() >= CRVId::nChanPerFEB) {

CRVReco/src/CrvCalibration_module.cc

@@ -46,11 +46,16 @@ namespace mu2e
       fhicl::Atom<double>      histMaxPulseHeight{Name("histMaxPulseHeight"), Comment("end range of pulseArea histogram"), 150.0};
       fhicl::Atom<double>      fitRangeStart{Name("fitRangeStart"), Comment("low end of the 1PE fit range as fraction of peak"), 0.8};
       fhicl::Atom<double>      fitRangeEnd{Name("fitRangeEnd"), Comment("high end of the 1PE fit range as fraction of peak"), 1.2};
+      fhicl::Atom<double>      minPeakPulseArea{Name("minPeakPulseArea"), Comment("minimum accepted SPE peak for pulseArea histogram"), 250.0};
+      fhicl::Atom<double>      minPeakPulseHeight{Name("minPeakPulseHeight"), Comment("minimum accepted SPE peak for pulseHeight histogram"), 10.0};
       fhicl::Atom<int>         minHistEntries{Name("minHistEntries"), Comment("minimum number of entries required for a fit"), 100};
-      fhicl::Atom<int>         spectrumNPeaks{Name("spectrumNPeaks"), Comment("maximum number of peaks searched by TSpectrum. numbers less then 4 result in warnings"), 6};
+      fhicl::Atom<int>         spectrumNPeaks{Name("spectrumNPeaks"), Comment("maximum number of peaks searched by TSpectrum"), 100};
       fhicl::Atom<double>      spectrumPeakSigma{Name("spectrumPeakSigma"), Comment("TSpectrum search parameter sigma"), 4.0};
       fhicl::Atom<double>      spectrumPeakThreshold{Name("spectrumPeakThreshold"), Comment("TSpectrum search parameter threshold"), 0.01};
-      fhicl::Atom<double>      peakRatioTolerance{Name("peakRatioTolerance"), Comment("allowed deviation of the ratio between 1PE peak and 2PE peak from 2.0"), 0.1};
+      fhicl::Atom<double>      maxFitDifferencePulseArea{Name("maxFitDifferencePulseArea"),
+                                   Comment("maximum difference between the TSpectrum peak and the fitted peak (for pulse areas). Indicates fit problems."), 100.0};
+      fhicl::Atom<double>      maxFitDifferencePulseHeight{Name("maxFitDifferencePulseHeight"),
+                                   Comment("maximum difference between the TSpectrum peak and the fitted peak (for pulse heights). Indicates fit problems."), 4.0};
       fhicl::Atom<std::string> tmpDBfileName{Name("tmpDBfileName"), Comment("name of the tmp. DB file name for the pedestals")};
     };
 
@@ -60,18 +65,19 @@ namespace mu2e
     void analyze(const art::Event& e);
     void beginRun(const art::Run&);
     void endJob();
-    bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak);
+    bool FindSPEpeak(TH1F *hist, TSpectrum &spectrum, TF1 &function, double &SPEpeak, double minPeak, double maxFitDifference);
 
     private:
     std::string        _crvRecoPulsesModuleLabel;
     int                _histBinsPulseArea, _histBinsPulseHeight;
     double             _histMaxPulseArea, _histMaxPulseHeight;
     double             _fitRangeStart, _fitRangeEnd;
     int                _minHistEntries;
+    double             _minPeakPulseArea, _minPeakPulseHeight;
     int                _spectrumNPeaks;
     double             _spectrumPeakSigma;
     double             _spectrumPeakThreshold;
-    double             _peakRatioTolerance;
+    double             _maxFitDifferencePulseArea, _maxFitDifferencePulseHeight;
     std::string        _tmpDBfileName;
     std::vector<TH1F*> _calibHistsPulseArea;
     std::vector<TH1F*> _calibHistsPulseHeight;
@@ -80,9 +86,6 @@ namespace mu2e
 
     std::vector<double> _pedestals;
     std::vector<double> _timeOffsets;
-
-    std::pair<int,int>  _firstRunSubrun;
-    std::pair<int,int>  _lastRunSubrun;
   };
 
 
@@ -96,10 +99,13 @@ 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()),
-    _peakRatioTolerance(conf().peakRatioTolerance()),
+    _maxFitDifferencePulseArea(conf().maxFitDifferencePulseArea()),
+    _maxFitDifferencePulseHeight(conf().maxFitDifferencePulseHeight()),
     _tmpDBfileName(conf().tmpDBfileName())
   {
   }
@@ -143,39 +149,27 @@ namespace mu2e
     treePedestal->Branch("channel", &channel);
     treePedestal->Branch("pedestal", &pedestal);
 
-    TF1 funcCalib("SPEpeak", "gaus");
-    TSpectrum spectrum(_spectrumNPeaks); //any value of 3 or less results in a "Peak buffer full" warning.
+    TF1 function("calibPeak","gaus");
+    TSpectrum spectrum(_spectrumNPeaks);
 
     std::ofstream outputFile;
     outputFile.open(_tmpDBfileName);
-    outputFile<<"TABLE CRVSiPM "<<_firstRunSubrun.first<<":"<<_firstRunSubrun.second<<"-"<<_lastRunSubrun.first<<":"<<_lastRunSubrun.second<<std::endl;
+    if(!outputFile.is_open()) throw cet::exception("CRVCALIB") << "Couldn't create output txt file " << _tmpDBfileName << ".";
+    outputFile<<"TABLE CRVSiPM"<<std::endl;
     outputFile<<"#channel, pedestal, calibPulseHeight, calibPulseArea"<<std::endl;
 
     for(channel=0; channel<_pedestals.size(); ++channel)
     {
       TH1F *hist;
-      double calibValue[2];
+      double calibValue[2]={-1,-1};
       for(int i=0; i<2; ++i) //loop over hisograms with pulse areas and pulse heights
       {
         if(i==1) hist=_calibHistsPulseArea.at(channel);
         else hist=_calibHistsPulseHeight.at(channel);
 
-        double peakCalib=0;
-        if(!FindSPEpeak(hist, spectrum, peakCalib))
-        {
-          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);
+        double SPEpeak=-1;
+        if(!FindSPEpeak(hist, spectrum, function, SPEpeak, (i==0?_minPeakPulseHeight:_minPeakPulseArea), (i==0?_maxFitDifferencePulseHeight:_maxFitDifferencePulseArea))) continue;
+        calibValue[i]=SPEpeak;
       }
 
       pedestal=_pedestals.at(channel);
@@ -220,10 +214,7 @@ namespace mu2e
         _pedestals[channelIndex] = calib.pedestal(channelIndex);
         _timeOffsets[channelIndex] = calib.timeOffset(channelIndex);
       }
-
-      _firstRunSubrun=std::pair<int,int>(event.run(),event.subRun());
     }
-    _lastRunSubrun=std::pair<int,int>(event.run(),event.subRun());
 
     for(auto iter=crvRecoPulseCollection->begin(); iter!=crvRecoPulseCollection->end(); ++iter)
     {
@@ -241,29 +232,36 @@ namespace mu2e
     }
   }
 
-  bool CrvCalibration::FindSPEpeak(TH1F *hist, TSpectrum &spectrum, double &SPEpeak)
+  bool CrvCalibration::FindSPEpeak(TH1F *hist, TSpectrum &spectrum, TF1 &function, double &SPEpeak, double minPeak, double maxFitDifference)
   {
     if(hist->GetEntries()<_minHistEntries) return false; //not enough data
 
     int nPeaks = spectrum.Search(hist,_spectrumPeakSigma,"nodraw",_spectrumPeakThreshold);
-    if(nPeaks==0) return false;
-
-    //peaks are not returned sorted
+    if(nPeaks<=0) return false;
+
+    //peaks are returned sorted by Y
+    //from our long-time experience:
+    //-SPE peak is either the highest peak or second highest peak
+    //-if the SPE peak is the second highest, then the highest peak comes from the baseline
+    //-the peak from the baseline is always below the minPeak threshold, while the SPE peak is not
+    //-the minPeak threshold may have to be adjusted for non-standard bias voltages
     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;});
-
-    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
-                                       //never seen peaks at 0, but still checking to avoid division by 0.
+    double x=peaksX[0];
+    if(x<minPeak)
     {
-      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(nPeaks==1) return false;
+      x=peaksX[1];
+      if(x<minPeak) return false;
     }
-    SPEpeak = peaks[peakToUse].first;
+
+    if(hist->FindBin(x*_fitRangeStart)==hist->FindBin(x*_fitRangeEnd)) return false; //fit range start/end are in the same bin
+    function.SetRange(x*_fitRangeStart,x*_fitRangeEnd);
+    function.SetParameter(1,x);
+    hist->Fit(&function, "QR");
+    SPEpeak = function.GetParameter(1);
+
+    if(fabs(SPEpeak-x)>maxFitDifference) return false;
+    if(SPEpeak<minPeak) return false;
     return true;
   }
 

CRVReco/src/CrvDQMcollector_module.cc

@@ -411,10 +411,10 @@ namespace mu2e
                                             Form("crvDigiRatesNZS_ROC%zu",ROC),
                                             CRVId::nFEBPerROC*CRVId::nChanPerFEB,0,CRVId::nFEBPerROC*CRVId::nChanPerFEB));
     }
-    _hist2DPEsMPVROC=tfs->make<TH2F>("crvPEsMPV","crvPEsMPV", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC);
-    _hist2DDigiRatesROC=tfs->make<TH2F>("crvDigiRates","crvDigiRates", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC);
-    _hist2DDigiRatesROCNZS=tfs->make<TH2F>("crvDigiRatesNZS","crvDigiRatesNZS", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC);
-    _histCoincidenceClusters=tfs->make<TH1I>("crvCoincidencesClusters","crvCoincidenceClusters",10,0,10);
+    _hist2DPEsMPVROC=tfs->make<TH2F>("crvPEsMPV","crvPEsMPV:FEBchannel:FEB", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC);
+    _hist2DDigiRatesROC=tfs->make<TH2F>("crvDigiRates","crvDigiRates:FEBchannel:FEB", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC);
+    _hist2DDigiRatesROCNZS=tfs->make<TH2F>("crvDigiRatesNZS","crvDigiRatesNZS:FEBchannel:FEB", CRVId::nChanPerFEB,0,CRVId::nChanPerFEB, CRVId::nROC*CRVId::nFEBPerROC,0,CRVId::nROC*CRVId::nFEBPerROC);
+    _histCoincidenceClusters=tfs->make<TH1I>("crvCoincidencesClusters","crvCoincidenceClusters:sectorType",10,0,10);
 
     _treeMetaData=tfs->make<TTree>("crvMetaData","crvMetaData");
     _treeMetaData->Branch("runNumberStart",&_firstRunSubrun.first);

CRVReco/src/CrvRecoPulsesFinder_module.cc

@@ -115,6 +115,10 @@ namespace mu2e
     _useTimeOffsetDB(conf().useTimeOffsetDB()),
     _ignoreChannels(conf().ignoreChannels())
   {
+    if(conf().pulseAreaThreshold()>conf().minADCdifference())
+      throw cet::exception("CRVRECO_BAD_CONFIG")
+      << "CrvRecoPulseFinder pulseAreaThreshold " << conf().pulseAreaThreshold()
+      << "  larger than minADCdifference " << conf().minADCdifference();
     produces<CrvRecoPulseCollection>(_NZSdata?"NZS":"");
     _makeCrvRecoPulses=boost::shared_ptr<mu2eCrv::MakeCrvRecoPulses>(new mu2eCrv::MakeCrvRecoPulses(conf().minADCdifference(),
                                                                                                     conf().defaultBeta(),

CRVReco/src/CrvTimingStudies_module.cc

@@ -10,6 +10,7 @@
 #include "Offline/DataProducts/inc/CRVId.hh"
 #include "Offline/ProditionsService/inc/ProditionsHandle.hh"
 #include "Offline/CRVConditions/inc/CRVOrdinal.hh"
+#include "Offline/CRVConditions/inc/CRVCalib.hh"
 
 #include "Offline/GeometryService/inc/DetectorSystem.hh"
 #include "Offline/GeometryService/inc/GeomHandle.hh"
@@ -49,6 +50,7 @@ namespace mu2e
     {
       fhicl::Atom<std::string> crvRecoPulsesModuleLabel{ Name("crvRecoPulseModuleLabel"), Comment("CrvRecoPulse Label")};
       fhicl::Atom<double>      PEthreshold{ Name("PEthreshold"), Comment("PE threshold")};
+      fhicl::Atom<bool>        removeTimeOffsets{ Name("removeTimeOffsets"), Comment("remove time offsets added by reco")};
     };
 
     using Parameters = art::EDAnalyzer::Table<Config>;
@@ -61,16 +63,21 @@ namespace mu2e
     private:
     std::string   _crvRecoPulsesModuleLabel;
     double        _PEthreshold;
+    bool          _removeTimeOffsets;
+    bool          _firstEvent;
 
     std::map<std::pair<int,int>,TH1F*>  _histTimeDiffs;
 
     ProditionsHandle<CRVOrdinal> _crvChannelMap_h;
+    ProditionsHandle<CRVCalib>   _calib_h;
   };
 
   CrvTimingStudies::CrvTimingStudies(const Parameters& conf) :
     art::EDAnalyzer{conf},
     _crvRecoPulsesModuleLabel(conf().crvRecoPulsesModuleLabel()),
-    _PEthreshold(conf().PEthreshold())
+    _PEthreshold(conf().PEthreshold()),
+    _removeTimeOffsets(conf().removeTimeOffsets()),
+    _firstEvent(true)
   {
   }
 
@@ -92,6 +99,53 @@ namespace mu2e
     if(!event.getByLabel(_crvRecoPulsesModuleLabel,"",crvRecoPulseCollection)) return;
 
     auto const& crvChannelMap = _crvChannelMap_h.get(event.id());
+    auto const& calib = _calib_h.get(event.id());
+
+    if(_firstEvent)
+    {
+      _firstEvent=false;
+
+      //store channel map, pedestals and calibration constants in the file,
+      //so that it can later be used to write a full calibration set
+      //of pedestals, calib consts, and time offsets
+
+      GeomHandle<CosmicRayShield> CRS;
+      const std::vector<std::shared_ptr<CRSScintillatorBar> > &counters = CRS->getAllCRSScintillatorBars();
+      art::ServiceHandle<art::TFileService> tfs;
+
+      TTree *treeChannelMap = tfs->make<TTree>("channelMap","channelMap");
+      size_t channel;
+      int    roc, feb, febChannel;
+      treeChannelMap->Branch("channel", &channel);
+      treeChannelMap->Branch("roc", &roc);
+      treeChannelMap->Branch("feb", &feb);
+      treeChannelMap->Branch("febChannel", &febChannel);
+
+      for(channel=0; channel<counters.size()*CRVId::nChanPerBar; ++channel)
+      {
+        if(!crvChannelMap.onlineExists(channel)) continue;
+        CRVROC onlineChannel  = crvChannelMap.online(channel);
+        roc         = onlineChannel.ROC();
+        feb         = onlineChannel.FEB();
+        febChannel  = onlineChannel.FEBchannel();
+        treeChannelMap->Fill();
+      }
+
+      TTree *treeCalib = tfs->make<TTree>("crvCalib","crvCalib");
+      double pedestal, calibPulseHeight, calibPulseArea;
+      treeCalib->Branch("channel", &channel);
+      treeCalib->Branch("pedestal", &pedestal);
+      treeCalib->Branch("calibPulseHeight", &calibPulseHeight);
+      treeCalib->Branch("calibPulseArea", &calibPulseArea);
+
+      for(channel=0; channel<counters.size()*CRVId::nChanPerBar; ++channel)
+      {
+        pedestal = calib.pedestal(channel);
+        calibPulseHeight = calib.pulseHeight(channel);
+        calibPulseArea = calib.pulseArea(channel);
+        treeCalib->Fill();
+      }
+    }
 
     int previousOfflineChannel=-1;
     std::map<uint16_t,std::vector<double> > fpgaTimes;
@@ -111,11 +165,14 @@ 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);
+
+      if(_removeTimeOffsets) //remove time offsets introduced during reconstruction to see the "pure" time differences, e.g. to generate new time calibration tables
+      {
+        double timeOffset = calib.timeOffset(offlineChannel);
+        recoPulseTime-=timeOffset;
+      }
 
       fpgaTimes[fpgaIndex].push_back(recoPulseTime);
     }
@@ -141,11 +198,8 @@ feb--;
       {
         art::ServiceHandle<art::TFileService> tfs;
         _histTimeDiffs[histIndex] = tfs->make<TH1F>(Form("fpgaTimeDiff_%i_%i",fpga1->first,fpga2->first),
-                                                    Form("Time Diffs between FGPAs %i and %i;time difference [ns];Counts",fpga1->first,fpga2->first),
-                                                    100,-50,50);
-//        _histTimeDiffs[histIndex] = tfs->make<TH1F>(Form("AFETimeDiff_%i_%i",fpga1->first,fpga2->first),
-//                                                    Form("Time Diffs between AFEs %i and %i;time difference [ns];Counts",fpga1->first,fpga2->first),
-//                                                    100,-50,50);
+                                                    Form("Time Diffs between FPGAs %i and %i;time difference [ns];Counts",fpga1->first,fpga2->first),
+                                                    300,-150,150);
       }
 
       double timeDiff=fpga1->second-fpga2->second;

CRVReco/src/MakeCrvRecoPulses.cc

@@ -31,7 +31,12 @@ MakeCrvRecoPulses::MakeCrvRecoPulses(float minADCdifference, float defaultBeta,
                                      _pulseThreshold(pulseThreshold),
                                      _pulseAreaThreshold(pulseAreaThreshold),
                                      _doublePulseSeparation(doublePulseSeparation)
-{}
+{
+  if(_pulseAreaThreshold>_minADCdifference)
+      throw cet::exception("CRVRECO_BAD_CONFIG")
+      << "MakeCrvRecoPulses pulseAreaThreshold " << _pulseAreaThreshold
+      << "  larger than minADCdifference " << _minADCdifference;
+}
 
 void MakeCrvRecoPulses::FillGraphAndFindPeaks(const std::vector<int16_t> &waveform, uint16_t startTDC,
                                               float digitizationPeriod, float pedestal,

EventDisplay/fcl/EventDisplayExtractedCrvOnly.fcl

@@ -2,6 +2,7 @@
 
 physics.analyzers.eventdisplay.extractedCrvOnly      : true
 services.GeometryService.inputFile                   : "Offline/Mu2eG4/geom/geom_common_extracted.txt"
-services.ProditionsService.crvStatus.useDb           : false
-services.ProditionsService.crvCalib.useDb            : false
-services.ProditionsService.crvCalib.pedestal         : 2047
+services.ProditionsService.crvCalib.useDb            : true
+services.ProditionsService.crvStatus.useDb           : true
+#comment out this line, if displaying simulated events
+services.DbService.textFile                          : ["Offline/CRVConditions/data/status_extracted_20260214.txt","Offline/CRVConditions/data/calib_extracted_20260214.txt"]