From df1a16e2dc2dfee7b7e3d17c5b547db302a44a92 Mon Sep 17 00:00:00 2001 From: echenard Date: Tue, 4 Aug 2026 22:13:37 -0500 Subject: [PATCH 1/7] Improved CaloMC modules --- CaloMC/inc/ShowerStepUtil.hh | 4 +- CaloMC/src/CaloHitTruthMatch_module.cc | 247 ++++++-------- CaloMC/src/CaloShowerROMaker_module.cc | 204 +++++------- CaloMC/src/CaloShowerStepMaker_module.cc | 408 ++++++++--------------- 4 files changed, 320 insertions(+), 543 deletions(-) diff --git a/CaloMC/inc/ShowerStepUtil.hh b/CaloMC/inc/ShowerStepUtil.hh index 23a23f0555..6641a6c859 100644 --- a/CaloMC/inc/ShowerStepUtil.hh +++ b/CaloMC/inc/ShowerStepUtil.hh @@ -17,7 +17,7 @@ namespace mu2e { imax_(imax),type_(type),n_(imax,0),eDepG4_(imax,0),eDepVis_(imax,0), pIn_(imax,0),time_(imax,0),t0_(imax,0),x_(imax,0),y_(imax,0),z_(imax,0), w_(imax,0),pos_(0,0,0) - {}; + {} void add(unsigned i, double eDepG4, double eDepVis, double time, double momentum, const CLHEP::Hep3Vector& pos); void reset(unsigned i); @@ -34,6 +34,8 @@ namespace mu2e { private: + void rangeCheck(unsigned i) const; + unsigned imax_; weight_type type_; std::vector n_; diff --git a/CaloMC/src/CaloHitTruthMatch_module.cc b/CaloMC/src/CaloHitTruthMatch_module.cc index 850d2ae178..79867c1004 100644 --- a/CaloMC/src/CaloHitTruthMatch_module.cc +++ b/CaloMC/src/CaloHitTruthMatch_module.cc @@ -3,8 +3,6 @@ // #include "art/Framework/Core/EDProducer.h" #include "art/Framework/Principal/Event.h" -#include "art_root_io/TFileService.h" -#include "art_root_io/TFileDirectory.h" #include "fhiclcpp/types/Atom.h" #include "Offline/MCDataProducts/inc/CaloEDepMC.hh" @@ -17,13 +15,13 @@ #include "Offline/RecoDataProducts/inc/CaloHit.hh" #include "Offline/Mu2eUtilities/inc/CaloPulseShape.hh" -#include "TH2F.h" -#include "TFile.h" - -#include -#include +#include #include +#include #include +#include +#include +#include #include @@ -36,26 +34,28 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - fhicl::Atom caloShowerSimCollection { Name("caloShowerSimCollection"), Comment("Name of caloShowerSim Collection") }; - fhicl::Atom caloHitCollection { Name("caloHitCollection"), Comment("Name of CaloHit collection") }; - fhicl::Atom primaryParticle { Name("primaryParticle"), Comment("PrimaryParticle producer")}; - fhicl::Atom pulseFileName { Name("pulseFileName"), Comment("Calo pulse file name") }; - fhicl::Atom pulseHistName { Name("pulseHistName"), Comment("Calo pulse hist name") }; - fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; - fhicl::Atom minAmplitude { Name("minAmplitude"), Comment("Minimum amplitude of waveform to define hit length") }; - fhicl::Atom fillDetailedMC { Name("fillDetailedMC"), Comment("Fill SimParticle - SimShower Assn map")}; - fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; + fhicl::Atom caloShowerSimCollection { Name("caloShowerSimCollection"), Comment("Name of caloShowerSim Collection") }; + fhicl::Atom caloHitCollection { Name("caloHitCollection"), Comment("Name of CaloHit collection") }; + fhicl::Atom primaryParticle { Name("primaryParticle"), Comment("PrimaryParticle producer")}; + fhicl::Atom pulseFileName { Name("pulseFileName"), Comment("Calo pulse file name") }; + fhicl::Atom pulseHistName { Name("pulseHistName"), Comment("Calo pulse hist name") }; + fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; + fhicl::Atom deltaTimeMinus { Name("deltaTimeMinus"), Comment("Max time (ns) a MC hit may precede the reco hit to be matched"), 100.0 }; + fhicl::Atom minAmplitude { Name("minAmplitude"), Comment("Minimum amplitude of waveform to define hit length") }; + fhicl::Atom fillDetailedMC { Name("fillDetailedMC"), Comment("Fill SimParticle - SimShower Assn map")}; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; }; explicit CaloHitTruthMatch(const art::EDProducer::Table& config) : EDProducer{config}, - caloShowerSimToken_ {consumes (config().caloShowerSimCollection())}, + caloShowerSimToken_ {consumes(config().caloShowerSimCollection())}, caloHitToken_ {consumes(config().caloHitCollection())}, ppToken_ {consumes(config().primaryParticle())}, pulseFileName_ (config().pulseFileName()), pulseHistName_ (config().pulseHistName()), digiSampling_ (config().digiSampling()), + deltaTimeMinus_ (config().deltaTimeMinus()), minAmplitude_ (config().minAmplitude()), fillDetailedMC_ (config().fillDetailedMC()), diagLevel_ (config().diagLevel()) @@ -65,71 +65,40 @@ namespace mu2e { if (fillDetailedMC_) produces(); } - void beginJob() override; - void produce(art::Event& e) override; - void beginRun(art::Run& aRun) override; + void beginRun(art::Run&) override; + void produce (art::Event& e) override; private: using SimParticlePtr = art::Ptr; - void makeTruthMatch (art::Event&, CaloHitMCCollection&, CaloHitMCTruthAssn&, CaloShowerMCTruthAssn&, const PrimaryParticle&); - void fillEdeps (const PrimaryParticle& primaryParticle, std::vector& edeps, const CaloShowerSim* showerSim); - void diag (const CaloShowerSim*, const CaloHit& ); - + void makeTruthMatch(art::Event&, CaloHitMCCollection&, CaloHitMCTruthAssn&, CaloShowerMCTruthAssn&, const PrimaryParticle&); + void fillEdeps (const PrimaryParticle& primaryParticle, std::vector& edeps, const CaloShowerSim* showerSim); const art::ProductToken caloShowerSimToken_; const art::ProductToken caloHitToken_; const art::ProductToken ppToken_; std::string pulseFileName_; std::string pulseHistName_; - double deltaTimeMinus_; double digiSampling_; + double deltaTimeMinus_; double minAmplitude_; bool fillDetailedMC_; std::vector wf_; - size_t wfBinMax_; + std::size_t wfBinMax_{0}; int diagLevel_; - - //some diagnostic histograms - TH1F* hTime_; - TH1F* hTime2_; - TH2F* hTime2d_; - TH2F* hEner2d_; - TH2F* hEnerTime_; - TH2F* hdEdT_; - TH1F* hChi2_; }; - - - - //-------------------------------------------------------------------- - void CaloHitTruthMatch::beginJob() - { - if ( diagLevel_ > 2) - { - art::ServiceHandle tfs; - hTime_ = tfs->make("hTime", "delta Time", 2000, -20., 180); - hTime2_ = tfs->make("hTime2", "delta Time", 2000, -20., 180); - hTime2d_ = tfs->make("hTime2d", "Reco vs Gen time", 200,500,1700, 200,500,1700); - hEner2d_ = tfs->make("hEner2d", "Reco vs gen Ener", 200,0,40, 200,0.,40); - hEnerTime_ = tfs->make("hTimeEner","delta Time vs Ener",500,0,100, 300,-50.,250); - hdEdT_ = tfs->make("hdEdt", "delta Time vs delta Ener",100,-10,70, 170,-10.,160); - hChi2_ = tfs->make("hChi2", "chi2 large dE", 50, 0., 10); - } - } - //----------------------------------------------------------------------------- - void CaloHitTruthMatch::beginRun(art::Run& aRun) + void CaloHitTruthMatch::beginRun(art::Run&) { CaloPulseShape cps(pulseFileName_,pulseHistName_,digiSampling_); cps.buildShapes(); - wf_ = cps.digitizedPulse(0); + wf_ = cps.digitizedPulse(0); wfBinMax_ = std::distance(wf_.begin(),std::max_element(wf_.begin(),wf_.end())); } @@ -137,12 +106,11 @@ namespace mu2e { //-------------------------------------------------------------------- void CaloHitTruthMatch::produce(art::Event& event) { - auto pph = event.getValidHandle(ppToken_); - auto const& primaryParticles = *pph; + const auto& primaryParticles = *event.getValidHandle(ppToken_); - std::unique_ptr caloHitMCs(new CaloHitMCCollection); - std::unique_ptr caloHitMCTruth(new CaloHitMCTruthAssn); - std::unique_ptr caloShowerMCTruth(new CaloShowerMCTruthAssn); + auto caloHitMCs = std::make_unique(); + auto caloHitMCTruth = std::make_unique(); + auto caloShowerMCTruth = std::make_unique(); makeTruthMatch(event, *caloHitMCs, *caloHitMCTruth, *caloShowerMCTruth, primaryParticles); @@ -153,74 +121,77 @@ namespace mu2e { - // perform the association with the following rules: - // MCtime must be inside the window [recoTime-deltaTimeMinus, recoTime+deltaTimePlus] to be associated to RecoHit _unless_ - // MCtime is already in the window of the next hit, in which case it is associate to this one. + // Association rules: + // MCtime must be inside [recoTime-deltaTimeMinus, recoTime+deltaTimePlus] to be associated to a RecoHit, + // _unless_ MCtime is already inside the window of the next hit, in which case it goes to that one. //-------------------------------------------------------------------- void CaloHitTruthMatch::makeTruthMatch(art::Event& event, CaloHitMCCollection& caloHitMCs, - CaloHitMCTruthAssn& CaloHitTruthMatch, CaloShowerMCTruthAssn& caloShowerTruthMatch, + CaloHitMCTruthAssn& caloHitTruthMatch, CaloShowerMCTruthAssn& caloShowerTruthMatch, const PrimaryParticle& primaryParticle) { - - int nMatched(0); + int nMatched(0); double totalEnergyMatched(0); - // access collections of CaloHits and caloShowerHits, need art::ProductID for creating art::Ptr - art::ProductID hitMCProductID(event.getProductID()); + const art::ProductID hitMCProductID(event.getProductID()); const art::EDProductGetter* hitMCProductGetter = event.productGetter(hitMCProductID); + const auto caloHitHandle = event.getValidHandle(caloHitToken_); const auto caloShowerSimHandle = event.getValidHandle(caloShowerSimToken_); const auto& caloHits(*caloHitHandle); const auto& caloShowerSims(*caloShowerSimHandle); - // sort the caloHits and caloShowerSim per crystal and then per time for each crystal to help with the matching algorithm. - std::map> caloHitMap; - std::map> caloShowerSimsMap; - for (const auto& caloHit: caloHits) caloHitMap[caloHit.crystalID()].push_back(&caloHit); - for (const auto& caloShowerSim: caloShowerSims) caloShowerSimsMap[caloShowerSim.crystalID()].push_back(&caloShowerSim); + // sort caloHits and caloShowerSims per crystal, then per time, to help the matching + std::unordered_map> caloHitMap; + std::unordered_map> caloShowerSimsMap; + for (const auto& caloHit : caloHits) caloHitMap[caloHit.crystalID()].push_back(&caloHit); + for (const auto& caloShowerSim : caloShowerSims) caloShowerSimsMap[caloShowerSim.crystalID()].push_back(&caloShowerSim); + + const auto byTime = [](const auto* a, const auto* b){return a->time() < b->time();}; + for (auto& [id, v] : caloHitMap) std::sort(v.begin(),v.end(),byTime); + for (auto& [id, v] : caloShowerSimsMap) std::sort(v.begin(),v.end(),byTime); - for (auto &kv : caloHitMap) std::sort(kv.second.begin(),kv.second.end(), [](auto const a, auto const b){return a->time() < b->time();}); - for (auto &kv : caloShowerSimsMap) std::sort(kv.second.begin(),kv.second.end(), [](auto const a, auto const b){return a->time() < b->time();}); + // O(1) shower-pointer -> index lookup for the detailed-MC association (was an O(n) scan per match) + std::unordered_map showerIndex; + if (fillDetailedMC_) + { + showerIndex.reserve(caloShowerSims.size()); + for (std::size_t i=0; i(caloHitHandle,ihit); - auto showerIt = sortedSims.begin(); - auto showerItEnd = sortedSims.end(); + const CaloHit& hit = caloHits[ihit]; + const auto& sortedHits = caloHitMap.at(hit.crystalID()); + const auto& sortedSims = caloShowerSimsMap.at(hit.crystalID()); + const auto hitIt = std::find(sortedHits.begin(),sortedHits.end(),&hit); + const auto hitNextIt = std::next(hitIt); + const auto hitPtr = art::Ptr(caloHitHandle,ihit); if (diagLevel_ > 2) - for (const auto& shower : sortedSims) std::cout<<"[CaloHitTruthMatch] Sim shower id/time/energy="<crystalID() - <<" / "<time()<<" / "<energyDep()<time()<<" / "<energyDep()<time()- (*hitIt)->time() - 2*digiSampling_ < deltaTimePlus ) - { - deltaTimePlus = (*hitNextIt)->time()- (*hitIt)->time() - 2*digiSampling_; - } - if (diagLevel_ > 2) std::cout<<"[CaloHitTruthMatch] inspect hit id/time/energy/length "<time() - hit.time() - 2*digiSampling_ < deltaTimePlus) + deltaTimePlus = (*hitNextIt)->time() - hit.time() - 2*digiSampling_; - // FIXME - // deltaTimePlus = 200.; - double deltaTimeMin = 100.; // 2*digiSampling_ + if (diagLevel_ > 2) std::cout<<"[CaloHitTruthMatch] inspect hit id/time/energy/length "< edeps; - while (showerIt != showerItEnd && ( (*showerIt)->time() < (*hitIt)->time() - deltaTimeMin ) ) ++showerIt; - while (showerIt != showerItEnd && ( (*showerIt)->time() < (*hitIt)->time() + deltaTimePlus) ) + auto showerIt = sortedSims.begin(); + while (showerIt != sortedSims.end() && (*showerIt)->time() < hit.time() - deltaTimeMinus_) ++showerIt; + while (showerIt != sortedSims.end() && (*showerIt)->time() < hit.time() + deltaTimePlus) { hitIsMatched = true; const CaloShowerSim* showerSim = *showerIt; @@ -228,50 +199,43 @@ namespace mu2e { if (fillDetailedMC_) { - size_t idxShower(0); - while (idxShower < caloShowerSims.size()) {if (&caloShowerSims[idxShower]==*showerIt) break; ++idxShower;} - auto ShowerSimPtr = art::Ptr(caloShowerSimHandle,idxShower); - caloShowerTruthMatch.addSingle(hitPtr, showerSim->sim(), ShowerSimPtr); + const auto showerSimPtr = art::Ptr(caloShowerSimHandle, showerIndex.at(showerSim)); + caloShowerTruthMatch.addSingle(hitPtr, showerSim->sim(), showerSimPtr); } - if (diagLevel_ > 1) diag(showerSim,hit); - if (diagLevel_ > 2) std::cout<<"[CaloHitTruthMatch] matched shower id/time/energyDep()= "<crystalID() + if (diagLevel_ > 2) std::cout<<"[CaloHitTruthMatch] matched shower id/time/energyDep= "<crystalID() <<" / "<time()<<" / "<energyDep()< MChit association std::sort(edeps.begin(),edeps.end(),[](const auto& a, const auto& b){return a.energyDep() > b.energyDep();}); - caloHitMCs.emplace_back(CaloHitMC(std::move(edeps),hit.crystalID())); - - art::Ptr hitMCPtr = art::Ptr(hitMCProductID, caloHitMCs.size()-1, hitMCProductGetter); - CaloHitTruthMatch.addSingle(hitPtr,hitMCPtr); - - if (hitIsMatched) {totalEnergyMatched += (*hitIt)->energyDep();++nMatched;} - ++hitIt; + caloHitMCs.emplace_back(std::move(edeps),hit.crystalID()); - if (diagLevel_ > 2 && !hitIsMatched) std::cout<<"[CaloHitTruthMatch] hit not matched"<(hitMCProductID, caloHitMCs.size()-1, hitMCProductGetter); + caloHitTruthMatch.addSingle(hitPtr,hitMCPtr); + if (hitIsMatched) {totalEnergyMatched += hit.energyDep(); ++nMatched;} + else if (diagLevel_ > 2) std::cout<<"[CaloHitTruthMatch] hit not matched"< 0) std::cout<<"[CaloHitTruthMatch] total particles / energy matched = "< 0) std::cout<<"[CaloHitTruthMatch] total particles / energy matched = "<& edeps, const CaloShowerSim* showerSim) { - // check if there is already a caloEdep object with same caloShowerSim's SimParticle - auto it = edeps.begin(); - while (it != edeps.end()) {if (it->sim() == showerSim->sim()) break; ++it;} + // if there is already a CaloEDepMC for this SimParticle, merge; otherwise create a new one + auto it = std::find_if(edeps.begin(), edeps.end(), + [&](const CaloEDepMC& e){return e.sim() == showerSim->sim();}); - // if found the add caloShowerSim to matching caloEdep object, otherwise create new caloEdep - if (it!= edeps.end()) + if (it != edeps.end()) { - it->addEDep(showerSim->energyDep()); + it->addEDep (showerSim->energyDep()); it->addEDepG4(showerSim->energyDepG4()); - it->addTime(showerSim->time()); - it->addMom(showerSim->momentumIn()); + it->addTime (showerSim->time()); + it->addMom (showerSim->momentumIn()); } else { @@ -281,34 +245,11 @@ namespace mu2e { MCRelationship mcr(spp,showerSim->sim()); if (mcr > mcrel) mcrel = mcr; } - edeps.emplace_back(CaloEDepMC(showerSim->sim(),showerSim->energyDep(),showerSim->energyDepG4(), - showerSim->time(),showerSim->momentumIn(),mcrel)); + edeps.emplace_back(showerSim->sim(),showerSim->energyDep(),showerSim->energyDepG4(), + showerSim->time(),showerSim->momentumIn(),mcrel); } } - - - //-------------------------------------------------------------------- - void CaloHitTruthMatch::diag(const CaloShowerSim* shower, const CaloHit& hit) - { - hTime_->Fill(shower->time()-hit.time()); - hEnerTime_->Fill(shower->energyDep(),shower->time()-hit.time()); - hTime2d_->Fill(shower->time(),hit.time()); - hEner2d_->Fill(shower->energyDep(),hit.energyDep()); - hdEdT_->Fill(hit.energyDep()-shower->energyDep(),shower->time()-hit.time()); - - if (shower->energyDep() > 5) hTime2_->Fill(shower->time()-hit.time()); - - double deltaE = std::abs(shower->energyDep()-hit.energyDep()); - if (deltaE > 5 && shower->energyDep() > 5) - if (!hit.recoCaloDigis().empty()) hChi2_->Fill(hit.recoCaloDigis().at(0)->chi2()/hit.recoCaloDigis().at(0)->ndf()); - } - - } -using mu2e::CaloHitTruthMatch; -DEFINE_ART_MODULE(CaloHitTruthMatch) - - - +DEFINE_ART_MODULE(mu2e::CaloHitTruthMatch) diff --git a/CaloMC/src/CaloShowerROMaker_module.cc b/CaloMC/src/CaloShowerROMaker_module.cc index 29b48d4138..09d3a55bb5 100644 --- a/CaloMC/src/CaloShowerROMaker_module.cc +++ b/CaloMC/src/CaloShowerROMaker_module.cc @@ -1,11 +1,12 @@ // // Transform the energy deposited in the scintillator into photo-electrons (PE) seen by the photosensor. -// Includes corrections from Birks law, longitudinal response uniformity and photo-statistcs fluctuations. +// Includes corrections from Birks law, longitudinal response uniformity and photo-statistics fluctuations. // The PE are generated individually and corrected for transit time. // #include "art/Framework/Core/EDProducer.h" #include "art/Framework/Principal/Event.h" #include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" #include "fhiclcpp/types/Sequence.h" #include "Offline/CaloMC/inc/CaloPhotonPropagation.hh" @@ -25,14 +26,14 @@ #include "Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh" #include "Offline/SeedService/inc/SeedService.hh" #include "CLHEP/Random/RandPoissonQ.h" -#include "CLHEP/Random/RandFlat.h" -#include -#include +#include #include +#include #include -#include +#include #include +#include namespace { @@ -42,7 +43,6 @@ namespace { StepEntry(const art::Ptr& step, float edepCorr, float timeCorr) : step_(step),edepCorr_(edepCorr),timeCorr_(timeCorr) {} - art::Ptr step_; float edepCorr_,timeCorr_; }; @@ -52,29 +52,25 @@ namespace { SimParticleSummary(const art::Ptr& step, float edepCorr, float timeCorr) : steps_{step},edepCorr_(edepCorr),timeCorr_(timeCorr) {} - void add(const art::Ptr& step, float edepCorr, float timeCorr) { steps_.push_back(step); edepCorr_ += edepCorr; - timeCorr_ = std::min(timeCorr,timeCorr_); + timeCorr_ = std::min(timeCorr,timeCorr_); } - std::vector> steps_; float edepCorr_,timeCorr_; }; - struct diagSummary + struct DiagSummary { - diagSummary() : totSteps(0),totNPE(0),totEdep(0.),totEdepCorr(0.),totEdepNPE(0.) {}; - int totSteps,totNPE; - float totEdep,totEdepCorr,totEdepNPE; + int totSteps{0}, totNPE{0}; + float totEdep{0.f}, totEdepCorr{0.f}, totEdepNPE{0.f}; }; } - namespace mu2e { class CaloShowerROMaker : public art::EDProducer @@ -101,8 +97,8 @@ namespace mu2e { explicit CaloShowerROMaker(const art::EDProducer::Table& config) : EDProducer{config}, - ewMarkerTag_ (config().ewMarkerTag()), - pbtmcTag_ (config().pbtmcTag()), + ewMarkerToken_ {consumes(config().ewMarkerTag())}, + pbtmcToken_ {consumes(config().pbtmcTag())}, digitizationStart_ (config().digitizationStart()), digitizationEnd_ (config().digitizationEnd()), digitizationBuffer_ (config().digitizationBuffer()), @@ -115,31 +111,29 @@ namespace mu2e { randPoisson_ (engine_), photonProp_ (config().propagationFileName(),config().propagationHistName(),engine_) { - - for (auto const& tag : config().caloShowerStepCollection()) crystalShowerTokens_.push_back(consumes(tag)); - consumes(ewMarkerTag_); - consumes(pbtmcTag_); + for (const auto& tag : config().caloShowerStepCollection()) + crystalShowerTokens_.push_back(consumes(tag)); produces(); produces(); } void beginRun(art::Run& aRun) override; - void produce(art::Event& e) override; + void produce(art::Event& e) override; private: using StepHandles = std::vector>; - void makeReadoutHits (const StepHandles&, CaloShowerROCollection&, CaloShowerSimCollection&, const CalSimParams&, - const EventWindowMarker&, const ProtonBunchTimeMC&, float timeFromProtonsToDRMarker); - float LRUCorrection (float normalizedPosZ, float edepInit, float lru); - void dumpCaloShowerSim (const CaloShowerSimCollection& caloShowerSims); + void makeReadoutHits (const StepHandles&, CaloShowerROCollection&, CaloShowerSimCollection&, const CalSimParams&, + const EventWindowMarker&, const ProtonBunchTimeMC&, float timeFromProtonsToDRMarker); + float LRUCorrection (float normalizedPosZ, float edepInit, float lru) const; + void dumpCaloShowerSim(const CaloShowerSimCollection& caloShowerSims) const; - ProditionsHandle calCrystalConds_; + ProditionsHandle calCrystalConds_; std::vector> crystalShowerTokens_; - art::InputTag ewMarkerTag_; - art::InputTag pbtmcTag_; + art::ProductToken ewMarkerToken_; + art::ProductToken pbtmcToken_; float digitizationStart_; float digitizationEnd_; float digitizationBuffer_; @@ -156,7 +150,7 @@ namespace mu2e { //----------------------------------------------- - void CaloShowerROMaker::beginRun(art::Run& aRun) + void CaloShowerROMaker::beginRun(art::Run&) { photonProp_.buildTable(); } @@ -167,34 +161,26 @@ namespace mu2e { { if (diagLevel_ > 0) std::cout << "[CaloShowerROMaker::produce] begin" << std::endl; - //get Event window and bunch timing info - art::Handle ewMarkerHandle; - event.getByLabel(ewMarkerTag_, ewMarkerHandle); - const EventWindowMarker& ewMarker(*ewMarkerHandle); - - art::Handle pbtmcHandle; - event.getByLabel(pbtmcTag_, pbtmcHandle); - const ProtonBunchTimeMC& pbtmc(*pbtmcHandle); + const EventWindowMarker& ewMarker = *event.getValidHandle(ewMarkerToken_); + const ProtonBunchTimeMC& pbtmc = *event.getValidHandle(pbtmcToken_); ProditionsHandle eventTimingHandle; - const EventTiming &eventTiming = eventTimingHandle.get(event.id()); - float timeFromProtonsToDRMarker = eventTiming.timeFromProtonsToDRMarker(); //fixed time between CFO first tick and event start + const EventTiming& eventTiming = eventTimingHandle.get(event.id()); + const float timeFromProtonsToDRMarker = eventTiming.timeFromProtonsToDRMarker(); const auto& calCrystalConds = calCrystalConds_.get(event.id()); - // Containers to hold the output hits. - auto CaloShowerROs = std::make_unique(); + auto caloShowerROs = std::make_unique(); auto caloShowerSims = std::make_unique(); - StepHandles newCrystalShowerTokens; - std::transform(std::begin(crystalShowerTokens_), std::end(crystalShowerTokens_), - back_inserter(newCrystalShowerTokens), - [&event](const auto& token) {return event.getValidHandle(token);}); + StepHandles crystalShowerHandles; + crystalShowerHandles.reserve(crystalShowerTokens_.size()); + std::transform(crystalShowerTokens_.begin(), crystalShowerTokens_.end(), std::back_inserter(crystalShowerHandles), + [&event](const auto& token){return event.getValidHandle(token);}); - makeReadoutHits(newCrystalShowerTokens, *CaloShowerROs, *caloShowerSims, calCrystalConds, ewMarker, pbtmc, timeFromProtonsToDRMarker ); + makeReadoutHits(crystalShowerHandles, *caloShowerROs, *caloShowerSims, calCrystalConds, ewMarker, pbtmc, timeFromProtonsToDRMarker); - // Add the output hit collection to the event - event.put(std::move(CaloShowerROs)); + event.put(std::move(caloShowerROs)); event.put(std::move(caloShowerSims)); if (diagLevel_ > 0) std::cout << "[CaloShowerROMaker::produce] end" << std::endl; @@ -202,54 +188,49 @@ namespace mu2e { //----------------------------------------------------------------------------------------------------- - void CaloShowerROMaker::makeReadoutHits(const StepHandles& crystalShowerHandles, CaloShowerROCollection& CaloShowerROs, + void CaloShowerROMaker::makeReadoutHits(const StepHandles& crystalShowerHandles, CaloShowerROCollection& caloShowerROs, CaloShowerSimCollection& caloShowerSims, const CalSimParams& calCrystalConds, const EventWindowMarker& ewMarker, const ProtonBunchTimeMC& pbtmc, float timeFromProtonsToDRMarker) { - GlobalConstantsHandle pdt; - - float mbtime = GlobalConstantsHandle()->getNominalDRPeriod(); + GlobalConstantsHandle pdt; + const float mbtime = GlobalConstantsHandle()->getNominalDRPeriod(); const Calorimeter& cal = *(GeomHandle()); const float crystalLength = cal.G4Info().get("crystalZLength"); std::map> simEntriesMap; - diagSummary diagSum; + DiagSummary diagSum; - // Digitization start / end from accelerator DR marker with PB jitter - float correctedDigitizeStart = digitizationStart_ - pbtmc.pbtime_ - timeFromProtonsToDRMarker - digitizationBuffer_; - float correctedDigitizeEnd = digitizationEnd_ - pbtmc.pbtime_ - timeFromProtonsToDRMarker ; + // Digitization start / end from accelerator DR marker with PB jitter + const float correctedDigitizeStart = digitizationStart_ - pbtmc.pbtime_ - timeFromProtonsToDRMarker - digitizationBuffer_; + const float correctedDigitizeEnd = digitizationEnd_ - pbtmc.pbtime_ - timeFromProtonsToDRMarker; //----------------------------------------------------------------------- - //store corrected energy deposits for each redouts - for (const auto& showerHandle: crystalShowerHandles) + // store corrected energy deposits for each readout + for (const auto& showerHandle : crystalShowerHandles) { const CaloShowerStepCollection& caloShowerSteps(*showerHandle); - for (auto istep = caloShowerSteps.begin(); istep !=caloShowerSteps.end(); ++istep) + for (auto istep = caloShowerSteps.begin(); istep != caloShowerSteps.end(); ++istep) { const CaloShowerStep& step = *istep; - // see doc-db for calo folding description, this is non-trivial. Note pbtmc.pbtime_ is NEGATIVE! - // fold hits into the DR -> maxHitTime window - // then move early hits (in the PB) to the end of the digi window spilling over the next pulse - // finally filter hits to match the digitization window - double mbLength = (ewMarker.spillType() == EventWindowMarker::SpillType::onspill) ? mbtime : ewMarker.eventLength(); - double maxHitTime = std::max(mbtime,correctedDigitizeEnd); + // see doc-db for calo folding description. Note pbtmc.pbtime_ is NEGATIVE! + const double mbLength = (ewMarker.spillType() == EventWindowMarker::SpillType::onspill) ? mbtime : ewMarker.eventLength(); + const double maxHitTime = std::max(mbtime,correctedDigitizeEnd); - double hitTime = fmod(istep->time()+pbtmc.pbtime_,mbLength) - pbtmc.pbtime_; + double hitTime = std::fmod(step.time()+pbtmc.pbtime_, mbLength) - pbtmc.pbtime_; if (hitTime < maxHitTime-mbLength) hitTime += mbLength; - if (hitTime < correctedDigitizeStart) continue; - size_t idx = std::distance(caloShowerSteps.begin(), istep); - art::Ptr stepPtr = art::Ptr(showerHandle,idx); + const std::size_t idx = std::distance(caloShowerSteps.begin(), istep); + const auto stepPtr = art::Ptr(showerHandle,idx); - auto crystalID = CrystalId(step.volumeG4ID()); - int SiPMIDBase = crystalID.SiPMId(CaloConst::SiPM0); - float posZ = step.position().z(); - float lru = calCrystalConds.LRU(crystalID); - auto pePerMeVs = calCrystalConds.pePerMeVs(crystalID); + const auto crystalID = CrystalId(step.volumeG4ID()); + const int SiPMIDBase = crystalID.SiPMId(CaloConst::SiPM0); + const float posZ = step.position().z(); + const float lru = calCrystalConds.LRU(crystalID); + const auto pePerMeVs = calCrystalConds.pePerMeVs(crystalID); float edep_corr(step.energyDepG4()); if (BirksCorrection_) edep_corr = step.energyDepBirks(); @@ -258,60 +239,58 @@ namespace mu2e { // Generate individual PEs and their arrival times for (int i=0; i PETime(NPE,hitTime); if (addTravelTime_) + for (auto& t : PETime) t += photonProp_.propTimeSimu(crystalLength-posZ); + + if (diagLevel_ > 2) { - for (auto& time : PETime) time += photonProp_.propTimeSimu(crystalLength-posZ); + std::cout<<"[CaloShowerROMaker] SiPMID:"< 2) std::cout<<"[CaloShowerROMaker::generatePE] SiPMID:"< 2) {std::cout<<"Time hit "< newSteps = kv.second; - - // fill the summary map for each simPtr for a given crystalID std::map,SimParticleSummary> summaryMap; - for (auto& newStep : newSteps) + for (const auto& newStep : newSteps) { const art::Ptr& sim = newStep.step_->simParticle(); auto mfind = summaryMap.find(sim); if (mfind==summaryMap.end()) - summaryMap.insert(std::make_pair(sim,SimParticleSummary(newStep.step_,newStep.edepCorr_,newStep.timeCorr_))); + summaryMap.emplace(sim, SimParticleSummary(newStep.step_,newStep.edepCorr_,newStep.timeCorr_)); else mfind->second.add(newStep.step_,newStep.edepCorr_,newStep.timeCorr_); } - // create the CaloShowerSim (MC truth) objects for a given crystalID - for (auto& kvsumm : summaryMap) caloShowerSims.push_back(CaloShowerSim(kvsumm.second.steps_, kvsumm.second.edepCorr_, kvsumm.second.timeCorr_)); + for (auto& [sim, summ] : summaryMap) + caloShowerSims.emplace_back(summ.steps_, summ.edepCorr_, summ.timeCorr_); } - - //-------------------------------------------------- // Diag if (diagLevel_ > 2) dumpCaloShowerSim(caloShowerSims); @@ -321,49 +300,44 @@ namespace mu2e { std::set crIds; for (const auto& css : caloShowerSims) crIds.insert(css.crystalID()); - for (auto crId : crIds) + for (int crId : crIds) { - std::map, double> simMap; + std::map, double> simMap; for (const auto& css : caloShowerSims) if (css.crystalID()==crId) simMap[css.sim()] += css.energyDep(); - for (auto& kv : simMap) std::cout<<"CrId: "< 0) std::cout<<"[CaloShowerROMaker] found energy (energy corr) (edep_npe) / nStepsMC / nPE " < 2) std::cout<<"[CaloShowerROMaker::LRUCorrection] before / after LRU -> edep_corr = " << edepInit<<" / "< #include -#include -#include #include -#include +#include #include +#include #include -#include -#include - +#include namespace { + // Accumulates the steps of one ancestor SimParticle together with the set of + // (secondary) SimParticles that were folded into it. class CaloCompressUtil { public: - CaloCompressUtil() : steps_(), sims_() {} + using SimPtr = art::Ptr; const std::vector& steps() const {return steps_;} - const std::set>& sims() const {return sims_;} + const std::set& sims() const {return sims_;} - void fill(const mu2e::StepPointMC* step, std::vector> sims) + void fill(const mu2e::StepPointMC* step, const std::vector& sims) { steps_.push_back(step); - for (const auto& sim: sims) sims_.insert(sim); + sims_.insert(sims.begin(), sims.end()); } private: std::vector steps_; - std::set> sims_; + std::set sims_; }; - struct diagSummary + struct DiagSummary { - diagSummary() : totalEdep_(0.0),totalStep_(0),totalSim_(0),totalChk_(0),nCompress_(0),ncompressInfo_(0) {}; - void reset() {totalEdep_=0.0;totalStep_=totalSim_=totalChk_=nCompress_=ncompressInfo_=0;} - - float totalEdep_; - unsigned totalStep_,totalSim_,totalChk_,nCompress_,ncompressInfo_; + float totalEdep_{0.f}; + unsigned totalStep_{0}, totalSim_{0}, totalChk_{0}, nCompress_{0}, nCompressInfo_{0}; + void reset() {*this = DiagSummary{};} }; } - namespace mu2e { class CaloShowerStepMaker : public art::EDProducer @@ -93,142 +79,72 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - fhicl::Sequence caloStepPointCollection { Name("caloStepPointCollection"), Comment("Calo crystal stepPointMC collection name") }; - fhicl::Atom physVolInfoInput { Name("physVolInfoInput"), Comment("Physics volume token names") }; - fhicl::Atom numZSlices { Name("numZSlices"), Comment("Number of crystal longitudinal slices ") }; - fhicl::Atom deltaTime { Name("deltaTime"), Comment("Max time difference to be inside a ShowerStep") }; - fhicl::Sequence caloMaterial { Name("caloMaterial"), Comment("List of calo material names") }; - fhicl::Atom compressData { Name("compressData"), Comment("Compress stepPointMC and SimParticles in crystal") }; - fhicl::Atom eDepThreshold { Name("eDepThreshold"), Comment("Threshold on energy deposited by SimParticle to keep it") }; - fhicl::Atom diagLevel { Name("diagLevel"), Comment("Debug"),0 }; + fhicl::Sequence caloStepPointCollection { Name("caloStepPointCollection"), Comment("Calo crystal stepPointMC collection name") }; + fhicl::Atom numZSlices { Name("numZSlices"), Comment("Number of crystal longitudinal slices") }; + fhicl::Atom deltaTime { Name("deltaTime"), Comment("Max time difference to be inside a ShowerStep") }; + fhicl::Atom compressData { Name("compressData"), Comment("Compress stepPointMC and SimParticles in crystal") }; + fhicl::Atom eDepThreshold { Name("eDepThreshold"), Comment("Threshold on energy deposited by SimParticle to keep it") }; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Debug"),0 }; }; explicit CaloShowerStepMaker(const art::EDProducer::Table& config); - - void beginJob() override; - void beginSubRun(art::SubRun& sr) override; - void produce( art::Event& e) override; + void produce(art::Event& e) override; private: - using HandleVector = std::vector>; - using SimPtr = art::Ptr; - using SimStepMap = std::map>; - - void makeCompressedHits (const HandleVector&, CaloShowerStepCollection&, SimParticlePtrCollection&); - void collectStepBySimAncestor (const Calorimeter&, const PhysicalVolumeMultiHelper&, const HandleVector&, std::map&); - void collectStepBySim (const HandleVector&, SimStepMap&); - void compressSteps (const Calorimeter&, CaloShowerStepCollection&, int, const SimPtr&, std::vector&); - void fillHisto1 (const Calorimeter&, const SimPtr&, const std::set&); - void fillHisto2 (int, float, const SimPtr&); - void dumpAllInfo (const HandleVector&, const Calorimeter&); - - - std::vector calorimeterStepPoints_; - art::InputTag physVolInfoInput_; - std::set mapPhysVol_; - std::vector caloMaterial_; - int numZSlices_; - double deltaTime_; - bool compressData_; - double eDepThreshold_; - int diagLevel_; - const PhysicalVolumeInfoMultiCollection* vols_ = nullptr; - double zSliceSize_; - - diagSummary diagSummary_; - TH2F* hStartPos_; - TH2F* hStopPos_; - TH1F* hStopPos2_; - TH1F* hStartPos2_; - TH1F* hZpos_; - TH1F* hEtot_; - TH1F* hStot_; - TH2F* hZpos2_; - TH1F* hGenId_; + using SimPtr = art::Ptr; + using StepHandles = std::vector>; + + void makeCompressedHits (const StepHandles&, CaloShowerStepCollection&, SimParticlePtrCollection&); + void collectStepBySimAncestor(const Calorimeter&, const StepHandles&, std::map&); + void compressSteps (const Calorimeter&, CaloShowerStepCollection&, int volId, const SimPtr&, std::vector& steps); + void dumpAllInfo (const StepHandles&, const Calorimeter&) const; + + std::vector> stepTokens_; + unsigned numZSlices_; + double deltaTime_; + bool compressData_; + double eDepThreshold_; + int diagLevel_; + double zSliceSize_{0}; + DiagSummary diag_; }; CaloShowerStepMaker::CaloShowerStepMaker(const art::EDProducer::Table& config) : art::EDProducer{config}, - calorimeterStepPoints_(config().caloStepPointCollection()), - physVolInfoInput_ (config().physVolInfoInput()), - caloMaterial_ (config().caloMaterial()), - numZSlices_ (config().numZSlices()), - deltaTime_ (config().deltaTime()), - compressData_ (config().compressData()), - eDepThreshold_ (config().eDepThreshold()), - diagLevel_ (config().diagLevel()), - vols_(), - zSliceSize_(0), - diagSummary_() - { - consumesMany(); - produces(); - produces(); - } - - - //-------------------------------------------------------------------- - void CaloShowerStepMaker::beginJob() - { - if (diagLevel_ > 1) - { - art::ServiceHandle tfs; - hStartPos_ = tfs->make("hStartPos", "Sim start position", 1000, 5000, 15000, 200, 0, 1000); - hStopPos_ = tfs->make("hStopPos", "Sim stop position", 1000, 5000, 15000, 200, 0, 1000); - hStartPos2_ = tfs->make("hStartPos2","Sim start position", 1000, 10000, 13000); - hStopPos2_ = tfs->make("hStopPos2", "Sim stop position", 1000, 10000, 13000); - hZpos_ = tfs->make("hZpos", "Step z pos", 20, 0, 20); - hZpos2_ = tfs->make("hZpos2", "Step z pos", 20, 0, 20, 100, 0, 5); - hEtot_ = tfs->make("hEtot", "Total E dep", 150, 0, 150); - hStot_ = tfs->make("hStot", "Total numebr steps", 100, 0, 10000); - hGenId_ = tfs->make("hSimId", "Gen Id", 150, -10, 140); - } - } - - - - void CaloShowerStepMaker::beginSubRun(art::SubRun& sr) + numZSlices_ (config().numZSlices()), + deltaTime_ (config().deltaTime()), + compressData_ (config().compressData()), + eDepThreshold_(config().eDepThreshold()), + diagLevel_ (config().diagLevel()) { - mapPhysVol_.clear(); + if (numZSlices_ == 0) + throw cet::exception("CONFIG") << "[CaloShowerStepMaker] numZSlices must be > 0\n"; - art::Handle volh; - sr.getByLabel(physVolInfoInput_, volh); - if (!volh.isValid()) return; + for (const auto& tag : config().caloStepPointCollection()) + stepTokens_.push_back(consumes(art::InputTag(tag))); - vols_ = volh.product(); - for (const auto& vol : *volh) - { - for (const auto& mv : vol) - { - if (std::find(caloMaterial_.begin(),caloMaterial_.end(), mv.second.materialName()) != caloMaterial_.end()) - mapPhysVol_.insert(&mv.second); - } - } + produces(); + produces(); } //------------------------------------------------------------------------------------------------------------ - // void CaloShowerStepMaker::produce(art::Event& event) { - diagSummary_.reset(); + diag_.reset(); if (diagLevel_ > 0) std::cout << "[CaloShowerStepMaker::produce] begin" << std::endl; auto caloShowerStepMCs = std::make_unique(); auto simsToKeep = std::make_unique(); - HandleVector crystalStepsHandles; - for (const auto& stepPts : calorimeterStepPoints_) - { - art::Handle hc; - event.getByLabel(art::InputTag(stepPts), hc); - crystalStepsHandles.push_back(hc); - } + StepHandles crystalStepsHandles; + crystalStepsHandles.reserve(stepTokens_.size()); + for (const auto& token : stepTokens_) crystalStepsHandles.push_back(event.getValidHandle(token)); - makeCompressedHits(crystalStepsHandles,*caloShowerStepMCs,*simsToKeep); + makeCompressedHits(crystalStepsHandles, *caloShowerStepMCs, *simsToKeep); event.put(std::move(caloShowerStepMCs)); event.put(std::move(simsToKeep)); @@ -238,173 +154,146 @@ namespace mu2e { //------------------------------------------------------------------------------------------------------------------ - void CaloShowerStepMaker::makeCompressedHits(const HandleVector& crystalStepsHandle, - CaloShowerStepCollection& caloShowerStepMCs,SimParticlePtrCollection& simsToKeep) + void CaloShowerStepMaker::makeCompressedHits(const StepHandles& crystalStepsHandles, + CaloShowerStepCollection& caloShowerStepMCs, SimParticlePtrCollection& simsToKeep) { - PhysicalVolumeMultiHelper vi(vols_); - const Calorimeter& cal = *(GeomHandle()); - zSliceSize_ = cal.G4Info().get("crystalZLength")/float(numZSlices_)+1e-5; - + zSliceSize_ = cal.G4Info().get("crystalZLength")/float(numZSlices_) + 1e-5; - //----------------------------------------------------------------- - // Collect the StepPointMC's produced by each SimParticle Ancestor + // Collect the StepPointMCs produced by each ancestor SimParticle std::map crystalAncestorsMap; - collectStepBySimAncestor(cal,vi,crystalStepsHandle,crystalAncestorsMap); + collectStepBySimAncestor(cal, crystalStepsHandles, crystalAncestorsMap); - if (diagLevel_ > 2) dumpAllInfo(crystalStepsHandle,cal); + if (diagLevel_ > 2) dumpAllInfo(crystalStepsHandles, cal); - - //--------------------------------------------------------------------------------------------------------------- - //Loop over ancestor simParticles, check if they are compressible, and produce the corresponding caloShowerStepMC - - std::set SimsToKeepUnique; - for (const auto& iter : crystalAncestorsMap ) + // Loop over ancestor SimParticles, check if they are compressible, and produce the corresponding CaloShowerStep + std::set simsToKeepUnique; + for (const auto& [sim, info] : crystalAncestorsMap) { - const SimPtr& sim = iter.first; - const CaloCompressUtil& info = iter.second; - - diagSummary_.totalSim_ += info.sims().size(); + diag_.totalSim_ += info.sims().size(); std::map> crystalMap; for (const StepPointMC* step : info.steps()) crystalMap[step->volumeId()].push_back(step); - for (const auto& iterCrystal : crystalMap) + for (auto& [crid, steps] : crystalMap) { - unsigned crid = iterCrystal.first; - std::vector steps = iterCrystal.second; - //Filter very small energy deposits at this stage double eDep(0); - for (const auto& step : steps) eDep += step->totalEDep(); + for (const auto* step : steps) eDep += step->totalEDep(); if (eDep < eDepThreshold_) continue; if (compressData_) { - SimsToKeepUnique.insert(sim); + simsToKeepUnique.insert(sim); compressSteps(cal, caloShowerStepMCs, crid, sim, steps); - if (diagLevel_ > 1) fillHisto1(cal,sim,info.sims()); } else { std::map> newSimStepMap; for (const StepPointMC* step : steps) newSimStepMap[step->simParticle()].push_back(step); - for (auto& iter : newSimStepMap) + for (auto& [stepSim, stepVec] : newSimStepMap) { - compressSteps(cal, caloShowerStepMCs, crid, iter.first, iter.second); - SimsToKeepUnique.insert(iter.first); + compressSteps(cal, caloShowerStepMCs, crid, stepSim, stepVec); + simsToKeepUnique.insert(stepSim); } } } - ++diagSummary_.ncompressInfo_; - if (compressData_) ++diagSummary_.nCompress_; + ++diag_.nCompressInfo_; + if (compressData_) ++diag_.nCompress_; } - //dump the unique set of SimParticles to keep into final vector - simsToKeep.assign(SimsToKeepUnique.begin(),SimsToKeepUnique.end()); + simsToKeep.assign(simsToKeepUnique.begin(), simsToKeepUnique.end()); - //--------------------------------------------------------------------------------------------------------------- // Final diag info if (diagLevel_ > 1) { - hEtot_->Fill(diagSummary_.totalEdep_); - hStot_->Fill(diagSummary_.totalStep_); std::cout<<"CaloShowerStepMaker summary"< volIds{}; - for (auto caloShowerStepMC : caloShowerStepMCs) volIds.insert(caloShowerStepMC.volumeG4ID()); + std::set volIds; + for (const auto& css : caloShowerStepMCs) volIds.insert(css.volumeG4ID()); - for (auto volId: volIds) + for (int volId : volIds) { - std::map, double> simMap; - for (const auto& caloShowerStepMC : caloShowerStepMCs) - if (caloShowerStepMC.volumeG4ID()==volId) simMap[caloShowerStepMC.simParticle()] += caloShowerStepMC.energyDepG4(); + std::map, double> simMap; + for (const auto& css : caloShowerStepMCs) + if (css.volumeG4ID()==volId) simMap[css.simParticle()] += css.energyDepG4(); - for (auto& kv : simMap) std::cout<<"Vol id: "<& ancestorsMap) + void CaloShowerStepMaker::collectStepBySimAncestor(const Calorimeter& cal, + const StepHandles& stepsHandles, + std::map& ancestorsMap) { + std::vector inspectedSims; std::unordered_map simToAncestorMap; - for (HandleVector::const_iterator i=stepsHandles.begin(), e=stepsHandles.end(); i != e; ++i ) + + for (const auto& handle : stepsHandles) { - const art::Handle& handle(*i); const StepPointMCCollection& steps(*handle); - for (const auto& step : steps ) + + for (const auto& step : steps) { SimPtr sim = step.simParticle(); - SimParticlePtrCollection inspectedSims; - while (sim->hasParent() && cal.isInsideAnyCrystal(sim->startPosition()) ) + inspectedSims.clear(); + while (sim->hasParent()) { - //simparticle starting in one section and ending in another one see note above - if (!cal.isInsideSameDisk(sim->startPosition(),sim->endPosition()) ) break; - const auto alreadyInspected = simToAncestorMap.find(sim); if (alreadyInspected != simToAncestorMap.end()) {sim = alreadyInspected->second; break;} - inspectedSims.push_back(sim); + + if (!cal.isInsideAnyCrystal(sim->startPosition())) break; + if (!cal.isInsideSameDisk(sim->startPosition(),sim->endPosition())) break; + sim = sim->parent(); } for (const SimPtr& inspectedSim : inspectedSims) simToAncestorMap[inspectedSim] = sim; - ancestorsMap[sim].fill(&step,inspectedSims); + ancestorsMap[sim].fill(&step, inspectedSims); - diagSummary_.totalEdep_ += step.totalEDep(); + diag_.totalEdep_ += step.totalEDep(); } - diagSummary_.totalStep_ += steps.size(); - } - } - - - //----------------------------------------------------------------------------------------------------------------------------------------------- - void CaloShowerStepMaker::collectStepBySim(const HandleVector& stepsHandles, - std::map>& simStepMap) - { - for (HandleVector::const_iterator i=stepsHandles.begin(), e=stepsHandles.end(); i != e; ++i) - { - const art::Handle& handle(*i); - const StepPointMCCollection& steps(*handle); - for (const auto& step : steps ) simStepMap[step.simParticle()].push_back(&step); + diag_.totalStep_ += steps.size(); } } //------------------------------------------------------------------------------------------------------------------------------- - void CaloShowerStepMaker::compressSteps(const Calorimeter& cal, CaloShowerStepCollection &caloShowerStepMCs, + void CaloShowerStepMaker::compressSteps(const Calorimeter& cal, CaloShowerStepCollection& caloShowerStepMCs, int volId, const SimPtr& sim, std::vector& steps) { - auto sortFunctor = [](const StepPointMC* a, const StepPointMC* b) {return a->time() < b->time();}; - std::sort(steps.begin(), steps.end(), sortFunctor); + std::sort(steps.begin(), steps.end(), [](const StepPointMC* a, const StepPointMC* b){return a->time() < b->time();}); - ShowerStepUtil buffer(numZSlices_, ShowerStepUtil::weight_type::energy ); + ShowerStepUtil buffer(numZSlices_, ShowerStepUtil::weight_type::energy); for (const StepPointMC* step : steps) { - CLHEP::Hep3Vector pos = cal.mu2eToCrystal(volId,step->position()); - int idx = int(std::max(1e-6,pos.z())/zSliceSize_); + const CLHEP::Hep3Vector pos = cal.mu2eToCrystal(volId, step->position()); - if (buffer.entries(idx)>0 && (step->time()-buffer.t0(idx) > deltaTime_) ) + // clamp: a step at/beyond the back face would otherwise index past the last slice + const unsigned idx = std::min(unsigned(std::max(1e-6, pos.z())/zSliceSize_), numZSlices_-1); + + if (buffer.entries(idx)>0 && (step->time()-buffer.t0(idx) > deltaTime_)) { - if (diagLevel_ > 1) {fillHisto2(idx,buffer.energyG4(idx),sim);} if (diagLevel_ > 2) {std::cout<<"[CaloShowerStepMaker::compressSteps] inserted "; buffer.printBucket(idx);} - diagSummary_.totalChk_ += buffer.entries(idx); + diag_.totalChk_ += buffer.entries(idx); - caloShowerStepMCs.push_back(CaloShowerStep(volId, sim, buffer.entries(idx), buffer.time(idx), buffer.energyG4(idx), - buffer.energyVis(idx),buffer.pIn(idx),buffer.pos(idx))); + caloShowerStepMCs.emplace_back(volId, sim, buffer.entries(idx), buffer.time(idx), buffer.energyG4(idx), + buffer.energyVis(idx), buffer.pIn(idx), buffer.pos(idx)); buffer.reset(idx); } @@ -412,55 +301,27 @@ namespace mu2e { } //do not forget to flush the final buffer(s) :-) - for (unsigned i=0;i 1) {fillHisto2(i,buffer.energyG4(i),sim);} - if (diagLevel_ > 2) {std::cout<<"[CaloShowerStepMaker::compressSteps] inserted "; buffer.printBucket(i);} - diagSummary_.totalChk_ += buffer.entries(i); + if (diagLevel_ > 2) {std::cout<<"[CaloShowerStepMaker::compressSteps] inserted "; buffer.printBucket(i);} + diag_.totalChk_ += buffer.entries(i); - caloShowerStepMCs.push_back(CaloShowerStep(volId, sim, buffer.entries(i), buffer.time(i), buffer.energyG4(i), - buffer.energyVis(i),buffer.pIn(i),buffer.pos(i))); + caloShowerStepMCs.emplace_back(volId, sim, buffer.entries(i), buffer.time(i), buffer.energyG4(i), + buffer.energyVis(i), buffer.pIn(i), buffer.pos(i)); } } //------------------------------------------------------------------------------------------------------------- - void CaloShowerStepMaker::fillHisto1(const Calorimeter& cal, const art::Ptr& sim, const std::set>& infoSims) - { - CLHEP::Hep3Vector startSection = cal.mu2eToDisk(0,sim->startPosition()); - CLHEP::Hep3Vector endSection = cal.mu2eToDisk(0,sim->endPosition()); - double rStart = sqrt(startSection.x()*startSection.x()+startSection.y()*startSection.y()); - double rEnd = sqrt(endSection.x()*endSection.x()+endSection.y()*endSection.y()); - - hStartPos_->Fill(sim->startPosition().z(),rStart); - hStopPos_->Fill( sim->endPosition().z(), rEnd); - for (const auto& simD: infoSims) - { - hStartPos2_->Fill(simD->startPosition().z()); - hStopPos2_->Fill(simD->endPosition().z()); - } - } - - //------------------------------------------------------------------------------------------------------------- - void CaloShowerStepMaker::fillHisto2(int idx, float edep, const SimPtr& sim) - { - hZpos_->Fill(idx); - hZpos2_->Fill(idx,edep); - if (sim->genParticle()) hGenId_->Fill(sim->genParticle()->generatorId().id()); - } - - //------------------------------------------------------------------------------------------------------------- - void CaloShowerStepMaker::dumpAllInfo(const HandleVector& stepsHandles, const Calorimeter& cal) + void CaloShowerStepMaker::dumpAllInfo(const StepHandles& stepsHandles, const Calorimeter& cal) const { std::cout<<"Dumping StepPointMCs Mu2e / crystal / disk / diskFF frames"<& handle(*i); const StepPointMCCollection& steps(*handle); - std::cout< Date: Wed, 5 Aug 2026 00:12:45 -0500 Subject: [PATCH 2/7] More refactoring --- CaloMC/inc/CaloNoiseSimGenerator.hh | 82 ---------- CaloMC/inc/CaloPhotonPropagation.hh | 41 ++--- CaloMC/inc/CaloWFExtractor.hh | 28 ---- CaloMC/src/CaloClusterTruthMatch_module.cc | 85 +++++----- CaloMC/src/CaloDigiMaker_module.cc | 155 +++++++++--------- CaloMC/src/CaloNoiseSimGenerator.cc | 178 --------------------- CaloMC/src/CaloPhotonPropagation.cc | 99 ++++++------ CaloMC/src/CaloWFExtractor.cc | 47 ------ 8 files changed, 198 insertions(+), 517 deletions(-) delete mode 100644 CaloMC/inc/CaloNoiseSimGenerator.hh delete mode 100644 CaloMC/inc/CaloWFExtractor.hh delete mode 100644 CaloMC/src/CaloNoiseSimGenerator.cc delete mode 100644 CaloMC/src/CaloWFExtractor.cc diff --git a/CaloMC/inc/CaloNoiseSimGenerator.hh b/CaloMC/inc/CaloNoiseSimGenerator.hh deleted file mode 100644 index 4d66cd59cb..0000000000 --- a/CaloMC/inc/CaloNoiseSimGenerator.hh +++ /dev/null @@ -1,82 +0,0 @@ -#ifndef CaloNoiseSimGenerator_HH -#define CaloNoiseSimGenerator_HH -// -// Generate long noise waveform to use for calorimeter digitization -// -#include "fhiclcpp/types/Atom.h" -#include "fhiclcpp/types/Sequence.h" -#include "art/Framework/Services/Optional/RandomNumberGenerator.h" -#include "Offline/SeedService/inc/SeedService.hh" - -#include "Offline/CaloMC/inc/CaloWFExtractor.hh" -#include "Offline/Mu2eUtilities/inc/CaloPulseShape.hh" - -#include "CLHEP/Random/RandPoissonQ.h" -#include "CLHEP/Random/RandGaussQ.h" -#include "CLHEP/Random/RandFlat.h" - - -namespace mu2e { - - class CaloNoiseSimGenerator - { - public: - struct Config - { - using Name = fhicl::Name; - using Comment = fhicl::Comment; - fhicl::Atom pulseFileName { Name("pulseFileName"), Comment("Calo pulse file name") }; - fhicl::Atom pulseHistName { Name("pulseHistName"), Comment("Calo pulse hist name") }; - fhicl::Atom elecNphotPerNs { Name("elecNphotPerNs"), Comment("Electronics noise number of PE / ns ") }; - fhicl::Atom rinNphotPerNs { Name("rinNphotPerNs"), Comment("RIN noise number of PE / ns ") }; - fhicl::Atom darkNphotPerNs { Name("darkNphotPerNs"), Comment("SiPM Dark noise number of PE / ns ") }; - fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; - fhicl::Atom pePerMeV { Name("readoutPEPerMeV"),Comment("Number of pe / MeV for Readout") }; - fhicl::Atom MeVToADC { Name("MeVToADC"), Comment("MeV to ADC conversion factor") }; - fhicl::Atom noiseWFSize { Name("noiseWFSize"), Comment("Noise WF size") }; - fhicl::Atom nMaxFragment { Name("nMaxFragment"), Comment("maximum number of wf generated for extracting noise fragments ") }; - fhicl::Atom minPeakADC { Name("minPeakADC"), Comment("Minimum ADC hits of local peak to digitize") }; - fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; - }; - - - CaloNoiseSimGenerator(const Config& config, CLHEP::HepRandomEngine& engine, int iRO); - - void initialize(const CaloWFExtractor& wfExtractor); - void refresh(); - - void addSampleNoise(std::vector& wfVector, unsigned istart, unsigned ilength); - void addSaltAndPepper(std::vector& wfVector); - void plotNoise(const std::string& name); - - const std::vector& noise() const {return waveform_;} - double pedestal() const {return pedestal_;} - - - private: - using vvd = std::vector>; - - void generateWF(std::vector& wfVector); - void generateFragments(const CaloWFExtractor& wfExtractor); - - unsigned iRO_; - std::vector waveform_; - int pedestal_; - vvd digiNoise_; - double digiNoiseProb_; - double digiSampling_; - double noiseRinDark_; - double noiseElec_; - double minPeakADC_; - double pePerMeV_; - double MeVToADC_; - CLHEP::RandPoissonQ randPoisson_; - CLHEP::RandGaussQ randGauss_; - CLHEP::RandFlat randFlat_; - unsigned nMaxFragment_; - CaloPulseShape pulseShape_; - int diagLevel_; - }; - -} -#endif diff --git a/CaloMC/inc/CaloPhotonPropagation.hh b/CaloMC/inc/CaloPhotonPropagation.hh index 8bd4c55077..104fe26f86 100644 --- a/CaloMC/inc/CaloPhotonPropagation.hh +++ b/CaloMC/inc/CaloPhotonPropagation.hh @@ -1,11 +1,14 @@ -#ifndef CaloPhotonPropagation_HH -#define CaloPhotonPropagation_HH - -// Calculate the propagation time from the location in the crystal -// Input based on detail Geant4 simulation of crystal - +#ifndef CaloMC_CaloPhotonPropagation_hh +#define CaloMC_CaloPhotonPropagation_hh +// +// Sample the scintillation-photon propagation time from a location in the crystal. +// The time distribution vs. depth is taken from a detailed Geant4 simulation, stored +// as a 2D histogram (z on X, propagation time on Y) and turned into a per-depth CDF. +// #include "CLHEP/Random/RandomEngine.h" #include "CLHEP/Random/RandFlat.h" + +#include #include namespace mu2e { @@ -15,20 +18,20 @@ namespace mu2e { public: CaloPhotonPropagation(const std::string& fileName, const std::string& histName, CLHEP::HepRandomEngine& engine); - void buildTable (); - float propTimeSimu(float z); - float propTimeLine(float z); + void buildTable(); + float propTimeSimu(float z); // sampled from the CDF (consumes a random number) + float propTimeLine(float z) const; // straight-line n*z/c estimate - private: - std::vector timeProp_; - std::vector cdf_; - unsigned nTimeDiv_; - unsigned nZDiv_; - float dzTime_; - CLHEP::RandFlat randFlat_; - std::string fileName_; - std::string histName_; - float lightSpeed_; + private: + std::vector timeProp_; // Y-bin centers: candidate propagation times + std::vector cdf_; // per-depth cumulative distribution, row-major [iz][itime] + unsigned nTimeDiv_{0}; // number of time bins (histogram Y) + unsigned nZDiv_{0}; // number of depth bins (histogram X) + float dz_{0.f}; // depth-bin width + CLHEP::RandFlat randFlat_; + std::string fileName_; + std::string histName_; + float lightSpeed_{300.f}; // mm/ns, overwritten with c/n in buildTable }; } diff --git a/CaloMC/inc/CaloWFExtractor.hh b/CaloMC/inc/CaloWFExtractor.hh deleted file mode 100644 index ad65606b4b..0000000000 --- a/CaloMC/inc/CaloWFExtractor.hh +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef CaloMC_CaloWFExtractor_hh -#define CaloMC_CaloWFExtractor_hh -// -// Utility to simulate waveform hit extraction in FPGA -// -#include -#include - -namespace mu2e { - - class CaloWFExtractor - { - public: - CaloWFExtractor(unsigned bufferDigi, unsigned nBinsPeak, int minPeakADC, unsigned startOffset) : - bufferDigi_(bufferDigi),nBinsPeak_(nBinsPeak),minPeakADC_(minPeakADC), startOffset_(startOffset) - {}; - - void extract(const std::vector& wf, std::vector& starts, std::vector& stops) const; - - private: - unsigned bufferDigi_; - unsigned nBinsPeak_; - int minPeakADC_; - unsigned startOffset_; - }; -} - -#endif diff --git a/CaloMC/src/CaloClusterTruthMatch_module.cc b/CaloMC/src/CaloClusterTruthMatch_module.cc index a9ad2b339a..599265b011 100644 --- a/CaloMC/src/CaloClusterTruthMatch_module.cc +++ b/CaloMC/src/CaloClusterTruthMatch_module.cc @@ -1,9 +1,9 @@ // // An EDProducer Module to match calo clusters to MC info // -// #include "art/Framework/Core/EDProducer.h" #include "art/Framework/Principal/Event.h" +#include "fhiclcpp/types/Atom.h" #include "Offline/MCDataProducts/inc/CaloEDepMC.hh" #include "Offline/MCDataProducts/inc/SimParticle.hh" @@ -11,16 +11,17 @@ #include "Offline/MCDataProducts/inc/CaloHitMC.hh" #include "Offline/MCDataProducts/inc/CaloClusterMC.hh" #include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/RecoDataProducts/inc/CaloHit.hh" +#include #include -#include #include +#include #include namespace mu2e { - class CaloClusterTruthMatch : public art::EDProducer { public: @@ -47,85 +48,81 @@ namespace mu2e { private: - void makeTruthMatch(art::Event&, CaloClusterMCCollection&,CaloClusterMCTruthAssn&); + void makeTruthMatch(art::Event&, CaloClusterMCCollection&, CaloClusterMCTruthAssn&) const; const art::ProductToken caloClusterToken_; const art::ProductToken caloHitMCTruthToken_; - int diagLevel_; + const int diagLevel_; }; - - //-------------------------------------------------------------------- void CaloClusterTruthMatch::produce(art::Event& event) { - std::unique_ptr caloClusterMCs(new CaloClusterMCCollection); - std::unique_ptr CaloClusterMCTruth(new CaloClusterMCTruthAssn); + auto caloClusterMCs = std::make_unique(); + auto caloClusterMCTruth = std::make_unique(); - makeTruthMatch(event, *caloClusterMCs, *CaloClusterMCTruth); + makeTruthMatch(event, *caloClusterMCs, *caloClusterMCTruth); - event.put(std::move(CaloClusterMCTruth)); + event.put(std::move(caloClusterMCTruth)); event.put(std::move(caloClusterMCs)); } //-------------------------------------------------------------------- void CaloClusterTruthMatch::makeTruthMatch(art::Event& event, CaloClusterMCCollection& caloClusterMCs, - CaloClusterMCTruthAssn& caloClusterTruthMatch) + CaloClusterMCTruthAssn& caloClusterTruthMatch) const { - - art::ProductID clusterMCProductID(event.getProductID()); + const art::ProductID clusterMCProductID(event.getProductID()); const art::EDProductGetter* clusterMCProductGetter = event.productGetter(clusterMCProductID); const auto caloClusterHandle = event.getValidHandle(caloClusterToken_); const auto& caloClusters(*caloClusterHandle); - const auto* caloClusterBase = caloClusters.data(); + const auto& caloHitTruth = *event.getValidHandle(caloHitMCTruthToken_); - const auto CaloHitMCHandle = event.getValidHandle(caloHitMCTruthToken_); - const auto& caloHitTruth(*CaloHitMCHandle); + // build the hit -> hitMC lookup once, instead of re-scanning the whole association per cluster + std::map, art::Ptr> hitToMC; + for (const auto& assn : caloHitTruth) hitToMC.emplace(assn.first, assn.second); double totalEnergyMatched(0); - int nMatched(0); + int nMatched(0); + caloClusterMCs.reserve(caloClusters.size()); - for (const auto& cluster : caloClusters) + for (std::size_t idx=0; idx clusterPtr = art::Ptr(caloClusterHandle,idx); - const auto& hits = cluster.caloHitsPtrVector(); + const CaloCluster& cluster = caloClusters[idx]; + const auto clusterPtr = art::Ptr(caloClusterHandle,idx); - if (diagLevel_ > 1) std::cout<<"[CaloClusterTruthMatch] Inspect cluster diskId/energy/time "<> digis; + if (diagLevel_ > 1) std::cout<<"[CaloClusterTruthMatch] Inspect cluster diskId/energy/time " + <first) == hits.end()) continue; - const auto& digiMC = i->second; - digis.push_back(digiMC); + // gather the CaloHitMC of every hit in this cluster that carries MC truth + std::vector> digis; + for (const auto& hitPtr : cluster.caloHitsPtrVector()) + { + const auto it = hitToMC.find(hitPtr); + if (it == hitToMC.end()) continue; - if (diagLevel_ > 1 && digiMC->nParticles()>0) std::cout<<"[CaloClusterTruthMatch] found hit in map "<nParticles()<<" "<time()<second); + if (diagLevel_ > 1 && it->second->nParticles()>0) + std::cout<<"[CaloClusterTruthMatch] found hit in map "<second->nParticles()<<" "<second->time()<totalEnergyDep() > b->totalEnergyDep();}); - caloClusterMCs.emplace_back(CaloClusterMC(std::move(digis))); + std::sort(digis.begin(),digis.end(),[](const auto& a, const auto& b){return a->totalEnergyDep() > b->totalEnergyDep();}); + caloClusterMCs.emplace_back(std::move(digis)); - art::Ptr clusterMCPtr = art::Ptr(clusterMCProductID, caloClusterMCs.size()-1, clusterMCProductGetter); - caloClusterTruthMatch.addSingle(clusterPtr,clusterMCPtr); + const auto clusterMCPtr = art::Ptr(clusterMCProductID, caloClusterMCs.size()-1, clusterMCProductGetter); + caloClusterTruthMatch.addSingle(clusterPtr,clusterMCPtr); - totalEnergyMatched += clusterPtr->energyDep(); - ++nMatched; + totalEnergyMatched += cluster.energyDep(); + ++nMatched; } if (diagLevel_ > 0) std::cout<<"[CaloClusterTruthMatch] total clusters / energy matched = "< #include #include @@ -55,7 +47,7 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - using CNG = mu2e::CaloNoiseSimGenerator::Config; + using CNG = mu2e::CaloNoiseGenerator::Config; fhicl::Table noise_gen_conf { Name("NoiseGenerator"), Comment("Noise generator config") }; fhicl::Atom caloShowerCollection { Name("caloShowerROCollection"), Comment("CaloShowerRO collection name") }; fhicl::Atom ewMarkerTag { Name("eventWindowMarker"), Comment("EventWindowMarker producer") }; @@ -67,7 +59,9 @@ namespace mu2e { fhicl::Atom addNoise { Name("addNoise"), Comment("Add noise to waveform") }; fhicl::Atom addRandomNoise { Name("addRandomNoise"), Comment("Add random salt and pepper noise") }; fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; - fhicl::Atom minNoiseAmplitude { Name("minNoiseAmplitude"), Comment("minimum amplitude to start adding noise") }; + fhicl::Atom pePerMeVCsI { Name("readoutPEPerMeVCsI"), Comment("Number of pe / MeV for Readout for CsI") }; + fhicl::Atom pePerMeVLyso { Name("readoutPEPerMeVLyso"), Comment("Number of pe / MeV for Readout for LYSO") }; + fhicl::Atom MeVToADC { Name("MeVToADC"), Comment("MeV to ADC conversion factor") }; fhicl::Atom nBits { Name("nBits"), Comment("ADC Number of bits") }; fhicl::Atom nBinsPeak { Name("nBinsPeak"), Comment("Window size for finding local maximum to digitize wf") }; fhicl::Atom minPeakADC { Name("minPeakADC"), Comment("Minimum ADC hits of local peak to digitize") }; @@ -86,12 +80,15 @@ namespace mu2e { bufferDigi_ (config().bufferDigi()), startTimeBuffer_ (config().digiSampling()*config().bufferDigi()), maxADCCounts_ ((1 << config().nBits()) -1), - minNoiseAmplitude_ (config().minNoiseAmplitude()), + pePerMeVCsI_ (config().pePerMeVCsI()), + pePerMeVLyso_ (config().pePerMeVLyso()), + MeVToADC_ (config().MeVToADC()), pulseShape_ (CaloPulseShape(config().pulseFileName(),config().pulseHistName(),config().digiSampling())), - wfExtractor_ (config().bufferDigi(),config().nBinsPeak(),config().minPeakADC(),config().bufferDigi()), + nBinsPeak_ (config().nBinsPeak()), + minPeakADC_ (config().minPeakADC()), engine_ (createEngine(art::ServiceHandle()->getSeed())), addNoise_ (config().addNoise()), - noiseGenerator_ (config().noise_gen_conf(), engine_, 0), + noiseSampler_ (config().noise_gen_conf(), engine_), addRandomNoise_ (config().addRandomNoise()), diagLevel_ (config().diagLevel()) { @@ -110,16 +107,14 @@ namespace mu2e { private: - void makeDigitization (const CaloShowerROCollection&, CaloDigiCollection&, const CalSimParams&, const EventWindowMarker&, const ProtonBunchTimeMC&); - bool fillROHits (unsigned iRO, std::vector& waveform, const CaloShowerROCollection&, const ProtonBunchTimeMC&, const CalSimParams&); + void makeDigitization (const CaloShowerROCollection&, CaloDigiCollection&, const EventWindowMarker&, const ProtonBunchTimeMC&); + bool fillROHits (unsigned iRO, std::vector& waveform, const CaloShowerROCollection&, const ProtonBunchTimeMC&); void generateSpotNoise (std::vector& waveform); void buildOutputDigi (unsigned iRO, std::vector& waveform, double pedestal, CaloDigiCollection&); + void extract (const std::vector& wf, std::vector& starts, std::vector& stops) const; void diag0 (unsigned, const std::vector&); void diag1 (unsigned, double, size_t, const std::vector&, int); - void plotWF (const std::vector& waveform, const std::string& pname, int pedestal); - void plotWF (const std::vector& waveform, const std::string& pname, int pedestal); - ProditionsHandle calCrystalConds_; const art::ProductToken caloShowerToken_; art::InputTag ewMarkerTag_; art::InputTag pbtmcTag_; @@ -130,12 +125,15 @@ namespace mu2e { unsigned bufferDigi_; float startTimeBuffer_; int maxADCCounts_; - float minNoiseAmplitude_; + float pePerMeVCsI_; + float pePerMeVLyso_; + float MeVToADC_; CaloPulseShape pulseShape_; - CaloWFExtractor wfExtractor_; + unsigned nBinsPeak_; + int minPeakADC_; CLHEP::HepRandomEngine& engine_; bool addNoise_; - CaloNoiseSimGenerator noiseGenerator_; + CaloNoiseGenerator noiseSampler_; bool addRandomNoise_; const Calorimeter* calorimeter_; int diagLevel_; @@ -146,7 +144,8 @@ namespace mu2e { void CaloDigiMaker::beginRun(art::Run& aRun) { pulseShape_.buildShapes(); - if (addNoise_) noiseGenerator_.initialize(wfExtractor_); + //if (addNoise_) noiseSampler_.initialize(); + //if (addNoise_) noiseSampler_.dumpNoise("noise_0.root"); } @@ -166,8 +165,6 @@ namespace mu2e { event.getByLabel(pbtmcTag_, pbtmcHandle); const ProtonBunchTimeMC& pbtmc(*pbtmcHandle); - const auto& calCrystalConds = calCrystalConds_.get(event.id()); - ProditionsHandle eventTimingHandle; const EventTiming &eventTiming = eventTimingHandle.get(event.id()); timeFromProtonsToDRMarker_ = eventTiming.timeFromProtonsToDRMarker(); @@ -176,7 +173,7 @@ namespace mu2e { const auto& CaloShowerROs = *caloShowerStepHandle; auto caloDigiColl = std::make_unique(); - makeDigitization(CaloShowerROs, *caloDigiColl, calCrystalConds, ewMarker, pbtmc); + makeDigitization(CaloShowerROs, *caloDigiColl,ewMarker, pbtmc); event.put(std::move(caloDigiColl)); if ( diagLevel_ > 0 ) std::cout<<"[CaloDigiMaker::produce] end" << std::endl; @@ -186,8 +183,7 @@ namespace mu2e { //----------------------------------------------------------------------------------------------------------------------------- // Note: DigitizationStart include the fixed delay from timeFromProtonsToDRMarker, need to subtract it to be in the digitizer frame void CaloDigiMaker::makeDigitization(const CaloShowerROCollection& CaloShowerROs, CaloDigiCollection& caloDigiColl, - const CalSimParams& calCrystalConds, const EventWindowMarker& ewMarker, - const ProtonBunchTimeMC& pbtmc) + const EventWindowMarker& ewMarker, const ProtonBunchTimeMC& pbtmc) { mu2e::GeomHandle ch; calorimeter_ = ch.get(); @@ -207,20 +203,19 @@ namespace mu2e { for (int iRO=0;iRO& waveform, const CaloShowerROCollection& CaloShowerROs, - const ProtonBunchTimeMC& pbtmc, const CalSimParams& calCrystalConds) + const ProtonBunchTimeMC& pbtmc) { - bool isEmpty = true; - auto SiPMID = CaloSiPMId(iRO); - auto crystalID = SiPMID.crystal(); - auto pePerMeV = calCrystalConds.pePerMeVs(crystalID).at(SiPMID.SiPMLocalId()); - auto ADCPerMeV = calCrystalConds.ADCPerMeVs(crystalID).at(SiPMID.SiPMLocalId()); - auto scaleFactor = ADCPerMeV/pePerMeV; + bool isEmpty = true; + bool isCaphri = CaloSiPMId(iRO).crystal().isCaphri(); + auto pePerMeV = isCaphri ? pePerMeVLyso_ : pePerMeVCsI_; + auto scaleFactor = MeVToADC_/pePerMeV; for (const auto& CaloShowerRO : CaloShowerROs) { @@ -262,6 +255,8 @@ namespace mu2e { //---------------------------------------------------------------------------------------------------------- void CaloDigiMaker::generateSpotNoise(std::vector& waveform) { + double minAmplitude = 0.1*MeVToADC_; + size_t timeSample(0); std::vector hitStarts{}, hitStops{}; hitStarts.reserve(16);hitStops.reserve(16); @@ -269,11 +264,11 @@ namespace mu2e { // First, find the ranges in the waveform with non-zero bins. while (timeSample < waveform.size()) { - if (waveform[timeSample] < minNoiseAmplitude_) {++timeSample; continue;} + if (waveform[timeSample] < minAmplitude) {++timeSample; continue;} size_t sampleStart = (timeSample > bufferDigi_) ? timeSample - bufferDigi_ : 0; size_t sampleStop(timeSample); - while (sampleStop < waveform.size() && waveform[sampleStop] > minNoiseAmplitude_) ++sampleStop; + while (sampleStop < waveform.size() && waveform[sampleStop] > minAmplitude) ++sampleStop; hitStarts.push_back(sampleStart); hitStops.push_back(sampleStop); @@ -295,7 +290,10 @@ namespace mu2e { //Now take a random part of the noise waveform and add it to the waveform content for (size_t ihit=0; ihit hitStarts, hitStops; hitStarts.reserve(16);hitStops.reserve(16); - wfExtractor_.extract(wf,hitStarts,hitStops); + extract(wf,hitStarts,hitStops); // Build digi for concatenated hits for (size_t ihit=0;ihit& wf, std::vector& starts, std::vector& stops) const + { + size_t timeSample(nBinsPeak_+bufferDigi_); + while (timeSample+nBinsPeak_ < wf.size()) + { + // find starting point + if (wf[timeSample] < minPeakADC_) {++timeSample; continue;} + + size_t imax(timeSample-nBinsPeak_); + for (auto i = timeSample-nBinsPeak_; i<=timeSample+nBinsPeak_;++i) {if (wf[i]>wf[imax]) imax=i;} + if (timeSample != imax) {++timeSample; continue;} + + // find the starting / stopping point of the peak (stop = first value under threshold) + size_t sampleStart = (timeSample > bufferDigi_) ? timeSample - bufferDigi_ : 0; + size_t sampleStop(timeSample); + ++sampleStop; + while (sampleStop < wf.size() && wf[sampleStop] >= minPeakADC_) ++sampleStop; + + starts.push_back(sampleStart); + stops.push_back(sampleStop); + + //fast forward to end of waveform to search for next one + timeSample = sampleStop+1; + } + + + // Concatenate peaks and remove unused values (flag value to remove past wf.size() since the latter is a legitimate value) + size_t iprev(0), icurrent(1); + while (icurrent < starts.size()) + { + if (stops[iprev] >= starts[icurrent]) {stops[iprev]=stops[icurrent]; starts[icurrent]=stops[icurrent]=wf.size()+1;} + else {iprev = icurrent;} + ++icurrent; + } + + auto pred = [&wf](const auto a) {return a>wf.size();}; + starts.erase(std::remove_if(starts.begin(),starts.end(),pred),starts.end()); + stops.erase(std::remove_if(stops.begin(), stops.end(), pred),stops.end()); + } + @@ -359,36 +398,6 @@ namespace mu2e { std::cout<& waveform, const std::string& pname, int pedestal) - { - double startTime = digitizationStart_ - timeFromProtonsToDRMarker_ - startTimeBuffer_; - TH1F h("h","Waveform",waveform.size(),startTime,waveform.size()*digiSampling_+startTime); - for (size_t i=1;i<=waveform.size();++i) h.SetBinContent(i,waveform[i-1]); - TLine line; - line.SetLineStyle(2); - TLine line2; - line2.SetLineStyle(3); - - gStyle->SetOptStat(0); - TCanvas c1("c1","c1"); - h.Draw(); - line.DrawLine(startTime,pedestal,waveform.size()*digiSampling_+startTime,pedestal); - line2.DrawLine(startTime,pedestal+16,waveform.size()*digiSampling_+startTime,pedestal+16); - c1.SaveAs(pname.c_str()); - } - - void CaloDigiMaker::plotWF(const std::vector& waveform, const std::string& pname, int pedestal) - { - std::vector v; - for (size_t i=0;i -#include -#include -#include -#include - - - - -namespace mu2e { - - CaloNoiseSimGenerator::CaloNoiseSimGenerator(const Config& config, CLHEP::HepRandomEngine& engine, int iRO) : - iRO_ (iRO), - waveform_ (config.noiseWFSize(),0.0), - pedestal_ (0.0), - digiNoise_ (), - digiNoiseProb_ (), - digiSampling_ (config.digiSampling()), - noiseRinDark_ (config.rinNphotPerNs() + config.darkNphotPerNs()), - noiseElec_ (config.elecNphotPerNs()), - minPeakADC_ (config.minPeakADC()), - pePerMeV_ (config.pePerMeV()), - MeVToADC_ (config.MeVToADC()), - randPoisson_ (engine), - randGauss_ (engine), - randFlat_ (engine), - nMaxFragment_ (config.nMaxFragment()), - pulseShape_ (config.pulseFileName(),config.pulseHistName(),digiSampling_), - diagLevel_ (config.diagLevel()) - {} - - - //------------------------------------------------------------------------------------------------------------------ - void CaloNoiseSimGenerator::initialize(const CaloWFExtractor& wfExtractor) - { - pulseShape_.buildShapes(); - generateWF(waveform_); - generateFragments(wfExtractor); - } - - //------------------------------------------------------------------------------------------------------------------ - void CaloNoiseSimGenerator::generateWF(std::vector& wfVector) - { - float scaleFactor(MeVToADC_/pePerMeV_); - - std::fill(wfVector.begin(),wfVector.end(),0); - - const auto& pulse = pulseShape_.digitizedPulse(0.0); - const unsigned pulseSize = pulse.size(); - const unsigned bufferSize = int(0.75*pulseSize); - const unsigned noiseSize = wfVector.size(); - const double totalTime = (noiseSize+bufferSize)*digiSampling_; - const int noiseLevelPE = int(totalTime*noiseRinDark_); - - //Generate the radiation induced noise (RIN) - const int nPh = randPoisson_(noiseLevelPE); - for (int i=0;i temp(length,0.0); - - std::vector wf; - wf.reserve(temp.size()); - for (const auto& val : temp) wf.emplace_back(val - pedestal_); - - std::vector starts, stops; - starts.reserve(16); stops.reserve(16); - wfExtractor.extract(wf,starts,stops); - if (starts.empty()) continue; - - std::vector fragment; - fragment.reserve(stops[0]-starts[0]); - std::copy(temp.begin()+starts[0], temp.begin()+stops[0]+1, std::back_inserter(fragment)); - digiNoise_.push_back(fragment); - - ++nfound; - if (nfound==enoughFragments) break; - } - - digiNoiseProb_= float(nfound)/float(nwf)/float(length); - } - - //------------------------------------------------------------------------------------------------------------------ - void CaloNoiseSimGenerator::refresh() {generateWF(waveform_);} - - - - //------------------------------------------------------------------------------------------------------------------ - void CaloNoiseSimGenerator::addSampleNoise(std::vector& wfVector, unsigned istart, unsigned ilength) - { - if (ilength >=waveform_.size()) - throw cet::exception("CATEGORY")<<"[CaloNoiseSimGenerator] noise length request too long"; - - unsigned irandom = unsigned(randFlat_.fire(0.,waveform_.size()-ilength)); - for (unsigned i=0;i& wfVector) - { - double muNoise = waveform_.size()*digiNoiseProb_; - int nNoise = randPoisson_(muNoise); - for (int in=0;in& digi = digiNoise_[idigi]; - if (wfVector.size() < digi.size()) continue; - - unsigned istart = unsigned(randFlat_.fire(0.,wfVector.size()-digi.size())); - for (unsigned i=0;i x(waveform_.size()),y; - std::iota(x.begin(),x.end(),0); - TGraph gr(x.size(),x.data(),waveform_.data()); - gr.SetTitle("Original waveform"); - - TH1F h1("h1","Projection waveform",100,-50,50); - for (const auto& val: waveform_) h1.Fill(val-pedestal_); - - TCanvas c1("c1","c1"); - c1.Divide(2,2); - c1.cd(1); - gr.Draw("AL"); - c1.cd(2); - h1.Draw(); - c1.SaveAs(name.c_str()); - } - -} - - diff --git a/CaloMC/src/CaloPhotonPropagation.cc b/CaloMC/src/CaloPhotonPropagation.cc index d0165a128f..e92cd2d34e 100644 --- a/CaloMC/src/CaloPhotonPropagation.cc +++ b/CaloMC/src/CaloPhotonPropagation.cc @@ -1,86 +1,93 @@ -#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" #include "Offline/CaloMC/inc/CaloPhotonPropagation.hh" +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" #include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" -#include "Offline/SeedService/inc/SeedService.hh" +#include "cetlib_except/exception.h" -#include "CLHEP/Random/RandFlat.h" #include "TFile.h" #include "TH2F.h" -#include +#include +#include +#include namespace mu2e { - CaloPhotonPropagation::CaloPhotonPropagation(const std::string& fileName, const std::string& histName, CLHEP::HepRandomEngine& engine) : - timeProp_ (), - cdf_ (), - nTimeDiv_ (0), - nZDiv_ (0), - dzTime_ (0), - randFlat_ (engine), - fileName_(fileName), - histName_(histName), - lightSpeed_(300) + CaloPhotonPropagation::CaloPhotonPropagation(const std::string& fileName, const std::string& histName, + CLHEP::HepRandomEngine& engine) : + randFlat_(engine), fileName_(fileName), histName_(histName) {} + //---------------------------------------------------------------------------------------------------------------------- void CaloPhotonPropagation::buildTable() { + constexpr float kCdfFloor = 1e-6f; ConfigFileLookupPolicy resolveFullPath; - std::string fullFileName = resolveFullPath(fileName_); + const std::string fullFileName = resolveFullPath(fileName_); - TH2F *hist(0); TFile file(fullFileName.c_str()); - if (file.IsOpen()) hist = (TH2F*) file.Get(histName_.c_str()); - if (!hist) throw cet::exception("CATEGORY")<<"CaloROStepMaker:: Histogram "<SetDirectory(0); + if (!file.IsOpen()) + throw cet::exception("CaloPhotonPropagation") + << "cannot open propagation file " << fullFileName << "\n"; + + // take ownership of the histogram so it survives the file and is freed on return + std::unique_ptr hist(dynamic_cast(file.Get(histName_.c_str()))); + if (!hist) + throw cet::exception("CaloPhotonPropagation") << "histogram " << histName_ + << " not found in " << fullFileName << "\n"; + hist->SetDirectory(nullptr); file.Close(); - dzTime_ = hist->GetXaxis()->GetBinWidth(1); + dz_ = hist->GetXaxis()->GetBinWidth(1); + nZDiv_ = hist->GetNbinsX(); nTimeDiv_ = hist->GetNbinsY(); - nZDiv_ = hist->GetNbinsX(); - for (unsigned iy=1;iy<=nTimeDiv_;++iy) timeProp_.push_back(hist->GetYaxis()->GetBinCenter(iy)); - cdf_.reserve(hist->GetNbinsX()*hist->GetNbinsY()); - for (int ix=1;ix<=hist->GetNbinsX();++ix) + timeProp_.reserve(nTimeDiv_); + for (unsigned iy=1; iy<=nTimeDiv_; ++iy) timeProp_.push_back(hist->GetYaxis()->GetBinCenter(iy)); + + // build one normalized cumulative distribution per depth slice + cdf_.reserve(nZDiv_*nTimeDiv_); + for (unsigned ix=1; ix<=nZDiv_; ++ix) { - float sum(1e-6); - std::vector temp; - for (int iy=1;iy<=hist->GetNbinsY();++iy) - { - sum += hist->GetBinContent(ix,iy); - temp.push_back(sum); - } - for (auto& val: temp) val /= sum; - std::copy(temp.begin(),temp.end(),std::back_inserter(cdf_)); + std::vector column; + column.reserve(nTimeDiv_); + + float sum = kCdfFloor; + for (unsigned iy=1; iy<=nTimeDiv_; ++iy) { sum += hist->GetBinContent(ix,iy); column.push_back(sum); } + for (float& v : column) v /= sum; + + cdf_.insert(cdf_.end(), column.begin(), column.end()); } const Calorimeter& cal = *(GeomHandle()); - lightSpeed_ = 300.0 / cal.G4Info().get("refractiveIndex"); //in mm/ns + lightSpeed_ = 300.0f / cal.G4Info().get("refractiveIndex"); // mm/ns } + //---------------------------------------------------------------------------- float CaloPhotonPropagation::propTimeSimu(float z) { - unsigned iz = z/dzTime_; - if(iz>=nZDiv_) iz = nZDiv_ - 1; - float test = randFlat_.fire(0.0,1.0); - unsigned ibin = nTimeDiv_*iz; - unsigned iend = ibin + nTimeDiv_ - 1; - - while (cdf_[ibin](z/dz_); + const unsigned iz = (izSigned <= 0) ? 0u : std::min(izSigned, nZDiv_-1); + + // inverse-CDF sampling: first time bin whose cumulative probability reaches the random draw + const float test = randFlat_.fire(0.0f,1.0f); + const auto first = cdf_.begin() + iz*nTimeDiv_; + const auto last = first + nTimeDiv_; + const auto hit = std::lower_bound(first, last, test); + + const std::size_t timeBin = std::min(std::distance(first,hit), nTimeDiv_-1); + return timeProp_[timeBin]; } + //---------------------------------------------------------------------------- - float CaloPhotonPropagation::propTimeLine(float z) + float CaloPhotonPropagation::propTimeLine(float z) const { return z/lightSpeed_; } - } - diff --git a/CaloMC/src/CaloWFExtractor.cc b/CaloMC/src/CaloWFExtractor.cc deleted file mode 100644 index 9b21efd3d9..0000000000 --- a/CaloMC/src/CaloWFExtractor.cc +++ /dev/null @@ -1,47 +0,0 @@ -#include "Offline/CaloMC/inc/CaloWFExtractor.hh" -#include -#include - -namespace mu2e { - - void CaloWFExtractor::extract(const std::vector& wf, std::vector& starts, std::vector& stops) const - { - size_t timeSample(nBinsPeak_+startOffset_); - while (timeSample+nBinsPeak_ < wf.size()) - { - // find starting point - if (wf[timeSample] < minPeakADC_) {++timeSample; continue;} - - size_t imax(timeSample-nBinsPeak_); - for (auto i = timeSample-nBinsPeak_; i<=timeSample+nBinsPeak_;++i) {if (wf[i]>wf[imax]) imax=i;} - if (timeSample != imax) {++timeSample; continue;} - - // find the starting / stopping point of the peak (stop = first value under threshold) - size_t sampleStart = (timeSample > bufferDigi_) ? timeSample - bufferDigi_ : 0; - size_t sampleStop(timeSample); - ++sampleStop; - while (sampleStop < wf.size() && wf[sampleStop] >= minPeakADC_) ++sampleStop; - - starts.push_back(sampleStart); - stops.push_back(sampleStop); - - //fast forward to end of waveform to search for next one - timeSample = sampleStop+1; - } - - - // Concatenate peaks and remove unused values (flag value to remove past wf.size() since the latter is a legitimate value) - size_t iprev(0), icurrent(1); - while (icurrent < starts.size()) - { - if (stops[iprev] >= starts[icurrent]) {stops[iprev]=stops[icurrent]; starts[icurrent]=stops[icurrent]=wf.size()+1;} - else {iprev = icurrent;} - ++icurrent; - } - - auto pred = [&wf](const auto a) {return a>wf.size();}; - starts.erase(std::remove_if(starts.begin(),starts.end(),pred),starts.end()); - stops.erase(std::remove_if(stops.begin(), stops.end(), pred),stops.end()); - } -} - From f9560d4a58ffa0a5679c2041244d3bed7a8ac77b Mon Sep 17 00:00:00 2001 From: echenard Date: Thu, 6 Aug 2026 16:00:33 -0500 Subject: [PATCH 3/7] Further Noise generator fixes --- CaloMC/CMakeLists.txt | 2 - CaloMC/fcl/prolog.fcl | 41 +----- CaloMC/src/CaloDigiMaker_module.cc | 159 ++++++++++---------- CaloMC/src/CaloHitTruthMatch_module.cc | 32 ++-- CaloReco/fcl/common.fcl | 24 ++- CaloReco/fcl/prolog.fcl | 3 +- CaloReco/inc/CaloTemplateWFProcessor.hh | 10 +- CaloReco/inc/CaloTemplateWFUtil.hh | 8 +- CaloReco/src/CaloTemplateWFProcessor.cc | 4 +- CaloReco/src/CaloTemplateWFUtil.cc | 12 +- Mu2eUtilities/CMakeLists.txt | 3 +- Mu2eUtilities/inc/CaloNoiseUtil.hh | 77 ++++++++++ Mu2eUtilities/inc/CaloPulseUtil.hh | 62 ++++++++ Mu2eUtilities/src/CaloNoiseUtil.cc | 188 ++++++++++++++++++++++++ Mu2eUtilities/src/CaloPulseUtil.cc | 113 ++++++++++++++ 15 files changed, 581 insertions(+), 157 deletions(-) create mode 100644 Mu2eUtilities/inc/CaloNoiseUtil.hh create mode 100644 Mu2eUtilities/inc/CaloPulseUtil.hh create mode 100644 Mu2eUtilities/src/CaloNoiseUtil.cc create mode 100644 Mu2eUtilities/src/CaloPulseUtil.cc diff --git a/CaloMC/CMakeLists.txt b/CaloMC/CMakeLists.txt index 9f9468b0a3..fdac25c801 100644 --- a/CaloMC/CMakeLists.txt +++ b/CaloMC/CMakeLists.txt @@ -1,8 +1,6 @@ cet_make_library( SOURCE - src/CaloNoiseSimGenerator.cc src/CaloPhotonPropagation.cc - src/CaloWFExtractor.cc src/ShowerStepUtil.cc src/CaloDigiWrapper.cc src/CaloDigiWrapperCollection.cc diff --git a/CaloMC/fcl/prolog.fcl b/CaloMC/fcl/prolog.fcl index 8db3b1935d..698424b61c 100644 --- a/CaloMC/fcl/prolog.fcl +++ b/CaloMC/fcl/prolog.fcl @@ -7,45 +7,20 @@ BEGIN_PROLOG -CaloMC : { - - #See doc-db 35519 for noise level discussion - NoiseGenerator : - { - noiseWFSize : 10000 - pulseFileName : @local::pulseFileName - pulseHistName : @local::pulseHistName - elecNphotPerNs : 0.18 - rinNphotPerNs : 0.25 - darkNphotPerNs : 0.00 -# Run2 values -# rinNphotPerNs : 0.36 -# darkNphotPerNs : 0.60 - digiSampling : @local::HitMakerDigiSampling - nMaxFragment : 1000 - minPeakADC : @local::HitMakerMinPeakADC - readoutPEPerMeV : 30 - MeVToADC : 16 - diagLevel : 0 - } -} +CaloMC : { -CaloMC : { @table::CaloMC CaloShowerStepMaker : { module_type : CaloShowerStepMaker numZSlices : 20 deltaTime : 0.2 caloStepPointCollection : ["g4run:calorimeter"] - physVolInfoInput : "g4run" - caloMaterial : ["G4_CESIUM_IODIDE", "Polyethylene092","CarbonFiber"] compressData : true eDepThreshold : 0 diagLevel : 0 } - CaloShowerROMaker : { module_type : CaloShowerROMaker @@ -70,19 +45,16 @@ CaloMC : { @table::CaloMC caloShowerROCollection : CaloShowerROMaker eventWindowMarker : EWMProducer protonBunchTimeMC : EWMProducer - pulseFileName : @local::pulseFileName - pulseHistName : @local::pulseHistName + pulseCache : { @table::CaloPulseCache } + noiseCache : { @table::CaloNoiseCache } digitizationStart : @local::HitMakerDigitizationStart digitizationEnd : @local::HitMakerDigitizationEnd - addNoise : true - addRandomNoise : false - NoiseGenerator : { @table::CaloMC.NoiseGenerator } digiSampling : @local::HitMakerDigiSampling + bufferDigi : 16 + addNoise : true nBits : 12 minPeakADC : @local::HitMakerMinPeakADC nBinsPeak : 2 - bufferDigi : 16 - minNoiseAmplitude : 2 diagLevel : 0 } @@ -92,8 +64,7 @@ CaloMC : { @table::CaloMC caloShowerSimCollection : CaloShowerROMaker caloHitCollection : CaloHitMaker primaryParticle : FindMCPrimary - pulseFileName : @local::pulseFileName - pulseHistName : @local::pulseHistName + pulseCache : { @table::CaloPulseCache } digiSampling : @local::HitMakerDigiSampling minAmplitude : 2 fillDetailedMC : false diff --git a/CaloMC/src/CaloDigiMaker_module.cc b/CaloMC/src/CaloDigiMaker_module.cc index d55700bfb7..01376eecab 100644 --- a/CaloMC/src/CaloDigiMaker_module.cc +++ b/CaloMC/src/CaloDigiMaker_module.cc @@ -13,9 +13,10 @@ #include "art_root_io/TFileService.h" #include "art_root_io/TFileDirectory.h" -#include "Offline/Mu2eUtilities/inc/CaloPulseShape.hh" +#include "Offline/CaloConditions/inc/CalSimParams.hh" +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" +#include "Offline/Mu2eUtilities/inc/CaloNoiseUtil.hh" #include "Offline/CalorimeterGeom/inc/Calorimeter.hh" -#include "Offline/CaloMC/inc/CaloNoiseGenerator.hh" #include "Offline/ProditionsService/inc/ProditionsHandle.hh" #include "Offline/DAQConditions/inc/EventTiming.hh" #include "Offline/DataProducts/inc/EventWindowMarker.hh" @@ -47,21 +48,17 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - using CNG = mu2e::CaloNoiseGenerator::Config; - fhicl::Table noise_gen_conf { Name("NoiseGenerator"), Comment("Noise generator config") }; + using CNG = CaloNoiseUtil::Config; + using CPG = CaloPulseUtil::Config; + fhicl::Table noiseCache { Name("noiseCache"), Comment("Noise cache maker config") }; + fhicl::Table pulseCache { Name("pulseCache"), Comment("Pulse cache maker config") }; fhicl::Atom caloShowerCollection { Name("caloShowerROCollection"), Comment("CaloShowerRO collection name") }; fhicl::Atom ewMarkerTag { Name("eventWindowMarker"), Comment("EventWindowMarker producer") }; fhicl::Atom pbtmcTag { Name("protonBunchTimeMC"), Comment("ProtonBunchTimeMC producer") }; - fhicl::Atom pulseFileName { Name("pulseFileName"), Comment("Calo pulse file name") }; - fhicl::Atom pulseHistName { Name("pulseHistName"), Comment("Calo pulse hist name") }; fhicl::Atom digitizationStart { Name("digitizationStart"), Comment("Start of digitization window relative to nominal pb time") }; fhicl::Atom digitizationEnd { Name("digitizationEnd"), Comment("End of digitization window relative to nominal pb time")}; fhicl::Atom addNoise { Name("addNoise"), Comment("Add noise to waveform") }; - fhicl::Atom addRandomNoise { Name("addRandomNoise"), Comment("Add random salt and pepper noise") }; fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; - fhicl::Atom pePerMeVCsI { Name("readoutPEPerMeVCsI"), Comment("Number of pe / MeV for Readout for CsI") }; - fhicl::Atom pePerMeVLyso { Name("readoutPEPerMeVLyso"), Comment("Number of pe / MeV for Readout for LYSO") }; - fhicl::Atom MeVToADC { Name("MeVToADC"), Comment("MeV to ADC conversion factor") }; fhicl::Atom nBits { Name("nBits"), Comment("ADC Number of bits") }; fhicl::Atom nBinsPeak { Name("nBinsPeak"), Comment("Window size for finding local maximum to digitize wf") }; fhicl::Atom minPeakADC { Name("minPeakADC"), Comment("Minimum ADC hits of local peak to digitize") }; @@ -80,16 +77,12 @@ namespace mu2e { bufferDigi_ (config().bufferDigi()), startTimeBuffer_ (config().digiSampling()*config().bufferDigi()), maxADCCounts_ ((1 << config().nBits()) -1), - pePerMeVCsI_ (config().pePerMeVCsI()), - pePerMeVLyso_ (config().pePerMeVLyso()), - MeVToADC_ (config().MeVToADC()), - pulseShape_ (CaloPulseShape(config().pulseFileName(),config().pulseHistName(),config().digiSampling())), + pulseCache_ (CaloPulseUtil(config().pulseCache())), nBinsPeak_ (config().nBinsPeak()), minPeakADC_ (config().minPeakADC()), engine_ (createEngine(art::ServiceHandle()->getSeed())), addNoise_ (config().addNoise()), - noiseSampler_ (config().noise_gen_conf(), engine_), - addRandomNoise_ (config().addRandomNoise()), + noiseSampler_ (config().noiseCache(), engine_), diagLevel_ (config().diagLevel()) { consumes(ewMarkerTag_); @@ -98,7 +91,8 @@ namespace mu2e { produces(); //check that StartTimeBuffer is shorter than BlindTime_ - if (startTimeBuffer_ > digitizationStart_) throw cet::exception("CATEGORY")<< "CaloDigiMaker: digitizationStart is too small to accommodate start time buffer"; + if (startTimeBuffer_ > digitizationStart_) throw cet::exception("CALODIGIMAKER") + << "digitizationStart too small to accommodate start time buffer"; } void produce(art::Event& e) override; @@ -107,14 +101,16 @@ namespace mu2e { private: - void makeDigitization (const CaloShowerROCollection&, CaloDigiCollection&, const EventWindowMarker&, const ProtonBunchTimeMC&); - bool fillROHits (unsigned iRO, std::vector& waveform, const CaloShowerROCollection&, const ProtonBunchTimeMC&); - void generateSpotNoise (std::vector& waveform); - void buildOutputDigi (unsigned iRO, std::vector& waveform, double pedestal, CaloDigiCollection&); - void extract (const std::vector& wf, std::vector& starts, std::vector& stops) const; - void diag0 (unsigned, const std::vector&); - void diag1 (unsigned, double, size_t, const std::vector&, int); + void makeDigitization (const CaloShowerROCollection&, CaloDigiCollection&, const EventWindowMarker&, const ProtonBunchTimeMC&, const CalSimParams&); + bool fillROHits (unsigned iRO, std::vector& waveform, const CaloShowerROCollection&, const ProtonBunchTimeMC&, const CalSimParams&); + void generateSpotNoise (std::vector& waveform, double scaleFactor); + void buildOutputDigi (unsigned iRO, std::vector& waveform, double pedestal, CaloDigiCollection&); + void extract (const std::vector& wf, std::vector& starts, std::vector& stops) const; + double readoutScaleFactor(unsigned iRO, const CalSimParams& conds) const; + void diag0 (unsigned, const std::vector&); + void diag1 (unsigned, double, size_t, const std::vector&, int); + ProditionsHandle calCrystalConds_; const art::ProductToken caloShowerToken_; art::InputTag ewMarkerTag_; art::InputTag pbtmcTag_; @@ -125,16 +121,12 @@ namespace mu2e { unsigned bufferDigi_; float startTimeBuffer_; int maxADCCounts_; - float pePerMeVCsI_; - float pePerMeVLyso_; - float MeVToADC_; - CaloPulseShape pulseShape_; + CaloPulseUtil pulseCache_; unsigned nBinsPeak_; int minPeakADC_; CLHEP::HepRandomEngine& engine_; bool addNoise_; - CaloNoiseGenerator noiseSampler_; - bool addRandomNoise_; + CaloNoiseUtil noiseSampler_; const Calorimeter* calorimeter_; int diagLevel_; }; @@ -143,9 +135,7 @@ namespace mu2e { //----------------------------------------------------------------------------- void CaloDigiMaker::beginRun(art::Run& aRun) { - pulseShape_.buildShapes(); - //if (addNoise_) noiseSampler_.initialize(); - //if (addNoise_) noiseSampler_.dumpNoise("noise_0.root"); + pulseCache_.buildCache(); } @@ -172,8 +162,10 @@ namespace mu2e { auto caloShowerStepHandle = event.getValidHandle(caloShowerToken_); const auto& CaloShowerROs = *caloShowerStepHandle; + const auto& calCrystalConds = calCrystalConds_.get(event.id()); + auto caloDigiColl = std::make_unique(); - makeDigitization(CaloShowerROs, *caloDigiColl,ewMarker, pbtmc); + makeDigitization(CaloShowerROs, *caloDigiColl,ewMarker, pbtmc,calCrystalConds); event.put(std::move(caloDigiColl)); if ( diagLevel_ > 0 ) std::cout<<"[CaloDigiMaker::produce] end" << std::endl; @@ -183,7 +175,7 @@ namespace mu2e { //----------------------------------------------------------------------------------------------------------------------------- // Note: DigitizationStart include the fixed delay from timeFromProtonsToDRMarker, need to subtract it to be in the digitizer frame void CaloDigiMaker::makeDigitization(const CaloShowerROCollection& CaloShowerROs, CaloDigiCollection& caloDigiColl, - const EventWindowMarker& ewMarker, const ProtonBunchTimeMC& pbtmc) + const EventWindowMarker& ewMarker, const ProtonBunchTimeMC& pbtmc, const CalSimParams& calCrystalConds) { mu2e::GeomHandle ch; calorimeter_ = ch.get(); @@ -203,59 +195,25 @@ namespace mu2e { for (int iRO=0;iRO& waveform, const CaloShowerROCollection& CaloShowerROs, - const ProtonBunchTimeMC& pbtmc) - { - bool isEmpty = true; - bool isCaphri = CaloSiPMId(iRO).crystal().isCaphri(); - auto pePerMeV = isCaphri ? pePerMeVLyso_ : pePerMeVCsI_; - auto scaleFactor = MeVToADC_/pePerMeV; - - for (const auto& CaloShowerRO : CaloShowerROs) - { - unsigned SiPMID = CaloShowerRO.SiPMID(); - if (SiPMID != iRO) continue; - - isEmpty = false; - for (const auto PEtime : CaloShowerRO.PETime()) - { - //PE time is given in DR frame, we need to subtract the event window start and the digi Start time - float time = PEtime + pbtmc.pbtime_- digitizationStart_ + timeFromProtonsToDRMarker_ + startTimeBuffer_; - unsigned startSample = std::max(0u,unsigned(time/digiSampling_)); - const auto& pulse = pulseShape_.digitizedPulse(time); - unsigned stopSample = std::min(startSample+pulse.size(), waveform.size()); - - for (size_t timeSample = startSample; timeSample < stopSample; ++timeSample) - waveform.at(timeSample) += pulse.at(timeSample - startSample)*scaleFactor; - } - } - return isEmpty; - } - - //---------------------------------------------------------------------------------------------------------- - void CaloDigiMaker::generateSpotNoise(std::vector& waveform) + void CaloDigiMaker::generateSpotNoise(std::vector& waveform, double scaleFactor) { - double minAmplitude = 0.1*MeVToADC_; + const int NoiseWFID(0); // will need this from proditions later; + + const double minAmplitude(0.01); + noiseSampler_.prepare(NoiseWFID, scaleFactor); size_t timeSample(0); std::vector hitStarts{}, hitStops{}; @@ -292,11 +250,39 @@ namespace mu2e { { unsigned istart = hitStarts[ihit]; unsigned ilength = hitStops[ihit]-hitStarts[ihit]; - const auto& noiseWF = noiseSampler_.noiseSegment(0,istart,ilength); + const auto& noiseWF = noiseSampler_.noiseSegment(NoiseWFID,istart,ilength); for (unsigned i=0;i& waveform, const CaloShowerROCollection& CaloShowerROs, + const ProtonBunchTimeMC& pbtmc, const CalSimParams& calCrystalConds) + { + bool isEmpty = true; + const double scaleFactor = readoutScaleFactor(iRO, calCrystalConds); + + for (const auto& CaloShowerRO : CaloShowerROs) + { + unsigned SiPMID = CaloShowerRO.SiPMID(); + if (SiPMID != iRO) continue; + + isEmpty = false; + for (const auto PEtime : CaloShowerRO.PETime()) + { + //PE time is given in DR frame, we need to subtract the event window start and the digi Start time + float time = PEtime + pbtmc.pbtime_- digitizationStart_ + timeFromProtonsToDRMarker_ + startTimeBuffer_; + unsigned startSample = std::max(0u,unsigned(time/digiSampling_)); + const auto& pulse = pulseCache_.digitizedPulse(time); + unsigned stopSample = std::min(startSample+pulse.size(), waveform.size()); + + for (size_t timeSample = startSample; timeSample < stopSample; ++timeSample) + waveform.at(timeSample) += pulse.at(timeSample - startSample)*scaleFactor; + } + } + return isEmpty; + } + //------------------------------------------------------------------------------------------------------------------- void CaloDigiMaker::buildOutputDigi(unsigned iRO, std::vector& waveform, double pedestal, CaloDigiCollection& caloDigiColl) @@ -380,6 +366,15 @@ namespace mu2e { stops.erase(std::remove_if(stops.begin(), stops.end(), pred),stops.end()); } + //------------------------------------------------------------------------------------------------------------------- + double CaloDigiMaker::readoutScaleFactor(unsigned iRO, const CalSimParams& conds) const + { + const auto SiPMID = CaloSiPMId(iRO); + const auto crystalID = SiPMID.crystal(); + const auto pePerMeV = conds.pePerMeVs (crystalID).at(SiPMID.SiPMLocalId()); + const auto ADCPerMeV = conds.ADCPerMeVs(crystalID).at(SiPMID.SiPMLocalId()); + return ADCPerMeV / pePerMeV; + } diff --git a/CaloMC/src/CaloHitTruthMatch_module.cc b/CaloMC/src/CaloHitTruthMatch_module.cc index 79867c1004..9a9745598e 100644 --- a/CaloMC/src/CaloHitTruthMatch_module.cc +++ b/CaloMC/src/CaloHitTruthMatch_module.cc @@ -13,7 +13,7 @@ #include "Offline/MCDataProducts/inc/SimParticle.hh" #include "Offline/MCDataProducts/inc/MCRelationship.hh" #include "Offline/RecoDataProducts/inc/CaloHit.hh" -#include "Offline/Mu2eUtilities/inc/CaloPulseShape.hh" +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" #include #include @@ -34,16 +34,16 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - fhicl::Atom caloShowerSimCollection { Name("caloShowerSimCollection"), Comment("Name of caloShowerSim Collection") }; - fhicl::Atom caloHitCollection { Name("caloHitCollection"), Comment("Name of CaloHit collection") }; - fhicl::Atom primaryParticle { Name("primaryParticle"), Comment("PrimaryParticle producer")}; - fhicl::Atom pulseFileName { Name("pulseFileName"), Comment("Calo pulse file name") }; - fhicl::Atom pulseHistName { Name("pulseHistName"), Comment("Calo pulse hist name") }; - fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; - fhicl::Atom deltaTimeMinus { Name("deltaTimeMinus"), Comment("Max time (ns) a MC hit may precede the reco hit to be matched"), 100.0 }; - fhicl::Atom minAmplitude { Name("minAmplitude"), Comment("Minimum amplitude of waveform to define hit length") }; - fhicl::Atom fillDetailedMC { Name("fillDetailedMC"), Comment("Fill SimParticle - SimShower Assn map")}; - fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; + using CPG = CaloPulseUtil::Config; + fhicl::Table pulseCache { Name("pulseCache"), Comment("Pulse cache maker config") }; + fhicl::Atom caloShowerSimCollection { Name("caloShowerSimCollection"), Comment("Name of caloShowerSim Collection") }; + fhicl::Atom caloHitCollection { Name("caloHitCollection"), Comment("Name of CaloHit collection") }; + fhicl::Atom primaryParticle { Name("primaryParticle"), Comment("PrimaryParticle producer")}; + fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitizer sampling time (ns) ") }; + fhicl::Atom deltaTimeMinus { Name("deltaTimeMinus"), Comment("Max time (ns) a MC hit may precede the reco hit to be matched"), 100.0 }; + fhicl::Atom minAmplitude { Name("minAmplitude"), Comment("Minimum amplitude of waveform to define hit length") }; + fhicl::Atom fillDetailedMC { Name("fillDetailedMC"), Comment("Fill SimParticle - SimShower Assn map")}; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; }; @@ -52,8 +52,7 @@ namespace mu2e { caloShowerSimToken_ {consumes(config().caloShowerSimCollection())}, caloHitToken_ {consumes(config().caloHitCollection())}, ppToken_ {consumes(config().primaryParticle())}, - pulseFileName_ (config().pulseFileName()), - pulseHistName_ (config().pulseHistName()), + pulseCache_ (config().pulseCache()), digiSampling_ (config().digiSampling()), deltaTimeMinus_ (config().deltaTimeMinus()), minAmplitude_ (config().minAmplitude()), @@ -78,8 +77,7 @@ namespace mu2e { const art::ProductToken caloShowerSimToken_; const art::ProductToken caloHitToken_; const art::ProductToken ppToken_; - std::string pulseFileName_; - std::string pulseHistName_; + CaloPulseUtil::Config pulseCache_; double digiSampling_; double deltaTimeMinus_; double minAmplitude_; @@ -95,8 +93,8 @@ namespace mu2e { //----------------------------------------------------------------------------- void CaloHitTruthMatch::beginRun(art::Run&) { - CaloPulseShape cps(pulseFileName_,pulseHistName_,digiSampling_); - cps.buildShapes(); + CaloPulseUtil cps(pulseCache_); + cps.buildCache(); wf_ = cps.digitizedPulse(0); wfBinMax_ = std::distance(wf_.begin(),std::max_element(wf_.begin(),wf_.end())); diff --git a/CaloReco/fcl/common.fcl b/CaloReco/fcl/common.fcl index 5efd0185b1..cff9ee882f 100644 --- a/CaloReco/fcl/common.fcl +++ b/CaloReco/fcl/common.fcl @@ -5,9 +5,29 @@ BEGIN_PROLOG HitMakerMinPeakADC : 16 pulseFileName : "OfflineData/ConditionsService/CsI-waveform-2020-02-12.root" pulseHistName : "h_waveform" - readoutPEPerMeVCsI : 30.0 - readoutPEPerMeVLyso : 500.0 MeVToADC : 16.0 ADCToMeV : 0.0625 + + CaloPulseCache:{ + fileName : "OfflineData/ConditionsService/CsI-waveform-2020-02-12.root" + histName : "h_waveform" + digiSampling : @local::HitMakerDigiSampling + } + + #See doc-db 35519 for noise level discussion + CaloNoiseCache : + { + generate : false + dumpGenerated : false + histoFileName : "CaloMC/data/ReadoutNoise.root" + histoPrefix : "histo_" + pulseCache : @local::CaloPulseCache + elecNphotPerNs : 0.18 + rinNphotPerNs : 0.25 + darkNphotPerNs : 0.00 + digiSampling : @local::HitMakerDigiSampling + } + + END_PROLOG diff --git a/CaloReco/fcl/prolog.fcl b/CaloReco/fcl/prolog.fcl index d075acb7db..4c99555010 100644 --- a/CaloReco/fcl/prolog.fcl +++ b/CaloReco/fcl/prolog.fcl @@ -18,8 +18,7 @@ CaloReco : { TemplateProcessor : { - pulseFileName : @local::pulseFileName - pulseHistName : @local::pulseHistName + pulseCache : {@table::CaloPulseCache} windowPeak : 3 minPeakAmplitude : 24 numNoiseBins : 5 diff --git a/CaloReco/inc/CaloTemplateWFProcessor.hh b/CaloReco/inc/CaloTemplateWFProcessor.hh index 0f9abd1541..d44c222724 100644 --- a/CaloReco/inc/CaloTemplateWFProcessor.hh +++ b/CaloReco/inc/CaloTemplateWFProcessor.hh @@ -15,10 +15,12 @@ // --> the first peak is never a pile-up // -#include "Offline/CaloReco/inc/CaloWaveformProcessor.hh" -#include "Offline/CaloReco/inc/CaloTemplateWFUtil.hh" #include "fhiclcpp/types/Atom.h" #include "fhiclcpp/types/Sequence.h" +#include "fhiclcpp/types/Table.h" +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" +#include "Offline/CaloReco/inc/CaloWaveformProcessor.hh" +#include "Offline/CaloReco/inc/CaloTemplateWFUtil.hh" #include "TH2.h" #include @@ -33,8 +35,8 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - fhicl::Atom pulseFileName { Name("pulseFileName"), Comment("Calo pulse file name") }; - fhicl::Atom pulseHistName { Name("pulseHistName"), Comment("Calo pulse hist name") }; + using CPG = CaloPulseUtil::Config; + fhicl::Table pulseCache { Name("pulseCache"), Comment("Pulse cache maker config") }; fhicl::Atom windowPeak { Name("windowPeak"), Comment("Number of bins around central value to inspect") }; fhicl::Atom minPeakAmplitude { Name("minPeakAmplitude"), Comment("Minimum peak amplitude") }; fhicl::Atom minDTPeaks { Name("minDTPeaks"), Comment("Minimum time difference between consecutive peaks") }; diff --git a/CaloReco/inc/CaloTemplateWFUtil.hh b/CaloReco/inc/CaloTemplateWFUtil.hh index 01687c829f..b6b2dcfaa0 100644 --- a/CaloReco/inc/CaloTemplateWFUtil.hh +++ b/CaloReco/inc/CaloTemplateWFUtil.hh @@ -1,7 +1,7 @@ #ifndef CaloTemplateWFUtil_HH #define CaloTemplateWFUtil_HH -#include "Offline/Mu2eUtilities/inc/CaloPulseShape.hh" +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" #include #include @@ -11,8 +11,8 @@ namespace mu2e { class CaloTemplateWFUtil { public: - CaloTemplateWFUtil(const std::string& pulseFileName, const std::string& pulseHistName, - double minPeakAmplitude, double digiSampling, double minDTPeaks, int printLevel=-1); + CaloTemplateWFUtil(const CaloPulseUtil::Config& configPulseCache, double minPeakAmplitude, + double minDTPeaks, int printLevel=-1); void initialize (); void setXYVector (const std::vector& xvec, const std::vector& yvec); @@ -48,7 +48,7 @@ namespace mu2e { private: bool selectComponent(const std::vector& tempPar, const std::vector& tempErr, unsigned ip); - CaloPulseShape pulseCache_; + CaloPulseUtil pulseCache_; double minPeakAmplitude_; double minDTPeaks_; int fitStrategy_; diff --git a/CaloReco/src/CaloTemplateWFProcessor.cc b/CaloReco/src/CaloTemplateWFProcessor.cc index fc23dd57a3..ff72e1ab16 100644 --- a/CaloReco/src/CaloTemplateWFProcessor.cc +++ b/CaloReco/src/CaloTemplateWFProcessor.cc @@ -1,3 +1,4 @@ +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" #include "Offline/CaloReco/inc/CaloTemplateWFProcessor.hh" #include "Offline/CaloReco/inc/CaloTemplateWFUtil.hh" #include "art/Framework/Services/Registry/ServiceHandle.h" @@ -25,8 +26,7 @@ namespace mu2e { chiThreshold_ (config.chiThreshold()), refitLeadingEdge_(config.refitLeadingEdge()), diagLevel_ (config.diagLevel()), - fmutil_ (config.pulseFileName(),config.pulseHistName(),minPeakAmplitude_, - config.digiSampling(),minDTPeaks_,config.fitPrintLevel()), + fmutil_ (config.pulseCache(),minDTPeaks_,config.fitPrintLevel()), chi2_ (999.), ndf_ (-1), resAmp_ (), diff --git a/CaloReco/src/CaloTemplateWFUtil.cc b/CaloReco/src/CaloTemplateWFUtil.cc index 183770fc56..cf882f2049 100644 --- a/CaloReco/src/CaloTemplateWFUtil.cc +++ b/CaloReco/src/CaloTemplateWFUtil.cc @@ -1,5 +1,5 @@ #include "Offline/CaloReco/inc/CaloTemplateWFUtil.hh" -#include "Offline/Mu2eUtilities/inc/CaloPulseShape.hh" +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" #include "TMinuit.h" #include "TF1.h" @@ -22,7 +22,7 @@ namespace { unsigned npTot_(0),npFcn_(0),npBkg_(0),x0_(0),x1_(0); std::vector xvec_{},yvec_{}; - mu2e::CaloPulseShape* pulseCachePtr_=(nullptr); + mu2e::CaloPulseUtil* pulseCachePtr_=(nullptr); double logn(double x, double *par) {return par[0]*pulseCachePtr_->evaluate(x-par[1]); } @@ -57,9 +57,9 @@ namespace namespace mu2e { - CaloTemplateWFUtil::CaloTemplateWFUtil(const std::string& pulseFileName, const std::string& pulseHistName, - double minPeakAmplitude, double digiSampling, double minDTPeaks, int printLevel) : - pulseCache_(CaloPulseShape(pulseFileName, pulseHistName, digiSampling)), + CaloTemplateWFUtil::CaloTemplateWFUtil(const CaloPulseUtil::Config& configPulseCache, double minPeakAmplitude, + double minDTPeaks, int printLevel) : + pulseCache_(CaloPulseUtil(configPulseCache)), minPeakAmplitude_(minPeakAmplitude), minDTPeaks_(minDTPeaks), fitStrategy_(1), @@ -80,7 +80,7 @@ namespace mu2e { //----------------------------------------------------------------------------------------------------- - void CaloTemplateWFUtil::initialize () {pulseCache_.buildShapes();} + void CaloTemplateWFUtil::initialize () {pulseCache_.buildCache();} void CaloTemplateWFUtil::reset () {param_.clear(); paramErr_.clear(); nParTot_=0; npTot_ = 0;} void CaloTemplateWFUtil::setXYVector(const std::vector& xvec, const std::vector& yvec) {xvec_ = xvec; yvec_ = yvec; x0_=0; x1_ = xvec_.size();} void CaloTemplateWFUtil::setPar (const std::vector& par) {param_ = par; nParTot_ = npTot_ = par.size();} diff --git a/Mu2eUtilities/CMakeLists.txt b/Mu2eUtilities/CMakeLists.txt index 2131d1ec21..e3c1dd1579 100644 --- a/Mu2eUtilities/CMakeLists.txt +++ b/Mu2eUtilities/CMakeLists.txt @@ -2,7 +2,8 @@ cet_make_library( SOURCE src/BinnedSpectrum.cc src/BuildLinearFitMatrixSums.cc - src/CaloPulseShape.cc + src/CaloPulseUtil.cc + src/CaloNoiseUtil.cc src/checkSimParticleCollection.cc src/compressPdgId.cc src/ConversionSpectrum.cc diff --git a/Mu2eUtilities/inc/CaloNoiseUtil.hh b/Mu2eUtilities/inc/CaloNoiseUtil.hh new file mode 100644 index 0000000000..81cce59d9b --- /dev/null +++ b/Mu2eUtilities/inc/CaloNoiseUtil.hh @@ -0,0 +1,77 @@ +#ifndef CaloNoiseUtil_HH +#define CaloNoiseUtil_HH +// +// Cache and provide noise waveforms for readouts +// +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" +#include "art/Framework/Services/Optional/RandomNumberGenerator.h" +#include "Offline/SeedService/inc/SeedService.hh" + +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" + +#include "CLHEP/Random/RandPoissonQ.h" +#include "CLHEP/Random/RandGaussQ.h" +#include "CLHEP/Random/RandFlat.h" + +#include +#include +#include + + +namespace mu2e { + + class CaloNoiseUtil + { + public: + struct Config + { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + using CPG = CaloPulseUtil::Config; + fhicl::Table pulseCache { Name("pulseCache"), Comment("Pulse cache maker config") }; + fhicl::Atom generate { Name("generate"), Comment("Regenerate waveform (true) or use histogram (false)") }; + fhicl::Atom dumpGenerated { Name("dumpGenerated"), Comment("Dump generated waveform") }; + fhicl::Atom histoFileName { Name("histoFileName"), Comment("Calo noise histo file name") }; + fhicl::Atom histoPrefix { Name("histoPrefix"), Comment("Noise histogram prefix") }; + fhicl::Atom elecNphotPerNs { Name("elecNphotPerNs"), Comment("Electronics noise number of PE / ns ") }; + fhicl::Atom rinNphotPerNs { Name("rinNphotPerNs"), Comment("RIN noise number of PE / ns ") }; + fhicl::Atom darkNphotPerNs { Name("darkNphotPerNs"), Comment("SiPM Dark noise number of PE / ns ") }; + fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; + }; + + + CaloNoiseUtil(const Config& config, CLHEP::HepRandomEngine& engine); + + void prepare(int histoID, double peToADC); + std::span noiseSegment(int histoID, size_t istart, size_t ilength); + int pedestal(); + void printCache(); + void dumpNoise(const std::string& name, const std::vector& wave); + + + private: + void fillCache(int histoBaseID); + void generateCache(int histoID, double peToADC); + + bool generate_; + std::string fileName_; + std::string prefix_; + double digiSampling_; + double noiseRinDark_; + double noiseElec_; + double minPeakADC_; + CLHEP::RandPoissonQ randPoisson_; + CLHEP::RandGaussQ randGauss_; + CLHEP::RandFlat randFlat_; + CaloPulseUtil pulseCache_; + bool dumpGenerated_; + int histoBaseID_; + int pedestal_; + std::map> noiseMap_; + + static constexpr int base = 10000; + }; + +} +#endif diff --git a/Mu2eUtilities/inc/CaloPulseUtil.hh b/Mu2eUtilities/inc/CaloPulseUtil.hh new file mode 100644 index 0000000000..7575cbee4c --- /dev/null +++ b/Mu2eUtilities/inc/CaloPulseUtil.hh @@ -0,0 +1,62 @@ +#ifndef CaloPulseUtil_HH +#define CaloPulseUtil_HH + +// Calculate the values of the digitized pulse shape as a function of the hit time. +// The value stored are the integral of the waveform over the digitization bin width. +// +// The waveform "starting point" corresponds to the time the PE hit the readout. +// The definition is arbitrary, but it must be internally consistent! +// +// Shifting the waveform forward by a time dt is equivalent to shifting the time origin backward by dt +// - value at nbin-i correspond to waveform shifted by time nbin+i, +// - t0 value is located at bin nSteps_ +// - the t0 value correspond to the content of the digitized bin whose LOW EDGE is at time t0 +// +// 1) digitizedPulse(hitTime) returns a waveform with hitTime corresponding to low edge of first bin +// 2) evaluate(deltaTime) return value of digitized bin at a given time difference with peak time value +// + +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" + +#include +#include + +namespace mu2e { + + class CaloPulseUtil + { + public: + struct Config + { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom fileName { Name("fileName"), Comment("Pulse file name") }; + fhicl::Atom histName { Name("histName"), Comment("Pulse histogram name") }; + fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitizer sampling time (ns) ") }; + }; + + CaloPulseUtil(const Config& config); + CaloPulseUtil(const std::string& fileName, const std::string& histName, double digiSampling); + ~CaloPulseUtil() = default; + + void buildCache(); + + const std::vector& digitizedPulse (double hitTime) const; + double evaluate (double timeDifference) const; + double fromPeakToT0 (double timePeak) const; + void diag (bool fullDiag=false) const; + + private: + std::string fileName_; + std::string histName_; + int nSteps_; + double digiStep_; + int nBinShape_; + std::vector pulseVec_; + double deltaT_; + mutable std::vector digitizedPulse_; + }; + +} +#endif diff --git a/Mu2eUtilities/src/CaloNoiseUtil.cc b/Mu2eUtilities/src/CaloNoiseUtil.cc new file mode 100644 index 0000000000..fc35c64224 --- /dev/null +++ b/Mu2eUtilities/src/CaloNoiseUtil.cc @@ -0,0 +1,188 @@ +#include "cetlib_except/exception.h" +#include "art/Framework/Services/Optional/RandomNumberGenerator.h" +#include "Offline/SeedService/inc/SeedService.hh" +#include "Offline/Mu2eUtilities/inc/CaloNoiseUtil.hh" +#include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" + +#include "TFile.h" +#include "TH1F.h" +#include "TKey.h" +#include "TDirectory.h" + +#include +#include +#include +#include +#include + + +namespace mu2e { + + CaloNoiseUtil::CaloNoiseUtil(const Config& config, CLHEP::HepRandomEngine& engine) : + generate_ {config.generate()}, + fileName_ {config.histoFileName()}, + prefix_ {config.histoPrefix()}, + digiSampling_ {config.digiSampling()}, + noiseRinDark_ {config.rinNphotPerNs() + config.darkNphotPerNs()}, + noiseElec_ {config.elecNphotPerNs()}, + randPoisson_ {engine}, + randGauss_ {engine}, + randFlat_ {engine}, + pulseCache_ {config.pulseCache()}, + dumpGenerated_ {config.dumpGenerated()}, + histoBaseID_ {-1}, + pedestal_ {0}, + noiseMap_ {} + {} + + + //---------------------------------------------------------------------------------------------------------------------- + void CaloNoiseUtil::prepare(int histoID, double peToADC) + { + int baseID = histoID / base; + if (baseID == histoBaseID_) return; + if (generate_) generateCache(baseID, peToADC); + else fillCache(baseID); + } + + //---------------------------------------------------------------------------------------------------------------------- + void CaloNoiseUtil::fillCache(int histoBaseID) + { + // Cache is for all baseID, clear it + noiseMap_.clear(); + histoBaseID_ = histoBaseID; + + // Refill the cache with all histos sharing the same baseID + ConfigFileLookupPolicy resolveFullPath; + std::string fullFileName = resolveFullPath(fileName_); + + TFile file(fullFileName.c_str()); + if (!file.IsOpen()) throw cet::exception("NOISEREADER")<<"Filename"<ReadObj(); + if (!obj->InheritsFrom(TH1F::Class())) + continue; + + TH1F* histo = static_cast(obj); + std::string name = histo->GetName(); + + // Parse the integer after the prefix + std::string suffix = name.substr(prefix_.size()); + std::istringstream iss(suffix); + + int hid; + if (!(iss >> hid) || !iss.eof()) + throw cet::exception("NOISEREADER")<<"Hitsogram "<GetArray(); + noiseMap_[hid].assign(array + 1, array + histo->GetNbinsX() + 1); + + // estimate pedestal, take a single value for every waveform + // (pedestals will be stored somewhere else later) + + float sum = 0.0; + for (int i = 1; i <= histo->GetNbinsX(); ++i) sum += histo->GetBinContent(i); + pedestal_ = std::trunc(sum /histo->GetNbinsX() ); + } + } + + + //---------------------------------------------------------------------------------------------------------------------- + void CaloNoiseUtil::generateCache(int histoBaseID, double peToADC) + { + // Cache is for all baseID, clear it + noiseMap_.clear(); + histoBaseID_ = histoBaseID; + constexpr unsigned noiseSize{10000}; + std::vector waveform(noiseSize,0.0); + + pulseCache_.buildCache(); + const auto& pulse = pulseCache_.digitizedPulse(0.0); + const unsigned pulseSize = pulse.size(); + const unsigned bufferSize = int(0.75*pulseSize); + const double totalTime = (noiseSize+bufferSize)*digiSampling_; + const int noiseLevelPE = int(totalTime*noiseRinDark_); + + //Generate the radiation induced noise (RIN) + const int nPh = randPoisson_(noiseLevelPE); + for (int i=0;i CaloNoiseUtil::noiseSegment(int histoID, size_t istart, size_t ilength) + { + int baseID = histoID/base; + + if (baseID != histoBaseID_) { + if (generate_) throw cet::exception("CaloNoiseUtil") + << "noiseSegment called before prepare() for baseID " << histoID/base << "\n"; + else fillCache(baseID); + } + + auto iter = noiseMap_.find(histoID); + if (iter == noiseMap_.end()) + throw cet::exception("CALONOISEUTIL")<<"histoID "<second; + + if (ilength >= vec.size()) + throw cet::exception("CALONOISEUTIL")<<"noise length request too long\n"; + + size_t irandom = size_t(randFlat_.fire(0.,vec.size()-ilength)); + return std::span(vec.data() + irandom, ilength); + } + + //---------------------------------------------------------------------------------------------------------------------- + int CaloNoiseUtil::pedestal() {return pedestal_;} + + //---------------------------------------------------------------------------------------------------------------------- + void CaloNoiseUtil::printCache() + { + std::cout<<"CaloNoiseUtil cache\n"; + for (const auto& kv : noiseMap_) std::cout<<"Histo id "<& wave) + { + TFile outfile(fname.c_str(), "RECREATE"); + + TH1F h("histo_0","histo_0", wave.size(), 0, wave.size()); + for (size_t i = 0; i < wave.size(); ++i) h.SetBinContent(i+1, wave[i]); + + h.Write(); + outfile.Close(); + std::cout<<"CaloNoiseUtil written waveform in "< +#include +#include + + +namespace mu2e { + + CaloPulseUtil::CaloPulseUtil(const Config& config) : + fileName_(config.fileName()), + histName_(config.histName()), + nSteps_(100), + digiStep_(config.digiSampling()/double(nSteps_)), + nBinShape_(0), + pulseVec_(), + deltaT_(0.), + digitizedPulse_() + {} + CaloPulseUtil::CaloPulseUtil(const std::string& fileName, const std::string& histName, double digiSampling) : + fileName_(fileName), + histName_(histName), + nSteps_(100), + digiStep_(digiSampling/double(nSteps_)), + nBinShape_(0), + pulseVec_(), + deltaT_(0.), + digitizedPulse_() + {} + + //---------------------------------------------------------------------------------------------------------------------- + void CaloPulseUtil::buildCache() + { + pulseVec_.clear(); + + ConfigFileLookupPolicy resolveFullPath; + std::string fullFileName = resolveFullPath(fileName_); + + std::unique_ptr pshape(nullptr); + TFile pulseFile(fullFileName.c_str()); + if (pulseFile.IsOpen()) pshape.reset((TH1F*) pulseFile.Get(histName_.c_str())); + if (!pshape) throw cet::exception("CATEGORY")<<"CaloPulseUtil:: Hitsogram "<SetDirectory(0); + pulseFile.Close(); + + + // Adjust binning to match digitizer sampling period, shift to zero and normalize + int nbins = int((pshape->GetXaxis()->GetXmax()-pshape->GetXaxis()->GetXmin())/digiStep_); + TH1F pulseShape("ps","ps", nbins, 0.0, pshape->GetXaxis()->GetXmax()-pshape->GetXaxis()->GetXmin()); + for (int i=1;i<=nbins;++i) pulseShape.SetBinContent(i,pshape->Interpolate(pulseShape.GetBinCenter(i))); + pulseShape.Scale(1.0/pulseShape.GetMaximum(),"nosw2"); + + // Cache histogram content into vector and shift waveform (see note), + // calculate the number of bins for the digitized waveform + for (int j=1;j<=(nbins+nSteps_);++j)pulseVec_.push_back((j>nSteps_) ? pulseShape.GetBinContent(j-nSteps_) : 0.0); + nBinShape_ = int(nbins/nSteps_); + digitizedPulse_ = std::vector(nBinShape_,0); + + deltaT_ = 0.0; + // find difference between peak time and t0 for digitized waveform. + for (int i=1;i& CaloPulseUtil::digitizedPulse(double hitTime) const + { + int shiftBin = nSteps_ - int(hitTime/digiStep_)%nSteps_; + for (int i=0;i= int(pulseVec_.size()-1)) return 0.0; + double t0bin = (ibin-nSteps_)*digiStep_; //t0 is located at nSteps_ + return (pulseVec_[ibin+1]-pulseVec_[ibin])/digiStep_*(t-t0bin)+pulseVec_[ibin]; + } + + //---------------------------------------------------------------------------- + double CaloPulseUtil::fromPeakToT0(double timePeak) const + { + return timePeak-deltaT_-0.5*digiStep_*nSteps_; + } + + //-------------------------------- + void CaloPulseUtil::diag(bool fullDiag) const + { + std::cout<<"Number of digi bins "< Date: Thu, 6 Aug 2026 22:53:57 -0500 Subject: [PATCH 4/7] More refactorign and CMake fix --- CaloMC/CMakeLists.txt | 3 + CaloMC/fcl/prolog.fcl | 6 +- CaloMC/inc/CaloPhotonPropagation.hh | 3 +- CaloMC/src/CaloClusterTruthMatch_module.cc | 6 +- CaloMC/src/CaloDigiMaker_module.cc | 55 +++++------ CaloMC/src/CaloHitTruthMatch_module.cc | 24 ++--- CaloMC/src/CaloPhotonPropagation.cc | 2 +- CaloMC/src/CaloShowerROMaker_module.cc | 19 ++-- CaloMC/src/CaloShowerStepMaker_module.cc | 8 +- CaloMC/src/ShowerStepUtil.cc | 26 +++-- CaloReco/fcl/common.fcl | 2 +- Mu2eUtilities/inc/CaloPulseShape.hh | 51 ---------- Mu2eUtilities/src/CaloPulseShape.cc | 105 --------------------- Mu2eUtilities/src/CaloPulseUtil.cc | 9 +- 14 files changed, 75 insertions(+), 244 deletions(-) delete mode 100644 Mu2eUtilities/inc/CaloPulseShape.hh delete mode 100644 Mu2eUtilities/src/CaloPulseShape.cc diff --git a/CaloMC/CMakeLists.txt b/CaloMC/CMakeLists.txt index fdac25c801..fa7dd8373f 100644 --- a/CaloMC/CMakeLists.txt +++ b/CaloMC/CMakeLists.txt @@ -6,6 +6,7 @@ cet_make_library( src/CaloDigiWrapperCollection.cc LIBRARIES PUBLIC + Offline::ConfigTools Offline::CalorimeterGeom Offline::GeometryService Offline::Mu2eUtilities @@ -85,6 +86,8 @@ cet_build_plugin(CaloShowerUpdater art::module configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/prolog.fcl ${CURRENT_BINARY_DIR} fcl/prolog.fcl) +install(DIRECTORY data DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/Offline/CaloReco) + install_source(SUBDIRS src) install_headers(USE_PROJECT_NAME SUBDIRS inc) install_fhicl(SUBDIRS fcl SUBDIRNAME Offline/CaloMC/fcl) diff --git a/CaloMC/fcl/prolog.fcl b/CaloMC/fcl/prolog.fcl index 698424b61c..4574661b67 100644 --- a/CaloMC/fcl/prolog.fcl +++ b/CaloMC/fcl/prolog.fcl @@ -45,13 +45,13 @@ CaloMC : { caloShowerROCollection : CaloShowerROMaker eventWindowMarker : EWMProducer protonBunchTimeMC : EWMProducer - pulseCache : { @table::CaloPulseCache } - noiseCache : { @table::CaloNoiseCache } + pulseCache : @local::CaloPulseCache + noiseCache : @local::CaloNoiseCache + addNoise : true digitizationStart : @local::HitMakerDigitizationStart digitizationEnd : @local::HitMakerDigitizationEnd digiSampling : @local::HitMakerDigiSampling bufferDigi : 16 - addNoise : true nBits : 12 minPeakADC : @local::HitMakerMinPeakADC nBinsPeak : 2 diff --git a/CaloMC/inc/CaloPhotonPropagation.hh b/CaloMC/inc/CaloPhotonPropagation.hh index 104fe26f86..9b3e4ab322 100644 --- a/CaloMC/inc/CaloPhotonPropagation.hh +++ b/CaloMC/inc/CaloPhotonPropagation.hh @@ -16,7 +16,8 @@ namespace mu2e { class CaloPhotonPropagation { public: - CaloPhotonPropagation(const std::string& fileName, const std::string& histName, CLHEP::HepRandomEngine& engine); + CaloPhotonPropagation(const std::string& fileName, const std::string& histName, + CLHEP::HepRandomEngine& engine); void buildTable(); float propTimeSimu(float z); // sampled from the CDF (consumes a random number) diff --git a/CaloMC/src/CaloClusterTruthMatch_module.cc b/CaloMC/src/CaloClusterTruthMatch_module.cc index 599265b011..b8057dcabb 100644 --- a/CaloMC/src/CaloClusterTruthMatch_module.cc +++ b/CaloMC/src/CaloClusterTruthMatch_module.cc @@ -90,8 +90,7 @@ namespace mu2e { caloClusterMCs.reserve(caloClusters.size()); - for (std::size_t idx=0; idx(caloClusterHandle,idx); @@ -100,8 +99,7 @@ namespace mu2e { // gather the CaloHitMC of every hit in this cluster that carries MC truth std::vector> digis; - for (const auto& hitPtr : cluster.caloHitsPtrVector()) - { + for (const auto& hitPtr : cluster.caloHitsPtrVector()) { const auto it = hitToMC.find(hitPtr); if (it == hitToMC.end()) continue; diff --git a/CaloMC/src/CaloDigiMaker_module.cc b/CaloMC/src/CaloDigiMaker_module.cc index 01376eecab..8a1092c1b0 100644 --- a/CaloMC/src/CaloDigiMaker_module.cc +++ b/CaloMC/src/CaloDigiMaker_module.cc @@ -13,19 +13,19 @@ #include "art_root_io/TFileService.h" #include "art_root_io/TFileDirectory.h" -#include "Offline/CaloConditions/inc/CalSimParams.hh" -#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" -#include "Offline/Mu2eUtilities/inc/CaloNoiseUtil.hh" #include "Offline/CalorimeterGeom/inc/Calorimeter.hh" -#include "Offline/ProditionsService/inc/ProditionsHandle.hh" -#include "Offline/DAQConditions/inc/EventTiming.hh" +#include "Offline/CaloConditions/inc/CalSimParams.hh" +#include "Offline/DataProducts/inc/CaloSiPMId.hh" #include "Offline/DataProducts/inc/EventWindowMarker.hh" +#include "Offline/DAQConditions/inc/EventTiming.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/MCDataProducts/inc/CaloShowerRO.hh" +#include "Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh" +#include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" +#include "Offline/Mu2eUtilities/inc/CaloNoiseUtil.hh" +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" #include "Offline/RecoDataProducts/inc/CaloDigi.hh" #include "Offline/SeedService/inc/SeedService.hh" -#include "Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh" -#include "Offline/DataProducts/inc/CaloSiPMId.hh" #include "CLHEP/Vector/ThreeVector.h" #include "CLHEP/Random/RandPoissonQ.h" @@ -182,18 +182,16 @@ namespace mu2e { if (calorimeter_->nCrystals()<1 || calorimeter_->G4Info().get("nSiPMPerCrystal")<1) return; int waveformSize = (digitizationEnd_ - digitizationStart_ + startTimeBuffer_) / digiSampling_; - if (ewMarker.spillType() != EventWindowMarker::SpillType::onspill) - { + if (ewMarker.spillType() != EventWindowMarker::SpillType::onspill) { waveformSize = (ewMarker.eventLength() - digitizationStart_ + startTimeBuffer_) / digiSampling_; } int nWaveforms = calorimeter_->nCrystals()*calorimeter_->G4Info().get("nSiPMPerCrystal"); - if (waveformSize<1) throw cet::exception("Rethrow")<< "[CaloMC/CaloDigiMaker] digitization size too short " << std::endl; + if (waveformSize<1) throw cet::exception("CALODIGIMAKER")<< "Digitization size too short " << std::endl; bool resetWaveform(false); std::vector waveform(waveformSize,0.0); - for (int iRO=0;iRO& waveform, double scaleFactor) { - const int NoiseWFID(0); // will need this from proditions later; + const int NoiseWFID(0); // will get this from proditions later; - const double minAmplitude(0.01); + const double minAmplitude(2); noiseSampler_.prepare(NoiseWFID, scaleFactor); size_t timeSample(0); @@ -220,8 +218,7 @@ namespace mu2e { hitStarts.reserve(16);hitStops.reserve(16); // First, find the ranges in the waveform with non-zero bins. - while (timeSample < waveform.size()) - { + while (timeSample < waveform.size()) { if (waveform[timeSample] < minAmplitude) {++timeSample; continue;} size_t sampleStart = (timeSample > bufferDigi_) ? timeSample - bufferDigi_ : 0; @@ -235,8 +232,7 @@ namespace mu2e { // ranges might overlap and need to be concatenated if this is the case size_t iprev(0),ic(1); - while (ic < hitStarts.size()) - { + while (ic < hitStarts.size()) { if (hitStops[iprev] >= hitStarts[ic]) {hitStops[iprev]=hitStops[ic]; hitStarts[ic]=hitStops[ic]=waveform.size()+1;} else {iprev = ic;} ++ic; @@ -246,8 +242,7 @@ namespace mu2e { hitStops.erase( std::remove_if(hitStops.begin(), hitStops.end(), pred),hitStops.end()); //Now take a random part of the noise waveform and add it to the waveform content - for (size_t ihit=0; ihit wfsample{}; wfsample.reserve(sampleStop-sampleStart); @@ -329,8 +321,7 @@ namespace mu2e { void CaloDigiMaker::extract(const std::vector& wf, std::vector& starts, std::vector& stops) const { size_t timeSample(nBinsPeak_+bufferDigi_); - while (timeSample+nBinsPeak_ < wf.size()) - { + while (timeSample+nBinsPeak_ < wf.size()){ // find starting point if (wf[timeSample] < minPeakADC_) {++timeSample; continue;} @@ -351,11 +342,9 @@ namespace mu2e { timeSample = sampleStop+1; } - // Concatenate peaks and remove unused values (flag value to remove past wf.size() since the latter is a legitimate value) size_t iprev(0), icurrent(1); - while (icurrent < starts.size()) - { + while (icurrent < starts.size()){ if (stops[iprev] >= starts[icurrent]) {stops[iprev]=stops[icurrent]; starts[icurrent]=stops[icurrent]=wf.size()+1;} else {iprev = icurrent;} ++icurrent; @@ -376,8 +365,6 @@ namespace mu2e { return ADCPerMeV / pePerMeV; } - - //------------------------------------------------------------------------------------------------------------------- void CaloDigiMaker::diag0(unsigned iSiPM, const std::vector& wf) { diff --git a/CaloMC/src/CaloHitTruthMatch_module.cc b/CaloMC/src/CaloHitTruthMatch_module.cc index 9a9745598e..a78c1cec58 100644 --- a/CaloMC/src/CaloHitTruthMatch_module.cc +++ b/CaloMC/src/CaloHitTruthMatch_module.cc @@ -150,17 +150,14 @@ namespace mu2e { // O(1) shower-pointer -> index lookup for the detailed-MC association (was an O(n) scan per match) std::unordered_map showerIndex; - if (fillDetailedMC_) - { + if (fillDetailedMC_) { showerIndex.reserve(caloShowerSims.size()); for (std::size_t i=0; i edeps; auto showerIt = sortedSims.begin(); while (showerIt != sortedSims.end() && (*showerIt)->time() < hit.time() - deltaTimeMinus_) ++showerIt; - while (showerIt != sortedSims.end() && (*showerIt)->time() < hit.time() + deltaTimePlus) - { + while (showerIt != sortedSims.end() && (*showerIt)->time() < hit.time() + deltaTimePlus) { hitIsMatched = true; const CaloShowerSim* showerSim = *showerIt; fillEdeps(primaryParticle, edeps, showerSim); - if (fillDetailedMC_) - { + if (fillDetailedMC_) { const auto showerSimPtr = art::Ptr(caloShowerSimHandle, showerIndex.at(showerSim)); caloShowerTruthMatch.addSingle(hitPtr, showerSim->sim(), showerSimPtr); } @@ -228,18 +223,15 @@ namespace mu2e { auto it = std::find_if(edeps.begin(), edeps.end(), [&](const CaloEDepMC& e){return e.sim() == showerSim->sim();}); - if (it != edeps.end()) - { + if (it != edeps.end()) { it->addEDep (showerSim->energyDep()); it->addEDepG4(showerSim->energyDepG4()); it->addTime (showerSim->time()); it->addMom (showerSim->momentumIn()); } - else - { + else{ MCRelationship mcrel; - for (const auto& spp : primaryParticle.primarySimParticles()) - { + for (const auto& spp : primaryParticle.primarySimParticles()) { MCRelationship mcr(spp,showerSim->sim()); if (mcr > mcrel) mcrel = mcr; } diff --git a/CaloMC/src/CaloPhotonPropagation.cc b/CaloMC/src/CaloPhotonPropagation.cc index e92cd2d34e..93e949736d 100644 --- a/CaloMC/src/CaloPhotonPropagation.cc +++ b/CaloMC/src/CaloPhotonPropagation.cc @@ -62,7 +62,7 @@ namespace mu2e { } const Calorimeter& cal = *(GeomHandle()); - lightSpeed_ = 300.0f / cal.G4Info().get("refractiveIndex"); // mm/ns + lightSpeed_ = 300.0f / cal.G4Info().get("refractiveIndex"); // mm/ns } diff --git a/CaloMC/src/CaloShowerROMaker_module.cc b/CaloMC/src/CaloShowerROMaker_module.cc index 09d3a55bb5..ac7ec79420 100644 --- a/CaloMC/src/CaloShowerROMaker_module.cc +++ b/CaloMC/src/CaloShowerROMaker_module.cc @@ -208,11 +208,10 @@ namespace mu2e { //----------------------------------------------------------------------- // store corrected energy deposits for each readout - for (const auto& showerHandle : crystalShowerHandles) - { + for (const auto& showerHandle : crystalShowerHandles) { const CaloShowerStepCollection& caloShowerSteps(*showerHandle); - for (auto istep = caloShowerSteps.begin(); istep != caloShowerSteps.end(); ++istep) - { + + for (auto istep = caloShowerSteps.begin(); istep != caloShowerSteps.end(); ++istep) { const CaloShowerStep& step = *istep; // see doc-db for calo folding description. Note pbtmc.pbtime_ is NEGATIVE! @@ -237,10 +236,9 @@ namespace mu2e { if (LRUCorrection_) edep_corr = LRUCorrection(posZ/crystalLength, edep_corr, lru); // Generate individual PEs and their arrival times - for (int i=0; i 2) - { + if (diagLevel_ > 2) { std::cout<<"[CaloShowerROMaker] SiPMID:"< 1) { std::cout<<"CaloShowerStepMaker summary"<hasParent()) { const auto alreadyInspected = simToAncestorMap.find(sim); @@ -283,10 +282,7 @@ namespace mu2e { for (const StepPointMC* step : steps) { const CLHEP::Hep3Vector pos = cal.mu2eToCrystal(volId, step->position()); - - // clamp: a step at/beyond the back face would otherwise index past the last slice const unsigned idx = std::min(unsigned(std::max(1e-6, pos.z())/zSliceSize_), numZSlices_-1); - if (buffer.entries(idx)>0 && (step->time()-buffer.t0(idx) > deltaTime_)) { if (diagLevel_ > 2) {std::cout<<"[CaloShowerStepMaker::compressSteps] inserted "; buffer.printBucket(idx);} diff --git a/CaloMC/src/ShowerStepUtil.cc b/CaloMC/src/ShowerStepUtil.cc index 6904249e83..7d8168436a 100644 --- a/CaloMC/src/ShowerStepUtil.cc +++ b/CaloMC/src/ShowerStepUtil.cc @@ -11,9 +11,12 @@ namespace mu2e { - void ShowerStepUtil::add(unsigned i, double eDepG4, double eDepVis, double time, double momentum, const CLHEP::Hep3Vector& pos) + void ShowerStepUtil::add(unsigned i, double eDepG4, double eDepVis, double time, + double momentum, const CLHEP::Hep3Vector& pos) { - if (i > imax_) throw cet::exception("Rethrow")<< "[CaloMC/ShowerStepUtil] Index out of bound " << i << std::endl; + if (i > imax_) throw cet::exception("Rethrow") + << "[CaloMC/ShowerStepUtil] Index out of bound " + << i << "\n"; //init buffer if needed if (n_[i]==0) {pIn_[i] = momentum; t0_[i] = time;} @@ -34,7 +37,10 @@ namespace mu2e { void ShowerStepUtil::reset(unsigned i) { - if (i > imax_) throw cet::exception("Rethrow")<< "[CaloMC/ShowerStepUtil] Index out of bound " << i << std::endl; + if (i > imax_) throw cet::exception("Rethrow") + << "[CaloMC/ShowerStepUtil] Index out of bound " + << i << "\n"; + n_[i] = 0; eDepG4_[i] = 0; eDepVis_[i] = 0; @@ -45,7 +51,9 @@ namespace mu2e { CLHEP::Hep3Vector& ShowerStepUtil::pos(unsigned i) { - if (i > imax_) throw cet::exception("Rethrow")<< "[CaloMC/ShowerStepUtil] Index out of bound " << i << std::endl; + if (i > imax_) throw cet::exception("Rethrow") + << "[CaloMC/ShowerStepUtil] Index out of bound " + << i << "\n"; pos_[0] = x_[i]/w_[i]; pos_[1] = y_[i]/w_[i]; @@ -56,9 +64,13 @@ namespace mu2e { void ShowerStepUtil::printBucket(unsigned i) { - if (i > imax_) throw cet::exception("Rethrow")<< "[CaloMC/ShowerStepUtil] Index out of bound " << i << std::endl; - std::cout<<"Entries= "< -#include - -namespace mu2e { - - class CaloPulseShape - { - public: - CaloPulseShape(const std::string& fileName, const std::string& histName, double digiSampling); - ~CaloPulseShape() {}; - - void buildShapes(); - - const std::vector& digitizedPulse (double hitTime) const; - double evaluate (double timeDifference) const; - double fromPeakToT0 (double timePeak) const; - void diag (bool fullDiag=false) const; - - private: - std::string fileName_; - std::string histName_; - int nSteps_; - double digiStep_; - int nBinShape_; - std::vector pulseVec_; - double deltaT_; - mutable std::vector digitizedPulse_; - }; - -} -#endif diff --git a/Mu2eUtilities/src/CaloPulseShape.cc b/Mu2eUtilities/src/CaloPulseShape.cc deleted file mode 100644 index 118b665357..0000000000 --- a/Mu2eUtilities/src/CaloPulseShape.cc +++ /dev/null @@ -1,105 +0,0 @@ -#include "cetlib_except/exception.h" -#include "Offline/Mu2eUtilities/inc/CaloPulseShape.hh" -#include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" - -#include "TFile.h" -#include "TH2F.h" - -#include -#include -#include - - -namespace mu2e { - - CaloPulseShape::CaloPulseShape(const std::string& fileName, const std::string& histName, double digiSampling) : - fileName_(fileName), - histName_(histName), - nSteps_(100), - digiStep_(digiSampling/double(nSteps_)), - nBinShape_(0), - pulseVec_(), - deltaT_(0.), - digitizedPulse_() - {} - - //---------------------------------------------------------------------------------------------------------------------- - void CaloPulseShape::buildShapes() - { - - pulseVec_.clear(); - - ConfigFileLookupPolicy resolveFullPath; - std::string fullFileName = resolveFullPath(fileName_); - - std::unique_ptr pshape(nullptr); - TFile pulseFile(fullFileName.c_str()); - if (pulseFile.IsOpen()) pshape.reset((TH1F*) pulseFile.Get(histName_.c_str())); - if (!pshape) throw cet::exception("CATEGORY")<<"CaloPulseShape:: Hitsogram "<SetDirectory(0); - pulseFile.Close(); - - - // Adjust binning to match digitizer sampling period, shift to zero and normalize - int nbins = int((pshape->GetXaxis()->GetXmax()-pshape->GetXaxis()->GetXmin())/digiStep_); - TH1F pulseShape("ps","ps", nbins, 0.0, pshape->GetXaxis()->GetXmax()-pshape->GetXaxis()->GetXmin()); - for (int i=1;i<=nbins;++i) pulseShape.SetBinContent(i,pshape->Interpolate(pulseShape.GetBinCenter(i))); - pulseShape.Scale(1.0/pulseShape.GetMaximum(),"nosw2"); - - // Cache histogram content into vector and shift waveform (see note), - // calculate the number of bins for the digitized waveform - for (int j=1;j<=(nbins+nSteps_);++j)pulseVec_.push_back((j>nSteps_) ? pulseShape.GetBinContent(j-nSteps_) : 0.0); - nBinShape_ = int(nbins/nSteps_); - digitizedPulse_ = std::vector(nBinShape_,0); - - deltaT_ = 0.0; - // find difference between peak time and t0 for digitized waveform. - for (int i=1;i& CaloPulseShape::digitizedPulse(double hitTime) const - { - int shiftBin = nSteps_ - int(hitTime/digiStep_)%nSteps_; - for (int i=0;i= int(pulseVec_.size()-1)) return 0.0; - double t0bin = (ibin-nSteps_)*digiStep_; //t0 is located at nSteps_ - return (pulseVec_[ibin+1]-pulseVec_[ibin])/digiStep_*(t-t0bin)+pulseVec_[ibin]; - } - - //---------------------------------------------------------------------------- - double CaloPulseShape::fromPeakToT0(double timePeak) const - { - return timePeak-deltaT_-0.5*digiStep_*nSteps_; - } - - //-------------------------------- - void CaloPulseShape::diag(bool fullDiag) const - { - std::cout<<"Number of digi bins "< Date: Fri, 7 Aug 2026 00:22:25 -0500 Subject: [PATCH 5/7] Added missing dir --- CaloReco/data/ReadoutNoise.root | Bin 0 -> 38082 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 CaloReco/data/ReadoutNoise.root diff --git a/CaloReco/data/ReadoutNoise.root b/CaloReco/data/ReadoutNoise.root new file mode 100644 index 0000000000000000000000000000000000000000..57298869e7ccedb14cb747d4647c2bc475335f68 GIT binary patch literal 38082 zcma&MRZtvEv@M(v+}+*XU4uIWC%8j!f_rdx2=4Cg4ub`Ech?zw(7E3^RsX~PbZ=L! z?%mbB_FmPGyVq|gC$}$Oo=(1e`C{?qOX|>ni2g6;{|{gPgD~`eljWB$kz`-KAR~VH zx}5v*6dHNo6|EhiY0;4JO%ia>%(dn0~ zr7=6}@Bb0=|KELq`5(@Iko)ht>3iZ`D(RM~e7jmiNW$Yitk958~_g8^x_4c&W&}o0};#g01as!w{N5 z3WfiIZLx!!r;+XLkG{91&wN*>gvPDtRgctJ%v zk#^4ddHn_nXjeIoI?f#7JTP_7jOXOM{(P^ebHDy;aCt+HS8NXHtbOd?+3Rl`;p#ew zJF@(IqZL@X%b$ka8r_Kj8x9~{!w{c#NV8<)CZtcG_ic6s;>Y%q@2UF|xV;Z~Ci;Zz z7V^){emcqSi;+F8cxn$OTmx)TnLcvvzcP1T00-|ik4P=PhWMd zvc&OI-I6v1^h{jXsPD{>Ht2i%)z8?Ie6Njm;YHTDvy~v_Qc&&~nca{ybZPqASe}35 z&sU_)G96un_a}cV@|aKqOnE-6mk%L1L0<6cl70t|&7mVN{PW~C zXxVKfP|etn=#?@p0brhr4MDmGb@?PVEIY|h35oAK1EsFAq-~uHF|`7d=j($@g4xwe z4wK#1(;Z@-~7O#gjc;m!n2Ao6Db5FxQ-zn3{5G`|{YMEUgC8J23? zJ^q?19;JZsxSn9M)}4DH&dOaDP{Pq19o0U=$8pseHmR*(8c5<=wcQ4+xbH$N)%n|#iS(z(e!0KxpK7^_<6_59Sr1!L&BW$i@ z&b38A2Vk!bzx1u0Ko2I}tT94j%DaR=x2F>ONYGecz3f-XJ%)u05)3Hy4kW>8L;$cm z`T4jmrRxnAS!eLV3(~|rL~-fKY(c|zqQ^HTRW=JeWu&C^L~S-nmm=y3yYH)=sAs5f z{;FigX_dU~B5-@b6y@M0|1cQ>ZjB$YlOp8cBf9k}j6O{!p-v;AJy_@E%zr_DM_*W+ zReP@vPc0lNsYp*UNRwFVehZBr8-D+b0s03w+f)ke@2#?A#TQJ(ga9rJYd_|qJeX>< zzlwLtV%R@~9Bj5>KIfEobQ{tEn&T4*?MiN7`Db5(x@8#3wp&~X_(KZbruqF=WCRD1 zCvR{0M~_S|8CM2VefjZFM;Pj{W}4ESeHUVlYw>1&d4i@?2Rc`obI5G=Jqk%(%V)iq z@DDs#;M;rFy%Z>;a?HYZYdCvkcPMSPbku$MlSyl3x@4eWU`Q&q45Z^`_+>mczH7Y* z*|ygd`#^jyDE(}NEo4X5UMvLsEx^^VRd~5nAHP^$^t0r&@|E>AF&4~JF>;)>iHEPL zQEF+=kG@OL8PlB93pJ{{H|fPp1h}VoJr`T7W?+eMLG7mLdF4h)V!GfDiy=b4ISx~t z<61f5LOTE>?qhDK`+gOIcCDh9uO~zc6t9uR?wC*rUX%<%+VttCGNfEo=n(tJ3imF) zFCclOaBwGW02M2*$_MBr`D9prl)i?txpS!6RCp8&O$PnDM5|=R+Lg{Z70&Dj6e%2J zdLncL+uzCh`AX=e ztF_SGYrvR?k;;pluSExP{fqAPNtM)06&2$v`Nb)U z2_gwt(KZA4-1{2Ktd;L{81;vATa)=chB*85!=10u$hqd3i}`)!wH`I)rNJ!cTIHzG z7^`f-6YQ4Zcz7UBC8a6a6h|9A3M6Q{#vN`0O5F6oHvPPLc)28opv_bf$;>w@ke+J z^lk-r)arx!n$U*aZ`V0d*of?e&B_k8716S=_QIF6XG5PT( z4QPTa%q}xN(`(d}8at`2&@WAO_7v3+^r9yR@{(WN`}?+`th@5x=l$P_h$B_nDY}59 zL%Yv3!*=+Et=3VUrQ!o|z6bAUqxFIbXCI9HU)VvZv4vdc^X+U(Cc% zui4YyU%IMa3MG4dATH~7Ypg$SyWi32$zPsMcBRv zj_&E%L_gdM30J*4qeTx2cp?=d3GLtO5*Ir z>B=gV_%^Dag^Fu%y?#`4*xHKMdm~wWBscmwhCFtntUV?VZph>G@2i0<8BXrGT{403 z&rFtnT1D_GsH2c~ItdGX_+@WtsnbU=gWmp0bvzei2CHzw^r)T9gq!u;%3y`6s9e7W z*V$78RdirF*BXD^*Uce(uUa%j$c)7RV3S@+L&hDf4(kvwG$-irYM5$TVa{snKSsVf z^0juJ*M>f|;A1%J^A5LyO)Oqp9@Ca)5&yj&Gn@g|hJLK*v=?&`$lxWs{b#q!$IsD- z#!m`PeP)LI`~dnO;F0F00oJJXF*IZR(EFJZHBjtuzv0WNI9{Xi>sMs7J+aBK0kCSv zb?X4+&NaA6558|9?QF-B{)d6m(0CFxcVJ0}?>Fxe^YI^jG68=mI>kqiUP-9lUfh~i z9cbJPIgLE+BXwN~)5Cvd#CZ<8{1)6mWUbf{P!{X`95?0{AH#cn*s7v=@n|kT8vZ9d zv3ervw)&dRx>`byEF5CLJNTgcoEcyT;OJ^5;CXZA&YEI(%nU!=g34h@IP8S5CnFjv zjl}}&iX}z z=|1jgb|>FWEQ$Q3rx=gqb6XBb^r!2xIQjd8s7l1HU8f?$vTLKqa5> zW6}7)Z=0&b0oRmC8&hK;WN$-Qf+{bm$&zDY0oe>sgtF9zW+oCYEsC)$R(LVE&wFq_Uje>0gPhm*B{6=_BrF zLJZf`;0Fv%=6B`e8%Jz?fy92^yp7Cu`kXq6?QJU=n@lx7#k8i%U#XdG9RChU8Qd!^ zh;yV?5^62HcWFFi2);^&AY|~s&ORf6K`@@hn8B1*hwW|{YQ1rhtzyy#MT!7YQ5EpK z%>kS(W@n|CXs-`upETAiaS%_y#s1xS>A=d&nLGH~rv^JoW`bootF9%^o(~uN?480^ zRzHT8DCur+AK)E7V$~LFD5uKhF+bGe6S0w$(jfR;ZUGAbeT^=-xhUxj<{ z?3F3@!M9jQg)HQJlii;+Ca8^|YL?fe2!c3Gg<8Ry6~8k@>*kR&)UuDQCor)_+dE>} zj-8%jVqN&h*LU9f)Of@)r-Fbk`L^C)eyY*HfNdE=|0(7wQ8YB~j`}04=Bzp*%9}TW zy3`tAXio`@893pyN~}%$9vnp4c_R|FA);?GFKuDlrIzTG)^pfQAG;;lDuG31zMX9cFXIJ!#xFO;MuDakpD`cX~bCxog$LkDR9U)XogQTCY5ZhIX$FD*cfyK#`=j_-dB(*dV4QtD*Am0!Al zd_#LePZca!y+=ePy>Sve%i-I@CPvvHOR|uS)^PrUxTz-VoXRD!+pk3I$u~8t#cBVK zNRyRV)@AKGwj-MCHOS#8T6Pg~X#pERG`m-mjf@fymo=;eE9zjfHbER)UBk|K6=Jza z7z^|pu>DXJzV^VDXleV0mae*e-UangPZd<6b2-Mk zd*Vx?_RP5on4$LBrk2|d@+yr8b^D|?vX=NO?McapTug19GZ2&u@AONx%b2h*E^vb1 zRKZyaVLlq~)-{V_rjVGTtKRw1mg4F{0~3l4<6e78YV3Nv-z)4vWtDzbK<;qQKh}FG z#l?zfMEHDQiU3hblXQ~#{(K5wBt)KM(WVi9G9#Y7%#%%z`Pzp4zQ}PMAuD^qhII;) zUoJhoAPEX-yIjRg`K4{7h*kTC@6SwP$d97<8|vFrw<0HF%Ex=SoNf{7r)?|0cQ^Q) z57>#wWiNvIQHl$cm2Ngu2h+Ig0;h~E3=l^UjKPp7m}(~qk4)aG8Ictz7>qfN&IO+^Iqs>PD`-zCil`5k=(vXkHO#1`aL`R zReoufvD8tY@a3fSh#)4Rl_OgsCH*KzhoapWTl(8?x!Fw>Js=-88%z0K&p=9bT9k93 z%=El5p7QIZ>Y4~=d_fwE4`gKu{66rGY%khQNgOrB2hW+g(jr3q~G<4sQ^JW zvGvgPK4ks5!2?YT#83Bl>Y@lme0qxv?sJ=cdLmkF-?4j%VQfzZ(WlOq=5K3Z$Niqa zXpNz;J)G0!L;to+{lMoh^VmVC2ef1FDNHy|`mqG#L6C|*QgFmTpXo2#{D8f^lwP>^ zihr!tWb@!19saDCT83NDOb=IVZrIqc9sY7?tE6)UQ0R~T0wM64sT;QZ@Ol?&-B%a; zJ_X&_F+GDV;>ahyTI8}r2QFor=MQ)ZJATkp8=f)c?WB6e^C@s-u_}~vkHvi$yzs*m zbnl+q5kkBVgbW7g5ch=_p@!uQ`%1F(Ro9sxsVEa`#5UD8tg8zNL(x`+u zhDJ=&MCQfMR-N`ozZF_*R3HJ6H+IDs!hW-et%p;cNe@X`5lYC=s~it?lC;@w%DT#_C*^tByr^}@B_*Ye*oJ?*;fL#mx&|9zy~nNwo-~*$9C^6X;HL}; zzdQJ#K_z)Yx6HY((SWrXN=I7)aUtUIrAE5XNP@z&TPgmP1Jgb%RevgBLP7+7jnL)8 z7yO945!QRer=S#$oK`5}-_AgS*p_X9;Ct3Xk@399N!f|Bj);o-`iWehq-gdz}3WGmxe?vo=Bh60rD?ujc| zBrMteRI0gG3+mLn%?HK>z7f8Qk`67?5$6LOGjKkoyR4|%K1)rogD%0lx-&lUN2g{= z<-3CyODy3Y8Pn{&H1s9M#jc-%wySgtc3M1lmk|Z*%}96W^)$vOgyi+>2?0BSXGE4k z`(`xUyze6&t4j8Mwj5nT>4_jQa6~?amj!N}<491>hK=&t#9VL-Ioq$S(kzp4i3VH(4 zK&4L<>C;=yXpG4OY6wnd3| z*hZPWgDMtD!H@nd9XRH~C6fSyK)CvSA^1eJLAx_sC?Np44IO`PSn%Z!#j(CcXW@|5 z{Y9Uz<#1xRk8a%d}-Vqxl#u3y{q z$l0>NITTSaBvFAj&@j!BpF=Spk<{t_>X3#g4xWrlQdl~(S zJcZ#`2T@_VEkWo{3Ga<*(!K9-9J z3|r5MsHNddN@sIK#&vw653N3Ep`e&wtb2snD9Wjn zZpXN9UyM0=jjgcA+JlORA_R~gZV}ic^;u7isa}rBH=mN#7&uo@2Hjup4-{>FLzSYOSJ-!m?=of4CxFTm^iceDPARi>l~lvFA6TNwvo)%ZG!7j$ zCRmgG;%|%F`TG!Y?(wO8^zdpIkuOsINZiP15)#j65X9WQL4{{A6j~?8s<}X&_-1qH zF9Fd-mdlSo009?N~XVUvH zm&CHZFQ0}|kE6G}qXn=KhRP@)Z1kjNjG;jHHjK(PC;oo(#I}4@Z}7y|lCF)_8ZRFa z3tzd;&wNrymzTep%BdSas;$O%l?~buAGQgAb%#X${Mf)=cxIusR-sY5e8ya;OawDI zA#4hNi!hya_E-C3IEJhwA%+ysBL#Z+2PQFyO4}ku#(pa1V)!=6qhR9@FV);3B6F+b zESwSEX-VWW6Z@UcKMH&j@bR+;mWQQnw4ZW9475WZc}{o6R{GS zzNk@_mspbCj4Volizgg)-}|CH%$Wthg}kgXL=K&YLi*!x(gFVTJG&AG3$Cf9MC!yJ zF{wB{N0t+b?rvXi;V9+$j+MGYIXKT|Yna^>7)_B2=g)KKp zx{7%!i8vb50{_iZA-osBww#<*;e zeK064L*`lK{@mwEz<+baq$D@|uCx5i2!DLwlL=q11-16dL~|L{)QQO}qRl)bvUA37 zp0CC8%A&v@GlbEpke=UQq^cyeu$?Z<(oN;P;UW8btiAZCwMO&K zTv}%$LQIy%skv?5MWy#>WIZgoe3s;#^5i65_3xtMg`^Gy;h}F^Ml};}TqJMg4~Y&w zyDi@E0?xcm3r<%twr^>?m%>OaVG~jUR}^$}Ji3U6i%1B74HM`KL@?M?_%~^Aca)z} zpecU)V14ATD`oVHN%}q}YQs!1wn$j*6wq^O8k=p!vlI1*FH>!n@k(@f^Y5&$E3$5s zCY;5=G9y^1E53L{*c^qHLif~u==jmQX4vUa9Y5p+s?kHFes4KBPW+RSBYk$%(TQ2T zJ4Ux8VBjj5*odwfA$JEsELSmDh!YSv0-sDefl=Yaa^m;12SpWqN7|%icM>6EQKsVx zYwNM)cHA%AbD9(?qU}J)5OnDOAlGFlXP^mdDU&^sOeWR$4?Y@}w6Md4L|O#)$AdW; zK@YX|w)?-6It^FBU<>UXX9ipuyJe=pLkeT}a zu_+2W*&|*(E6hzQP<*RRgcTxy@7OB!oP-9uE;JWI9F{an0qe&pvB~#kxlU|;9XIl` z9qM+=(=8rsj_D06$;QHrI%3AtPB#-^o>U=YLYPkB1|;=Dd5t~!aMgWfGT~FNsiaAO zA$^PXzObtiBRbgPmRjn%m;0QDx$FnEA>0oi=mlW?c`3Edkd|l9q2MH4EyvO&pLbze zBvwapYe7&tF-mW_jAoUCet*VF9>ZElk)_*QLPG)R+h9Nd8QL#oQRDty_Tm= zxt$E(#R#6r;gOo?rCYHS{3Q`F<*Bb!G;NCL8aL5pE?cGy!=ba9;$ew{5dn+N&xPMU z4d?9fJf`MBbrPr&?X!z<%_!v7{rfBH3-gbekya)IjQM72o>#Nv zDKgl;om>T{s@M;UYDBvUg(=YWmRfO+SuAB6&Ft>);VVI<9W(tgt;g6`GrotHk1)Fl zXt##wDqbfd2b2)>eyX&gGKA1w!8qe7VyLqG0X`|*-N8ILPouH>7TWnHpj@6neoGdg zg&(2#P8u*3pg@1iFR6{@8eMj*(deEDLVc&u=uC_qg~2A~rMsMx5Qa0P`);FlqfdLl zqpt>XggpxoPS9Gjbt6SGJpMa?EEceLZ4pWRz<}%WGsF9(n9B&%;0trQP>J`d{mq$( zB3t7;uFysMC*t7QsA2(tU{KPc&xD=OV&|EOcCnZyy5y9qMt5?MT*(BP=sD zA!8F~Ctq{HDxteO+`Fd&cO{JW46a2I=5FRlu!53W85!rDb|PfB?0V%wp(Ji$^) zXMFN~ae(bL)irXVUxnJ;Sjy6i#sfM!QB~K*$x=w-9Inhw|9%?%qFm0XL%6qH@^VxEcFF@4nmXA(DLU4#z$!Ssmv zk6a!2Z22#V4`I}GIp#&{3VHZiTd|HuXPs!PpI_j!b9M)7(FEz0d4Nzad2AqI}5IJjgDi55>R zJ7DOU1u{NQQLr`ssi>B>80U}FD_6aea$TTvr06>!_b!HhR(s29%9vU5%a=;a<;!HM z`$6U&o3-YG9S~3L?<+cNoW-v04AI^Ec(e^H`S5 ze{0gPvA7!s=BjiTn7bTZF*t`Qh#6bfSSh~k<%3wG-DqWPhFp(A^$TmIWFwgj*tOuH zis{wye+!fXM?96QWf+DP3Tw@+0En&v;IGY)t16jR@xI^yiDFN5^d;_UMV zx=GAufz1KFsg&sq{75#j3wN1#xkHj!hem+j^yOzO=c~>}@wN8g-aig{mfnAJ6?n)` zzYE)6L|^gUi)OV+!acBOC(b$SPl+^jY>RpBxcBLe785 z5+=I!dfKJUMNR4A5s&%xqq4dypz^`H9dLZ6ZK5bb1X}ZNkdc5Q3f$mbcSPf;atZb+ zEx@TQyXndXv<%#yzeCg+Fdn+rD{eL0Az0@vq}iKR)H8lhtT&jw+fj<_jzR{UT;u^y z+O7hcdHxu8Wue^gCo8mSmv`5pb_c*5u+uX7$L#e%&|N?HOf>?~%YMi1Z5AU4Gg(E} z20^po8kG=n${BqOO{oSiyF-zXcMcOx9=7ABMK4yq18>3u@Ss9RSp343!nu3mhVWbJswSs1UinS9mpv%*-OR)* zS>(MZ@&Ftf2{M{}9W#oHc`O6gEpl5Q9JTzi#TbC)hv8PKJ=Ec*!Y*z3w!i-nh^caV5LLc#gM#-P)VZetcR);h?aWxvINCEyfV1%tS)R380gv;V-4ne+`#3*WJ6 z&_RXYXTxj;ClAIq3cDF}xh9}(@{59hnl!yJlN0)UGTy6$@}06Dba9t4z`+UKyB!_4 z#p4pdi{Zi?o(Nn73C!g_z}-=+aA6w3Q&F~G)A~I6CHK$Lt(Eo<`0GKBA_g95>E!s= znOt`fZ;Bl2f@?-MDG?(m;2K>#gZH@jm3QSqn-mi1+u_wwFr*m4)+`u)n|poGog_5h zooD-04nCsCqJ7A{CJR$Q{Y|*QbJ>ymIxiHA`>a#+_ZeFZ`y=Ha_wZ29r-)>K_h4~7 zC;`6+)T&!R=|3YOEqz121#+67Izv}@`xKS5HZ4^9@l9R=21B1+Am)IEQ65*QHlMJy zy7b|?du=y~#d$(^sKYl)-RmF29X+rKzeTP#XwldE@i#%xgR4FRHJJNaS#itj%4=&| z61HL(BoWneHA07!eoP(=i|hpXp(BRhVB@)ldu2)7 z7#nv!zc3A>KiGc@YyX7{n|pV@`PQ@tO-aD{2yGi9SYVeMm;5fRDf2k$c4ifO_c4dR zjbWaTa`ukfWB8#}L&b`u;)`5|H@naO-T=;y&A+lNLVKWO&`p89q+W1HnOLo&yPrk5 zOBGAmrN{$J*imK@X_yX3!bwSQ+n(z^haT;8rVa}#Nv6Xy&ZA~a5a$8o&s zE$xtr`+xdn1&Z&lBxphqETkSqo3N`G4C`cnXQr}S3Sh3PV*0OqJ%t*L4%a8IYC0Sj@@gfpr3ltqdj5= z-iBO@tJ~o&&uw*y1WZ}a!(UI7&-CK(^!&t_^}@(ZCR~+brO-CBf);lM4_fsb^vNUK zFVOwm-CbV%6G5c^hK8FXJFd=$j9+5F1B29LdO^>3zR;IQvxDIn9gi8T;Ar(L?lBCn zR+Mfu*INg*p*k>rxc`coYbnlQ|0T345ApD7mQDIg!KM$+v^cQwmE*e?ks!ywemmHp zd0zK#h^>_vb4cslBFm7ySr{mB*|e`5HHI{eK6`gqSdz4v)WIf@B=SX4AeBu~jl?gv zkfSPtOlWjKZTm-V0GVOo&$BZ6g1#cu-U_AEQ3T`dq3FoiVGZkh>=8=DvD0wm6@Mt5BIX5^IYWgAOt3PUhUi)p%#2!Yi8Qd49zWOEyR}O1i z8LS-(+%0T?Voop_>~lKjk()0pvSx90D;y(T%~`Jn!AAn3YJ@*^NWeB$&bIt3b~xvM zZQ&5iu?4$nCe42v^Ousy(PC)uf!V?#p4Icm!!}F2E$l5Xgp#vP0QSvrgdSa|CBE&61y;BN% z69A&-k@nj%rsi6U&z--u#Po9H87 zL*qd_3qXf8xOUE*jk~cE20lp{Zwy01?1}p$a1&9+%C~?BO_~H6u@joiY+{D@#@LO7 zNql3H;WghU&3xRQc&ue^&ov2UknP2h%zVyjUK8k8)AFIv=Na|lN_|4d345p2TeFio z#ptKZb%m^^)nggs*q?tiUo{P;=^i8ZQmja6yqpBZ7pWp0I4K`wX@dTeL`FW1v3+;W zHCJ!QXXeB=#zrrp@c8`2-)TIy*+Tl}`zX*Qx=yL5DyWpy`xd;Y#>$@0BvV4`IO`Fg2u( zEiQuHjQ)JMs)13i8=@3OS>7>rDtG^DqEOk(`EH#)j=bzx8ie*$RzG8&2YVk7?t~?k z5fctYREfUvUMSYgT~XQUb>COyW9ts9*ib%Ug<9eKc@*}&DqG!uN!VM4MKbCX-`mUV zlGB+UvUK`8CvJ@H)eTj>KFHtEJLIW>c_U4EzyE?FJH~mn^J0366*|TtG?;`s0CfE} z2-6+98X$r(F)O7q90}GL(r3OX6py5dk84-n&rKp0MIpH!ru|ZD~-I2u|Mep+_r$E)0RxOTfQ8! zx$W85AAUd+2KRIM*+$P#7M`KhjU9Y4 zA0Op3T86|P3>zsG#j1sq$JMJlM}_7l0|db3uYY) zr)4)^3jJz!!2r2gdDL46UTAn^jr_2Iy&kRU^C_x&=FX5YEhk3hXhmO7*wYyi@+Q-J zzYj??!dBA3-v+r^gWn~YT6oTfM=(@MVrg_Z9BXB|g}OGL?+FsIZ` zP8ae4B>b%26nx5e0n$!_+xMjr)dNhk+~=EyU^%N^iS`Ga;UYLw0S9iGKFLITkw2TA zrn=rUgYPzQ1`<-DHieuw!EYWN?MBX}#Ut-*UY`5=zd8H=xXr8h;kEBW6fNqWgvTIz z*apft2wYxLCkg9-Yp_g3J|9DR?|5P>bs2Xy{o*yCs*j@SRF+bT4Z+X#2be=CjO}L) zq))zM`u9-f(AN+RN+His;VIn&!}JSvut|Z{5)QV)husz#%@}Z5+~6Wiyogy21BjX?y?B{^J*~ zDc(p@PqnGCUL0ro8EHK+_{}H@rw;rm!x#A;Nob=Oyp54sG{9*DW{a#do$akf-lgho zl@w#Wf#uRMx6Iq?VhfOuB?_)mCtE?|j^6f)iC&FNq7cMhsu@k3YbcsW7r3|9gLWAu z%R%sfY1ln@#$P}9-K!r~zXp#W{yqafU_<;`bTmcJ*^3uH|7W@MH4#%w`3L{6QIi zPe?u`)waYQg#pJoG1s%sU?pTC4K%AV#|e}vA6#hFB=aynmBHgq@C=sKi2II}p?%m~ zi*I;qI|3n-(FFXeYw?e?`|FL^ zh_*6O(0z=n2-4X$q9J^kNt04#UGZs48KO}<^zjFh97gf>N!Fy3Vpw%T8z@kg)(gAr zvP2o&!2rMZ%kf`-kImpaj2V{GY^1nUhr@g@jN$}-Ut6>TiyBjRkS6sqWQVRwiQC`k zN&MKcdcQC}XtW_T#mj@AHXuEp`T#uFODweMTHJ=7R!bc%ujFA8)G3vqO~X54i1O(S zx^C2vdyNPmpiPb~JU^*=WiFg6MLtSwed1`B}J?VsSJ9YmIgb3y}_bS#&pls!dCt`-aC zpz60oH-B~t(=e}p(ZnMRB9u5tk?X{=Q=s3R?_uTc;c&-T+UgDd@aw~}2)mnNsgL}U zWg2vd&viu|KBUJt`Sh6t)a-81FY>1hq~z?02Q+4)wY2-=xJjxgG}n3+pzKvFu~VeX z(qnATrs5iXD&=gS1znPpRqMQJ{Nin)I(&3kMqR+4hnh=vN_x`SptZ9nAuM0@D{o(!7ZEEmJ- zQDY$e7zS}-Wo;cJznkGqbHg$15sJ91|nzw^$BH<=cRLeqd`a)i}>rc&ys$4c#5cMNs zVVohIvW??>tC=zI4V~+n7Xu+=;>lm)>X@FD zvfqMF7L5W(`O&6gQw|l#^vgVeGTRhSwqRAqtBprE2|jo0JQAvH6S!x?12gF)zV_rI zL%*u~uST-VARB#O9Nxw_8dwqn9W5LbsP@I<;|hjsC%o8^SMT|+dFx&msLacxFq5`n zn;CiJ)@2JbEN|Z0O)q;c05kUYjyN*4hMIDMxO*jPRyoXR@qDQEA2`XeYuSEP6v;6^ z2`4i+H{2V`$S6TUE~&EP;(=p@X&e2|LCOpz6&xCwr``tO321iprW zIe&9U{5RO_wJS}mU$G}PY-?7=O&2I1KD+E|bo5OAmc(B#;ns+b7hBW?y^_K$8N-V* zHVgf$mr^2FRpViSg4_bKHVgbyMZ2DkOm?!_=QV)-O8KC~6Vv}z*8 z6;#>sQ62Pg8h|=sghn+!{0Rt-bEWc*1b*r)$`grYS3P2Aifl*^-~S%n@MdqUrg@}d z`Gvg|ko(WTylsLFjNs>Kii{jzwSJAYHm^;5Io{AmT;w(amJm zazTeHvDzNDBHlYIJ4LEv>TOG+ZJi}G`qN_CvHZ{)lKkg8l0;bgmajNWl2KEr-vYw6fLO7902;V&M5J1Q0A@ZAif-a#X%7vc(mbO6$^wo za<+eI@5SLUy9(Be7XrDs`|fo}z5^WvG%b-Yb)T;~ zSyH67S7bUe$C(R$0bkN^4y+w%)wx_LihWp0D3?=}P`)K;iBdGZ1RrYM(VWf8mptd~ zdqYI)no0LT37}lz*P(_L^9+?#*_S8ZhrTIvgUC4*QUtccd=1K3orLNk$6Fx`;Yf(S zY8yzOhA5m$tPKfJ z+&j+e`4Vj{&bzNodGFcw=L|VX%Oz6+{MTP99>eFf$arXy5{zH zq&*skR`V(kemDJ)TX#@yYKIZvprmRmi7wfdP}BiS)%p1r)$1l7_e7`48$rincmADZ z03ecSOl_KjowDhCB8$P>*smERQP?{o{w+|KKT;rr1p%{a$sE!j7v$Jm((#}SVjt@; zwzIs9PJe(4?nXNk9=u*)(i}Ta(tmF9u1Dstc%tMYj*?50cA*yVPDgT1O4!){x9cgN@&9lu zftWKr8wKb&18Q%IcIr*|>9>W20(y8#MoL%@)1YX>Yz=GA)h)ND96?t>ee3C+s~fM`-M93y z#>%~6|G@rLccasnPL(*F{A}@oy4}Z9#)`i5leuC8(##F$TsA=M9x=r#cgz}}eD4wb z`^@>w`EglhEBEi!4X=c=VAC979Ml!?(t6DBZy^c)&%w(=xvu<}U0^wU{ccXl zFQM@)ds+1pMD=*lLc<=qb&?`^qWkFh19$4rHr>Ytw+2{gYBbJekqv_KdJ#9cjQzhf zPTKa=y!|sv_U1GwaVLu=$3jZq!sn#6(i9dWG=L@D`-{8V$UqWKQ2CeB`k-xet zRS!{m{rZR}FJ2Y+0sn3-aGj4AbFUU-GhG6O2*s{cEuZQlL^)}T<=jX>*QBL#8UMxH=HX}VN@)}a%oU@k@8(>7r7_^ zW-8c+p{<&}V83Ui$D{Mtd&=6rkLD2F1rf;hcr+GJ^R|bMF>$~Rs1Eck_YMtc%Y25C zW8$YW6#pzTHG`9?r>T`CQCF2p-nrs!99uYZzN2;7Rl({2el1-796DX1)=SkK6VGzG z7lD}aI+0n)zE|uHq7gtK+NRVHpNIv>CL6CyiG@Fh^OoAtsWCLw)#hPc_*$__GM)%(R?20slU~a?43(>OUp)PXyOdx zFquhc&0hCN_&xecwJK5Rc~`hyUI6@)0`Kndi^|(5%pf>JggOAb9k?8$_cGyOl&xfA z^acof-&krT$iLoTNMj)oL1Sj)`$dR2pcJLZZu;w8sXrAhZT3-&S*22y^4llgdCLnr z3FY6U6}Fs;zTDe)q4*Ls&)dMl`tlbUuL+%G*CS>AjtuHeC_P`2!s{dJlDzoE8f$8S zd&voXou~&qW%}QN-)Cou@E2yHD3(vzH{M7sr~eV$NfMjGNwM3?^MsSom)M7=`5oRO z_uwk_Q|7krQ{n&MPmC>LwI~%WO0fzc>FQ3lIhL}WHf#|vOa~3=Ty*C|o-;XgEm-LB z(m!mb_crn_K{gbjCr1C6Cu}(q@KJ?I{SZsBh6Xv4uJB1 zyb|gE2Q@&-zgrP>@ti!T=>)$17vF{QUF44^V)kNMPY5D<->-O|=`PklDP9Ez&2dzJ z7S1QQ@(`&G&o4Q$sTMw z%)slLd^dA{5Nktv<)lo))^4RunIoPnP7~cUmc;*>=oxn=vUn41hCd|!4{4iW%J@AU zKt%h6WDV(0*6(|1>*_%AQ*j~O^VIX1kk|JU2|M@6FFip13r7?iqH}r=I;>hu%Ce5~ zJk_P^TlvDHq!n8{B|g`fj0qVeRkqSRK^fF^d50~glh(781Zh7hBOcSv*_3vxRCi=o zC2O%Mss4lH6D=iWpk@fcmE_}?BhPR*>?8H=Zk~gW$9dfEMe2Tu(d?l!ZJL&=uW75E z@v?e|AMx$n4c`e%aW0d;^Uh7Z@78$ytGl+FvJS(=MtUiR$Dc#RSkW06AuTYv7701Z zo3@ZoAgwR$#`==;pV(SqqsR`BFWlv_bon5{msD5AXoWInW16PFjqg@C_s6adtan)n_|7I1I>#%VA);wSjkE1gugR< z+}ev8R!^xj=Y-+akZDwRdTzL%u*7id*w2O+{r4Ik+O#%&Y{tTRUvc4;Q1v$`|Yx(20*=NWrHobhIkjM`SfNdIirfF;UNpI03m#o+TQV-i+V zR^<&7f`gb*>oVm!&q_D^!PrXPOiTGlscQlK+5}Lz5PC#iqigJIQmvm8e(DCv#^NDp zUg`92Auc=Qbv$N5D=Qb=kILh&`DfVlqnc+a&u=33Owe#KMpU1aJtS)TI6`lS;Q!x4 zeCo-|@w_H}nM*KCT8Ptj3qr?g);4+!?hiiUcEAl=<$e9?yWn>FFL~G35n!*3?p-VO za{1!Bg{c1dmx$h~Kc@#_>8`WT;v=m(cEsb)`8aDfZ};0zIQ{vE7Dr7niO!`(u(f8+ zWi*~OoK~tQTpRdek^G33S3gkO>lQ}kuc)x{6?N}b#MJRD)lVyH8C{izZu@C|w3Nnw zyrSVZFB;r$NYn36VAt6Q^TCNU9=ZB&$kx>hp0>S1hnw|SG{nUAaC&7X&< zVDL9Qb8c<;=C>D<)S1-ZRh1?&bFiNmLDR*tIIWV`%IO=n0f%v0Ca%ke@9}tjP&%_4 zA!||yeq=)EPeruO`bO&?nv-H&NLHh0`TgqaS#u*nnfd62CrD0BqC-d%>6dj#UV2ZO z{TC88*^#B$d{W)+WYii+qRTtls6UK7Glw8|D-3q;2yOOR{oFgf-z z`X>m>{1urQF-$mtptr5HYDs+aqa;vx>#vg`73d(&~8@~T!=1lH|OqDcvc zO@lBD`9kda*Ydw55pm!H;gxOVQ?E_9%Pf*MWvOnHKX;i29a=P^`|;1TJ$aG1a?Onj z8j`j&ku2#=3Fp1ZEN~|F^go(UbW#T4J5qM4&+#8iiq6vXDl5Q7-(TJ%Y48WJlNQv% zb^!J@>rj2#BEx6f{|s+D8&GM|Dyp40Lxt?_hObrqsj#J-YRB7C`M_7IHFKfDV)du; z^)yRrf@NNTJj6rr5BY}QGR-^Q>?h^nJW>y9Ruj8g45K zl3&-Mz3;!ImA4~xz)wVfy-wVfvvjB`2Gt93Ll(cHyYn0J4unydT1fIV<$tGNA!*%u z!aK(kzFoC^`|0Y{uMv^85k|Zq=9E844Yrc7-Hy~Vjmd42Af0Cp#qYwzo~X*mTQT%= zE@6P4|6wot(&zWT8PeIC{?%5}w@x31oHwFi(Z6(^e3kr_n`vvIY{Q&hXc9(L#(i3w z%GdYJmgqQdX)g2dZTSNC3UN5LX^H2R-uShPAw+#%%!`l272hLywKeJOrAyA1w!d&b z`4RQW{QMtnn{6Wb&>iBNP9vhE0VIig+cdUIrtWn z!T#7Dy^VDb8{BitX>qE6mgk^a=z|Np-iCkNdTODEFmyDd2W`4$(GKe5?p ziN&osSPw8$9@-SUNB>}vI~}99Td+Imh1HH9upD{;udA=I4iAu*S@lVV9j!$Xy+ed-f+H-uBg`!dECEvPekAGLI^ z*R|h4)7tVMwL2%DYp zu8oaDLu%9-O5Ku*s^`Ssxv&40&=A{adhg=bz^Wr-aZ|@?dt~O?U_KgNFYAL)1&}7(t8nm1!nMVbj55Fy0@HQSGD zcq@)uU(&}Fk{S1v=nfy{yHLL|`mSb)!|C+Ng!V5@>0YBTMbQnEJ^qJwJ;dodHkrl#?dsn*M(;9Zc=mnHR^Yar$$z&@;NiGI6W2fx=(5FM?abt9;5PK2QbQO zM5C_axRw4$&7m(atMLO(#072c^A4M?Q)pBgsVvMpjLH(RyekjzCP(?NZLvCdToJsaec8~^GWSwPg2EiMEAcz%1!xZlXO138A{}~o8nsL;lDurK)M!OgPmp&y$N&rLdqJKc!;-~X*y5!vN#!sZ^&bnM9j9i>UVT6*#1k$@9U7+*-iSN=7)o4=nTC^rv|e% z*Skc|s-fg2*Q3i2KZ?U{)8(kRI&F^7Ek|DZymFE~^XQPVi3AVvz1)V;)+3NMBhAIm z-B0VZ5t{pGK3^eBebib)%nEQ>IT^hMYsN-wW5RDo81wsjMjwh~VD|*ZRA@xmi*N??K1;!= z4ivmBBz^N}Iyc-;x1N=xnOV}lWQ2IH(R2!1O8P#{$6V^rS@V=U6MuP~6X~i9Om4tS z+E4jJdL^-mHi;Rxx()@=^2a_ZCA*ojVqJ8%?$-S+oj2(ECJbBuQeH5Js5E)>t~nBW zUHbX7s<>7gNZg05#F(i!J~M!%{yhlJ2p0#W7P%{>&BTr(-7=Q$ourkX(A;vy0Vb3_ zr8M_86O4bP)N>V+oCYxN*SCz3pJL?KN%Z~DlfpqI^l03Jj`?D;mG7go?>}Up&|Mn$ zjqJ#M@|X@MdG}dbSC&sXy`@+$;;G)y>`rxU>&xdb9Eca!I0BdKU_uPza5oRaYpXKC zql|H!ycw%U{V_ZCFP3w9t2cF^h0lB}gU4X&cS7Ec21I65!}p4`!qpSRdtIq>Z7tqK z^RQpKPFX<*nw&mF?eA7&6xEc<_lpg$8>CUe?6l!k_yxmT$2NwK?t7^<=9Av-dssZ( zh}qq>U7?0Qesa%`s#wnQeUVz!s9W>Zd3G-%0 zsT%z&7PVUHT)IZR;}faBsV(&<@77u6L_P5ujLVDgOs|ed!vXkeejhURA)Pwtd9LvX zIf0Mqv$-Q9m$p~G^ampk#VF^V$k5zb4B!7vnI$X6TD@S*gd)b2u4Kg1ILhvLB zA)#x;c&)mup%|A+?3R@^%{wrg@o5+N2i41S>f^i#n zFm8M<6D!9u@$eQV^|(&ih((Mak<8es7xZr*KneZ~ZdZ+hn$PLH;b(G8RqxGHpRoHL zIom!Eubg=5z=;GNSC;5~lyaZabb?PJ77IgiCpp12>gz&WWuG?_w#%Ufc@ zvPm=kMCzyeWOZGl_uwd<7M-Qz6a$%Eqe=8%{r|kY#A$z#u|--&-=B$DsF{y?{G3P0 zWL+yK+s%h=ot80Z|9pmdD(m*!QF^C;q}#y(6kXm))|^wM6b9?hKP0WNp?a5tq+OB@ z6+N9!54($@^&7(nR%F ziMM0upQgFWn^TM)v6o>5mni<}G-LMvM=$BBUB(*IH+hQu{8uPgFQ%jQRMHo1pP_Ykxg?JXR9;fPt3|iTIr@lpN9K@yuqi2xl|{^TBVv7X0=_7t z^m~d}R~rd_>4ST|HC~U>3GT9y5a}(R%k}&zC(x$uef5OORR7+aHXVn<)$TeI69{i& zMM(Dr7`9!+Gi#Up_nK`Te@)N_BU~+Y-Dk9vSJB%{+4O7Q;o8gw&nj(bxoadIrU^J*7=%mY65M`ztQm}SzsxjRB|pdJ zbvyh99mT0Y^PA^h*iFgA@>WG!SoXnai?VZ9ffcuiydnov4>W%J#fzT z$1TYflM)-{v!CMJu{-X@mgvA6jsQ#DC*9;R*r9$ht6241 zKRQhGAobfRvO;E(TBDGTf2WeSRe99(_H^!-LrzJ3I^Ees#&_1rWvkYEl1SlS0raZ8 zfbPm%6t+A-{{j2yUHv+pmpsu-yca_=7c%TrJf;1cQ*yj5W!=;pXBpG)`8h_kS7t*! za_?hX=~JVtys+x^M$95*aJ2l=CpAyKh1VVN33@cax&Iql^;143R6c#53_|Np#!Y@{ z-(_T9+;21nj+=5Nmj%_M0ULAZsos(WV6R)M)lwp0KOq7`8 zbxbkT>_+N-Wy$uurF|<0dPm-+`|6<#+-=XWonm2KJ4v5?Vj8}T)*XAE!PA1s4{u4Q zRUgS;ET&wSzld-DD_vfWBgNtb345-RVQ)n2<)26&Fq-5hkMYhIH|xR^JT~gNja04q zx)hIF+wm;Eg;h{<-6@BuUnhZ9mk-m*Z6D45{u8@)zhKo@Iq<86>a)Q8LMguX4e(y# zr_7o(&kecwe(8$)Go8z~U*cxn4R?cl^zZf%R6c~z8On{V8%1mDABdTuykvMLv2`oa z)>{3)a@I-8-esIStt^%WSv6~rQC}RID#~T3KBTO*&nR5xbOmX^UN_uHux4`MYH|4I13R?r;$1`JZugyiT)> z(bNgpPu-uz9GaMlncHdc^n$VJQUS{^;n+I;hUszxj%zpI6Vr{>t>zQ|!hytA;#bYp zoPNV|V&$I+pDE92-XVMsE~52?o$}|&b7WLmXZBqRynV@hI7fbJSF)Rw)64mWo|`8O zUA~Bd0mtcaJc$7>*E9UTc8m-iN6GBj^juSmG7C>i!ml!NmLhpQ zy8qCFPQRaLK)yV5-^VE@@PPQge$>ydFY(6{NEkI8-(BC6YCj*t+k1q!*G%r$s?r5Y z)U)aPORhj%zs=+v4%FSXn9#*9Ng5^Qxce{irf7CGQMu1$;-**6AWB)GFuzZPR{R-^ zeMlMAlD2~)32S*)Jbg0~PmCiW-~`c&B5A!VlgQ8Y=yc;2*+UxBHhwByyL=Lpr-N!> z7|?bGeUiq=1F?sK={cnN7$|T~qvOmJ;^!_T!)!VE6aC5745;YXM&((QzrWa>f>+UW zJG-2$5865L%2&EvoOH6W7Pa zG>f{lNIh{!`E45G*-x33em~-7Twk@6}NS`N0;H}C_n$JzX)v*5Bh9TcBzL&Szf_$>Yx?|{?r*CG7hAEwRbv*dnvgU%yHlDB?@o;S^6?6xq%^dN&)JY!5+ zb-L$lpwQWiu9Lo}=icG8FMdVN2<;=NrEEx)JjUN{;lJkuJ|ly0{d^JQ1NX6V5F_Kg z{NVN$c#W42A*2R2RaNuPSc3IBy>n+z=3)SLdX3V-ERgL;?X z{Ph|>I&=KLPE=m33`3QDIBe-j@UKsB)VUk?!!*3BE~L5dVw%lchS|cZ)YxGsE`O2X z^}sA@51dMEzY5ez`+%w6Q5rl+c zwW&f#8!>}yC*jhlAMSsIYOb*i`{c1W?^{UYkSuPqk;ig-T$7N zE2?1}wv(nm4WrJ)%9!jwfo05OnzTKN%lBEdkY?_EE&$IAKZ5th5Z-7DxzkGN>i-iR zZXF@lFpK`p0_j~}%#z0k=+&ni1w&KmBM(;3Grn~CQ%vQNv*=ty^UR^*RF~V(=9Twx8b8NTv1gPwfP z(3R)u{#9{) zcDauIR&&jM&#PuD*4*(&_+6eC$;ydQhosRbk$Sbd=F8IJ4WmgPwUFcyjp$IfIf-qy z(0lw&X)xv5H87TpX0ap>C@0mgjLb(7q#ErcEkiSl4qG+LQs(rwX6Q?15aIZ8!>pY!te$V;g z=J(7PF^TW5Z)f@g_%3EUKi+M}4};&b!22=tCK|J7eP!mZ@n+_ti_ASUk{^TnFuO?_ zvoc(na@>}2JLYO0Vk=JWK}L;OLZ4Qetp`q~=ceB%{{0KxyO_x%uMGCF46;k_lX_e; zh*>pBmlT?DznG+ZVth%!7HqYcZI369;+JZCuM(sJ0N!RDLj<> zbYGf9$RA$v`E(~KP&|Y~;usa)5p%gC1&bdtv|lV^tLuAwQH?%J#JcO)mOj0wQrz-S z1{cY{Ur9aMiYgTSno3DXCOv~U(Yd4#J+^gbc&xnr*5kC>rd)XoXiB>b z*7B&mT!YV_Q}z81$GB9ikw?{OGO+`-Hr}GzvK$)uMNqd~tn-MDG&YvTzV>HK!;P`A zY>9Q1c{KVOM3YOFG|pK>!?yjXm0Z>EanUInUHXXC!N0|UJtfc2FSI)F9>4feVn@2; z`_DPpr2Q|Wc97vK-#|Uhubh4&Iry-=>e`voJ&NexjFEow-8`XAk!|-&fg743(|cwZhe0N zzXKPsp6^J@gX3tm@(9kps)gnqA=Fi$h5vo!+P)H;X`;+=4Ke&Z!~*<;*`dA|DH~zF z;3SQ2L}3y847cmbw|cyw+0XTH=_qgIrWE;IBIUU_qq$84Z9N|ncH$YakKBlQoCKBB zx6Xc~?{gplyXAYIc%AT_@_$9%AvmBV26szB-fpG!_78;aF(xeiX9B-SUookP-Hg$A zIf+R>RE*+lp@emuPV|PGgbgwwdfHFKuTbw+OnY z$Du$zp}arvJl>qhSKo_id|Y~>W>+1Rl{j*hPIa!66Vj9($Bf9;ylSM|b^3juDevYG z@^@^Z7xx+1Q2Qm$7Sm_F2Zi01XdbKj#3z~D&(S0;&myVuW91Cw&1e@&TxEl_(YpAL zbReLe9lYvK?MJfd3ijs6ZXJ8!)w zC8Zi|yWAuC_;#X}k0oI8D?+4jzo|?`!nJ)MduST)0HsO2ZZdrrHt>N2na@N*S$m!8iqkNaop<*v|b=? z@6iN8o{1}6Lvwc@8w?4Ll)Kmdhl>+&%J;|7{{}6ajmG!dChV+Ya8KM%aGxOre)Yg5 z&l1n?N^u+fFCM#pA!xH|0JqulR*8MP4?04^Q?W0qe@D+3HOZTnO}CTcDy&^j z-{@nE_npY-<1Oj@eF**ke9s{F-za&qlc6!T3>>1&*X{^Pn~k9Fnmh)cc&Z$v72Qrx zp{t&^^VaZBrkUYk>L$ajI_s(W%89yFhhn_p z17?@k;PFd+f|spV2B9tii((1g9ZGVH3Az74j=pATN<{|jv}4#S@zG8zH}L5Zqay3j z_oX?bG=J}3djVYr{X+NZeJGsUgibqL#n&B2pYKxX|4=z*lfiU9D<4$oQ{7SH$+!4I zYvqIu>%EALm?0nhN@DYr)4nM!Lvuykq?T2Q ze$|E00x=>gimSX@EFaI|1m)c(uxl8RjUHn#DI?1IC-Fo-6U-51?dlM*HDBzN{v?^G zW|MvwBP)V`X*+y2Eg;xB9>10+rK1#U@MezpqTTpSUZ^u@D6S)=$&C69@BR5WnS1NZ z9F4=uW!SA)jqQOy@V8$oZ_Pq6rMnQhdNXZ0>#VM<8DwfLQhj4decX+ZpmVb(*sh78rOje&OY31icZ+;0XH?^tQ0-xDnmmbb*6Obn85ufI{VXhRji&AETbG8ZJXWi>*+I?+yQ>fRh2GIFt8H=5eoIcHL_Ki0Q@p~+HdrK1nv?6{Igqp)@A3JukoYsf z@_HEXUa#}%r{%P)YJrc9d`MrH)8f%Wyu4}O2CbaQ=!56!`RpkD$DEjAiks(35_j@-*5!K&)wwL6@#)_GQRP>akA81FtZ#7 z`B>dgeO6|^6#<9l;_fW3!ov;tP19$zeJ@&rc?B&51La*O*Z^Z%TXyQ)8|q`$?= zGqq_hA=^E1-X#V=uXA|+9xFXPjo=@|hu-!+b5q*%L&kJlMDee080OGP zKbt5Dh8(5GfsNXg*pt3P>rhxDj4sM<6kJnYJVT6??B?|NQTt^6*h1P?>1<`Qh}Hcb zTI)BGhFKFf_X`;nmBoLwneN>Nkn?bmve3!oZn0Lf>666{QU<^Z6_GOUtx4sYRFsQ|0G2qY8@67KEafqd;)_CpqP!4tdpV}Q5 zK*!_SzqCL((=2Dw>KoJj(w_`U8O9j(!Xxj@*F2yB;~jGuS66win88e0HHh)H_DrZ* zjVb;1Q<_tg(z5xCPjF>?{XI-tJXbxcMvX|J*{(ikieGYj4!W z3FNx}P1d5j%4IE}^9M19245l1rl~Z@ab(mqA@apPfBHa&znRq&*cGnkoMv4|e8hajfM@v|j~2Wh7=E(teHAgw(bn zwe%LrLCUb*X+(_l>0B%Q+#CEwZi0G;PIhEFULa$TvNd-HOEZ~>?+vkYCZ*uBYoUCn z$%GB>MQ|^BF=qakGnWQ4*qoG?j~QZ+r(}D7dUi8ou#ffx88m<2n8NU~L5wx|jj?6g znPk+EVZC;0{-wO0^s~YT)yXHGuF4e`pV`mQItv(BWKVgIUQEhr&$LZ{vp9JKv!-5R z)vYKNO^#r0uMifL?qS|PQ&~83HM4fUU`jztrgW*!xJ^Eci_B*H6Gz6GY7X9b1Y>J% zXY9k_j2##!y}+Hm!O|V4_aLvGB}q9gi70PK!0mKW)5Lz=GgzM=1bVyRAF+?%$|nfe zEUvEEGJLYMNA9J(|3>PAK8mB4a7?__5oA2vPW#6V$q$=Eo?l109@2f*+@5Yh+Gq2% zy|j{CdSo`!?6`oQbrx#3(^0z29wA@OG_vnYd$CTZeQ9SpJp7M*heLE`E!lY|=)`EU zGFy>ZLz?gIlcc&0l5YN$_C0KA^EN>of>rp8)y(JZILyz9lV}-(=a%u}js1&rk@VYe zXZ(+7zjSC8tqMlrTd5m`F?0_UB=zl%JlCzU_?PWvy*D$;yie?p}X zg_E_%Ke6p6?_cBoxJuLYzw!Z(A(e5z zI33S>iP|GHioo1WxF+1k+t>^5fi(%z9$N3=`^0UVBCfnLl?!H&l&(GFPn|F{eMVg8 zM@0OmzUPlAV!_%HGp2#&dbbGc86=L6yub70BS$&5_Z{T3YeHDuHj)Ef=-@hsHghY` z_KpwvzS4|M`bmfOWpLtc3hTNtAl`%A^V8_Dw-!C`|4x@NVjNn&ql;NN9S^5dRA2j* zZ(EU)b4vbYUuivib)SpbFutL9nfAoL5lhSTXWE>*Ncvmx0-nDlKGTIXX-RoTYjo%6 zGu(QKoM+#}FPJ5t^CYr9YqoSw{Qpl2iC2CnvX}C};woo9xJf&sPIQ?op6i8k{=4Hrhrn6}FFw$BLC;jO&x@2`HYoNGOr}N1ep|kPB3GJ=Ytm@frf(Fc$5AA?5 zb-L5OXin%a4pOf6LiuTz#wyMH=8Yl$ofFyP%}IFHoLI-h#0T1An0$kf(~k*x{8jIn zScl(VQtkP#eC(6(tXod#=I6Ly6YsKp22M?laDU{A=cg}t4(NdIbp!6VoyEYsj{8X$ zLV{!T4sRhqpNB#F`NB14PjFPfw8e#tk&%R7ED)=7m+IU~VVej4Gk#KSqco1`=U zk(K^R?~^@+X4lEvaGzZ3(R2us9&eRGLWuTBZ(c~4`FuP&w8mGz*Ym_E+#3~Z7mO($ znFBRnwALOr?c<%5FK)Z?N4l5nuihi&@@+eIw&yWI%!ABJD4=6JQ$$A8n@`{e>c+LaJ%9!+>l-E$in5x46oQH4n)*V38b zs(K@FEwP1lh;3d>hIGrQF9S(-?J58KUxYmvMog=VqY>ZM~>s@6XeMYzo!e8qGM*_bG={4IQcEG6Pi5%CKm2rm-z%DEz~ z#TH4ZzFGUK<;kfqj_@7z@b4_&ZjBu6qmrMj3?iyUV3)fFi}g1!HA=+%$VBOdZD^6y z9M7&@@Sii5klaW3H|dRQPger{>Zo?Ur0nlqqRrZpQPPgK`?cGsT@+d0{ii!wyZ$~5 zqWzUkOZTM zYL3+!`?>`fnKTwJJD8?(#kmaVL6g;+X!fuo_G_17v#teAs=cLIq8Merl+$Q=43^tt zrFT}M*2(~?8lR(P+#qVr8AR1ey$uhF9~sX79c}n$?`*@P6AulSYla){_6{`MUyx%s z@#2-?wS^Zo=1(zvw0>=P`&xOfKb}za@gti2Y){44?`Ymneqv=4EEc)r89f8*5S2EO zUNrA7FH-M;(#I?7{ENgr=tn{Z-o)>nr}iqWCTzF-9*gCjswW*yK1Jt>uL%g7f#C;n zH;!2&DqA#$HJX{ANJF^&MO8q(Z0@{ zrpn`l(D7IRxyLt&?-NHrVm~_8P9(SE6_WePAH7z)hhAtlbWn_}CGx)Z3#X{{4u*ti zUY=1$zM^6ZfBQ!M{xmTTHH*@@-22H5I_1l2WPOs34mu}G%#^!zCGX1!I`;aDu9@4( zpAB7)&Z2G4y`F$SC!1GUN~Wk%$^3+P7ARh`Kmxkvjm zpOV)_y8M#@@@uBhr;0WCN$2DVJx}qKdJOKeh*55l472=&5ue2xx_6MWk@YAm9Ye{^ zUdIDWzDreqZe?PJ#Y|{B zp9z|8mQL?W`IJSJ^hu}emJ5T5>o7hpjlRvc(apIgZB6V*kDaUfPRtjN2Sm*X5X0rL zX4gju{N;NBDzwFCp9?Khv{P$}bi^qOa9H{YPvgcoWo*PF{19Qq$_WG|;Q5bsaBNg2 zOf!0S_5aRQZqVZGPFjvu8FhLK#{X=?X4W`t`kGTA@vU}j9>vmPG0odvRNnDV9IVB= z*d+}~^+Bs6ZWvbApheY5+8dpR?^`n*?@hqnSNp_&Y%4ZoHN2l?>)HEBTBwJ-P8YNb z?~O8n+P^(#JxST>L+&V-WuY3mrS?y?e@pAfif~C=2+ncCr$SS`zsh+W?m&xH()8w; zVQSx3JKue%v@Msqw$3zer+vHbfo`!cc$2oiwQ;$Dz`6M=gdUmU!-s1gSk)N#} z4z<4FqI{V9;Jpe0;>jgR|rl)!wb-%HntAGkjqOdX65gcjzI7ku~TZr`=2kTaf&>yfvBE zNSmQu+4YZ!Nt%v9J~HQ5n*FyMh^u_u-c8lF?VGP%r*|=q??5N zU?~kXQw)i_88`>-!G6#0xWAL%ZCVFB#^}u5r#wVMZ|TF9xW)cS%bDM>w6w$H<8qt_ zx5VaPAkN&z+d@n_)eWJ=+9CJlGOjmW@Vc}Bm*!P*91*U$$Zov6YvJ|j5Uqz@C+y>& zVvc?xveP=El601}Z$xYh8?}o{cWYRWly5wo8Y;jF@Eza%gcFchYU&khu%$x(a(67d=J`-+BbsYw?X9ma*O=sr?qSQ5Sb5hh}KL$%DWlKS7SBT z(rj?OSUDfPiOg9cUW^#_EyQI|#xQb!1ENb*yEi{19*}yw;_<4%X41A}YuY~AOPgQB zsB>2y{mfe;wXvxbbMO z;daqi!>fQ>RQ`CzaQ_!OYE(?3_T)6G_->)m-0m3fX+(pQeW*RugeKwwH3|F|iYMoitX@|8W$5CbjBivJBnDzI;Q$1+mLUi)VBrPaL;ev@^{k<(oJLZ$=Z@_Ir|& z#N+BZhv<2l6Mv2*DE$-u$L``I9+>A(f8*v^ERJ~qA&ssP9<5(Hs!lF<#$%H1pw$nt z-Fgz+#28E-{y^PJYBz*sOf9JRxI$rv0X+B(3?E&R$~QPR!M=^>{kJ*SWLU zg!Cy*XnTJuiqi7t@4qJv=?$ehCU~y6K|7dqwI0x zZMZ7)6)}_v%H48zg6nIhXJs5fiGB{GE8qOXBGk zU`g)lE9Atyrt9XXbY7ZH$8~|^f0n<{MLTZoF4O&_b`h?qPQTOu1}qDZ7kNIt=GCKM zqosVj9T~iMFQwU!wd?k!b`XwdRB&DTOgcrM+g&Lh-JRZdt}tL=1$sIcle=mO9lOe_ z(sv~3+wYUOK)mhMKJv$pBLB)dF(`hZ>t*F`2NsZSQHi8k*NE4QBeKCoqAvK*p{wrA z82Q#FQKt>!LD-t2=SB^}sbP4ZrVPvCk1pR3)v6 zav^qQcW5+InI4man2oi;GIg=M!SZlC+=%&`_cT8_7V~b}i)T2Bk@ip2ZtqFs1%=c- zGMZYoJ7ZGRO8Xw?;51ab8QM42PTvifuDFEbymal*%fhFU_{@8ytt8GQ?RXW!OGgr# z&w_b&~Q_%cy=q zJ5g<3Qm1}P!{fRMhWGm>7;cU6GMv8`V%Rx+fMI!dxna(yH-`O_FBtZ2FEp%PJ<@Qx zNrB;NM6uytqa%io(>EJl%`-Q=+SQio>eFj}zJSqw?W%kwUAKJ-jVc$ZAF;qBDqTC? zMqs|z0=J2e@$er{%S~4JNt<){qX3&#LueKFQrd1gwkIoNJJSWbI>Yh$#hn&^T$5*} z7EQmEV|jft&3?R%QJ!|Cq-%HDgQr+6=!nH%qcN@CTXjnvG2zuGoBTwRV_k69ywNJu z0t@*g{Z{lM=(1Q@CYozaFeAt`TdV|gf_o^R?XA7RUw^@H%YvA!2%^Nuit@WfyOP@4 z2YiK;8Fh&tw1+gbqvT0Hl81gJubq#rgLl}@|EgJg=?X$t_|ZCF`&Z64BzC(p11)UuIPsVWAa7GZJ{C~7i&U)yOoUc zk7C%P4Eopcq{#Lb!&WtBkfmlk^ZsOD`WHI>l|z2X98wNFrjzLw`G1t9lztPQZi#QE z<{5RaYF1ty=j+{ZnfwgjhwA4e<*1c06x4N{r0~^Yjy^ECFK&T(|PStqcFHau8eF-A{rtu2%~`D3f_RqE{Le?n}7?5g5ZX+goB6+ z?h!fxiq@!I=91=$N#&a3I%!fOnK~(HrZp}&hfgmTUwF6w|M@-Voco+3;oJ+&M2DFA zDwDCdjPN+_hW}pqqun0h)#As9(rj@jJ+_eWqA%$#oy2uc zC%;Fg`d3zyHQs>|vu7058c~>fk)o5=C^id_w(Z2SzJE~BOZ{q=fvkBQNqJK+6&@WZ z|Lrklb+=jdpLtYX31a<6>N}n_g%4joWUaCq>swk`eR?aaoPVXlIgQd!J=J%mZ0?m@ zis#mmce7DB^ET4ON1g%j6A$bjrmVo3WB4ys57^ls7$0{}v*CGlFZ_g^G;w?R^F8b3 z5os{sYpZ=!gA>7Tj^O{9c9coVJ${#_y<;7L--^M|?Xa%HOSRt{jq~Y!81BnI;`~6I zipyeakHoUE0J9b9)(HuwYmj_OhZkY!&=bql0xa$yqw6{+QLY(`vc$S%HmJMLwhKqY(p4- zThOM_QCSbTddm|V<@N2Xeda%8@GxmrG3q!R=*k4$r-PKm30l=&+^_?T zy7!7;6Y)N_ttHO-8qr_xBd%edx^%_WFyA5O%Te`GiA~wMLUWxM-@fXZUsp{`RTJq! zV@OMVre4wgWIvnDp>~FJmhU;)jrJZb6RM zE%F!sliao^4YvKkj=OYUd}}Pxvj}&KyXeT+c>=eV-$LO5rh4a-pxbH5(_velHdUwX>l=`q{N`rE~rGA+A<$K9y z-WK(Fo#{^T%3*vUUADMqBTKSJkRR(y+GBkl#DKGiP?ZyxN5ro;J~c<`@RF^u)385`(1M z4DNRa*D`S@%bIa{?S{|jDEzv~*SYB+V+>CiwMOUYu5THeHI=}Zdk9=C|3&pnaS(?w zwZ=|-I^`^)gP4By2I+}n?Q~hm`2Iy=!;X;tTsoXr4^r1X zB{#^PEP3IRquY^k#h8@Izbg+nOdhT}@jfHf`_vxyMV|6m>KfZ;3&Z74FxWj1f9oGO zJo!^u*26em*nrI&v9Xg~aO$O9;MiUmzFSTI6}I&I))PY~Ze?GIk6UPgrCe*f26gtn2)9FN*ny&H>BP5XHxs2 z#2)&)vI-qZv=_@I=7QKpKNG#-CuWF4p4?*2#7nQ3X_O+LixC`sDlWxPLg#xBY1_ci zvDfg7`x37u-{avc=HMT#49nSukJvpPGpZO>w@>$8^#@tp)UHnak|&j7#>(Ti=@k(d z{vy6xve@pL#fyBz{&3J9eh(SOex%O+jd`ZKNNv8L*Sj6b#me>XkjGJisQy4`Y&(8E_4LeYhvip zaXMDEg;;g&Lzg`c*gAReZfSvXK0A4L`+017`JL}_zT2_r zz!DXyK$$jdQQ^b8Xr;L%a_pKDoeSKM@b&w&)gl_jaPz52j#5)L+7^P9+(^NUvhP zI6GT1L-~Rqr*P%kHy{4#vFwN5smZxB_<%K1k^=^vUa-K>uabfL%L08M3te=X&Jl%X zdQ~S#JlAs~BAUeY5UzDwKQj^l7SW9c6`nyqPX_KU3XYk@>2ph7;yXhcDB&0C0s}-8 z&NhgyptnR^JJabQt|kXbGl%812gARMexy}O27s9feNV1O%h98ZCW zXsi2D;N7Y_=cXUhq}1g-&$>@xWG}9D+fY4-y@kN3nqH7tM70JPi|8FPaKA3NhqAau zcgjoRon^Ab_9#B*++4(nxq)%3mKCg}GW zrcLvH`STxTKR30!$ZPqip%(vRS^SO^1WQpMr2k&T;P)|@7yhmF7QTW3cQ9Z{7OsT; zR5~XS4VUJeTc<34H->NG~*)Bfc@^HKI=lN90FZySn~ zUkNFXq_7~R1(urgSPz1x!O|K*pILh}?Q)z%b|o>(uv%utx8pdw_U7{|fBE$X&1deL zH1$9{h3_@gMtxUk+LNUz3rl(us7qlHf8ot~qg6;fJiA4-fB|>V{|}6)*HFgHCIf(1 zM3E1tL=bdLjOMNYmeUMfGg-`P0e+vv)v4tS_0C@)!y>W;>|`En&U@W4gLcG4zG5eJjR`VZAI84T(^^FTXej&j*H-Q`U6_y@&Hbxn z;O!Mrp~v&w=Q)M@+xW?}5&vvL1RVm22r@uj%7JZWN z?7@tNMf7MxDeKNwF_1Ki`c7g5qPV^zGU1x<#)QSv z#}O7U?`;gYgMOKMaIcB(KbMuN4&=v5Fw$^oUTnHyWZEVh@2ZZn%|rZl)Rd_7FcF=` zm1{Z;NWV4RUBN<>*icJYH|`Oe%9L&6n`tMAn#=XXY3o>D4evge!r-Ba)C;6*TI4U2 zfwzl7Y|hN~xfjY_5D#BSodv&Ahpqv0E@*`7WZ>4(dvQ(X}E=Xz0$ zaKQ@*?n7Vb+PA(cM}Lk*1dlHv?+r5Wc3GfWyzG3Qnt}?MJb}C^|g=Qoz z*~$=TGH|WihST(@>HWW^LD5-6uM6n&SWFc}?>}>>Vw{5=; zK|-vFV`KNb99B7zSK&Rf9ajO@x~;byhQ3AAXsGybg$%q!vKDxrwLcHaT4utiW+#tN z+lg~f&|!G+;i&1XoH9>5eSCc4_|s3D;}fS&K6&cI3F*T}$q zUY7Q`?9%$70cFRJhdC8v2Ke$yVoYWG;<0@jNYlbGi^y-#$N}_KGH|~pNXN2BTih4I zD!}61O{iK}<=O&rzK${bBxM`}NyXuY_>E%UI$@Kz4sAVQEPMZm@o*izziV|JgM?;@ z7{vf(#Wo#^=^AHQoi`enn_K+9U4D@HWr*TkC=EI3nuUojI#sj!L^}-$M?;FFS=Wu7g=}0= zuv5GnU1K5ZpNg8g!n?1?1{5B1aC?)Do5Xrev$9Nv#dT1a20EJEVa!UEK8fO;g4Sza!Q9&47_bL)c>+N)KHXlQ@#%z>D>i&W_U~Pg%bp1Z5a< z3PWKHhbCIB;a&{DT>n3TG?x#yNeyA+Kw8%?m^^9CK-#c~u)pXI`U>^nuEN}$&5)3t zQ{~aMe0gIyov>50W#OQM0#_C|<&|@xZYHN4(V5!DXgcw2f~8r4VG+Hlbub}*3q<2{ z;Yq)Hq2l;{)~9i#duV+6@KIBKF%L~lA2yG?bbR)<6Co8`?k0zb+k_^pVOL~5>i8p9 z7e=JbJ78p9{FzyG;!5N#B~l-HmWY0jYu%=DMhmiR(4L>5-vH*4z+C!o!(0}a%by75 z4K2i&4A7+V+Go=ZPFzu&vLj`QACo|2EEETKSo~<&j?BuU9id8yY}a=d!tUY(R0ze6 z-4ghlqC?MQ;kPop)WHrYbSZ>A7;smEXM>bnJ!^}5CO{7D#3O4_87iXu1)SNUFLdqO zS`NCPk=9LjacmXWYs;ghFtL4eu3lp=yX<%ii%FMIv4*-yhQ;SR?1!Gra`&{h6<%#_Y7Fz;)lXphjP0~$)1{Oz-ud@~>7uWZ@vrEj-k+Vpf!-dEOu6$q` zCvn}U6qFv6zoXOUw}i%Jp>ZZlo$Kj54kuJHf@wx z>ai1X+W941ykFOCI-j<1(xmOtB{J|fCkx4W~WP4ade7j{h|W z-X;TYtHObtiycNCE11JM`1fx0r>+o*1ZtOI5}s5;T{uo<** ztukF@N2}YRR+j@q>{i&{kW8bnS2B(5>C%vC^x`iw^yhcKGff~L%l{Z-k!sN2TT=}K z127pEJtQ2i%5YWr%_T`JpuOYAroEcEnC`cY&!cs!!wAxv;Wu(0J|7djtAeLgJ$F*l zUlNNlvn@3&a{R#@yQf_q{P>^Qk4>Ng{RbYy-2(ss0AT4T-Zf|LGWnVBZa&K>R000310F=UaHIxA8`B_r{ literal 0 HcmV?d00001 From be443da37a9e7b312fe2a8aa828265d6acd513dc Mon Sep 17 00:00:00 2001 From: echenard Date: Fri, 7 Aug 2026 16:21:57 -0500 Subject: [PATCH 6/7] Fixed issues and improve structure of noise sampler --- CaloMC/CMakeLists.txt | 2 - CaloMC/fcl/prolog.fcl | 2 +- CaloMC/src/CaloDigiMaker_module.cc | 10 +- CaloMC/src/CaloHitTruthMatch_module.cc | 13 ++- CaloMC/src/ShowerStepUtil.cc | 8 +- CaloReco/CMakeLists.txt | 1 + CaloReco/fcl/prolog.fcl | 3 +- CaloReco/inc/CaloTemplateWFProcessor.hh | 4 +- CaloReco/inc/CaloTemplateWFUtil.hh | 2 +- CaloReco/src/CaloTemplateWFProcessor.cc | 2 +- Mu2eUtilities/CMakeLists.txt | 2 +- Mu2eUtilities/inc/CaloNoiseUtil.hh | 25 ++--- Mu2eUtilities/inc/CaloPulseUtil.hh | 4 +- Mu2eUtilities/src/CaloNoiseUtil.cc | 139 ++++++++++++------------ Mu2eUtilities/src/CaloPulseUtil.cc | 2 +- 15 files changed, 110 insertions(+), 109 deletions(-) diff --git a/CaloMC/CMakeLists.txt b/CaloMC/CMakeLists.txt index fa7dd8373f..13d1012da4 100644 --- a/CaloMC/CMakeLists.txt +++ b/CaloMC/CMakeLists.txt @@ -86,8 +86,6 @@ cet_build_plugin(CaloShowerUpdater art::module configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/prolog.fcl ${CURRENT_BINARY_DIR} fcl/prolog.fcl) -install(DIRECTORY data DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/Offline/CaloReco) - install_source(SUBDIRS src) install_headers(USE_PROJECT_NAME SUBDIRS inc) install_fhicl(SUBDIRS fcl SUBDIRNAME Offline/CaloMC/fcl) diff --git a/CaloMC/fcl/prolog.fcl b/CaloMC/fcl/prolog.fcl index 4574661b67..96c7ecadb4 100644 --- a/CaloMC/fcl/prolog.fcl +++ b/CaloMC/fcl/prolog.fcl @@ -64,7 +64,7 @@ CaloMC : { caloShowerSimCollection : CaloShowerROMaker caloHitCollection : CaloHitMaker primaryParticle : FindMCPrimary - pulseCache : { @table::CaloPulseCache } + pulseCache : @local::CaloPulseCache digiSampling : @local::HitMakerDigiSampling minAmplitude : 2 fillDetailedMC : false diff --git a/CaloMC/src/CaloDigiMaker_module.cc b/CaloMC/src/CaloDigiMaker_module.cc index 8a1092c1b0..f8181761ae 100644 --- a/CaloMC/src/CaloDigiMaker_module.cc +++ b/CaloMC/src/CaloDigiMaker_module.cc @@ -10,8 +10,6 @@ #include "fhiclcpp/types/Atom.h" #include "fhiclcpp/types/Sequence.h" #include "fhiclcpp/types/Table.h" -#include "art_root_io/TFileService.h" -#include "art_root_io/TFileDirectory.h" #include "Offline/CalorimeterGeom/inc/Calorimeter.hh" #include "Offline/CaloConditions/inc/CalSimParams.hh" @@ -73,7 +71,7 @@ namespace mu2e { pbtmcTag_ (config().pbtmcTag()), digitizationStart_ (config().digitizationStart()), digitizationEnd_ (config().digitizationEnd()), - digiSampling_ (config().digiSampling()), + digiSampling_ (config().pulseCache().digiSampling()), bufferDigi_ (config().bufferDigi()), startTimeBuffer_ (config().digiSampling()*config().bufferDigi()), maxADCCounts_ ((1 << config().nBits()) -1), @@ -93,7 +91,7 @@ namespace mu2e { //check that StartTimeBuffer is shorter than BlindTime_ if (startTimeBuffer_ > digitizationStart_) throw cet::exception("CALODIGIMAKER") << "digitizationStart too small to accommodate start time buffer"; - } + } void produce(art::Event& e) override; void beginRun(art::Run& aRun) override; @@ -180,6 +178,8 @@ namespace mu2e { mu2e::GeomHandle ch; calorimeter_ = ch.get(); + const int NoiseWFID(0); // will get this from proditions later; + if (calorimeter_->nCrystals()<1 || calorimeter_->G4Info().get("nSiPMPerCrystal")<1) return; int waveformSize = (digitizationEnd_ - digitizationStart_ + startTimeBuffer_) / digiSampling_; if (ewMarker.spillType() != EventWindowMarker::SpillType::onspill) { @@ -201,7 +201,7 @@ namespace mu2e { const double scaleFactor = readoutScaleFactor(iRO, calCrystalConds); generateSpotNoise(waveform,scaleFactor); } - buildOutputDigi(iRO, waveform, noiseSampler_.pedestal(), caloDigiColl); + buildOutputDigi(iRO, waveform, noiseSampler_.pedestal(NoiseWFID), caloDigiColl); } } diff --git a/CaloMC/src/CaloHitTruthMatch_module.cc b/CaloMC/src/CaloHitTruthMatch_module.cc index a78c1cec58..a8ca9a5fde 100644 --- a/CaloMC/src/CaloHitTruthMatch_module.cc +++ b/CaloMC/src/CaloHitTruthMatch_module.cc @@ -53,7 +53,7 @@ namespace mu2e { caloHitToken_ {consumes(config().caloHitCollection())}, ppToken_ {consumes(config().primaryParticle())}, pulseCache_ (config().pulseCache()), - digiSampling_ (config().digiSampling()), + digiSampling_ (config().pulseCache().digiSampling()), deltaTimeMinus_ (config().deltaTimeMinus()), minAmplitude_ (config().minAmplitude()), fillDetailedMC_ (config().fillDetailedMC()), @@ -160,15 +160,20 @@ namespace mu2e { for (std::size_t ihit=0; ihit < caloHits.size(); ++ihit) { const CaloHit& hit = caloHits[ihit]; const auto& sortedHits = caloHitMap.at(hit.crystalID()); - const auto& sortedSims = caloShowerSimsMap.at(hit.crystalID()); const auto hitIt = std::find(sortedHits.begin(),sortedHits.end(),&hit); const auto hitNextIt = std::next(hitIt); const auto hitPtr = art::Ptr(caloHitHandle,ihit); - if (diagLevel_ > 2) - for (const auto* shower : sortedSims) + const auto& sortedSimsIt = caloShowerSimsMap.find(hit.crystalID()); + const auto& sortedSims = sortedSimsIt == caloShowerSimsMap.end() ? + std::vector{} : sortedSimsIt->second; + + if (diagLevel_ > 2){ + if (sortedSims.empty()) std::cout<<"No shower sims for "<time()<<" / "<energyDep()< imax_) throw cet::exception("Rethrow") + if (i > imax_) throw cet::exception("SHOWERSTEPUTIL") << "[CaloMC/ShowerStepUtil] Index out of bound " << i << "\n"; @@ -37,7 +37,7 @@ namespace mu2e { void ShowerStepUtil::reset(unsigned i) { - if (i > imax_) throw cet::exception("Rethrow") + if (i > imax_) throw cet::exception("SHOWERSTEPUTIL") << "[CaloMC/ShowerStepUtil] Index out of bound " << i << "\n"; @@ -51,7 +51,7 @@ namespace mu2e { CLHEP::Hep3Vector& ShowerStepUtil::pos(unsigned i) { - if (i > imax_) throw cet::exception("Rethrow") + if (i > imax_) throw cet::exception("SHOWERSTEPUTIL") << "[CaloMC/ShowerStepUtil] Index out of bound " << i << "\n"; @@ -64,7 +64,7 @@ namespace mu2e { void ShowerStepUtil::printBucket(unsigned i) { - if (i > imax_) throw cet::exception("Rethrow") + if (i > imax_) throw cet::exception("SHOWERSTEPUTIL") << "[CaloMC/ShowerStepUtil] Index out of bound " << i << "\n"; std::cout<<"Entries= "< chiThreshold { Name("chiThreshold"), Comment("Min chi2 for refit strategy") }; fhicl::Atom refitLeadingEdge { Name("refitLeadingEdge"), Comment("Refit the leading edge to extract peak time") }; - fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; fhicl::Atom fitPrintLevel { Name("fitPrintLevel"), Comment("minuit fit print level") }; fhicl::Atom fitStrategy { Name("fitStrategy"), Comment("Minuit fit strategy") }; fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diagnosis level") }; @@ -72,7 +71,8 @@ namespace mu2e { void initHistos (); void setPrimaryPeakPar1 (const std::vector& xvec, const std::vector& yvec); void setPrimaryPeakPar2 (const std::vector& xvec, const std::vector& yvec); - void findRisingPeak (int ipeak, std::vector& parInit, const std::vector& xvec, const std::vector& yvec, std::vector& ywork); + void findRisingPeak (int ipeak, std::vector& parInit, const std::vector& xvec, + const std::vector& yvec, std::vector& ywork); void setSecondaryPeakPar(const std::vector& xvec, const std::vector& yvec); double estimatePeakTime (const std::vector& xvec, const std::vector& ywork, int ic); bool checkPeakDist (double x0); diff --git a/CaloReco/inc/CaloTemplateWFUtil.hh b/CaloReco/inc/CaloTemplateWFUtil.hh index b6b2dcfaa0..b4ae06c3f2 100644 --- a/CaloReco/inc/CaloTemplateWFUtil.hh +++ b/CaloReco/inc/CaloTemplateWFUtil.hh @@ -12,7 +12,7 @@ namespace mu2e { public: CaloTemplateWFUtil(const CaloPulseUtil::Config& configPulseCache, double minPeakAmplitude, - double minDTPeaks, int printLevel=-1); + double minDTPeaks, int printLevel); void initialize (); void setXYVector (const std::vector& xvec, const std::vector& yvec); diff --git a/CaloReco/src/CaloTemplateWFProcessor.cc b/CaloReco/src/CaloTemplateWFProcessor.cc index ff72e1ab16..11b15835d6 100644 --- a/CaloReco/src/CaloTemplateWFProcessor.cc +++ b/CaloReco/src/CaloTemplateWFProcessor.cc @@ -26,7 +26,7 @@ namespace mu2e { chiThreshold_ (config.chiThreshold()), refitLeadingEdge_(config.refitLeadingEdge()), diagLevel_ (config.diagLevel()), - fmutil_ (config.pulseCache(),minDTPeaks_,config.fitPrintLevel()), + fmutil_ (config.pulseCache(),minPeakAmplitude_,minDTPeaks_,config.fitPrintLevel()), chi2_ (999.), ndf_ (-1), resAmp_ (), diff --git a/Mu2eUtilities/CMakeLists.txt b/Mu2eUtilities/CMakeLists.txt index e3c1dd1579..d93cad8744 100644 --- a/Mu2eUtilities/CMakeLists.txt +++ b/Mu2eUtilities/CMakeLists.txt @@ -2,8 +2,8 @@ cet_make_library( SOURCE src/BinnedSpectrum.cc src/BuildLinearFitMatrixSums.cc - src/CaloPulseUtil.cc src/CaloNoiseUtil.cc + src/CaloPulseUtil.cc src/checkSimParticleCollection.cc src/compressPdgId.cc src/ConversionSpectrum.cc diff --git a/Mu2eUtilities/inc/CaloNoiseUtil.hh b/Mu2eUtilities/inc/CaloNoiseUtil.hh index 81cce59d9b..3674481d8d 100644 --- a/Mu2eUtilities/inc/CaloNoiseUtil.hh +++ b/Mu2eUtilities/inc/CaloNoiseUtil.hh @@ -1,15 +1,10 @@ -#ifndef CaloNoiseUtil_HH -#define CaloNoiseUtil_HH +#ifndef Mu2eUtilities_CaloNoiseUtil_hh +#define Mu2eUtilities_CaloNoiseUtil_hh // // Cache and provide noise waveforms for readouts // #include "fhiclcpp/types/Atom.h" -#include "fhiclcpp/types/Sequence.h" -#include "art/Framework/Services/Optional/RandomNumberGenerator.h" -#include "Offline/SeedService/inc/SeedService.hh" - #include "Offline/Mu2eUtilities/inc/CaloPulseUtil.hh" - #include "CLHEP/Random/RandPoissonQ.h" #include "CLHEP/Random/RandGaussQ.h" #include "CLHEP/Random/RandFlat.h" @@ -43,11 +38,11 @@ namespace mu2e { CaloNoiseUtil(const Config& config, CLHEP::HepRandomEngine& engine); - void prepare(int histoID, double peToADC); - std::span noiseSegment(int histoID, size_t istart, size_t ilength); - int pedestal(); - void printCache(); - void dumpNoise(const std::string& name, const std::vector& wave); + void prepare(int histoID, double peToADC); + std::span noiseSegment(int histoID, size_t istart, size_t ilength); + int pedestal(int histoID) const; + void printCache() const; + void dumpNoise(const std::string& name, const std::vector& wave); private: @@ -60,15 +55,13 @@ namespace mu2e { double digiSampling_; double noiseRinDark_; double noiseElec_; - double minPeakADC_; CLHEP::RandPoissonQ randPoisson_; CLHEP::RandGaussQ randGauss_; CLHEP::RandFlat randFlat_; CaloPulseUtil pulseCache_; bool dumpGenerated_; - int histoBaseID_; - int pedestal_; - std::map> noiseMap_; + std::map pedestal_; + std::map> noiseMap_; static constexpr int base = 10000; }; diff --git a/Mu2eUtilities/inc/CaloPulseUtil.hh b/Mu2eUtilities/inc/CaloPulseUtil.hh index 7575cbee4c..c38eae9ec5 100644 --- a/Mu2eUtilities/inc/CaloPulseUtil.hh +++ b/Mu2eUtilities/inc/CaloPulseUtil.hh @@ -1,5 +1,5 @@ -#ifndef CaloPulseUtil_HH -#define CaloPulseUtil_HH +#ifndef Mu2eUtilities_CaloPulseUtil_hh +#define Mu2eUtilities_CaloPulseUtil_hh // Calculate the values of the digitized pulse shape as a function of the hit time. // The value stored are the integral of the waveform over the digitization bin width. diff --git a/Mu2eUtilities/src/CaloNoiseUtil.cc b/Mu2eUtilities/src/CaloNoiseUtil.cc index fc35c64224..8e2f94d729 100644 --- a/Mu2eUtilities/src/CaloNoiseUtil.cc +++ b/Mu2eUtilities/src/CaloNoiseUtil.cc @@ -14,6 +14,8 @@ #include #include #include +#include +#include namespace mu2e { @@ -30,8 +32,7 @@ namespace mu2e { randFlat_ {engine}, pulseCache_ {config.pulseCache()}, dumpGenerated_ {config.dumpGenerated()}, - histoBaseID_ {-1}, - pedestal_ {0}, + pedestal_ {}, noiseMap_ {} {} @@ -39,25 +40,23 @@ namespace mu2e { //---------------------------------------------------------------------------------------------------------------------- void CaloNoiseUtil::prepare(int histoID, double peToADC) { - int baseID = histoID / base; - if (baseID == histoBaseID_) return; - if (generate_) generateCache(baseID, peToADC); - else fillCache(baseID); + if (generate_) generateCache(histoID, peToADC); } //---------------------------------------------------------------------------------------------------------------------- - void CaloNoiseUtil::fillCache(int histoBaseID) + void CaloNoiseUtil::fillCache(int histoID) { + int histoBaseID = histoID/base; + // Cache is for all baseID, clear it noiseMap_.clear(); - histoBaseID_ = histoBaseID; // Refill the cache with all histos sharing the same baseID ConfigFileLookupPolicy resolveFullPath; std::string fullFileName = resolveFullPath(fileName_); TFile file(fullFileName.c_str()); - if (!file.IsOpen()) throw cet::exception("NOISEREADER")<<"Filename"<> hid) || !iss.eof()) - throw cet::exception("NOISEREADER")<<"Hitsogram "<GetArray(); noiseMap_[hid].assign(array + 1, array + histo->GetNbinsX() + 1); @@ -89,92 +88,98 @@ namespace mu2e { // estimate pedestal, take a single value for every waveform // (pedestals will be stored somewhere else later) - float sum = 0.0; + double sum = 0.0; for (int i = 1; i <= histo->GetNbinsX(); ++i) sum += histo->GetBinContent(i); - pedestal_ = std::trunc(sum /histo->GetNbinsX() ); + pedestal_[hid] = std::trunc(sum /histo->GetNbinsX() ); } } //---------------------------------------------------------------------------------------------------------------------- - void CaloNoiseUtil::generateCache(int histoBaseID, double peToADC) + void CaloNoiseUtil::generateCache(int histoID, double peToADC) { - // Cache is for all baseID, clear it - noiseMap_.clear(); - histoBaseID_ = histoBaseID; - constexpr unsigned noiseSize{10000}; - std::vector waveform(noiseSize,0.0); - - pulseCache_.buildCache(); - const auto& pulse = pulseCache_.digitizedPulse(0.0); - const unsigned pulseSize = pulse.size(); - const unsigned bufferSize = int(0.75*pulseSize); - const double totalTime = (noiseSize+bufferSize)*digiSampling_; - const int noiseLevelPE = int(totalTime*noiseRinDark_); - - //Generate the radiation induced noise (RIN) - const int nPh = randPoisson_(noiseLevelPE); - for (int i=0;i waveform(noiseSize,0.0); + + pulseCache_.buildCache(); + const auto& pulse = pulseCache_.digitizedPulse(0.0); + const unsigned pulseSize = pulse.size(); + const unsigned bufferSize = int(0.75*pulseSize); + const double totalTime = (noiseSize+bufferSize)*digiSampling_; + const int noiseLevelPE = int(totalTime*noiseRinDark_); + + //Generate the radiation induced noise (RIN) + const int nPh = randPoisson_(noiseLevelPE); + for (int i=0;i CaloNoiseUtil::noiseSegment(int histoID, size_t istart, size_t ilength) + std::span CaloNoiseUtil::noiseSegment(int histoID, size_t istart, size_t ilength) { - int baseID = histoID/base; + auto iter = noiseMap_.find(histoID); + if (iter == noiseMap_.end()) { + + if (generate_) + throw cet::exception("CALONOISEUTIL") + << "noiseSegment called before prepare() for histoID " << histoID << "\n"; - if (baseID != histoBaseID_) { - if (generate_) throw cet::exception("CaloNoiseUtil") - << "noiseSegment called before prepare() for baseID " << histoID/base << "\n"; - else fillCache(baseID); + fillCache(histoID); + iter = noiseMap_.find(histoID); + if (iter == noiseMap_.end()) + throw cet::exception("CALONOISEUTIL")<<"histoID "<second; - if (ilength >= vec.size()) throw cet::exception("CALONOISEUTIL")<<"noise length request too long\n"; size_t irandom = size_t(randFlat_.fire(0.,vec.size()-ilength)); - return std::span(vec.data() + irandom, ilength); + return std::span(vec.data() + irandom, ilength); } //---------------------------------------------------------------------------------------------------------------------- - int CaloNoiseUtil::pedestal() {return pedestal_;} + int CaloNoiseUtil::pedestal(int histoID) const { + auto iter = pedestal_.find(histoID); + if (iter == pedestal_.end()) + throw cet::exception("CALONOISEUTIL")<<"phistoID "<second;; + } //---------------------------------------------------------------------------------------------------------------------- - void CaloNoiseUtil::printCache() + void CaloNoiseUtil::printCache() const { std::cout<<"CaloNoiseUtil cache\n"; for (const auto& kv : noiseMap_) std::cout<<"Histo id "<& wave) + void CaloNoiseUtil::dumpNoise(const std::string& fname, const std::vector& wave) { TFile outfile(fname.c_str(), "RECREATE"); diff --git a/Mu2eUtilities/src/CaloPulseUtil.cc b/Mu2eUtilities/src/CaloPulseUtil.cc index 4bd9c2c9ef..1eb546c01d 100644 --- a/Mu2eUtilities/src/CaloPulseUtil.cc +++ b/Mu2eUtilities/src/CaloPulseUtil.cc @@ -44,7 +44,7 @@ namespace mu2e { std::unique_ptr pshape(nullptr); TFile pulseFile(fullFileName.c_str()); if (pulseFile.IsOpen()) pshape.reset((TH1F*) pulseFile.Get(histName_.c_str())); - if (!pshape) throw cet::exception("CATEGORY")<<"CaloPulseUtil:: Hitsogram "<SetDirectory(0); pulseFile.Close(); From 8d9ef253c85aade17b0f5f6e02ace60f7851984c Mon Sep 17 00:00:00 2001 From: echenard Date: Sat, 8 Aug 2026 19:38:22 -0500 Subject: [PATCH 7/7] Further bug fixes --- CaloMC/fcl/prolog.fcl | 2 +- CaloMC/src/CaloDigiMaker_module.cc | 33 +++++++++++----------- CaloMC/src/CaloHitTruthMatch_module.cc | 12 ++++---- CaloMC/test/RunCaloCalibGun.fcl | 2 -- Mu2eUtilities/inc/CaloNoiseUtil.hh | 1 + Mu2eUtilities/src/CaloNoiseUtil.cc | 39 +++++++++++++++----------- 6 files changed, 48 insertions(+), 41 deletions(-) diff --git a/CaloMC/fcl/prolog.fcl b/CaloMC/fcl/prolog.fcl index 96c7ecadb4..52712956bd 100644 --- a/CaloMC/fcl/prolog.fcl +++ b/CaloMC/fcl/prolog.fcl @@ -50,10 +50,10 @@ CaloMC : { addNoise : true digitizationStart : @local::HitMakerDigitizationStart digitizationEnd : @local::HitMakerDigitizationEnd - digiSampling : @local::HitMakerDigiSampling bufferDigi : 16 nBits : 12 minPeakADC : @local::HitMakerMinPeakADC + minNoiseAmp : 2 nBinsPeak : 2 diagLevel : 0 } diff --git a/CaloMC/src/CaloDigiMaker_module.cc b/CaloMC/src/CaloDigiMaker_module.cc index f8181761ae..8d2e9443e0 100644 --- a/CaloMC/src/CaloDigiMaker_module.cc +++ b/CaloMC/src/CaloDigiMaker_module.cc @@ -56,10 +56,10 @@ namespace mu2e { fhicl::Atom digitizationStart { Name("digitizationStart"), Comment("Start of digitization window relative to nominal pb time") }; fhicl::Atom digitizationEnd { Name("digitizationEnd"), Comment("End of digitization window relative to nominal pb time")}; fhicl::Atom addNoise { Name("addNoise"), Comment("Add noise to waveform") }; - fhicl::Atom digiSampling { Name("digiSampling"), Comment("Digitization time sampling") }; fhicl::Atom nBits { Name("nBits"), Comment("ADC Number of bits") }; fhicl::Atom nBinsPeak { Name("nBinsPeak"), Comment("Window size for finding local maximum to digitize wf") }; fhicl::Atom minPeakADC { Name("minPeakADC"), Comment("Minimum ADC hits of local peak to digitize") }; + fhicl::Atom minNoiseAmp { Name("minNoiseAmp"), Comment("Minimum noise amplitude to add noise") }; fhicl::Atom bufferDigi { Name("bufferDigi"), Comment("Number of timeStamps for the buffer digi") }; fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level"),0 }; }; @@ -73,11 +73,12 @@ namespace mu2e { digitizationEnd_ (config().digitizationEnd()), digiSampling_ (config().pulseCache().digiSampling()), bufferDigi_ (config().bufferDigi()), - startTimeBuffer_ (config().digiSampling()*config().bufferDigi()), + startTimeBuffer_ (config().pulseCache().digiSampling()*config().bufferDigi()), maxADCCounts_ ((1 << config().nBits()) -1), pulseCache_ (CaloPulseUtil(config().pulseCache())), nBinsPeak_ (config().nBinsPeak()), minPeakADC_ (config().minPeakADC()), + minNoiseAmp_ (config().minNoiseAmp()), engine_ (createEngine(art::ServiceHandle()->getSeed())), addNoise_ (config().addNoise()), noiseSampler_ (config().noiseCache(), engine_), @@ -101,7 +102,7 @@ namespace mu2e { void makeDigitization (const CaloShowerROCollection&, CaloDigiCollection&, const EventWindowMarker&, const ProtonBunchTimeMC&, const CalSimParams&); bool fillROHits (unsigned iRO, std::vector& waveform, const CaloShowerROCollection&, const ProtonBunchTimeMC&, const CalSimParams&); - void generateSpotNoise (std::vector& waveform, double scaleFactor); + void AddWFNoise (std::vector& waveform, int noiseWFID); void buildOutputDigi (unsigned iRO, std::vector& waveform, double pedestal, CaloDigiCollection&); void extract (const std::vector& wf, std::vector& starts, std::vector& stops) const; double readoutScaleFactor(unsigned iRO, const CalSimParams& conds) const; @@ -122,6 +123,7 @@ namespace mu2e { CaloPulseUtil pulseCache_; unsigned nBinsPeak_; int minPeakADC_; + float minNoiseAmp_; CLHEP::HepRandomEngine& engine_; bool addNoise_; CaloNoiseUtil noiseSampler_; @@ -178,8 +180,6 @@ namespace mu2e { mu2e::GeomHandle ch; calorimeter_ = ch.get(); - const int NoiseWFID(0); // will get this from proditions later; - if (calorimeter_->nCrystals()<1 || calorimeter_->G4Info().get("nSiPMPerCrystal")<1) return; int waveformSize = (digitizationEnd_ - digitizationStart_ + startTimeBuffer_) / digiSampling_; if (ewMarker.spillType() != EventWindowMarker::SpillType::onspill) { @@ -195,35 +195,36 @@ namespace mu2e { if (resetWaveform) std::fill(waveform.begin(), waveform.end(), 0.0); bool isEmpty = fillROHits(iRO, waveform, CaloShowerROs, pbtmc, calCrystalConds ); resetWaveform = !isEmpty; - if (isEmpty) continue; + + double pedestal(0); if (addNoise_) { + const int NoiseWFID(0); // will get this from proditions later; const double scaleFactor = readoutScaleFactor(iRO, calCrystalConds); - generateSpotNoise(waveform,scaleFactor); + noiseSampler_.prepare(NoiseWFID, scaleFactor); + + AddWFNoise(waveform,NoiseWFID); + pedestal = noiseSampler_.pedestal(NoiseWFID); } - buildOutputDigi(iRO, waveform, noiseSampler_.pedestal(NoiseWFID), caloDigiColl); + + buildOutputDigi(iRO, waveform, pedestal, caloDigiColl); } } //---------------------------------------------------------------------------------------------------------- - void CaloDigiMaker::generateSpotNoise(std::vector& waveform, double scaleFactor) + void CaloDigiMaker::AddWFNoise(std::vector& waveform, int NoiseWFID) { - const int NoiseWFID(0); // will get this from proditions later; - - const double minAmplitude(2); - noiseSampler_.prepare(NoiseWFID, scaleFactor); - size_t timeSample(0); std::vector hitStarts{}, hitStops{}; hitStarts.reserve(16);hitStops.reserve(16); // First, find the ranges in the waveform with non-zero bins. while (timeSample < waveform.size()) { - if (waveform[timeSample] < minAmplitude) {++timeSample; continue;} + if (waveform[timeSample] < minNoiseAmp_) {++timeSample; continue;} size_t sampleStart = (timeSample > bufferDigi_) ? timeSample - bufferDigi_ : 0; size_t sampleStop(timeSample); - while (sampleStop < waveform.size() && waveform[sampleStop] > minAmplitude) ++sampleStop; + while (sampleStop < waveform.size() && waveform[sampleStop] > minNoiseAmp_) ++sampleStop; hitStarts.push_back(sampleStart); hitStops.push_back(sampleStop); diff --git a/CaloMC/src/CaloHitTruthMatch_module.cc b/CaloMC/src/CaloHitTruthMatch_module.cc index a8ca9a5fde..732ccc91f4 100644 --- a/CaloMC/src/CaloHitTruthMatch_module.cc +++ b/CaloMC/src/CaloHitTruthMatch_module.cc @@ -164,15 +164,15 @@ namespace mu2e { const auto hitNextIt = std::next(hitIt); const auto hitPtr = art::Ptr(caloHitHandle,ihit); - const auto& sortedSimsIt = caloShowerSimsMap.find(hit.crystalID()); - const auto& sortedSims = sortedSimsIt == caloShowerSimsMap.end() ? - std::vector{} : sortedSimsIt->second; + static const std::vector noSims; + const auto it = caloShowerSimsMap.find(hit.crystalID()); + const auto& sortedSims = (it == caloShowerSimsMap.end()) ? noSims : it->second; if (diagLevel_ > 2){ if (sortedSims.empty()) std::cout<<"No shower sims for "<time()<<" / "<energyDep()<time()<<" / "<energyDep()< pedestal_; + std::map peToADC_; std::map> noiseMap_; static constexpr int base = 10000; diff --git a/Mu2eUtilities/src/CaloNoiseUtil.cc b/Mu2eUtilities/src/CaloNoiseUtil.cc index 8e2f94d729..6bfd98c361 100644 --- a/Mu2eUtilities/src/CaloNoiseUtil.cc +++ b/Mu2eUtilities/src/CaloNoiseUtil.cc @@ -33,6 +33,7 @@ namespace mu2e { pulseCache_ {config.pulseCache()}, dumpGenerated_ {config.dumpGenerated()}, pedestal_ {}, + peToADC_ {}, noiseMap_ {} {} @@ -41,14 +42,19 @@ namespace mu2e { void CaloNoiseUtil::prepare(int histoID, double peToADC) { if (generate_) generateCache(histoID, peToADC); + else fillCache(histoID); } //---------------------------------------------------------------------------------------------------------------------- void CaloNoiseUtil::fillCache(int histoID) { - int histoBaseID = histoID/base; + //need to check that the cache is not already constructed for a histoID + if (noiseMap_.find(histoID) != noiseMap_.end()) return; + +std::cout<<"Fill Cache noise\n"; // Cache is for all baseID, clear it + int histoBaseID = histoID/base; noiseMap_.clear(); // Refill the cache with all histos sharing the same baseID @@ -98,6 +104,15 @@ namespace mu2e { //---------------------------------------------------------------------------------------------------------------------- void CaloNoiseUtil::generateCache(int histoID, double peToADC) { + // If a histo is alredy there with the correct peToADC, then return. If the + // histo is there but the peToADC is missing or different this is an error + if (noiseMap_.find(histoID) != noiseMap_.end()){ + auto it = peToADC_.find(histoID); + if (it != peToADC_.end() && std::abs(it->second - peToADC)<0.01) return; + throw cet::exception("CALONOISEUTIL")<<"Same histoID can't have two differnt peToADC\n"; + } + +std::cout<<"Generate Cache noise\n"; // Clear from cache all non basedID entries int histoBaseID = histoID/base; std::erase_if(noiseMap_, [&](const auto& pair) {return pair.first/base != histoBaseID;}); @@ -127,12 +142,13 @@ namespace mu2e { //add electronics noise double noiseADC = noiseElec_*digiSampling_*peToADC; for (auto& val : waveform) val += randGauss_.fire(0.0,noiseADC); - - noiseMap_.emplace(histoBaseID, waveform); + noiseMap_[histoID] = std::move(waveform); //estimate pedestal for this waveform - set it to theoretical value for the time being pedestal_[histoID] = std::trunc(noiseRinDark_*digiSampling_*std::accumulate(pulse.begin(),pulse.end(),0.0)*peToADC); + peToADC_[histoID] = peToADC; + if (dumpGenerated_){ dumpNoise("caloNoise.root",waveform); dumpGenerated_ = false; @@ -143,17 +159,8 @@ namespace mu2e { std::span CaloNoiseUtil::noiseSegment(int histoID, size_t istart, size_t ilength) { auto iter = noiseMap_.find(histoID); - if (iter == noiseMap_.end()) { - - if (generate_) - throw cet::exception("CALONOISEUTIL") - << "noiseSegment called before prepare() for histoID " << histoID << "\n"; - - fillCache(histoID); - iter = noiseMap_.find(histoID); - if (iter == noiseMap_.end()) - throw cet::exception("CALONOISEUTIL")<<"histoID "<second; if (ilength >= vec.size()) @@ -167,8 +174,8 @@ namespace mu2e { int CaloNoiseUtil::pedestal(int histoID) const { auto iter = pedestal_.find(histoID); if (iter == pedestal_.end()) - throw cet::exception("CALONOISEUTIL")<<"phistoID "<second;; + throw cet::exception("CALONOISEUTIL")<<"histoID "<second; } //----------------------------------------------------------------------------------------------------------------------