small bug fixes in CRV simulation

CRVResponse/src/CrvSiPMChargeGenerator_module.cc

@@ -215,7 +215,8 @@ namespace mu2e
 
         //time wrapping happened in the photon generator
         std::vector<std::pair<double,size_t> > photonTimesNew;   //pair of photon time and index in the original photon vector
-        CrvPhotonsCollection::const_iterator crvPhotons;
+        CrvPhotonsCollection::const_iterator crvPhotons;   //after this loop, it will point to the correct index in the photon collection, if barIndex and SiPM is found
+                                                           //otherwise it will be crvPhotonsCollection->end()
         for(crvPhotons=crvPhotonsCollection->begin(); crvPhotons!=crvPhotonsCollection->end(); crvPhotons++)
         {
           if(crvPhotons->GetScintillatorBarIndex()==barIndex && crvPhotons->GetSiPMNumber()==(int)SiPM)
@@ -249,6 +250,9 @@ namespace mu2e
             bool darkNoise=responseIter->_darkNoise;
             if(!darkNoise)
             {
+              //crvPhotons should always point to the correct index in the photon collection, if not dark noise
+              //this exception line protects against future unintended code changes
+              if(crvPhotons==crvPhotonsCollection->end()) throw cet::exception("CrvSiPMChargeGenerator") << "invalid iterator of the photon collection";
               const std::vector<CrvPhotons::SinglePhoton> &photonTimes = crvPhotons->GetPhotons();
               charges.emplace_back(time, charge, chargeInPEs, photonTimes[photonIndex]._step);
             }

CRVResponse/src/CrvStepsFromStepPointMCs_module.cc

@@ -306,7 +306,9 @@ namespace mu2e
         //calculating the end time of the last step
         double avgEnergy = 0.5*(endEnergyBefore+endEnergyAfter);
         double avgGamma  = avgEnergy/mass;
-        double avgBeta   = sqrt(1.0-1.0/(avgGamma*avgGamma));
+        double radicand  = 1.0-1.0/(avgGamma*avgGamma);
+        if(radicand<0) radicand=0;  //shouldn't happen except due to numerical precision
+        double avgBeta   = sqrt(radicand);
         double velocity  = avgBeta*CLHEP::c_light;
         endTime         += last->stepLength()/velocity;
 

CRVResponse/src/CrvWaveformsGenerator_module.cc

@@ -103,12 +103,11 @@ namespace mu2e
 
     static constexpr int nDigiPeriods = 4; //number of digitization periods to check for charges
 
-    struct ChargeCluster
-    {
-      std::vector<std::pair<double,double> > timesAndCharges;
-    };
+    typedef std::vector<std::pair<double,double> > ChargeCluster;  //first: time, second: charge
+                                                                   //this structure is required by the MakeCrvWaveforms class,
+                                                                   //which cannot be modified, since it is used in other projects.
 
-    void FindChargeClusters(const std::vector<CrvSiPMCharges::SingleCharge> &timesAndCharges,
+    void FindChargeClusters(const std::vector<CrvSiPMCharges::SingleCharge> &singleCharges,
                             std::vector<ChargeCluster> &chargeClusters, double timeOffset);
     bool SingleWaveformStart(std::vector<double> &fullWaveform, size_t i);
   };
@@ -221,26 +220,25 @@ namespace mu2e
       double TDC0time=eventWindowStart+digitizationPointShiftFEB;  //that's the time when TDC=0 for this FEB
 
       //charges and times
-      const std::vector<CrvSiPMCharges::SingleCharge> &timesAndCharges = iter->GetCharges();
+      const std::vector<CrvSiPMCharges::SingleCharge> &singleCharges = iter->GetCharges();
 
       //zero suppressed data
       std::vector<ChargeCluster> chargeClusters;
-      FindChargeClusters(timesAndCharges, chargeClusters, timeOffset);
+      FindChargeClusters(singleCharges, chargeClusters, timeOffset);
 
-      for(size_t iCluster=0; iCluster<chargeClusters.size(); ++iCluster)
+      for(const ChargeCluster &thisChargeCluster : chargeClusters)
       {
         //if the number of charges in this cluster cannot achieve the minimum voltage, skip this cluster
-        if(chargeClusters[iCluster].timesAndCharges.size()*_makeCrvWaveforms->GetSinglePEMaxVoltage()<_minVoltage) continue;
+        if(thisChargeCluster.size()*_makeCrvWaveforms->GetSinglePEMaxVoltage()<_minVoltage) continue;
 
         //find the TDC time when the first charge occurs (adjusted for this FEB)
-        double firstChargeTime=chargeClusters[iCluster].timesAndCharges.front().first;
+        double firstChargeTime=thisChargeCluster.front().first;
         firstChargeTime-=1.0*CRVDigitizationPeriod;  //start somewhere before the first charge
         double TDCstartTime=ceil((firstChargeTime-TDC0time)/CRVDigitizationPeriod) * CRVDigitizationPeriod + TDC0time;
 
         //first create the full waveform
         std::vector<double> fullWaveform;
-        _makeCrvWaveforms->MakeWaveform(chargeClusters[iCluster].timesAndCharges,
-                                        fullWaveform, TDCstartTime, CRVDigitizationPeriod);
+        _makeCrvWaveforms->MakeWaveform(thisChargeCluster, fullWaveform, TDCstartTime, CRVDigitizationPeriod);
         _makeCrvWaveforms->AddElectronicNoise(fullWaveform, _noise, _randGaussQ);
 
         //break the waveform apart into short pieces (_numberSampleZS) and apply the zero suppression
@@ -268,13 +266,13 @@ namespace mu2e
             //collect CrvSteps and SimParticles responsible for this single waveform
             std::set<art::Ptr<CrvStep> > steps;  //use a set to remove duplicate steppoints
             std::map<art::Ptr<SimParticle>, int> simparticles;
-            for(size_t j=0; j<timesAndCharges.size(); ++j)
+            for(size_t j=0; j<singleCharges.size(); ++j)
             {
-              if(timesAndCharges[j]._time>=digiStartTime-_singlePEWaveformMaxTime &&
-                 timesAndCharges[j]._time<=digiStartTime+_numberSamplesZS*CRVDigitizationPeriod)
+              if(singleCharges[j]._time+timeOffset>=digiStartTime-_singlePEWaveformMaxTime &&
+                 singleCharges[j]._time+timeOffset<=digiStartTime+_numberSamplesZS*CRVDigitizationPeriod)
               {
-                steps.insert(timesAndCharges[j]._step);
-                if(timesAndCharges[j]._step.isNonnull()) simparticles[timesAndCharges[j]._step->simParticle()]++;
+                steps.insert(singleCharges[j]._step);
+                if(singleCharges[j]._step.isNonnull()) simparticles[singleCharges[j]._step->simParticle()]++;
               }
             }
 
@@ -310,19 +308,19 @@ namespace mu2e
       {
         //non-zero suppressed data
         //use only charges for the first 134 samples (_numberSamplesNZS)
-        std::vector<std::pair<double,double> > timesAndChargesNZS;
-        for(size_t iCharge=0; iCharge<timesAndCharges.size(); ++iCharge)
+        ChargeCluster chargeClusterNZS;
+        for(size_t iCharge=0; iCharge<singleCharges.size(); ++iCharge)
         {
-          if(timesAndCharges[iCharge]._time>=TDC0time-_singlePEWaveformMaxTime &&
-             timesAndCharges[iCharge]._time<=TDC0time+CRVDigitizationPeriod*_numberSamplesNZS)
+          if(singleCharges[iCharge]._time+timeOffset>=TDC0time-_singlePEWaveformMaxTime &&
+             singleCharges[iCharge]._time+timeOffset<=TDC0time+CRVDigitizationPeriod*_numberSamplesNZS)
           {
-            timesAndChargesNZS.emplace_back(timesAndCharges[iCharge]._time,timesAndCharges[iCharge]._charge);
+            chargeClusterNZS.emplace_back(singleCharges[iCharge]._time+timeOffset,singleCharges[iCharge]._charge);
           }
         }
 
         //create NZS waveform
         std::vector<double> fullWaveformNZS;
-        _makeCrvWaveforms->MakeWaveform(timesAndChargesNZS, fullWaveformNZS, TDC0time, CRVDigitizationPeriod);
+        _makeCrvWaveforms->MakeWaveform(chargeClusterNZS, fullWaveformNZS, TDC0time, CRVDigitizationPeriod);
         //record first 134 samples (_numberSamplesNZS)
         //assumed to be always within eventWindowStart and eventWindowEnd
         fullWaveformNZS.resize(_numberSamplesNZS);
@@ -352,20 +350,20 @@ namespace mu2e
       if(chargeClusters.empty())
       {
         chargeClusters.resize(1);
-        chargeClusters.back().timesAndCharges.reserve(timesAndCharges.size());
+        chargeClusters.back().reserve(timesAndCharges.size());
       }
       else
       {
-        if(timeTmp-chargeClusters.back().timesAndCharges.back().first>_singlePEWaveformMaxTime+nDigiPeriods*CRVDigitizationPeriod)
+        if(timeTmp-chargeClusters.back().back().first>_singlePEWaveformMaxTime+nDigiPeriods*CRVDigitizationPeriod)
         //if the difference b/w the time of the next charge and the time of the last charge
         //is greater than a full single PE waveform plus four additional digitization periods
         //-->start a new charge cluster
         {
           chargeClusters.resize(chargeClusters.size()+1);
-          chargeClusters.back().timesAndCharges.reserve(timesAndCharges.size());
+          chargeClusters.back().reserve(timesAndCharges.size());
         }
       }
-      chargeClusters.back().timesAndCharges.emplace_back(timeTmp,timesAndCharges[i]._charge);
+      chargeClusters.back().emplace_back(timeTmp,timesAndCharges[i]._charge);
     }
   }
 

CRVResponse/src/MakeCrvPhotons.cc

@@ -515,6 +515,7 @@ int MakeCrvPhotons::GetNumberOfPhotonsFromAverage(double average, int nSteps)  /
   {
     double sigma = std::sqrt(average);
     nPhotons = lrint(_randGaussQ.fire(average,sigma)/nSteps);
+    if(nPhotons<0) nPhotons=0;
   }
   else
   {

CRVResponse/src/MakeCrvSiPMCharges.cc

@@ -45,7 +45,7 @@ std::pair<int,int> MakeCrvSiPMCharges::FindFiberPhotonsPixelId()
 {
   double x,y;
   _photonMap->GetRandom2(x,y);
-  return std::pair<int,int>(x,y);
+  return std::pair<int,int>(lrint(x),lrint(y));
 }
 
 bool MakeCrvSiPMCharges::IsInactivePixelId(const std::pair<int,int> &pixelId)