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/EventGenerator/src/PhotonGun_module.cc b/EventGenerator/src/PhotonGun_module.cc index 6b9e4359e2..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" @@ -38,7 +39,11 @@ 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")}; + fhicl::OptionalAtom verbose{ Name("verbose"), Comment("Print verbose messages")}; }; using Parameters = art::EDProducer::Table; explicit PhotonGun(const Parameters& conf); @@ -46,7 +51,9 @@ 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; + bool verbose = false; }; PhotonGun::PhotonGun(const Parameters& conf): @@ -56,17 +63,33 @@ 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 ((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) { 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()) { + if (verbose) + std::cout << "PhotonGun: Using delta position 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/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/Absorber.fcl b/STMMC/fcl/Absorber.fcl index 14c1785f8d..bf07469418 100644 --- a/STMMC/fcl/Absorber.fcl +++ b/STMMC/fcl/Absorber.fcl @@ -25,18 +25,18 @@ physics: { producers : { generate : { module_type : PhotonGun - x : @local::Efficiency.Absorber.x - y : @local::Efficiency.Absorber.y - z : @local::Efficiency.Absorber.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 } analyzers : { EDep : { - module_type : MakeVirtualDetectorTree - VirtualDetectorId : 101 + module_type : VirtualDetectorTree StepPointMCsTag: "g4run:virtualdetector" SimParticlemvTag: "g4run:" } @@ -50,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/AbsorberFromSTLaBr.fcl b/STMMC/fcl/AbsorberFromSTLaBr.fcl new file mode 100644 index 0000000000..a128c021ce --- /dev/null +++ b/STMMC/fcl/AbsorberFromSTLaBr.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.LaBr.x + deltay : @local::Efficiency.FromSTToDet.LaBr.y + deltaz : @local::Efficiency.FromSTToDet.LaBr.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/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/HPGeReco.fcl b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl similarity index 58% rename from STMMC/fcl/HPGeReco.fcl rename to STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl index b6dc9d5362..e574407859 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,27 +26,32 @@ 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 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 + EnergyPerADCBin : 10 # @local::HPGeDigitization.EnergyPerADCBin + NoiseSD : 0.0 # @local::HPGeDigitization.PreamplifierNoiseSD risingEdgeDecayConstant : @local::HPGeDigitization.DecayConstant microspillBufferLengthCount : @local::HPGeDigitization.microspillBufferLengthCount makeTTree : 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 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 @@ -45,13 +59,15 @@ 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 } 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 + # stmWaveformDigisTag: "DigiHPGe:" # DEBUGGING TODO - delete me! tau : @local::STMMCAnalysis.MWD.HPGe.tau M : @local::STMMCAnalysis.MWD.HPGe.M L : @local::STMMCAnalysis.MWD.HPGe.L @@ -59,10 +75,10 @@ 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 : 0 + verbosityLevel : 1 xAxis : "" } } @@ -72,8 +88,10 @@ physics : { STMWaveformDigisTag : @local::HPGeDigitization.concatenation.filterTag } } - digitization_path : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] - trigger_paths : [digitization_path] + 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_no_ZS"] # Populate me! TODO - make me @nil before git push output_path : [compressedOutput] end_paths : [output_path] } @@ -81,8 +99,8 @@ physics : { outputs : { compressedOutput : { module_type : RootOutput - fileName : "dts.owner.HPGeReco.version.sequencer.art" - SelectEvents: ["digitization_path"] + 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" diff --git a/STMMC/fcl/MakeTree.fcl b/STMMC/fcl/MakeTree.fcl deleted file mode 100644 index 1bcf442025..0000000000 --- a/STMMC/fcl/MakeTree.fcl +++ /dev/null @@ -1,56 +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.StepPointMCsTag - 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 : @nil - end_paths: [o1] -} - -services.TFileService.fileName : @nil diff --git a/STMMC/fcl/ROOTAnalysisDump.fcl b/STMMC/fcl/ROOTAnalysisDump.fcl new file mode 100644 index 0000000000..f42073a298 --- /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 energy depositions in crystal - MC truth + HPGeEDeps : { + module_type : HPGeTree + 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 +# } + 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 : ["HPGeMWDspectra"] # Populate me! + end_paths: [o1] +} + +services.TFileService.fileName : "tree.root" # Populate me! diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 1b56b4f7d1..3033337323 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 : { @@ -13,7 +75,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 @@ -34,76 +96,39 @@ ComponentPositions : { # In beam order LaBr : { # aperture centre x : -3863.4 y : 0.0 + z : 40404 } -} - -# 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 + ST : { # ST centre + x : -3904.0 y : 0.0 - z : 40300.9 + z : 627.0 } - 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 @@ -114,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 @@ -130,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 : { + StepPointMCsTag : { + 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 } } @@ -174,25 +255,80 @@ STMDAQParameters : { LaBr : 370.370370370 # MSPS - copied from Offline/STMConditions/fcl/prolog.fcl } } + +# HPGe waveform generation parameters HPGeDigitization : { - EnergyPerADCBin : 0.5 # keV/bin + EnergyPerADCBin : 1 # keV/bin PreamplifierNoiseSD : 0.32 # mV DecayConstant : 50 # us - resetEventNumber : 31858 - microspillBufferLengthCount : 1000 + resetEventNumber : @local::DataProducts.nMicrospillsPerSpill + microspillBufferLengthCount : 2 concatenation : { - STMWaveformDigisTag : "DigiHPGe:" - nMerge : 31858 - filterTag : "concatenateWaveformsHPGe:" + STMWaveformDigisTag : @local::DataProducts.Waveforms.HPGe.uSpill + nMerge : @local::DataProducts.nMicrospillsPerSpill + filterTag : @local::DataProducts.Waveforms.HPGe.concatenated } } -# Normalization parameter determination studies -# Detector efficiency simulation +# 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.concatenated + } +} + +# Detector efficiency simulation parameters Efficiency : { NPhotons : 1e7 - PhotonEnergy : 1.809 # MeV - OutputFilename : "Efficiency.root" + PhotonEnergy : @nil # 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 + } + } +} +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/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"; diff --git a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc index f084862a4e..770f0124d8 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 @@ -91,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] @@ -99,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 @@ -137,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 @@ -155,6 +167,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); @@ -172,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.); @@ -205,7 +221,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,13 +232,22 @@ 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; }; 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; @@ -230,12 +255,14 @@ 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; 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; @@ -243,42 +270,80 @@ 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 }; }; 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 = nADCs_init; + + // Assign the variables for time and number of ADCs + eventTimeBuffer = eventId % 5; + if (!(eventTimeBuffer == 0 || eventTimeBuffer == 3)) + nADCs++; + + // Update the last event decayed charge before the noise so the noise isn't added twice + if (resetEventNumber != 0 && eventId == resetEventNumber) { + lastEventEndDecayedCharge = 0.0; + lastEventEndIntegratedCharge = 0.0; + std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); + } + + // 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); - // Add a collection of charge depositions to _charge - for(const StepPointMC& step : StepsEle){ - if (step.ionizingEdep() != 0) - depositCharge(step); - }; - for(const StepPointMC& step : StepsMu){ - if (step.ionizingEdep() != 0) - depositCharge(step); - }; - for(const StepPointMC& step : Steps1809){ - if (step.ionizingEdep() != 0) - depositCharge(step); - }; + // If there are no hits, produce an empty STMWaveformDigiCollection and 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(). + } // 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,21 +356,13 @@ 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); 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]; @@ -315,18 +372,27 @@ 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)); + + // Update the event time for the next waveform + eventTime += nADCs; return; }; @@ -436,7 +502,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 @@ -453,10 +520,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() { @@ -476,10 +561,19 @@ 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("HPGeWaevfrmsFromStepPointMCs"); + log << "\n"; + log << "=====HPGeWaveformsFromStepPointMCs summary=====\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 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 6156d1fc6b..a976156ce1 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,17 +56,19 @@ 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, 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) : 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 +77,36 @@ 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 + ttree->Branch("SimParticleId", &simParticleId, "SimParticleId/i"); }; 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 +119,16 @@ 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 + simParticleId = particle.id().asInt(); + if (Ekin < 0) throw cet::exception("LogicError", "Energy is negative"); ttree->Fill(); diff --git a/STMReco/src/STMMovingWindowDeconvolution_module.cc b/STMReco/src/STMMovingWindowDeconvolution_module.cc index 6da8c05fc4..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,10 @@ 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; }; STMMovingWindowDeconvolution::STMMovingWindowDeconvolution(const Parameters& conf) : @@ -191,7 +202,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 +213,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 }; }; @@ -230,6 +243,7 @@ namespace mu2e { nADCs = ADCs.size(); processedWaveforms++; waveformID++; + processedADCs+=nADCs; if (M > nADCs) M = nADCs; @@ -240,11 +254,14 @@ namespace mu2e { differentiate(); average(); calculate_baseline(); + // CHANGE HERE + // 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 @@ -286,112 +303,190 @@ 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; - }; - }; - - 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::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; - }; - }; - - void STMMovingWindowDeconvolution::calculate_baseline() { + _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() { + double subtrahend = 0.0; + + 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; 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; - lowest_height = 0; - lowest_height_time = -1; // in clock ticks - - 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 + // 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 - 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 { + 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); // NOTE - THIS WAS FOR DEBUGGING ONLY + + baseline_stddev = calculated_sd; + } 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; + 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; + }; }; }; }; 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 << std::left << std::setw(25) << "\tnADCs: " << nADCs << "\n"; + log << std::left << std::setw(25) << "\tprocessedADCs: " << processedADCs << "\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; };