From f4545fe3ecf8a9a922aeb343c9acae954fa4e9cc Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 13 Jan 2026 12:32:15 -0500 Subject: [PATCH 01/36] tentative debug tree --- .../PMT/Algorithms/PMTsimulationAlg.cxx | 54 +++++++- icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 41 +++++- icaruscode/PMT/SimPMTIcarus_module.cc | 121 +++++++++++++++++- 3 files changed, 200 insertions(+), 16 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 5e3ad9d3f..7f3762963 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -156,17 +156,24 @@ icarus::opdet::PMTsimulationAlg::PMTsimulationAlg // ----------------------------------------------------------------------------- -std::tuple, std::optional> +std::tuple, std::optional, + std::optional> icarus::opdet::PMTsimulationAlg::simulate(sim::SimPhotons const& photons, - sim::SimPhotonsLite const& lite_photons) + sim::SimPhotonsLite const& lite_photons, + bool enableDebug) { std::optional photons_used; + std::optional debug; + + if (enableDebug) debug.emplace(); - Waveform_t const waveform = CreateFullWaveform(photons, lite_photons, photons_used); + Waveform_t const waveform = CreateFullWaveform(photons, lite_photons, photons_used, + enableDebug ? &(*debug) : nullptr); return { CreateFixedSizeOpDetWaveforms(photons.OpChannel(), waveform), - std::move(photons_used) + std::move(photons_used), + std::move(debug) }; } // icarus::opdet::PMTsimulationAlg::simulate() @@ -191,7 +198,8 @@ auto icarus::opdet::PMTsimulationAlg::makeGainFluctuator() const { auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform (sim::SimPhotons const& photons, sim::SimPhotonsLite const& lite_photons, - std::optional& photons_used) + std::optional& photons_used, + icarus::opdet::DebugInfo* debug) const -> Waveform_t { @@ -217,6 +225,19 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform } microsecond const timeDelay { totalDelay }; + if (debug) { + debug->opChannel = channel; + debug->QE = fQE; + debug->sampling_MHz = fSampling.value(); + debug->readoutEnablePeriod_us = fParams.readoutEnablePeriod.value(); + debug->timeDelay_us = timeDelay.value(); + debug->triggerOffsetPMT_us = fParams.triggerOffsetPMT.value(); + debug->nSamples = fNsamples; + debug->nSubsamples = wsp.nSubsamples(); + debug->photons.clear(); + debug->peDeposits.clear(); + } + // // collect the amount of photoelectrons arriving at each subtick; // the waveform is split in groups of photons at the same relative subtick @@ -264,6 +285,16 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform if (tick >= endSample) continue; ++peMaps[subtick][tick]; + + if (debug) { + icarus::opdet::DebugPhoton dbgPh; + dbgPh.simTime_ns = ph.Time.value(); + dbgPh.trigTime_us = mytime.value(); + dbgPh.tick = tick.value(); + dbgPh.subtick = subtick.value(); + debug->photons.push_back(std::move(dbgPh)); + } + } // for photons // auto end = std::chrono::high_resolution_clock::now(); @@ -327,8 +358,18 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform ); // std::cout << "Channel=" << channel << ", subsample=" << iSubsample << ", tick=" << startTick << ", nPE=" << nPE << ", ePE=" << nEffectivePE << std::endl; + if (debug) { + DebugPEDeposit dep; + dep.tick = startTick.value(); + dep.subtick = iSubsample; + dep.nPE = nPE; + dep.nEffectivePE = nEffectivePE; + debug->peDeposits.push_back(std::move(dep)); + } + } // for sample } // for subsamples + MF_LOG_TRACE("PMTsimulationAlg") << nTotalPE << " photoelectrons at " << std::accumulate( @@ -361,6 +402,9 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform // end=std::chrono::high_resolution_clock::now(); diff = end-start; // std::cout << "\tadded saturation... " << channel << " " << diff.count() << std::endl; + if(debug) + debug.waveform = waveform; + return waveform; } // CreateFullWaveform() diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index a7e49a61c..19963fc05 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -80,9 +80,36 @@ namespace icarus::opdet { class PMTsimulationAlg; class PMTsimulationAlgMaker; - -} // namespace icarus::opdet + struct DebugPhoton { + float simTime_ns = 0.f; // input photon time (simulation) + float trigTime_us = 0.f; // timings.toTriggerTime(...) - offset + delay + int32_t tick = -1; // tick index (optical clock tick) + uint16_t subtick = 0; // subsample index + }; + + struct DebugPEDeposit { + int32_t tick = -1; // where the PE got deposited + uint16_t subtick = 0; + uint16_t nPE = 0; // integer PE count + float nEffectivePE = 0.f; // after gain fluctuation (what you actually used) + float gainFactor() const { return (nPE > 0) ? (nEffectivePE / float(nPE)) : 0.f; } + }; + + struct DebugInfo { + uint32_t opChannel = 0; + float QE = 0.f; + float sampling_MHz = 0.f; + float readoutEnablePeriod_us = 0.f; + float timeDelay_us = 0.f; + float triggerOffsetPMT_us = 0.f; + uint32_t nSamples = 0; + uint16_t nSubsamples = 0; + std::vector photons; // per-photon bookkeeping + std::vector peDeposits; // aggregated PE deposits actually used to build waveform + OpDetWaveformMakerClass::WaveformData_t waveform; // waveform data (ADC counts) + }; +} // namespace icarus::opdet // ----------------------------------------------------------------------------- /// Helper class to cut a `raw::OpDetWaveform` from a longer waveform data. @@ -559,9 +586,11 @@ class icarus::opdet::PMTsimulationAlg { * In that case, the returned `sim::SimPhotons` contains a copy of each of * the `photons` contributing to any of the waveforms. */ - std::tuple, std::optional> + std::tuple, std::optional, + std::optional> simulate(sim::SimPhotons const& photons, - sim::SimPhotonsLite const& lite_photons); + sim::SimPhotonsLite const& lite_photons, + bool enableDebug = false); /// Prints the configuration into the specified output stream. template @@ -673,8 +702,8 @@ class icarus::opdet::PMTsimulationAlg { Waveform_t CreateFullWaveform( sim::SimPhotons const& photons, sim::SimPhotonsLite const& lite_photons, - std::optional& photons_used - ) const; + std::optional& photons_used, + icarus::opdet::DebugInfo* debug) const; /** * @brief Creates `raw::OpDetWaveform` objects from a waveform data. diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 6da16d173..bd735abfc 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -30,6 +30,7 @@ #include "nurandom/RandomUtils/NuRandomService.h" // framework libraries +#include "art_root_io/TFileService.h" #include "art/Framework/Core/EDProducer.h" #include "art/Framework/Core/ModuleMacros.h" #include "art/Framework/Principal/Event.h" @@ -240,6 +241,12 @@ namespace icarus::opdet { "HepJamesRandom" }; + fhicl::Atom DebugTree { + Name("DebugTree"), + Comment("enable debug tree output"), + false + }; + }; // struct Config using Parameters = art::EDProducer::Table; @@ -254,6 +261,7 @@ namespace icarus::opdet { SimPMTIcarus & operator = (SimPMTIcarus &&) = delete; // Required functions. + virtual void beginJob() override; void produce(art::Event & e) override; private: @@ -271,7 +279,7 @@ namespace icarus::opdet { bool fMakeMetadata; ///< Whether to produce waveform metadata. bool fWritePhotons { false }; ///< Whether to save contributing photons. - + CLHEP::HepRandomEngine& fEfficiencyEngine; CLHEP::HepRandomEngine& fDarkNoiseEngine; CLHEP::HepRandomEngine& fElectronicsNoiseEngine; @@ -285,7 +293,6 @@ namespace icarus::opdet { /// The actual simulation algorithm. icarus::opdet::PMTsimulationAlgMaker makePMTsimulator; - /// True if `firstTime()` has already been called. std::atomic_flag fNotFirstTime; @@ -302,7 +309,35 @@ namespace icarus::opdet { /// Returns whether no other event has been processed yet. bool firstTime() { return !fNotFirstTime.test_and_set(); } - + + bool fDebugTree { false }; ///< Whether to enable debug tree output. + TTree* fTree { nullptr }; ///< Debug tree pointer. + // debug tree variables + int run, event, channel; + float QE; + float sampling_MHz; + float readoutEnablePeriod_us; + float timeDelay_us; + float triggerOffsetPMT_us; + int nSamples; + int nSubsamples; + // info per waveform + int nsize; + std::vector wf; + // info per photon + int nPhotons; + std::vector photon_simTime_ns; + std::vector photon_trigTime_us; + std::vector photon_tick; + std::vector photon_subtick; + // info per PE deposit + int nPEDeposits; + std::vector pedeposit_tick; + std::vector pedeposit_subtick; + std::vector pedeposit_nPE; + std::vector pedeposit_nEffectivePE; + std::vector pedeposit_gainFactor; + }; // class SimPMTIcarus @@ -346,6 +381,8 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) ->makeGenerator(fElectronicsNoiseEngine) } , makePMTsimulator(config().algoConfig()) + , fDebugTree(config().DebugTree()) + { // Call appropriate produces<>() functions here. produces>(); @@ -358,6 +395,38 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) fNotFirstTime.clear(); // superfluous in C++20 } // SimPMTIcarus::SimPMTIcarus() + // --------------------------------------------------------------------------- + void SimPMTIcarus::beginJob() + { + if (fDebugTree) { + art::ServiceHandle tfs; + fTree = tfs->make("SimPMTDebug","Debug tree for SimPMTIcarus"); + fTree->Branch("run", &run, "run/I"); + fTree->Branch("event", &event, "event/I"); + fTree->Branch("channel", &channel, "channel/I"); + fTree->Branch("nSamples", &nSamples, "nSamples/I"); + fTree->Branch("nSubsamples", &nSubsamples, "nSubsamples/I"); + fTree->Branch("QE", &QE, "QE/F"); + fTree->Branch("sampling_MHz", &sampling_MHz, "sampling_MHz/F"); + fTree->Branch("readoutEnablePeriod_us", &readoutEnablePeriod_us, "readoutEnablePeriod_us/F"); + fTree->Branch("timeDelay_us", &timeDelay_us, "timeDelay_us/F"); + fTree->Branch("triggerOffsetPMT_us", &triggerOffsetPMT_us, "triggerOffsetPMT_us/F"); + fTree->Branch("nsize", &nsize, "nsize/I"); + fTree->Branch("wf", &wf); + fTree->Branch("nPhotons", &nPhotons, "nPhotons/I"); + fTree->Branch("photon_simTime_ns", &photon_simTime_ns); + fTree->Branch("photon_trigTime_us", &photon_trigTime_us); + fTree->Branch("photon_tick", &photon_tick); + fTree->Branch("photon_subtick", &photon_subtick); + fTree->Branch("nPEDeposits", &nPEDeposits, "nPEDeposits/I"); + fTree->Branch("pedeposit_tick", &pedeposit_tick); + fTree->Branch("pedeposit_subtick", &pedeposit_subtick); + fTree->Branch("pedeposit_nPE", &pedeposit_nPE); + fTree->Branch("pedeposit_nEffectivePE", &pedeposit_nEffectivePE); + fTree->Branch("pedeposit_gainFactor", &pedeposit_gainFactor); + } + } // SimPMTIcarus::beginJob() + // --------------------------------------------------------------------------- void SimPMTIcarus::produce(art::Event & e) @@ -437,13 +506,55 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) sim::SimPhotons photons(lite_photons.OpChannel); - auto const& [ channelWaveforms, photons_used ] - = PMTsimulator->simulate(photons, lite_photons); + auto const& [ channelWaveforms, photons_used, debug] + = PMTsimulator->simulate(photons, lite_photons, fDebugTree); std::move( channelWaveforms.cbegin(), channelWaveforms.cend(), std::back_inserter(*pulseVecPtr) ); + if(fDebugTree && debug) { + // fill debug tree variables + channel = debug->channel; + QE = debug->QE; + sampling_MHz = debug->sampling_MHz; + readoutEnablePeriod_us = debug->readoutEnablePeriod_us; + timeDelay_us = debug->timeDelay_us; + triggerOffsetPMT_us = debug->triggerOffsetPMT_us; + nSamples = debug->nSamples; + nSubsamples = debug->nSubsamples; + //fill per waveform info + nsize = debug->waveform.size(); + wf = debug->waveform; + // fill per photon info + nPhotons = debug->photons.size(); + photon_simTime_ns.clear(); + photon_trigTime_us.clear(); + photon_tick.clear(); + photon_subtick.clear(); + for(auto const& phot : debug->photons) { + photon_simTime_ns.push_back(phot.simTime_ns); + photon_trigTime_us.push_back(phot.trigTime_us); + photon_tick.push_back(phot.tick); + photon_subtick.push_back(phot.subtick); + } + // fill per PE deposit info + nPEDeposits = debug->PEdeposits.size(); + pedeposit_tick.clear(); + pedeposit_subtick.clear(); + pedeposit_nPE.clear(); + pedeposit_nEffectivePE.clear(); + pedeposit_gainFactor.clear(); + for(auto const& pedep : debug->PEdeposits) { + pedeposit_tick.push_back(pedep.tick); + pedeposit_subtick.push_back(pedep.subtick); + pedeposit_nPE.push_back(pedep.nPE); + pedeposit_nEffectivePE.push_back(pedep.nEffectivePE); + pedeposit_gainFactor.push_back(pedep.gainFactor); + } + + fTree->Fill(); + } // if debug tree } // for } From 64a330409bc9d0d92906ec6f07e97e565320357c Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 13 Jan 2026 12:33:37 -0600 Subject: [PATCH 02/36] tweaks for success build --- icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx | 11 +++++++---- icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 2 +- icaruscode/PMT/SimPMTIcarus_module.cc | 16 +++++++++------- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 7f3762963..58fdb4091 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -288,10 +288,10 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform if (debug) { icarus::opdet::DebugPhoton dbgPh; - dbgPh.simTime_ns = ph.Time.value(); + dbgPh.simTime_ns = ph.Time; dbgPh.trigTime_us = mytime.value(); dbgPh.tick = tick.value(); - dbgPh.subtick = subtick.value(); + dbgPh.subtick = subtick; debug->photons.push_back(std::move(dbgPh)); } @@ -402,8 +402,11 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform // end=std::chrono::high_resolution_clock::now(); diff = end-start; // std::cout << "\tadded saturation... " << channel << " " << diff.count() << std::endl; - if(debug) - debug.waveform = waveform; + if(debug){ + for(std::size_t i=0; iwaveform.push_back(waveform[i].value()); + } + } return waveform; } // CreateFullWaveform() diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index 19963fc05..2f085efbf 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -107,7 +107,7 @@ namespace icarus::opdet { uint16_t nSubsamples = 0; std::vector photons; // per-photon bookkeeping std::vector peDeposits; // aggregated PE deposits actually used to build waveform - OpDetWaveformMakerClass::WaveformData_t waveform; // waveform data (ADC counts) + std::vector> waveform; // waveform data (ADC counts) }; } // namespace icarus::opdet diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index bd735abfc..8566cc93f 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -48,6 +48,9 @@ #include "fhiclcpp/types/Atom.h" #include "fhiclcpp/ParameterSet.h" +// ROOT libraries +#include "TTree.h" + // CLHEP libraries #include "CLHEP/Random/RandEngine.h" // CLHEP::HepRandomEngine @@ -60,7 +63,6 @@ #include // std::move() #include - namespace icarus::opdet { /** @@ -323,7 +325,7 @@ namespace icarus::opdet { int nSubsamples; // info per waveform int nsize; - std::vector wf; + std::vector wf; // info per photon int nPhotons; std::vector photon_simTime_ns; @@ -487,7 +489,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) sim::SimPhotonsLite lite_photons(photons.OpChannel()); - auto const& [ channelWaveforms, photons_used ] + auto const& [ channelWaveforms, photons_used, debug ] = PMTsimulator->simulate(photons, lite_photons); std::move( channelWaveforms.cbegin(), channelWaveforms.cend(), @@ -515,7 +517,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) if(fDebugTree && debug) { // fill debug tree variables - channel = debug->channel; + channel = debug->opChannel; QE = debug->QE; sampling_MHz = debug->sampling_MHz; readoutEnablePeriod_us = debug->readoutEnablePeriod_us; @@ -539,18 +541,18 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) photon_subtick.push_back(phot.subtick); } // fill per PE deposit info - nPEDeposits = debug->PEdeposits.size(); + nPEDeposits = debug->peDeposits.size(); pedeposit_tick.clear(); pedeposit_subtick.clear(); pedeposit_nPE.clear(); pedeposit_nEffectivePE.clear(); pedeposit_gainFactor.clear(); - for(auto const& pedep : debug->PEdeposits) { + for(auto const& pedep : debug->peDeposits) { pedeposit_tick.push_back(pedep.tick); pedeposit_subtick.push_back(pedep.subtick); pedeposit_nPE.push_back(pedep.nPE); pedeposit_nEffectivePE.push_back(pedep.nEffectivePE); - pedeposit_gainFactor.push_back(pedep.gainFactor); + pedeposit_gainFactor.push_back(pedep.gainFactor()); } fTree->Fill(); From 2dda045d2710bbc8a4d94c082a0f54477cd33dff Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 13 Jan 2026 14:16:33 -0500 Subject: [PATCH 03/36] move to other loop --- icaruscode/PMT/SimPMTIcarus_module.cc | 37 +++++++++++++-------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 8566cc93f..432820a21 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -490,31 +490,14 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) sim::SimPhotonsLite lite_photons(photons.OpChannel()); auto const& [ channelWaveforms, photons_used, debug ] - = PMTsimulator->simulate(photons, lite_photons); + = PMTsimulator->simulate(photons, lite_photons, fDebugTree); std::move( channelWaveforms.cbegin(), channelWaveforms.cend(), std::back_inserter(*pulseVecPtr) ); if (simphVecPtr && photons_used) simphVecPtr->emplace_back(std::move(photons_used.value())); - - } // for - } - else if(pmtLiteVector.isValid()) { - nopch = pmtLiteVector->size(); - for(auto const& lite_photons : *pmtLiteVector) { - - // Make an empty SimPhotons with the same channel number. - - sim::SimPhotons photons(lite_photons.OpChannel); - - auto const& [ channelWaveforms, photons_used, debug] - = PMTsimulator->simulate(photons, lite_photons, fDebugTree); - std::move( - channelWaveforms.cbegin(), channelWaveforms.cend(), - std::back_inserter(*pulseVecPtr) - ); - + if(fDebugTree && debug) { // fill debug tree variables channel = debug->opChannel; @@ -559,6 +542,22 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) } // if debug tree } // for } + else if(pmtLiteVector.isValid()) { + nopch = pmtLiteVector->size(); + for(auto const& lite_photons : *pmtLiteVector) { + + // Make an empty SimPhotons with the same channel number. + + sim::SimPhotons photons(lite_photons.OpChannel); + + auto const& [ channelWaveforms, photons_used, debug] + = PMTsimulator->simulate(photons, lite_photons); + std::move( + channelWaveforms.cbegin(), channelWaveforms.cend(), + std::back_inserter(*pulseVecPtr) + ); + } // for + } mf::LogInfo("SimPMTIcarus") << "Generated " << pulseVecPtr->size() << " waveforms out of " << nopch << " optical channels."; From 0bd7e7497c234507749b705b8c76880c6e463b77 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 13 Jan 2026 13:45:36 -0600 Subject: [PATCH 04/36] save debug only with photons --- icaruscode/PMT/SimPMTIcarus_module.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 432820a21..ed3d7aed3 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -498,7 +498,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) if (simphVecPtr && photons_used) simphVecPtr->emplace_back(std::move(photons_used.value())); - if(fDebugTree && debug) { + if(fDebugTree && debug && debug->photons.size()>0) { // fill debug tree variables channel = debug->opChannel; QE = debug->QE; From fc5e966e8d5080b7ae8876330477c269ac6bbf9b Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 13 Jan 2026 14:58:03 -0500 Subject: [PATCH 05/36] store photon start position in debug --- .../PMT/Algorithms/PMTsimulationAlg.cxx | 3 +++ icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 6 ++--- icaruscode/PMT/SimPMTIcarus_module.cc | 25 +++++++++++-------- 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 58fdb4091..83982ecbe 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -288,6 +288,9 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform if (debug) { icarus::opdet::DebugPhoton dbgPh; + dbgPh.startX = ph.InitialPosition.X(); + dbgPh.startY = ph.InitialPosition.Y(); + dbgPh.startZ = ph.InitialPosition.Z(); dbgPh.simTime_ns = ph.Time; dbgPh.trigTime_us = mytime.value(); dbgPh.tick = tick.value(); diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index 2f085efbf..e6068ca02 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -82,6 +82,9 @@ namespace icarus::opdet { class PMTsimulationAlgMaker; struct DebugPhoton { + float startX = 0.f; // photon start X position [cm] + float startY = 0.f; // photon start Y position [cm] + float startZ = 0.f; // photon start Z position [cm] float simTime_ns = 0.f; // input photon time (simulation) float trigTime_us = 0.f; // timings.toTriggerTime(...) - offset + delay int32_t tick = -1; // tick index (optical clock tick) @@ -98,9 +101,6 @@ namespace icarus::opdet { struct DebugInfo { uint32_t opChannel = 0; - float QE = 0.f; - float sampling_MHz = 0.f; - float readoutEnablePeriod_us = 0.f; float timeDelay_us = 0.f; float triggerOffsetPMT_us = 0.f; uint32_t nSamples = 0; diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index ed3d7aed3..78012c2d3 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -316,20 +316,19 @@ namespace icarus::opdet { TTree* fTree { nullptr }; ///< Debug tree pointer. // debug tree variables int run, event, channel; - float QE; - float sampling_MHz; - float readoutEnablePeriod_us; float timeDelay_us; float triggerOffsetPMT_us; int nSamples; int nSubsamples; // info per waveform - int nsize; std::vector wf; // info per photon int nPhotons; std::vector photon_simTime_ns; std::vector photon_trigTime_us; + std::vector photon_start_x; + std::vector photon_start_y; + std::vector photon_start_z; std::vector photon_tick; std::vector photon_subtick; // info per PE deposit @@ -408,14 +407,13 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) fTree->Branch("channel", &channel, "channel/I"); fTree->Branch("nSamples", &nSamples, "nSamples/I"); fTree->Branch("nSubsamples", &nSubsamples, "nSubsamples/I"); - fTree->Branch("QE", &QE, "QE/F"); - fTree->Branch("sampling_MHz", &sampling_MHz, "sampling_MHz/F"); - fTree->Branch("readoutEnablePeriod_us", &readoutEnablePeriod_us, "readoutEnablePeriod_us/F"); fTree->Branch("timeDelay_us", &timeDelay_us, "timeDelay_us/F"); fTree->Branch("triggerOffsetPMT_us", &triggerOffsetPMT_us, "triggerOffsetPMT_us/F"); - fTree->Branch("nsize", &nsize, "nsize/I"); fTree->Branch("wf", &wf); fTree->Branch("nPhotons", &nPhotons, "nPhotons/I"); + fTree->Branch("photon_start_x", &photon_start_x); + fTree->Branch("photon_start_y", &photon_start_y); + fTree->Branch("photon_start_z", &photon_start_z); fTree->Branch("photon_simTime_ns", &photon_simTime_ns); fTree->Branch("photon_trigTime_us", &photon_trigTime_us); fTree->Branch("photon_tick", &photon_tick); @@ -501,9 +499,8 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) if(fDebugTree && debug && debug->photons.size()>0) { // fill debug tree variables channel = debug->opChannel; - QE = debug->QE; - sampling_MHz = debug->sampling_MHz; - readoutEnablePeriod_us = debug->readoutEnablePeriod_us; + run = e.id().run(); + event = e.id().event(); timeDelay_us = debug->timeDelay_us; triggerOffsetPMT_us = debug->triggerOffsetPMT_us; nSamples = debug->nSamples; @@ -517,7 +514,13 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) photon_trigTime_us.clear(); photon_tick.clear(); photon_subtick.clear(); + photon_start_x.clear(); + photon_start_y.clear(); + photon_start_z.clear(); for(auto const& phot : debug->photons) { + photon_start_x.push_back(phot.startX); + photon_start_y.push_back(phot.startY); + photon_start_z.push_back(phot.startZ); photon_simTime_ns.push_back(phot.simTime_ns); photon_trigTime_us.push_back(phot.trigTime_us); photon_tick.push_back(phot.tick); From 69c8e178c06194df303b51f17680436172d83c20 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 13 Jan 2026 15:01:08 -0500 Subject: [PATCH 06/36] remove old variables --- icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 83982ecbe..435158ce2 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -227,15 +227,13 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform if (debug) { debug->opChannel = channel; - debug->QE = fQE; - debug->sampling_MHz = fSampling.value(); - debug->readoutEnablePeriod_us = fParams.readoutEnablePeriod.value(); debug->timeDelay_us = timeDelay.value(); debug->triggerOffsetPMT_us = fParams.triggerOffsetPMT.value(); debug->nSamples = fNsamples; debug->nSubsamples = wsp.nSubsamples(); debug->photons.clear(); debug->peDeposits.clear(); + debug->waveform.clear(); } // From e490f278485a1407285d0a0c60258d16a316eea2 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 13 Jan 2026 15:22:49 -0500 Subject: [PATCH 07/36] remove forgotten variable --- icaruscode/PMT/SimPMTIcarus_module.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 78012c2d3..baea67355 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -506,7 +506,6 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) nSamples = debug->nSamples; nSubsamples = debug->nSubsamples; //fill per waveform info - nsize = debug->waveform.size(); wf = debug->waveform; // fill per photon info nPhotons = debug->photons.size(); From f2ddbc94aeb4e8c29967daaf797e9846c87736b1 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Thu, 12 Mar 2026 12:05:50 -0400 Subject: [PATCH 08/36] augment gain db service --- .../ICARUSPhotonCalibratorServiceFromDB.h | 4 ++ .../Calibration/PhotonCalibratorFromDB.cxx | 39 ++++++++++++++++--- .../PMT/Calibration/PhotonCalibratorFromDB.h | 19 ++++++--- 3 files changed, 51 insertions(+), 11 deletions(-) diff --git a/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h b/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h index e9c62a189..c4ef076be 100644 --- a/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h +++ b/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h @@ -58,6 +58,10 @@ namespace calib { std::string getAreaDatabaseTag() const& { return fProvider.getAreaDatabaseTag(); } + double getSPEArea(int channel) const { return fProvider.getSPEArea(channel); } + + double getSPEFitWidth(int channel) const { return fProvider.getSPEFitWidth(channel); } + private: provider_type const* provider() const override { return &fProvider; } diff --git a/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx b/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx index fd4b3c1bc..054dbe59f 100644 --- a/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx +++ b/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx @@ -71,16 +71,43 @@ void icarusDB::PhotonCalibratorFromDB::readCalibrationFromDB(unsigned int run) mf::LogDebug log(fLogCategory); log << "Loading calibration for " << channelList.size() << " channels in run " << run; for( auto channel : channelList ) { - - // SPE area from database + + // SPE area [ADC x tick] double area = 0; int error = fDB.GetNamedChannelData( channel, "area", area ); - if( error ) throw cet::exception( "PhotonCalibratorFromDB" ) << "Encountered error (code " << error << ") while trying to access 'area' from the table\n"; - + if( error ) throw cet::exception( fLogCategory ) << "Error (code " << error << ") reading 'area' for channel " << channel << "\n"; fDatabaseSPECalibrations[channel].speArea = area; // ADC x tick + + // SPE area uncertainty from fit + double areaErr = 0; + error = fDB.GetNamedChannelData( channel, "area_err", areaErr ); + if( error ) throw cet::exception( fLogCategory ) << "Error (code " << error << ") reading 'area_err' for channel " << channel << "\n"; + fDatabaseSPECalibrations[channel].speAreaErr = areaErr; // ADC x tick + + // SPE fit width (sigma) [ADC x tick] + double width = 0; + error = fDB.GetNamedChannelData( channel, "width", width ); + if( error ) throw cet::exception( fLogCategory ) << "Error (code " << error << ") reading 'width' for channel " << channel << "\n"; + fDatabaseSPECalibrations[channel].speFitWidth = width; // ADC x tick + + // SPE fit width uncertainty + double widthErr = 0; + error = fDB.GetNamedChannelData( channel, "width_err", widthErr ); + if( error ) throw cet::exception( fLogCategory ) << "Error (code " << error << ") reading 'width_err' for channel " << channel << "\n"; + fDatabaseSPECalibrations[channel].speFitWidthErr = widthErr; // ADC x tick + + // Fit quality + double fitQuality = 0; + error = fDB.GetNamedChannelData( channel, "fit_quality", fitQuality ); + if( error ) throw cet::exception( fLogCategory ) << "Error (code " << error << ") reading 'fit_quality' for channel " << channel << "\n"; + fDatabaseSPECalibrations[channel].speFitQuality = fitQuality; + if (fVerbose) - log << "\n channel " << std::setw(3) << channel << " SPE area " << area; - } + log << "\n channel " << std::setw(3) << channel + << " SPE area " << area << " +/- " << areaErr + << " sigma " << width << " +/- " << widthErr + << " fit quality " << fitQuality; + } } diff --git a/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.h b/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.h index c5a6026cf..38805be69 100644 --- a/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.h +++ b/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.h @@ -29,11 +29,12 @@ namespace icarusDB::details { struct PhotonCalibratorInfo { /// Area [positive, ADC count sum] of response to single photoelectron. - double speArea = -1.0; - double speAreaErr = -1.0; /// Uncertainty on `speArea` [ADC count sum] - double speFitWidth = -1.0; - double speFitWidthErr = -1.0; - + double speArea = -1.0; /// `speArea` from fit [ADC x tick] + double speAreaErr = -1.0; /// fit uncertainty on `speArea` [ADC x tick] + double speFitWidth = -1.0; /// width (sigma) from fit [ADC x tick] + double speFitWidthErr = -1.0; /// fit uncertainty on `speFitWidth` [ADC x tick] + double speFitQuality = -1.0; /// Fit quality + }; } // icarusDB::details @@ -151,6 +152,14 @@ class icarusDB::PhotonCalibratorFromDB: public calib::IPhotonCalibrator { /// Get current area database tag std::string getAreaDatabaseTag() const& { return fAreaTag; } + /// Returns the SPE area [ADC x tick] for the given channel (default if not in DB). + double getSPEArea(int channel) const + { return getChannelCalibOrDefault(channel).speArea; } + + /// Returns the SPE fit width (sigma) [ADC x tick] for the given channel (default if not in DB). + double getSPEFitWidth(int channel) const + { return getChannelCalibOrDefault(channel).speFitWidth; } + private: using PhotonCalibratorInfo = details::PhotonCalibratorInfo; From a7bb214c20ce32668d8263ca12d7089f71d2ae23 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Thu, 12 Mar 2026 12:40:01 -0400 Subject: [PATCH 09/36] add gain service provider interface --- .../PMT/Algorithms/PMTsimulationAlg.cxx | 11 +++++-- icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 29 +++++++++++++++++++ .../PMT/Algorithms/pmtsimulation_icarus.fcl | 6 ++-- icaruscode/PMT/SimPMTIcarus_module.cc | 4 +++ 4 files changed, 45 insertions(+), 5 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index c8e6f7e8d..c03e6ffe1 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -906,6 +906,7 @@ icarus::opdet::PMTsimulationAlgMaker::PMTsimulationAlgMaker (PMTspecs.VoltageDistribution()); fBaseConfig.PMTspecs.gain = PMTspecs.Gain(); fBaseConfig.doGainFluctuations = config.FluctuateGain(); + fBaseConfig.useGainCalibDB = config.UseGainDatabase(); fBaseConfig.doTimingDelays = config.ApplyTimingDelays(); // @@ -962,6 +963,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, + calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -972,7 +974,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( { return std::make_unique(makeParams( beamGateTimestamp, - larProp, clockData, timingDelays, + larProp, clockData, timingDelays, gainCalibService, SPRfunction, pedestalGenerator, mainRandomEngine, darkNoiseRandomEngine, elecNoiseRandomEngine, trackSelectedPhotons @@ -986,7 +988,8 @@ auto icarus::opdet::PMTsimulationAlgMaker::makeParams( std::uint64_t beamGateTimestamp, detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, - icarusDB::PMTTimingCorrections const* timingDelays, + icarusDB::PMTTimingCorrections const* timingDelays, + calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -1010,8 +1013,10 @@ auto icarus::opdet::PMTsimulationAlgMaker::makeParams( params.clockData = &clockData; params.detTimings = detinfo::makeDetectorTimings(params.clockData); params.timingDelays = timingDelays; - + params.gainCalibService = gainCalibService; + assert(!params.doTimingDelays || params.timingDelays); + assert(!params.useGainCalibDB || params.gainCalibService); params.pulseFunction = &SPRfunction; params.pedestalGen = &pedestalGenerator; diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index d24c19623..8e9bac83f 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -22,6 +22,7 @@ #include "icaruscode/Utilities/quantities_utils.h" // util::value_t #include "icarusalg/Utilities/SampledFunction.h" #include "icaruscode/Timing/PMTTimingCorrections.h" +#include "icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h" // LArSoft libraries #include "lardataobj/RawData/OpDetWaveform.h" @@ -224,6 +225,10 @@ class icarus::opdet::OpDetWaveformMakerClass { * to a nominal gain (`PMTspecs.gain` configuration parameter), which is * then fluctuated to obtain the effective gain. This feature can be * disabled by setting configuration parameter `FluctuateGain` to `false`. + * + * Two gain fluctuation models are available: + * + * **Legacy Poisson model** (default, `UseGainDatabase: false`): * The approximation used here is that the fluctuation is entirely due to * the first stage of multiplication. The gain on the first stage is * described as a random variable with Poisson distribution around the mean @@ -244,6 +249,19 @@ class icarus::opdet::OpDetWaveformMakerClass { * The first stage gain is computed by * `icarus::opdet::PMTsimulationAlg::ConfigurationParameters_t::PMTspecs_t::multiplicationStageGain()`. * + * **Database gain model** (`UseGainDatabase: true`): + * When this mode is enabled the algorithm reads, for each PMT channel, + * the measured SPE area and SPE area width from the calibration database via + * the `calib::ICARUSPhotonCalibratorServiceFromDB` service. + * The number of effective photoelectrons added per true photoelectron is + * drawn from a Gaussian whose mean is the ratio of the channel SPE area + * to the nominal SPR template area @f$ q_0 @f$, and whose standard + * deviation is the ratio of the channel SPE width to @f$ q_0 @f$. + * This model preserves the mean deposited charge per photoelectron to + * match the measured per-channel gain, while the nominal gain set in + * `PMTspecs` and in `SinglePhotonResponse` cancels out and can be left + * at any convenient value across run periods. + * * * Dark noise * ----------- @@ -483,6 +501,7 @@ class icarus::opdet::PMTsimulationAlg { float saturation; //equivalent to the number of p.e. that saturates the electronic signal PMTspecs_t PMTspecs; ///< PMT specifications. bool doGainFluctuations; ///< Whether to simulate gain fluctuations. + bool useGainCalibDB; ///< Whether to use per-channel DB gain for fluctuations. bool doTimingDelays; ///< Whether to simulate timing delays. /// @} @@ -499,6 +518,9 @@ class icarus::opdet::PMTsimulationAlg { /// Timing delays service interfacing with database icarusDB::PMTTimingCorrections const* timingDelays = nullptr; + /// SPE calibration service for per-channel gain (nullptr = disabled). + calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService = nullptr; + // detTimings is not really "optional" but it needs delayed construction. /// Detector clocks data wrapper. std::optional detTimings; @@ -1014,6 +1036,11 @@ class icarus::opdet::PMTsimulationAlgMaker { Comment("add timing delays (cable, transit time) to photon times"), true }; + fhicl::Atom UseGainDatabase { + Name("UseGainDatabase"), + Comment("use per-channel SPE area and width from calibration DB for gain fluctuations"), + true + }; // // single photoelectron response @@ -1101,6 +1128,7 @@ class icarus::opdet::PMTsimulationAlgMaker { detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& detClocks, icarusDB::PMTTimingCorrections const* timingDelays, + calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -1133,6 +1161,7 @@ class icarus::opdet::PMTsimulationAlgMaker { detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, + calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, diff --git a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl index 0082653a5..1957be5ca 100644 --- a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl +++ b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl @@ -386,7 +386,8 @@ icarus_pmtsimulationalg_202202_noise: { QE: @local::icarus_opticalproperties.ScintPreScale # from opticalproperties_icarus.fcl FluctuateGain: true # apply per-photoelectron gain fluctuations ApplyTimingDelays: false # legacy value - + UseGainDatabase: false # legacy value + PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage # is proportional to dV^k @@ -534,7 +535,8 @@ icarus_pmtsimulationalg_2018: { QE: @local::icarus_opticalproperties.ScintPreScale # from opticalproperties_icarus.fcl FluctuateGain: true # apply per-photoelectron gain fluctuations ApplyTimingDelays: false # legacy value - + UseGainDatabase: false # legacy value + PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage # is proportional to dV^k diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 9bf70fefe..03bb95b1b 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -17,6 +17,7 @@ #include "icaruscode/PMT/Algorithms/PhotoelectronPulseFunction.h" #include "icaruscode/IcarusObj/OpDetWaveformMeta.h" #include "icaruscode/Timing/IPMTTimingCorrectionService.h" +#include "icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h" // LArSoft libraries #include "larcore/CoreUtils/ServiceUtil.h" @@ -272,6 +273,7 @@ namespace icarus::opdet { bool fMakeMetadata; ///< Whether to produce waveform metadata. bool fWritePhotons { false }; ///< Whether to save contributing photons. bool fDoTimingDelays; ///< Whether timing delay corrections are applied. + bool fUseGainCalibDB; ///< Whether per-channel gain calibration from DB is applied. CLHEP::HepRandomEngine& fEfficiencyEngine; CLHEP::HepRandomEngine& fDarkNoiseEngine; @@ -317,6 +319,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) , fMakeMetadata(config().MakeMetadata()) , fWritePhotons(config().writePhotons()) , fDoTimingDelays(config().algoConfig().ApplyTimingDelays()) + , fUseGainCalibDB(config().algoConfig().UseGainDatabase()) // random engines , fEfficiencyEngine(art::ServiceHandle()->registerAndSeedEngine( createEngine(0, "HepJamesRandom", "Efficiencies"), @@ -394,6 +397,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) *(lar::providerFrom()), clockData, fDoTimingDelays ? lar::providerFrom() : nullptr, + fUseGainCalibDB ? &(*art::ServiceHandle()) : nullptr, *fSinglePhotonResponseFunc, *fPedestalGen, fEfficiencyEngine, From 60fc0359a961d7ad179aab4d666a1ebd6ef40e18 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Thu, 12 Mar 2026 14:26:14 -0500 Subject: [PATCH 10/36] fix linking for db interface --- icaruscode/PMT/Algorithms/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/icaruscode/PMT/Algorithms/CMakeLists.txt b/icaruscode/PMT/Algorithms/CMakeLists.txt index 01a8b7f25..e2ac0ec24 100644 --- a/icaruscode/PMT/Algorithms/CMakeLists.txt +++ b/icaruscode/PMT/Algorithms/CMakeLists.txt @@ -6,6 +6,7 @@ art_make_library( SOURCE ${lib_srcs} LIBRARIES icaruscode::Decode_DecoderTools + icaruscode::PMT_Calibration icarusalg::Utilities larcorealg::Geometry lardataobj::RawData From 384a5ef7ba5f498203c723453f54f924e511a1e6 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 16 Mar 2026 11:44:55 -0500 Subject: [PATCH 11/36] update speareas database tag --- icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl b/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl index 39a2d2efe..b2b4a1476 100644 --- a/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl +++ b/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl @@ -48,7 +48,7 @@ pmt_bkgphotons_calibration: icarus_photon_calibration: { service_provider: ICARUSPhotonCalibratorServiceFromDB - AreaTag: "v1r0" + AreaTag: "v1r1" Defaults: { SPEArea: @local::SPERun2.Area # ADC x tick, Run2 default } From 213b9e003a6063ba8abc2e8e2828c5a1762fdfcc Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 16 Mar 2026 15:37:12 -0400 Subject: [PATCH 12/36] add pmt channel status service --- icaruscode/PMT/CMakeLists.txt | 1 + icaruscode/PMT/Status/CMakeLists.txt | 26 ++++ .../PMT/Status/IPMTChannelStatusService.h | 32 +++++ icaruscode/PMT/Status/PMTChannelStatus.h | 86 +++++++++++ .../PMT/Status/PMTChannelStatusProvider.cxx | 133 ++++++++++++++++++ .../PMT/Status/PMTChannelStatusProvider.h | 112 +++++++++++++++ .../Status/PMTChannelStatusService_service.cc | 60 ++++++++ .../PMT/Status/pmt_channel_status_icarus.fcl | 18 +++ 8 files changed, 468 insertions(+) create mode 100644 icaruscode/PMT/Status/CMakeLists.txt create mode 100644 icaruscode/PMT/Status/IPMTChannelStatusService.h create mode 100644 icaruscode/PMT/Status/PMTChannelStatus.h create mode 100644 icaruscode/PMT/Status/PMTChannelStatusProvider.cxx create mode 100644 icaruscode/PMT/Status/PMTChannelStatusProvider.h create mode 100644 icaruscode/PMT/Status/PMTChannelStatusService_service.cc create mode 100644 icaruscode/PMT/Status/pmt_channel_status_icarus.fcl diff --git a/icaruscode/PMT/CMakeLists.txt b/icaruscode/PMT/CMakeLists.txt index a34d31b95..21084fab3 100644 --- a/icaruscode/PMT/CMakeLists.txt +++ b/icaruscode/PMT/CMakeLists.txt @@ -5,6 +5,7 @@ add_subdirectory(Algorithms) add_subdirectory(LibraryMappingTools) add_subdirectory(Trigger) add_subdirectory(Calibration) +add_subdirectory(Status) # Removing AVX instructions because assuming the grid nodes # being less than 6 year old proved to be pretentious. diff --git a/icaruscode/PMT/Status/CMakeLists.txt b/icaruscode/PMT/Status/CMakeLists.txt new file mode 100644 index 000000000..0b628976d --- /dev/null +++ b/icaruscode/PMT/Status/CMakeLists.txt @@ -0,0 +1,26 @@ +cet_enable_asserts() + +set( LIB_LIBRARIES + art::Framework_Services_Registry + messagefacility::MF_MessageLogger + larcorealg::CoreUtils + larevt::CalibrationDBI_IOVData + larevt::CalibrationDBI_Providers + fhiclcpp::fhiclcpp + cetlib_except::cetlib_except +) + +set( SERVICE_LIBRARIES + icaruscode_PMT_Status + larcore::CoreUtils +) + +file(GLOB lib_srcs *.cxx) + +art_make_library( SOURCE ${lib_srcs} LIBRARIES PUBLIC ${LIB_LIBRARIES}) + +cet_build_plugin( PMTChannelStatusService art::service LIBRARIES PUBLIC ${SERVICE_LIBRARIES}) + +install_headers() +install_fhicl() +install_source() diff --git a/icaruscode/PMT/Status/IPMTChannelStatusService.h b/icaruscode/PMT/Status/IPMTChannelStatusService.h new file mode 100644 index 000000000..e429278a7 --- /dev/null +++ b/icaruscode/PMT/Status/IPMTChannelStatusService.h @@ -0,0 +1,32 @@ +/** + * @file icaruscode/PMT/Status/IPMTChannelStatusService.h + * @brief Art service interface for PMT channel status. + * @author Matteo Vicenzi (mvicenzi@bnl.gov) + */ + +#ifndef ICARUSCODE_PMT_STATUS_IPMTCHANNELSTATUSSERVICE_H +#define ICARUSCODE_PMT_STATUS_IPMTCHANNELSTATUSSERVICE_H + +// ICARUS libraries +#include "icaruscode/PMT/Status/PMTChannelStatus.h" + +// LArSoft libraries +#include "larcore/CoreUtils/ServiceProviderWrappers.h" + + +// ----------------------------------------------------------------------------- +namespace icarusDB { + /// The only thing this service does is to return its service provider of type + /// `icarusDB::PMTChannelStatus`. + using IPMTChannelStatusService + = lar::ServiceProviderInterfaceWrapper; +} + + +// ----------------------------------------------------------------------------- +DECLARE_ART_SERVICE_INTERFACE(icarusDB::IPMTChannelStatusService, SHARED) + + +// ----------------------------------------------------------------------------- + +#endif // ICARUSCODE_PMT_STATUS_IPMTCHANNELSTATUSSERVICE_H diff --git a/icaruscode/PMT/Status/PMTChannelStatus.h b/icaruscode/PMT/Status/PMTChannelStatus.h new file mode 100644 index 000000000..8708b4959 --- /dev/null +++ b/icaruscode/PMT/Status/PMTChannelStatus.h @@ -0,0 +1,86 @@ +/** + * @file icaruscode/PMT/Status/PMTChannelStatus.h + * @brief Interface for PMT channel status provider. + * @author Matteo Vicenzi (mvicenzi@bnl.gov) + */ + +#ifndef ICARUSCODE_PMT_STATUS_PMTCHANNELSTATUS_H +#define ICARUSCODE_PMT_STATUS_PMTCHANNELSTATUS_H + +// LArSoft libraries +#include "larcorealg/CoreUtils/UncopiableAndUnmovableClass.h" + +// Framework libraries +#include "art/Framework/Services/Registry/ServiceDeclarationMacros.h" + +// C++ standard libraries +#include +#include + + +namespace icarusDB { + + /// Possible status values for a PMT channel, matching the database integers. + enum PMTChannelStatusValue : int { + kOFF = 0, ///< Channel is off (not powered). + kON = 1, ///< Channel is on and good. + kBad = 2 ///< Channel is powered but flagged as bad. + }; + + + /** + * @brief Interface for PMT channel status information. + * + * Currently, the class provides interface for the following information: + * - channel status: kON, kOFF, or kBad + * - convenience predicates: isOn(), isOff(), isBad() + * - channel sets by status: getOnChannels(), getOffChannels(), getBadChannels() + * + */ + class PMTChannelStatus : private lar::UncopiableAndUnmovableClass { + + public: + + using ChannelSet_t = std::set; + + virtual ~PMTChannelStatus() noexcept = default; + + /// Returns the status of the specified channel. + virtual PMTChannelStatusValue getChannelStatus(unsigned int channel) const = 0; + + /// Returns whether the specified channel is on and good. + virtual bool isOn(unsigned int channel) const + { return getChannelStatus(channel) == kON; } + + /// Returns whether the specified channel is off (not powered). + virtual bool isOff(unsigned int channel) const + { return getChannelStatus(channel) == kOFF; } + + /// Returns whether the specified channel is powered but bad. + virtual bool isBad(unsigned int channel) const + { return getChannelStatus(channel) == kBad; } + + /// Returns the set of channels with status kON. + virtual ChannelSet_t getOnChannels() const = 0; + + /// Returns the set of channels with status kOFF. + virtual ChannelSet_t getOffChannels() const = 0; + + /// Returns the set of channels with status kBad. + virtual ChannelSet_t getBadChannels() const = 0; + + /// Returns the nominal voltage [V] of the specified channel. + virtual double getChannelVoltage(unsigned int channel) const = 0; + + /// Returns the database tag currently in use. + virtual std::string getDatabaseTag() const = 0; + + }; // class PMTChannelStatus + +} // namespace icarusDB + + +DECLARE_ART_SERVICE_INTERFACE(icarusDB::PMTChannelStatus, SHARED) + + +#endif // ICARUSCODE_PMT_STATUS_PMTCHANNELSTATUS_H diff --git a/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx b/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx new file mode 100644 index 000000000..a46c3eea0 --- /dev/null +++ b/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx @@ -0,0 +1,133 @@ +/** + * @file icaruscode/PMT/Status/PMTChannelStatusProvider.cxx + * @brief Implementation of PMT channel status provider. + * @author Matteo Vicenzi (mvicenzi@bnl.gov) + * @see icaruscode/PMT/Status/PMTChannelStatusProvider.h + */ + +#include "icaruscode/PMT/Status/PMTChannelStatusProvider.h" + +// LArSoft libraries +#include "larevt/CalibrationDBI/IOVData/TimeStampDecoder.h" + +// Framework libraries +#include "messagefacility/MessageLogger/MessageLogger.h" + +// C++ standard libraries +#include // std::setw() + + +// ----------------------------------------------------------------------------- +icarusDB::PMTChannelStatusProvider::PMTChannelStatusProvider + (const fhicl::ParameterSet& pset) + : fVerbose ( pset.get ("Verbose", false ) ) + , fLogCategory( pset.get("LogCategory", "PMTChannelStatusProvider") ) + , fStatusTag ( pset.get("StatusTag", "" ) ) + , fDB ( pset.get("DBname", "pmt_voltage_data"), + "", "", fStatusTag, true, false ) +{ + int const defaultStatusInt = pset.get("DefaultStatus", static_cast(kON)); + fDefault.status = static_cast(defaultStatusInt); + fDefault.voltage = pset.get("DefaultVoltage", 1500.0); + + mf::LogInfo(fLogCategory) + << "PMTChannelStatusProvider connected to '" + << pset.get("DBname", "pmt_voltage_data") + << "' DB, tag '" << fStatusTag << "'" + << ", default status: " << defaultStatusInt + << ", default voltage: " << fDefault.voltage << " V"; +} + + +// ----------------------------------------------------------------------------- +void icarusDB::PMTChannelStatusProvider::readStatusFromDB(unsigned int run) +{ + mf::LogInfo(fLogCategory) + << "Reading PMT channel statuses from database for run " << run; + + bool const ret = fDB.UpdateData( RunToDatabaseTimestamp(run) ); + mf::LogTrace(fLogCategory) + << "Status" << (ret ? "" : " not") << " updated for run " << run + << "\nFetched IoV [ " << fDB.CachedStart().DBStamp() + << " ; " << fDB.CachedEnd().DBStamp() + << " ] to cover t=" << RunToDatabaseTimestamp(run) + << " [=" << lariov::TimeStampDecoder::DecodeTimeStamp(RunToDatabaseTimestamp(run)).DBStamp() << "]"; + + std::vector channelList; + if (int const res = fDB.GetChannelList(channelList); res != 0) { + throw cet::exception(fLogCategory) + << "GetChannelList() returned " << res << " on run " << run << " query\n"; + } + + if (channelList.empty()) { + throw cet::exception(fLogCategory) + << "Got an empty channel list for run " << run << "\n"; + } + + fDatabaseStatus.clear(); + + mf::LogDebug log(fLogCategory); + log << "Loading status for " << channelList.size() << " channels in run " << run; + + for (auto const channel : channelList) { + + long statusInt = 0; + int error = fDB.GetNamedChannelData(channel, "status", statusInt); + if (error) { + throw cet::exception(fLogCategory) + << "Error (code " << error << ") reading 'status' for channel " + << channel << "\n"; + } + + double voltage = 0.0; + error = fDB.GetNamedChannelData(channel, "voltage", voltage); + if (error) { + throw cet::exception(fLogCategory) + << "Error (code " << error << ") reading 'voltage' for channel " + << channel << "\n"; + } + + fDatabaseStatus[channel].status = static_cast(statusInt); + fDatabaseStatus[channel].voltage = voltage; + + if (fVerbose) + log << "\n channel " << std::setw(3) << channel + << " status " << statusInt + << " voltage " << voltage << " V"; + + } // for channel + +} // readStatusFromDB() + + +// ----------------------------------------------------------------------------- +std::uint64_t +icarusDB::PMTChannelStatusProvider::RunToDatabaseTimestamp(unsigned int run) const +{ + // Same convention as other ICARUS DB services: + // run XXXXX -> (XXXXX + 1'000'000'000) * 1'000'000'000 [nanoseconds] + // e.g. run XXXXX -> 10000XXXXX * 1e9 + std::uint64_t timestamp = std::uint64_t(run) + 1'000'000'000ULL; + timestamp *= 1'000'000'000ULL; + + if (fVerbose) + mf::LogInfo(fLogCategory) + << "Run " << run << " status from DB timestamp " << timestamp; + + return timestamp; +} + + +// ----------------------------------------------------------------------------- +icarusDB::PMTChannelStatusProvider::ChannelSet_t +icarusDB::PMTChannelStatusProvider::getChannelsWithStatus + (PMTChannelStatusValue status) const +{ + ChannelSet_t result; + for (auto const& [channel, info] : fDatabaseStatus) + if (info.status == status) result.insert(channel); + return result; +} + + +// ----------------------------------------------------------------------------- diff --git a/icaruscode/PMT/Status/PMTChannelStatusProvider.h b/icaruscode/PMT/Status/PMTChannelStatusProvider.h new file mode 100644 index 000000000..9b917ff46 --- /dev/null +++ b/icaruscode/PMT/Status/PMTChannelStatusProvider.h @@ -0,0 +1,112 @@ +/** + * @file icaruscode/PMT/Status/PMTChannelStatusProvider.h + * @brief Provider for PMT channel status from the conditions database. + * @author Matteo Vicenzi (mvicenzi@bnl.gov) + * @see icaruscode/PMT/Status/PMTChannelStatusProvider.cxx + */ + +#ifndef ICARUSCODE_PMT_STATUS_PMTCHANNELSTATUSPROVIDER_H +#define ICARUSCODE_PMT_STATUS_PMTCHANNELSTATUSPROVIDER_H + +// ICARUS libraries +#include "icaruscode/PMT/Status/PMTChannelStatus.h" + +// LArSoft libraries +#include "larevt/CalibrationDBI/Providers/DBFolder.h" + +// Framework libraries +#include "fhiclcpp/ParameterSet.h" +#include "cetlib_except/exception.h" + +// C++ standard libraries +#include +#include +#include + + +namespace icarusDB::details { + + /// Internal structure holding the status and voltage for a single channel. + struct PMTChannelStatusDB { + PMTChannelStatusValue status = kOFF; + double voltage = 0.0; ///< Nominal voltage [V]. + }; + +} // namespace icarusDB::details + + +// ----------------------------------------------------------------------------- +namespace icarusDB { class PMTChannelStatusProvider; } +/** + * @brief Provider for PMT channel status from the conditions database. + * + * Reads the channel status (kOFF=0, kON=1, kBad=2) from the + * `pmt_voltage_data` database table, indexed by run number. + * + * The database interface is accessed only on `readStatusFromDB()` calls, + * and the relevant information is cached. + * + * Configuration parameters + * ------------------------- + * * `DBname` (default: `"pmt_voltage_data"`): database table name. + * * `StatusTag` (default: `""`): database tag to select. + * * `DefaultStatus` (default: `0` = kOFF): status for channels absent from DB. + * * `Verbose` (default: `false`): print channel statuses when loading. + * * `LogCategory` (default: `"PMTChannelStatusProvider"`). + * + */ +class icarusDB::PMTChannelStatusProvider : public PMTChannelStatus { + +public: + + /// Constructor from FHiCL configuration (no database access yet). + PMTChannelStatusProvider(const fhicl::ParameterSet& pset); + + /// Loads channel statuses for the given run from the database. + void readStatusFromDB(unsigned int run); + + /// Converts a run number to the internal 19-digit nanosecond timestamp. + std::uint64_t RunToDatabaseTimestamp(unsigned int run) const; + + // --- PMTChannelStatus interface --- + + PMTChannelStatusValue getChannelStatus(unsigned int channel) const override + { return getChannelStatusOrDefault(channel).status; } + + ChannelSet_t getOnChannels() const override { return getChannelsWithStatus(kON); } + ChannelSet_t getOffChannels() const override { return getChannelsWithStatus(kOFF); } + ChannelSet_t getBadChannels() const override { return getChannelsWithStatus(kBad); } + + double getChannelVoltage(unsigned int channel) const override + { return getChannelStatusOrDefault(channel).voltage; } + + std::string getDatabaseTag() const override { return fStatusTag; } + +private: + + using PMTChannelStatusDB = details::PMTChannelStatusDB; + + bool fVerbose; + std::string fLogCategory; + std::string fStatusTag; + + PMTChannelStatusDB fDefault; ///< Status used for channels not present in DB. + + lariov::DBFolder fDB; ///< Cached database interface. + + /// Map of channel ID to status information. + std::map fDatabaseStatus; + + /// Returns the channel status record, or the default if channel is absent. + PMTChannelStatusDB const& getChannelStatusOrDefault(unsigned int channel) const { + auto const it = fDatabaseStatus.find(channel); + return (it == fDatabaseStatus.end()) ? fDefault : it->second; + } + + /// Returns the set of all channels with the given status. + ChannelSet_t getChannelsWithStatus(PMTChannelStatusValue status) const; + +}; // class icarusDB::PMTChannelStatusProvider + + +#endif // ICARUSCODE_PMT_STATUS_PMTCHANNELSTATUSPROVIDER_H diff --git a/icaruscode/PMT/Status/PMTChannelStatusService_service.cc b/icaruscode/PMT/Status/PMTChannelStatusService_service.cc new file mode 100644 index 000000000..ba031690c --- /dev/null +++ b/icaruscode/PMT/Status/PMTChannelStatusService_service.cc @@ -0,0 +1,60 @@ +/** + * @file icaruscode/PMT/Status/PMTChannelStatusService_service.cc + * @brief Art service for PMT channel status from the conditions database. + * @author Matteo Vicenzi (mvicenzi@bnl.gov) + */ + +#include "icaruscode/PMT/Status/IPMTChannelStatusService.h" +#include "icaruscode/PMT/Status/PMTChannelStatusProvider.h" + +// Framework libraries +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Services/Registry/ActivityRegistry.h" +#include "art/Framework/Services/Registry/ServiceDefinitionMacros.h" +#include "art/Framework/Services/Registry/ServiceDeclarationMacros.h" +#include "fhiclcpp/ParameterSet.h" +#include "cetlib_except/exception.h" + + +// ----------------------------------------------------------------------------- +namespace icarusDB { class PMTChannelStatusService; } +class icarusDB::PMTChannelStatusService + : public IPMTChannelStatusService, private PMTChannelStatusProvider { + + void preBeginRun(const art::Run& run); + + /// Returns a constant pointer to the service provider. + virtual PMTChannelStatusProvider const* do_provider() const override + { return this; } + + public: + + PMTChannelStatusService(const fhicl::ParameterSet& pset, art::ActivityRegistry& reg); + +}; // class icarusDB::PMTChannelStatusService + + +// ----------------------------------------------------------------------------- +// --- Implementation +// ----------------------------------------------------------------------------- +icarusDB::PMTChannelStatusService::PMTChannelStatusService + (const fhicl::ParameterSet& pset, art::ActivityRegistry& reg) + : PMTChannelStatusProvider(pset) +{ + reg.sPreBeginRun.watch(this, &PMTChannelStatusService::preBeginRun); +} + + +// ----------------------------------------------------------------------------- +void icarusDB::PMTChannelStatusService::preBeginRun(const art::Run& run) +{ + readStatusFromDB(run.run()); +} + + +// ----------------------------------------------------------------------------- +DECLARE_ART_SERVICE_INTERFACE_IMPL(icarusDB::PMTChannelStatusService, icarusDB::IPMTChannelStatusService, SHARED) +DEFINE_ART_SERVICE_INTERFACE_IMPL(icarusDB::PMTChannelStatusService, icarusDB::IPMTChannelStatusService) + + +// ----------------------------------------------------------------------------- diff --git a/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl b/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl new file mode 100644 index 000000000..542eedf2b --- /dev/null +++ b/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl @@ -0,0 +1,18 @@ +# +# @file pmt_channel_status_icarus.fcl +# @brief Standard configuration for the ICARUS PMT channel status service. +# + +BEGIN_PROLOG + +icarus_pmt_channel_status: { + service_provider: PMTChannelStatusService + DBname: "pmt_voltage_data" + StatusTag: "v1r3" + DefaultStatus: 1 # kON + DefaultVoltage: 1500.0 # [V] + Verbose: false + LogCategory: "PMTChannelStatusProvider" +} + +END_PROLOG From ee0f6af7ca577fd89b028c3dd7365003703eddbf Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 16 Mar 2026 15:44:16 -0400 Subject: [PATCH 13/36] tie up the service --- fcl/reco/Definitions/stage0_icarus_driver_common.fcl | 2 ++ fcl/reco/Stage0/data/partial/decodePMT_icarus.fcl | 4 +++- fcl/services/services_icarus_simulation.fcl | 3 ++- 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/fcl/reco/Definitions/stage0_icarus_driver_common.fcl b/fcl/reco/Definitions/stage0_icarus_driver_common.fcl index 942cd67c1..b674ce281 100644 --- a/fcl/reco/Definitions/stage0_icarus_driver_common.fcl +++ b/fcl/reco/Definitions/stage0_icarus_driver_common.fcl @@ -4,6 +4,7 @@ #include "stage0_icarus_defs.fcl" #include "services_common_icarus.fcl" #include "channelmapping_icarus.fcl" +#include "pmt_channel_status_icarus.fcl" process_name: Stage0 @@ -13,6 +14,7 @@ services: IICARUSChannelMap: @local::icarus_channelmappinggservice IPMTTimingCorrectionService: @local::icarus_pmttimingservice IPhotonCalibratorService: @local::icarus_photon_calibration + IPMTChannelStatusService: @local::icarus_pmt_channel_status @table::icarus_wirecalibration_minimum_services } diff --git a/fcl/reco/Stage0/data/partial/decodePMT_icarus.fcl b/fcl/reco/Stage0/data/partial/decodePMT_icarus.fcl index d22be3e79..b70b737a1 100644 --- a/fcl/reco/Stage0/data/partial/decodePMT_icarus.fcl +++ b/fcl/reco/Stage0/data/partial/decodePMT_icarus.fcl @@ -53,6 +53,7 @@ #include "services_common_icarus.fcl" #include "channelmapping_icarus.fcl" #include "timing_icarus.fcl" +#include "pmt_channel_status_icarus.fcl" #include "rootoutput_icarus.fcl" #include "decoderdefs_icarus.fcl" @@ -71,7 +72,8 @@ services: { DetectorClocksService: @local::icarus_detectorclocks IICARUSChannelMap: @local::icarus_channelmappinggservice IPMTTimingCorrectionService: @local::icarus_pmttimingservice - + IPMTChannelStatusService: @local::icarus_pmt_channel_status + TFileService: { fileName: "Trees-%ifb_%tc-%p.root" } } diff --git a/fcl/services/services_icarus_simulation.fcl b/fcl/services/services_icarus_simulation.fcl index bbdc9eeee..b3605c5bd 100644 --- a/fcl/services/services_icarus_simulation.fcl +++ b/fcl/services/services_icarus_simulation.fcl @@ -40,6 +40,7 @@ #include "signalservices_icarus.fcl" #include "photpropservices_icarus.fcl" #include "timing_icarus.fcl" +#include "pmt_channel_status_icarus.fcl" BEGIN_PROLOG @@ -118,7 +119,7 @@ icarus_detsim_services: { @table::icarus_detsim_dark_services IPMTTimingCorrectionService: @local::icarus_pmttimingservice - # PmtGainService: @local::icarus_pmtgain_service + IPMTChannelStatusService: @local::icarus_pmt_channel_status } # icarus_detsim_services From be29f4160680a724fc9a9e1dd6f643db7848d821 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 16 Mar 2026 15:19:13 -0500 Subject: [PATCH 14/36] fix linking --- icaruscode/PMT/Status/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/icaruscode/PMT/Status/CMakeLists.txt b/icaruscode/PMT/Status/CMakeLists.txt index 0b628976d..73581e3de 100644 --- a/icaruscode/PMT/Status/CMakeLists.txt +++ b/icaruscode/PMT/Status/CMakeLists.txt @@ -12,7 +12,8 @@ set( LIB_LIBRARIES set( SERVICE_LIBRARIES icaruscode_PMT_Status - larcore::CoreUtils + art::Framework_Principal + larcore::Geometry_Geometry_service ) file(GLOB lib_srcs *.cxx) From a47faadd2f7eb762f2c4690544fec58ecd64d395 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 16 Mar 2026 16:51:07 -0400 Subject: [PATCH 15/36] change kOn to kGood --- icaruscode/PMT/Status/PMTChannelStatus.h | 20 +++++++++---------- .../PMT/Status/PMTChannelStatusProvider.cxx | 2 +- .../PMT/Status/PMTChannelStatusProvider.h | 10 +++++----- .../PMT/Status/pmt_channel_status_icarus.fcl | 2 +- 4 files changed, 17 insertions(+), 17 deletions(-) diff --git a/icaruscode/PMT/Status/PMTChannelStatus.h b/icaruscode/PMT/Status/PMTChannelStatus.h index 8708b4959..86241ce2d 100644 --- a/icaruscode/PMT/Status/PMTChannelStatus.h +++ b/icaruscode/PMT/Status/PMTChannelStatus.h @@ -22,8 +22,8 @@ namespace icarusDB { /// Possible status values for a PMT channel, matching the database integers. enum PMTChannelStatusValue : int { - kOFF = 0, ///< Channel is off (not powered). - kON = 1, ///< Channel is on and good. + kOff = 0, ///< Channel is off (not powered). + kGood = 1, ///< Channel is on and good. kBad = 2 ///< Channel is powered but flagged as bad. }; @@ -32,8 +32,8 @@ namespace icarusDB { * @brief Interface for PMT channel status information. * * Currently, the class provides interface for the following information: - * - channel status: kON, kOFF, or kBad - * - convenience predicates: isOn(), isOff(), isBad() + * - channel status: kGood, kOff, or kBad + * - convenience predicates: isGood(), isOff(), isBad() * - channel sets by status: getOnChannels(), getOffChannels(), getBadChannels() * */ @@ -48,22 +48,22 @@ namespace icarusDB { /// Returns the status of the specified channel. virtual PMTChannelStatusValue getChannelStatus(unsigned int channel) const = 0; - /// Returns whether the specified channel is on and good. - virtual bool isOn(unsigned int channel) const - { return getChannelStatus(channel) == kON; } + /// Returns whether the specified channel is on and good (status kGood). + virtual bool isGood(unsigned int channel) const + { return getChannelStatus(channel) == kGood; } /// Returns whether the specified channel is off (not powered). virtual bool isOff(unsigned int channel) const - { return getChannelStatus(channel) == kOFF; } + { return getChannelStatus(channel) == kOff; } /// Returns whether the specified channel is powered but bad. virtual bool isBad(unsigned int channel) const { return getChannelStatus(channel) == kBad; } - /// Returns the set of channels with status kON. + /// Returns the set of channels with status kGood. virtual ChannelSet_t getOnChannels() const = 0; - /// Returns the set of channels with status kOFF. + /// Returns the set of channels with status kOff. virtual ChannelSet_t getOffChannels() const = 0; /// Returns the set of channels with status kBad. diff --git a/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx b/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx index a46c3eea0..0c1d7a7d8 100644 --- a/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx +++ b/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx @@ -26,7 +26,7 @@ icarusDB::PMTChannelStatusProvider::PMTChannelStatusProvider , fDB ( pset.get("DBname", "pmt_voltage_data"), "", "", fStatusTag, true, false ) { - int const defaultStatusInt = pset.get("DefaultStatus", static_cast(kON)); + int const defaultStatusInt = pset.get("DefaultStatus", static_cast(kGood)); fDefault.status = static_cast(defaultStatusInt); fDefault.voltage = pset.get("DefaultVoltage", 1500.0); diff --git a/icaruscode/PMT/Status/PMTChannelStatusProvider.h b/icaruscode/PMT/Status/PMTChannelStatusProvider.h index 9b917ff46..913a1c7df 100644 --- a/icaruscode/PMT/Status/PMTChannelStatusProvider.h +++ b/icaruscode/PMT/Status/PMTChannelStatusProvider.h @@ -28,7 +28,7 @@ namespace icarusDB::details { /// Internal structure holding the status and voltage for a single channel. struct PMTChannelStatusDB { - PMTChannelStatusValue status = kOFF; + PMTChannelStatusValue status = kOff; double voltage = 0.0; ///< Nominal voltage [V]. }; @@ -40,7 +40,7 @@ namespace icarusDB { class PMTChannelStatusProvider; } /** * @brief Provider for PMT channel status from the conditions database. * - * Reads the channel status (kOFF=0, kON=1, kBad=2) from the + * Reads the channel status (kOff=0, kGood=1, kBad=2) from the * `pmt_voltage_data` database table, indexed by run number. * * The database interface is accessed only on `readStatusFromDB()` calls, @@ -50,7 +50,7 @@ namespace icarusDB { class PMTChannelStatusProvider; } * ------------------------- * * `DBname` (default: `"pmt_voltage_data"`): database table name. * * `StatusTag` (default: `""`): database tag to select. - * * `DefaultStatus` (default: `0` = kOFF): status for channels absent from DB. + * * `DefaultStatus` (default: `0` = kOff): status for channels absent from DB. * * `Verbose` (default: `false`): print channel statuses when loading. * * `LogCategory` (default: `"PMTChannelStatusProvider"`). * @@ -73,8 +73,8 @@ class icarusDB::PMTChannelStatusProvider : public PMTChannelStatus { PMTChannelStatusValue getChannelStatus(unsigned int channel) const override { return getChannelStatusOrDefault(channel).status; } - ChannelSet_t getOnChannels() const override { return getChannelsWithStatus(kON); } - ChannelSet_t getOffChannels() const override { return getChannelsWithStatus(kOFF); } + ChannelSet_t getOnChannels() const override { return getChannelsWithStatus(kGood); } + ChannelSet_t getOffChannels() const override { return getChannelsWithStatus(kOff); } ChannelSet_t getBadChannels() const override { return getChannelsWithStatus(kBad); } double getChannelVoltage(unsigned int channel) const override diff --git a/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl b/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl index 542eedf2b..b1cd61dec 100644 --- a/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl +++ b/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl @@ -9,7 +9,7 @@ icarus_pmt_channel_status: { service_provider: PMTChannelStatusService DBname: "pmt_voltage_data" StatusTag: "v1r3" - DefaultStatus: 1 # kON + DefaultStatus: 1 # kGood DefaultVoltage: 1500.0 # [V] Verbose: false LogCategory: "PMTChannelStatusProvider" From 27c3dea88b331d953324654a096b353e2e33ce73 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 16 Mar 2026 16:51:53 -0400 Subject: [PATCH 16/36] use channel status service in pmtsim --- icaruscode/PMT/CMakeLists.txt | 1 + icaruscode/PMT/SimPMTIcarus_module.cc | 19 ++++++++++++++++++- 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/icaruscode/PMT/CMakeLists.txt b/icaruscode/PMT/CMakeLists.txt index 21084fab3..0fba647c5 100644 --- a/icaruscode/PMT/CMakeLists.txt +++ b/icaruscode/PMT/CMakeLists.txt @@ -15,6 +15,7 @@ add_subdirectory(Status) cet_build_plugin(SimPMTIcarus art::module LIBRARIES icaruscode_PMT_Algorithms + icaruscode_PMT_Status lardataobj::RawData lardataobj::Simulation larcore::Geometry_Geometry_service diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 03bb95b1b..7e047f2b4 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -18,6 +18,7 @@ #include "icaruscode/IcarusObj/OpDetWaveformMeta.h" #include "icaruscode/Timing/IPMTTimingCorrectionService.h" #include "icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h" +#include "icaruscode/PMT/Status/IPMTChannelStatusService.h" // LArSoft libraries #include "larcore/CoreUtils/ServiceUtil.h" @@ -201,6 +202,12 @@ namespace icarus::opdet { fhicl::TableFragment algoConfig; + fhicl::Atom UseChannelStatusDatabase { + Name("UseChannelStatusDatabase"), + Comment("skip simulation of channels flagged in the PMT channel status DB"), + false + }; + fhicl::Atom MakeMetadata { Name("MakeMetadata"), Comment("writes a metadata object for each waveform"), @@ -274,6 +281,7 @@ namespace icarus::opdet { bool fWritePhotons { false }; ///< Whether to save contributing photons. bool fDoTimingDelays; ///< Whether timing delay corrections are applied. bool fUseGainCalibDB; ///< Whether per-channel gain calibration from DB is applied. + bool fUseChannelStatusDB; ///< Whether to skip non-ON channels using the status DB. CLHEP::HepRandomEngine& fEfficiencyEngine; CLHEP::HepRandomEngine& fDarkNoiseEngine; @@ -320,6 +328,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) , fWritePhotons(config().writePhotons()) , fDoTimingDelays(config().algoConfig().ApplyTimingDelays()) , fUseGainCalibDB(config().algoConfig().UseGainDatabase()) + , fUseChannelStatusDB(config().UseChannelStatusDatabase()) // random engines , fEfficiencyEngine(art::ServiceHandle()->registerAndSeedEngine( createEngine(0, "HepJamesRandom", "Efficiencies"), @@ -419,7 +428,11 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) if(pmtVector.isValid()) { nopch = pmtVector->size(); for(auto const& photons : *pmtVector) { - + + if (fUseChannelStatusDB + && !lar::providerFrom()->isGood(photons.OpChannel())) + continue; + // Make an empty SimPhotonsLite with the same channel number. sim::SimPhotonsLite lite_photons(photons.OpChannel()); @@ -439,6 +452,10 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) nopch = pmtLiteVector->size(); for(auto const& lite_photons : *pmtLiteVector) { + if (fUseChannelStatusDB + && !lar::providerFrom()->isGood(lite_photons.OpChannel)) + continue; + // Make an empty SimPhotons with the same channel number. sim::SimPhotons photons(lite_photons.OpChannel); From 4245aa26f73141972b86e0d0aada58e9428d4920 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 16 Mar 2026 17:09:53 -0400 Subject: [PATCH 17/36] update fcls for legacy vs overlay flow --- fcl/detsim/detsim_1d_icarus.fcl | 4 +++- fcl/detsim/detsim_2d_icarus.fcl | 4 +++- fcl/detsim/detsim_2d_icarus_refactored.fcl | 4 +++- .../detsim_2d_icarus_refactored_Run4.fcl | 4 +++- .../detsim_2d_icarus_refactored_overlay.fcl | 4 +++- ...etsim_2d_icarus_refactored_overlay_Run4.fcl | 4 +++- .../detsim_2d_icarus_refactored_yzsim.fcl | 4 +++- .../detsim_2d_icarus_refactored_yzsim_Run4.fcl | 4 +++- ...im_2d_icarus_refactored_yzsim_notrigger.fcl | 4 +++- ...arus_refactored_yzsim_notrigger_overlay.fcl | 4 +++- ...tsim_2d_icarus_refactored_yzsim_overlay.fcl | 4 +++- ...2d_icarus_refactored_yzsim_overlay_Run4.fcl | 4 +++- ...etonly_icarus_overlay_reprocessing_run2.fcl | 4 +++- ...etonly_icarus_overlay_reprocessing_run4.fcl | 4 +++- .../PMT/Algorithms/pmtsimulation_icarus.fcl | 18 ++++++++++++------ 15 files changed, 54 insertions(+), 20 deletions(-) diff --git a/fcl/detsim/detsim_1d_icarus.fcl b/fcl/detsim/detsim_1d_icarus.fcl index 9a6acf4b4..81e364794 100644 --- a/fcl/detsim/detsim_1d_icarus.fcl +++ b/fcl/detsim/detsim_1d_icarus.fcl @@ -99,4 +99,6 @@ services.Geometry.GDML: "icarus_complete_20220518_overburden.gdml" services.Geometry.ROOT: "icarus_complete_20220518_overburden.gdml" physics.producers.crtdaq.G4ModuleLabel: "largeant" physics.producers.opdaq.InputModule: "largeant" -physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.UseGainDatabase: false # keep false for non-overlay MC +physics.producers.opdaq.UseChannelStatusDatabase: false # keep false for non-overlay MC diff --git a/fcl/detsim/detsim_2d_icarus.fcl b/fcl/detsim/detsim_2d_icarus.fcl index 55ca76db5..4fe9e6a6b 100644 --- a/fcl/detsim/detsim_2d_icarus.fcl +++ b/fcl/detsim/detsim_2d_icarus.fcl @@ -55,4 +55,6 @@ outputs: { services.Geometry.GDML: "icarus_complete_20220518_overburden.gdml" services.Geometry.ROOT: "icarus_complete_20220518_overburden.gdml" physics.producers.crtdaq.G4ModuleLabel: "shifted" -physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.UseGainDatabase: false # keep false for non-overlay MC +physics.producers.opdaq.UseChannelStatusDatabase: false # keep false for non-overlay MC diff --git a/fcl/detsim/detsim_2d_icarus_refactored.fcl b/fcl/detsim/detsim_2d_icarus_refactored.fcl index 5bc3c7f6e..0e026fb70 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored.fcl @@ -53,6 +53,8 @@ outputs: { physics.producers.crtdaq.G4ModuleLabel: "shifted" physics.producers.opdaq.InputModule: "pdfastsim" -physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.UseGainDatabase: false # keep false for non-overlay MC +physics.producers.opdaq.UseChannelStatusDatabase: false # keep false for non-overlay MC physics.producers.shifted.InitAuxDetSimChannelLabel: "genericcrt" physics.producers.shifted.InitSimPhotonsLabel: "pdfastsim" diff --git a/fcl/detsim/detsim_2d_icarus_refactored_Run4.fcl b/fcl/detsim/detsim_2d_icarus_refactored_Run4.fcl index 9c2b91bb6..e1bbabd1c 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_Run4.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_Run4.fcl @@ -2,4 +2,6 @@ # Run3/4 optical tune physics.producers.opdaq: @local::icarus_simpmt_run4 -physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.UseGainDatabase: false # keep false for non-overlay MC +physics.producers.opdaq.UseChannelStatusDatabase: false # keep false for non-overlay MC diff --git a/fcl/detsim/detsim_2d_icarus_refactored_overlay.fcl b/fcl/detsim/detsim_2d_icarus_refactored_overlay.fcl index b9e8ab646..cb2b2b80a 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_overlay.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_overlay.fcl @@ -4,5 +4,7 @@ physics.producers.daq: @local::icarus_simwire_wirecell_shifted_overlay # turn off mc noise on pmt waveforms physics.producers.opdaq: @local::icarus_simpmt_nonoise -physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.UseGainDatabase: true +physics.producers.opdaq.UseChannelStatusDatabase: true diff --git a/fcl/detsim/detsim_2d_icarus_refactored_overlay_Run4.fcl b/fcl/detsim/detsim_2d_icarus_refactored_overlay_Run4.fcl index 4fa4a4d9b..26f416910 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_overlay_Run4.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_overlay_Run4.fcl @@ -2,4 +2,6 @@ # Run3/4 optical tune (with no noise) physics.producers.opdaq: @local::icarus_simpmt_run4_nonoise -physics.producers.opdaq.ApplyTimingDelays: true \ No newline at end of file +physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.UseGainDatabase: true +physics.producers.opdaq.UseChannelStatusDatabase: true \ No newline at end of file diff --git a/fcl/detsim/detsim_2d_icarus_refactored_yzsim.fcl b/fcl/detsim/detsim_2d_icarus_refactored_yzsim.fcl index 3a341343e..23ce1d8f6 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_yzsim.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_yzsim.fcl @@ -62,7 +62,9 @@ outputs: { physics.producers.crtdaq.G4ModuleLabel: "shifted" physics.producers.opdaq.InputModule: "pdfastsim" -physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.UseGainDatabase: false # keep false for non-overlay MC +physics.producers.opdaq.UseChannelStatusDatabase: false # keep false for non-overlay MC physics.producers.shifted.InitAuxDetSimChannelLabel: "genericcrt" physics.producers.shifted.InitSimPhotonsLabel: "pdfastsim" physics.producers.filtersed.InitSimEnergyDepositLabel: "shifted" diff --git a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_Run4.fcl b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_Run4.fcl index 4c3280574..8980ce32a 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_Run4.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_Run4.fcl @@ -4,4 +4,6 @@ physics.producers.daq.wcls_main.params.YZScaleMapJson: "yzmap_gain_icarus_v4_run # Run3/4 optical tune physics.producers.opdaq: @local::icarus_simpmt_run4 -physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.UseGainDatabase: false # keep false for non-overlay MC +physics.producers.opdaq.UseChannelStatusDatabase: false # keep false for non-overlay MC diff --git a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger.fcl b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger.fcl index a58a8a221..d31d057c3 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger.fcl @@ -57,4 +57,6 @@ outputs: { } physics.producers.opdaq.InputModule: "pdfastsim" -physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.ApplyTimingDelays: false # keep false for non-overlay MC +physics.producers.opdaq.UseGainDatabase: false # keep false for non-overlay MC +physics.producers.opdaq.UseChannelStatusDatabase: false # keep false for non-overlay MC diff --git a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger_overlay.fcl b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger_overlay.fcl index 639eaf030..4725fb200 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger_overlay.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_notrigger_overlay.fcl @@ -4,4 +4,6 @@ physics.producers.daq: @local::icarus_simwire_wirecell_yz_overlay # turn off mc noise on pmt waveforms physics.producers.opdaq: @local::icarus_simpmt_nonoise -physics.producers.opdaq.ApplyTimingDelays: true \ No newline at end of file +physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.UseGainDatabase: true +physics.producers.opdaq.UseChannelStatusDatabase: true \ No newline at end of file diff --git a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay.fcl b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay.fcl index f4e935a36..d621e63d7 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay.fcl @@ -4,4 +4,6 @@ physics.producers.daq: @local::icarus_simwire_wirecell_yz_overlay # turn off mc noise on pmt waveforms physics.producers.opdaq: @local::icarus_simpmt_nonoise -physics.producers.opdaq.ApplyTimingDelays: true \ No newline at end of file +physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.UseGainDatabase: true +physics.producers.opdaq.UseChannelStatusDatabase: true \ No newline at end of file diff --git a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay_Run4.fcl b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay_Run4.fcl index 152c26da6..6673bbc90 100644 --- a/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay_Run4.fcl +++ b/fcl/detsim/detsim_2d_icarus_refactored_yzsim_overlay_Run4.fcl @@ -4,4 +4,6 @@ physics.producers.daq.wcls_main.params.YZScaleMapJson: "yzmap_gain_icarus_v4_run # Run3/4 optical tune (with no noise) physics.producers.opdaq: @local::icarus_simpmt_run4_nonoise -physics.producers.opdaq.ApplyTimingDelays: true \ No newline at end of file +physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.UseGainDatabase: true +physics.producers.opdaq.UseChannelStatusDatabase: true \ No newline at end of file diff --git a/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run2.fcl b/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run2.fcl index 9611436c8..e9a395661 100644 --- a/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run2.fcl +++ b/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run2.fcl @@ -72,7 +72,9 @@ physics.producers: { # build waveforms from simPhotons (already shifted) physics.producers.opdaq.InputModule: "shifted::DetSim" -physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.UseGainDatabase: true +physics.producers.opdaq.UseChannelStatusDatabase: true # make sure the triggersim chain used the "new" products # instead of the "old" ones which I cannot drop on input diff --git a/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run4.fcl b/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run4.fcl index 5c104fcb5..4c8e0ac01 100644 --- a/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run4.fcl +++ b/fcl/detsim/partial/detsim_opdetonly_icarus_overlay_reprocessing_run4.fcl @@ -11,4 +11,6 @@ # switch to Run3/Run-4 tune physics.producers.opdaq: @local::icarus_simpmt_run4_nonoise physics.producers.opdaq.InputModule: "shifted::DetSim" -physics.producers.opdaq.ApplyTimingDelays: true \ No newline at end of file +physics.producers.opdaq.ApplyTimingDelays: true +physics.producers.opdaq.UseGainDatabase: true +physics.producers.opdaq.UseChannelStatusDatabase: true \ No newline at end of file diff --git a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl index 1957be5ca..668c251a2 100644 --- a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl +++ b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl @@ -387,6 +387,7 @@ icarus_pmtsimulationalg_202202_noise: { FluctuateGain: true # apply per-photoelectron gain fluctuations ApplyTimingDelays: false # legacy value UseGainDatabase: false # legacy value + UseChannelStatusDatabase: false # legacy value PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage @@ -450,17 +451,19 @@ icarus_pmtsimulationalg_run2_2025: { # Apply timing delays (cables, transit time, ..) to the simulated photons # These delays are extracted from the calibration database - ApplyTimingDelays: true + ApplyTimingDelays: true + # Use channel status databaset to skip dead or bad channels + UseChannelStatusDatabase: true # Threshold in ADC counts for a PMT self-trigger ThresholdADC: 15 # ADC counts - + # arranged ad hoc to get ~40% Poisson fluctuation on gain (effectively: 40.8%) PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage VoltageDistribution: [ 3.0, 3.4, 5.0, 3.33, 1.67, 1.0, 1.2, 1.5, 2.2, 3.0 ] Gain: @local::icarus_settings_opdet_run2.nominalPMTgain # total PMT gain - } + } } # icarus_pmtsimulationalg_run2_2025 @@ -481,17 +484,19 @@ icarus_pmtsimulationalg_run4_2025: { # Apply timing delays (cables, transit time, ..) to the simulated photons # These delays are extracted from the calibration database - ApplyTimingDelays: true + ApplyTimingDelays: true + # Use channel status databaset to skip dead or bad channels + UseChannelStatusDatabase: true # Threshold in ADC counts for a PMT self-trigger ThresholdADC: 15 # ADC counts - + # arranged ad hoc to get ~40% Poisson fluctuation on gain (effectively: 40.8%) PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage VoltageDistribution: [ 3.0, 3.4, 5.0, 3.33, 1.67, 1.0, 1.2, 1.5, 2.2, 3.0 ] Gain: @local::icarus_settings_opdet_run4.nominalPMTgain # total PMT gain - } + } } # icarus_pmtsimulationalg_run4_2025 @@ -536,6 +541,7 @@ icarus_pmtsimulationalg_2018: { FluctuateGain: true # apply per-photoelectron gain fluctuations ApplyTimingDelays: false # legacy value UseGainDatabase: false # legacy value + UseChannelStatusDatabase: false # legacy value PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage From 8203c273cc04f37f2402eb33770f78fd8152873f Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 17 Mar 2026 09:35:25 -0400 Subject: [PATCH 18/36] switch reco to local entry point --- .../OpReco/Algorithms/OpRecoFactoryStuff.h | 61 +++++++++++++++++-- .../PMT/OpReco/ICARUSOpHitFinder_module.cc | 32 +++++++++- .../PMT/OpReco/fcl/icarus_ophitfinder.fcl | 9 ++- 3 files changed, 91 insertions(+), 11 deletions(-) diff --git a/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h b/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h index 8a1173f76..ac5343c7f 100644 --- a/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h +++ b/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h @@ -30,6 +30,7 @@ // ----------------------------------------------------------------------------- // forward declarations namespace art { class Event; } +namespace pmtana { class RiseTimeCalculatorBase; } // ----------------------------------------------------------------------------- @@ -291,15 +292,21 @@ class opdet::factory::AlgorithmFactory { */ std::unique_ptr create (std::string const& name, fhicl::ParameterSet const& pset) const; - + + /// As `create(name, pset)`, also injecting a rise time calculator into the algorithm. + std::unique_ptr create( + std::string const& name, + fhicl::ParameterSet const& pset, + std::unique_ptr calc) const; + /** * @brief Creates an instance of the algorithm constructed with `pset`. * @param pset the configuration of the algorithm * @return the newly created algorithm object - * + * * The type of algorithm is discovered from `pset` itself with the mechanism * documented in `setAlgorithmConfigurationKey()`. - * + * * For example: * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} * factory.setAlgorithmConfigurationKey("Name"); @@ -311,6 +318,11 @@ class opdet::factory::AlgorithmFactory { * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ std::unique_ptr create(fhicl::ParameterSet const& pset) const; + + /// As `create(pset)`, also injecting a rise time calculator into the algorithm. + std::unique_ptr create( + fhicl::ParameterSet const& pset, + std::unique_ptr calc) const; /// Returns a list with the names of supported algorithms. std::vector names() const; @@ -449,7 +461,12 @@ struct opdet::factory::AlgorithmFactory::AlgoMaker { */ virtual std::unique_ptr makeAlgo (fhicl::ParameterSet const& pset) const = 0; - + + /// Creates an algorithm of the wrapped type, injecting a rise time calculator. + virtual std::unique_ptr makeAlgo( + fhicl::ParameterSet const& pset, + std::unique_ptr calc) const = 0; + /// Comparison: sort by name in lexicographic order. bool operator< (AlgoMaker const& other) const { return name < other.name; } @@ -468,7 +485,12 @@ struct opdet::factory::AlgorithmFactory::AlgoMakerFor std::unique_ptr makeAlgo (fhicl::ParameterSet const& pset) const override { return std::make_unique(pset); } - + + std::unique_ptr makeAlgo( + fhicl::ParameterSet const& pset, + std::unique_ptr calc) const override + { return std::make_unique(pset, std::move(calc)); } + }; // opdet::factory::AlgorithmFactory::AlgoMakerFor @@ -534,6 +556,35 @@ auto opdet::factory::AlgorithmFactory::create } // opdet::factory::AlgorithmFactory::create() +// ----------------------------------------------------------------------------- +template +auto opdet::factory::AlgorithmFactory::create( + std::string const& name, + fhicl::ParameterSet const& pset, + std::unique_ptr calc) const + -> std::unique_ptr +{ + AlgoMaker const* maker = getMaker(name); + if (maker) return maker->makeAlgo(pset, std::move(calc)); + + throw cet::exception{ "AlgorithmFactory" } + << "Unknown algorithm: '" << name << "'.\nSupported algorithms:\n - '" + << names("'\n - '") << "'\n"; + +} // opdet::factory::AlgorithmFactory::create(name, pset, calc) + + +// ----------------------------------------------------------------------------- +template +auto opdet::factory::AlgorithmFactory::create( + fhicl::ParameterSet const& pset, + std::unique_ptr calc) const + -> std::unique_ptr +{ + return create(pset.get(fAlgoNameKey), pset, std::move(calc)); +} // opdet::factory::AlgorithmFactory::create(pset, calc) + + // ----------------------------------------------------------------------------- template std::vector opdet::factory::AlgorithmFactory::names() const { diff --git a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc index f9ae98c29..db5f15067 100644 --- a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc +++ b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc @@ -13,6 +13,7 @@ #include "icaruscode/PMT/Data/WaveformRMS.h" #include "sbnobj/ICARUS/PMT/Data/WaveformBaseline.h" // LArSoft libraries +#include "larana/OpticalDetector/OpHitFinder/RiseTimeTools/RiseTimeCalculatorBase.h" #include "larana/OpticalDetector/OpHitFinder/AlgoCFD.h" #include "larana/OpticalDetector/OpHitFinder/AlgoFixedWindow.h" #include "larana/OpticalDetector/OpHitFinder/AlgoSiPM.h" @@ -45,9 +46,11 @@ #include "messagefacility/MessageLogger/MessageLogger.h" #include "fhiclcpp/types/Sequence.h" #include "fhiclcpp/types/OptionalAtom.h" +#include "fhiclcpp/types/OptionalDelegatedParameter.h" #include "fhiclcpp/types/Atom.h" #include "fhiclcpp/types/DelegatedParameter.h" #include "fhiclcpp/ParameterSet.h" +#include "art/Utilities/make_tool.h" // C++ standard libraries #include // std::copy_if(), std::binary_search() @@ -244,13 +247,21 @@ class opdet::ICARUSOpHitFinder: public art::ReplicatedProducer { fhicl::DelegatedParameter PedAlgoPset { Name{ "PedAlgoPset" }, - Comment{ + Comment{ "parameters of the pedestal extraction algorithm." " The parameters must include the algorithm \"Name\", one of: " + PedAlgoFactory.names(", ") + "." } }; - + + fhicl::OptionalDelegatedParameter RiseTimeCalculator { + Name{ "RiseTimeCalculator" }, + Comment{ + "Configuration of the rise time calculator art tool (optional)." + " If omitted, rise time is not computed and the OpHit RiseTime field is left at zero." + } + }; + fhicl::Atom UseCalibrator { Name{ "UseCalibrator" }, Comment{ "Use the photon calibration service configured in the job" }, @@ -325,6 +336,10 @@ class opdet::ICARUSOpHitFinder: public art::ReplicatedProducer { std::vector selectWaveforms (std::vector const& waveforms) const; + /// Creates a rise time calculator tool from config, or returns nullptr if not configured. + static std::unique_ptr + makeRiseTimeCalculator(Config const& config); + }; // opdet::ICARUSOpHitFinder @@ -494,7 +509,8 @@ opdet::ICARUSOpHitFinder::ICARUSOpHitFinder // algorithms , fPulseRecoMgr{} , fThreshAlg - { HitAlgoFactory.create(params().HitAlgoPset.get()) } + { HitAlgoFactory.create(params().HitAlgoPset.get(), + makeRiseTimeCalculator(params())) } , fPedAlg { PedAlgoFactory.create(params().PedAlgoPset.get()) } { @@ -666,6 +682,16 @@ std::vector opdet::ICARUSOpHitFinder::selectWaveforms } // selectWaveforms() +// ----------------------------------------------------------------------------- +std::unique_ptr +opdet::ICARUSOpHitFinder::makeRiseTimeCalculator(Config const& config) +{ + fhicl::ParameterSet calcPset; + if (!config.RiseTimeCalculator.get_if_present(calcPset)) return nullptr; + return art::make_tool(calcPset); +} // opdet::ICARUSOpHitFinder::makeRiseTimeCalculator() + + // ----------------------------------------------------------------------------- DEFINE_ART_MODULE(opdet::ICARUSOpHitFinder) diff --git a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl index 650351356..7f3a2f7e1 100644 --- a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl +++ b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl @@ -186,10 +186,13 @@ icarus_opreco_hit_cfd: { icarus_ophit: # some basic configuration { - module_type: "OpHitFinder" + ## in the past we were using directly module: "OpHitFinder" from larana + ## we are now using our own entrypoint to the reconstruction "ICARUSOpHitFinder" + ## this allows interfacing with our PMT channel status service + module_type: "ICARUSOpHitFinder" GenModule: "generator" InputModule: "opdaq" - InputLabels: [] + # InputLabels: [] #not supported by ICARUSOpHitFinder ChannelMasks: [] UseCalibrator: true # use calibrator service (SPEArea database) # if true, overrides SPEArea / SPEShift below @@ -200,7 +203,7 @@ icarus_ophit: # some basic configuration SPEShift: @local::SPE.Shift # Used by PhotonCalibratorStandard to compute # PE = area / SPEArea + SPEShift UseStartTime: false # start and peak times each in its own data member - reco_man: @local::standard_preco_manager + # reco_man: @local::standard_preco_manager #now managed internally by ICARUSOpHitFinder HitAlgoPset: @local::icarus_opreco_hit_slidingwindow_201910 PedAlgoPset: @local::icarus_opreco_pedestal_DocDB24969 RiseTimeCalculator: { From ca85cc39e263fd4cb2db3c0bf10275fd66e025f0 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 17 Mar 2026 09:52:27 -0400 Subject: [PATCH 19/36] fix sharing factory with ped algos --- icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h b/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h index ac5343c7f..053481960 100644 --- a/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h +++ b/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h @@ -489,7 +489,14 @@ struct opdet::factory::AlgorithmFactory::AlgoMakerFor std::unique_ptr makeAlgo( fhicl::ParameterSet const& pset, std::unique_ptr calc) const override - { return std::make_unique(pset, std::move(calc)); } + { + if constexpr (std::is_constructible_v>) + return std::make_unique(pset, std::move(calc)); + else + return std::make_unique(pset); // algorithm does not use a rise time calculator + } }; // opdet::factory::AlgorithmFactory::AlgoMakerFor From f457297bf06706a05a193112c3d20662c6535a06 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 17 Mar 2026 16:16:44 -0400 Subject: [PATCH 20/36] fix standard calibrator + add channel status service --- .../PMT/OpReco/ICARUSOpHitFinder_module.cc | 77 ++++++++++++------- .../PMT/OpReco/fcl/icarus_ophitfinder.fcl | 11 ++- 2 files changed, 56 insertions(+), 32 deletions(-) diff --git a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc index db5f15067..660496fc5 100644 --- a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc +++ b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc @@ -8,6 +8,7 @@ */ // ICARUS libraries +#include "icaruscode/PMT/Status/IPMTChannelStatusService.h" #include "icaruscode/PMT/OpReco/Algorithms/PedAlgoFixed.h" #include "icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h" #include "icaruscode/PMT/Data/WaveformRMS.h" @@ -256,10 +257,7 @@ class opdet::ICARUSOpHitFinder: public art::ReplicatedProducer { fhicl::OptionalDelegatedParameter RiseTimeCalculator { Name{ "RiseTimeCalculator" }, - Comment{ - "Configuration of the rise time calculator art tool (optional)." - " If omitted, rise time is not computed and the OpHit RiseTime field is left at zero." - } + Comment{ "Configuration of the rise time calculator art tool (optional)."} }; fhicl::Atom UseCalibrator { @@ -285,7 +283,13 @@ class opdet::ICARUSOpHitFinder: public art::ReplicatedProducer { Comment{ "shift on the single photoelectron response" }, [this](){ return !UseCalibrator(); } }; - + + fhicl::Atom UseChannelStatus { + Name{ "UseChannelStatusDatabase" }, + Comment{"ignore waveforms on channels flagged in PMT channel status database."}, + false + }; + }; // Config using Parameters = art::ReplicatedProducer::Table; @@ -309,12 +313,15 @@ class opdet::ICARUSOpHitFinder: public art::ReplicatedProducer { // --- BEGIN -- Cached service values ---------------------------------------- /// Storage for our own calibration algorithm, if any. std::unique_ptr fMyCalib; - + unsigned int const fMaxOpChannel; ///< Number of channels in the detector. - + /// The calibration algorithm to be used (not owned). calib::IPhotonCalibrator const* fCalib = nullptr; - + + /// Optional PMT channel status provider (not owned); null if not in use. + icarusDB::PMTChannelStatus const* fChannelStatus = nullptr; + // --- END ---- Cached service values ---------------------------------------- @@ -527,14 +534,15 @@ opdet::ICARUSOpHitFinder::ICARUSOpHitFinder // if (params().UseCalibrator()) { fCalib = frame.serviceHandle()->provider(); + mf::LogInfo{ "ICARUSOpHitFinder" } << "Using 'calib::IPhotonCalibratorService' for gain calibration."; } else { - fMyCalib = std::make_unique( - params().AreaToPE() - , params().SPEArea() - , params().SPEShift() - ); + fMyCalib = std::make_unique(params().SPEArea(), params().SPEShift(), params().AreaToPE()); fCalib = fMyCalib.get(); + mf::LogInfo{ "ICARUSOpHitFinder" } << "Using 'AreaToPE' = " << params().AreaToPE() + << ", 'SPEArea' = " << params().SPEArea() + << ", 'SPEShift' = " << params().SPEShift() + << " for gain calibration."; } // if ... else // @@ -543,8 +551,16 @@ opdet::ICARUSOpHitFinder::ICARUSOpHitFinder fPulseRecoMgr.AddRecoAlgo(fThreshAlg.get()); fPulseRecoMgr.SetDefaultPedAlgo(&(fPedAlg->algo())); + // + // channel status service + // + if (params().UseChannelStatus()){ + fChannelStatus = frame.serviceHandle()->provider(); + mf::LogInfo{ "ICARUSOpHitFinder" } << "Using 'icarusDB::IPMTChannelStatusService' for channel status."; + } + // framework hooks to the algorithms - if (!fChannelMasks.empty() && (fPedAlg->algo().Name() == "Fixed")) { + if ((!fChannelMasks.empty() || fChannelStatus) && (fPedAlg->algo().Name() == "Fixed")) { /* TODO * The reason why this is not working is complicate. * The interface of ophit mini-framework is quite minimal and tight, @@ -566,7 +582,7 @@ opdet::ICARUSOpHitFinder::ICARUSOpHitFinder */ throw art::Exception{ art::errors::Configuration } << "Pedestal algorithm \"Fixed\" will not work" - " since a channel mask list is specified.\n"; + " since a channel mask list or channel status database is specified.\n"; } fPedAlg->initialize(consumesCollector()); @@ -593,17 +609,18 @@ void opdet::ICARUSOpHitFinder::produce // MaybeOwnerPtr const> waveforms // = fChannelMasks.empty()? &allWaveforms: selectWaveforms(allWaveforms); using WaveformsPtr_t = MaybeOwnerPtr const>; - WaveformsPtr_t waveforms = fChannelMasks.empty() - ? WaveformsPtr_t{ &allWaveforms } - : WaveformsPtr_t{ selectWaveforms(allWaveforms) } + bool const needsSelection = !fChannelMasks.empty() || fChannelStatus; + WaveformsPtr_t waveforms = needsSelection + ? WaveformsPtr_t{ selectWaveforms(allWaveforms) } + : WaveformsPtr_t{ &allWaveforms } ; - - if (fChannelMasks.empty() && (&allWaveforms != &*waveforms)) { + + if (!needsSelection && (&allWaveforms != &*waveforms)) { // if this happens, contact the author throw art::Exception{ art::errors::LogicError } << "Bug in MaybeOwnerPtr!\n"; } - assert(fChannelMasks.empty() || (&allWaveforms == &*waveforms)); + assert(needsSelection || (&allWaveforms == &*waveforms)); // // run the algorithm @@ -668,16 +685,20 @@ std::vector opdet::ICARUSOpHitFinder::selectWaveforms (std::vector const& waveforms) const { std::vector selected; - - auto const isNotMasked = [this](raw::OpDetWaveform const& waveform) + + auto const isGood = [this](raw::OpDetWaveform const& waveform) { - return !std::binary_search - (fChannelMasks.begin(), fChannelMasks.end(), waveform.ChannelNumber()); + unsigned int const ch = waveform.ChannelNumber(); + if (std::binary_search(fChannelMasks.begin(), fChannelMasks.end(), ch)) + return false; + if (fChannelStatus && !fChannelStatus->isGood(ch)) + return false; + return true; }; - + std::copy_if - (waveforms.begin(), waveforms.end(), back_inserter(selected), isNotMasked); - + (waveforms.begin(), waveforms.end(), back_inserter(selected), isGood); + return selected; } // selectWaveforms() diff --git a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl index 7f3a2f7e1..726a66c3b 100644 --- a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl +++ b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl @@ -193,8 +193,9 @@ icarus_ophit: # some basic configuration GenModule: "generator" InputModule: "opdaq" # InputLabels: [] #not supported by ICARUSOpHitFinder - ChannelMasks: [] - UseCalibrator: true # use calibrator service (SPEArea database) + ChannelMasks: [] + UseChannelStatusDatabase: true # skip channels flagged in the status database + UseCalibrator: true # use calibrator service (SPEArea database) # if true, overrides SPEArea / SPEShift below HitThreshold: 0.2 # PE AreaToPE: true # Use area to calculate number of PEs @@ -207,8 +208,9 @@ icarus_ophit: # some basic configuration HitAlgoPset: @local::icarus_opreco_hit_slidingwindow_201910 PedAlgoPset: @local::icarus_opreco_pedestal_DocDB24969 RiseTimeCalculator: { - tool_type: RiseTimeThreshold - PeakRatio: 0.15 # at 15% of the peak amplitude -- not tuned + tool_type: RiseTimeThreshold + PeakRatio: 0.15 # at 15% of the peak amplitude -- not tuned + InterpolateSamples: true # sub-sample linear interpolation } } @@ -221,6 +223,7 @@ icarus_ophitdebugger.OutputFile: "ophit_debug.root" icarus_ophit_MC: { @table::icarus_ophit UseCalibrator: false + UseChannelStatusDatabase: false InputModule: "opdaq" } From 617b74c9931c3b8d3369870cdda8a9fe20fee774 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 17 Mar 2026 17:17:14 -0400 Subject: [PATCH 21/36] remove conflicting defaults --- icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc index 660496fc5..94f63c16d 100644 --- a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc +++ b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc @@ -269,19 +269,16 @@ class opdet::ICARUSOpHitFinder: public art::ReplicatedProducer { fhicl::Atom AreaToPE { Name{ "AreaToPE" }, Comment{ "Whether the `SPEArea` parameter refers to area or amplitude" }, - [this](){ return !UseCalibrator(); } }; fhicl::Atom SPEArea { Name{ "SPEArea" }, Comment{ "area or amplitude of PMT response to single photoelectron" }, - [this](){ return !UseCalibrator(); } }; fhicl::Atom SPEShift { Name{ "SPEShift" }, Comment{ "shift on the single photoelectron response" }, - [this](){ return !UseCalibrator(); } }; fhicl::Atom UseChannelStatus { From 66fda19c09f3c6a006a08f4c06ffed9230f12e94 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 17 Mar 2026 16:18:15 -0500 Subject: [PATCH 22/36] remove mcophit from mc-overlays --- .../Definitions/stage0_icarus_mc_defs.fcl | 1 - .../overlay/stage0_run3_wc_icarus_overlay.fcl | 7 ------ .../stage0_run3_wcdnn_icarus_overlay.fcl | 7 ------ .../overlay/stage0_run4_wc_icarus_overlay.fcl | 7 ------ .../stage0_run4_wcdnn_icarus_overlay.fcl | 7 ------ ..._opdetonly_icarus_overlay_reprocessing.fcl | 9 ------- ..._opdetonly_icarus_overlay_reprocessing.fcl | 12 --------- ..._opdetonly_icarus_overlay_reprocessing.fcl | 12 --------- .../PMT/OpReco/fcl/icarus_opana_modules.fcl | 4 +-- .../PMT/OpReco/fcl/icarus_ophitfinder.fcl | 25 ++++++++++++------- 10 files changed, 18 insertions(+), 73 deletions(-) delete mode 100644 fcl/reco/Stage0/overlay/stage0_run3_wc_icarus_overlay.fcl delete mode 100644 fcl/reco/Stage0/overlay/stage0_run3_wcdnn_icarus_overlay.fcl delete mode 100644 fcl/reco/Stage0/overlay/stage0_run4_wc_icarus_overlay.fcl delete mode 100644 fcl/reco/Stage0/overlay/stage0_run4_wcdnn_icarus_overlay.fcl delete mode 100644 fcl/reco/Stage0/partial/stage0_run3_opdetonly_icarus_overlay_reprocessing.fcl delete mode 100644 fcl/reco/Stage0/partial/stage0_run4_opdetonly_icarus_overlay_reprocessing.fcl diff --git a/fcl/reco/Definitions/stage0_icarus_mc_defs.fcl b/fcl/reco/Definitions/stage0_icarus_mc_defs.fcl index 22cc2e7ea..dfe65af81 100644 --- a/fcl/reco/Definitions/stage0_icarus_mc_defs.fcl +++ b/fcl/reco/Definitions/stage0_icarus_mc_defs.fcl @@ -76,7 +76,6 @@ icarus_stage0_overlay_PMT: [ ophituncorrected, # @local::icarus_ophit_data ophit, # ovveride to @local::icarus_ophit_timing_correction # this happens in enable_overlay.fcl - mcophit, opflashCryoE, opflashCryoW ] diff --git a/fcl/reco/Stage0/overlay/stage0_run3_wc_icarus_overlay.fcl b/fcl/reco/Stage0/overlay/stage0_run3_wc_icarus_overlay.fcl deleted file mode 100644 index 605db37ef..000000000 --- a/fcl/reco/Stage0/overlay/stage0_run3_wc_icarus_overlay.fcl +++ /dev/null @@ -1,7 +0,0 @@ -#include "stage0_run2_wc_icarus_overlay.fcl" - -# set RUN-3 tune for optical reconstruction -physics.producers.ophit.SPEArea: @local::SPERun3.Area -physics.producers.ophit.SPEShift: @local::SPERun3.Shift -physics.producers.mcophit.SPEArea: @local::SPERun3.Area -physics.producers.mcophit.SPEAmplitude: @local::SPERun3.Amplitude diff --git a/fcl/reco/Stage0/overlay/stage0_run3_wcdnn_icarus_overlay.fcl b/fcl/reco/Stage0/overlay/stage0_run3_wcdnn_icarus_overlay.fcl deleted file mode 100644 index 60a56c23b..000000000 --- a/fcl/reco/Stage0/overlay/stage0_run3_wcdnn_icarus_overlay.fcl +++ /dev/null @@ -1,7 +0,0 @@ -#include "stage0_run2_wcdnn_icarus_overlay.fcl" - -# set RUN-3 tune for optical reconstruction -physics.producers.ophit.SPEArea: @local::SPERun3.Area -physics.producers.ophit.SPEShift: @local::SPERun3.Shift -physics.producers.mcophit.SPEArea: @local::SPERun3.Area -physics.producers.mcophit.SPEAmplitude: @local::SPERun3.Amplitude diff --git a/fcl/reco/Stage0/overlay/stage0_run4_wc_icarus_overlay.fcl b/fcl/reco/Stage0/overlay/stage0_run4_wc_icarus_overlay.fcl deleted file mode 100644 index 79b8d7ac0..000000000 --- a/fcl/reco/Stage0/overlay/stage0_run4_wc_icarus_overlay.fcl +++ /dev/null @@ -1,7 +0,0 @@ -#include "stage0_run2_wc_icarus_overlay.fcl" - -# set RUN-4 tune for optical reconstruction -physics.producers.ophit.SPEArea: @local::SPERun4.Area -physics.producers.ophit.SPEShift: @local::SPERun4.Shift -physics.producers.mcophit.SPEArea: @local::SPERun4.Area -physics.producers.mcophit.SPEAmplitude: @local::SPERun4.Amplitude diff --git a/fcl/reco/Stage0/overlay/stage0_run4_wcdnn_icarus_overlay.fcl b/fcl/reco/Stage0/overlay/stage0_run4_wcdnn_icarus_overlay.fcl deleted file mode 100644 index b1dfd535c..000000000 --- a/fcl/reco/Stage0/overlay/stage0_run4_wcdnn_icarus_overlay.fcl +++ /dev/null @@ -1,7 +0,0 @@ -#include "stage0_run2_wcdnn_icarus_overlay.fcl" - -# set RUN-4 tune for optical reconstruction -physics.producers.ophit.SPEArea: @local::SPERun4.Area -physics.producers.ophit.SPEShift: @local::SPERun4.Shift -physics.producers.mcophit.SPEArea: @local::SPERun4.Area -physics.producers.mcophit.SPEAmplitude: @local::SPERun4.Amplitude diff --git a/fcl/reco/Stage0/partial/stage0_run2_opdetonly_icarus_overlay_reprocessing.fcl b/fcl/reco/Stage0/partial/stage0_run2_opdetonly_icarus_overlay_reprocessing.fcl index eafb4a00d..c8eb75138 100644 --- a/fcl/reco/Stage0/partial/stage0_run2_opdetonly_icarus_overlay_reprocessing.fcl +++ b/fcl/reco/Stage0/partial/stage0_run2_opdetonly_icarus_overlay_reprocessing.fcl @@ -16,10 +16,6 @@ ## This would make the data component consistent and the MC signal biased. ## Such a choice makes the configuration independent of the run period... -## ... excpet for the `mcophit` module (fake hit reconstruction) -## This module requires a fixed SPE tune (no db interface) so it must be manually set -## to the run-period SPE tune used in DetSim. - process_name: OpMCstage0 ## only perform PMT/CRT MCstage0 on top of reprocessed MC overlay files @@ -33,10 +29,5 @@ physics.path: [ @sequence::icarus_stage0_overlay_PMT, @sequence::icarus_stage0_mc_crtreco # only reco, hits from overlay ] -## update mcophit producer to pick up Run-2 SPE tune -## (ophit picks up SPEArea from the database automatically) -physics.producers.mcophit.SPEArea: @local::SPERun2.Area -physics.producers.mcophit.SPEAmplitude: @local::SPERun2.Amplitude - diff --git a/fcl/reco/Stage0/partial/stage0_run3_opdetonly_icarus_overlay_reprocessing.fcl b/fcl/reco/Stage0/partial/stage0_run3_opdetonly_icarus_overlay_reprocessing.fcl deleted file mode 100644 index c4fc39427..000000000 --- a/fcl/reco/Stage0/partial/stage0_run3_opdetonly_icarus_overlay_reprocessing.fcl +++ /dev/null @@ -1,12 +0,0 @@ -## File: stage0_run3_opdetonly_icarus_overlay_reprocessing.fcl -## Purpose: Reprocessing of the optical products for the 2025 SBN "spring" MC OVERLAY production - -#include "stage0_run2_opdetonly_icarus_overlay_reprocessing.fcl" - -## This job selectively re-runs only the PMT/CRT reconstruction path. -## This is valid only for the post-processing of the 2025 SBN MC OVERLAY productions. - -## update mcophit to pick up Run-3 SPE tune; see Run-2 base file for the rationale -## (ophit picks up SPEArea from the database automatically) -physics.producers.mcophit.SPEArea: @local::SPERun3.Area -physics.producers.mcophit.SPEAmplitude: @local::SPERun3.Amplitude diff --git a/fcl/reco/Stage0/partial/stage0_run4_opdetonly_icarus_overlay_reprocessing.fcl b/fcl/reco/Stage0/partial/stage0_run4_opdetonly_icarus_overlay_reprocessing.fcl deleted file mode 100644 index c5d7c4fa2..000000000 --- a/fcl/reco/Stage0/partial/stage0_run4_opdetonly_icarus_overlay_reprocessing.fcl +++ /dev/null @@ -1,12 +0,0 @@ -## File: stage0_run4_opdetonly_icarus_overlay_reprocessing.fcl -## Purpose: Reprocessing of the optical products for the 2025 SBN "spring" MC OVERLAY production - -#include "stage0_run2_opdetonly_icarus_overlay_reprocessing.fcl" - -## This job selectively re-runs only the PMT/CRT reconstruction path. -## This is valid only for the post-processing of the 2025 SBN MC OVERLAY productions. - -## update mcophit to pick up Run-4 SPE tune; see Run-2 base file for the rationale -## (ophit picks up SPEArea from the database automatically) -physics.producers.mcophit.SPEArea: @local::SPERun4.Area -physics.producers.mcophit.SPEAmplitude: @local::SPERun4.Amplitude diff --git a/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl b/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl index 628ef5c5d..6d39330ca 100644 --- a/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl +++ b/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl @@ -25,8 +25,8 @@ ICARUSMCOpHit: { module_type: "ICARUSMCOpHit" MergePeriod: 0.01 SimPhotonsProducer: "pdfastsim" - SPEArea: @local::SPE.Area - SPEAmplitude: @local::SPE.Amplitude + SPEArea: @local::SPERun2.Area + SPEAmplitude: @local::SPERun2.Amplitude } ICARUSMCOpFlash: { diff --git a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl index 726a66c3b..d01da3070 100644 --- a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl +++ b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl @@ -186,31 +186,38 @@ icarus_opreco_hit_cfd: { icarus_ophit: # some basic configuration { - ## in the past we were using directly module: "OpHitFinder" from larana - ## we are now using our own entrypoint to the reconstruction "ICARUSOpHitFinder" + ## in the past we were using directly module "OpHitFinder" from larana + ## we are now using our own entry point to the reconstruction "ICARUSOpHitFinder" ## this allows interfacing with our PMT channel status service module_type: "ICARUSOpHitFinder" GenModule: "generator" InputModule: "opdaq" - # InputLabels: [] #not supported by ICARUSOpHitFinder + #InputLabels: [] # not supported by ICARUSOpHitFinder + ChannelMasks: [] UseChannelStatusDatabase: true # skip channels flagged in the status database + UseCalibrator: true # use calibrator service (SPEArea database) # if true, overrides SPEArea / SPEShift below - HitThreshold: 0.2 # PE AreaToPE: true # Use area to calculate number of PEs - SPEArea: @local::SPE.Area # If AreaToPE is true, this number is - # used as single PE area (in ADC counts) - SPEShift: @local::SPE.Shift # Used by PhotonCalibratorStandard to compute + SPEArea: @local::SPERun2.Area # If AreaToPE is true, this number is + # used as single PE area (in ADC counts) + SPEShift: @local::SPERun2.Shift # Used by PhotonCalibratorStandard to compute # PE = area / SPEArea + SPEShift + UseStartTime: false # start and peak times each in its own data member - # reco_man: @local::standard_preco_manager #now managed internally by ICARUSOpHitFinder + + #reco_man: @local::standard_preco_manager # now managed internally by ICARUSOpHitFinder + + HitThreshold: 0.2 # PE HitAlgoPset: @local::icarus_opreco_hit_slidingwindow_201910 + PedAlgoPset: @local::icarus_opreco_pedestal_DocDB24969 + RiseTimeCalculator: { tool_type: RiseTimeThreshold PeakRatio: 0.15 # at 15% of the peak amplitude -- not tuned - InterpolateSamples: true # sub-sample linear interpolation + #InterpolateSamples: true # sub-sample linear interpolation } } From a6a4bcf68030c4faf8a9a8e352f77970e5b60687 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Tue, 17 Mar 2026 22:19:16 -0400 Subject: [PATCH 23/36] never hold service pointers, you fool --- icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx | 10 +++++----- icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 10 +++++----- .../Calibration/ICARUSPhotonCalibratorServiceFromDB.h | 3 ++- icaruscode/PMT/SimPMTIcarus_module.cc | 2 +- 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index c03e6ffe1..3d7a7b98f 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -963,7 +963,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, - calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -974,7 +974,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( { return std::make_unique(makeParams( beamGateTimestamp, - larProp, clockData, timingDelays, gainCalibService, + larProp, clockData, timingDelays, gainCalibProvider, SPRfunction, pedestalGenerator, mainRandomEngine, darkNoiseRandomEngine, elecNoiseRandomEngine, trackSelectedPhotons @@ -989,7 +989,7 @@ auto icarus::opdet::PMTsimulationAlgMaker::makeParams( detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, - calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -1013,10 +1013,10 @@ auto icarus::opdet::PMTsimulationAlgMaker::makeParams( params.clockData = &clockData; params.detTimings = detinfo::makeDetectorTimings(params.clockData); params.timingDelays = timingDelays; - params.gainCalibService = gainCalibService; + params.gainCalibProvider = gainCalibProvider; assert(!params.doTimingDelays || params.timingDelays); - assert(!params.useGainCalibDB || params.gainCalibService); + assert(!params.useGainCalibDB || params.gainCalibProvider); params.pulseFunction = &SPRfunction; params.pedestalGen = &pedestalGenerator; diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index 8e9bac83f..023c326c6 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -22,7 +22,7 @@ #include "icaruscode/Utilities/quantities_utils.h" // util::value_t #include "icarusalg/Utilities/SampledFunction.h" #include "icaruscode/Timing/PMTTimingCorrections.h" -#include "icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h" +#include "icaruscode/PMT/Calibration/PhotonCalibratorFromDB.h" // LArSoft libraries #include "lardataobj/RawData/OpDetWaveform.h" @@ -518,8 +518,8 @@ class icarus::opdet::PMTsimulationAlg { /// Timing delays service interfacing with database icarusDB::PMTTimingCorrections const* timingDelays = nullptr; - /// SPE calibration service for per-channel gain (nullptr = disabled). - calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService = nullptr; + /// SPE calibration provider for per-channel gain (nullptr = disabled). + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider = nullptr; // detTimings is not really "optional" but it needs delayed construction. /// Detector clocks data wrapper. @@ -1128,7 +1128,7 @@ class icarus::opdet::PMTsimulationAlgMaker { detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& detClocks, icarusDB::PMTTimingCorrections const* timingDelays, - calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -1161,7 +1161,7 @@ class icarus::opdet::PMTsimulationAlgMaker { detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, - calib::ICARUSPhotonCalibratorServiceFromDB const* gainCalibService, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, diff --git a/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h b/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h index c4ef076be..993f2155d 100644 --- a/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h +++ b/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h @@ -62,9 +62,10 @@ namespace calib { double getSPEFitWidth(int channel) const { return fProvider.getSPEFitWidth(channel); } - private: provider_type const* provider() const override { return &fProvider; } + private: + void preBeginRun(art::Run const& run); icarusDB::PhotonCalibratorFromDB fProvider; diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 7e047f2b4..038b5307c 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -406,7 +406,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) *(lar::providerFrom()), clockData, fDoTimingDelays ? lar::providerFrom() : nullptr, - fUseGainCalibDB ? &(*art::ServiceHandle()) : nullptr, + fUseGainCalibDB ? art::ServiceHandle()->provider() : nullptr, *fSinglePhotonResponseFunc, *fPedestalGen, fEfficiencyEngine, From bd7dd28d8eaa0e4abb8afb2ccf7303748fc2edd0 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 10:06:00 -0400 Subject: [PATCH 24/36] expose integral and bias from sampled waveform tool --- .../Algorithms/PhotoelectronPulseFunction.h | 24 ++++++++++++ .../PMT/Algorithms/SampledWaveformFunction.h | 37 +++++++++++-------- .../PMT/Algorithms/pmtsimulation_icarus.fcl | 28 +++++++------- .../PMT/SampledWaveformFunctionTool_tool.cc | 12 +++++- 4 files changed, 72 insertions(+), 29 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PhotoelectronPulseFunction.h b/icaruscode/PMT/Algorithms/PhotoelectronPulseFunction.h index 78c6f9036..2f3139477 100644 --- a/icaruscode/PMT/Algorithms/PhotoelectronPulseFunction.h +++ b/icaruscode/PMT/Algorithms/PhotoelectronPulseFunction.h @@ -75,6 +75,24 @@ class icarus::opdet::PhotoelectronPulseFunction { /// Returns the polarity of the pulse (`+1`: positive, or `-1`: negative). int polarity() const { return doPolarity(); } + /** + * @brief Returns the integral of the discretised pulse [ADC × tick]. + * + * For sampled implementations this is the exact sum of all stored samples. + * The default implementation returns zero; concrete classes that support + * integration should override this method. + */ + ADCcount integral() const { return doIntegral(); } + + /** + * @brief Returns a bias constant associated with this pulse. + * + * The physical interpretation of this constant is defined by the concrete + * implementation and the FHiCL configuration. Returns `1.0` if the + * implementation does not define one. + */ + double biasConstant() const { return doBiasConstant(); } + // @{ /** @@ -137,6 +155,12 @@ class icarus::opdet::PhotoelectronPulseFunction { virtual int doPolarity() const { return ((peakAmplitude() - baseline()) >= ADCcount{ 0 })? +1: -1; } + /// Returns the integral of the discretised pulse [ADC × tick]. + virtual ADCcount doIntegral() const { return ADCcount{ 0 }; } + + /// Returns the bias constant associated with this pulse (1.0 if not set). + virtual double doBiasConstant() const { return 1.0; } + /** * @brief Prints into the stream the parameters of this shape. diff --git a/icaruscode/PMT/Algorithms/SampledWaveformFunction.h b/icaruscode/PMT/Algorithms/SampledWaveformFunction.h index a232c4e96..c1675c013 100644 --- a/icaruscode/PMT/Algorithms/SampledWaveformFunction.h +++ b/icaruscode/PMT/Algorithms/SampledWaveformFunction.h @@ -63,7 +63,7 @@ class icarus::opdet::SampledWaveformFunction struct WaveformSpecs_t { std::string name { "" }; ///< Name of this waveform (short). std::string description; ///< Description of this waveform. - std::string date { "n/a" }; ///< An indication of the date of the waveform. + std::string date { "n/a" }; ///< An indication of the date of the waveform. unsigned int version { 1 }; ///< Version number. std::vector samples; ///< Samples [mV] Time sampleDuration; ///< Time extension of one sample. @@ -86,7 +86,9 @@ class icarus::opdet::SampledWaveformFunction * The polarity of the waveform is deduced by the value at peak. * */ - SampledWaveformFunction(WaveformSpecs_t specs, Time peakTime, float gain); + SampledWaveformFunction( + WaveformSpecs_t specs, Time peakTime, float gain, + double biasRatio = 1.0); /// @{ /// @name Parameter accessors. @@ -94,6 +96,12 @@ class icarus::opdet::SampledWaveformFunction /// Returns the gain the waveform is representing. float gain() const { return fGain; } + /// Returns the integral of the waveform [ADC × tick]. + ADCcount doIntegral() const override; + + /// Returns the bias constant as set in the FHiCL configuration. + double doBiasConstant() const override { return fBiasRatio; } + /// @} private: @@ -103,6 +111,8 @@ class icarus::opdet::SampledWaveformFunction std::vector const fSamples; ///< All samples. float const fGain; ///< The gain this waveform represents. + + double const fBiasRatio; ///< User-configurable bias ratio from FHiCL. std::size_t const fPeakSample; ///< The sample with the absolute peak. @@ -146,11 +156,7 @@ class icarus::opdet::SampledWaveformFunction /// Returns whether a sample with the specified index is within the range. bool hasSample(std::ptrdiff_t index) const { return (index >= 0) && (std::size_t(index) < fSamples.size()); } - - - /// Returns the integral of the waveform. - ADCcount integral() const; - + /** * @brief Transforms the input waveform. * @param waveform input waveform (in millivolt and for a known gain) @@ -174,12 +180,13 @@ class icarus::opdet::SampledWaveformFunction // ----------------------------------------------------------------------------- template icarus::opdet::SampledWaveformFunction::SampledWaveformFunction - (WaveformSpecs_t waveformSpecs, Time peakTime, float gain) - : fSource { std::move(waveformSpecs) } - , fSamples { buildSamples(gain) } - , fGain { gain } - , fPeakSample{ findPeak(fSamples) } - , fRefTime { peakTime - fPeakSample * sampleDuration() } + (WaveformSpecs_t waveformSpecs, Time peakTime, float gain, double biasRatio) + : fSource { std::move(waveformSpecs) } + , fSamples { buildSamples(gain) } + , fGain { gain } + , fBiasRatio { biasRatio } + , fPeakSample { findPeak(fSamples) } + , fRefTime { peakTime - fPeakSample * sampleDuration() } {} @@ -212,7 +219,7 @@ void icarus::opdet::SampledWaveformFunction::doDump( << " with amplitude " << Base_t::peakAmplitude() << "\n" << indent << " start at " << fRefTime << ", gain " << fGain - << " (integral: " << integral() << ")" + << " (integral: " << doIntegral() << ")" << '\n' ; @@ -221,7 +228,7 @@ void icarus::opdet::SampledWaveformFunction::doDump( // ----------------------------------------------------------------------------- template -auto icarus::opdet::SampledWaveformFunction::integral() const -> ADCcount +auto icarus::opdet::SampledWaveformFunction::doIntegral() const -> ADCcount { return std::reduce(fSamples.begin(), fSamples.end()); } diff --git a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl index 668c251a2..d2cec5b5a 100644 --- a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl +++ b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl @@ -206,28 +206,30 @@ icarus_photoelectronresponse_202202: { # This is valid for Run1/2 cables. # icarus_photoelectronresponse_run2_2025: { - + tool_type: SampledWaveformFunctionTool - - WaveformData: "Responses/SPRisolHitRun9271.txt" - TransitTime: "55.1 ns" - Gain: @local::icarus_settings_opdet_run2.nominalPMTgain - + + WaveformData: "Responses/SPRisolHitRun9271.txt" + TransitTime: "55.1 ns" + Gain: @local::icarus_settings_opdet_run2.nominalPMTgain + BiasRatio: 1.5332 + } # icarus_photoelectronresponse_run2_2025 # # Single photoelectron response from single-PE data. # This is valid for Run3/4/5+ cables. -# +# icarus_photoelectronresponse_run4_2025: { - + tool_type: SampledWaveformFunctionTool - - WaveformData: "Responses/SPRisolHitRun12715.txt" - TransitTime: "55.1 ns" - Gain: @local::icarus_settings_opdet_run4.nominalPMTgain - + + WaveformData: "Responses/SPRisolHitRun12715.txt" + TransitTime: "55.1 ns" + Gain: @local::icarus_settings_opdet_run4.nominalPMTgain + BiasRatio: 1.2196 + } # icarus_photoelectronresponse_run4_2025 diff --git a/icaruscode/PMT/SampledWaveformFunctionTool_tool.cc b/icaruscode/PMT/SampledWaveformFunctionTool_tool.cc index 8f2059ed7..6c0f3a049 100644 --- a/icaruscode/PMT/SampledWaveformFunctionTool_tool.cc +++ b/icaruscode/PMT/SampledWaveformFunctionTool_tool.cc @@ -93,6 +93,8 @@ namespace icarus::opdet { struct SampledWaveformFunctionTool; } * gain of the response, and that response will be rescaled from its gain * value to the one specified with this parameter. If not specified, * the response is used as is. + * * `BiasRatio` (real, optional): a user-configurable ratio associated with + * the waveform gain; it affects the bias of the response (default: 1.0). * */ struct icarus::opdet::SampledWaveformFunctionTool @@ -119,7 +121,12 @@ struct icarus::opdet::SampledWaveformFunctionTool Name("Gain"), Comment("PMT amplification gain factor") }; - + + fhicl::OptionalAtom BiasRatio { + Name("BiasRatio"), + Comment("user-configurable bias of the response gain (default: 1.0)") + }; + }; // struct Config @@ -351,10 +358,13 @@ auto icarus::opdet::SampledWaveformFunctionTool::makePulseFunction // // create the algorithm // + double const biasRatio = config.BiasRatio().value_or(1.0); + return std::make_unique( std::move(waveformSpecs) // waveformSpecs , config.TransitTime() // peakTime , reqGain // gain + , biasRatio // bias ); } // icarus::opdet::SampledWaveformFunctionTool::makePulseFunction() From 2192b3b224c865b0504fab6bdb8c464fac9b63a8 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 10:06:24 -0400 Subject: [PATCH 25/36] plug in db gain and spread in simulation --- .../PMT/Algorithms/PMTsimulationAlg.cxx | 67 +++++++++++++------ icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 51 +++++++++----- 2 files changed, 80 insertions(+), 38 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 3d7a7b98f..663622fe6 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -24,6 +24,7 @@ // CLHEP libraries #include "CLHEP/Random/RandFlat.h" +#include "CLHEP/Random/RandGaussQ.h" #include "CLHEP/Random/RandPoisson.h" #include "CLHEP/Random/RandExponential.h" @@ -115,6 +116,8 @@ icarus::opdet::PMTsimulationAlg::PMTsimulationAlg 1.0e-4_ADCf, // stop sampling when signal is closer to baseline than this... 20_ns // ... for at least this long ) + , fNominalSPEArea(std::abs(fParams.pulseFunction->integral().value())) + , fBiasRatio(fParams.pulseFunction->biasConstant()) , fPedestalGen(fParams.pedestalGen) , fDiscrAlgo(selectDiscriminationAlgo(fParams.discrimAlgo)) { @@ -173,16 +176,48 @@ std::tuple, std::optional> //------------------------------------------------------------------------------ -auto icarus::opdet::PMTsimulationAlg::makeGainFluctuator() const { +icarus::opdet::PMTsimulationAlg::GainFluctuator::GainFluctuator + (double const refGain, CLHEP::RandPoisson rng) +{ + fFluctuate = [refGain, rng = std::move(rng)](double n) mutable -> double { + return rng.fire(n * refGain) / refGain; + }; +} // GainFluctuator::GainFluctuator(Poisson) + + +//------------------------------------------------------------------------------ +icarus::opdet::PMTsimulationAlg::GainFluctuator::GainFluctuator + (double const gainRatio, double const relSigma, CLHEP::RandGaussQ rng) +{ + fFluctuate = [gainRatio, relSigma, rng = std::move(rng)](double n) mutable -> double { + double const result = rng.fire(n * gainRatio, std::sqrt(n) * relSigma); + return (result > 0.0)? result: 0.0; + }; +} // GainFluctuator::GainFluctuator(Gaussian) - using Fluctuator_t = GainFluctuator; - if (fParams.doGainFluctuations) { - double const refGain = fParams.PMTspecs.firstStageGain(); - return Fluctuator_t - { refGain, CLHEP::RandPoisson{ *fParams.gainRandomEngine, refGain } }; +//------------------------------------------------------------------------------ +auto icarus::opdet::PMTsimulationAlg::makeGainFluctuator(int channel) const + -> GainFluctuator +{ + if (!fParams.doGainFluctuations) return GainFluctuator{}; // identity + + if (fParams.useGainCalibDB && fParams.gainCalibProvider) { + // DB Gaussian: per-channel SPE area and width from database + // fNominalSPEArea is the integral of the SPR template, i.e. the mean area per PE + double const speArea = fParams.gainCalibProvider->getSPEArea(channel); + double const speFitWidth = fParams.gainCalibProvider->getSPEFitWidth(channel); + // gainRatio = speArea / fNominalSPEArea: mean effective PEs per true PE + // relSigma = speFitWidth / fNominalSPEArea: sigma per sqrt(PE) + double const gainRatio = (fNominalSPEArea > 0.0) ? (speArea / fNominalSPEArea): 1.0; + double const relSigma = (fNominalSPEArea > 0.0) ? (speFitWidth / fNominalSPEArea): 0.0; + + return GainFluctuator{gainRatio, relSigma, CLHEP::RandGaussQ{ *fParams.gainRandomEngine, gainRatio, relSigma }}; } - else return Fluctuator_t{}; // default-constructed does not fluctuate anything + + // Legacy Poisson mode: fluctuate around first-dynode gain + double const refGain = fParams.PMTspecs.firstStageGain(); + return GainFluctuator{refGain, CLHEP::RandPoisson{ *fParams.gainRandomEngine, refGain }}; } // icarus::opdet::PMTsimulationAlg::makeGainFluctuator() @@ -306,7 +341,7 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform unsigned int nTotalPE [[maybe_unused]] = 0U; // unused if not in `debug` mode double nTotalEffectivePE [[maybe_unused]] = 0U; // unused if not in `debug` mode - auto gainFluctuation = makeGainFluctuator(); + auto gainFluctuation = makeGainFluctuator(channel); // go though all subsamples (starting each at a fraction of a tick) for (auto const& [ iSubsample, peMap ]: util::enumerate(peMaps)) { @@ -342,7 +377,7 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform // start=std::chrono::high_resolution_clock::now(); AddPedestal(channel, waveformStartTS, waveform); - if(fParams.darkNoiseRate > 0.0_Hz) AddDarkNoise(waveform); + if(fParams.darkNoiseRate > 0.0_Hz) AddDarkNoise(waveform, channel); // end=std::chrono::high_resolution_clock::now(); diff = end-start; // std::cout << "\tadded noise... " << channel << " " << diff.count() << std::endl; @@ -670,7 +705,7 @@ void icarus::opdet::PMTsimulationAlg::AddPedestal // ----------------------------------------------------------------------------- -void icarus::opdet::PMTsimulationAlg::AddDarkNoise(Waveform_t& wave) const { +void icarus::opdet::PMTsimulationAlg::AddDarkNoise(Waveform_t& wave, int channel) const { /* * We assume leakage current ("dark noise") is completely stochastic and * distributed uniformly in time with a fixed and known rate. @@ -701,7 +736,7 @@ void icarus::opdet::PMTsimulationAlg::AddDarkNoise(Waveform_t& wave) const { TimeToTickAndSubtickConverter const toTickAndSubtick(wsp.nSubsamples()); - auto gainFluctuation = makeGainFluctuator(); + auto gainFluctuation = makeGainFluctuator(channel); MF_LOG_TRACE("PMTsimulationAlg") << "Adding dark noise (" << fParams.darkNoiseRate << ") up to " << maxTime; @@ -845,16 +880,6 @@ auto icarus::opdet::PMTsimulationAlg::TimeToTickAndSubtickConverter::operator() } // icarus::opdet::PMTsimulationAlg::TimeToTickAndSubtickConverter::operator() -// ----------------------------------------------------------------------------- -template -double icarus::opdet::PMTsimulationAlg::GainFluctuator::operator() - (double const n) -{ - return fRandomGain - ? (fRandomGain->fire(n * fReferenceGain) / fReferenceGain): n; -} - - // ----------------------------------------------------------------------------- // --- icarus::opdet::PMTsimulationAlgMaker // ----------------------------------------------------------------------------- diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index 023c326c6..b8e8b140d 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -47,8 +47,11 @@ // CLHEP libraries #include "CLHEP/Random/RandEngine.h" // CLHEP::HepRandomEngine +#include "CLHEP/Random/RandGaussQ.h" +#include "CLHEP/Random/RandPoisson.h" // C++ standard library +#include // std::function #include #include #include @@ -255,8 +258,8 @@ class icarus::opdet::OpDetWaveformMakerClass { * the `calib::ICARUSPhotonCalibratorServiceFromDB` service. * The number of effective photoelectrons added per true photoelectron is * drawn from a Gaussian whose mean is the ratio of the channel SPE area - * to the nominal SPR template area @f$ q_0 @f$, and whose standard - * deviation is the ratio of the channel SPE width to @f$ q_0 @f$. + * to the nominal SPR template area, and whose standard + * deviation is the ratio of the channel SPE width to that nominal area. * This model preserves the mean deposited charge per photoelectron to * match the measured per-channel gain, while the nominal gain set in * `PMTspecs` and in `SinglePhotonResponse` cancels out and can be left @@ -634,28 +637,39 @@ class icarus::opdet::PMTsimulationAlg { }; // TimeToTickAndSubtickConverter - /// Applies a random gain fluctuation to the specified number of - /// photoelectrons. - template + /** + * @brief Applies a random gain fluctuation to a number of photoelectrons. + * + * Two modes are supported, selected at construction time: + * - **Legacy Poisson** (doGainFluctuations = true, useGainCalibDB = false): + * models first-dynode statistics via Poisson(n*refGain)/refGain.. + * - **DB Gaussian** (doGainFluctuations = true, useGainCalibDB = true): + * draws from Normal(n*gainRatio, sqrt(n)*relSigma) using per-channel DB values. + * - **Identity** (default-constructed): returns n unchanged. + */ class GainFluctuator { - std::optional fRandomGain; ///< Random gain extractor (optional). - double const fReferenceGain = 0.0; ///< Reference (average) gain. + /// The fluctuation callable; empty means identity (no fluctuation). + std::function fFluctuate; public: + /// Identity fluctuator: returns `n` unchanged. GainFluctuator() = default; - GainFluctuator(double const refGain, Rand&& randomGain) - : fRandomGain(std::move(randomGain)) - , fReferenceGain(refGain) - {} - /// Returns the new number of photoelectrons after fluctuation from `n`. - double operator() (double const n); + /// Legacy Poisson mode: Poisson(n*refGain) / refGain. + GainFluctuator(double refGain, CLHEP::RandPoisson rng); + + /// DB Gaussian mode: Normal(n*gainRatio, sqrt(n)*relSigma), clamped positive. + GainFluctuator(double gainRatio, double relSigma, CLHEP::RandGaussQ rng); + + /// Returns the fluctuated number of effective photoelectrons from n. + double operator()(double n) const + { return fFluctuate? fFluctuate(n): n; } }; // GainFluctuator - /// Returns a configured gain fluctuator object. - auto makeGainFluctuator() const; + /// Returns a gain fluctuator configured for the given optical channel. + GainFluctuator makeGainFluctuator(int channel) const; // --- END -- Helper functors ------------------------------------------------ @@ -666,7 +680,10 @@ class icarus::opdet::PMTsimulationAlg { megahertz fSampling; ///< Wave sampling frequency [MHz]. std::size_t fNsamples; ///< Samples per waveform. - DiscretePhotoelectronPulse wsp; /// Single photon pulse (sampled). + DiscretePhotoelectronPulse wsp; ///< Single photon pulse (sampled). + + double fNominalSPEArea = 0.0; ///< Integral of the SPR template q₀ [ADC×tick]. + double fBiasConstant = 1.0; ///< Bias constant from the FHiCL configuration. /// Pedestal and electronics noise generator algorithm. PedestalGenerator_t* fPedestalGen = nullptr; @@ -804,7 +821,7 @@ class icarus::opdet::PMTsimulationAlg { (raw::Channel_t channel, std::uint64_t time, Waveform_t& wave) const; // Add "dark" noise to baseline. - void AddDarkNoise(Waveform_t& wave) const; + void AddDarkNoise(Waveform_t& wave, int channel) const; /** From 46f9ef68959ada5d5a54f1370aa18a778cca2719 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 09:32:03 -0500 Subject: [PATCH 26/36] fix variable name --- icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 663622fe6..2040951ee 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -117,7 +117,7 @@ icarus::opdet::PMTsimulationAlg::PMTsimulationAlg 20_ns // ... for at least this long ) , fNominalSPEArea(std::abs(fParams.pulseFunction->integral().value())) - , fBiasRatio(fParams.pulseFunction->biasConstant()) + , fBiasConstant(fParams.pulseFunction->biasConstant()) , fPedestalGen(fParams.pedestalGen) , fDiscrAlgo(selectDiscriminationAlgo(fParams.discrimAlgo)) { From 206283ccafea13aa6b3863b583afcfa75da50c8e Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 10:48:08 -0400 Subject: [PATCH 27/36] keep database use false for default non-overlay MC --- .../PMT/Algorithms/pmtsimulation_icarus.fcl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl index d2cec5b5a..c407af071 100644 --- a/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl +++ b/icaruscode/PMT/Algorithms/pmtsimulation_icarus.fcl @@ -453,9 +453,11 @@ icarus_pmtsimulationalg_run2_2025: { # Apply timing delays (cables, transit time, ..) to the simulated photons # These delays are extracted from the calibration database - ApplyTimingDelays: true - # Use channel status databaset to skip dead or bad channels - UseChannelStatusDatabase: true + ApplyTimingDelays: false + # Use channel status database to skip dead or bad channels + UseChannelStatusDatabase: false + # Use per-channel gain from database + UseGainDatabase: false # Threshold in ADC counts for a PMT self-trigger ThresholdADC: 15 # ADC counts @@ -486,9 +488,11 @@ icarus_pmtsimulationalg_run4_2025: { # Apply timing delays (cables, transit time, ..) to the simulated photons # These delays are extracted from the calibration database - ApplyTimingDelays: true - # Use channel status databaset to skip dead or bad channels - UseChannelStatusDatabase: true + ApplyTimingDelays: false + # Use channel status database to skip dead or bad channels + UseChannelStatusDatabase: false + # Use per-channel gain from database + UseGainDatabase: false # Threshold in ADC counts for a PMT self-trigger ThresholdADC: 15 # ADC counts From 09cb33cf81e9e783078d232d6d536f9525abcd80 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 11:12:14 -0400 Subject: [PATCH 28/36] add forgotten service --- fcl/services/services_icarus_simulation.fcl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/fcl/services/services_icarus_simulation.fcl b/fcl/services/services_icarus_simulation.fcl index b3605c5bd..51d9e9d60 100644 --- a/fcl/services/services_icarus_simulation.fcl +++ b/fcl/services/services_icarus_simulation.fcl @@ -41,6 +41,7 @@ #include "photpropservices_icarus.fcl" #include "timing_icarus.fcl" #include "pmt_channel_status_icarus.fcl" +#include "pmt_calibration_icarus.fcl" BEGIN_PROLOG @@ -116,11 +117,12 @@ icarus_g4_services: { # Define icarus_detsim_services ... (2*) icarus_detsim_services: { - + @table::icarus_detsim_dark_services IPMTTimingCorrectionService: @local::icarus_pmttimingservice IPMTChannelStatusService: @local::icarus_pmt_channel_status - + IPhotonCalibratorService: @local::icarus_photon_calibration + } # icarus_detsim_services From cb07257b781035825e17138d06d4da6d316841d6 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 11:22:33 -0400 Subject: [PATCH 29/36] improve prints --- icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 6 +++++- icaruscode/PMT/SimPMTIcarus_module.cc | 1 + icaruscode/PMT/Status/PMTChannelStatusProvider.cxx | 4 +--- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index b8e8b140d..a407fbf44 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -1262,7 +1262,9 @@ void icarus::opdet::PMTsimulationAlg::printConfiguration << '\n' << indent << "Saturation: " << fParams.saturation << " p.e." << '\n' << indent << "doGainFluctuations: " << std::boolalpha << fParams.doGainFluctuations - << '\n' << indent << "doTimingDelays: " + << '\n' << indent << "useGainCalibDB: " + << std::boolalpha << fParams.useGainCalibDB + << '\n' << indent << "doTimingDelays: " << std::boolalpha << fParams.doTimingDelays << '\n' << indent << "PulsePolarity: " << ((fParams.pulsePolarity == 1)? "positive": "negative") << " (=" << fParams.pulsePolarity << ")" << '\n' << indent << "Sampling: " << fSampling; @@ -1271,6 +1273,8 @@ void icarus::opdet::PMTsimulationAlg::printConfiguration out << '\n' << indent << "Samples/waveform: " << fNsamples << " ticks" << '\n' << indent << "Gain at first stage: " << fParams.PMTspecs.firstStageGain() + << '\n' << indent << "SPR nominal area: " << fNominalSPEArea << " ADC x tick" + << '\n' << indent << "Bias ratio: " << fBiasConstant ; out << '\n' << indent << "Pedestal: " << fPedestalGen->toString(indent + " ", ""); diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 038b5307c..f7bb55ac9 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -419,6 +419,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) mf::LogInfo log { "SimPMTIcarus" }; log << "PMT simulation configuration (first event):\n"; PMTsimulator->printConfiguration(log); + log << "\nuseChannelStatusDB: " << std::boolalpha << fUseChannelStatusDB; } // if first time // diff --git a/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx b/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx index 0c1d7a7d8..1ecc539e0 100644 --- a/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx +++ b/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx @@ -33,9 +33,7 @@ icarusDB::PMTChannelStatusProvider::PMTChannelStatusProvider mf::LogInfo(fLogCategory) << "PMTChannelStatusProvider connected to '" << pset.get("DBname", "pmt_voltage_data") - << "' DB, tag '" << fStatusTag << "'" - << ", default status: " << defaultStatusInt - << ", default voltage: " << fDefault.voltage << " V"; + << "' DB, tag '" << fStatusTag << "'."; } From 37ffe21527b698af9021feb62606f10e40e04088 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 11:21:30 -0500 Subject: [PATCH 30/36] better printing --- icaruscode/PMT/Algorithms/PMTsimulationAlg.h | 10 ++++++---- icaruscode/PMT/SimPMTIcarus_module.cc | 2 +- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h index a407fbf44..3814b2304 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.h +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.h @@ -1251,7 +1251,7 @@ void icarus::opdet::PMTsimulationAlg::printConfiguration using namespace util::quantities::electronics_literals; out - << indent << "ADC bits: " << fParams.ADCbits + << indent << "ADC bits: " << fParams.ADCbits << " (" << fParams.ADCrange().first << " -- " << fParams.ADCrange().second << ")" << '\n' << indent << "ReadoutWindowSize: " << fParams.readoutWindowSize << " ticks" @@ -1266,7 +1266,9 @@ void icarus::opdet::PMTsimulationAlg::printConfiguration << std::boolalpha << fParams.useGainCalibDB << '\n' << indent << "doTimingDelays: " << std::boolalpha << fParams.doTimingDelays - << '\n' << indent << "PulsePolarity: " << ((fParams.pulsePolarity == 1)? "positive": "negative") << " (=" << fParams.pulsePolarity << ")" + << '\n' << indent << "PulsePolarity: " + << ((fParams.pulsePolarity == 1)? "positive": "negative") + << " (=" << fParams.pulsePolarity << ")" << '\n' << indent << "Sampling: " << fSampling; if (fParams.pulseSubsamples > 1U) out << " (subsampling: x" << fParams.pulseSubsamples << ")"; @@ -1281,14 +1283,14 @@ void icarus::opdet::PMTsimulationAlg::printConfiguration if (fParams.createBeamGateTriggers) { out << '\n' << indent << "Create " << fParams.beamGateTriggerNReps - << " beam gate triggers, one every " << fParams.beamGateTriggerRepPeriod << "."; + << " beam gate triggers, one every " << fParams.beamGateTriggerRepPeriod << "."; } else out << '\n' << indent << "Do not create beam gate triggers."; out << '\n' << indent << "... and more."; out << '\n' << indent << "Template photoelectron waveform settings:" - << '\n'; + << '\n'; wsp.dump(std::forward(out), indent + " "); out << '\n' << indent << "Track used photons: " diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index f7bb55ac9..08ed7a8af 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -418,8 +418,8 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) if (firstTime()) { mf::LogInfo log { "SimPMTIcarus" }; log << "PMT simulation configuration (first event):\n"; + log << "useChannelStatusDB: " << std::boolalpha << fUseChannelStatusDB << "\n"; PMTsimulator->printConfiguration(log); - log << "\nuseChannelStatusDB: " << std::boolalpha << fUseChannelStatusDB; } // if first time // From def9fd7ae2de83cf53b20d53d820e673f1165c3f Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 11:53:18 -0500 Subject: [PATCH 31/36] apply integral definition bias --- icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 2040951ee..3b46c3904 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -205,13 +205,14 @@ auto icarus::opdet::PMTsimulationAlg::makeGainFluctuator(int channel) const if (fParams.useGainCalibDB && fParams.gainCalibProvider) { // DB Gaussian: per-channel SPE area and width from database // fNominalSPEArea is the integral of the SPR template, i.e. the mean area per PE + // fBiasConstant covers the bias in the integral definitions btw SPR template and official reco double const speArea = fParams.gainCalibProvider->getSPEArea(channel); double const speFitWidth = fParams.gainCalibProvider->getSPEFitWidth(channel); // gainRatio = speArea / fNominalSPEArea: mean effective PEs per true PE // relSigma = speFitWidth / fNominalSPEArea: sigma per sqrt(PE) - double const gainRatio = (fNominalSPEArea > 0.0) ? (speArea / fNominalSPEArea): 1.0; + double const gainRatio = (fNominalSPEArea > 0.0) ? (speArea * fBiasConstant / fNominalSPEArea): 1.0; double const relSigma = (fNominalSPEArea > 0.0) ? (speFitWidth / fNominalSPEArea): 0.0; - + return GainFluctuator{gainRatio, relSigma, CLHEP::RandGaussQ{ *fParams.gainRandomEngine, gainRatio, relSigma }}; } From b4298cf75db35c1a5f7bb5e427c918867c858a5a Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 18 Mar 2026 14:24:31 -0500 Subject: [PATCH 32/36] update db tag --- icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl b/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl index b2b4a1476..f34699f25 100644 --- a/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl +++ b/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl @@ -48,7 +48,7 @@ pmt_bkgphotons_calibration: icarus_photon_calibration: { service_provider: ICARUSPhotonCalibratorServiceFromDB - AreaTag: "v1r1" + AreaTag: "v1r2" Defaults: { SPEArea: @local::SPERun2.Area # ADC x tick, Run2 default } From 72f6fe661f458846f642abfd6ab7e831b617b58f Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 23 Mar 2026 11:52:58 -0400 Subject: [PATCH 33/36] update variable names --- .../PMT/Algorithms/PMTsimulationAlg.cxx | 5 +- icaruscode/PMT/SimPMTIcarus_module.cc | 85 +++++++++++-------- 2 files changed, 51 insertions(+), 39 deletions(-) diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index 273fdde2d..facb22b30 100644 --- a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +++ b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx @@ -392,7 +392,8 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform subsample, waveform, startTick, static_cast(nEffectivePE) ); - // std::cout << "Channel=" << channel << ", subsample=" << iSubsample << ", tick=" << startTick << ", nPE=" << nPE << ", ePE=" << nEffectivePE << std::endl; + //std::cout << "Channel=" << channel << ", subsample=" << iSubsample << ", tick=" << startTick + // << ", nPE=" << nPE << ", ePE=" << nEffectivePE << std::endl; if (debug) { DebugPEDeposit dep; @@ -440,7 +441,7 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform if(debug){ for(std::size_t i=0; iwaveform.push_back(waveform[i].value()); + debug->waveform.push_back(waveform[i].value()); } } diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 86092ea02..e5cab132b 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -324,32 +324,41 @@ namespace icarus::opdet { bool firstTime() { return !fNotFirstTime.test_and_set(); } bool fDebugTree { false }; ///< Whether to enable debug tree output. - TTree* fTree { nullptr }; ///< Debug tree pointer. - - // debug tree variables - int run, event, channel; - float timeDelay_us; - float triggerOffsetPMT_us; - int nSamples; - int nSubsamples; - // info per waveform - std::vector wf; - // info per photon - int nPhotons; - std::vector photon_simTime_ns; - std::vector photon_trigTime_us; - std::vector photon_start_x; - std::vector photon_start_y; - std::vector photon_start_z; - std::vector photon_tick; - std::vector photon_subtick; - // info per PE deposit - int nPEDeposits; - std::vector pedeposit_tick; - std::vector pedeposit_subtick; - std::vector pedeposit_nPE; - std::vector pedeposit_nEffectivePE; - std::vector pedeposit_gainFactor; + TTree* fTree { nullptr }; ///< Debug tree pointer (owned by TFileService). + + // debug tree branch variables (one entry per optical channel per event) + // event identification + int run; ///< Run number. + int event; ///< Event number. + int channel; ///< Optical channel number. + // per-channel timing parameters + float timeDelay_us; ///< Per-channel timing delay applied [us]: laser+cosmics DB corrections. + float triggerOffsetPMT_us; ///< Global PMT readout start offset relative to trigger [us] + // waveform geometry + int nSamples; ///< number of ticks. + int nSubsamples; ///< number of subsamples (sub-tick interpolation steps) + + // full-window waveform (before zero-suppression splitting) + std::vector waveform_adc; ///< full waveform ADC values [ADC counts]. + + // per-photon info + int nDetectedPhotons; ///< Number of photons accepted after QE cut. + int nInputPhotons; ///< Total number of photons in sim::SimPhotons. + std::vector photon_simTime_ns; ///< G4 simulation time of each detected photon [ns]. + std::vector photon_waveformTime_us; ///< Photon time in waveform coordinates [us]: toTriggerTime - triggerOffsetPMT + timeDelay. + std::vector photon_start_x; ///< X position of the photon emission point [cm]. + std::vector photon_start_y; ///< Y position of the photon emission point [cm]. + std::vector photon_start_z; ///< Z position of the photon emission point [cm]. + std::vector photon_tick; ///< Waveform tick index where the photon was placed. + std::vector photon_subtick; ///< Sub-tick index (subsample bin) where the photon was placed. + + // per-PE-deposit info + int nPEDeposits; ///< Number of distinct PE deposit entries. + std::vector pedeposit_tick; ///< Waveform tick of the PE deposit. + std::vector pedeposit_subtick; ///< Sub-tick index of the PE deposit. + std::vector pedeposit_nPE; ///< Integer number of photoelectrons at this deposit (after QE, before gain fluctuation). + std::vector pedeposit_nEffectivePE; ///< Effective PE count after gain fluctuation (what was actually added to the waveform). + std::vector pedeposit_gainFactor; ///< Ratio nEffectivePE/nPE: encodes the per-deposit gain fluctuation (DB-calibrated or Poisson). }; // class SimPMTIcarus @@ -419,18 +428,19 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) fTree = tfs->make("SimPMTDebug","Debug tree for SimPMTIcarus"); fTree->Branch("run", &run, "run/I"); fTree->Branch("event", &event, "event/I"); - fTree->Branch("channel", &channel, "channel/I"); + fTree->Branch("channel", &channel, "channel/I"); fTree->Branch("nSamples", &nSamples, "nSamples/I"); fTree->Branch("nSubsamples", &nSubsamples, "nSubsamples/I"); fTree->Branch("timeDelay_us", &timeDelay_us, "timeDelay_us/F"); fTree->Branch("triggerOffsetPMT_us", &triggerOffsetPMT_us, "triggerOffsetPMT_us/F"); - fTree->Branch("wf", &wf); - fTree->Branch("nPhotons", &nPhotons, "nPhotons/I"); + fTree->Branch("waveform_adc", &waveform_adc); + fTree->Branch("nInputPhotons", &nInputPhotons, "nInputPhotons/I"); + fTree->Branch("nDetectedPhotons", &nDetectedPhotons, "nDetectedPhotons/I"); fTree->Branch("photon_start_x", &photon_start_x); fTree->Branch("photon_start_y", &photon_start_y); fTree->Branch("photon_start_z", &photon_start_z); fTree->Branch("photon_simTime_ns", &photon_simTime_ns); - fTree->Branch("photon_trigTime_us", &photon_trigTime_us); + fTree->Branch("photon_waveformTime_us", &photon_waveformTime_us); fTree->Branch("photon_tick", &photon_tick); fTree->Branch("photon_subtick", &photon_subtick); fTree->Branch("nPEDeposits", &nPEDeposits, "nPEDeposits/I"); @@ -519,19 +529,20 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) if(fDebugTree && debug && debug->photons.size()>0) { // fill debug tree variables - channel = debug->opChannel; run = e.id().run(); event = e.id().event(); + channel = debug->opChannel; timeDelay_us = debug->timeDelay_us; triggerOffsetPMT_us = debug->triggerOffsetPMT_us; nSamples = debug->nSamples; nSubsamples = debug->nSubsamples; - //fill per waveform info - wf = debug->waveform; + // fill full-window waveform + waveform_adc = debug->waveform; // fill per photon info - nPhotons = debug->photons.size(); + nInputPhotons = photons.size(); + nDetectedPhotons = debug->photons.size(); photon_simTime_ns.clear(); - photon_trigTime_us.clear(); + photon_waveformTime_us.clear(); photon_tick.clear(); photon_subtick.clear(); photon_start_x.clear(); @@ -542,7 +553,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) photon_start_y.push_back(phot.startY); photon_start_z.push_back(phot.startZ); photon_simTime_ns.push_back(phot.simTime_ns); - photon_trigTime_us.push_back(phot.trigTime_us); + photon_waveformTime_us.push_back(phot.trigTime_us); photon_tick.push_back(phot.tick); photon_subtick.push_back(phot.subtick); } @@ -560,7 +571,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) pedeposit_nEffectivePE.push_back(pedep.nEffectivePE); pedeposit_gainFactor.push_back(pedep.gainFactor()); } - + fTree->Fill(); } // if debug tree } // for From bda70244d73afa6b18bfde3773cbba7821bb15d3 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 25 Mar 2026 11:44:15 -0400 Subject: [PATCH 34/36] TailNSigmaThreshold is wrong, should be TailNSigma --- icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl index d01da3070..18240ac18 100644 --- a/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl +++ b/icaruscode/PMT/OpReco/fcl/icarus_ophitfinder.fcl @@ -146,7 +146,7 @@ icarus_opreco_hit_slidingwindow: { ADCThreshold: 5 # ADC threshold (absolute) above pedestal mean to fire a pulse NSigmaThreshold: 3 # ADC threshold (N*pedestal sigma) above pedestal mean to fire a pulse TailADCThreshold: 2 # ADC threshold (absolute) below which next pulse is allowed to fire - TailNSigmaThreshold: 2 # ADC threshold (N*pedestal sigma) below which next pulse is allowed to fire + TailNSigma: 2 # ADC threshold (N*pedestal sigma) below which next pulse is allowed to fire EndADCThreshold: 1 # ADC threshold (absolute) at which the pulse ends EndNSigmaThreshold: 1 # ADC threshold (N*pedetal sigma) at which the pulse ends MinPulseWidth: 2 # The width of a pulse needs to be equal or larger than this to be recorded @@ -162,7 +162,7 @@ icarus_opreco_hit_slidingwindow_201910: { # based on icarus_opreco_hit_slidingwi ADCThreshold: 10 # ADC threshold (absolute) above pedestal mean to fire a pulse NSigmaThreshold: 3 # ADC threshold (N*pedestal sigma) above pedestal mean to fire a pulse TailADCThreshold: 6 # ADC threshold (absolute) below which next pulse is allowed to fire - TailNSigmaThreshold: 2 # ADC threshold (N*pedestal sigma) below which next pulse is allowed to fire + TailNSigma: 2 # ADC threshold (N*pedestal sigma) below which next pulse is allowed to fire EndADCThreshold: 2 # ADC threshold (absolute) at which the pulse ends EndNSigmaThreshold: 1 # ADC threshold (N*pedetal sigma) at which the pulse ends MinPulseWidth: 5 # The width of a pulse needs to be equal or larger than this to be recorded From cec91e5f6fc4b611ae5dbf56d9bd2ec810b4ff66 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 25 Mar 2026 11:52:40 -0400 Subject: [PATCH 35/36] re-tune flash algo for new SPEAreas to match past performance --- icaruscode/PMT/OpReco/fcl/icarus_flashalgo.fcl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/icaruscode/PMT/OpReco/fcl/icarus_flashalgo.fcl b/icaruscode/PMT/OpReco/fcl/icarus_flashalgo.fcl index a85a3f5cd..f8f69fd82 100644 --- a/icaruscode/PMT/OpReco/fcl/icarus_flashalgo.fcl +++ b/icaruscode/PMT/OpReco/fcl/icarus_flashalgo.fcl @@ -25,9 +25,9 @@ SimpleFlashCryo1.OpChannelRange: [180,359] SimpleFlashDataCryo0: @local::SimpleFlashCryo0 -SimpleFlashDataCryo0.PEThreshold: 100 +SimpleFlashDataCryo0.PEThreshold: 200 SimpleFlashDataCryo0.PEThresholdHit: 1.0 -SimpleFlashDataCryo0.MinPECoinc: 100 +SimpleFlashDataCryo0.MinPECoinc: 200 SimpleFlashDataCryo0.MinMultCoinc: 5 SimpleFlashDataCryo0.IntegralTime: 1. SimpleFlashDataCryo0.PreSample: 0.02 @@ -35,9 +35,9 @@ SimpleFlashDataCryo0.VetoSize: 1. SimpleFlashDataCryo0.TimeResolution: 0.01 SimpleFlashDataCryo1: @local::SimpleFlashCryo1 -SimpleFlashDataCryo1.PEThreshold: 100 +SimpleFlashDataCryo1.PEThreshold: 200 SimpleFlashDataCryo1.PEThresholdHit: 1.0 -SimpleFlashDataCryo1.MinPECoinc: 100 +SimpleFlashDataCryo1.MinPECoinc: 200 SimpleFlashDataCryo1.MinMultCoinc: 5 SimpleFlashDataCryo1.IntegralTime: 1. SimpleFlashDataCryo1.PreSample: 0.02 From 7f04d5b278e4f7d1562f8bf63f94d7896e107086 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 25 Mar 2026 18:00:25 -0400 Subject: [PATCH 36/36] add ophit integral to flashana tree --- icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc b/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc index acd96a95e..1fc5a4c75 100644 --- a/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc +++ b/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc @@ -144,6 +144,7 @@ class opana::ICARUSFlashAssAna : public art::EDAnalyzer std::vector &pmt_start_rwm_time, std::vector &pmt_pe, std::vector &pmt_amplitude, + std::vector &pmt_integral, TTree *ophittree); /// Return RWM-relative time from a trigger-relative time @@ -219,6 +220,7 @@ class opana::ICARUSFlashAssAna : public art::EDAnalyzer std::vector m_pmt_start_time_rwm; std::vector m_pmt_pe; std::vector m_pmt_amplitude; + std::vector m_pmt_integral; // Ophit trees int m_channel_id; @@ -409,6 +411,7 @@ void opana::ICARUSFlashAssAna::beginJob() ttree->Branch("time_pmt_rwm", &m_pmt_start_time_rwm); ttree->Branch("pe_pmt", &m_pmt_pe); ttree->Branch("amplitude_pmt", &m_pmt_amplitude); + ttree->Branch("integral_pmt", &m_pmt_integral); fOpFlashTrees.push_back(ttree); @@ -627,6 +630,7 @@ void opana::ICARUSFlashAssAna::processOpHitsFlash(std::vector &pmt_start_time_rwm, std::vector &pmt_pe, std::vector &pmt_amplitude, + std::vector &pmt_integral, TTree *ophittree) { @@ -671,6 +675,7 @@ void opana::ICARUSFlashAssAna::processOpHitsFlash(std::vectorFill(); @@ -869,6 +874,7 @@ void opana::ICARUSFlashAssAna::analyze(art::Event const &e) m_pmt_rise_time.resize(360); m_pmt_start_time_rwm.resize(360); m_pmt_amplitude.resize(360); + m_pmt_integral.resize(360); m_flash_id = idx; m_flash_time = flash.Time(); @@ -894,7 +900,7 @@ void opana::ICARUSFlashAssAna::analyze(art::Event const &e) m_multiplicity_left, m_multiplicity_right, m_sum_pe_left, m_sum_pe_right, m_pmt_start_time, m_pmt_rise_time, m_pmt_start_time_rwm, m_pmt_pe, m_pmt_amplitude, - fOpHitFlashTrees[iFlashLabel]); + m_pmt_integral, fOpHitFlashTrees[iFlashLabel]); m_multiplicity = m_multiplicity_left + m_multiplicity_right; @@ -908,6 +914,8 @@ void opana::ICARUSFlashAssAna::analyze(art::Event const &e) m_pmt_start_time.clear(); m_pmt_rise_time.clear(); m_pmt_start_time_rwm.clear(); + m_pmt_amplitude.clear(); + m_pmt_integral.clear(); idx++; }