From 53a09c47d9a218b4217e18f11178c80f9237a9e7 Mon Sep 17 00:00:00 2001 From: ehrlich-uva Date: Fri, 24 Apr 2026 21:23:28 -0500 Subject: [PATCH] small bug fixes in CRV simulation --- .../src/CrvSiPMChargeGenerator_module.cc | 6 ++- .../src/CrvStepsFromStepPointMCs_module.cc | 4 +- .../src/CrvWaveformsGenerator_module.cc | 52 +++++++++---------- CRVResponse/src/MakeCrvPhotons.cc | 1 + CRVResponse/src/MakeCrvSiPMCharges.cc | 2 +- 5 files changed, 35 insertions(+), 30 deletions(-) diff --git a/CRVResponse/src/CrvSiPMChargeGenerator_module.cc b/CRVResponse/src/CrvSiPMChargeGenerator_module.cc index 1aeaf21be2..a94c39e10c 100644 --- a/CRVResponse/src/CrvSiPMChargeGenerator_module.cc +++ b/CRVResponse/src/CrvSiPMChargeGenerator_module.cc @@ -215,7 +215,8 @@ namespace mu2e //time wrapping happened in the photon generator std::vector > 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 &photonTimes = crvPhotons->GetPhotons(); charges.emplace_back(time, charge, chargeInPEs, photonTimes[photonIndex]._step); } diff --git a/CRVResponse/src/CrvStepsFromStepPointMCs_module.cc b/CRVResponse/src/CrvStepsFromStepPointMCs_module.cc index 2981c53d9f..bfa66cdf9d 100644 --- a/CRVResponse/src/CrvStepsFromStepPointMCs_module.cc +++ b/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; diff --git a/CRVResponse/src/CrvWaveformsGenerator_module.cc b/CRVResponse/src/CrvWaveformsGenerator_module.cc index a6a484f4b4..3ecdb65dd8 100644 --- a/CRVResponse/src/CrvWaveformsGenerator_module.cc +++ b/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 > timesAndCharges; - }; + typedef std::vector > 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 ×AndCharges, + void FindChargeClusters(const std::vector &singleCharges, std::vector &chargeClusters, double timeOffset); bool SingleWaveformStart(std::vector &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 ×AndCharges = iter->GetCharges(); + const std::vector &singleCharges = iter->GetCharges(); //zero suppressed data std::vector chargeClusters; - FindChargeClusters(timesAndCharges, chargeClusters, timeOffset); + FindChargeClusters(singleCharges, chargeClusters, timeOffset); - for(size_t iCluster=0; iClusterGetSinglePEMaxVoltage()<_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 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 > steps; //use a set to remove duplicate steppoints std::map, int> simparticles; - for(size_t j=0; j=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 > timesAndChargesNZS; - for(size_t iCharge=0; iCharge=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 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); } } diff --git a/CRVResponse/src/MakeCrvPhotons.cc b/CRVResponse/src/MakeCrvPhotons.cc index 790d0fb239..a28844e61a 100644 --- a/CRVResponse/src/MakeCrvPhotons.cc +++ b/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 { diff --git a/CRVResponse/src/MakeCrvSiPMCharges.cc b/CRVResponse/src/MakeCrvSiPMCharges.cc index 92c8a1840c..921a539265 100644 --- a/CRVResponse/src/MakeCrvSiPMCharges.cc +++ b/CRVResponse/src/MakeCrvSiPMCharges.cc @@ -45,7 +45,7 @@ std::pair MakeCrvSiPMCharges::FindFiberPhotonsPixelId() { double x,y; _photonMap->GetRandom2(x,y); - return std::pair(x,y); + return std::pair(lrint(x),lrint(y)); } bool MakeCrvSiPMCharges::IsInactivePixelId(const std::pair &pixelId)