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/HPGeReco.fcl b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl similarity index 64% rename from STMMC/fcl/HPGeReco.fcl rename to STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl index b6dc9d5362..3c0058ef24 100644 --- a/STMMC/fcl/HPGeReco.fcl +++ b/STMMC/fcl/HPGeWaveformGenerationAndAnalysis.fcl @@ -1,4 +1,13 @@ - +# Executes the full HPGe waveform generation and analysis chain on STMMC data +# Outline: +# 1. DigiHPGe - generates HPGe waveforms from StepPointMCs +# 2. concatenateWaveformsHPGe - concatenates microspill waveforms into spill waveforms +# 3. (optional) ZSHPGe - zero-suppresses the concatenated waveforms +# 4. MWDHPGe - performs moving window deconvolution on the concatenated (or zero-suppressed) waveforms +# 5. Ou1tputs the MWD digis to a ROOT file +# Note - each module is documented in the _module.cc file. This is the clearest way to understand what each module does and what parameters it takes. +# Usage: +# 1. Select which data you want to run with or without zero suppression by replacing "@nil" in trigger_paths (keeping the [BRAKCETS]!) with the single analyer you want to run, e.g. trigger_paths: [digitization_path_ZS] # Original author: Pawel Plesniak #include "Offline/fcl/standardServices.fcl" @@ -17,6 +26,7 @@ services : @local::Services.Sim physics : { producers : { DigiHPGe : { + # Digitize the StepPointMCs from the HPGe crystals into waveforms module_type : HPGeWaveformsFromStepPointMCs StepPointMCsTagEle : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagEle StepPointMCsTagMu : @local::STMMCAnalysis.MixedEventsTags.StepPointMCsTagMu @@ -27,17 +37,20 @@ physics : { risingEdgeDecayConstant : @local::HPGeDigitization.DecayConstant microspillBufferLengthCount : @local::HPGeDigitization.microspillBufferLengthCount makeTTree : false + makeADCPlot : false timeOffset : 0 resetEventNumber : @local::HPGeDigitization.resetEventNumber verbosityLevel : 0 } concatenateWaveformsHPGe : { + # Concatenate the microspill waveforms into macrospill waveforms module_type : ConcatenateDigitizedWaveforms STMWaveformDigisTag : @local::HPGeDigitization.concatenation.STMWaveformDigisTag nMerge : @local::HPGeDigitization.concatenation.nMerge makeTTree : false } ZSHPGe : { + # Apply zero-suppression to the concatenated waveforms module_type : STMZeroSuppression stmWaveformDigisTag : @local::STMMCAnalysis.ZS.HPGe.stmWaveformDigisTag.concatenated tbefore : @local::STMMCAnalysis.ZS.HPGe.tbefore @@ -50,6 +63,7 @@ physics : { makeTTreeWaveforms: false } MWDHPGe : { + # Perform moving window deconvolution on the concatenated (or zero-suppressed) waveforms to determine the pulse energies module_type : STMMovingWindowDeconvolution stmWaveformDigisTag : @local::STMMCAnalysis.MWD.HPGe.stmWaveformDigisTag.concatenated tau : @local::STMMCAnalysis.MWD.HPGe.tau @@ -60,7 +74,7 @@ physics : { defaultBaselineMean : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineMean.suppressed defaultBaselineSD : @local::STMMCAnalysis.MWD.HPGe.defaultBaselineSD.suppressed makeTTreeMWD: false - makeTTreeEnergies: false + makeTTreeEnergies: true TTreeEnergyCalib : @local::HPGeDigitization.EnergyPerADCBin verbosityLevel : 0 xAxis : "" @@ -72,21 +86,22 @@ physics : { STMWaveformDigisTag : @local::HPGeDigitization.concatenation.filterTag } } - digitization_path : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] - trigger_paths : [digitization_path] - output_path : [compressedOutput] - end_paths : [output_path] + digitization_path_no_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, MWDHPGe] + digitization_path_ZS : [DigiHPGe, concatenateWaveformsHPGe, concatenationFilterHPGe, ZSHPGe, MWDHPGe] + trigger_paths : ["digitization_path_ZS"] # Populate me! TODO - make me @nil before git push + # output_path : [compressedOutput] + # end_paths : [output_path] } outputs : { compressedOutput : { module_type : RootOutput fileName : "dts.owner.HPGeReco.version.sequencer.art" - SelectEvents: ["digitization_path"] - outputCommands: [ - "drop *_*_*_*", - "keep mu2e::STMMWDDigis_MWDHPGe_*_HPGeReco" - ] + SelectEvents: [@nil] # Populate me! + # outputCommands: [ + # "drop *_*_*_*", + # "keep mu2e::STMMWDDigis_MWDHPGe_*_HPGeReco" + # ] } } diff --git a/STMMC/fcl/MakeTree.fcl b/STMMC/fcl/MakeTree.fcl deleted file mode 100644 index 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..41ffa70db8 --- /dev/null +++ b/STMMC/fcl/ROOTAnalysisDump.fcl @@ -0,0 +1,107 @@ +# Extracts parameters from STM simulation scripts into ROOT TTrees for use in the scripts defined it Mu2e/STMAnalysis repository +# This generates ROOT files containing different things depending on which analyzer is selected in the o1 path: +# - stage1virtualdetectorStepPointMCs : Dumps StepPointMCs from the virtual detector in Stage 1 simulation +# - stage2virtualdetectorStepPointMCs : Dumps StepPointMCs from the virtual detector in Stage 2 simulation +# - HPGeWaveforms : Dumps HPGe waveforms from Stage 2 simulation +# - LaBrWaveforms : Dumps LaBr3 waveforms from Stage 2 simulation (not implemented yet) +# - ZSwaveforms : Dumps zero-suppressed waveforms from HPGe digitisers in Stage 2 simulation (not implemented yet) +# - MWDspectra : Dumps MWD spectra from HPGe digitisers in Stage 2 simulation +# Usage: +# 1. Select which data you want to dump to a ROOT TTree by replacing "@nil" in the list o1 (keeping the [BRAKCETS]!) with the single analyer you want to run, e.g. o1: [stage1virtualdetectorStepPointMCs] +# 2. Set the output ROOT file name in services.TFileService.fileName +# 3. Prepare a text file containing a list of input STMMC data files (one file name per line) +# 3. Run as: `mu2e -c ROOTAnalysisDump.fcl -S FileWithListOfDataFiles.txt` +# +# Original author: Pawel Plesniak + +#include "Offline/STMMC/fcl/prolog.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name: ROOTAnalysisDump + +source : { + module_type : RootInput + fileNames: @nil +} +services : { + message : @local::default_message + GlobalConstantsService : { + inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt" + } +} +physics: { + analyzers : { + # Simulation stage 1 ROOT files + Stage1virtualdetectorStepPointMCs : { + module_type : VirtualDetectorTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage1.StepPointMCsTag # 1809 + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + } + Stage1SimParticeandVirtualDetectorBacktrace : { + module_type : SimParticleAndVDBacktrace + StepPointMCsTag : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + StartVirtualDetectorId : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + TraceVirtualdetectorIds : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage1 + } + + # Simulation stage 2 ROOT files + Stage2virtualdetectorStepPointMCs : { + module_type : VirtualDetectorTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage2.StepPointMCsTagVirtualdetector + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + } + Stage2SimParticeandVirtualDetectorBacktrace : { + module_type : SimParticleAndVDBacktrace + StepPointMCsTag : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + StartVirtualDetectorId : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + TraceVirtualdetectorIds : @local::ROOTAnalysisDump.SimParticleVirtualDetectorBacktrace.StepPointMCsTag.Stage2 + } + # TODO - see if this can be removed. + # ROOTDump : { + # module_type : SignalParticleBacktraceVDs + # # StepPointMCsTag : "compressDetStepMCsSTM:virtualdetector" #Stage 1 + # # # StepPointMCsTag: "compressDetStepMCsSTM:" #Stage 1 1809 as Stage 1 Mu3 + # # StartVirtualDetectorId : 101 + # # TrackVirtualDetectorIDs : [101, 100, 10, 9, 8, 2, 1] + # StepPointMCsTag : "compressSTMDet:virtualdetector" #Stage 2 + # StartVirtualDetectorId : 90 + # TrackVirtualDetectorIDs : [90, 100, 10, 9, 8, 2, 1] + # EnergyWindowSize : 0.1 + # } + + # HPGe waveform and analysis ROOT files + HPGeWaveforms : { + module_type : HPGeTree + StepPointMCsTag : @local::ROOTAnalysisDump.Stage2.StepPointMCsTagSTMDet + SimParticlemvTag : @local::ROOTAnalysisDump.Stage2.SimParticlemvTag + Detector: @local::ROOTAnalysisDump.Stage2.DetectorName.HPGe + } + HPGeZSWaveforms : { + module_type : @nil # Doesn't exist yet! + StepPointMCsTag : @local::SimplifyStage2Data.StepPointMCsTagSTMDet + SimParticlemvTag : @local::SimplifyStage2Data.SimParticlemvTag # TODO: Remove the SimParticle dependency + Detector: @local::SimplifyStage2Data.DetectorName.HPGe + } + HPGeMWDspectra : { + module_type : MWDTree + STMMWDDigiTag : @local::STMMCAnalysis.MWD.HPGe.STMMWDDigiTag + EnergyCalib : @local::HPGeDigitization.EnergyPerADCBin # Once the proditions are implemneted, this should be removed. + } + + # LaBr3 waveform ROOT files - not implemented yet + LaBrWaveforms : { + module_type : @nil + } + LaBrZSWaveforms : { + module_type : @nil + } + LaBrMWDspectra : { + module_type : @nil + } + } + + o1 : [@nil] # Populate me! + end_paths: [o1] +} + +services.TFileService.fileName : @nil # Populate me! diff --git a/STMMC/fcl/prolog.fcl b/STMMC/fcl/prolog.fcl index 1b56b4f7d1..9ba94f8829 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 : { + StepPointMCsTags : { + Stage1 : @local::DataProducts.Stage1.StepPointMCs + Stage2 : @local::DataProducts.Stage2.StepPointMCs + } + StartVirtualDetectorID : { + Stage1 : @local::DataProducts.Stage1.EndVirtualdetectorID + Stage2 : { + HPGe : @local::DataProducts.Stage2.EndVirtualdetectorID.HPGe + LaBr : @local::DataProducts.Stage2.EndVirtualdetectorID.LaBr + } + } + TraceVirtualdetectorIds : { + Stage1: [101, 100, 10, 9, 8] + Stage2: [90, 101, 100, 10, 9, 8] + } + } } -SimplifyStage2Data : { - StepPointMCsTagSTMDet : "compressSTMDet:STMDet" - StepPointMCsTagVirtualdetector : "compressSTMDet:virtualdetector" - SimParticlemvTag : "compressSTMDet:" - VirtualDetectorID : { - HPGe : 90 - LaBr : 89 - } - DetectorName : { - HPGe : "HPGe" - LaBr : "LaBr" + +# For analysis studies, you should only need to edit lines between the "edit lines" comments +######################################## Start of edit lines ######################################## + +# Mixing parameters +MixSTMEvents : { + extendedMean2BB : 3.93e7 # Copied from Production/JobConfig/mixing/TwoBB.fcl + cutMax2BB : 2.36e8 # Copied from Production/JobConfig/mixing/TwoBB.fcl + extendedMean1BB : 1.58e7 # Copied from Production/JobConfig/mixing/OneBB.fcl + cutMax1BB : 9.48e7 # Copied from Production/JobConfig/mixing/OneBB.fcl + SDF : 0.6 # Copied from Production/JobConfig/mixing/OneBB.fcl + nPOTs : 26919873604557116 # Effective number of POTs, correcting for the resampling factor + nMicroSpills : 690253169 + meanEventsPerPOTFactors : { + EleBeamCat : 5.465876332164181e-11 + MuBeamCat : 7.887481312841522e-13 + TargetStopsCat1809 : 3.114418308224206e-16 } } @@ -174,25 +255,70 @@ STMDAQParameters : { LaBr : 370.370370370 # MSPS - copied from Offline/STMConditions/fcl/prolog.fcl } } + +# HPGe waveform generation parameters HPGeDigitization : { - EnergyPerADCBin : 0.5 # keV/bin - PreamplifierNoiseSD : 0.32 # mV + EnergyPerADCBin : 1 # keV/bin + PreamplifierNoiseSD : 0.0 # mV DecayConstant : 50 # us - resetEventNumber : 31858 + resetEventNumber : @local::DataProducts.nMicrospillsPerSpill microspillBufferLengthCount : 1000 concatenation : { - STMWaveformDigisTag : "DigiHPGe:" - nMerge : 31858 - filterTag : "concatenateWaveformsHPGe:" + STMWaveformDigisTag : @local::DataProducts.Waveforms.HPGe.uSpill + nMerge : @local::DataProducts.nMicrospillsPerSpill + filterTag : @local::DataProducts.Waveforms.HPGe.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 + } + } } +######################################### End of edit lines ######################################### END_PROLOG diff --git a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc index f084862a4e..413151ad00 100644 --- a/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc +++ b/STMMC/src/HPGeWaveformsFromStepPointMCs_module.cc @@ -1,6 +1,11 @@ // Simulates the electronics response of the HPGe detector. Simulates the pulse height, decay tail, and ADC digitization. Generates one STMWaveformDigi per micropulse. // Model based heavily on example provided in docDb 43617 // See docDb 51487 for full documentation +// Remaining TODOs for future development: +// - Include the commented out includes correctly +// - Implement the makeADCPlot functionality correctly +// - Implement different sampling frequencies correctly +// - Remove the implementation of stepPositionTolerance, it is a hack to address a position resolution I did not have time to fully address // Original author: Pawel Plesniak // stdlib includes @@ -45,6 +50,7 @@ // ROOT includes #include "art_root_io/TFileService.h" #include "TTree.h" +#include "TGraph.h" namespace mu2e { @@ -62,6 +68,7 @@ namespace mu2e { fhicl::Atom risingEdgeDecayConstant{ Name("risingEdgeDecayConstant"), Comment("Rising edge decay time [us]")}; fhicl::OptionalAtom microspillBufferLengthCount{ Name("microspillBufferLengthCount"), Comment("Number of microspills to buffer ahead for, in number of microspills")}; fhicl::OptionalAtom makeTTree{ Name("makeTTree"), Comment("Controls whether to make the TTree with branches chargeCollected, chargeDecayed, ADC, eventId, time")}; + fhicl::OptionalAtom makeADCPlot{ Name("makeADCPlot"), Comment("Controls whether to make a plot of the generated ADC values")}; fhicl::OptionalAtom timeOffset{ Name("timeOffset"), Comment("For debugging, adds the named time offset in [ns], used for testing analysis algorithms")}; fhicl::OptionalAtom resetEventNumber{ Name("resetEventNumber"), Comment("Simulates the off-spill period by resetting the inter-event last decayed charge to zero")}; fhicl::OptionalAtom verbosityLevel{ Name("verbosityLevel"), Comment("Controls verbosity")}; @@ -75,6 +82,7 @@ namespace mu2e { void decayCharge(); void addNoise(); void digitize(); + void endJob(); // fhicl variables art::ProductToken StepPointMCsTokenEle, StepPointMCsTokenMu, StepPointMCsToken1809; // Token of StepPointMCs in STMDet @@ -83,6 +91,7 @@ namespace mu2e { double noiseSD = 0; // Standard deviation of ADC noise [mV] double risingEdgeDecayConstant = 0; // [us] bool makeTTree = false; // Controls whether an analysis TTree is made + bool makeADCPlot = false; // Controls whether the ADC graph is made double timeOffset = 0.0; // Used for debugging [ns] uint resetEventNumber = 0; // Event ID at which to reset the charge amplitude int verbosityLevel = 0; // How much output to generate @@ -155,6 +164,10 @@ namespace mu2e { std::vector _chargeCarryOver; // Temporary buffer that will store _chargeCollected over the course of the next event std::vector _adcs; // Buffer for storing the ADC values to put into the STMWaveformDigi + // Debugging variables + double waveform_time = 0.0; // Time of start of waveform for debugging. TODO - REMOVE ME! + int32_t kept_events = 0, dropped_events = 0; // Counters for debugging + // Offline utilities // TODO: include the prodition to get the sampling frequency // mu2e::STMChannel::enum_type _HPGeChannel = static_cast(1); @@ -205,7 +218,7 @@ namespace mu2e { // Assign the approrpiate time offset timeOffset = conf().timeOffset() ? *(conf().timeOffset()) : 0.0; - // Assign TTrees + // Assign optional ROOT parameters makeTTree = conf().makeTTree() ? *(conf().makeTTree()) : false; if (makeTTree) { art::ServiceHandle tfs; @@ -216,7 +229,15 @@ namespace mu2e { ttree->Branch("eventId", &eventId, "eventId/i"); ttree->Branch("time", &time, "time/i"); }; + makeADCPlot = conf().makeADCPlot() ? *(conf().makeADCPlot()) : false; + if (makeADCPlot) { + throw cet::exception("IMPLEMENT", "makeADCPlot functionality not yet implemented correctly\n"); + art::ServiceHandle tfs; + TTree* adcTree = tfs->make("adcTree", "HPGeWaveformsFromStepPointMCs ADC Tree"); + adcTree->Branch("ADC", &ADC, "ADC/D"); + }; + // Assign optional variables resetEventNumber = conf().resetEventNumber() ? *(conf().resetEventNumber()) : 0; }; @@ -230,6 +251,7 @@ namespace mu2e { std::cout << std::left << "\t\t" << std::setw(60) << "risingEdgeDecayConstant [us]" << risingEdgeDecayConstant << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "microspillBufferLengthCount" << microspillBufferLengthCount << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "makeTTree" << makeTTree << std::endl; + std::cout << std::left << "\t\t" << std::setw(60) << "makeADCPlot" << makeADCPlot << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "timeOffset [ns]" << timeOffset << std::endl; std::cout << std::left << "\t\t" << std::setw(60) << "resetEventNumber" << resetEventNumber << std::endl; std::cout << "\tDerived parameters: " << std::endl; @@ -249,36 +271,92 @@ namespace mu2e { void HPGeWaveformsFromStepPointMCs::produce(art::Event& event) { eventId = event.id().event(); + // Simulation takes the POT time as t = 0, and has sequential microspills (events). The trigger time offset is not used here, left as a TODO + // Create the STMWaveformDigi and insert all the relevant attributes + // TODO - this only keeps accurate time if the sampling is 320MHz. Needs to be rewritten to work for other times + if (std::abs(fADC - 320) > 1e-12) + throw cet::exception("RANGE") << "Currently only fADC of 320 MHz is supported (got " << fADC << " MHz)\n"; + nADCs = 542; + + // Assign the variables for time and number of ADCs + eventTimeBuffer = eventId % 5; + if (!(eventTimeBuffer == 0 || eventTimeBuffer == 3)) + nADCs++; + eventTime += nADCs; + + // Update the last event decayed charge before the noise so the noise isn't added twice + if (resetEventNumber != 0 && eventId == resetEventNumber) { + lastEventEndDecayedCharge = 0; + std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); + } + else + lastEventEndDecayedCharge = _chargeDecayed.back(); + + // Update the parameters to carry over to the next event + _chargeCollected.clear(); + _chargeCollected.assign(_chargeCarryOver.begin(), _chargeCarryOver.end()); + _chargeCollected.insert(_chargeCollected.end(), nADCs, 0); + _chargeCarryOver.clear(); + _chargeCarryOver.assign(_chargeCollected.begin() + nADCs, _chargeCollected.end()); + + // Clear previous buffer vectors + _chargeDecayed.clear(); + _chargeDecayed.insert(_chargeDecayed.end(), nADCs, 0); + _adcs.clear(); + _adcs.insert(_adcs.end(), nADCs, 0); + + // Reset waveform time for debugging + waveform_time = 0.0; + // Get the hits in the detector std::vector StepsEle = event.getProduct(StepPointMCsTokenEle); std::vector StepsMu = event.getProduct(StepPointMCsTokenMu); std::vector Steps1809 = event.getProduct(StepPointMCsToken1809); + // If there are no hits, produce an empty STMWaveformDigiCollection and return + if (StepsEle.size() == 0 && StepsMu.size() == 0 && Steps1809.size() == 0) { + std::vector emptyAdcs(nADCs, static_cast(0)); + STMWaveformDigi _waveformDigi(0.0, emptyAdcs); + std::unique_ptr outputDigis(new STMWaveformDigiCollection); + outputDigis->emplace_back(_waveformDigi); + event.put(std::move(outputDigis)); + dropped_events++; + return; + } + kept_events++; + // Add a collection of charge depositions to _charge for(const StepPointMC& step : StepsEle){ - if (step.ionizingEdep() != 0) - depositCharge(step); + if (step.ionizingEdep() != 0) { + depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); + }; }; + // Process the collected charge for each step + for (const StepPointMC& step : StepsEle) { + if (step.ionizingEdep() != 0) { + depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); + } + } for(const StepPointMC& step : StepsMu){ if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); }; for(const StepPointMC& step : Steps1809){ if (step.ionizingEdep() != 0) depositCharge(step); + if (step.time() < waveform_time) + waveform_time = step.time(); }; // Decay all of the collected charges decayCharge(); - // Update the last event decayed charge before the noise so the noise isn't added twice - if (resetEventNumber != 0 && eventId == resetEventNumber) { - lastEventEndDecayedCharge = 0; - std::fill(_chargeCarryOver.begin(), _chargeCarryOver.end(), 0); - } - else - lastEventEndDecayedCharge = _chargeDecayed.back(); - // Add preamplifier electronics noise with SD defined in noiseSD addNoise(); @@ -291,15 +369,8 @@ namespace mu2e { throw cet::exception("LogicError", "ADC values too high!"); }; - // Simulation takes the POT time as t = 0, and has sequential microspills (events). The trigger time offset is not used here, left as a TODO - // Create the STMWaveformDigi and insert all the relevant attributes - // TODO - this only keeps accurate time if the sampling is 320MHz. Needs to be rewritten to work for other times - eventTimeBuffer = eventId % 5; - if (eventTimeBuffer == 0 || (eventId % 3) == 0) - eventTime += nADCs + 1; - else - eventTime += nADCs; - STMWaveformDigi _waveformDigi(eventTime, _adcs); + // STMWaveformDigi _waveformDigi(eventTime, _adcs); + STMWaveformDigi _waveformDigi(waveform_time, _adcs); std::unique_ptr outputDigis(new STMWaveformDigiCollection); outputDigis->emplace_back(_waveformDigi); @@ -315,15 +386,21 @@ namespace mu2e { }; }; - // Update the parameters to carry over to the next event - _chargeCarryOver.clear(); - _chargeCarryOver.assign(_chargeCollected.begin() + nADCs, _chargeCollected.end()); - _chargeCollected.clear(); - _chargeCollected.assign(_chargeCarryOver.begin(), _chargeCarryOver.end()); - _chargeCollected.insert(_chargeCollected.end(), nADCs, 0); - // Clear previous buffer vectors - std::fill(_chargeDecayed.begin(), _chargeDecayed.end(), 0); - std::fill(_adcs.begin(), _adcs.end(), 0); + // Make the ADC plot if appropriate + if (makeADCPlot) { + // NOTE - THIS DOES NOT CURRENTLY FUNCTION, AND I DO NOT HAVE TIME TO FIX IT RIGHT NOW + std::stringstream histsuffix; + histsuffix.str(""); + histsuffix << "_evt" << event.event(); + + art::ServiceHandle tfs; + std::vector _adc_double = std::vector(_adcs.begin(), _adcs.end()); + TGraph* g_adcs = tfs->make(nADCs, _adc_double.data()); + g_adcs->SetName(("g_adcs"+histsuffix.str()).c_str()); + // g_adcs->SetTitle("ADC Values"); + // for (uint i = 0; i < nADCs; i++) + // g_adcs->SetPoint(i, i, _adcs[i]); + }; // Add the STMWaveformDigi to the event event.put(std::move(outputDigis)); @@ -480,6 +557,14 @@ namespace mu2e { }; return; }; + + void HPGeWaveformsFromStepPointMCs::endJob() { + mf::LogInfo log("STMResamplingFilter summary"); + log << "=====HPGeWaveformsFromStepPointMCs summary=====\n"; + log << "No. kept events: " << kept_events << "\n"; + log << "No. discarded events: " << dropped_events << "\n"; + log << "===============================================\n"; + }; }; // namespace mu2e DEFINE_ART_MODULE(mu2e::HPGeWaveformsFromStepPointMCs) diff --git a/STMMC/src/STMResamplingProducer_module.cc b/STMMC/src/STMResamplingProducer_module.cc index 02c83d003f..30137f61df 100644 --- a/STMMC/src/STMResamplingProducer_module.cc +++ b/STMMC/src/STMResamplingProducer_module.cc @@ -44,7 +44,7 @@ namespace mu2e { art::EDProducer{conf}, StepPointMCsToken(consumes(conf().stepPointMCsTag())), virtualDetectorID(conf().virtualDetectorID()) { - produces(); + produces("virtualdetector"); }; void STMResamplingProducer::produce(art::Event& event) { @@ -63,7 +63,7 @@ namespace mu2e { // Update counter includedStepPointMCs += outputStepPointMCs->size(); // Add the new data products to the event - event.put(std::move(outputStepPointMCs)); + event.put(std::move(outputStepPointMCs), "virtualdetector"); return; }; diff --git a/STMMC/src/SimParticleAndVDBacktrace_module.cc b/STMMC/src/SimParticleAndVDBacktrace_module.cc new file mode 100644 index 0000000000..82f4538588 --- /dev/null +++ b/STMMC/src/SimParticleAndVDBacktrace_module.cc @@ -0,0 +1,263 @@ +// Adapted from ReadVirtualDetector_module.cc +// Generates a large ROOT TTree with SimParticle and StepPointMC information. For each SimParticle that has a StepPointMC in the specified virtual detector, gets its +// - parentCounter - 0 for the original particle, 1 for its parent, 2 for its grandparent, etc. This is traced all the way back to the POT +// - particleId - SimParticle ID (from GEANT) +// - pdgId - PDG ID +// - creationCode - creation code. See Offline/MCDataProducts/inc/ProcessCode.hh for details +// - x [mm] - start x position of the SimParticle +// - y [mm] - start y position of the SimParticle +// - z [mm] - start z position of the SimParticle +// - px [MeV/c] - start x momentum of the SimParticle +// - py [MeV/c] - start y momentum of the SimParticle +// - pz [MeV/c] - start z momentum of the SimParticle +// - time [ns] - start global time of the SimParticle +// - Ekin [MeV] - end kinetic energy of the SimParticle +// For each chosen virtual detector ID to trace back to, also gets the following StepPointMC information (0 if the SimParticle did not have a StepPointMC in that virtual detector): +// - VD__x [mm] - x position of the StepPointMC in that virtual detector +// - VD__y [mm] - y position of the StepPointMC in that virtual detector +// - VD__z [mm] - z position of the StepPointMC in that virtual detector +// - VD__px [MeV/c] - x momentum of the StepPointMC in that virtual detector +// - VD__py [MeV/c] - y momentum of the StepPointMC in that virtual detector +// - VD__pz [MeV/c] - z momentum of the StepPointMC in that virtual detector +// - VD__time [ns] - time of the StepPointMC in that virtual detector +// - VD__Ekin [MeV] - kinetic energy of the StepPointMC in that virtual detector +// - VD__pdgId - PDG ID of the SimParticle at that StepPointMC +// Original author: Ivan Logashenko +// Adapted by: Pawel Plesniak + +// stdlib includes +#include +#include + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" + +// exception handling +#include "cetlib_except/exception.h" + +// fhicl includes +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" + +// message handling +#include "messagefacility/MessageLogger/MessageLogger.h" + +// Offline includes +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/StepPointMC.hh" + +// ROOT includes +#include "art_root_io/TFileService.h" +#include "TTree.h" + +typedef cet::map_vector_key key_type; + +namespace mu2e { + class SimParticleAndVDBacktrace : public art::EDAnalyzer { + public: + using Name=fhicl::Name; + using Comment=fhicl::Comment; + struct Config { + fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; + fhicl::Atom StartVirtualDetectorId{Name("StartVirtualDetectorId"), Comment("ID of the virtual detector to filter")}; + fhicl::Sequence TraceVirtualdetectorIds{Name("TraceVirtualdetectorIds"), Comment("IDs of the virtual detectors to trace back to the origin")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit SimParticleAndVDBacktrace(const Parameters& conf); + void analyze(const art::Event& e); + void addToTree(const art::Ptr &particle, const std::vector steps); + void addParentToTree(const art::Ptr &particle, const std::vector steps); + void traceSteps(const art::Ptr &particle, const std::vector steps); + void clear_virtualdetector_vectors(); + void populate_virtualdetector_vector(const StepPointMC & step_it, const unsigned index); + private: + art::ProductToken StepPointMCsToken; + unsigned StartVirtualDetectorId = 0; + std::vector TraceVirtualdetectorIds; + int consecutiveEmptyFileThreshold = 0; + + GlobalConstantsHandle pdt; + int pdgId = 0, creationCode = 0, consecutiveEmptyFileCounter = 0, nTraceVirtualdetectorIds = 0; + double x = 0.0, y = 0.0, z = 0.0, px = 0.0, py = 0.0, pz = 0.0, time = 0.0, Ekin = 0.0; + std::vector VD_x, VD_y, VD_z, VD_px, VD_py, VD_pz, VD_time, VD_Ekin; + std::vector VD_pdgId; + uint16_t parentCounter = 0; + cet::map_vector_key particleId; + TTree* ttree; + CLHEP::Hep3Vector startPosition, startMomentum; + art::Ptr parent; + std::vector*> VD_vectors; + }; + + SimParticleAndVDBacktrace::SimParticleAndVDBacktrace(const Parameters& conf) : + art::EDAnalyzer(conf), + StepPointMCsToken(consumes(conf().StepPointMCsTag())), + StartVirtualDetectorId(conf().StartVirtualDetectorId()), + TraceVirtualdetectorIds(conf().TraceVirtualdetectorIds()) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; + if (std::find(TraceVirtualdetectorIds.begin(), TraceVirtualdetectorIds.end(), StartVirtualDetectorId) == TraceVirtualdetectorIds.end()) { + TraceVirtualdetectorIds.push_back(StartVirtualDetectorId); + }; + nTraceVirtualdetectorIds = TraceVirtualdetectorIds.size(); + VD_vectors = {&VD_x, &VD_y, &VD_z, &VD_px, &VD_py, &VD_pz, &VD_time, &VD_Ekin}; + for (auto vec : VD_vectors) { + vec->resize(nTraceVirtualdetectorIds, 0.0); + }; + VD_pdgId.resize(nTraceVirtualdetectorIds, 0); + + art::ServiceHandle tfs; + ttree = tfs->make( "ttree", "SimParticle ttree"); + ttree->Branch("parentCounter", &parentCounter, "parentCounter/s"); + ttree->Branch("particleId", &particleId, "particleId/l"); + ttree->Branch("pdgId", &pdgId, "pdgId/I"); + ttree->Branch("creationCode", &creationCode, "creationCode/I"); + ttree->Branch("x", &x, "x/D"); // mm + ttree->Branch("y", &y, "y/D"); // mm + ttree->Branch("z", &z, "z/D"); // mm + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("time", &time, "time/D"); // ns + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + for (size_t i=0; iBranch(Form("VD_%u_x", TraceVirtualdetectorIds[i]), &VD_x[i], Form("VD_%u_x/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_y", TraceVirtualdetectorIds[i]), &VD_y[i], Form("VD_%u_y/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_z", TraceVirtualdetectorIds[i]), &VD_z[i], Form("VD_%u_z/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_px", TraceVirtualdetectorIds[i]), &VD_px[i], Form("VD_%u_px/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_py", TraceVirtualdetectorIds[i]), &VD_py[i], Form("VD_%u_py/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_pz", TraceVirtualdetectorIds[i]), &VD_pz[i], Form("VD_%u_pz/D", TraceVirtualdetectorIds[i])); // mm + ttree->Branch(Form("VD_%u_time", TraceVirtualdetectorIds[i]), &VD_time[i], Form("VD_%u_time/D", TraceVirtualdetectorIds[i])); // ns + ttree->Branch(Form("VD_%u_Ekin", TraceVirtualdetectorIds[i]), &VD_Ekin[i], Form("VD_%u_Ekin/D", TraceVirtualdetectorIds[i])); // MeV + ttree->Branch(Form("VD_%u_pdgId", TraceVirtualdetectorIds[i]), &VD_pdgId[i], Form("VD_%u_pdgId/I", TraceVirtualdetectorIds[i])); + }; + }; + + void SimParticleAndVDBacktrace::populate_virtualdetector_vector(const StepPointMC & step_it, const unsigned index) { + const CLHEP::Hep3Vector momentum = step_it.momentum(); + const CLHEP::Hep3Vector position = step_it.position(); + const double mass = pdt->particle(step_it.simParticle()->pdgId()).mass(); + const double EKin = std::sqrt(momentum.mag2() + mass * mass) - mass; + VD_x[index] = position.x(); + VD_y[index] = position.y(); + VD_z[index] = position.z(); + VD_px[index] = momentum.x(); + VD_py[index] = momentum.y(); + VD_pz[index] = momentum.z(); + VD_time[index] = step_it.time(); + VD_Ekin[index] = EKin; + VD_pdgId[index] = step_it.simParticle()->pdgId(); + }; + + void SimParticleAndVDBacktrace::traceSteps(const art::Ptr &particle, const std::vector steps) { + unsigned step_virtualdetectorId = 0, index = 0; + for (const StepPointMC & step_it : steps) { + if (step_it.simParticle()->id() != particle->id()) continue; + step_virtualdetectorId = step_it.virtualDetectorId(); + auto it = std::find(TraceVirtualdetectorIds.begin(), TraceVirtualdetectorIds.end(), step_virtualdetectorId); + if (it == TraceVirtualdetectorIds.end()) continue; + index = std::distance(TraceVirtualdetectorIds.begin(), it); + populate_virtualdetector_vector(step_it, index); + }; + }; + + void SimParticleAndVDBacktrace::clear_virtualdetector_vectors() { + for (auto & vec : VD_vectors) { + std::fill(vec->begin(), vec->end(), 0.0); + }; + std::fill(VD_pdgId.begin(), VD_pdgId.end(), 0); + return; + }; + + void SimParticleAndVDBacktrace::addToTree(const art::Ptr &particle, const std::vector steps) { + // Assigne all the SimParticle variables to the tree variables + parentCounter = 0; + particleId = particle->id(); + pdgId = particle->pdgId(); + creationCode = particle->creationCode(); + startPosition = particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle->startGlobalTime(); + Ekin = particle->endKineticEnergy(); + + // Trace the StepPointMCs in the relevant virtualdetectors + clear_virtualdetector_vectors(); + traceSteps(particle, steps); + ttree->Fill(); + + // Now trace the parents + if (!particle->hasParent()) return; + + // Require separate definitions of addParentToTree because of data accessors + parent = particle->parent(); + addParentToTree(parent, steps); + + while (parent->hasParent()) { + parent = parent->parent(); + addParentToTree(parent, steps); + }; + }; + + void SimParticleAndVDBacktrace::addParentToTree(const art::Ptr &particle, const std::vector steps) { + // Increment the parent counter + parentCounter++; + + // Assign all the SimParticle variables to the tree variables + particleId = particle->id(); + pdgId = particle->pdgId(); + creationCode = particle->creationCode(); + startPosition = particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle->startGlobalTime(); + Ekin = particle->endKineticEnergy(); + + // Trace the StepPointMCs in the relevant virtualdetectors + clear_virtualdetector_vectors(); + traceSteps(particle, steps); + ttree->Fill(); + }; + + void SimParticleAndVDBacktrace::analyze(const art::Event& event) { + // Get the data products + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, are you sure you have the correct data product name?\n"); + }; + return; + }; + const std::vector StepPointMCs = *stepHandle; + consecutiveEmptyFileCounter = 0; + + for (const StepPointMC& step : StepPointMCs) { + // Filter by virtual detector ID if requested + if (step.virtualDetectorId() != StartVirtualDetectorId) continue; + + // Get the associated particle + const art::Ptr particle = step.simParticle(); + addToTree(particle, StepPointMCs); + }; + return; + }; +}; // end namespace mu2e + +DEFINE_ART_MODULE(mu2e::SimParticleAndVDBacktrace) diff --git a/STMMC/src/SimParticleDump_module.cc b/STMMC/src/SimParticleDump_module.cc new file mode 100644 index 0000000000..48cf595cf6 --- /dev/null +++ b/STMMC/src/SimParticleDump_module.cc @@ -0,0 +1,182 @@ +// Adapted from ReadVirtualDetector_module.cc +// For StepPointMCs in virtualdetectors, generates a TTree showing the provenance of the generated SimParticles with the following branches +// - pdgId - PDG ID +// - creationCode - creation code. See Offline/MCDataProducts/inc/ProcessCode.hh for details +// - x [mm] - start x position of the SimParticle +// - y [mm] - start y position of the SimParticle +// - z [mm] - start z position of the SimParticle +// - px [MeV/c] - start x momentum of the SimParticle +// - py [MeV/c] - start y momentum of the SimParticle +// - pz [MeV/c] - start z momentum of the SimParticle +// - time [ns] - start global time of the SimParticle +// - Ekin [MeV] - end kinetic energy of the SimParticle +// Original author: Ivan Logashenko +// Adapted by: Pawel Plesniak + +// TODO before upload - validate that the code works and generates the same particle trace as the SimParticleAndVDBacktrace module + +// stdlib includes +#include +#include + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" + +// exception handling +#include "cetlib_except/exception.h" + +// fhicl includes +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" + +// message handling +#include "messagefacility/MessageLogger/MessageLogger.h" + +// Offline includes +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" +#include "Offline/GlobalConstantsService/inc/ParticleDataList.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/StepPointMC.hh" + +// ROOT includes +#include "art_root_io/TFileService.h" +#include "TTree.h" + +typedef cet::map_vector_key key_type; + +namespace mu2e { + class SimParticleDump : public art::EDAnalyzer { + public: + using Name=fhicl::Name; + using Comment=fhicl::Comment; + struct Config { + fhicl::Atom StepPointMCsTag{Name("StepPointMCsTag"), Comment("Tag identifying the StepPointMCs")}; + fhicl::Atom SimParticlemvTag{Name("SimParticlemvTag"), Comment("Tag identifying the SimParticlemv")}; + fhicl::Atom FilterVirtualDetectorId{Name("FilterVirtualDetectorId"), Comment("ID of the virtual detector to filter")}; + fhicl::OptionalAtom consecutiveEmptyFileThreshold{Name("consecutiveEmptyFileThreshold"), Comment("Number of consecutive empty files before stopping the job")}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit SimParticleDump(const Parameters& conf); + void analyze(const art::Event& e); + void addToTree(const SimParticle& particle, TTree* ttree); + void addParentToTree(art::Ptr &_particle, TTree* ttree); + private: + art::ProductToken StepPointMCsToken; + art::ProductToken SimParticlemvToken; + GlobalConstantsHandle pdt; + int pdgId = 0, creationCode = 0, consecutiveEmptyFileCounter = 0, consecutiveEmptyFileThreshold = 0; + double x = 0.0, y = 0.0, z = 0.0, px = 0.0, py = 0.0, pz = 0.0, time = 0.0, Ekin = 0.0; + unsigned virtualdetectorId = 0, FilterVirtualDetectorId = 0; + bool forward = false; + uint16_t parentCounter = 0; + cet::map_vector_key particleId; + TTree* ttree; + CLHEP::Hep3Vector startPosition, startMomentum; + art::Ptr parent; + }; + + SimParticleDump::SimParticleDump(const Parameters& conf) : + art::EDAnalyzer(conf), + StepPointMCsToken(consumes(conf().StepPointMCsTag())), + SimParticlemvToken(consumes(conf().SimParticlemvTag())), + FilterVirtualDetectorId(conf().FilterVirtualDetectorId()) { + consecutiveEmptyFileThreshold = conf().consecutiveEmptyFileThreshold() ? *(conf().consecutiveEmptyFileThreshold()) : 10; + art::ServiceHandle tfs; + ttree = tfs->make( "ttree", "SimParticle ttree"); + ttree->Branch("parentCounter", &parentCounter, "parentCounter/s"); + ttree->Branch("particleId", &particleId, "particleId/l"); + ttree->Branch("pdgId", &pdgId, "pdgId/I"); + ttree->Branch("creationCode", &creationCode, "creationCode/I"); + ttree->Branch("x", &x, "x/D"); // mm + ttree->Branch("y", &y, "y/D"); // mm + ttree->Branch("z", &z, "z/D"); // mm + ttree->Branch("px", &px, "px/D"); // mm + ttree->Branch("py", &py, "py/D"); // mm + ttree->Branch("pz", &pz, "pz/D"); // mm + ttree->Branch("time", &time, "time/D"); // ns + ttree->Branch("Ekin", &Ekin, "Ekin/D"); // MeV + }; + + void SimParticleDump::addToTree(const SimParticle& particle, TTree* ttree) { + parentCounter = 0; + particleId = particle.id(); + pdgId = particle.pdgId(); + creationCode = particle.creationCode(); + startPosition = particle.startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = particle.startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = particle.startGlobalTime(); + Ekin = particle.endKineticEnergy(); + ttree->Fill(); + + if (!particle.hasParent()) { + return; + }; + parentCounter++; + parent = particle.parent(); + addParentToTree(parent, ttree); + while (parent->hasParent()) { + parentCounter++; + parent = parent->parent(); + addParentToTree(parent, ttree); + }; + }; + + void SimParticleDump::addParentToTree(art::Ptr &_particle, TTree* ttree) { + particleId = _particle->id(); + pdgId = _particle->pdgId(); + creationCode = _particle->creationCode(); + startPosition = _particle->startPosition(); + x = startPosition.x(); + y = startPosition.y(); + z = startPosition.z(); + startMomentum = _particle->startMomentum(); + px = startMomentum.x(); + py = startMomentum.y(); + pz = startMomentum.z(); + time = _particle->startGlobalTime(); + Ekin = _particle->endKineticEnergy(); + ttree->Fill(); + }; + + void SimParticleDump::analyze(const art::Event& event) { + auto stepHandle = event.getHandle< std::vector >(StepPointMCsToken); + if (!stepHandle || stepHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + }; + auto simHandle = event.getHandle< SimParticleCollection >(SimParticlemvToken); + if (!simHandle || simHandle->empty()) { + consecutiveEmptyFileCounter++; + return; + }; + if (consecutiveEmptyFileCounter > consecutiveEmptyFileThreshold) { + throw cet::exception("LogicError", "Too many consecutive empty files, stopping the job"); + }; + + auto const& StepPointMCs = *stepHandle; + auto const& SimParticles = *simHandle; + consecutiveEmptyFileCounter = 0; + + for (const StepPointMC& step : StepPointMCs) { + // Filter by virtual detector ID if requested + if (step.virtualDetectorId() != FilterVirtualDetectorId) continue; + + // Get the associated particle + const SimParticle& particle = SimParticles.at(step.trackId()); + addToTree(particle, ttree); + }; + return; + }; +}; // end namespace mu2e + +DEFINE_ART_MODULE(mu2e::SimParticleDump) diff --git a/STMMC/src/VirtualDetectorTree_module.cc b/STMMC/src/VirtualDetectorTree_module.cc index 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();