CRVResponse: fix defects found by a static sweep of the package

CRVResponse/CMakeLists.txt

@@ -66,6 +66,7 @@ cet_build_plugin(CrvPlot art::module
       Offline::CRVConditions
       Offline::DataProducts
       Offline::GeometryService
+      Offline::GlobalConstantsService
       Offline::MCDataProducts
       Offline::ProditionsService
       Offline::RecoDataProducts

CRVResponse/src/CrvPhotonGenerator_module.cc

@@ -198,7 +198,10 @@ namespace mu2e
       bool tableLoaded=false;
       for(size_t j=0; j<i; ++j)
       {
-        if(_lookupTableFileNames[i]==_lookupTableFileNames[j])
+        // the scintillation yield is state held on the shared MakeCrvPhotons instance, so two
+        // sectors may only share an instance if they agree on the yield as well as the table
+        if(_lookupTableFileNames[i]==_lookupTableFileNames[j] &&
+           _scintillationYields[i]==_scintillationYields[j])
         {
            tableLoaded=true;
            _makeCrvPhotons.emplace_back(_makeCrvPhotons[j]);

CRVResponse/src/CrvPlot_module.cc

@@ -7,6 +7,7 @@
 #include "Offline/CosmicRayShieldGeom/inc/CosmicRayShield.hh"
 #include "Offline/DataProducts/inc/CRSScintillatorBarIndex.hh"
 
+#include "cetlib_except/exception.h"
 #include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh"
 #include "Offline/GlobalConstantsService/inc/PhysicsParams.hh"
 #include "Offline/CRVConditions/inc/CRVDigitizationPeriod.hh"
@@ -122,6 +123,18 @@ namespace mu2e
     art::Handle<CrvRecoPulseCollection> crvRecoPulseCollection;
     event.getByLabel(_crvRecoPulsesModuleLabel,"",crvRecoPulseCollection);
 
+    // any of these products can legitimately be absent from a file whose intermediate CRV
+    // collections were dropped, so say which one is missing instead of dereferencing an
+    // invalid handle (CrvMCHelper.cc and CrvCoincidenceClusterMatchMC_module.cc guard the same way)
+    if(!crvPhotonsCollection.isValid())
+      throw cet::exception("SIM")<<"mu2e::CrvPlot: no CrvPhotonsCollection with label "<<_crvPhotonsModuleLabel<<std::endl;
+    if(!crvSiPMChargesCollection.isValid())
+      throw cet::exception("SIM")<<"mu2e::CrvPlot: no CrvSiPMChargesCollection with label "<<_crvSiPMChargesModuleLabel<<std::endl;
+    if(!crvDigiCollection.isValid())
+      throw cet::exception("SIM")<<"mu2e::CrvPlot: no CrvDigiCollection with label "<<_crvDigiModuleLabel<<std::endl;
+    if(!crvRecoPulseCollection.isValid())
+      throw cet::exception("SIM")<<"mu2e::CrvPlot: no CrvRecoPulseCollection with label "<<_crvRecoPulsesModuleLabel<<std::endl;
+
     auto const& calib = _calib_h.get(event.id());
 
     GeomHandle<CosmicRayShield> CRS;
@@ -233,8 +246,10 @@ namespace mu2e
 //waveforms
       std::vector<std::vector<double> > ADCs[4];  //there can be multiple disconnected ADC waveforms per SiPM
       std::vector<std::vector<double> > times[4];
-      double scale[4]={NAN};
-      double maxADC[4]={NAN};
+      // NB: `double x[4]={NAN};` would set only element 0 to NAN and value-initialize 1-3 to 0.0,
+      // which silently defeats the isnan() "no digis" guard below for SiPMs 1-3
+      double scale[4]={NAN,NAN,NAN,NAN};
+      double maxADC[4]={NAN,NAN,NAN,NAN};
       CrvDigiCollection::const_iterator digis;
       for(digis=crvDigiCollection->begin(); digis!=crvDigiCollection->end(); ++digis)
       {

CRVResponse/src/CrvWidebandTest_module.cc

@@ -241,8 +241,23 @@ namespace mu2e
   CrvWidebandTest::~CrvWidebandTest()
   {
     delete[] _recoPEs;
+    delete[] _recoTime;
+    delete[] _fitStatus;
     delete[] _depositedEnergy;
     delete[] _coincidencePDGid;
+    delete[] _coincidenceTime;
+    delete[] _coincidencePosX;
+    delete[] _coincidencePosY;
+    delete[] _coincidencePosZ;
+    delete[] _trackSlope;
+    delete[] _trackIntercept;
+    delete[] _trackPoints;
+    delete[] _trackPEs;
+    delete[] _trackChi2;
+    delete[] _summaryPEs;
+    delete[] _summaryFWHMs;
+    delete[] _summarySignals;
+    delete[] _summaryChi2s;
   }
 
   void CrvWidebandTest::beginJob()
@@ -396,18 +411,22 @@ namespace mu2e
     //fits for the entire stack of modules/sectors (sectorType=0) and for individual modules/sectors (sectorType=1,...)
     for(int iSectorType=0; iSectorType<_nSectorTypes; ++iSectorType)
     {
+     //reset the per-event track variables for EVERY sector type, including the trigger sectors
+     //skipped below: these arrays are branched into the tree at dimension _nSectorTypes and are
+     //filled every event, so an un-reset trigger-sector entry would persist whatever new[] left there
+     _trackSlope[iSectorType]    =0;
+     _trackIntercept[iSectorType]=0;
+     _trackPEs[iSectorType]      =0;
+     _trackPoints[iSectorType]   =0;
+     _trackChi2[iSectorType]     =-1;
+
      if(std::find(_triggerSectorTypes.begin(),_triggerSectorTypes.end(), iSectorType) != _triggerSectorTypes.end()) continue;
 
      //initialize track variables
      float sumX     =0;
      float sumY     =0;
      float sumXY    =0;
      float sumYY    =0;
-     _trackSlope[iSectorType]    =0;
-     _trackIntercept[iSectorType]=0;
-     _trackPEs[iSectorType]      =0;
-     _trackPoints[iSectorType]   =0;
-     _trackChi2[iSectorType]     =-1;
 
       double counterWidth = counters.at(0)->getBarDetail().getHalfWidth()*2.0;
       int widthDirection = counters.at(0)->getBarDetail().getWidthDirection();  //assumes that all counters are oriented in the same way

CRVResponse/src/MakeCrvDigis.cc

@@ -9,10 +9,12 @@ void MakeCrvDigis::SetWaveform(const std::vector<double> &waveform, double ADCco
   _ADCs.resize(waveform.size());
   for(size_t i=0; i<waveform.size(); i++)
   {
-    int16_t ADC = static_cast<int16_t>(waveform[i]*ADCconversionFactor+pedestal+0.5);
-    if(ADC<minADC) ADC=minADC;
-    if(ADC>maxADC) ADC=maxADC;
-    _ADCs.at(i)=ADC;
+    // clamp before narrowing: an extreme sample cast to int16_t first would wrap around and
+    // land back inside [minADC,maxADC], so the saturation clamp would never see it
+    double ADCd = waveform[i]*ADCconversionFactor+pedestal+0.5;
+    if(ADCd<minADC) ADCd=minADC;
+    if(ADCd>maxADC) ADCd=maxADC;
+    _ADCs.at(i)=static_cast<int16_t>(ADCd);
   }
 
   int TDCtmp=lrint(startTime/digitizationPrecision);

CRVResponse/src/MakeCrvPhotons.cc

@@ -275,6 +275,10 @@ void LookupBin::Read(std::ifstream &lookupfile, const unsigned int &i)
   for(unsigned int i=0; i<nScintillatorScintillationBins; i++) _bins[0][i].Read(lookupfile,i);
   for(unsigned int i=0; i<nScintillatorCerenkovBins; i++)      _bins[1][i].Read(lookupfile,i);
   for(unsigned int i=0; i<nFiberCerenkovBins; i++)             _bins[2][i].Read(lookupfile,i);
+  // a truncated table would otherwise leave the bin boundaries and yield curves holding
+  // whatever was in the buffers, and the job would continue with silently wrong photon yields
+  if(!lookupfile) throw std::logic_error("Error while reading CRV lookup table file "+filename+
+                                         " (truncated or corrupt).");
   if(debug>0) std::cout<<"Done."<<std::endl;
 
   lookupfile.close();

CRVResponse/src/MakeCrvSiPMCharges.cc

@@ -45,7 +45,11 @@ std::pair<int,int> MakeCrvSiPMCharges::FindFiberPhotonsPixelId()
 {
   double x,y;
   _photonMap->GetRandom2(x,y);
-  return std::pair<int,int>(lrint(x),lrint(y));
+  // GetRandom2 samples the full width of the outermost bin, so lrint can return _nPixelsX/Y,
+  // one past the last pixel; that phantom pixel would carry its own discharge state
+  int ix=std::min(std::max(static_cast<int>(lrint(x)),0),_nPixelsX-1);
+  int iy=std::min(std::max(static_cast<int>(lrint(y)),0),_nPixelsY-1);
+  return std::pair<int,int>(ix,iy);
 }
 
 bool MakeCrvSiPMCharges::IsInactivePixelId(const std::pair<int,int> &pixelId)

CRVResponse/src/MakeCrvWaveforms.cc

@@ -48,6 +48,8 @@ void MakeCrvWaveforms::LoadSinglePEWaveform(const std::string &filename, double
   }
   f.close();
 
+  if(_singlePEWaveform.empty())
+    throw std::logic_error("Could not read any single PE waveform points from "+filename);
   _singlePEMaxVoltage = *std::max_element(_singlePEWaveform.begin(), _singlePEWaveform.end());
 }