Updates to the STM simulation framework

Analyses/src/CountVirtualDetectorHits_module.cc

@@ -51,10 +51,17 @@ namespace mu2e {
     : art::EDAnalyzer(conf),
       StepPointMCsToken(consumes<StepPointMCCollection>(conf().stepPointMCsTag())),
       enabledVDs(conf().enabledVDs()) {
-        // Create list of unique enabled virtual detectors
-        std::set<int> enabledVDsSet(enabledVDs.begin(), enabledVDs.end());
+        // Create list of unique enabled virtual detectors, preserving the order
+        std::vector<int> unqiueEnabledVDsVec;;
+        std::unordered_set<int> 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();

EventGenerator/src/PhotonGun_module.cc

@@ -4,6 +4,7 @@
 
 // stdlib includes
 #include <cmath>
+#include <iostream>
 
 // art includes
 #include "art/Framework/Core/EDProducer.h"
@@ -38,15 +39,21 @@ namespace mu2e {
       fhicl::Atom<double> z{ Name("z"), Comment("z position of generated photon")};
       fhicl::OptionalAtom<double> px{ Name("px"), Comment("x momentum of generated photon")};
       fhicl::OptionalAtom<double> py{ Name("py"), Comment("y momentum of generated photon")};
+      fhicl::OptionalAtom<double> deltax{ Name("deltax"), Comment("Difference in x position of generated photon, use to calculate targeted position")};
+      fhicl::OptionalAtom<double> deltay{ Name("deltay"), Comment("Difference in y position of generated photon, use to calculate targeted position")};
+      fhicl::OptionalAtom<double> deltaz{ Name("deltaz"), Comment("Difference in z position of generated photon, use to calculate targeted position")};
       fhicl::Atom<double> E{ Name("E"), Comment("Energy of generated photon")};
+      fhicl::OptionalAtom<bool> verbose{ Name("verbose"), Comment("Print verbose messages")};
     };
     using Parameters = art::EDProducer::Table<Config>;
     explicit PhotonGun(const Parameters& conf);
     virtual void produce(art::Event& event);
   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<GenParticleCollection>();
+      if (E < std::numeric_limits<double>::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<double>::epsilon() || std::abs(py) > std::numeric_limits<double>::epsilon()) && (std::abs(deltax) > std::numeric_limits<double>::epsilon() || std::abs(deltay) > std::numeric_limits<double>::epsilon() || std::abs(deltaz) > std::numeric_limits<double>::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<GenParticleCollection> 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<double>::epsilon() ||
+        std::abs(deltay) > std::numeric_limits<double>::epsilon() ||
+        std::abs(deltaz) > std::numeric_limits<double>::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));

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<int16_t>()) {};
+    STMWaveformDigi() : _DetID(0), _EWT(0), _DTCtime(0), _ADCtime(0), _trigTimeOffset(0), _peak_fitTime1(0), _peak_fitTime2(0), _peak_sep(0), _adcs(std::vector<int16_t>()), _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<int16_t> &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<int16_t> &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<int16_t> &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<int16_t> &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<int16_t> &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<int16_t> &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<int16_t> &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<int16_t>& 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<int16_t> _adcs; // vector of ADC values for the waveform
+    double   _tmp_sim_time; // temporary storage of sim time for validation (ns)
   };
 
   typedef std::vector<STMWaveformDigi> STMWaveformDigiCollection;

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

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

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

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

STMMC/fcl/HPGeReco.fcl → 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"
+    # ]
   }
 }