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/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