From 7cc503fbd9943e9bd5292c1f14fd0a99018b922f Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Wed, 9 Jul 2025 04:13:19 -0500 Subject: [PATCH 1/8] setting up siom --- EventGenerator/src/PhotonGun_module.cc | 22 +++++++++- STMMC/fcl/Absorber.fcl | 11 +++-- STMMC/fcl/AbsorberFromST.fcl | 59 ++++++++++++++++++++++++++ STMMC/fcl/prolog.fcl | 39 ++++++++++++++++- 4 files changed, 121 insertions(+), 10 deletions(-) create mode 100644 STMMC/fcl/AbsorberFromST.fcl diff --git a/EventGenerator/src/PhotonGun_module.cc b/EventGenerator/src/PhotonGun_module.cc index 6b9e4359e2..a20da245be 100644 --- a/EventGenerator/src/PhotonGun_module.cc +++ b/EventGenerator/src/PhotonGun_module.cc @@ -38,6 +38,9 @@ namespace mu2e { fhicl::Atom z{ Name("z"), Comment("z position of generated photon")}; fhicl::OptionalAtom px{ Name("px"), Comment("x momentum of generated photon")}; fhicl::OptionalAtom py{ Name("py"), Comment("y momentum of generated photon")}; + fhicl::OptionalAtom deltax{ Name("deltax"), Comment("Difference in x position of generated photon, use to calculate targeted position")}; + fhicl::OptionalAtom deltay{ Name("deltay"), Comment("Difference in y position of generated photon, use to calculate targeted position")}; + fhicl::OptionalAtom deltaz{ Name("deltaz"), Comment("Difference in z position of generated photon, use to calculate targeted position")}; fhicl::Atom E{ Name("E"), Comment("Energy of generated photon")}; }; using Parameters = art::EDProducer::Table; @@ -46,6 +49,7 @@ namespace mu2e { private: double x = 0.0, y = 0.0, z = 0.0; double px = 0.0, py = 0.0, pz = 0.0; + double deltax = 0.0, deltay = 0.0, deltaz = 0.0; double E = 0.0; }; @@ -56,17 +60,31 @@ namespace mu2e { z(conf().z()), E(conf().E()) { produces(); + if (E < std::numeric_limits::epsilon()) + throw cet::exception("RANGE") << "Energy must be greater than zero, exiting."; px = conf().px() ? *conf().px() : 0; - px = conf().py() ? *conf().py() : 0; + py = conf().py() ? *conf().py() : 0; if ((px*px + py*py) > (E*E)) throw cet::exception("RANGE") << "magnitude of px and py is greater than E, exiting."; pz = std::sqrt(E*E - px*px - py*py); + deltax = conf().deltax() ? *conf().deltax() : 0; + deltay = conf().deltay() ? *conf().deltay() : 0; + deltaz = conf().deltaz() ? *conf().deltaz() : 0; + if ((px > std::numeric_limits::epsilon() || py > std::numeric_limits::epsilon()) && (deltax > std::numeric_limits::epsilon() || deltay > std::numeric_limits::epsilon() || deltaz > std::numeric_limits::epsilon())) + throw cet::exception("RANGE") << "Cannot specify both momentum and delta position, exiting."; }; void PhotonGun::produce(art::Event& event) { std::unique_ptr output(new GenParticleCollection); const CLHEP::Hep3Vector pos(x, y, z); - const CLHEP::Hep3Vector p(px, py, pz); + CLHEP::Hep3Vector p(px, py, pz); + if (std::abs(deltax) > std::numeric_limits::epsilon() || + std::abs(deltay) > std::numeric_limits::epsilon() || + std::abs(deltaz) > std::numeric_limits::epsilon()) { + std::cout << "PhotonGun: Using deltax, deltay, deltaz to calculate momentum." << std::endl; + const CLHEP::Hep3Vector delta(deltax, deltay, deltaz); + p = delta.unit() * E; + } CLHEP::HepLorentzVector mom(p, E); output->push_back(GenParticle(PDGCode::gamma, GenId::particleGun, pos, mom, 0.)); event.put(std::move(output)); diff --git a/STMMC/fcl/Absorber.fcl b/STMMC/fcl/Absorber.fcl index 14c1785f8d..654cdd1372 100644 --- a/STMMC/fcl/Absorber.fcl +++ b/STMMC/fcl/Absorber.fcl @@ -10,7 +10,7 @@ process_name: HPGeAbsorberSpectrumShift source : { module_type : EmptyEvent - maxEvents : @local::Efficiency.NPhotons + maxEvents : 10 # @local::Efficiency.NPhotons } services : { @@ -25,9 +25,9 @@ physics: { producers : { generate : { module_type : PhotonGun - x : @local::Efficiency.Absorber.x - y : @local::Efficiency.Absorber.y - z : @local::Efficiency.Absorber.z + x : @local::ComponentPositions.LaBr.x + y : @local::ComponentPositions.LaBr.y + z : @local::ComponentPositions.LaBr.z E : @local::Efficiency.PhotonEnergy } g4run : @local::g4run @@ -35,8 +35,7 @@ physics: { analyzers : { EDep : { - module_type : MakeVirtualDetectorTree - VirtualDetectorId : 101 + module_type : VirtualDetectorTree StepPointMCsTag: "g4run:virtualdetector" SimParticlemvTag: "g4run:" } diff --git a/STMMC/fcl/AbsorberFromST.fcl b/STMMC/fcl/AbsorberFromST.fcl new file mode 100644 index 0000000000..03026040df --- /dev/null +++ b/STMMC/fcl/AbsorberFromST.fcl @@ -0,0 +1,59 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/STMMC/fcl/prolog.fcl" + +# This module simulates photons being fired from the ST incident on the STM detectors to investigate the relative shift in spectrum +# Original author : Pawel Plesniak +# Note - edit the parameters in this prolog.fcl, aside from which detector is being used which is stored in the generate module + +process_name: HPGeAbsorberSpectrumShift +source : { + module_type : EmptyEvent + maxEvents : @local::Efficiency.NPhotons +} + +services : { + # @local::Services.Sim + @table::Services.Core + @table::Services.SimOnly + message: @local::default_message + GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } +} + +physics: { + producers : { + generate : { + module_type : PhotonGun + x : @local::Efficiency.ST.x + y : @local::Efficiency.ST.y + z : @local::Efficiency.ST.z + deltax : @local::Efficiency.FromSTToDet.LaBr.x + deltay : @local::Efficiency.FromSTToDet.LaBr.y + deltaz : @local::Efficiency.FromSTToDet.LaBr.z + E : @local::Efficiency.PhotonEnergy + } + g4run : @local::g4run + } + + analyzers : { + EDep : { + module_type : VirtualDetectorTree + StepPointMCsTag: "g4run:virtualdetector" + SimParticlemvTag: "g4run:" + } + } + + p1 : [ generate, g4run ] + trigger_paths : [p1] + o1 : [ EDep ] + end_paths: [ o1 ] +} + +physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" +physics.producers.g4run.SDConfig.enableSD: [STMDet, virtualdetector] +services.SeedService.baseSeed : 8 +services.SeedService.maxUniqueEngines : 20 +services.TFileService.fileName : @local::Efficiency.OutputFilename +physics.producers.g4run.debug.trackingVerbosityLevel : 1 +physics.producers.g4run.debug.steppingVerbosityLevel : 1 diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 1b56b4f7d1..93895d3d16 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -34,6 +34,12 @@ ComponentPositions : { # In beam order LaBr : { # aperture centre x : -3863.4 y : 0.0 + z : 40404 + } + ST : { # ST centre + x : -3904.0 + y : 0.0 + z : 627.0 } } @@ -191,8 +197,37 @@ HPGeDigitization : { # Detector efficiency simulation Efficiency : { NPhotons : 1e7 - PhotonEnergy : 1.809 # MeV - OutputFilename : "Efficiency.root" + PhotonEnergy : @nil #1.809 # MeV + OutputFilename : "nts.owner.Absorber.Spectrum.sequencer.root" + ST : { + x : -3904.0 # mm + y : 0.0 # mm + z : 627.0 # mm + } + FromSTToSSCAperture: { + HPGe: { + x : -40.6 # mm + y : 0.0 # mm + z : 39777 # mm + } + LaBr: { + x : 40.6 # mm + y : 0.0 # mm + z : 39777 # mm + } + } + FromSTToDet: { # ST at (-3904, 0, 627) + HPGe: { + x : -40.6 # mm + y : 0.0 # mm + z : 40072.1 # mm + } + LaBr: { + x : 40.6 # mm + y : 0.0 # mm + z : 40264.5 # mm + } + } } END_PROLOG From 8b263a11dc95521cf5bb480de2bd8dd5dc366ce4 Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Wed, 30 Jul 2025 04:14:03 -0500 Subject: [PATCH 2/8] Updating the STM code to include momentum in the virtual detector trees --- EventGenerator/src/PhotonGun_module.cc | 9 ++- STMMC/fcl/Absorber.fcl | 11 ++-- STMMC/fcl/AbsorberFromSTHPGe.fcl | 60 +++++++++++++++++++ ...orberFromST.fcl => AbsorberFromSTLaBr.fcl} | 5 +- STMMC/fcl/MakeTree.fcl | 7 ++- STMMC/fcl/prolog.fcl | 2 +- STMMC/src/VirtualDetectorTree_module.cc | 48 +++++++++++---- 7 files changed, 118 insertions(+), 24 deletions(-) create mode 100644 STMMC/fcl/AbsorberFromSTHPGe.fcl rename STMMC/fcl/{AbsorberFromST.fcl => AbsorberFromSTLaBr.fcl} (92%) diff --git a/EventGenerator/src/PhotonGun_module.cc b/EventGenerator/src/PhotonGun_module.cc index a20da245be..202df23091 100644 --- a/EventGenerator/src/PhotonGun_module.cc +++ b/EventGenerator/src/PhotonGun_module.cc @@ -4,6 +4,7 @@ // stdlib includes #include +#include // art includes #include "art/Framework/Core/EDProducer.h" @@ -42,6 +43,7 @@ namespace mu2e { fhicl::OptionalAtom deltay{ Name("deltay"), Comment("Difference in y position of generated photon, use to calculate targeted position")}; fhicl::OptionalAtom deltaz{ Name("deltaz"), Comment("Difference in z position of generated photon, use to calculate targeted position")}; fhicl::Atom E{ Name("E"), Comment("Energy of generated photon")}; + fhicl::OptionalAtom verbose{ Name("verbose"), Comment("Print verbose messages")}; }; using Parameters = art::EDProducer::Table; explicit PhotonGun(const Parameters& conf); @@ -51,6 +53,7 @@ namespace mu2e { double px = 0.0, py = 0.0, pz = 0.0; double deltax = 0.0, deltay = 0.0, deltaz = 0.0; double E = 0.0; + bool verbose = false; }; PhotonGun::PhotonGun(const Parameters& conf): @@ -70,8 +73,9 @@ namespace mu2e { deltax = conf().deltax() ? *conf().deltax() : 0; deltay = conf().deltay() ? *conf().deltay() : 0; deltaz = conf().deltaz() ? *conf().deltaz() : 0; - if ((px > std::numeric_limits::epsilon() || py > std::numeric_limits::epsilon()) && (deltax > std::numeric_limits::epsilon() || deltay > std::numeric_limits::epsilon() || deltaz > std::numeric_limits::epsilon())) + if ((std::abs(px) > std::numeric_limits::epsilon() || std::abs(py) > std::numeric_limits::epsilon()) && (std::abs(deltax) > std::numeric_limits::epsilon() || std::abs(deltay) > std::numeric_limits::epsilon() || std::abs(deltaz) > std::numeric_limits::epsilon())) throw cet::exception("RANGE") << "Cannot specify both momentum and delta position, exiting."; + verbose = conf().verbose() ? *conf().verbose() : false; }; void PhotonGun::produce(art::Event& event) { @@ -81,7 +85,8 @@ namespace mu2e { if (std::abs(deltax) > std::numeric_limits::epsilon() || std::abs(deltay) > std::numeric_limits::epsilon() || std::abs(deltaz) > std::numeric_limits::epsilon()) { - std::cout << "PhotonGun: Using deltax, deltay, deltaz to calculate momentum." << std::endl; + if (verbose) + std::cout << "PhotonGun: Using delta position to calculate momentum." << std::endl; const CLHEP::Hep3Vector delta(deltax, deltay, deltaz); p = delta.unit() * E; } diff --git a/STMMC/fcl/Absorber.fcl b/STMMC/fcl/Absorber.fcl index 654cdd1372..bf07469418 100644 --- a/STMMC/fcl/Absorber.fcl +++ b/STMMC/fcl/Absorber.fcl @@ -10,7 +10,7 @@ process_name: HPGeAbsorberSpectrumShift source : { module_type : EmptyEvent - maxEvents : 10 # @local::Efficiency.NPhotons + maxEvents : @local::Efficiency.NPhotons } services : { @@ -25,10 +25,11 @@ physics: { producers : { generate : { module_type : PhotonGun - x : @local::ComponentPositions.LaBr.x - y : @local::ComponentPositions.LaBr.y - z : @local::ComponentPositions.LaBr.z + x : @local::ComponentPositions.Absorber.x + y : @local::ComponentPositions.Absorber.y + z : @local::ComponentPositions.Absorber.z E : @local::Efficiency.PhotonEnergy + verbose : false } g4run : @local::g4run } @@ -49,7 +50,7 @@ physics: { physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" physics.producers.g4run.SDConfig.enableSD: [STMDet, virtualdetector] -services.SeedService.baseSeed : 8 +services.SeedService.baseSeed : 18 services.SeedService.maxUniqueEngines : 20 services.TFileService.fileName : @local::Efficiency.OutputFilename # physics.producers.g4run.debug.trackingVerbosityLevel : 1 diff --git a/STMMC/fcl/AbsorberFromSTHPGe.fcl b/STMMC/fcl/AbsorberFromSTHPGe.fcl new file mode 100644 index 0000000000..db1cfd9f34 --- /dev/null +++ b/STMMC/fcl/AbsorberFromSTHPGe.fcl @@ -0,0 +1,60 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/STMMC/fcl/prolog.fcl" + +# This module simulates photons being fired from the ST incident on the STM detectors to investigate the relative shift in spectrum +# Original author : Pawel Plesniak +# Note - edit the parameters in this prolog.fcl, aside from which detector is being used which is stored in the generate module + +process_name: HPGeAbsorberSpectrumShift +source : { + module_type : EmptyEvent + maxEvents : @local::Efficiency.NPhotons +} + +services : { + # @local::Services.Sim + @table::Services.Core + @table::Services.SimOnly + message: @local::default_message + GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } +} + +physics: { + producers : { + generate : { + module_type : PhotonGun + x : @local::Efficiency.ST.x + y : @local::Efficiency.ST.y + z : @local::Efficiency.ST.z + deltax : @local::Efficiency.FromSTToDet.HPGe.x + deltay : @local::Efficiency.FromSTToDet.HPGe.y + deltaz : @local::Efficiency.FromSTToDet.HPGe.z + E : @local::Efficiency.PhotonEnergy + verbose: false + } + g4run : @local::g4run + } + + analyzers : { + EDep : { + module_type : VirtualDetectorTree + StepPointMCsTag: "g4run:virtualdetector" + SimParticlemvTag: "g4run:" + } + } + + p1 : [ generate, g4run ] + trigger_paths : [p1] + o1 : [ EDep ] + end_paths: [ o1 ] +} + +physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" +physics.producers.g4run.SDConfig.enableSD: [STMDet, virtualdetector] +services.SeedService.baseSeed : 8 +services.SeedService.maxUniqueEngines : 20 +services.TFileService.fileName : @local::Efficiency.OutputFilename +# physics.producers.g4run.debug.trackingVerbosityLevel : 1 +# physics.producers.g4run.debug.steppingVerbosityLevel : 1 diff --git a/STMMC/fcl/AbsorberFromST.fcl b/STMMC/fcl/AbsorberFromSTLaBr.fcl similarity index 92% rename from STMMC/fcl/AbsorberFromST.fcl rename to STMMC/fcl/AbsorberFromSTLaBr.fcl index 03026040df..a128c021ce 100644 --- a/STMMC/fcl/AbsorberFromST.fcl +++ b/STMMC/fcl/AbsorberFromSTLaBr.fcl @@ -32,6 +32,7 @@ physics: { deltay : @local::Efficiency.FromSTToDet.LaBr.y deltaz : @local::Efficiency.FromSTToDet.LaBr.z E : @local::Efficiency.PhotonEnergy + verbose: false } g4run : @local::g4run } @@ -55,5 +56,5 @@ physics.producers.g4run.SDConfig.enableSD: [STMDet, virtualdetector] services.SeedService.baseSeed : 8 services.SeedService.maxUniqueEngines : 20 services.TFileService.fileName : @local::Efficiency.OutputFilename -physics.producers.g4run.debug.trackingVerbosityLevel : 1 -physics.producers.g4run.debug.steppingVerbosityLevel : 1 +# physics.producers.g4run.debug.trackingVerbosityLevel : 1 +# physics.producers.g4run.debug.steppingVerbosityLevel : 1 diff --git a/STMMC/fcl/MakeTree.fcl b/STMMC/fcl/MakeTree.fcl index 1bcf442025..0861aa8b08 100644 --- a/STMMC/fcl/MakeTree.fcl +++ b/STMMC/fcl/MakeTree.fcl @@ -23,7 +23,7 @@ physics: { analyzers : { Stage1virtualdetector : { module_type : VirtualDetectorTree - StepPointMCsTag : @local::SimplifyStage1Data.StepPointMCsTag + StepPointMCsTag : @local::SimplifyStage1Data.StepPointMCsTag1809 SimParticlemvTag : @local::SimplifyStage1Data.SimParticlemvTag } Stage2virtualdetector : { @@ -49,8 +49,9 @@ physics: { EnergyCalib : @local::HPGeDigitization.EnergyPerADCBin } } - o1 : @nil + + o1 : [Stage1virtualdetector] end_paths: [o1] } -services.TFileService.fileName : @nil +services.TFileService.fileName : "test.root" diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 93895d3d16..8063b34a5b 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -13,7 +13,7 @@ ComponentPositions : { # In beam order } Absorber : { # 50mm x 50mm x 390mm x : -3944.6 - y : -24.9 # TODO - why is this not 0.0? + y : 0 z : 39900 } VD101 : { # Centre diff --git a/STMMC/src/VirtualDetectorTree_module.cc b/STMMC/src/VirtualDetectorTree_module.cc index 6156d1fc6b..6f6ed43048 100644 --- a/STMMC/src/VirtualDetectorTree_module.cc +++ b/STMMC/src/VirtualDetectorTree_module.cc @@ -19,6 +19,7 @@ // fhicl includes #include "canvas/Utilities/InputTag.h" #include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" // message handling #include "messagefacility/MessageLogger/MessageLogger.h" @@ -45,6 +46,7 @@ namespace mu2e { struct Config { fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; fhicl::Atom SimParticlemvTag{Name("SimParticlemvTag"), Comment("Tag identifying the SimParticlemv")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; }; using Parameters = art::EDAnalyzer::Table; explicit VirtualDetectorTree(const Parameters& conf); @@ -54,8 +56,8 @@ namespace mu2e { art::ProductToken StepPointMCsToken; art::ProductToken SimParticlemvToken; GlobalConstantsHandle pdt; - int pdgId = 0; - double x = 0.0, y = 0.0, z = 0.0, mass = 0.0, E = 0.0, time = 0.0; + int pdgId = 0, consecutiveEmptyFileCounter = 0, consecutiveEmptyFileThreshold = 0; + double x = 0.0, y = 0.0, z = 0.0, mass = 0.0, Ekin = 0.0, Etot = 0.0, time = 0.0, p = 0.0, p2 = 0.0, px = 0.0, py = 0.0, pz = 0.0; VolumeId_type virtualdetectorId = 0; TTree* ttree; std::map pdgIds; // @@ -65,6 +67,7 @@ namespace mu2e { art::EDAnalyzer(conf), StepPointMCsToken(consumes(conf().StepPointMCsTag())), SimParticlemvToken(consumes(conf().SimParticlemvTag())) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; art::ServiceHandle tfs; ttree = tfs->make( "ttree", "Virtual Detectors ttree"); ttree->Branch("time", &time, "time/D"); // ns @@ -73,18 +76,35 @@ namespace mu2e { ttree->Branch("x", &x, "x/D"); // mm ttree->Branch("y", &y, "y/D"); // mm ttree->Branch("z", &z, "z/D"); // mm - ttree->Branch("E", &E, "E/D"); // MeV + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("p", &p, "p/D"); // mm + ttree->Branch("p2", &p2, "p2/D"); // mm + ttree->Branch("mass", &mass, "mass/D"); // MeV/c^2 + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + ttree->Branch("Etot", &Etot, "Etot/D"); // MeV }; void VirtualDetectorTree::analyze(const art::Event& event) { - // Get the data products from the event - auto const& StepPointMCs = event.getProduct(StepPointMCsToken); - if (StepPointMCs.empty()) - return; - auto const& SimParticles = event.getProduct(SimParticlemvToken); - if (SimParticles.empty()) + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; return; + } + // Try to get SimParticles (assuming it's a map or vector) + auto simHandle = event.getHandle< SimParticleCollection >(SimParticlemvToken); + if (!simHandle || simHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + } + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, stopping the job"); + } + auto const& StepPointMCs = *stepHandle; + auto const& SimParticles = *simHandle; + consecutiveEmptyFileCounter = 0; // Loop over all VD hits for (const StepPointMC& step : StepPointMCs) { // Get the associated particle @@ -97,9 +117,15 @@ namespace mu2e { x = step.position().x(); y = step.position().y(); z = step.position().z(); + px = step.momentum().x(); + py = step.momentum().y(); + pz = step.momentum().z(); + p2 = step.momentum().mag2(); + p = std::sqrt(p2); mass = pdt->particle(pdgId).mass(); - E = std::sqrt(step.momentum().mag2()+mass*mass)-mass; // Subtract the rest mass - if (E < 0) + Etot = std::sqrt(p2 + mass * mass); // Total energy + Ekin = Etot - mass; // Subtract the rest mass + if (Ekin < 0) throw cet::exception("LogicError", "Energy is negative"); ttree->Fill(); From 1d11c8a8dc5f6bbcd5970ad6526a8987ad1eda17 Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Sun, 23 Nov 2025 18:26:48 -0600 Subject: [PATCH 3/8] Stage 1 simulation files corrected --- STMMC/fcl/prolog.fcl | 271 +++++++++++++++++++++++++++++-------------- 1 file changed, 181 insertions(+), 90 deletions(-) diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 8063b34a5b..cf1d75b682 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -1,10 +1,72 @@ -# -# prolog.fcl for fcl paramaters that will be used in STM simulation modules -# +# Defines standard paramaters used in STM simulation modules +# Adapted by: Pawel Plesniak + +# TODOs before Pawel uploads: +# - Review all parameters for correctness +# - Add comments where necessary +# - Verify that all necessary parameters are included +# - Simplify the ROOTAnalysisDump data product names by defining them separately here, once the simulation generats them correctly +# - SimParticlemv -> SimParticles +# - StepPointMCsTag -> StepPointMCs #include "Offline/STMReco/fcl/prolog.fcl" +#include "Production/JobConfig/pileup/STM/prolog.fcl" BEGIN_PROLOG +# Data product definitions +DataProducts : { + # Data products generated in Stage 1 of the simulation + Stage1 : { + StepPointMCs : @local::STMSimDataProducts.Stage1.CompressedOutput.StepPointMCs + SimParticles : @local::STMSimDataProducts.Stage1.CompressedOutput.SimParticles + EndVirtualdetectorID : @local::STMPileup.ResamplingProducer.VirtualDetectorID + } + + # Data products generated in Stage 2 of the simulation + Stage2 : { + StepPointMCs : { + Virtualdetector : @local::STMSimDataProducts.Stage2.CompressedOutput.StepPointMCs.Virtualdetector + Detector : @local::STMSimDataProducts.Stage2.CompressedOutput.StepPointMCs.Detector + } + SimParticles : @local::STMSimDataProducts.Stage2.CompressedOutput.SimParticles + EndVirtualdetectorID : { + HPGe : 90 + LaBr : 89 + } + } + + # Data products used in mixing from POT frame to microspill frame + MixingStepPointMCTags : { + StepPointMCsTagEle : "STMStepMixerEle:STMDet" + StepPointMCsTagMu : "STMStepMixerMu:STMDet" + StepPointMCsTag1809 : "STMStepMixer1809:STMDet" + } + + # Data products generated in STM waveform generation and analysis + Waveforms : { + HPGe : { + uSpill : "DigiHPGe:" + concatenated : "concatenateWaveformsHPGe:" + ZS : "ZSHPGe:" + } + LaBr : { + uSpill : "DigiLaBr:" + concatenated : "concatenateWaveformsLaBr:" + ZS : "ZSLaBr:" + } + } + + # Data products generated in STM MWD analysis + PeakEnergies : { + HPGe : "MWDHPGe:" + LaBr : "MWDLaBr:" + } + + # Miscellaneous parameters + nMicrospillsPerSpill : 31858 +} + +## Simulation parameters # Geometry, all in mm ComponentPositions : { # In beam order Beam : { @@ -43,73 +105,30 @@ ComponentPositions : { # In beam order } } -# Stage 1 propagation parameters -ResamplingProducer : { - StepPointMCsTag : "g4run:virtualdetector" - VirtualDetectorID : 101 -} -ResamplingFilter : { - StepPointMCsTag : "compressDetStepMCsSTM" -} -VirtualDetectorCounter : { - StepPointMCsTag : "g4run:virtualdetector" - virtualDetectorIDs : [88, 89, 90, 100, 101] -} - -# ShiftVirtualDetectorStepPointMCs_module.cc parameters -ShiftVD101Steps : { - StepPointMCsTag : "compressDetStepMCsSTM:virtualdetector" - VirtualDetectorID : 101 - InputRadius : 200 # VD101 aperture radius - OutputRadius : 3.98942280401 # SSC aperture radius - HPGeUpStr : { # Position just upstream of absorber centered on the HPGe SSC aperture - x : -3944.6 - y : 0.0 - z : 39906.0 - } - LaBrUpStr : { # Position just upstream of SSC centered on the Labr SSC aperture - x : -3863.4 - y : 0.0 - z : 40300.9 - } - pdgID : 0 -} - -# Mixing parameters -MixSTMEvents : { - extendedMean2BB : 3.93e7 # Copied from Production/JobConfig/mixing/TwoBB.fcl - cutMax2BB : 2.36e8 # Copied from Production/JobConfig/mixing/TwoBB.fcl - extendedMean1BB : 1.58e7 # Copied from Production/JobConfig/mixing/OneBB.fcl - cutMax1BB : 9.48e7 # Copied from Production/JobConfig/mixing/OneBB.fcl - SDF : 0.6 - nPOTs : 26919873604557116 - nMicroSpills : 690253169 - meanEventsPerPOTFactors : { - EleBeamCat : 5.465876332164181e-11 - MuBeamCat : 7.887481312841522e-13 - TargetStopsCat1809 : 3.114418308224206e-16 - } -} +# Waveform generation and digitization parameters STMMCAnalysis : { MixedEventsTags : { - StepPointMCsTagEle : "STMStepMixerEle:STMDet" - StepPointMCsTagMu : "STMStepMixerMu:STMDet" - StepPointMCsTag1809 : "STMStepMixer1809:STMDet" + StepPointMCsTagEle : @local::DataProducts.MixingStepPointMCTags.StepPointMCsTagEle + StepPointMCsTagMu : @local::DataProducts.MixingStepPointMCTags.StepPointMCsTagMu + StepPointMCsTag1809 : @local::DataProducts.MixingStepPointMCTags.StepPointMCsTag1809 } ZS : { HPGe : { stmWaveformDigisTag : { - concatenated : "concatenateWaveformsHPGe:" - uSpill : "DigiHPGe:" + uSpill : @local::DataProducts.Waveforms.HPGe.uSpill + concatenated : @local::DataProducts.Waveforms.HPGe.concatenated } tbefore : @local::STM.HPGe.tbefore tafter : @local::STM.HPGe.tafter threshold : @local::STM.HPGe.threshold - window : @local::STM.HPGe.window # TODO - when using the optimized number from STMReco, the gradients were fluctuating between the minimum band of an int16_t and 0, but this disappeared when set to 50. Figure out why. + window : @local::STM.HPGe.window naverage : @local::STM.HPGe.naverage } LaBr : { - stmWaveformDigisTag : "DigiLaBr:" + stmWaveformDigisTag : { + uSpill : @local::DataProducts.Waveforms.LaBr.uSpill + concatenated : @local::DataProducts.Waveforms.LaBr.concatenated + } tbefore : @local::STM.LaBr.tbefore tafter : @local::STM.LaBr.tafter threshold : @local::STM.HPGe.threshold @@ -120,8 +139,8 @@ STMMCAnalysis : { MWD : { HPGe : { stmWaveformDigisTag : { - ZS: "ZSHPGe:" - concatenated: "concatenateWaveformsHPGe:" + concatenated: @local::DataProducts.Waveforms.HPGe.concatenated + ZS: @local::DataProducts.Waveforms.HPGe.ZS } tau : @local::STM.HPGe.tau M : @local::STM.HPGe.M @@ -136,40 +155,96 @@ STMMCAnalysis : { full: @local::STM.HPGe.defaultBaselineSD suppressed: 0 } - STMMWDDigiTag : "MWDHPGe:" + STMMWDDigiTag : @local::DataProducts.PeakEnergies.HPGe } LaBr : { - stmWaveformDigisTag : "ZSLaBr:" + stmWaveformDigisTag : { + concatenated: @local::DataProducts.Waveforms.LaBr.concatenated + ZS: @local::DataProducts.Waveforms.LaBr.ZS + } tau : @local::STM.HPGe.tau M : @local::STM.HPGe.M L : @local::STM.HPGe.L nsigma_cut : @local::STM.HPGe.nsigma_cut thresholdgrad : @local::STM.HPGe.thresholdgrad - defaultBaselineMean : @local::STM.HPGe.defaultBaselineMean - defaultBaselineSD : @local::STM.HPGe.defaultBaselineSD - STMMWDDigiTag : "MWDLaBr:" + defaultBaselineMean : { # TODO - this should be the pedestal + full : @local::STM.HPGe.defaultBaselineMean + suppressed : 0 + } + defaultBaselineSD : { # TODO - this should be NoiseSD converted to ADC bins + full: @local::STM.HPGe.defaultBaselineSD + suppressed: 0 + } + STMMWDDigiTag : @local::DataProducts.PeakEnergies.LaBr } } } -# MakeTree_module.cc parameters -SimplifyStage1Data : { - StepPointMCsTag : "compressDetStepMCsSTM:" - StepPointMCsTag1809 : "compressDetStepMCsSTM:virtualdetector" - SimParticlemvTag : "compressDetStepMCsSTM:" - VirtualDetectorID : 101 +# ROOTAnalysisDump.fcl parameters +ROOTAnalysisDump : { + Stage1 : { + StepPointMCsTag : @local::DataProducts.Stage1.StepPointMCs + SimParticlemvTag : @local::DataProducts.Stage1.SimParticles + VirtualDetectorID : @local::DataProducts.Stage1.EndVirtualdetectorID + } + Stage2 : { + StepPointMCsTagVirtualdetector : @local::DataProducts.Stage2.StepPointMCs.Virtualdetector + StepPointMCsTagSTMDet : @local::DataProducts.Stage2.StepPointMCs.Detector + SimParticlemvTag : @local::DataProducts.Stage2.SimParticles + VirtualDetectorID : { + HPGe : @local::DataProducts.Stage2.EndVirtualdetectorID.HPGe + LaBr : @local::DataProducts.Stage2.EndVirtualdetectorID.LaBr + } + DetectorName : { + HPGe : "HPGe" + LaBr : "LaBr" + } + } + HPGeAnalysis: { + MicrospillHPGeWaveformTag : @local::DataProducts.Waveforms.HPGe.uSpill + ZSWaveformTag : @local::DataProducts.Waveforms.HPGe.ZS + MWDResultsTag : @local::DataProducts.PeakEnergies.HPGe + } + LaBrAnalysis: { + MicrospillHPGeWaveformTag : @local::DataProducts.Waveforms.LaBr.uSpill + ZSWaveformTag : @local::DataProducts.Waveforms.LaBr.ZS + MWDResultsTag : @local::DataProducts.PeakEnergies.LaBr + } + SimParticleVirtualDetectorBacktrace : { + StepPointMCsTags : { + Stage1 : @local::DataProducts.Stage1.StepPointMCs + Stage2 : @local::DataProducts.Stage2.StepPointMCs + } + StartVirtualDetectorID : { + Stage1 : @local::DataProducts.Stage1.EndVirtualdetectorID + Stage2 : { + HPGe : @local::DataProducts.Stage2.EndVirtualdetectorID.HPGe + LaBr : @local::DataProducts.Stage2.EndVirtualdetectorID.LaBr + } + } + TraceVirtualdetectorIds : { + Stage1: [101, 100, 10, 9, 8] + Stage2: [90, 101, 100, 10, 9, 8] + } + } } -SimplifyStage2Data : { - StepPointMCsTagSTMDet : "compressSTMDet:STMDet" - StepPointMCsTagVirtualdetector : "compressSTMDet:virtualdetector" - SimParticlemvTag : "compressSTMDet:" - VirtualDetectorID : { - HPGe : 90 - LaBr : 89 - } - DetectorName : { - HPGe : "HPGe" - LaBr : "LaBr" + +# For analysis studies, you should only need to edit lines between the "edit lines" comments +######################################## Start of edit lines ######################################## + +# Mixing parameters +MixSTMEvents : { + extendedMean2BB : 3.93e7 # Copied from Production/JobConfig/mixing/TwoBB.fcl + cutMax2BB : 2.36e8 # Copied from Production/JobConfig/mixing/TwoBB.fcl + extendedMean1BB : 1.58e7 # Copied from Production/JobConfig/mixing/OneBB.fcl + cutMax1BB : 9.48e7 # Copied from Production/JobConfig/mixing/OneBB.fcl + SDF : 0.6 # Copied from Production/JobConfig/mixing/OneBB.fcl + nPOTs : 26919873604557116 # Effective number of POTs, correcting for the resampling factor + nMicroSpills : 690253169 + meanEventsPerPOTFactors : { + EleBeamCat : 5.465876332164181e-11 + MuBeamCat : 7.887481312841522e-13 + TargetStopsCat1809 : 3.114418308224206e-16 } } @@ -180,24 +255,39 @@ STMDAQParameters : { LaBr : 370.370370370 # MSPS - copied from Offline/STMConditions/fcl/prolog.fcl } } + +# HPGe waveform generation parameters HPGeDigitization : { - EnergyPerADCBin : 0.5 # keV/bin - PreamplifierNoiseSD : 0.32 # mV + EnergyPerADCBin : 1 # keV/bin + PreamplifierNoiseSD : 0.0 # mV DecayConstant : 50 # us - resetEventNumber : 31858 + resetEventNumber : @local::DataProducts.nMicrospillsPerSpill microspillBufferLengthCount : 1000 concatenation : { - STMWaveformDigisTag : "DigiHPGe:" - nMerge : 31858 - filterTag : "concatenateWaveformsHPGe:" + STMWaveformDigisTag : @local::DataProducts.Waveforms.HPGe.uSpill + nMerge : @local::DataProducts.nMicrospillsPerSpill + filterTag : @local::DataProducts.Waveforms.HPGe.ZS + } +} + +# LaBr waveform generation parameters - not implemented yet +LaBrDigitization : { + EnergyPerADCBin : @nil # keV/bin + PreamplifierNoiseSD : @nil # mV + DecayConstant : @nil # us + resetEventNumber : @nil + microspillBufferLengthCount : @nil + concatenation : { + STMWaveformDigisTag : @local::DataProducts.Waveforms.LaBr.uSpill + nMerge : @nil + filterTag : @local::DataProducts.Waveforms.LaBr.ZS } } -# Normalization parameter determination studies -# Detector efficiency simulation +# Detector efficiency simulation parameters Efficiency : { NPhotons : 1e7 - PhotonEnergy : @nil #1.809 # MeV + PhotonEnergy : @nil # MeV OutputFilename : "nts.owner.Absorber.Spectrum.sequencer.root" ST : { x : -3904.0 # mm @@ -229,5 +319,6 @@ Efficiency : { } } } +######################################### End of edit lines ######################################### END_PROLOG From cdad575e1a528018750d0d67c70bbf612801e403 Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Thu, 15 Jan 2026 08:22:46 -0600 Subject: [PATCH 4/8] First push of a very large update --- .../src/CountVirtualDetectorHits_module.cc | 13 +- RecoDataProducts/inc/STMWaveformDigi.hh | 12 +- STMMC/CMakeLists.txt | 16 ++ ... => HPGeWaveformGenerationAndAnalysis.fcl} | 37 ++- STMMC/fcl/MakeTree.fcl | 57 ---- STMMC/fcl/ROOTAnalysisDump.fcl | 107 +++++++ STMMC/fcl/prolog.fcl | 4 +- .../HPGeWaveformsFromStepPointMCs_module.cc | 143 ++++++++-- STMMC/src/STMResamplingProducer_module.cc | 4 +- STMMC/src/SimParticleAndVDBacktrace_module.cc | 263 ++++++++++++++++++ STMMC/src/SimParticleDump_module.cc | 182 ++++++++++++ STMMC/src/VirtualDetectorTree_module.cc | 5 +- 12 files changed, 734 insertions(+), 109 deletions(-) rename STMMC/fcl/{HPGeReco.fcl => HPGeWaveformGenerationAndAnalysis.fcl} (64%) delete mode 100644 STMMC/fcl/MakeTree.fcl create mode 100644 STMMC/fcl/ROOTAnalysisDump.fcl create mode 100644 STMMC/src/SimParticleAndVDBacktrace_module.cc create mode 100644 STMMC/src/SimParticleDump_module.cc diff --git a/Analyses/src/CountVirtualDetectorHits_module.cc b/Analyses/src/CountVirtualDetectorHits_module.cc index 656f5db9a3..2f80b510c6 100644 --- a/Analyses/src/CountVirtualDetectorHits_module.cc +++ b/Analyses/src/CountVirtualDetectorHits_module.cc @@ -51,10 +51,17 @@ namespace mu2e { : art::EDAnalyzer(conf), StepPointMCsToken(consumes(conf().stepPointMCsTag())), enabledVDs(conf().enabledVDs()) { - // Create list of unique enabled virtual detectors - std::set enabledVDsSet(enabledVDs.begin(), enabledVDs.end()); + // Create list of unique enabled virtual detectors, preserving the order + std::vector unqiueEnabledVDsVec;; + std::unordered_set uniqueEnabledVDsSet; + for (const int & enabledVD : enabledVDs) { + if (uniqueEnabledVDsSet.find(enabledVD) == uniqueEnabledVDsSet.end()) { + uniqueEnabledVDsSet.insert(enabledVD); + unqiueEnabledVDsVec.push_back(enabledVD); + }; + }; enabledVDs.clear(); - enabledVDs.insert(enabledVDs.end(), enabledVDsSet.begin(), enabledVDsSet.end()); + enabledVDs = unqiueEnabledVDsVec; // Insert _enabledVDs.size() zeros into the counter vector nVDs = enabledVDs.size(); diff --git a/RecoDataProducts/inc/STMWaveformDigi.hh b/RecoDataProducts/inc/STMWaveformDigi.hh index 276b923df4..e212e4624e 100644 --- a/RecoDataProducts/inc/STMWaveformDigi.hh +++ b/RecoDataProducts/inc/STMWaveformDigi.hh @@ -18,13 +18,15 @@ namespace mu2e { class STMWaveformDigi { public: // Initialise all variables - STMWaveformDigi() : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(0), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(std::vector()) {}; + STMWaveformDigi() : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(0), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(std::vector()), _tmp_sim_time(0) {}; // Constructor for timing plus trig offset - STMWaveformDigi(int16_t DetID, uint64_t EWT, uint64_t DTCtime, uint64_t ADCtime, uint32_t trigTimeOffset, std::vector &adcs) : _DetID(DetID), _EWT(EWT), _DTCtime(DTCtime), _ADCtime(ADCtime), _trigTimeOffset(trigTimeOffset), _adcs(adcs) {}; + STMWaveformDigi(int16_t DetID, uint64_t EWT, uint64_t DTCtime, uint64_t ADCtime, uint32_t trigTimeOffset, std::vector &adcs) : _DetID(DetID), _EWT(EWT), _DTCtime(DTCtime), _ADCtime(ADCtime), _trigTimeOffset(trigTimeOffset), _adcs(adcs), _tmp_sim_time(0) {}; // Constructor for peak fitting - STMWaveformDigi(int16_t DetID, uint64_t EWT, uint64_t DTCtime, uint64_t ADCtime, uint32_t trigTimeOffset, double peak_fitTime1, double peak_fitTime2, double peak_sep, std::vector &adcs) : _DetID(DetID), _EWT(EWT), _DTCtime(DTCtime), _ADCtime(ADCtime), _trigTimeOffset(trigTimeOffset), _peak_fitTime1(peak_fitTime1), _peak_fitTime2(peak_fitTime2), _peak_sep(peak_sep), _adcs(adcs) {}; + STMWaveformDigi(int16_t DetID, uint64_t EWT, uint64_t DTCtime, uint64_t ADCtime, uint32_t trigTimeOffset, double peak_fitTime1, double peak_fitTime2, double peak_sep, std::vector &adcs) : _DetID(DetID), _EWT(EWT), _DTCtime(DTCtime), _ADCtime(ADCtime), _trigTimeOffset(trigTimeOffset), _peak_fitTime1(peak_fitTime1), _peak_fitTime2(peak_fitTime2), _peak_sep(peak_sep), _adcs(adcs), _tmp_sim_time(0) {}; // Basic constructor - STMWaveformDigi(uint32_t trigTimeOffset, std::vector &adcs) : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(trigTimeOffset), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(adcs) {}; + STMWaveformDigi(uint32_t trigTimeOffset, std::vector &adcs) : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(trigTimeOffset), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(adcs), _tmp_sim_time(0) {}; + // Constructor for simulation validation + STMWaveformDigi(double sim_time, std::vector &adcs) : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(0), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(adcs), _tmp_sim_time(sim_time) {}; int16_t DetID () const { return _DetID; } uint64_t EWT () const { return _EWT; } @@ -35,6 +37,7 @@ namespace mu2e { double peak_fitTime2 () const { return _peak_fitTime2; } double peak_sep () const { return _peak_sep; } const std::vector& adcs () const { return _adcs; } + double tmp_sim_time() const { return _tmp_sim_time; } private: int16_t _DetID; @@ -46,6 +49,7 @@ namespace mu2e { double _peak_fitTime2; // fit time of second rising edge (ns) double _peak_sep; // separation time (ns) std::vector _adcs; // vector of ADC values for the waveform + double _tmp_sim_time; // temporary storage of sim time for validation (ns) }; typedef std::vector STMWaveformDigiCollection; diff --git a/STMMC/CMakeLists.txt b/STMMC/CMakeLists.txt index 7e719143d7..1f9f7d2768 100644 --- a/STMMC/CMakeLists.txt +++ b/STMMC/CMakeLists.txt @@ -57,6 +57,22 @@ cet_build_plugin(ShiftVirtualDetectorStepPointMCs art::module Offline::MCDataProducts ) +cet_build_plugin(SimParticleDump art::module + REG_SOURCE src/SimParticleDump_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::GlobalConstantsService + Offline::MCDataProducts +) + +cet_build_plugin(SimParticleAndVDBacktrace art::module + REG_SOURCE src/SimParticleAndVDBacktrace_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::GlobalConstantsService + Offline::MCDataProducts +) + cet_build_plugin(STMResamplingProducer art::module REG_SOURCE src/STMResamplingProducer_module.cc LIBRARIES REG diff --git a/STMMC/fcl/HPGeReco.fcl b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl similarity index 64% rename from STMMC/fcl/HPGeReco.fcl rename to STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl index b6dc9d5362..3c0058ef24 100644 --- a/STMMC/fcl/HPGeReco.fcl +++ b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl @@ -1,4 +1,13 @@ - +# Executes the full HPGe waveform generation and analysis chain on STMMC data +# Outline: +# 1. DigiHPGe - generates HPGe waveforms from StepPointMCs +# 2. concatenateWaveformsHPGe - concatenates microspill waveforms into spill waveforms +# 3. (optional) ZSHPGe - zero-suppresses the concatenated waveforms +# 4. MWDHPGe - performs moving window deconvolution on the concatenated (or zero-suppressed) waveforms +# 5. Ou1tputs the MWD digis to a ROOT file +# Note - each module is documented in the _module.cc file. This is the clearest way to understand what each module does and what parameters it takes. +# Usage: +# 1. Select which data you want to run with or without zero suppression by replacing "@nil" in trigger_paths (keeping the [BRAKCETS]!) with the single analyer you want to run, e.g. trigger_paths: [digitization_path_ZS] # Original author: Pawel Plesniak #include "Offline/fcl/standardServices.fcl" @@ -17,6 +26,7 @@ services : @local::Services.Sim physics : { producers : { DigiHPGe : { + # Digitize the StepPointMCs from the HPGe crystals into waveforms module_type : HPGeWaveformsFromStepPointMCs StepPointMCsTagEle : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagEle StepPointMCsTagMu : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagMu @@ -27,17 +37,20 @@ physics : { risingEdgeDecayConstant : @local::HPGeDigitization.DecayConstant microspillBufferLengthCount : @local::HPGeDigitization.microspillBufferLengthCount makeTTree : false + makeADCPlot : false timeOffset : 0 resetEventNumber : @local::HPGeDigitization.resetEventNumber verbosityLevel : 0 } concatenateWaveformsHPGe : { + # Concatenate the microspill waveforms into macrospill waveforms module_type : ConcatenateDigitizedWaveforms STMWaveformDigisTag : @local::HPGeDigitization.concatenation.STMWaveformDigisTag nMerge : @local::HPGeDigitization.concatenation.nMerge makeTTree : false } ZSHPGe : { + # Apply zero-suppression to the concatenated waveforms module_type : STMZeroSuppression stmWaveformDigisTag : @local::STMMCAnalysis.ZS.HPGe.stmWaveformDigisTag.concatenated tbefore : @local::STMMCAnalysis.ZS.HPGe.tbefore @@ -50,6 +63,7 @@ physics : { makeTTreeWaveforms: false } MWDHPGe : { + # Perform moving window deconvolution on the concatenated (or zero-suppressed) waveforms to determine the pulse energies module_type : STMMovingWindowDeconvolution stmWaveformDigisTag : @local::STMMCAnalysis.MWD.HPGe.stmWaveformDigisTag.concatenated tau : @local::STMMCAnalysis.MWD.HPGe.tau @@ -60,7 +74,7 @@ physics : { defaultBaselineMean : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineMean.suppressed defaultBaselineSD : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineSD.suppressed makeTTreeMWD: false - makeTTreeEnergies: false + makeTTreeEnergies: true TTreeEnergyCalib : @local::HPGeDigitization.EnergyPerADCBin verbosityLevel : 0 xAxis : "" @@ -72,21 +86,22 @@ physics : { STMWaveformDigisTag : @local::HPGeDigitization.concatenation.filterTag } } - digitization_path : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] - trigger_paths : [digitization_path] - output_path : [compressedOutput] - end_paths : [output_path] + digitization_path_no_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, MWDHPGe] + digitization_path_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] + trigger_paths : ["digitization_path_ZS"] # Populate me! TODO - make me @nil before git push + # output_path : [compressedOutput] + # end_paths : [output_path] } outputs : { compressedOutput : { module_type : RootOutput fileName : "dts.owner.HPGeReco.version.sequencer.art" - SelectEvents: ["digitization_path"] - outputCommands: [ - "drop *_*_*_*", - "keep mu2e::STMMWDDigis_MWDHPGe_*_HPGeReco" - ] + SelectEvents: [@nil] # Populate me! + # outputCommands: [ + # "drop *_*_*_*", + # "keep mu2e::STMMWDDigis_MWDHPGe_*_HPGeReco" + # ] } } diff --git a/STMMC/fcl/MakeTree.fcl b/STMMC/fcl/MakeTree.fcl deleted file mode 100644 index 0861aa8b08..0000000000 --- a/STMMC/fcl/MakeTree.fcl +++ /dev/null @@ -1,57 +0,0 @@ -# Extracts data from STM propagation scripts into ROOT TTrees -# Edit services.TFileService.fileName before running! -# Run as: $ mu2e -c MakeTree.fcl -S FileWithListOfDataFiles.txt -# -# Original author: Pawel Plesniak - -#include "Offline/STMMC/fcl/prolog.fcl" -#include "Offline/fcl/standardServices.fcl" - -process_name: MakeTree - -source : { - module_type : RootInput - fileNames: @nil -} -services : { - message : @local::default_message - GlobalConstantsService : { - inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" - } -} -physics: { - analyzers : { - Stage1virtualdetector : { - module_type : VirtualDetectorTree - StepPointMCsTag : @local::SimplifyStage1Data.StepPointMCsTag1809 - SimParticlemvTag : @local::SimplifyStage1Data.SimParticlemvTag - } - Stage2virtualdetector : { - module_type : VirtualDetectorTree - StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagVirtualdetector - SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag - } - Stage2HPGe : { - module_type : HPGeTree - StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagSTMDet - SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag - Detector: @local::SimplifyStage2Data.DetectorName.HPGe - } - Stage2LaBr : { - module_type : HPGeTree - StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagSTMDet - SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag - Detector: @local::SimplifyStage2Data.DetectorName.LaBr - } - MWDSpectra : { - module_type : MWDTree - STMMWDDigiTag : @local::STMMCAnalysis.MWD.HPGe.STMMWDDigiTag - EnergyCalib : @local::HPGeDigitization.EnergyPerADCBin - } - } - - o1 : [Stage1virtualdetector] - end_paths: [o1] -} - -services.TFileService.fileName : "test.root" diff --git a/STMMC/fcl/ROOTAnalysisDump.fcl b/STMMC/fcl/ROOTAnalysisDump.fcl new file mode 100644 index 0000000000..41ffa70db8 --- /dev/null +++ b/STMMC/fcl/ROOTAnalysisDump.fcl @@ -0,0 +1,107 @@ +# Extracts parameters from STM simulation scripts into ROOT TTrees for use in the scripts defined it Mu2e/STMAnalysis repository +# This generates ROOT files containing different things depending on which analyzer is selected in the o1 path: +# - stage1virtualdetectorStepPointMCs : Dumps StepPointMCs from the virtual detector in Stage 1 simulation +# - stage2virtualdetectorStepPointMCs : Dumps StepPointMCs from the virtual detector in Stage 2 simulation +# - HPGeWaveforms : Dumps HPGe waveforms from Stage 2 simulation +# - LaBrWaveforms : Dumps LaBr3 waveforms from Stage 2 simulation (not implemented yet) +# - ZSwaveforms : Dumps zero-suppressed waveforms from HPGe digitisers in Stage 2 simulation (not implemented yet) +# - MWDspectra : Dumps MWD spectra from HPGe digitisers in Stage 2 simulation +# Usage: +# 1. Select which data you want to dump to a ROOT TTree by replacing "@nil" in the list o1 (keeping the [BRAKCETS]!) with the single analyer you want to run, e.g. o1: [stage1virtualdetectorStepPointMCs] +# 2. Set the output ROOT file name in services.TFileService.fileName +# 3. Prepare a text file containing a list of input STMMC data files (one file name per line) +# 3. Run as: `mu2e -c ROOTAnalysisDump.fcl -S FileWithListOfDataFiles.txt` +# +# Original author: Pawel Plesniak + +#include "Offline/STMMC/fcl/prolog.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name: ROOTAnalysisDump + +source : { + module_type : RootInput + fileNames: @nil +} +services : { + message : @local::default_message + GlobalConstantsService : { + inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" + } +} +physics: { + analyzers : { + # Simulation stage 1 ROOT files + Stage1virtualdetectorStepPointMCs : { + module_type : VirtualDetectorTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage1.StepPointMCsTag # 1809 + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + } + Stage1SimParticeandVirtualDetectorBacktrace : { + module_type : SimParticleAndVDBacktrace + StepPointMCsTag : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + StartVirtualDetectorId : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + TraceVirtualdetectorIds : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + } + + # Simulation stage 2 ROOT files + Stage2virtualdetectorStepPointMCs : { + module_type : VirtualDetectorTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage2.StepPointMCsTagVirtualdetector + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + } + Stage2SimParticeandVirtualDetectorBacktrace : { + module_type : SimParticleAndVDBacktrace + StepPointMCsTag : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + StartVirtualDetectorId : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + TraceVirtualdetectorIds : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + } + # TODO - see if this can be removed. + # ROOTDump : { + # module_type : SignalParticleBacktraceVDs + # # StepPointMCsTag : "compressDetStepMCsSTM:virtualdetector" #Stage 1 + # # # StepPointMCsTag: "compressDetStepMCsSTM:" #Stage 1 1809 as Stage 1 Mu3 + # # StartVirtualDetectorId : 101 + # # TrackVirtualDetectorIDs : [101, 100, 10, 9, 8, 2, 1] + # StepPointMCsTag : "compressSTMDet:virtualdetector" #Stage 2 + # StartVirtualDetectorId : 90 + # TrackVirtualDetectorIDs : [90, 100, 10, 9, 8, 2, 1] + # EnergyWindowSize : 0.1 + # } + + # HPGe waveform and analysis ROOT files + HPGeWaveforms : { + module_type : HPGeTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage2.StepPointMCsTagSTMDet + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + Detector: @local::ROOTAnalysisDump.Stage2.DetectorName.HPGe + } + HPGeZSWaveforms : { + module_type : @nil # Doesn't exist yet! + StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagSTMDet + SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag # TODO: Remove the SimParticle dependency + Detector: @local::SimplifyStage2Data.DetectorName.HPGe + } + HPGeMWDspectra : { + module_type : MWDTree + STMMWDDigiTag : @local::STMMCAnalysis.MWD.HPGe.STMMWDDigiTag + EnergyCalib : @local::HPGeDigitization.EnergyPerADCBin # Once the proditions are implemneted, this should be removed. + } + + # LaBr3 waveform ROOT files - not implemented yet + LaBrWaveforms : { + module_type : @nil + } + LaBrZSWaveforms : { + module_type : @nil + } + LaBrMWDspectra : { + module_type : @nil + } + } + + o1 : [@nil] # Populate me! + end_paths: [o1] +} + +services.TFileService.fileName : @nil # Populate me! diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index cf1d75b682..9ba94f8829 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -266,7 +266,7 @@ HPGeDigitization : { concatenation : { STMWaveformDigisTag : @local::DataProducts.Waveforms.HPGe.uSpill nMerge : @local::DataProducts.nMicrospillsPerSpill - filterTag : @local::DataProducts.Waveforms.HPGe.ZS + filterTag : @local::DataProducts.Waveforms.HPGe.concatenated } } @@ -280,7 +280,7 @@ LaBrDigitization : { concatenation : { STMWaveformDigisTag : @local::DataProducts.Waveforms.LaBr.uSpill nMerge : @nil - filterTag : @local::DataProducts.Waveforms.LaBr.ZS + filterTag : @local::DataProducts.Waveforms.LaBr.concatenated } } diff --git a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc index f084862a4e..413151ad00 100644 --- a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc +++ b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc @@ -1,6 +1,11 @@ // Simulates the electronics response of the HPGe detector. Simulates the pulse height, decay tail, and ADC digitization. Generates one STMWaveformDigi per micropulse. // Model based heavily on example provided in docDb 43617 // See docDb 51487 for full documentation +// Remaining TODOs for future development: +// - Include the commented out includes correctly +// - Implement the makeADCPlot functionality correctly +// - Implement different sampling frequencies correctly +// - Remove the implementation of stepPositionTolerance, it is a hack to address a position resolution I did not have time to fully address // Original author: Pawel Plesniak // stdlib includes @@ -45,6 +50,7 @@ // ROOT includes #include "art_root_io/TFileService.h" #include "TTree.h" +#include "TGraph.h" namespace mu2e { @@ -62,6 +68,7 @@ namespace mu2e { fhicl::Atom risingEdgeDecayConstant{ Name("risingEdgeDecayConstant"), Comment("Rising edge decay time [us]")}; fhicl::OptionalAtom microspillBufferLengthCount{ Name("microspillBufferLengthCount"), Comment("Number of microspills to buffer ahead for, in number of microspills")}; fhicl::OptionalAtom makeTTree{ Name("makeTTree"), Comment("Controls whether to make the TTree with branches chargeCollected, chargeDecayed, ADC, eventId, time")}; + fhicl::OptionalAtom makeADCPlot{ Name("makeADCPlot"), Comment("Controls whether to make a plot of the generated ADC values")}; fhicl::OptionalAtom timeOffset{ Name("timeOffset"), Comment("For debugging, adds the named time offset in [ns], used for testing analysis algorithms")}; fhicl::OptionalAtom resetEventNumber{ Name("resetEventNumber"), Comment("Simulates the off-spill period by resetting the inter-event last decayed charge to zero")}; fhicl::OptionalAtom verbosityLevel{ Name("verbosityLevel"), Comment("Controls verbosity")}; @@ -75,6 +82,7 @@ namespace mu2e { void decayCharge(); void addNoise(); void digitize(); + void endJob(); // fhicl variables art::ProductToken StepPointMCsTokenEle, StepPointMCsTokenMu, StepPointMCsToken1809; // Token of StepPointMCs in STMDet @@ -83,6 +91,7 @@ namespace mu2e { double noiseSD = 0; // Standard deviation of ADC noise [mV] double risingEdgeDecayConstant = 0; // [us] bool makeTTree = false; // Controls whether an analysis TTree is made + bool makeADCPlot = false; // Controls whether the ADC graph is made double timeOffset = 0.0; // Used for debugging [ns] uint resetEventNumber = 0; // Event ID at which to reset the charge amplitude int verbosityLevel = 0; // How much output to generate @@ -155,6 +164,10 @@ namespace mu2e { std::vector _chargeCarryOver; // Temporary buffer that will store _chargeCollected over the course of the next event std::vector _adcs; // Buffer for storing the ADC values to put into the STMWaveformDigi + // Debugging variables + double waveform_time = 0.0; // Time of start of waveform for debugging. TODO - REMOVE ME! + int32_t kept_events = 0, dropped_events = 0; // Counters for debugging + // Offline utilities // TODO: include the prodition to get the sampling frequency // mu2e::STMChannel::enum_type _HPGeChannel = static_cast(1); @@ -205,7 +218,7 @@ namespace mu2e { // Assign the approrpiate time offset timeOffset = conf().timeOffset() ? *(conf().timeOffset()) : 0.0; - // Assign TTrees + // Assign optional ROOT parameters makeTTree = conf().makeTTree() ? *(conf().makeTTree()) : false; if (makeTTree) { art::ServiceHandle tfs; @@ -216,7 +229,15 @@ namespace mu2e { ttree->Branch("eventId", &eventId, "eventId/i"); ttree->Branch("time", &time, "time/i"); }; + makeADCPlot = conf().makeADCPlot() ? *(conf().makeADCPlot()) : false; + if (makeADCPlot) { + throw cet::exception("IMPLEMENT", "makeADCPlot functionality not yet implemented correctly\n"); + art::ServiceHandle tfs; + TTree* adcTree = tfs->make("adcTree", "HPGeWaveformsFromStepPointMCs ADC Tree"); + adcTree->Branch("ADC", &ADC, "ADC/D"); + }; + // Assign optional variables resetEventNumber = conf().resetEventNumber() ? *(conf().resetEventNumber()) : 0; }; @@ -230,6 +251,7 @@ namespace mu2e { std::cout << std::left << "\t\t" << std::setw(60) << "risingEdgeDecayConstant [us]" << risingEdgeDecayConstant << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "microspillBufferLengthCount" << microspillBufferLengthCount << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "makeTTree" << makeTTree << std::endl; + std::cout << std::left << "\t\t" << std::setw(60) << "makeADCPlot" << makeADCPlot << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "timeOffset [ns]" << timeOffset << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "resetEventNumber" << resetEventNumber << std::endl; std::cout << "\tDerived parameters: " << std::endl; @@ -249,36 +271,92 @@ namespace mu2e { void HPGeWaveformsFromStepPointMCs::produce(art::Event& event) { eventId = event.id().event(); + // Simulation takes the POT time as t = 0, and has sequential microspills (events). The trigger time offset is not used here, left as a TODO + // Create the STMWaveformDigi and insert all the relevant attributes + // TODO - this only keeps accurate time if the sampling is 320MHz. Needs to be rewritten to work for other times + if (std::abs(fADC - 320) > 1e-12) + throw cet::exception("RANGE") << "Currently only fADC of 320 MHz is supported (got " << fADC << " MHz)\n"; + nADCs = 542; + + // Assign the variables for time and number of ADCs + eventTimeBuffer = eventId % 5; + if (!(eventTimeBuffer == 0 || eventTimeBuffer == 3)) + nADCs++; + eventTime += nADCs; + + // Update the last event decayed charge before the noise so the noise isn't added twice + if (resetEventNumber != 0 && eventId == resetEventNumber) { + lastEventEndDecayedCharge = 0; + std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); + } + else + lastEventEndDecayedCharge = _chargeDecayed.back(); + + // Update the parameters to carry over to the next event + _chargeCollected.clear(); + _chargeCollected.assign(_chargeCarryOver.begin(), _chargeCarryOver.end()); + _chargeCollected.insert(_chargeCollected.end(), nADCs, 0); + _chargeCarryOver.clear(); + _chargeCarryOver.assign(_chargeCollected.begin() + nADCs, _chargeCollected.end()); + + // Clear previous buffer vectors + _chargeDecayed.clear(); + _chargeDecayed.insert(_chargeDecayed.end(), nADCs, 0); + _adcs.clear(); + _adcs.insert(_adcs.end(), nADCs, 0); + + // Reset waveform time for debugging + waveform_time = 0.0; + // Get the hits in the detector std::vector StepsEle = event.getProduct(StepPointMCsTokenEle); std::vector StepsMu = event.getProduct(StepPointMCsTokenMu); std::vector Steps1809 = event.getProduct(StepPointMCsToken1809); + // If there are no hits, produce an empty STMWaveformDigiCollection and return + if (StepsEle.size() == 0 && StepsMu.size() == 0 && Steps1809.size() == 0) { + std::vector emptyAdcs(nADCs, static_cast(0)); + STMWaveformDigi _waveformDigi(0.0, emptyAdcs); + std::unique_ptr outputDigis(new STMWaveformDigiCollection); + outputDigis->emplace_back(_waveformDigi); + event.put(std::move(outputDigis)); + dropped_events++; + return; + } + kept_events++; + // Add a collection of charge depositions to _charge for(const StepPointMC& step : StepsEle){ - if (step.ionizingEdep() != 0) - depositCharge(step); + if (step.ionizingEdep() != 0) { + depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); + }; }; + // Process the collected charge for each step + for (const StepPointMC& step : StepsEle) { + if (step.ionizingEdep() != 0) { + depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); + } + } for(const StepPointMC& step : StepsMu){ if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); }; for(const StepPointMC& step : Steps1809){ if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); }; // Decay all of the collected charges decayCharge(); - // Update the last event decayed charge before the noise so the noise isn't added twice - if (resetEventNumber != 0 && eventId == resetEventNumber) { - lastEventEndDecayedCharge = 0; - std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); - } - else - lastEventEndDecayedCharge = _chargeDecayed.back(); - // Add preamplifier electronics noise with SD defined in noiseSD addNoise(); @@ -291,15 +369,8 @@ namespace mu2e { throw cet::exception("LogicError", "ADC values too high!"); }; - // Simulation takes the POT time as t = 0, and has sequential microspills (events). The trigger time offset is not used here, left as a TODO - // Create the STMWaveformDigi and insert all the relevant attributes - // TODO - this only keeps accurate time if the sampling is 320MHz. Needs to be rewritten to work for other times - eventTimeBuffer = eventId % 5; - if (eventTimeBuffer == 0 || (eventId % 3) == 0) - eventTime += nADCs + 1; - else - eventTime += nADCs; - STMWaveformDigi _waveformDigi(eventTime, _adcs); + // STMWaveformDigi _waveformDigi(eventTime, _adcs); + STMWaveformDigi _waveformDigi(waveform_time, _adcs); std::unique_ptr outputDigis(new STMWaveformDigiCollection); outputDigis->emplace_back(_waveformDigi); @@ -315,15 +386,21 @@ namespace mu2e { }; }; - // Update the parameters to carry over to the next event - _chargeCarryOver.clear(); - _chargeCarryOver.assign(_chargeCollected.begin() + nADCs, _chargeCollected.end()); - _chargeCollected.clear(); - _chargeCollected.assign(_chargeCarryOver.begin(), _chargeCarryOver.end()); - _chargeCollected.insert(_chargeCollected.end(), nADCs, 0); - // Clear previous buffer vectors - std::fill(_chargeDecayed.begin(), _chargeDecayed.end(), 0); - std::fill(_adcs.begin(), _adcs.end(), 0); + // Make the ADC plot if appropriate + if (makeADCPlot) { + // NOTE - THIS DOES NOT CURRENTLY FUNCTION, AND I DO NOT HAVE TIME TO FIX IT RIGHT NOW + std::stringstream histsuffix; + histsuffix.str(""); + histsuffix << "_evt" << event.event(); + + art::ServiceHandle tfs; + std::vector _adc_double = std::vector(_adcs.begin(), _adcs.end()); + TGraph* g_adcs = tfs->make(nADCs, _adc_double.data()); + g_adcs->SetName(("g_adcs"+histsuffix.str()).c_str()); + // g_adcs->SetTitle("ADC Values"); + // for (uint i = 0; i < nADCs; i++) + // g_adcs->SetPoint(i, i, _adcs[i]); + }; // Add the STMWaveformDigi to the event event.put(std::move(outputDigis)); @@ -480,6 +557,14 @@ namespace mu2e { }; return; }; + + void HPGeWaveformsFromStepPointMCs::endJob() { + mf::LogInfo log("STMResamplingFilter summary"); + log << "=====HPGeWaveformsFromStepPointMCs summary=====\n"; + log << "No. kept events: " << kept_events << "\n"; + log << "No. discarded events: " << dropped_events << "\n"; + log << "===============================================\n"; + }; }; // namespace mu2e DEFINE_ART_MODULE(mu2e::HPGeWaveformsFromStepPointMCs) diff --git a/STMMC/src/STMResamplingProducer_module.cc b/STMMC/src/STMResamplingProducer_module.cc index 02c83d003f..30137f61df 100644 --- a/STMMC/src/STMResamplingProducer_module.cc +++ b/STMMC/src/STMResamplingProducer_module.cc @@ -44,7 +44,7 @@ namespace mu2e { art::EDProducer{conf}, StepPointMCsToken(consumes(conf().stepPointMCsTag())), virtualDetectorID(conf().virtualDetectorID()) { - produces(); + produces("virtualdetector"); }; void STMResamplingProducer::produce(art::Event& event) { @@ -63,7 +63,7 @@ namespace mu2e { // Update counter includedStepPointMCs += outputStepPointMCs->size(); // Add the new data products to the event - event.put(std::move(outputStepPointMCs)); + event.put(std::move(outputStepPointMCs), "virtualdetector"); return; }; diff --git a/STMMC/src/SimParticleAndVDBacktrace_module.cc b/STMMC/src/SimParticleAndVDBacktrace_module.cc new file mode 100644 index 0000000000..82f4538588 --- /dev/null +++ b/STMMC/src/SimParticleAndVDBacktrace_module.cc @@ -0,0 +1,263 @@ +// Adapted from ReadVirtualDetector_module.cc +// Generates a large ROOT TTree with SimParticle and StepPointMC information. For each SimParticle that has a StepPointMC in the specified virtual detector, gets its +// - parentCounter - 0 for the original particle, 1 for its parent, 2 for its grandparent, etc. This is traced all the way back to the POT +// - particleId - SimParticle ID (from GEANT) +// - pdgId - PDG ID +// - creationCode - creation code. See Offline/MCDataProducts/inc/ProcessCode.hh for details +// - x [mm] - start x position of the SimParticle +// - y [mm] - start y position of the SimParticle +// - z [mm] - start z position of the SimParticle +// - px [MeV/c] - start x momentum of the SimParticle +// - py [MeV/c] - start y momentum of the SimParticle +// - pz [MeV/c] - start z momentum of the SimParticle +// - time [ns] - start global time of the SimParticle +// - Ekin [MeV] - end kinetic energy of the SimParticle +// For each chosen virtual detector ID to trace back to, also gets the following StepPointMC information (0 if the SimParticle did not have a StepPointMC in that virtual detector): +// - VD__x [mm] - x position of the StepPointMC in that virtual detector +// - VD__y [mm] - y position of the StepPointMC in that virtual detector +// - VD__z [mm] - z position of the StepPointMC in that virtual detector +// - VD__px [MeV/c] - x momentum of the StepPointMC in that virtual detector +// - VD__py [MeV/c] - y momentum of the StepPointMC in that virtual detector +// - VD__pz [MeV/c] - z momentum of the StepPointMC in that virtual detector +// - VD__time [ns] - time of the StepPointMC in that virtual detector +// - VD__Ekin [MeV] - kinetic energy of the StepPointMC in that virtual detector +// - VD__pdgId - PDG ID of the SimParticle at that StepPointMC +// Original author: Ivan Logashenko +// Adapted by: Pawel Plesniak + +// stdlib includes +#include +#include + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" + +// exception handling +#include "cetlib_except/exception.h" + +// fhicl includes +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" + +// message handling +#include "messagefacility/MessageLogger/MessageLogger.h" + +// Offline includes +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/StepPointMC.hh" + +// ROOT includes +#include "art_root_io/TFileService.h" +#include "TTree.h" + +typedef cet::map_vector_key key_type; + +namespace mu2e { + class SimParticleAndVDBacktrace : public art::EDAnalyzer { + public: + using Name=fhicl::Name; + using Comment=fhicl::Comment; + struct Config { + fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; + fhicl::Atom StartVirtualDetectorId{Name("StartVirtualDetectorId"), Comment("ID of the virtual detector to filter")}; + fhicl::Sequence TraceVirtualdetectorIds{Name("TraceVirtualdetectorIds"), Comment("IDs of the virtual detectors to trace back to the origin")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit SimParticleAndVDBacktrace(const Parameters& conf); + void analyze(const art::Event& e); + void addToTree(const art::Ptr &particle, const std::vector steps); + void addParentToTree(const art::Ptr &particle, const std::vector steps); + void traceSteps(const art::Ptr &particle, const std::vector steps); + void clear_virtualdetector_vectors(); + void populate_virtualdetector_vector(const StepPointMC & step_it, const unsigned index); + private: + art::ProductToken StepPointMCsToken; + unsigned StartVirtualDetectorId = 0; + std::vector TraceVirtualdetectorIds; + int consecutiveEmptyFileThreshold = 0; + + GlobalConstantsHandle pdt; + int pdgId = 0, creationCode = 0, consecutiveEmptyFileCounter = 0, nTraceVirtualdetectorIds = 0; + double x = 0.0, y = 0.0, z = 0.0, px = 0.0, py = 0.0, pz = 0.0, time = 0.0, Ekin = 0.0; + std::vector VD_x, VD_y, VD_z, VD_px, VD_py, VD_pz, VD_time, VD_Ekin; + std::vector VD_pdgId; + uint16_t parentCounter = 0; + cet::map_vector_key particleId; + TTree* ttree; + CLHEP::Hep3Vector startPosition, startMomentum; + art::Ptr parent; + std::vector*> VD_vectors; + }; + + SimParticleAndVDBacktrace::SimParticleAndVDBacktrace(const Parameters& conf) : + art::EDAnalyzer(conf), + StepPointMCsToken(consumes(conf().StepPointMCsTag())), + StartVirtualDetectorId(conf().StartVirtualDetectorId()), + TraceVirtualdetectorIds(conf().TraceVirtualdetectorIds()) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; + if (std::find(TraceVirtualdetectorIds.begin(), TraceVirtualdetectorIds.end(), StartVirtualDetectorId) == TraceVirtualdetectorIds.end()) { + TraceVirtualdetectorIds.push_back(StartVirtualDetectorId); + }; + nTraceVirtualdetectorIds = TraceVirtualdetectorIds.size(); + VD_vectors = {&VD_x, &VD_y, &VD_z, &VD_px, &VD_py, &VD_pz, &VD_time, &VD_Ekin}; + for (auto vec : VD_vectors) { + vec->resize(nTraceVirtualdetectorIds, 0.0); + }; + VD_pdgId.resize(nTraceVirtualdetectorIds, 0); + + art::ServiceHandle tfs; + ttree = tfs->make( "ttree", "SimParticle ttree"); + ttree->Branch("parentCounter", &parentCounter, "parentCounter/s"); + ttree->Branch("particleId", &particleId, "particleId/l"); + ttree->Branch("pdgId", &pdgId, "pdgId/I"); + ttree->Branch("creationCode", &creationCode, "creationCode/I"); + ttree->Branch("x", &x, "x/D"); // mm + ttree->Branch("y", &y, "y/D"); // mm + ttree->Branch("z", &z, "z/D"); // mm + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("time", &time, "time/D"); // ns + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + for (size_t i=0; iBranch(Form("VD_%u_x", TraceVirtualdetectorIds[i]), &VD_x[i], Form("VD_%u_x/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_y", TraceVirtualdetectorIds[i]), &VD_y[i], Form("VD_%u_y/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_z", TraceVirtualdetectorIds[i]), &VD_z[i], Form("VD_%u_z/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_px", TraceVirtualdetectorIds[i]), &VD_px[i], Form("VD_%u_px/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_py", TraceVirtualdetectorIds[i]), &VD_py[i], Form("VD_%u_py/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_pz", TraceVirtualdetectorIds[i]), &VD_pz[i], Form("VD_%u_pz/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_time", TraceVirtualdetectorIds[i]), &VD_time[i], Form("VD_%u_time/D", TraceVirtualdetectorIds[i])); // ns + ttree->Branch(Form("VD_%u_Ekin", TraceVirtualdetectorIds[i]), &VD_Ekin[i], Form("VD_%u_Ekin/D", TraceVirtualdetectorIds[i])); // MeV + ttree->Branch(Form("VD_%u_pdgId", TraceVirtualdetectorIds[i]), &VD_pdgId[i], Form("VD_%u_pdgId/I", TraceVirtualdetectorIds[i])); + }; + }; + + void SimParticleAndVDBacktrace::populate_virtualdetector_vector(const StepPointMC & step_it, const unsigned index) { + const CLHEP::Hep3Vector momentum = step_it.momentum(); + const CLHEP::Hep3Vector position = step_it.position(); + const double mass = pdt->particle(step_it.simParticle()->pdgId()).mass(); + const double EKin = std::sqrt(momentum.mag2() + mass * mass) - mass; + VD_x[index] = position.x(); + VD_y[index] = position.y(); + VD_z[index] = position.z(); + VD_px[index] = momentum.x(); + VD_py[index] = momentum.y(); + VD_pz[index] = momentum.z(); + VD_time[index] = step_it.time(); + VD_Ekin[index] = EKin; + VD_pdgId[index] = step_it.simParticle()->pdgId(); + }; + + void SimParticleAndVDBacktrace::traceSteps(const art::Ptr &particle, const std::vector steps) { + unsigned step_virtualdetectorId = 0, index = 0; + for (const StepPointMC & step_it : steps) { + if (step_it.simParticle()->id() != particle->id()) continue; + step_virtualdetectorId = step_it.virtualDetectorId(); + auto it = std::find(TraceVirtualdetectorIds.begin(), TraceVirtualdetectorIds.end(), step_virtualdetectorId); + if (it == TraceVirtualdetectorIds.end()) continue; + index = std::distance(TraceVirtualdetectorIds.begin(), it); + populate_virtualdetector_vector(step_it, index); + }; + }; + + void SimParticleAndVDBacktrace::clear_virtualdetector_vectors() { + for (auto & vec : VD_vectors) { + std::fill(vec->begin(), vec->end(), 0.0); + }; + std::fill(VD_pdgId.begin(), VD_pdgId.end(), 0); + return; + }; + + void SimParticleAndVDBacktrace::addToTree(const art::Ptr &particle, const std::vector steps) { + // Assigne all the SimParticle variables to the tree variables + parentCounter = 0; + particleId = particle->id(); + pdgId = particle->pdgId(); + creationCode = particle->creationCode(); + startPosition = particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle->startGlobalTime(); + Ekin = particle->endKineticEnergy(); + + // Trace the StepPointMCs in the relevant virtualdetectors + clear_virtualdetector_vectors(); + traceSteps(particle, steps); + ttree->Fill(); + + // Now trace the parents + if (!particle->hasParent()) return; + + // Require separate definitions of addParentToTree because of data accessors + parent = particle->parent(); + addParentToTree(parent, steps); + + while (parent->hasParent()) { + parent = parent->parent(); + addParentToTree(parent, steps); + }; + }; + + void SimParticleAndVDBacktrace::addParentToTree(const art::Ptr &particle, const std::vector steps) { + // Increment the parent counter + parentCounter++; + + // Assign all the SimParticle variables to the tree variables + particleId = particle->id(); + pdgId = particle->pdgId(); + creationCode = particle->creationCode(); + startPosition = particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle->startGlobalTime(); + Ekin = particle->endKineticEnergy(); + + // Trace the StepPointMCs in the relevant virtualdetectors + clear_virtualdetector_vectors(); + traceSteps(particle, steps); + ttree->Fill(); + }; + + void SimParticleAndVDBacktrace::analyze(const art::Event& event) { + // Get the data products + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, are you sure you have the correct data product name?\n"); + }; + return; + }; + const std::vector StepPointMCs = *stepHandle; + consecutiveEmptyFileCounter = 0; + + for (const StepPointMC& step : StepPointMCs) { + // Filter by virtual detector ID if requested + if (step.virtualDetectorId() != StartVirtualDetectorId) continue; + + // Get the associated particle + const art::Ptr particle = step.simParticle(); + addToTree(particle, StepPointMCs); + }; + return; + }; +}; // end namespace mu2e + +DEFINE_ART_MODULE(mu2e::SimParticleAndVDBacktrace) diff --git a/STMMC/src/SimParticleDump_module.cc b/STMMC/src/SimParticleDump_module.cc new file mode 100644 index 0000000000..48cf595cf6 --- /dev/null +++ b/STMMC/src/SimParticleDump_module.cc @@ -0,0 +1,182 @@ +// Adapted from ReadVirtualDetector_module.cc +// For StepPointMCs in virtualdetectors, generates a TTree showing the provenance of the generated SimParticles with the following branches +// - pdgId - PDG ID +// - creationCode - creation code. See Offline/MCDataProducts/inc/ProcessCode.hh for details +// - x [mm] - start x position of the SimParticle +// - y [mm] - start y position of the SimParticle +// - z [mm] - start z position of the SimParticle +// - px [MeV/c] - start x momentum of the SimParticle +// - py [MeV/c] - start y momentum of the SimParticle +// - pz [MeV/c] - start z momentum of the SimParticle +// - time [ns] - start global time of the SimParticle +// - Ekin [MeV] - end kinetic energy of the SimParticle +// Original author: Ivan Logashenko +// Adapted by: Pawel Plesniak + +// TODO before upload - validate that the code works and generates the same particle trace as the SimParticleAndVDBacktrace module + +// stdlib includes +#include +#include + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" + +// exception handling +#include "cetlib_except/exception.h" + +// fhicl includes +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" + +// message handling +#include "messagefacility/MessageLogger/MessageLogger.h" + +// Offline includes +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/StepPointMC.hh" + +// ROOT includes +#include "art_root_io/TFileService.h" +#include "TTree.h" + +typedef cet::map_vector_key key_type; + +namespace mu2e { + class SimParticleDump : public art::EDAnalyzer { + public: + using Name=fhicl::Name; + using Comment=fhicl::Comment; + struct Config { + fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; + fhicl::Atom SimParticlemvTag{Name("SimParticlemvTag"), Comment("Tag identifying the SimParticlemv")}; + fhicl::Atom FilterVirtualDetectorId{Name("FilterVirtualDetectorId"), Comment("ID of the virtual detector to filter")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit SimParticleDump(const Parameters& conf); + void analyze(const art::Event& e); + void addToTree(const SimParticle& particle, TTree* ttree); + void addParentToTree(art::Ptr &_particle, TTree* ttree); + private: + art::ProductToken StepPointMCsToken; + art::ProductToken SimParticlemvToken; + GlobalConstantsHandle pdt; + int pdgId = 0, creationCode = 0, consecutiveEmptyFileCounter = 0, consecutiveEmptyFileThreshold = 0; + double x = 0.0, y = 0.0, z = 0.0, px = 0.0, py = 0.0, pz = 0.0, time = 0.0, Ekin = 0.0; + unsigned virtualdetectorId = 0, FilterVirtualDetectorId = 0; + bool forward = false; + uint16_t parentCounter = 0; + cet::map_vector_key particleId; + TTree* ttree; + CLHEP::Hep3Vector startPosition, startMomentum; + art::Ptr parent; + }; + + SimParticleDump::SimParticleDump(const Parameters& conf) : + art::EDAnalyzer(conf), + StepPointMCsToken(consumes(conf().StepPointMCsTag())), + SimParticlemvToken(consumes(conf().SimParticlemvTag())), + FilterVirtualDetectorId(conf().FilterVirtualDetectorId()) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; + art::ServiceHandle tfs; + ttree = tfs->make( "ttree", "SimParticle ttree"); + ttree->Branch("parentCounter", &parentCounter, "parentCounter/s"); + ttree->Branch("particleId", &particleId, "particleId/l"); + ttree->Branch("pdgId", &pdgId, "pdgId/I"); + ttree->Branch("creationCode", &creationCode, "creationCode/I"); + ttree->Branch("x", &x, "x/D"); // mm + ttree->Branch("y", &y, "y/D"); // mm + ttree->Branch("z", &z, "z/D"); // mm + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("time", &time, "time/D"); // ns + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + }; + + void SimParticleDump::addToTree(const SimParticle& particle, TTree* ttree) { + parentCounter = 0; + particleId = particle.id(); + pdgId = particle.pdgId(); + creationCode = particle.creationCode(); + startPosition = particle.startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle.startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle.startGlobalTime(); + Ekin = particle.endKineticEnergy(); + ttree->Fill(); + + if (!particle.hasParent()) { + return; + }; + parentCounter++; + parent = particle.parent(); + addParentToTree(parent, ttree); + while (parent->hasParent()) { + parentCounter++; + parent = parent->parent(); + addParentToTree(parent, ttree); + }; + }; + + void SimParticleDump::addParentToTree(art::Ptr &_particle, TTree* ttree) { + particleId = _particle->id(); + pdgId = _particle->pdgId(); + creationCode = _particle->creationCode(); + startPosition = _particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = _particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = _particle->startGlobalTime(); + Ekin = _particle->endKineticEnergy(); + ttree->Fill(); + }; + + void SimParticleDump::analyze(const art::Event& event) { + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + }; + auto simHandle = event.getHandle< SimParticleCollection >(SimParticlemvToken); + if (!simHandle || simHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + }; + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, stopping the job"); + }; + + auto const& StepPointMCs = *stepHandle; + auto const& SimParticles = *simHandle; + consecutiveEmptyFileCounter = 0; + + for (const StepPointMC& step : StepPointMCs) { + // Filter by virtual detector ID if requested + if (step.virtualDetectorId() != FilterVirtualDetectorId) continue; + + // Get the associated particle + const SimParticle& particle = SimParticles.at(step.trackId()); + addToTree(particle, ttree); + }; + return; + }; +}; // end namespace mu2e + +DEFINE_ART_MODULE(mu2e::SimParticleDump) diff --git a/STMMC/src/VirtualDetectorTree_module.cc b/STMMC/src/VirtualDetectorTree_module.cc index 6f6ed43048..a976156ce1 100644 --- a/STMMC/src/VirtualDetectorTree_module.cc +++ b/STMMC/src/VirtualDetectorTree_module.cc @@ -56,11 +56,12 @@ namespace mu2e { art::ProductToken StepPointMCsToken; art::ProductToken SimParticlemvToken; GlobalConstantsHandle pdt; - int pdgId = 0, consecutiveEmptyFileCounter = 0, consecutiveEmptyFileThreshold = 0; + int pdgId = 0, consecutiveEmptyFileCounter = 0, consecutiveEmptyFileThreshold = 0, simParticleId=0; double x = 0.0, y = 0.0, z = 0.0, mass = 0.0, Ekin = 0.0, Etot = 0.0, time = 0.0, p = 0.0, p2 = 0.0, px = 0.0, py = 0.0, pz = 0.0; VolumeId_type virtualdetectorId = 0; TTree* ttree; std::map pdgIds; // + uint simParticleIdKey = 0; }; VirtualDetectorTree::VirtualDetectorTree(const Parameters& conf) : @@ -84,6 +85,7 @@ namespace mu2e { ttree->Branch("mass", &mass, "mass/D"); // MeV/c^2 ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV ttree->Branch("Etot", &Etot, "Etot/D"); // MeV + ttree->Branch("SimParticleId", &simParticleId, "SimParticleId/i"); }; void VirtualDetectorTree::analyze(const art::Event& event) { @@ -125,6 +127,7 @@ namespace mu2e { mass = pdt->particle(pdgId).mass(); Etot = std::sqrt(p2 + mass * mass); // Total energy Ekin = Etot - mass; // Subtract the rest mass + simParticleId = particle.id().asInt(); if (Ekin < 0) throw cet::exception("LogicError", "Energy is negative"); ttree->Fill(); From 8597da075347b3ada7e49297bdcdcd64f8c062ab Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Sun, 15 Feb 2026 05:30:23 -0600 Subject: [PATCH 5/8] Formatting fixes --- STMMC/fcl/ROOTAnalysisDump.fcl | 24 ++++++++++++------------ STMMC/fcl/prolog.fcl | 2 +- STMMC/src/ConcatenationFilter_module.cc | 1 + 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/STMMC/fcl/ROOTAnalysisDump.fcl b/STMMC/fcl/ROOTAnalysisDump.fcl index 41ffa70db8..43e7e88c9b 100644 --- a/STMMC/fcl/ROOTAnalysisDump.fcl +++ b/STMMC/fcl/ROOTAnalysisDump.fcl @@ -69,19 +69,19 @@ physics: { # EnergyWindowSize : 0.1 # } - # HPGe waveform and analysis ROOT files - HPGeWaveforms : { + # HPGe energy depositions in crystal - MC truth + HPGeEDeps : { module_type : HPGeTree - StepPointMCsTag : @local::ROOTAnalysisDump.Stage2.StepPointMCsTagSTMDet - SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + StepPointMCsTag : "compressSTMDet:STMDet" + SimParticlemvTag : "compressSTMDet:" Detector: @local::ROOTAnalysisDump.Stage2.DetectorName.HPGe } - HPGeZSWaveforms : { - module_type : @nil # Doesn't exist yet! - StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagSTMDet - SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag # TODO: Remove the SimParticle dependency - Detector: @local::SimplifyStage2Data.DetectorName.HPGe - } +# HPGeZSWaveforms : { +# module_type : @nil # Doesn't exist yet! +# StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagSTMDet +# SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag # TODO: Remove the SimParticle dependency +# Detector: @local::SimplifyStage2Data.DetectorName.HPGe +# } HPGeMWDspectra : { module_type : MWDTree STMMWDDigiTag : @local::STMMCAnalysis.MWD.HPGe.STMMWDDigiTag @@ -100,8 +100,8 @@ physics: { } } - o1 : [@nil] # Populate me! + o1 : ["HPGeEDeps"] # Populate me! end_paths: [o1] } -services.TFileService.fileName : @nil # Populate me! +services.TFileService.fileName : "tree.root" # Populate me! diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 9ba94f8829..3fe4cefd30 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -211,7 +211,7 @@ ROOTAnalysisDump : { MWDResultsTag : @local::DataProducts.PeakEnergies.LaBr } SimParticleVirtualDetectorBacktrace : { - StepPointMCsTags : { + StepPointMCsTag : { Stage1 : @local::DataProducts.Stage1.StepPointMCs Stage2 : @local::DataProducts.Stage2.StepPointMCs } diff --git a/STMMC/src/ConcatenationFilter_module.cc b/STMMC/src/ConcatenationFilter_module.cc index ed360dee97..10d6712377 100644 --- a/STMMC/src/ConcatenationFilter_module.cc +++ b/STMMC/src/ConcatenationFilter_module.cc @@ -58,6 +58,7 @@ namespace mu2e { void ConcatenationFilter::endJob() { mf::LogInfo log("ConcatenationFilter summary"); + log << "\n"; log << "==========ConcatenationFilter summary==========\n"; log << std::left << std::setw(25) << "\tNo. input events:" << inputEvents << "\n"; log << std::left << std::setw(25) << "\tNo. kept events:" << keptEvents << "\n"; From 40fdc272ca27d0d812cb6d10e0809a391448407a Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Mon, 16 Feb 2026 03:58:16 -0600 Subject: [PATCH 6/8] Next step debugging --- STMMC/fcl/EventFilter.fcl | 42 +++++++ STMMC/fcl/FirePhotonToDetector.fcl | 60 +++++++++ .../fcl/HPGeWaveformGenerationAndAnalysis.fcl | 35 +++--- STMMC/fcl/ROOTAnalysisDump.fcl | 2 +- STMMC/fcl/prolog.fcl | 14 ++- .../HPGeWaveformsFromStepPointMCs_module.cc | 78 +++++++----- .../STMMovingWindowDeconvolution_module.cc | 119 ++++++++++++------ 7 files changed, 263 insertions(+), 87 deletions(-) create mode 100644 STMMC/fcl/EventFilter.fcl create mode 100644 STMMC/fcl/FirePhotonToDetector.fcl diff --git a/STMMC/fcl/EventFilter.fcl b/STMMC/fcl/EventFilter.fcl new file mode 100644 index 0000000000..e17ab65e2b --- /dev/null +++ b/STMMC/fcl/EventFilter.fcl @@ -0,0 +1,42 @@ +# Script to generate a separate file associated with a single event, for testing purposes. This is not intended to be used in production, but it can be useful for debugging and development. +# The event ID is specified in the 'idsToMatch' parameter of the EventIDFilter module. In this example, it is set to "1:0:5", which corresponds to run 1, subrun 0, event 5. You can change this value to select a different event. +# Mainly for use with the output of FirePhotonToDetector.fcl, as input to HPGeWaveformGenerationAndAnalysis.fcl, but it can be used with any input file that contains the specified event. +# Original author: Pawel Plesniak + +#include "Offline/fcl/minimalMessageService.fcl" + +process_name: STMCat + +source: { + module_type: RootInput +} + +physics: { + filters: { + isolateEvent: { + module_type: EventIDFilter + # The framework specifically asked for 'idsToMatch' + idsToMatch: [ "1:0:5" ] + } + } + + # This is the path name that SelectEvents must reference + STMCompressedPath: [ isolateEvent ] + trigger_paths: [ STMCompressedPath ] + o1 : [ out ] + end_paths: [ o1 ] +} + +services: { + message: @local::default_message +} + +outputs: { + out: { + module_type: RootOutput + fileName: "Filtered.art" + SelectEvents: ["STMCompressedPath"] + } +} + +services.scheduler.wantSummary: true diff --git a/STMMC/fcl/FirePhotonToDetector.fcl b/STMMC/fcl/FirePhotonToDetector.fcl new file mode 100644 index 0000000000..19e968ca10 --- /dev/null +++ b/STMMC/fcl/FirePhotonToDetector.fcl @@ -0,0 +1,60 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/STMMC/fcl/prolog.fcl" + +# This module simulates N photon(s) being fired upstream of the HPGe detector to validate the behaviour of the digitizer and analysis tools +# Original author : Pawel Plesniak + +process_name: FirePhotonToDetector +source : { + module_type : EmptyEvent + maxEvents : 1 +} + +services : { + # @local::Services.Sim + @table::Services.Core + @table::Services.SimOnly + message: @local::default_message + GlobalConstantsService : { inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" } +} + +physics: { + producers : { + generate : { + module_type : PhotonGun + x : @local::DigitizationTester.HPGe.x + y : @local::DigitizationTester.HPGe.y + z : @local::DigitizationTester.HPGe.z + E : 1.809 #@nil + verbose : true + } + g4run : @local::g4run + } + + p1 : [ generate, g4run ] + trigger_paths : [ p1 ] + o1 : [ out ] + end_paths: [ o1 ] +} + +outputs : { + out : { + module_type : RootOutput + fileName : "FirePhotonToDetector.art" + outputCommands: [ + "drop *_*_*_*", + "keep mu2e::StepPointMCs_g4run_STMDet_FirePhotonToDetector", + "keep mu2e::SimParticlemv_g4run__FirePhotonToDetector" + ] + } +} + +# physics.producers.g4run.physics.physicsListName: "QGSP_BERT_EMZ" +physics.producers.g4run.SDConfig.enableSD: [STMDet] +services.SeedService.baseSeed : 18 +services.SeedService.maxUniqueEngines : 20 +# services.TFileService.fileName : @local::Efficiency.OutputFilename +physics.producers.g4run.debug.trackingVerbosityLevel : 1 +physics.producers.g4run.debug.steppingVerbosityLevel : 1 diff --git a/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl index 3c0058ef24..fa10611b00 100644 --- a/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl +++ b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl @@ -31,16 +31,17 @@ physics : { StepPointMCsTagEle : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagEle StepPointMCsTagMu : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagMu StepPointMCsTag1809 : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTag1809 + # StepPointMCsTagEle: "g4run:STMDet" # DEBUGGING TODO - remove me! fADC : @local::STMDAQParameters.samplingFrequencies.HPGe #TODO - make this module work with multiple frequencies EnergyPerADCBin : @local::HPGeDigitization.EnergyPerADCBin - NoiseSD : @local::HPGeDigitization.PreamplifierNoiseSD + NoiseSD : 0.0 # @local::HPGeDigitization.PreamplifierNoiseSD risingEdgeDecayConstant : @local::HPGeDigitization.DecayConstant microspillBufferLengthCount : @local::HPGeDigitization.microspillBufferLengthCount makeTTree : false - makeADCPlot : false + makeADCPlot : false # Note - this does not work, I don't have the time to fix this. timeOffset : 0 resetEventNumber : @local::HPGeDigitization.resetEventNumber - verbosityLevel : 0 + verbosityLevel : 1 } concatenateWaveformsHPGe : { # Concatenate the microspill waveforms into macrospill waveforms @@ -58,7 +59,7 @@ physics : { threshold : @local::STMMCAnalysis.ZS.HPGe.threshold window : @local::STMMCAnalysis.ZS.HPGe.window naverage : @local::STMMCAnalysis.ZS.HPGe.naverage - verbosityLevel : 0 + verbosityLevel : 1 makeTTreeGradients: false makeTTreeWaveforms: false } @@ -66,6 +67,7 @@ physics : { # Perform moving window deconvolution on the concatenated (or zero-suppressed) waveforms to determine the pulse energies module_type : STMMovingWindowDeconvolution stmWaveformDigisTag : @local::STMMCAnalysis.MWD.HPGe.stmWaveformDigisTag.concatenated + # stmWaveformDigisTag: "DigiHPGe:" # DEBUGGING TODO - delete me! tau : @local::STMMCAnalysis.MWD.HPGe.tau M : @local::STMMCAnalysis.MWD.HPGe.M L : @local::STMMCAnalysis.MWD.HPGe.L @@ -74,9 +76,9 @@ physics : { defaultBaselineMean : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineMean.suppressed defaultBaselineSD : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineSD.suppressed makeTTreeMWD: false - makeTTreeEnergies: true + makeTTreeEnergies: false TTreeEnergyCalib : @local::HPGeDigitization.EnergyPerADCBin - verbosityLevel : 0 + verbosityLevel : 1 xAxis : "" } } @@ -86,25 +88,26 @@ physics : { STMWaveformDigisTag : @local::HPGeDigitization.concatenation.filterTag } } + test: [DigiHPGe, MWDHPGe] # DEBUGGING TODO - remove me! digitization_path_no_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, MWDHPGe] digitization_path_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] - trigger_paths : ["digitization_path_ZS"] # Populate me! TODO - make me @nil before git push - # output_path : [compressedOutput] - # end_paths : [output_path] + trigger_paths : ["digitization_path_no_ZS"] # Populate me! TODO - make me @nil before git push + output_path : [compressedOutput] + end_paths : [output_path] } outputs : { compressedOutput : { module_type : RootOutput - fileName : "dts.owner.HPGeReco.version.sequencer.art" - SelectEvents: [@nil] # Populate me! - # outputCommands: [ - # "drop *_*_*_*", - # "keep mu2e::STMMWDDigis_MWDHPGe_*_HPGeReco" - # ] + fileName : "dts.owner.HPGeReco.version.sequencer.art" # DEBUGGING: "test.art" + SelectEvents: ["digitization_path_no_ZS"] # TODO - delete this! Populate me! Comment below lines! + outputCommands: [ + "drop *_*_*_*", + "keep mu2e::STMMWDDigis_MWDHPGe_*_HPGeReco" + ] } } services.scheduler.wantSummary: true services.SeedService.baseSeed : 1 -services.TFileService.fileName : "nts.owner.HPGeRecoTree.version.sequencer.root" +# services.TFileService.fileName : "nts.owner.HPGeRecoTree.version.sequencer.root" diff --git a/STMMC/fcl/ROOTAnalysisDump.fcl b/STMMC/fcl/ROOTAnalysisDump.fcl index 43e7e88c9b..f42073a298 100644 --- a/STMMC/fcl/ROOTAnalysisDump.fcl +++ b/STMMC/fcl/ROOTAnalysisDump.fcl @@ -100,7 +100,7 @@ physics: { } } - o1 : ["HPGeEDeps"] # Populate me! + o1 : ["HPGeMWDspectra"] # Populate me! end_paths: [o1] } diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 3fe4cefd30..3033337323 100644 --- a/STMMC/fcl/prolog.fcl +++ b/STMMC/fcl/prolog.fcl @@ -259,10 +259,10 @@ STMDAQParameters : { # HPGe waveform generation parameters HPGeDigitization : { EnergyPerADCBin : 1 # keV/bin - PreamplifierNoiseSD : 0.0 # mV + PreamplifierNoiseSD : 0.32 # mV DecayConstant : 50 # us resetEventNumber : @local::DataProducts.nMicrospillsPerSpill - microspillBufferLengthCount : 1000 + microspillBufferLengthCount : 2 concatenation : { STMWaveformDigisTag : @local::DataProducts.Waveforms.HPGe.uSpill nMerge : @local::DataProducts.nMicrospillsPerSpill @@ -319,6 +319,16 @@ Efficiency : { } } } +DigitizationTester: { # For firing photons into the detector to validate the behaviour of the relevant algorithms + HPGe: { + x : -3958.55 # mm + y : 0.0 # mm + z : 40584.95 # mm + deltax : -1.0 # MeV/c + deltay : 0.0 # MeV/c + deltaz : 1.0 # MeV/c + } +} ######################################### End of edit lines ######################################### END_PROLOG diff --git a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc index 413151ad00..7eae499927 100644 --- a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc +++ b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc @@ -100,6 +100,7 @@ namespace mu2e { const double feedbackCapacitance = 1e-12; // [Farads] const double epsilonGe = 2.96; // Energy required to generate an eh pair in Ge at 77K [eV] const double micropulseTime = 1695.0; // [ns] + const double preamplifier_rise_time = 80.0; // [ns] // Define physics constants const double _e = 1.602176634e-19; // Electric charge constant [Coulombs] @@ -108,7 +109,8 @@ namespace mu2e { // ADC variables double chargeToADC = 0; // Conversion factor from charge built in capacitor to ADC determined voltage, multiply by this value to get from charge built to ADC voltage. - uint nADCs = 0; // Number of ADC values in an event + uint nADCs = 0; // Number of ADC values in an event, can change to account for timing variations + uint nADCs_init = 0; // Number of ADC values in an event, is fixed const int16_t ADCMax = static_cast((-1 * std::pow(2, 15)) + 1); // Maximum ADC value, power is 15 not 16 as using int16_t not uint16_t double ADC = 0; // iterator variable uint32_t eventTimeBuffer = 0; // Multiple of event ids to store @@ -146,7 +148,8 @@ namespace mu2e { double electronTravelTime = 0, holeTravelTime = 0; // Drift times [ns] uint32_t electronTravelTimeSteps = 0, holeTravelTimeSteps = 0; // Drift times [steps] double decayExp = 0; // Amount of decay with each tADC - double lastEventEndDecayedCharge = 0; // Carry over for starting new microspill waveforms [charge carrier pairs] + double lastEventEndDecayedCharge = 0; // Carry over the charge stored in the crystal from the previous event + double lastEventEndIntegratedCharge = 0; // Carry over for integrated charge stored in the preamplifier int microspillBufferLengthCount = 0; // Buffer to store the charge deposits that are allocated to this event but happen after the microspill ends e.g. 844keV const int defaultMicrospillBufferLengthCount = 2; // Default value for the microspill buffer length @@ -185,19 +188,19 @@ namespace mu2e { noiseSD(conf().noiseSD()), risingEdgeDecayConstant(conf().risingEdgeDecayConstant()) { produces(); - if (defaultMicrospillBufferLengthCount < 2) - throw cet::exception("RANGE", "defaultMicrospillBufferLengthCount has to be more than 1\n"); - crystalCentrePosition.set(crystalCentreX, crystalCentreY, crystalCentreZ); tADC = 1e3/fADC; // 1e3 converts [us] to [ns] as fADC is in [MHz] // Assign optional variables microspillBufferLengthCount = conf().microspillBufferLengthCount() ? *(conf().microspillBufferLengthCount()) : defaultMicrospillBufferLengthCount; + if (defaultMicrospillBufferLengthCount < 2) + throw cet::exception("RANGE", "defaultMicrospillBufferLengthCount has to be at least 1!\n"); verbosityLevel = conf().verbosityLevel() ? *(conf().verbosityLevel()) : 0; // Determine the number of ADC values in each STMWaveformDigi. Increase the number by one due to truncation. At 320MHz, this will be 543 ADC values per microbunch double _nADCs = (micropulseTime/tADC) + 1; - nADCs = (int) _nADCs; + nADCs = (uint) _nADCs; + nADCs_init = nADCs; _charge.insert(_charge.begin(), nADCs * microspillBufferLengthCount, 0.); _chargeCollected.insert(_chargeCollected.begin(), nADCs * microspillBufferLengthCount, 0.); _chargeCarryOver.insert(_chargeCarryOver.begin(), nADCs * (microspillBufferLengthCount - 1), 0.); @@ -243,7 +246,8 @@ namespace mu2e { void HPGeWaveformsFromStepPointMCs::beginJob() { if (verbosityLevel) { - std::cout << "STM HPGe digitization parameters" << std::endl; + std::cout << std::endl; + std::cout << "=======================================STM HPGe digitization parameters=======================================" << std::endl; std::cout << "\tInput parameters" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "fAD [MHz]" << fADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "EnergyPerADCBin [keV/bin]" << ADCToEnergy << std::endl; @@ -258,6 +262,7 @@ namespace mu2e { std::cout << std::left << "\t\t" << std::setw(60) << "tADC [ns]" << tADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "nADCs" << nADCs << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "NoiseSD [charge carriers]" << noiseSD << std::endl; + std::cout << std::left << "\t\t" << std::setw(60) << "NoiseSD [ADC bins]" << noiseSD * chargeToADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "chargeToADC [bin/charge carrier]" << chargeToADC << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range [V]" << "[+1, -1]" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range used [V]" << "[0, -1]" << std::endl; @@ -265,6 +270,7 @@ namespace mu2e { std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range used [charge carriers]" << "[0, " << ADCMax/chargeToADC << "]" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Voltage range used [C]" << "[0, " << ADCMax * _e/chargeToADC << "]" << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "Energy range [keV]" << "[0, " << -1*ADCMax * ADCToEnergy << "]" << std::endl; + std::cout << "==============================================================================================================" << std::endl; std::cout << std::endl; // buffer line }; }; @@ -276,21 +282,19 @@ namespace mu2e { // TODO - this only keeps accurate time if the sampling is 320MHz. Needs to be rewritten to work for other times if (std::abs(fADC - 320) > 1e-12) throw cet::exception("RANGE") << "Currently only fADC of 320 MHz is supported (got " << fADC << " MHz)\n"; - nADCs = 542; + nADCs = nADCs_init; // Assign the variables for time and number of ADCs eventTimeBuffer = eventId % 5; if (!(eventTimeBuffer == 0 || eventTimeBuffer == 3)) nADCs++; - eventTime += nADCs; // Update the last event decayed charge before the noise so the noise isn't added twice if (resetEventNumber != 0 && eventId == resetEventNumber) { - lastEventEndDecayedCharge = 0; + lastEventEndDecayedCharge = 0.0; + lastEventEndIntegratedCharge = 0.0; std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); } - else - lastEventEndDecayedCharge = _chargeDecayed.back(); // Update the parameters to carry over to the next event _chargeCollected.clear(); @@ -333,14 +337,6 @@ namespace mu2e { waveform_time = step.time(); }; }; - // Process the collected charge for each step - for (const StepPointMC& step : StepsEle) { - if (step.ionizingEdep() != 0) { - depositCharge(step); - if (step.time() < waveform_time) - waveform_time = step.time(); - } - } for(const StepPointMC& step : StepsMu){ if (step.ionizingEdep() != 0) depositCharge(step); @@ -369,14 +365,13 @@ namespace mu2e { throw cet::exception("LogicError", "ADC values too high!"); }; - // STMWaveformDigi _waveformDigi(eventTime, _adcs); - STMWaveformDigi _waveformDigi(waveform_time, _adcs); + STMWaveformDigi _waveformDigi(eventTime, _adcs); std::unique_ptr outputDigis(new STMWaveformDigiCollection); outputDigis->emplace_back(_waveformDigi); // Make the ttree if appropriate if (makeTTree) { - time = _waveformDigi.trigTimeOffset() - 1; + time = eventTime; for (uint i = 0; i < nADCs; i++) { chargeCollected = _chargeCollected[i]; chargeDecayed = _chargeDecayed[i]; @@ -404,6 +399,9 @@ namespace mu2e { // Add the STMWaveformDigi to the event event.put(std::move(outputDigis)); + + // Update the event time for the next waveform + eventTime += nADCs; return; }; @@ -482,7 +480,7 @@ namespace mu2e { N_ehPairs = -1.0 * step.ionizingEdep() * 1e6 / epsilonGe; // 1e6 converts MeV to eV. -1.0 as this is a decreasing peak // Define parameters required for charge deposition. Constants A and B are defined here for code brevity - uint tIndex = (step.time() + timeOffset) / tADC, tIndexStart = tIndex; + uint tIndex = (step.time() + timeOffset) / tADC + 1000, tIndexStart = tIndex; const double A = N_ehPairs / log(R2/R1); const double Be = electronDriftVelocity / R0; const double Bh = holeDriftVelocity / R0; @@ -513,7 +511,8 @@ namespace mu2e { }; // Allocate the rest of the charge for one more entry for the continuity - _charge[tIndex] = N_ehPairs; + for (uint i = tIndex; i < _charge.size(); i++) + _charge[i] = _charge[tIndex - 1]; tIndex++; // Update _chargeCollected. First case is treated separately as there is no charge deposited in the previous step @@ -530,10 +529,28 @@ namespace mu2e { }; void HPGeWaveformsFromStepPointMCs::decayCharge() { - _chargeDecayed[0] = lastEventEndDecayedCharge * decayExp + _chargeCollected[0]; - for (uint t = 1; t < nADCs; t++) - _chargeDecayed[t] = _chargeDecayed[t-1] * decayExp + _chargeCollected[t]; - return; + // Convert time constants to samples based on 320MHz clock, 80ns is the approximate ORTEC GMX preamplifier rise time + const double tau_rise_samples = preamplifier_rise_time / tADC; // Convert this to ADC steps + const double riseExp = 1.0 / (1.0 + tau_rise_samples); + + // Carry over the current in the last time step + double currentDecayed = lastEventEndDecayedCharge; + + // Holiding variable + double currentIntegrated = lastEventEndIntegratedCharge; + + for (uint t = 0; t < nADCs; t++) { + // Discharge the capacitor, add the next contribution of the step + currentDecayed = currentDecayed * decayExp + _chargeCollected[t]; + + // Smooth out the waveform accounting for the finite bandwidth of the preamplifier + currentIntegrated = riseExp * currentDecayed + (1.0 - riseExp) * currentIntegrated; + + // Store the results + _chargeDecayed[t] = currentIntegrated; + } + lastEventEndDecayedCharge = currentDecayed; + lastEventEndIntegratedCharge = currentIntegrated; }; void HPGeWaveformsFromStepPointMCs::addNoise() { @@ -553,13 +570,14 @@ namespace mu2e { // Convert the charge deposition to ADC voltage output. for (uint i = 0; i < nADCs; i++) { ADC = _chargeDecayed[i] * chargeToADC; - _adcs[i] = ADC > ADCMax ? static_cast(std::round(ADC)) : ADCMax; + _adcs[i] = (ADC > ADCMax) ? static_cast(std::round(ADC)) : ADCMax; }; return; }; void HPGeWaveformsFromStepPointMCs::endJob() { mf::LogInfo log("STMResamplingFilter summary"); + log << "\n"; log << "=====HPGeWaveformsFromStepPointMCs summary=====\n"; log << "No. kept events: " << kept_events << "\n"; log << "No. discarded events: " << dropped_events << "\n"; diff --git a/STMReco/src/STMMovingWindowDeconvolution_module.cc b/STMReco/src/STMMovingWindowDeconvolution_module.cc index 6da8c05fc4..8efce6494a 100644 --- a/STMReco/src/STMMovingWindowDeconvolution_module.cc +++ b/STMReco/src/STMMovingWindowDeconvolution_module.cc @@ -140,6 +140,9 @@ namespace mu2e { // Summary data uint processedEvents = 0, processedWaveforms = 0, foundPeaks = 0; + + // Debugging variables + double waveform_time_microspill_frame = 0.0; }; STMMovingWindowDeconvolution::STMMovingWindowDeconvolution(const Parameters& conf) : @@ -191,7 +194,8 @@ namespace mu2e { void STMMovingWindowDeconvolution::beginJob() { if (verbosityLevel) { - std::cout << "STM Moving Window Deconvolution" << std::endl; + std::cout << std::endl; + std::cout << "=======================================STM HPGe digitization parameters=======================================" << std::endl; std::cout << "\tAlgorithm parameters" << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "tau" << tau << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "M" << M << std::endl; @@ -201,6 +205,7 @@ namespace mu2e { std::cout << "\tChannel: " << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "Name" << channel.name() << std::endl; std::cout << std::left << "\t\t" << std::setw(15) << "ID" << static_cast(channel.id()) << std::endl; + std::cout << "==============================================================================================================" << std::endl; std::cout << std::endl; // buffer line }; }; @@ -240,11 +245,13 @@ namespace mu2e { differentiate(); average(); calculate_baseline(); + baseline_mean = 0.0; + baseline_stddev = 2.0; find_peaks(); nPeaks = peak_heights.size(); foundPeaks += nPeaks; - if (nPeaks && verbosityLevel) // TODO - change to verbosityLevel + if (verbosityLevel > 0) std::cout << "MWD: found " << nPeaks << " peaks in event " << event.id() << std::endl; for (i = 0; i < nPeaks; ++i) { mwd_energy = (peak_heights[i] < ADCMax) ? ADCMax : static_cast(peak_heights[i]); // When saturating the int16_t limit, deconvolution goes below the int16_t limit so the energy turns negative. This clips the energy and the limit @@ -339,59 +346,95 @@ namespace mu2e { }; }; - void STMMovingWindowDeconvolution::calculate_baseline() { +void STMMovingWindowDeconvolution::calculate_baseline() { i = M; foundBaselineData = false; using namespace boost::accumulators; accumulator_set > acc_data_without_peaks; - // Remove peaks so that we can calculate the baseline of the averaged data - while (i < nADCs){ - gradient = averaged_data[i + 1] - averaged_data[i]; - if(gradient < thresholdgrad) // if the gradient is too sharp (i.e. we have hit a peak) - i += (M + 2 * L); // jump ahead a little bit - else { - acc_data_without_peaks(averaged_data[i]); - i++; - foundBaselineData = true; - }; - }; - baseline_mean = foundBaselineData ? extract_result(acc_data_without_peaks) : defaultBaselineMean; - baseline_stddev = foundBaselineData ? std::sqrt(extract_result(acc_data_without_peaks)) : defaultBaselineSD; - }; - void STMMovingWindowDeconvolution::find_peaks() { - threshold_cut = baseline_mean - nsigma_cut * baseline_stddev; + // Safety: stop one sample early to allow for the i+1 gradient check + while (i < nADCs - 1) { + gradient = averaged_data[i + 1] - averaged_data[i]; + + // Use -1.0 to skip pulses while keeping the baseline calculation stable + if (gradient < thresholdgrad) { + i += (static_cast(M) + 2 * static_cast(L)); + if (i >= nADCs) break; + } + else { + acc_data_without_peaks(averaged_data[i]); + i++; + foundBaselineData = true; + } + } + + if (foundBaselineData) { + baseline_mean = extract_result(acc_data_without_peaks); + double calculated_sd = std::sqrt(extract_result(acc_data_without_peaks)); + + // DEBUG: Force a tiny floor since NoiseSD is 0.0 + baseline_stddev = std::max(calculated_sd, 2.0); + } else { + baseline_mean = defaultBaselineMean; + baseline_stddev = defaultBaselineSD; + } +} + +void STMMovingWindowDeconvolution::find_peaks() { + // 1. Use a more conservative trigger for the recovery phase. + double trigger_threshold = baseline_mean - 100.0; + lowest_height = 0; - lowest_height_time = -1; // in clock ticks + lowest_height_time = -1; + bool is_armed = true; // NEW: Ensure we cross ABOVE threshold before re-arming for(i = M; i < nADCs; i++){ - if (averaged_data[i] < threshold_cut) { // the waveforms are negative so if we go below this threshold we have seen a peak - if (averaged_data[i] < averaged_data[i - 1] && averaged_data[i] < lowest_height){ // if the current value is lower than the previous value and lower than the lowest value we've seen so far - lowest_height = averaged_data[i]; // record the lowest height - if (lowest_height_time == -1) - lowest_height_time = i; // record the time we cross the threshold + // Only allow a new peak if we are 'armed' (meaning we've been above threshold) + if (is_armed && averaged_data[i] < trigger_threshold) { + if (averaged_data[i] < lowest_height) { + lowest_height = averaged_data[i]; + lowest_height_time = i; + } } - else - continue; - }; - if (lowest_height_time == -1) // this will be true if we haven't seen a peak yet - continue; - else if (averaged_data[i] > threshold_cut) { // if we have seen a peak and go above the cut - peak_heights.push_back(lowest_height - baseline_mean); - peak_times.push_back(lowest_height_time); // ct - lowest_height_time = -1; // reset to 0 so we can find a new peak - lowest_height = 0; - }; - }; - }; + else if (averaged_data[i] >= trigger_threshold) { + // If we are above threshold, we are officially 'armed' and ready for a new pulse + is_armed = true; + + // If we just finished a pulse, record it + if (lowest_height_time != -1) + double amplitude = std::abs(lowest_height - baseline_mean); + if (amplitude >= 200.0) { + peak_heights.push_back(lowest_height - baseline_mean); + peak_times.push_back(lowest_height_time); + // Move forward + i += (M + 2 * L); + is_armed = false; // Force the algorithm to wait for baseline crossing + if (i >= nADCs) break; + } + lowest_height = 0; + lowest_height_time = -1; + } + } + } +} void STMMovingWindowDeconvolution::endJob() { mf::LogInfo log("MWD summary"); + log << "\n"; log << "==================MWD summary==================\n"; log << std::left << std::setw(25) << "\tProcessed events: " << processedEvents << "\n"; log << std::left << std::setw(25) << "\tProcessed waveforms: " << processedWaveforms << "\n"; log << std::left << std::setw(25) << "\tFound peaks: " << foundPeaks << "\n"; log << "===============================================\n"; + if (verbosityLevel) { + log << "\n"; + log << "==================Parameters==================\n"; + log << std::left << std::setw(25) << "\tnsPerCt: " << nsPerCt << "\n"; + log << std::left << std::setw(25) << "\tPedestal: " << pedestal << "\n"; + log << std::left << std::setw(25) << "\tBaseline mean: " << baseline_mean << "\n"; + log << std::left << std::setw(25) << "\tBaseline std dev: " << baseline_stddev << "\n"; + log << "===============================================\n"; + } return; }; From 1d7da45015ac09997a9f73d3d848ae52c8b1bac2 Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Mon, 16 Feb 2026 05:15:32 -0600 Subject: [PATCH 7/8] Removing remaining debug script --- .../fcl/HPGeWaveformGenerationAndAnalysis.fcl | 4 +- .../HPGeWaveformsFromStepPointMCs_module.cc | 2 +- .../STMMovingWindowDeconvolution_module.cc | 58 +++++++++---------- 3 files changed, 32 insertions(+), 32 deletions(-) diff --git a/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl index fa10611b00..d932c3e31b 100644 --- a/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl +++ b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl @@ -75,7 +75,7 @@ physics : { thresholdgrad : @local::STMMCAnalysis.MWD.HPGe.thresholdgrad defaultBaselineMean : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineMean.suppressed defaultBaselineSD : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineSD.suppressed - makeTTreeMWD: false + makeTTreeMWD: true makeTTreeEnergies: false TTreeEnergyCalib : @local::HPGeDigitization.EnergyPerADCBin verbosityLevel : 1 @@ -110,4 +110,4 @@ outputs : { services.scheduler.wantSummary: true services.SeedService.baseSeed : 1 -# services.TFileService.fileName : "nts.owner.HPGeRecoTree.version.sequencer.root" +services.TFileService.fileName : "nts.owner.HPGeRecoTree.version.sequencer.root" diff --git a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc index 7eae499927..30f9d815f5 100644 --- a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc +++ b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc @@ -480,7 +480,7 @@ namespace mu2e { N_ehPairs = -1.0 * step.ionizingEdep() * 1e6 / epsilonGe; // 1e6 converts MeV to eV. -1.0 as this is a decreasing peak // Define parameters required for charge deposition. Constants A and B are defined here for code brevity - uint tIndex = (step.time() + timeOffset) / tADC + 1000, tIndexStart = tIndex; + uint tIndex = (step.time() + timeOffset) / tADC, tIndexStart = tIndex; const double A = N_ehPairs / log(R2/R1); const double Be = electronDriftVelocity / R0; const double Bh = holeDriftVelocity / R0; diff --git a/STMReco/src/STMMovingWindowDeconvolution_module.cc b/STMReco/src/STMMovingWindowDeconvolution_module.cc index 8efce6494a..550a349955 100644 --- a/STMReco/src/STMMovingWindowDeconvolution_module.cc +++ b/STMReco/src/STMMovingWindowDeconvolution_module.cc @@ -388,35 +388,35 @@ void STMMovingWindowDeconvolution::find_peaks() { lowest_height_time = -1; bool is_armed = true; // NEW: Ensure we cross ABOVE threshold before re-arming - for(i = M; i < nADCs; i++){ - // Only allow a new peak if we are 'armed' (meaning we've been above threshold) - if (is_armed && averaged_data[i] < trigger_threshold) { - if (averaged_data[i] < lowest_height) { - lowest_height = averaged_data[i]; - lowest_height_time = i; - } - } - else if (averaged_data[i] >= trigger_threshold) { - // If we are above threshold, we are officially 'armed' and ready for a new pulse - is_armed = true; - - // If we just finished a pulse, record it - if (lowest_height_time != -1) - double amplitude = std::abs(lowest_height - baseline_mean); - if (amplitude >= 200.0) { - peak_heights.push_back(lowest_height - baseline_mean); - peak_times.push_back(lowest_height_time); - // Move forward - i += (M + 2 * L); - is_armed = false; // Force the algorithm to wait for baseline crossing - if (i >= nADCs) break; - } - lowest_height = 0; - lowest_height_time = -1; - } - } - } -} + for(i = M; i < nADCs; i++) { + // Only allow a new peak if we are 'armed' (meaning we've been above threshold) + if (is_armed && averaged_data[i] < trigger_threshold) { + if (averaged_data[i] < lowest_height) { + lowest_height = averaged_data[i]; + lowest_height_time = i; + } + } + else if (averaged_data[i] >= trigger_threshold) { + // If we are above threshold, we are officially 'armed' and ready for a new pulse + is_armed = true; + + // If we just finished a pulse, record it + if (lowest_height_time != -1) { + double amplitude = std::abs(lowest_height - baseline_mean); + if (amplitude >= 200.0) { + peak_heights.push_back(lowest_height - baseline_mean); + peak_times.push_back(lowest_height_time); + // Move forward + i += (M + 2 * L); + is_armed = false; // Force the algorithm to wait for baseline crossing + if (i >= nADCs) break; + }; + lowest_height = 0; + lowest_height_time = -1; + }; + }; + }; + }; void STMMovingWindowDeconvolution::endJob() { mf::LogInfo log("MWD summary"); From 1265f1427f28559a3d69ec55381721b97b686143 Mon Sep 17 00:00:00 2001 From: PawelPlesniak Date: Wed, 18 Feb 2026 03:41:42 -0600 Subject: [PATCH 8/8] MWD now no longer complains about things --- .../fcl/HPGeWaveformGenerationAndAnalysis.fcl | 2 +- .../HPGeWaveformsFromStepPointMCs_module.cc | 57 +++---- .../STMMovingWindowDeconvolution_module.cc | 158 ++++++++++++------ 3 files changed, 130 insertions(+), 87 deletions(-) diff --git a/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl index d932c3e31b..e574407859 100644 --- a/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl +++ b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl @@ -33,7 +33,7 @@ physics : { StepPointMCsTag1809 : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTag1809 # StepPointMCsTagEle: "g4run:STMDet" # DEBUGGING TODO - remove me! fADC : @local::STMDAQParameters.samplingFrequencies.HPGe #TODO - make this module work with multiple frequencies - EnergyPerADCBin : @local::HPGeDigitization.EnergyPerADCBin + EnergyPerADCBin : 10 # @local::HPGeDigitization.EnergyPerADCBin NoiseSD : 0.0 # @local::HPGeDigitization.PreamplifierNoiseSD risingEdgeDecayConstant : @local::HPGeDigitization.DecayConstant microspillBufferLengthCount : @local::HPGeDigitization.microspillBufferLengthCount diff --git a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc index 30f9d815f5..770f0124d8 100644 --- a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc +++ b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc @@ -318,37 +318,28 @@ namespace mu2e { std::vector Steps1809 = event.getProduct(StepPointMCsToken1809); // If there are no hits, produce an empty STMWaveformDigiCollection and return - if (StepsEle.size() == 0 && StepsMu.size() == 0 && Steps1809.size() == 0) { - std::vector emptyAdcs(nADCs, static_cast(0)); - STMWaveformDigi _waveformDigi(0.0, emptyAdcs); - std::unique_ptr outputDigis(new STMWaveformDigiCollection); - outputDigis->emplace_back(_waveformDigi); - event.put(std::move(outputDigis)); - dropped_events++; - return; + bool hasHits = (!StepsEle.empty() || !StepsMu.empty() || !Steps1809.empty()); + + if (hasHits) { + kept_events++; + // Add a collection of charge depositions to _charge + for(const StepPointMC& step : StepsEle) { + if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) waveform_time = step.time(); + } + for(const StepPointMC& step : StepsMu) { + if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) waveform_time = step.time(); + } + for(const StepPointMC& step : Steps1809) { + if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) waveform_time = step.time(); + } + } else { + dropped_events++; + // We do nothing here. _chargeCollected is already 0s from initialization. + // We just proceed to decayCharge(). } - kept_events++; - - // Add a collection of charge depositions to _charge - for(const StepPointMC& step : StepsEle){ - if (step.ionizingEdep() != 0) { - depositCharge(step); - if (step.time() < waveform_time) - waveform_time = step.time(); - }; - }; - for(const StepPointMC& step : StepsMu){ - if (step.ionizingEdep() != 0) - depositCharge(step); - if (step.time() < waveform_time) - waveform_time = step.time(); - }; - for(const StepPointMC& step : Steps1809){ - if (step.ionizingEdep() != 0) - depositCharge(step); - if (step.time() < waveform_time) - waveform_time = step.time(); - }; // Decay all of the collected charges decayCharge(); @@ -576,11 +567,11 @@ namespace mu2e { }; void HPGeWaveformsFromStepPointMCs::endJob() { - mf::LogInfo log("STMResamplingFilter summary"); + mf::LogInfo log("HPGeWaevfrmsFromStepPointMCs"); log << "\n"; log << "=====HPGeWaveformsFromStepPointMCs summary=====\n"; - log << "No. kept events: " << kept_events << "\n"; - log << "No. discarded events: " << dropped_events << "\n"; + log << std::left << std::setw(25) << "\tNo. kept events: " << kept_events << "\n"; + log << std::left << std::setw(25) << "\tNo. discarded events: " << dropped_events << "\n"; log << "===============================================\n"; }; }; // namespace mu2e diff --git a/STMReco/src/STMMovingWindowDeconvolution_module.cc b/STMReco/src/STMMovingWindowDeconvolution_module.cc index 550a349955..4383ba51c1 100644 --- a/STMReco/src/STMMovingWindowDeconvolution_module.cc +++ b/STMReco/src/STMMovingWindowDeconvolution_module.cc @@ -115,6 +115,14 @@ namespace mu2e { double sum = 0.0; // variable part of avarege() int count = 0; // counter variable int16_t mwd_energy = 0; + bool _isFirstWaveform = true; + double _last_ADC = 0.0; + double _last_deconvolved = 0.0; + // Rolling history buffers to bridge the gap between chunks + std::deque _deconvolved_history; // Stores last M samples + std::deque _differentiated_history; // Stores last L samples + double _last_average = 0.0; // Stores the last average value from previous chunk + // Peak finding variables double baseline_mean = 0.0; @@ -139,7 +147,7 @@ namespace mu2e { uint eventId = 0, waveformID = 0; // Summary data - uint processedEvents = 0, processedWaveforms = 0, foundPeaks = 0; + uint processedEvents = 0, processedWaveforms = 0, foundPeaks = 0, processedADCs; // Debugging variables double waveform_time_microspill_frame = 0.0; @@ -235,6 +243,7 @@ namespace mu2e { nADCs = ADCs.size(); processedWaveforms++; waveformID++; + processedADCs+=nADCs; if (M > nADCs) M = nADCs; @@ -245,8 +254,9 @@ namespace mu2e { differentiate(); average(); calculate_baseline(); - baseline_mean = 0.0; - baseline_stddev = 2.0; + // CHANGE HERE + // baseline_mean = 0.0; + // baseline_stddev = 2.0; find_peaks(); nPeaks = peak_heights.size(); @@ -293,58 +303,96 @@ namespace mu2e { }; - void STMMovingWindowDeconvolution::deconvolve() { - if (verbosityLevel > 4) { - std::cout << "MWD: input ADCs (" << ADCs.size() << "): "; - for (int16_t data : ADCs) - std::cout << data << ", "; - std::cout << "\n" << std::endl; - }; +void STMMovingWindowDeconvolution::deconvolve() { + // Handle the boundary (Index 0) + if (_isFirstWaveform) { deconvolved_data.push_back(ADCs[0] - pedestal); - for(i = 1; i < nADCs; i++) - deconvolved_data.push_back((ADCs[i] - pedestal) - timeFactor * (ADCs[i - 1] - pedestal) + deconvolved_data[i - 1]); - if (verbosityLevel > 4) { - std::cout << "MWD: deconvoluted data (" << deconvolved_data.size() << "): "; - for (double data : deconvolved_data) - std::cout << data << ", "; - std::cout << "\n" << std::endl; - }; - }; + _isFirstWaveform = false; + } else { + // Stitch to the previous chunk's last known state + deconvolved_data.push_back((ADCs[0] - pedestal) - timeFactor * (_last_ADC - pedestal) + _last_deconvolved); + } + + // Process the rest of the vector + for(i = 1; i < nADCs; i++) + deconvolved_data.push_back((ADCs[i] - pedestal) - timeFactor * (ADCs[i - 1] - pedestal) + deconvolved_data[i - 1]); + + // Save state for the next chunk + _last_ADC = ADCs.back(); + _last_deconvolved = deconvolved_data.back(); +} - void STMMovingWindowDeconvolution::differentiate() { - for (i = 0; i < M; i++) - differentiated_data.push_back(deconvolved_data[i]); - for (i = M; i < nADCs; i++) - differentiated_data.push_back(deconvolved_data[i] - deconvolved_data[i - M]); - if (verbosityLevel > 4) { - std::cout << "MWD: differentiated data (" << differentiated_data.size() << "): "; - for (double data : differentiated_data) - std::cout << data << ", "; - std::cout << "\n" << std::endl; - }; - }; +void STMMovingWindowDeconvolution::differentiate() { + double subtrahend = 0.0; - void STMMovingWindowDeconvolution::average() { - // sum the first L-1 elements of differentiated data - // and set the first L-1 elements of averaged data - sum = 0.0; - for (i = 0; i < L - 1; i++) { - sum += differentiated_data[i]; - averaged_data.push_back(differentiated_data[i]); // TODO: should this be divided by sum/i? - }; - sum += differentiated_data[L - 1]; - averaged_data.push_back(sum/L); - for (i = L; i < nADCs; ++i) { - sum += differentiated_data[i] - differentiated_data[i - L]; // move the sum across one sample - averaged_data.push_back(sum/L); - }; - if (verbosityLevel > 4) { - std::cout << "MWD: averaged data (" << averaged_data.size() << "): "; - for (double data : averaged_data) - std::cout << data << ", "; - std::cout << "\n" << std::endl; - }; - }; + for (i = 0; i < nADCs; i++) { + if (i >= M) { + // Standard case: Look back inside current vector + subtrahend = deconvolved_data[i - M]; + } else { + // Boundary case: Look back into history buffer + // If history is empty (very first event), assume 0. Otherwise fetch from deque. + if (_deconvolved_history.size() >= M) { + // The deque stores the LAST M samples. + // If i=0, we need the sample from M steps ago (index 0 of deque) + subtrahend = _deconvolved_history[i]; + } else { + subtrahend = deconvolved_data[0]; // Fallback for very first start + } + } + differentiated_data.push_back(deconvolved_data[i] - subtrahend); + } + + // Update History Buffer for the NEXT chunk + // We need to keep the LAST M samples of the current deconvolved_data + for (double val : deconvolved_data) { + _deconvolved_history.push_back(val); + if (_deconvolved_history.size() > M) _deconvolved_history.pop_front(); + } +} + +void STMMovingWindowDeconvolution::average() { + // We use a recursive moving average: Avg[i] = Avg[i-1] + (New - Old)/L + + for(i = 0; i < nADCs; i++) { + double added_val = differentiated_data[i]; + double removed_val = 0.0; + + if (i >= L) { + // Standard case: remove the sample L steps back in the current vector + removed_val = differentiated_data[i - L]; + } else { + // Boundary case: Look into history buffer + if (_differentiated_history.size() >= L) { + // The deque stores the LAST L samples from the previous chunk. + // If i=0, we need the oldest sample in history (index 0). + // If i=1, we need index 1, etc. + removed_val = _differentiated_history[i]; + } else { + removed_val = 0.0; // Startup case (first event ever) + } + } + + // Get the previous average. + // If i=0, use the stored average from the end of the previous chunk. + double prev_avg = (i == 0) ? _last_average : averaged_data.back(); + + // Calculate new average + double new_avg = prev_avg + (added_val - removed_val) / L; + averaged_data.push_back(new_avg); + } + + // Update State for the NEXT chunk + if (!averaged_data.empty()) { + _last_average = averaged_data.back(); + } + + // Update History Buffer + for (double val : differentiated_data) { + _differentiated_history.push_back(val); + if (_differentiated_history.size() > L) _differentiated_history.pop_front(); + } +} void STMMovingWindowDeconvolution::calculate_baseline() { i = M; @@ -373,7 +421,9 @@ void STMMovingWindowDeconvolution::calculate_baseline() { double calculated_sd = std::sqrt(extract_result(acc_data_without_peaks)); // DEBUG: Force a tiny floor since NoiseSD is 0.0 - baseline_stddev = std::max(calculated_sd, 2.0); + // baseline_stddev = std::max(calculated_sd, 2.0); // NOTE - THIS WAS FOR DEBUGGING ONLY + + baseline_stddev = calculated_sd; } else { baseline_mean = defaultBaselineMean; baseline_stddev = defaultBaselineSD; @@ -425,6 +475,8 @@ void STMMovingWindowDeconvolution::find_peaks() { log << std::left << std::setw(25) << "\tProcessed events: " << processedEvents << "\n"; log << std::left << std::setw(25) << "\tProcessed waveforms: " << processedWaveforms << "\n"; log << std::left << std::setw(25) << "\tFound peaks: " << foundPeaks << "\n"; + log << std::left << std::setw(25) << "\tnADCs: " << nADCs << "\n"; + log << std::left << std::setw(25) << "\tprocessedADCs: " << processedADCs << "\n"; log << "===============================================\n"; if (verbosityLevel) { log << "\n";