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/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/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/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/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/fcl/services/services_icarus_simulation.fcl b/fcl/services/services_icarus_simulation.fcl index bbdc9eeee..51d9e9d60 100644 --- a/fcl/services/services_icarus_simulation.fcl +++ b/fcl/services/services_icarus_simulation.fcl @@ -40,6 +40,8 @@ #include "signalservices_icarus.fcl" #include "photpropservices_icarus.fcl" #include "timing_icarus.fcl" +#include "pmt_channel_status_icarus.fcl" +#include "pmt_calibration_icarus.fcl" BEGIN_PROLOG @@ -115,11 +117,12 @@ icarus_g4_services: { # Define icarus_detsim_services ... (2*) icarus_detsim_services: { - + @table::icarus_detsim_dark_services IPMTTimingCorrectionService: @local::icarus_pmttimingservice - # PmtGainService: @local::icarus_pmtgain_service - + IPMTChannelStatusService: @local::icarus_pmt_channel_status + IPhotonCalibratorService: @local::icarus_photon_calibration + } # icarus_detsim_services 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 diff --git a/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx b/icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx index c8e6f7e8d..facb22b30 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())) + , fBiasConstant(fParams.pulseFunction->biasConstant()) , fPedestalGen(fParams.pedestalGen) , fDiscrAlgo(selectDiscriminationAlgo(fParams.discrimAlgo)) { @@ -156,33 +159,73 @@ 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() //------------------------------------------------------------------------------ -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) - using Fluctuator_t = GainFluctuator; - if (fParams.doGainFluctuations) { - double const refGain = fParams.PMTspecs.firstStageGain(); - return Fluctuator_t - { refGain, CLHEP::RandPoisson{ *fParams.gainRandomEngine, refGain } }; +//------------------------------------------------------------------------------ +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) + + +//------------------------------------------------------------------------------ +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 + // 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 * fBiasConstant / 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() @@ -191,7 +234,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 { @@ -216,6 +260,17 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform timeDelay -= microseconds{ fParams.timingDelays->getCosmicsCorrections(channel) }; } + if (debug) { + debug->opChannel = channel; + 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(); + } + // // collect the amount of photoelectrons arriving at each subtick; // the waveform is split in groups of photons at the same relative subtick @@ -263,6 +318,19 @@ auto icarus::opdet::PMTsimulationAlg::CreateFullWaveform if (tick >= endSample) continue; ++peMaps[subtick][tick]; + + 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(); + dbgPh.subtick = subtick; + debug->photons.push_back(std::move(dbgPh)); + } + } // for photons // auto end = std::chrono::high_resolution_clock::now(); @@ -306,7 +374,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)) { @@ -324,10 +392,21 @@ 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; + 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( @@ -342,7 +421,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; @@ -360,6 +439,12 @@ 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){ + for(std::size_t i=0; iwaveform.push_back(waveform[i].value()); + } + } + return waveform; } // CreateFullWaveform() @@ -670,7 +755,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 +786,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 +930,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 // ----------------------------------------------------------------------------- @@ -906,6 +981,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 +1038,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -972,7 +1049,7 @@ icarus::opdet::PMTsimulationAlgMaker::operator()( { return std::make_unique(makeParams( beamGateTimestamp, - larProp, clockData, timingDelays, + larProp, clockData, timingDelays, gainCalibProvider, SPRfunction, pedestalGenerator, mainRandomEngine, darkNoiseRandomEngine, elecNoiseRandomEngine, trackSelectedPhotons @@ -986,7 +1063,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, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -1010,8 +1088,10 @@ auto icarus::opdet::PMTsimulationAlgMaker::makeParams( params.clockData = &clockData; params.detTimings = detinfo::makeDetectorTimings(params.clockData); params.timingDelays = timingDelays; - + params.gainCalibProvider = gainCalibProvider; + assert(!params.doTimingDelays || params.timingDelays); + 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 d24c19623..4c2771bc2 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/PhotonCalibratorFromDB.h" // LArSoft libraries #include "lardataobj/RawData/OpDetWaveform.h" @@ -46,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 @@ -80,9 +84,36 @@ namespace icarus::opdet { class PMTsimulationAlg; class PMTsimulationAlgMaker; - -} // namespace icarus::opdet + 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) + 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 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 + std::vector> waveform; // waveform data (ADC counts) + }; +} // namespace icarus::opdet // ----------------------------------------------------------------------------- /// Helper class to cut a `raw::OpDetWaveform` from a longer waveform data. @@ -224,6 +255,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 +279,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, 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 + * at any convenient value across run periods. + * * * Dark noise * ----------- @@ -483,6 +531,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 +548,9 @@ class icarus::opdet::PMTsimulationAlg { /// Timing delays service interfacing with database icarusDB::PMTTimingCorrections const* timingDelays = 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. std::optional detTimings; @@ -560,9 +612,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 @@ -612,28 +666,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 ------------------------------------------------ @@ -644,7 +709,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; @@ -674,8 +742,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. @@ -782,7 +850,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; /** @@ -1014,6 +1082,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 +1174,7 @@ class icarus::opdet::PMTsimulationAlgMaker { detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& detClocks, icarusDB::PMTTimingCorrections const* timingDelays, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -1133,6 +1207,7 @@ class icarus::opdet::PMTsimulationAlgMaker { detinfo::LArProperties const& larProp, detinfo::DetectorClocksData const& clockData, icarusDB::PMTTimingCorrections const* timingDelays, + icarusDB::PhotonCalibratorFromDB const* gainCalibProvider, SinglePhotonResponseFunc_t const& SPRfunction, PedestalGenerator_t& pedestalGenerator, CLHEP::HepRandomEngine& mainRandomEngine, @@ -1205,7 +1280,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" @@ -1216,29 +1291,35 @@ 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 << "PulsePolarity: " + << ((fParams.pulsePolarity == 1)? "positive": "negative") + << " (=" << fParams.pulsePolarity << ")" << '\n' << indent << "Sampling: " << fSampling; if (fParams.pulseSubsamples > 1U) out << " (subsampling: x" << fParams.pulseSubsamples << ")"; 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 + " ", ""); 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/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 0082653a5..c407af071 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 @@ -386,7 +388,9 @@ 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 + UseChannelStatusDatabase: false # legacy value + PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage # is proportional to dV^k @@ -449,17 +453,21 @@ 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: 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 - + # 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 @@ -480,17 +488,21 @@ 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: 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 - + # 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 @@ -534,7 +546,9 @@ 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 + UseChannelStatusDatabase: false # legacy value + PMTspecs: { DynodeK: 0.75 # gain on a PMT multiplication stage # is proportional to dV^k diff --git a/icaruscode/PMT/CMakeLists.txt b/icaruscode/PMT/CMakeLists.txt index a34d31b95..0fba647c5 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. @@ -14,6 +15,7 @@ add_subdirectory(Calibration) 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/Calibration/CMakeLists.txt b/icaruscode/PMT/Calibration/CMakeLists.txt index b0ad6f563..7484db3f6 100644 --- a/icaruscode/PMT/Calibration/CMakeLists.txt +++ b/icaruscode/PMT/Calibration/CMakeLists.txt @@ -83,6 +83,28 @@ cet_build_plugin(OpHitRecalibrator art::module icaruscode_Timing ) +cet_build_plugin(VerticalTrackFlashWaveformAna art::module + LIBRARIES + icaruscode::Timing + larcore::Geometry_Geometry_service + larcorealg::Geometry + lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + lardataalg::DetectorInfo + lardataobj::RawData + lardataobj::RecoBase + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::TFileService_service + art_root_io::tfile_support + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib::cetlib + cetlib_except::cetlib_except + ROOT::Tree + ) + art_make_exec(NAME fitPulseDistribution SOURCE fitPulseDistribution.cc diff --git a/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h b/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h index e9c62a189..993f2155d 100644 --- a/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h +++ b/icaruscode/PMT/Calibration/ICARUSPhotonCalibratorServiceFromDB.h @@ -58,9 +58,14 @@ namespace calib { std::string getAreaDatabaseTag() const& { return fProvider.getAreaDatabaseTag(); } - private: + double getSPEArea(int channel) const { return fProvider.getSPEArea(channel); } + + double getSPEFitWidth(int channel) const { return fProvider.getSPEFitWidth(channel); } + provider_type const* provider() const override { return &fProvider; } + private: + void preBeginRun(art::Run const& run); icarusDB::PhotonCalibratorFromDB fProvider; diff --git a/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx b/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx index fd4b3c1bc..9bb7a146c 100644 --- a/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx +++ b/icaruscode/PMT/Calibration/PhotonCalibratorFromDB.cxx @@ -34,11 +34,12 @@ namespace icarusDB { // ----------------------------------------------------------------------------- icarusDB::PhotonCalibratorFromDB::PhotonCalibratorFromDB(const Config& config) - : fCalibDefaults( config.Defaults() ) - , fVerbose ( config.Verbose() ) - , fLogCategory ( config.LogCategory() ) - , fAreaTag ( config.AreaTag() ) - , fDB ( config.DBname(), "", "", config.AreaTag(), true, false) + : fCalibDefaults ( config.Defaults() ) + , fVerbose ( config.Verbose() ) + , fLogCategory ( config.LogCategory() ) + , fAreaTag ( config.AreaTag() ) + , fOverrideRunNumber ( config.OverrideRunNumber() ) + , fDB ( config.DBname(), "", "", config.AreaTag(), true, false) { mf::LogInfo(fLogCategory) << "PhotonCalibratorFromDB connected to " << config.DBname() << " DB, tag '" << config.AreaTag() << "'"; @@ -49,6 +50,13 @@ icarusDB::PhotonCalibratorFromDB::PhotonCalibratorFromDB(const Config& config) void icarusDB::PhotonCalibratorFromDB::readCalibrationFromDB(unsigned int run) { + if (fOverrideRunNumber >= 0) { + mf::LogInfo(fLogCategory) + << "Overriding run number " << run << " with " << fOverrideRunNumber + << " for DB queries."; + run = static_cast(fOverrideRunNumber); + } + mf::LogInfo(fLogCategory) << "Reading SPE area calibrations from database for run " << run; bool ret = fDB.UpdateData( RunToDatabaseTimestamp(run) ); // select table based on run number @@ -71,16 +79,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..fc3ca19b7 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 @@ -115,7 +116,13 @@ class icarusDB::PhotonCalibratorFromDB: public calib::IPhotonCalibrator { Name{ "Defaults" }, Comment{ "values used for channels not present in the database" } }; - + + fhicl::Atom OverrideRunNumber { + Name{ "OverrideRunNumber" }, + Comment{ "if non-negative, use this run number instead of the actual one for DB queries" }, + -1 + }; + }; // Config @@ -151,6 +158,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; @@ -159,6 +174,7 @@ class icarusDB::PhotonCalibratorFromDB: public calib::IPhotonCalibrator { bool const fVerbose; std::string const fLogCategory; std::string const fAreaTag; + int const fOverrideRunNumber; ///< If non-negative, overrides the run number for DB queries. lariov::DBFolder fDB; ///< Cached database interface. diff --git a/icaruscode/PMT/Calibration/VerticalTrackFlashWaveformAna_module.cc b/icaruscode/PMT/Calibration/VerticalTrackFlashWaveformAna_module.cc new file mode 100644 index 000000000..2b4f40ecd --- /dev/null +++ b/icaruscode/PMT/Calibration/VerticalTrackFlashWaveformAna_module.cc @@ -0,0 +1,639 @@ +//////////////////////////////////////////////////////////////////////// +// Class: VerticalTrackFlashWaveformAna +// Plugin Type: analyzer +// File: VerticalTrackFlashWaveformAna_module.cc +// +// Select nearly-vertical cosmic muons inside a Z-slice, match them to a +// recob::OpFlash by YZ-barycenter proximity, and dump for every PMT in +// the flash the OpHit info (integral/amplitude/time/pe) together with +// the single raw::OpDetWaveform covering the flash time on that channel. +// +// Output tree feeds the PMT gain calibration on vertical cosmic muons. +// +// mailto:mvicenzi@bnl.gov +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art_root_io/TFileService.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "larcore/Geometry/WireReadout.h" +#include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom() +#include "larcorealg/Geometry/OpDetGeo.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" +#include "lardataalg/DetectorInfo/DetectorTimings.h" +#include "lardataobj/RawData/OpDetWaveform.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "lardataobj/RecoBase/OpHit.h" +#include "lardataobj/RecoBase/SpacePoint.h" +#include "lardataobj/RecoBase/Track.h" +#include "lardataobj/Simulation/SimPhotons.h" + +#include "icaruscode/Timing/IPMTTimingCorrectionService.h" +#include "icaruscode/Timing/PMTTimingCorrections.h" + +#include "TTree.h" + +#include +#include +#include +#include + +namespace pmtcalib { + class VerticalTrackFlashWaveformAna; +} + +class pmtcalib::VerticalTrackFlashWaveformAna : public art::EDAnalyzer { +public: + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Sequence TrackLabels{ + Name("TrackLabels"), + Comment("recob::Track labels, one per cryostat (E, W)")}; + + fhicl::Sequence OpFlashLabels{ + Name("OpFlashLabels"), + Comment("recob::OpFlash labels, one per cryostat (E, W) in the same order as TrackLabels")}; + + fhicl::Sequence PandoraLabels{ + Name("PandoraLabels"), + Comment("Pandora producer labels (one per cryostat, E/W), used to fetch the" + " hit->SpacePoint association for the track charge barycenter")}; + + fhicl::Atom OpDetWaveformLabel{ + Name("OpDetWaveformLabel"), + Comment("Single raw::OpDetWaveform collection covering both cryostats")}; + + fhicl::Atom SimPhotonsLabel{ + Name("SimPhotonsLabel"), + Comment("Single sim::SimPhoton collection covering both cryostats"), + ""}; + + fhicl::Atom MinAbsDirY{ + Name("MinAbsDirY"), + Comment("Minimum |dir.Y()| for the track to be considered vertical"), + 0.90}; + + fhicl::Atom ZMin{ + Name("ZMin"), + Comment("Minimum Z (cm): both track endpoints must have Z >= ZMin"), + -500.0}; + + fhicl::Atom ZMax{ + Name("ZMax"), + Comment("Maximum Z (cm): both track endpoints must have Z <= ZMax"), + 500.0}; + + fhicl::Atom MinTrackLength{ + Name("MinTrackLength"), + Comment("Minimum track length in cm"), + 100.0}; + + fhicl::Atom RequireDownwards{ + Name("RequireDownwards"), + Comment("If true, only keep tracks going downward (dir.Y() < 0 after optional flip)"), + true}; + + fhicl::Atom MaxDeltaX{ + Name("MaxDeltaX"), + Comment("Max |X_end - X_start| (cm) between track endpoints"), + 20.0}; + + fhicl::Atom MaxDeltaZ{ + Name("MaxDeltaZ"), + Comment("Max |Z_end - Z_start| (cm) between track endpoints"), + 20.0}; + + fhicl::Atom BarycenterMaxDist{ + Name("BarycenterMaxDist"), + Comment("Max YZ distance (cm) between track midpoint and flash barycenter"), + 30.0}; + }; + + using Parameters = art::EDAnalyzer::Table; + + explicit VerticalTrackFlashWaveformAna(Parameters const& config); + + VerticalTrackFlashWaveformAna(VerticalTrackFlashWaveformAna const&) = delete; + VerticalTrackFlashWaveformAna(VerticalTrackFlashWaveformAna&&) = delete; + VerticalTrackFlashWaveformAna& operator=(VerticalTrackFlashWaveformAna const&) = delete; + VerticalTrackFlashWaveformAna& operator=(VerticalTrackFlashWaveformAna&&) = delete; + + void beginJob() override; + void analyze(art::Event const& e) override; + +private: + // --- configuration --- + std::vector fTrackLabels; + std::vector fOpFlashLabels; + std::vector fPandoraLabels; + art::InputTag fOpDetWaveformLabel; + art::InputTag fSimPhotonsLabel; + double fMinAbsDirY; + double fZMin; + double fZMax; + double fMinTrackLength; + bool fRequireDownwards; + double fMaxDeltaX; + double fMaxDeltaZ; + double fBarycenterMaxDist; + + // --- detector timing (optical tick in us) --- + double fOpticalTickUs = 0.002; // ICARUS PMT digitizer: 500 MHz -> 2 ns + + // --- PMT timing corrections service (DB-provided) --- + // OpHit/OpFlash times already include these corrections, but raw + // OpDetWaveform::TimeStamp() does not. We add the same correction to the + // waveform timestamp when checking coverage of a flash time. + icarusDB::PMTTimingCorrections const* fTimingCorrections = nullptr; + + // --- output tree --- + TTree* fTree = nullptr; + int fNSelTracks = 0; + + // event + unsigned int fRun = 0; + unsigned int fEvent = 0; + int fCryo = -1; + + // track + double fTrkStartX = 0., fTrkStartY = 0., fTrkStartZ = 0.; + double fTrkEndX = 0., fTrkEndY = 0., fTrkEndZ = 0.; + double fTrkDirX = 0., fTrkDirY = 0., fTrkDirZ = 0.; + double fTrkLength = 0.; + double fTrkBaryX = 0., fTrkBaryY = 0., fTrkBaryZ = 0.; + + // flash + double fFlashTime = 0.; + double fFlashTotalPE = 0.; + double fFlashY = 0., fFlashZ = 0.; + double fFlashWidthY = 0., fFlashWidthZ = 0.; + int fFlashNPMTs = 0; + double fFlashBarycenterDist = 0.; + + // per-PMT scalars (one tree row per PMT/waveform) + int fPMTChannel = -1; + double fPMTOpHitIntegral = 0.; + double fPMTOpHitAmplitude = 0.; + double fPMTOpHitPeakTime = 0.; + double fPMTOpHitStartTime = 0.; + double fPMTOpHitPE = 0.; + std::vector fPMTWaveform; + double fPMTWaveformTimestamp = 0.; + double fPMTPosX = 0., fPMTPosY = 0., fPMTPosZ = 0.; + double fPMTBarycenterDist = 0.; + + // simphotons + int fNSimPhotons = 0; + std::vector fSimPhotonTimes; + std::vector fSimPhotonInitX; + std::vector fSimPhotonInitY; + std::vector fSimPhotonInitZ; + + // --- helpers (stubs; to be filled in next iteration) --- + + /// Initialise branches on the output tree. + void setupTree(); + + /// Clear all per-entry buffers. + void resetEntry(); + + /// Select a nearly-vertical track inside the Z-slice. + /// Returns true and fills start/end/dir/length into members if selected. + bool selectVerticalTrack(recob::Track const& track); + + /// Find the OpFlash in `flashes` whose YZ-barycenter is closest to + /// (trkY, trkZ); return its index (or -1 if none within BarycenterMaxDist). + /// On success sets fFlashBarycenterDist. + int matchFlashByBarycenter(std::vector> const& flashes, + double trkY, double trkZ); + + /// Total timing correction (us) to apply to a raw waveform TimeStamp so + /// that it lives in the same reference frame as OpHit/OpFlash times. + /// Corrections are additive; the sum matches what OpHitTimingCorrection + /// applies to OpHits. + double getTimingCorrection(raw::Channel_t channel) const; + + /// Scan the full waveform collection and return the raw::OpDetWaveform on + /// `channel` whose time window [TimeStamp+corr, TimeStamp+corr+size*tick) + /// contains `flashTime`. `correction` is the per-channel correction already + /// extracted once by the caller. Returns nullptr if none found. + raw::OpDetWaveform const* findCoveringWaveform( + std::vector const& wfs, + raw::Channel_t channel, + double flashTime, + double correction) const; + + /// Fill per-PMT vectors for the matched flash by looping the waveform + /// collection once per PMT channel in the flash. + void fillFlashPMTs( + recob::OpFlash const& flash, + std::vector> const& ophits, + std::vector const& wfs, + std::vector const& simphotonsCollection, + detinfo::DetectorTimings const& timings); +}; + +// ===================================================================== + +pmtcalib::VerticalTrackFlashWaveformAna::VerticalTrackFlashWaveformAna(Parameters const& config) + : art::EDAnalyzer{config} + , fTrackLabels{config().TrackLabels()} + , fOpFlashLabels{config().OpFlashLabels()} + , fPandoraLabels{config().PandoraLabels()} + , fOpDetWaveformLabel{config().OpDetWaveformLabel()} + , fSimPhotonsLabel{config().SimPhotonsLabel()} + , fMinAbsDirY{config().MinAbsDirY()} + , fZMin{config().ZMin()} + , fZMax{config().ZMax()} + , fMinTrackLength{config().MinTrackLength()} + , fRequireDownwards{config().RequireDownwards()} + , fMaxDeltaX{config().MaxDeltaX()} + , fMaxDeltaZ{config().MaxDeltaZ()} + , fBarycenterMaxDist{config().BarycenterMaxDist()} + , fTimingCorrections{lar::providerFrom()} +{ + if (fTrackLabels.size() != fOpFlashLabels.size() + || fTrackLabels.size() != fPandoraLabels.size()) { + throw art::Exception(art::errors::Configuration) + << "TrackLabels, OpFlashLabels and PandoraLabels must all have the same" + " size (one per cryostat)."; + } +} + +void pmtcalib::VerticalTrackFlashWaveformAna::beginJob() +{ + setupTree(); + fNSelTracks = 0; +} + +void pmtcalib::VerticalTrackFlashWaveformAna::analyze(art::Event const& e) +{ + fRun = e.run(); + fEvent = e.event(); + + // Refresh the optical tick period from DetectorClocksService each event + // (cheap and keeps data/MC config consistent). + auto const clocks = + art::ServiceHandle()->DataFor(e); + fOpticalTickUs = clocks.OpticalClock().TickPeriod(); + + // Single collections spanning both cryostats. + auto const& wfs = + *e.getValidHandle>(fOpDetWaveformLabel); + + auto const timings = detinfo::makeDetectorTimings(clocks); + + static const std::vector emptySimPhotons; + auto const & simphotonsCollection = (!fSimPhotonsLabel.empty()) + ? e.getProduct>(fSimPhotonsLabel) + : emptySimPhotons; + + for (std::size_t icryo = 0; icryo < fTrackLabels.size(); ++icryo) { + auto const& trackHandle = + e.getValidHandle>(fTrackLabels[icryo]); + auto const& tracks = *trackHandle; + auto const& flashHandle = + e.getValidHandle>(fOpFlashLabels[icryo]); + auto const& flashes = *flashHandle; + + if (flashes.empty()) continue; + + // Promote OpFlash vector to Ptrs (for the matcher signature) and fetch + // the flash->OpHit associations once per cryostat. + std::vector> flashPtrs; + flashPtrs.reserve(flashes.size()); + for (std::size_t i = 0; i < flashes.size(); ++i) + flashPtrs.emplace_back(flashHandle, i); + + art::FindManyP ophitsPtr(flashHandle, e, fOpFlashLabels[icryo]); + + // Track -> Hit association, from the track producer (pattern borrowed + // from sbncode TrackCaloSkimmer_module). + art::FindManyP fmtrkHits(trackHandle, e, fTrackLabels[icryo]); + + std::vector> const emptyHitVector; + + for (std::size_t itrk = 0; itrk < tracks.size(); ++itrk) { + resetEntry(); + if (!selectVerticalTrack(tracks[itrk])) continue; + + art::Ptr const trkPtr(trackHandle, itrk); + std::vector> const& trkHits = + fmtrkHits.isValid() ? fmtrkHits.at(trkPtr.key()) : emptyHitVector; + + // Charge-weighted track barycenter using hit Integral as weight and + // the associated SpacePoint(s) for 3D position — same chain as + // sbncode TrackCaloSkimmer: Hit -> SpacePoint assns live under the + // pandora (PFP) producer label. + art::FindManyP fmtrkHitSPs(trkHits, e, fPandoraLabels[icryo]); + + double sumCharge = 0.; + double sumX = 0., sumY = 0., sumZ = 0.; + for (std::size_t k = 0; k < trkHits.size(); ++k) { + auto const& sps = fmtrkHitSPs.at(k); + if (sps.empty()) continue; + double const q = trkHits[k]->Integral(); + auto const* xyz = sps.front()->XYZ(); + sumCharge += q; + sumX += xyz[0] * q; + sumY += xyz[1] * q; + sumZ += xyz[2] * q; + } + if (sumCharge <= 0.) continue; + fTrkBaryX = sumX / sumCharge; + fTrkBaryY = sumY / sumCharge; + fTrkBaryZ = sumZ / sumCharge; + + int const iflash = matchFlashByBarycenter(flashPtrs, fTrkBaryY, fTrkBaryZ); + if (iflash < 0) continue; + + auto const& flash = flashes[iflash]; + fFlashTime = flash.Time(); + fFlashTotalPE = flash.TotalPE(); + fFlashY = flash.YCenter(); + fFlashZ = flash.ZCenter(); + fFlashWidthY = flash.YWidth(); + fFlashWidthZ = flash.ZWidth(); + + mf::LogInfo("VerticalTrackFlashWaveformAna") << "Vertical track matched to flash: " + << "\n TrackYDir: " << fTrkDirY + << "\n TrackX: " << fTrkStartX << " " << fTrkBaryX << " " << fTrkEndX + << "\n TrackY: " << fTrkStartY << " " << fTrkBaryY << " " << fTrkEndY + << "\n TrackZ: " << fTrkStartZ << " " << fTrkBaryZ << " " << fTrkEndZ + << "\n TrkLength: " << fTrkLength + << "\n FlashTime: " << fFlashTime + << "\n FlashY: " << fFlashY + << "\n FlashZ: " << fFlashZ + << "\n FlashPE: " << fFlashTotalPE; + + fNSelTracks++; + fCryo = static_cast(icryo); + + auto const& ophits = ophitsPtr.at(iflash); + fillFlashPMTs(flash, ophits, wfs, simphotonsCollection, timings); + + } + } + + mf::LogInfo("VerticalTrackFlashWaveformAna") << "Total selected tracks: " << fNSelTracks; +} + +// --- helper stubs (bodies to be added in the next step) ----------------- + +void pmtcalib::VerticalTrackFlashWaveformAna::setupTree() +{ + art::ServiceHandle tfs; + fTree = tfs->make("vertical", "Vertical muon / flash / PMT waveform tree"); + + fTree->Branch("run", &fRun); + fTree->Branch("event", &fEvent); + fTree->Branch("cryo", &fCryo); + + fTree->Branch("trk_start_x", &fTrkStartX); + fTree->Branch("trk_start_y", &fTrkStartY); + fTree->Branch("trk_start_z", &fTrkStartZ); + fTree->Branch("trk_end_x", &fTrkEndX); + fTree->Branch("trk_end_y", &fTrkEndY); + fTree->Branch("trk_end_z", &fTrkEndZ); + fTree->Branch("trk_dir_x", &fTrkDirX); + fTree->Branch("trk_dir_y", &fTrkDirY); + fTree->Branch("trk_dir_z", &fTrkDirZ); + fTree->Branch("trk_length", &fTrkLength); + fTree->Branch("trk_barycenter_x", &fTrkBaryX); + fTree->Branch("trk_barycenter_y", &fTrkBaryY); + fTree->Branch("trk_barycenter_z", &fTrkBaryZ); + + fTree->Branch("flash_time", &fFlashTime); + fTree->Branch("flash_total_pe", &fFlashTotalPE); + fTree->Branch("flash_y", &fFlashY); + fTree->Branch("flash_z", &fFlashZ); + fTree->Branch("flash_width_y", &fFlashWidthY); + fTree->Branch("flash_width_z", &fFlashWidthZ); + fTree->Branch("flash_n_pmts", &fFlashNPMTs); + fTree->Branch("flash_barycenter_dist", &fFlashBarycenterDist); + + fTree->Branch("pmt_channel", &fPMTChannel, "pmt_channel/I"); + fTree->Branch("pmt_ophit_integral", &fPMTOpHitIntegral, "pmt_ophit_integral/D"); + fTree->Branch("pmt_ophit_amplitude", &fPMTOpHitAmplitude, "pmt_ophit_amplitude/D"); + fTree->Branch("pmt_ophit_peak_time", &fPMTOpHitPeakTime, "pmt_ophit_peak_time/D"); + fTree->Branch("pmt_ophit_start_time", &fPMTOpHitStartTime, "pmt_ophit_start_time/D"); + fTree->Branch("pmt_ophit_pe", &fPMTOpHitPE, "pmt_ophit_pe/D"); + fTree->Branch("pmt_waveform", &fPMTWaveform); + fTree->Branch("pmt_waveform_timestamp", &fPMTWaveformTimestamp, "pmt_waveform_timestamp/D"); + fTree->Branch("pmt_pos_x", &fPMTPosX, "pmt_pos_x/D"); + fTree->Branch("pmt_pos_y", &fPMTPosY, "pmt_pos_y/D"); + fTree->Branch("pmt_pos_z", &fPMTPosZ, "pmt_pos_z/D"); + fTree->Branch("pmt_barycenter_dist", &fPMTBarycenterDist, "pmt_barycenter_dist/D"); + + if(!fSimPhotonsLabel.empty()) + { + fTree->Branch("n_sim_photons", &fNSimPhotons, "n_sim_photons/I"); + fTree->Branch("sim_photon_times", &fSimPhotonTimes); + fTree->Branch("sim_photon_init_x", &fSimPhotonInitX); + fTree->Branch("sim_photon_init_y", &fSimPhotonInitY); + fTree->Branch("sim_photon_init_z", &fSimPhotonInitZ); + } + +} + +void pmtcalib::VerticalTrackFlashWaveformAna::resetEntry() +{ + fCryo = -1; + fTrkStartX = fTrkStartY = fTrkStartZ = 0.; + fTrkEndX = fTrkEndY = fTrkEndZ = 0.; + fTrkDirX = fTrkDirY = fTrkDirZ = 0.; + fTrkLength = 0.; + fTrkBaryX = fTrkBaryY = fTrkBaryZ = 0.; + + fFlashTime = fFlashTotalPE = 0.; + fFlashY = fFlashZ = fFlashWidthY = fFlashWidthZ = 0.; + fFlashNPMTs = 0; + fFlashBarycenterDist = 0.; + + fPMTChannel = -1; + fPMTOpHitIntegral = 0.; + fPMTOpHitAmplitude = 0.; + fPMTOpHitPeakTime = 0.; + fPMTOpHitStartTime = 0.; + fPMTOpHitPE = 0.; + fPMTWaveform.clear(); + fPMTWaveformTimestamp = 0.; + fPMTPosX = fPMTPosY = fPMTPosZ = 0.; + + fNSimPhotons = 0; + fSimPhotonTimes.clear(); + fSimPhotonInitX.clear(); + fSimPhotonInitY.clear(); + fSimPhotonInitZ.clear(); +} + +bool pmtcalib::VerticalTrackFlashWaveformAna::selectVerticalTrack(recob::Track const& track) +{ + + if (track.Length() < fMinTrackLength) return false; + + auto const start = track.Start(); + auto const end = track.End(); + + if (start.Z() < fZMin || start.Z() > fZMax) return false; + if (end.Z() < fZMin || end.Z() > fZMax) return false; + if (std::abs(end.X() - start.X()) > fMaxDeltaX) return false; + if (std::abs(end.Z() - start.Z()) > fMaxDeltaZ) return false; + + auto dir = track.StartDirection(); + // Flip to downgoing if requested (same idiom as TimeTrackTreeStorageCRT). + if (fRequireDownwards && dir.Y() > 0.0) dir = -track.EndDirection(); + if (fRequireDownwards && dir.Y() >= 0.0) return false; + if (std::abs(dir.Y()) < fMinAbsDirY) return false; + + fTrkStartX = start.X(); fTrkStartY = start.Y(); fTrkStartZ = start.Z(); + fTrkEndX = end.X(); fTrkEndY = end.Y(); fTrkEndZ = end.Z(); + fTrkDirX = dir.X(); fTrkDirY = dir.Y(); fTrkDirZ = dir.Z(); + fTrkLength = track.Length(); + + return true; +} + +int pmtcalib::VerticalTrackFlashWaveformAna::matchFlashByBarycenter( + std::vector> const& flashes, + double trkY, double trkZ) +{ + int best = -1; + double bestDist = std::numeric_limits::max(); + for (std::size_t i = 0; i < flashes.size(); ++i) { + double const dy = flashes[i]->YCenter() - trkY; + double const dz = flashes[i]->ZCenter() - trkZ; + double const d = std::sqrt(dy * dy + dz * dz); + if (d < bestDist) { bestDist = d; best = static_cast(i); } + } + if (best < 0 || bestDist > fBarycenterMaxDist) return -1; + fFlashBarycenterDist = bestDist; + return best; +} + +double pmtcalib::VerticalTrackFlashWaveformAna::getTimingCorrection( + raw::Channel_t channel) const +{ + // Same additive combination applied to OpHits in OpHitTimingCorrection. + return fTimingCorrections->getLaserCorrections(channel) + + fTimingCorrections->getCosmicsCorrections(channel); +} + +raw::OpDetWaveform const* +pmtcalib::VerticalTrackFlashWaveformAna::findCoveringWaveform( + std::vector const& wfs, + raw::Channel_t channel, + double flashTime, + double correction) const +{ + // Linear scan (performance is not a concern). Both OpDetWaveform::TimeStamp() + // and OpFlash::Time() are in microseconds. The raw timestamp is un-corrected, + // while flashTime is already corrected, so shift the waveform window by the + // per-channel correction before the containment check. + for (auto const& wf : wfs) { + if (wf.ChannelNumber() != channel) continue; + double const t0 = wf.TimeStamp() + correction; + double const t1 = t0 + wf.size() * fOpticalTickUs; + if (flashTime >= t0 && flashTime < t1) return &wf; + } + return nullptr; +} + +void pmtcalib::VerticalTrackFlashWaveformAna::fillFlashPMTs( + recob::OpFlash const& flash, + std::vector> const& ophits, + std::vector const& wfs, + std::vector const& simphotonsCollection, + detinfo::DetectorTimings const& timings) +{ + auto const& channelMap = art::ServiceHandle()->Get(); + + // A flash can have more than one OpHit on the same channel. Following + // ICARUSFlashAssAna, keep only the earliest-in-time ophit per channel and + // record its values for that single hit. + std::map firstHit; + for (auto const& hit : ophits) { + auto [it, inserted] = firstHit.try_emplace(hit->OpChannel(), hit.get()); + if (!inserted && hit->StartTime() < it->second->StartTime()) + it->second = hit.get(); + } + + fFlashNPMTs = static_cast(firstHit.size()); + + for (auto const& [ch, hitPtr] : firstHit) { + recob::OpHit const& hit = *hitPtr; + + // Extract the per-channel correction once, then use it to shift the + // raw waveform timestamp into the same frame as the flash time. + double const correction = getTimingCorrection(ch); + + raw::OpDetWaveform const* wf = findCoveringWaveform(wfs, ch, flash.AbsTime(), correction); + if (!wf) { + mf::LogWarning("VerticalTrackFlashWaveformAna") + << "No raw::OpDetWaveform covering flash AbsTime=" << flash.AbsTime() + << " us on channel " << ch; + continue; + } + + auto const pmtPos = channelMap.OpDetGeoFromOpChannel(ch).GetCenter(); + + fPMTChannel = static_cast(ch); + fPMTOpHitIntegral = hit.Area(); + fPMTOpHitAmplitude = hit.Amplitude(); + fPMTOpHitPeakTime = hit.PeakTime(); + fPMTOpHitStartTime = hit.StartTime(); + fPMTOpHitPE = hit.PE(); + fPMTWaveform = wf->Waveform(); + fPMTWaveformTimestamp = wf->TimeStamp() + getTimingCorrection(fPMTChannel); + fPMTPosX = pmtPos.X(); + fPMTPosY = pmtPos.Y(); + fPMTPosZ = pmtPos.Z(); + fPMTBarycenterDist = std::hypot( fPMTPosX - fTrkBaryX, fPMTPosY - fTrkBaryY, fPMTPosZ - fTrkBaryZ ); + + // SimPhotons belonging to this channel and falling inside the current + // waveform window. Collection is empty when no SimPhotons product is + // available (data, or MC without photon propagation). + fNSimPhotons = 0; + fSimPhotonTimes.clear(); + fSimPhotonInitX.clear(); + fSimPhotonInitY.clear(); + fSimPhotonInitZ.clear(); + + double const wfstart = fPMTWaveformTimestamp; + double const wfend = wfstart + wf->size() * fOpticalTickUs; + + for (auto const& simphotons : simphotonsCollection) { + + if (simphotons.OpChannel() != fPMTChannel) continue; + + for (auto const& ph : simphotons) { + detinfo::timescales::simulation_time const photonTime{ ph.Time }; + double const t = timings.toElectronicsTime(photonTime).value(); + + if (t < wfstart || t > wfend) continue; + + fSimPhotonTimes.push_back(t); + fSimPhotonInitX.push_back(ph.InitialPosition.X()); + fSimPhotonInitY.push_back(ph.InitialPosition.Y()); + fSimPhotonInitZ.push_back(ph.InitialPosition.Z()); + } + } + fNSimPhotons = static_cast(fSimPhotonTimes.size()); + + fTree->Fill(); + } +} + +DEFINE_ART_MODULE(pmtcalib::VerticalTrackFlashWaveformAna) diff --git a/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl b/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl index 39a2d2efe..e657c1f51 100644 --- a/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl +++ b/icaruscode/PMT/Calibration/fcl/pmt_calibration_icarus.fcl @@ -48,13 +48,44 @@ pmt_bkgphotons_calibration: icarus_photon_calibration: { service_provider: ICARUSPhotonCalibratorServiceFromDB - AreaTag: "v1r0" + AreaTag: "v1r2" Defaults: { SPEArea: @local::SPERun2.Area # ADC x tick, Run2 default } + OverrideRunNumber: 9271 } # icarus_photon_calibration +# ============================================================================== +# === Vertical-track flash / waveform ntupler (for calibration) +# === +verticaltrack_flashwf_ana: +{ + module_type: "VerticalTrackFlashWaveformAna" + + # per-cryostat track + flash collections (order must match: E, W) + TrackLabels: [ "pandoraTrackGausCryoE", "pandoraTrackGausCryoW" ] + OpFlashLabels: [ "opflashCryoE", "opflashCryoW" ] + PandoraLabels: [ "pandoraGausCryoE", "pandoraGausCryoW" ] + + # single merged collections covering both cryostats + OpDetWaveformLabel: "daqPMT" + SimPhotonsLabel: "" + + # vertical-track selection + MinAbsDirY: 0.98 + MaxDeltaX: 30.0 # cm + MaxDeltaZ: 30.0 # cm + ZMin: -700.0 # cm + ZMax: 700.0 # cm + MinTrackLength: 100.0 # cm + RequireDownwards: true + + # flash match + BarycenterMaxDist: 30.0 # cm (YZ distance) +} + + # ============================================================================== END_PROLOG diff --git a/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h b/icaruscode/PMT/OpReco/Algorithms/OpRecoFactoryStuff.h index 8a1173f76..053481960 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,19 @@ 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 + { + 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 @@ -534,6 +563,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/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++; } diff --git a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc index f9ae98c29..94f63c16d 100644 --- a/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc +++ b/icaruscode/PMT/OpReco/ICARUSOpHitFinder_module.cc @@ -8,11 +8,13 @@ */ // 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" #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 +47,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 +248,18 @@ 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)."} + }; + fhicl::Atom UseCalibrator { Name{ "UseCalibrator" }, Comment{ "Use the photon calibration service configured in the job" }, @@ -260,21 +269,24 @@ 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 { + Name{ "UseChannelStatusDatabase" }, + Comment{"ignore waveforms on channels flagged in PMT channel status database."}, + false + }; + }; // Config using Parameters = art::ReplicatedProducer::Table; @@ -298,12 +310,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 ---------------------------------------- @@ -325,6 +340,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 +513,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()) } { @@ -511,14 +531,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 // @@ -527,8 +548,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, @@ -550,7 +579,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()); @@ -577,17 +606,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 @@ -652,20 +682,34 @@ 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() +// ----------------------------------------------------------------------------- +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_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 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 650351356..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 @@ -186,26 +186,38 @@ 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 entry point to the reconstruction "ICARUSOpHitFinder" + ## this allows interfacing with our PMT channel status service + module_type: "ICARUSOpHitFinder" GenModule: "generator" InputModule: "opdaq" - InputLabels: [] - ChannelMasks: [] - UseCalibrator: true # use calibrator service (SPEArea database) + #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 + + #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 + tool_type: RiseTimeThreshold + PeakRatio: 0.15 # at 15% of the peak amplitude -- not tuned + #InterpolateSamples: true # sub-sample linear interpolation } } @@ -218,6 +230,7 @@ icarus_ophitdebugger.OutputFile: "ophit_debug.root" icarus_ophit_MC: { @table::icarus_ophit UseCalibrator: false + UseChannelStatusDatabase: false InputModule: "opdaq" } 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() diff --git a/icaruscode/PMT/SimPMTIcarus_module.cc b/icaruscode/PMT/SimPMTIcarus_module.cc index 9bf70fefe..e5cab132b 100644 --- a/icaruscode/PMT/SimPMTIcarus_module.cc +++ b/icaruscode/PMT/SimPMTIcarus_module.cc @@ -17,6 +17,8 @@ #include "icaruscode/PMT/Algorithms/PhotoelectronPulseFunction.h" #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" @@ -30,6 +32,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" @@ -47,6 +50,9 @@ #include "fhiclcpp/types/Atom.h" #include "fhiclcpp/ParameterSet.h" +// ROOT libraries +#include "TTree.h" + // CLHEP libraries #include "CLHEP/Random/RandEngine.h" // CLHEP::HepRandomEngine @@ -59,7 +65,6 @@ #include // std::move() #include - namespace icarus::opdet { /** @@ -200,6 +205,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"), @@ -240,6 +251,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 +271,7 @@ namespace icarus::opdet { SimPMTIcarus & operator = (SimPMTIcarus &&) = delete; // Required functions. + virtual void beginJob() override; void produce(art::Event & e) override; private: @@ -272,6 +290,8 @@ 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. + bool fUseChannelStatusDB; ///< Whether to skip non-ON channels using the status DB. CLHEP::HepRandomEngine& fEfficiencyEngine; CLHEP::HepRandomEngine& fDarkNoiseEngine; @@ -286,7 +306,6 @@ namespace icarus::opdet { /// The actual simulation algorithm. icarus::opdet::PMTsimulationAlgMaker makePMTsimulator; - /// True if `firstTime()` has already been called. std::atomic_flag fNotFirstTime; @@ -303,7 +322,44 @@ 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 (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 @@ -317,6 +373,8 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) , fMakeMetadata(config().MakeMetadata()) , fWritePhotons(config().writePhotons()) , fDoTimingDelays(config().algoConfig().ApplyTimingDelays()) + , fUseGainCalibDB(config().algoConfig().UseGainDatabase()) + , fUseChannelStatusDB(config().UseChannelStatusDatabase()) // random engines , fEfficiencyEngine(art::ServiceHandle()->registerAndSeedEngine( createEngine(0, "HepJamesRandom", "Efficiencies"), @@ -348,6 +406,8 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) ->makeGenerator(fElectronicsNoiseEngine) } , makePMTsimulator(config().algoConfig()) + , fDebugTree(config().DebugTree()) + { // Call appropriate produces<>() functions here. produces>(); @@ -360,6 +420,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("timeDelay_us", &timeDelay_us, "timeDelay_us/F"); + fTree->Branch("triggerOffsetPMT_us", &triggerOffsetPMT_us, "triggerOffsetPMT_us/F"); + 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_waveformTime_us", &photon_waveformTime_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) @@ -394,6 +486,7 @@ SimPMTIcarus::SimPMTIcarus(Parameters const& config) *(lar::providerFrom()), clockData, fDoTimingDelays ? lar::providerFrom() : nullptr, + fUseGainCalibDB ? art::ServiceHandle()->provider() : nullptr, *fSinglePhotonResponseFunc, *fPedestalGen, fEfficiencyEngine, @@ -405,6 +498,7 @@ 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); } // if first time @@ -415,37 +509,91 @@ 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()); - 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 (simphVecPtr && photons_used) simphVecPtr->emplace_back(std::move(photons_used.value())); + + if(fDebugTree && debug && debug->photons.size()>0) { + // fill debug tree variables + 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 full-window waveform + waveform_adc = debug->waveform; + // fill per photon info + nInputPhotons = photons.size(); + nDetectedPhotons = debug->photons.size(); + photon_simTime_ns.clear(); + photon_waveformTime_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_waveformTime_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 } else if(pmtLiteVector.isValid()) { 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); - auto const& [ channelWaveforms, photons_used ] + auto const& [ channelWaveforms, photons_used, debug] = PMTsimulator->simulate(photons, lite_photons); std::move( channelWaveforms.cbegin(), channelWaveforms.cend(), std::back_inserter(*pulseVecPtr) ); - } // for } diff --git a/icaruscode/PMT/Status/CMakeLists.txt b/icaruscode/PMT/Status/CMakeLists.txt new file mode 100644 index 000000000..73581e3de --- /dev/null +++ b/icaruscode/PMT/Status/CMakeLists.txt @@ -0,0 +1,27 @@ +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 + art::Framework_Principal + larcore::Geometry_Geometry_service +) + +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..86241ce2d --- /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). + kGood = 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: kGood, kOff, or kBad + * - convenience predicates: isGood(), 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 (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; } + + /// 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 kGood. + 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..ba4e3cc1a --- /dev/null +++ b/icaruscode/PMT/Status/PMTChannelStatusProvider.cxx @@ -0,0 +1,140 @@ +/** + * @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", "" ) ) + , fOverrideRunNumber( pset.get("OverrideRunNumber", -1) ) + , fDB ( pset.get("DBname", "pmt_voltage_data"), + "", "", fStatusTag, true, false ) +{ + int const defaultStatusInt = pset.get("DefaultStatus", static_cast(kGood)); + 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 << "'."; +} + + +// ----------------------------------------------------------------------------- +void icarusDB::PMTChannelStatusProvider::readStatusFromDB(unsigned int run) +{ + + if (fOverrideRunNumber >= 0) { + mf::LogInfo(fLogCategory) + << "Overriding run number " << run << " with " << fOverrideRunNumber + << " for DB queries."; + run = static_cast(fOverrideRunNumber); + } + + 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..85873b2a6 --- /dev/null +++ b/icaruscode/PMT/Status/PMTChannelStatusProvider.h @@ -0,0 +1,113 @@ +/** + * @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, kGood=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(kGood); } + 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; + int fOverrideRunNumber; ///< If non-negative, overrides the run number for DB queries. + + 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..cfaf53788 --- /dev/null +++ b/icaruscode/PMT/Status/pmt_channel_status_icarus.fcl @@ -0,0 +1,19 @@ +# +# @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 # kGood + DefaultVoltage: 1500.0 # [V] + Verbose: false + OverrideRunNumber: -1 + LogCategory: "PMTChannelStatusProvider" +} + +END_PROLOG diff --git a/icaruscode/Timing/PMTTimingCorrectionsProvider.cxx b/icaruscode/Timing/PMTTimingCorrectionsProvider.cxx index 17dcd12f3..ca618a2bb 100644 --- a/icaruscode/Timing/PMTTimingCorrectionsProvider.cxx +++ b/icaruscode/Timing/PMTTimingCorrectionsProvider.cxx @@ -25,10 +25,11 @@ //-------------------------------------------------------------------------------- icarusDB::PMTTimingCorrectionsProvider::PMTTimingCorrectionsProvider - (const fhicl::ParameterSet& pset) + (const fhicl::ParameterSet& pset) : fVerbose{ pset.get("Verbose", false) } , fLogCategory{ pset.get("LogCategory", "PMTTimingCorrection") } - { + , fOverrideRunNumber{ pset.get("OverrideRunNumber", -1) } + { fhicl::ParameterSet const tags{ pset.get("CorrectionTags") }; fCablesTag = tags.get("CablesTag"); fLaserTag = tags.get("LaserTag"); @@ -210,12 +211,21 @@ void icarusDB::PMTTimingCorrectionsProvider::ReadCosmicsCorrections( uint32_t ru /// is the PMT channel number void icarusDB::PMTTimingCorrectionsProvider::readTimeCorrectionDatabase(const art::Run& run){ + uint32_t runNumber = run.id().run(); + + if (fOverrideRunNumber >= 0) { + mf::LogInfo(fLogCategory) + << "Overriding run number " << runNumber << " with " << fOverrideRunNumber + << " for DB queries."; + runNumber = static_cast(fOverrideRunNumber); + } + // Clear before the run fDatabaseTimingCorrections.clear(); - ReadPMTCablesCorrections(run.id().run()); - ReadLaserCorrections(run.id().run()); - ReadCosmicsCorrections(run.id().run()); + ReadPMTCablesCorrections(runNumber); + ReadLaserCorrections(runNumber); + ReadCosmicsCorrections(runNumber); if( fVerbose ) { diff --git a/icaruscode/Timing/PMTTimingCorrectionsProvider.h b/icaruscode/Timing/PMTTimingCorrectionsProvider.h index 3ef53a007..ee8002bd1 100644 --- a/icaruscode/Timing/PMTTimingCorrectionsProvider.h +++ b/icaruscode/Timing/PMTTimingCorrectionsProvider.h @@ -100,9 +100,10 @@ class icarusDB::PMTTimingCorrectionsProvider : public PMTTimingCorrections { bool fVerbose = false; ///< Whether to print the configuration we read. std::string fLogCategory; ///< Category tag for messages. - std::string fCablesTag; ///< Tag for cable corrections database. + int fOverrideRunNumber = -1; ///< If non-negative, overrides the run number for DB queries. + std::string fCablesTag; ///< Tag for cable corrections database. std::string fLaserTag; ///< Tag for laser corrections database. - std::string fCosmicsTag; ///< Tag for cosmics corrections database. + std::string fCosmicsTag; ///< Tag for cosmics corrections database. /// Map of corrections by channel std::map fDatabaseTimingCorrections; diff --git a/icaruscode/Timing/timing_icarus.fcl b/icaruscode/Timing/timing_icarus.fcl index c09df28b3..b0291232d 100644 --- a/icaruscode/Timing/timing_icarus.fcl +++ b/icaruscode/Timing/timing_icarus.fcl @@ -12,6 +12,7 @@ icarus_pmttimingservice: CosmicsTag: @local::ICARUS_Calibration_GlobalTags.pmt_cosmics_timing_data } Verbose: false + OverrideRunNumber: -1 } icarus_ophit_timing_correction: