diff --git a/CaloMC/CMakeLists.txt b/CaloMC/CMakeLists.txt index 9f9468b0a3..13d1012da4 100644 --- a/CaloMC/CMakeLists.txt +++ b/CaloMC/CMakeLists.txt @@ -1,13 +1,12 @@ cet_make_library( SOURCE - src/CaloNoiseSimGenerator.cc src/CaloPhotonPropagation.cc - src/CaloWFExtractor.cc src/ShowerStepUtil.cc src/CaloDigiWrapper.cc src/CaloDigiWrapperCollection.cc LIBRARIES PUBLIC + Offline::ConfigTools Offline::CalorimeterGeom Offline::GeometryService Offline::Mu2eUtilities diff --git a/CaloMC/fcl/prolog.fcl b/CaloMC/fcl/prolog.fcl index 8db3b1935d..52712956bd 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 : @local::CaloPulseCache + noiseCache : @local::CaloNoiseCache + addNoise : true digitizationStart : @local::HitMakerDigitizationStart digitizationEnd : @local::HitMakerDigitizationEnd - addNoise : true - addRandomNoise : false - NoiseGenerator : { @table::CaloMC.NoiseGenerator } - digiSampling : @local::HitMakerDigiSampling + bufferDigi : 16 nBits : 12 minPeakADC : @local::HitMakerMinPeakADC + minNoiseAmp : 2 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 : @local::CaloPulseCache digiSampling : @local::HitMakerDigiSampling minAmplitude : 2 fillDetailedMC : false 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..9b3e4ab322 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 { @@ -13,22 +16,23 @@ 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); - 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/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/CaloClusterTruthMatch_module.cc b/CaloMC/src/CaloClusterTruthMatch_module.cc index a9ad2b339a..b8057dcabb 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,79 @@ 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) - { - const CaloCluster* thisCaloCluster = &cluster; - size_t idx = (thisCaloCluster - caloClusterBase); - art::Ptr clusterPtr = art::Ptr(caloClusterHandle,idx); - const auto& hits = cluster.caloHitsPtrVector(); + for (std::size_t idx=0; idx(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,22 +46,20 @@ namespace mu2e { { using Name = fhicl::Name; using Comment = fhicl::Comment; - using CNG = mu2e::CaloNoiseSimGenerator::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 minNoiseAmplitude { Name("minNoiseAmplitude"), Comment("minimum amplitude to start adding noise") }; 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 }; }; @@ -82,17 +71,17 @@ 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()), + startTimeBuffer_ (config().pulseCache().digiSampling()*config().bufferDigi()), maxADCCounts_ ((1 << config().nBits()) -1), - minNoiseAmplitude_ (config().minNoiseAmplitude()), - pulseShape_ (CaloPulseShape(config().pulseFileName(),config().pulseHistName(),config().digiSampling())), - wfExtractor_ (config().bufferDigi(),config().nBinsPeak(),config().minPeakADC(),config().bufferDigi()), + pulseCache_ (CaloPulseUtil(config().pulseCache())), + nBinsPeak_ (config().nBinsPeak()), + minPeakADC_ (config().minPeakADC()), + minNoiseAmp_ (config().minNoiseAmp()), engine_ (createEngine(art::ServiceHandle()->getSeed())), addNoise_ (config().addNoise()), - noiseGenerator_ (config().noise_gen_conf(), engine_, 0), - addRandomNoise_ (config().addRandomNoise()), + noiseSampler_ (config().noiseCache(), engine_), diagLevel_ (config().diagLevel()) { consumes(ewMarkerTag_); @@ -101,8 +90,9 @@ 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; void beginRun(art::Run& aRun) override; @@ -110,16 +100,16 @@ 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 generateSpotNoise (std::vector& waveform); - void buildOutputDigi (unsigned iRO, std::vector& waveform, double pedestal, CaloDigiCollection&); - 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); + 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 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; + void diag0 (unsigned, const std::vector&); + void diag1 (unsigned, double, size_t, const std::vector&, int); - ProditionsHandle calCrystalConds_; + ProditionsHandle calCrystalConds_; const art::ProductToken caloShowerToken_; art::InputTag ewMarkerTag_; art::InputTag pbtmcTag_; @@ -130,13 +120,13 @@ namespace mu2e { unsigned bufferDigi_; float startTimeBuffer_; int maxADCCounts_; - float minNoiseAmplitude_; - CaloPulseShape pulseShape_; - CaloWFExtractor wfExtractor_; + CaloPulseUtil pulseCache_; + unsigned nBinsPeak_; + int minPeakADC_; + float minNoiseAmp_; CLHEP::HepRandomEngine& engine_; bool addNoise_; - CaloNoiseSimGenerator noiseGenerator_; - bool addRandomNoise_; + CaloNoiseUtil noiseSampler_; const Calorimeter* calorimeter_; int diagLevel_; }; @@ -145,8 +135,7 @@ namespace mu2e { //----------------------------------------------------------------------------- void CaloDigiMaker::beginRun(art::Run& aRun) { - pulseShape_.buildShapes(); - if (addNoise_) noiseGenerator_.initialize(wfExtractor_); + pulseCache_.buildCache(); } @@ -166,8 +155,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(); @@ -175,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, calCrystalConds, ewMarker, pbtmc); + makeDigitization(CaloShowerROs, *caloDigiColl,ewMarker, pbtmc,calCrystalConds); event.put(std::move(caloDigiColl)); if ( diagLevel_ > 0 ) std::cout<<"[CaloDigiMaker::produce] end" << std::endl; @@ -186,94 +175,56 @@ 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, const CalSimParams& calCrystalConds) { mu2e::GeomHandle ch; calorimeter_ = ch.get(); 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, const CaloShowerROCollection& CaloShowerROs, - const ProtonBunchTimeMC& pbtmc, const CalSimParams& calCrystalConds) - { - 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; - - 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; + buildOutputDigi(iRO, waveform, pedestal, caloDigiColl); + } } - //---------------------------------------------------------------------------------------------------------- - void CaloDigiMaker::generateSpotNoise(std::vector& waveform) + void CaloDigiMaker::AddWFNoise(std::vector& waveform, int NoiseWFID) { 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] < minNoiseAmplitude_) {++timeSample; continue;} + while (timeSample < waveform.size()) { + if (waveform[timeSample] < minNoiseAmp_) {++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] > minNoiseAmp_) ++sampleStop; hitStarts.push_back(sampleStart); hitStops.push_back(sampleStop); @@ -282,8 +233,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; @@ -293,12 +243,40 @@ 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& 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) @@ -313,15 +291,14 @@ namespace mu2e { //extract hits start / stop times std::vector 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 wfsample{}; wfsample.reserve(sampleStop-sampleStart); @@ -341,8 +318,53 @@ 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()){ + // 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()); + } + + //------------------------------------------------------------------------------------------------------------------- + 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; + } //------------------------------------------------------------------------------------------------------------------- void CaloDigiMaker::diag0(unsigned iSiPM, const std::vector& wf) @@ -359,36 +381,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 #include +#include +#include +#include #include @@ -36,26 +34,27 @@ 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 }; + 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 }; }; 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()), + pulseCache_ (config().pulseCache()), + digiSampling_ (config().pulseCache().digiSampling()), + deltaTimeMinus_ (config().deltaTimeMinus()), minAmplitude_ (config().minAmplitude()), fillDetailedMC_ (config().fillDetailedMC()), diagLevel_ (config().diagLevel()) @@ -65,71 +64,39 @@ 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_; + CaloPulseUtil::Config pulseCache_; 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(); + CaloPulseUtil cps(pulseCache_); + cps.buildCache(); - wf_ = cps.digitizedPulse(0); + wf_ = cps.digitizedPulse(0); wfBinMax_ = std::distance(wf_.begin(),std::max_element(wf_.begin(),wf_.end())); } @@ -137,12 +104,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,162 +119,132 @@ 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); - - 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();}); - - - // do the matching for a given calo hit - for (size_t ihit=0; ihit < caloHits.size(); ++ihit) - { - const auto& hit = caloHits[ihit]; - const auto sortedHits = caloHitMap[hit.crystalID()]; - const auto sortedSims = caloShowerSimsMap[hit.crystalID()]; - auto hitIt = std::find(sortedHits.begin(),sortedHits.end(),&hit); - auto hitNextIt = std::next(hitIt); - auto hitPtr = art::Ptr(caloHitHandle,ihit); - auto showerIt = sortedSims.begin(); - auto showerItEnd = sortedSims.end(); - - if (diagLevel_ > 2) - for (const auto& shower : sortedSims) std::cout<<"[CaloHitTruthMatch] Sim shower id/time/energy="<crystalID() - <<" / "<time()<<" / "<energyDep()<> 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); + + // 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; itime()- (*hitIt)->time() - 2*digiSampling_ < deltaTimePlus ) - { - deltaTimePlus = (*hitNextIt)->time()- (*hitIt)->time() - 2*digiSampling_; + // matching per calo hit + caloHitMCs.reserve(caloHits.size()); + for (std::size_t ihit=0; ihit < caloHits.size(); ++ihit) { + const CaloHit& hit = caloHits[ihit]; + const auto& sortedHits = caloHitMap.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); + + 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()< 2) std::cout<<"[CaloHitTruthMatch] inspect hit id/time/energy/length "<time() - hit.time() - 2*digiSampling_ < deltaTimePlus) + deltaTimePlus = (*hitNextIt)->time() - hit.time() - 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; fillEdeps(primaryParticle, edeps, showerSim); - 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); + if (fillDetailedMC_) { + 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())); + caloHitMCs.emplace_back(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; - - 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 found the add caloShowerSim to matching caloEdep object, otherwise create new caloEdep - if (it!= edeps.end()) - { - it->addEDep(showerSim->energyDep()); + // 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 (it != edeps.end()) { + 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 - { + 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; } - 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/CaloNoiseSimGenerator.cc b/CaloMC/src/CaloNoiseSimGenerator.cc deleted file mode 100644 index 6b03fbe8ce..0000000000 --- a/CaloMC/src/CaloNoiseSimGenerator.cc +++ /dev/null @@ -1,178 +0,0 @@ -#include "Offline/CaloMC/inc/CaloNoiseSimGenerator.hh" -#include "art_root_io/TFileService.h" -#include "art_root_io/TFileDirectory.h" -#include "art/Framework/Services/Optional/RandomNumberGenerator.h" -#include "Offline/SeedService/inc/SeedService.hh" - -#include "TFile.h" -#include "TH2.h" -#include "TGraph.h" -#include "TCanvas.h" - -#include -#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..93e949736d 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/CaloShowerROMaker_module.cc b/CaloMC/src/CaloShowerROMaker_module.cc index 29b48d4138..ac7ec79420 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,116 +188,106 @@ 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(); if (LRUCorrection_) edep_corr = LRUCorrection(posZ/crystalLength, edep_corr, lru); // Generate individual PEs and their arrival times - for (int i=0; i PETime(NPE,hitTime); if (addTravelTime_) - { - for (auto& time : PETime) time += photonProp_.propTimeSimu(crystalLength-posZ); - } - CaloShowerROs.push_back(CaloShowerRO(SiPMID,stepPtr,PETime)); + for (auto& t : PETime) t += photonProp_.propTimeSimu(crystalLength-posZ); - if (diagLevel_ > 2) std::cout<<"[CaloShowerROMaker::generatePE] SiPMID:"< 2) {std::cout<<"Time hit "< 2) { + std::cout<<"[CaloShowerROMaker] SiPMID:"< 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 +297,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,142 @@ 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 + // Detailed 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) { + inspectedSims.clear(); SimPtr sim = step.simParticle(); - - SimParticlePtrCollection inspectedSims; - while (sim->hasParent() && cal.isInsideAnyCrystal(sim->startPosition()) ) + 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_); - - if (buffer.entries(idx)>0 && (step->time()-buffer.t0(idx) > deltaTime_) ) + const CLHEP::Hep3Vector pos = cal.mu2eToCrystal(volId, step->position()); + 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 +297,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< -#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()); - } -} - diff --git a/CaloMC/src/ShowerStepUtil.cc b/CaloMC/src/ShowerStepUtil.cc index 6904249e83..71ccd8e045 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("SHOWERSTEPUTIL") + << "[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("SHOWERSTEPUTIL") + << "[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("SHOWERSTEPUTIL") + << "[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= "< 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") }; @@ -42,7 +44,6 @@ namespace mu2e { fhicl::Atom psdThreshold { Name("psdThreshold"), Comment("Pulse shape discrimination threshold for secondary peaks") }; fhicl::Atom 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") }; @@ -70,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 01687c829f..b4ae06c3f2 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); 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..11b15835d6 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(),minPeakAmplitude_,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..d93cad8744 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/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 new file mode 100644 index 0000000000..072de97bab --- /dev/null +++ b/Mu2eUtilities/inc/CaloNoiseUtil.hh @@ -0,0 +1,71 @@ +#ifndef Mu2eUtilities_CaloNoiseUtil_hh +#define Mu2eUtilities_CaloNoiseUtil_hh +// +// Cache and provide noise waveforms for readouts +// +#include "fhiclcpp/types/Atom.h" +#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(int histoID) const; + void printCache() const; + 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_; + CLHEP::RandPoissonQ randPoisson_; + CLHEP::RandGaussQ randGauss_; + CLHEP::RandFlat randFlat_; + CaloPulseUtil pulseCache_; + bool dumpGenerated_; + std::map pedestal_; + std::map peToADC_; + std::map> noiseMap_; + + static constexpr int base = 10000; + }; + +} +#endif diff --git a/Mu2eUtilities/inc/CaloPulseShape.hh b/Mu2eUtilities/inc/CaloPulseUtil.hh similarity index 66% rename from Mu2eUtilities/inc/CaloPulseShape.hh rename to Mu2eUtilities/inc/CaloPulseUtil.hh index 394a4a30e9..c38eae9ec5 100644 --- a/Mu2eUtilities/inc/CaloPulseShape.hh +++ b/Mu2eUtilities/inc/CaloPulseUtil.hh @@ -1,5 +1,5 @@ -#ifndef CaloPulseShape_HH -#define CaloPulseShape_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. @@ -15,21 +15,32 @@ // 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 // -// NOTE: uncomment the pline creation if the discontinuities in the second order derivative arising from the -// linear piecewise approxmiation are problematic for the minimization + +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" #include #include namespace mu2e { - class CaloPulseShape + class CaloPulseUtil { public: - CaloPulseShape(const std::string& fileName, const std::string& histName, double digiSampling); - ~CaloPulseShape() {}; - - void buildShapes(); + 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; diff --git a/Mu2eUtilities/src/CaloNoiseUtil.cc b/Mu2eUtilities/src/CaloNoiseUtil.cc new file mode 100644 index 0000000000..6bfd98c361 --- /dev/null +++ b/Mu2eUtilities/src/CaloNoiseUtil.cc @@ -0,0 +1,200 @@ +#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 +#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()}, + pedestal_ {}, + peToADC_ {}, + noiseMap_ {} + {} + + + //---------------------------------------------------------------------------------------------------------------------- + void CaloNoiseUtil::prepare(int histoID, double peToADC) + { + if (generate_) generateCache(histoID, peToADC); + else fillCache(histoID); + } + + //---------------------------------------------------------------------------------------------------------------------- + void CaloNoiseUtil::fillCache(int histoID) + { + //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 + ConfigFileLookupPolicy resolveFullPath; + std::string fullFileName = resolveFullPath(fileName_); + + TFile file(fullFileName.c_str()); + if (!file.IsOpen()) throw cet::exception("CALONOISEUTIL")<<"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("CALONOISEUTIL")<<"Histogram "<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) + + double sum = 0.0; + for (int i = 1; i <= histo->GetNbinsX(); ++i) sum += histo->GetBinContent(i); + pedestal_[hid] = std::trunc(sum /histo->GetNbinsX() ); + } + } + + + //---------------------------------------------------------------------------------------------------------------------- + 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;}); + + 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) + { + 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(int histoID) const { + auto iter = pedestal_.find(histoID); + if (iter == pedestal_.end()) + throw cet::exception("CALONOISEUTIL")<<"histoID "<second; + } + + //---------------------------------------------------------------------------------------------------------------------- + void CaloNoiseUtil::printCache() const + { + 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 @@ -12,7 +11,18 @@ namespace mu2e { - CaloPulseShape::CaloPulseShape(const std::string& fileName, const std::string& histName, double digiSampling) : + 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), @@ -24,9 +34,8 @@ namespace mu2e { {} //---------------------------------------------------------------------------------------------------------------------- - void CaloPulseShape::buildShapes() + void CaloPulseUtil::buildCache() { - pulseVec_.clear(); ConfigFileLookupPolicy resolveFullPath; @@ -35,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")<<"CaloPulseShape:: Hitsogram "<SetDirectory(0); pulseFile.Close(); @@ -55,13 +64,15 @@ namespace mu2e { deltaT_ = 0.0; // find difference between peak time and t0 for digitized waveform. - for (int i=1;i& CaloPulseShape::digitizedPulse(double hitTime) const + const std::vector& CaloPulseUtil::digitizedPulse(double hitTime) const { int shiftBin = nSteps_ - int(hitTime/digiStep_)%nSteps_; for (int i=0;i