Merge pull request #243 from tdixon97/main

TrackID in Germanium output scheme
legend-exp · Jan 29, 2025 · 115a62b · 115a62b
2 parents 1981b4b + e13e1e8
commit 115a62b
Show file tree

Hide file tree

Showing 9 changed files with 244 additions and 7 deletions.
diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md
diff --git a/include/RMGGermaniumDetector.hh b/include/RMGGermaniumDetector.hh
@@ -25,6 +25,7 @@
 #include "G4VHit.hh"
 #include "G4VSensitiveDetector.hh"
 
+/** @brief Class to store hits in the Germanium detectors, extends @c G4VHit */
 class RMGGermaniumDetectorHit : public G4VHit {
 
   public:
@@ -47,10 +48,13 @@ class RMGGermaniumDetectorHit : public G4VHit {
 
     int detector_uid = -1;
     int particle_type = -1;
+
     double energy_deposition = -1;
     double distance_to_surface = -1;
     G4ThreeVector global_position;
     double global_time = -1;
+    int track_id = -1;
+    int parent_track_id = -1;
 };
 
 using RMGGermaniumDetectorHitsCollection = G4THitsCollection<RMGGermaniumDetectorHit>;

diff --git a/include/RMGGermaniumOutputScheme.hh b/include/RMGGermaniumOutputScheme.hh
@@ -26,20 +26,61 @@
 #include "RMGVOutputScheme.hh"
 
 class G4Event;
+/** @brief Output scheme for Germanium detectors.
+ *
+ *  @details This output scheme records the hits in the Germanium detectors.
+ *  The properties of each @c RMGGermaniumDetectorHit are recorded:
+ *  - event index,
+ *  - particle type,
+ *  - time,
+ *  - position,
+ *  - energy deposition,
+ *  - distance to detector surface.
+ *
+ * Optionally, the track ID and parent track ID can be stored as well as the
+ * detector ID for the single table output mode.
+ *
+ * In addition, this class can be used for stacking tracks associated with optical
+ * photons when no energy was deposited in Germanium.
+ */
 class RMGGermaniumOutputScheme : public RMGVOutputScheme {
 
   public:
 
     RMGGermaniumOutputScheme();
 
+    /** @brief Sets the names of the output columns, invoked in @c RMGRunAction::SetupAnalysisManager */
     void AssignOutputNames(G4AnalysisManager* ana_man) override;
-    void StoreEvent(const G4Event*) override;
-    bool ShouldDiscardEvent(const G4Event*) override;
+
+    /** @brief Store the information from the event, invoked in @c RMGEventAction::EndOfEventAction
+     */
+    void StoreEvent(const G4Event* event) override;
+
+    /** @brief Decide whether to store the event, invoked in @c RMGEventAction::EndOfEventAction
+     *  @details @c true if the event should be discarded, else @c false .
+     *  The event is discarded if there is no hit in Germanium or the energy range
+     *  condition is not met.
+     */
+    bool ShouldDiscardEvent(const G4Event* event) override;
+
+    /** @brief Wraps @c G4UserStackingAction::StackingActionNewStage
+     *  @details discard all waiting events, if @c ShouldDiscardEvent() is true.
+     */
     std::optional<bool> StackingActionNewStage(int) override;
+
+    /** @brief Wraps @c G4UserStackingAction::StackingActionClassify
+     *  @details This is used to classify all optical photon tracks as @c fWaiting if
+     * @c fDiscardPhotonsIfNoGermaniumEdep is true.
+     */
     std::optional<G4ClassificationOfNewTrack> StackingActionClassify(const G4Track*, int) override;
 
+    /** @brief Set a lower cut on the energy deposited in the event to store it. */
     inline void SetEdepCutLow(double threshold) { fEdepCutLow = threshold; }
+
+    /** @brief Set a lower cut on the energy deposited in the event to store it. */
     inline void SetEdepCutHigh(double threshold) { fEdepCutHigh = threshold; }
+
+    /** @brief Add a detector uid to the list of detectors to apply the energy cut for. */
     inline void AddEdepCutDetector(int det_uid) { fEdepCutDetectors.insert(det_uid); }
 
   protected:
@@ -61,6 +102,8 @@ class RMGGermaniumOutputScheme : public RMGVOutputScheme {
 
     bool fStoreSinglePrecisionEnergy = false;
     bool fStoreSinglePrecisionPosition = false;
+
+    bool fStoreTrackID = false;
 };
 
 #endif

diff --git a/src/RMGGermaniumDetector.cc b/src/RMGGermaniumDetector.cc
@@ -126,10 +126,11 @@ bool RMGGermaniumDetector::ProcessHits(G4Step* step, G4TouchableHistory* /*histo
   hit->energy_deposition = step->GetTotalEnergyDeposit();
   hit->global_position = position;
   hit->global_time = prestep->GetGlobalTime();
+  hit->track_id = step->GetTrack()->GetTrackID();
+  hit->parent_track_id = step->GetTrack()->GetParentID();
 
   // Get distance to surface.
   // Check distance to surfaces of Mother volume
-  //
 
   // First transform coordinates into local system
   G4AffineTransform tf(pv->GetRotation(), pv->GetTranslation());

diff --git a/src/RMGGermaniumOutputScheme.cc b/src/RMGGermaniumOutputScheme.cc
@@ -33,7 +33,7 @@ namespace u = CLHEP;
 
 RMGGermaniumOutputScheme::RMGGermaniumOutputScheme() { this->DefineCommands(); }
 
-// invoked in RMGRunAction::SetupAnalysisManager()
+
 void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
 
   auto rmg_man = RMGManager::Instance();
@@ -61,9 +61,17 @@ void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
         rmg_man->RegisterNtuple(det.second.uid, ana_man->CreateNtuple(ntuple_name, "Event data"));
     registered_ntuples.emplace(ntuple_name, id);
 
+    // store the indices
     ana_man->CreateNtupleIColumn(id, "evtid");
     if (!fNtuplePerDetector) { ana_man->CreateNtupleIColumn(id, "det_uid"); }
     ana_man->CreateNtupleIColumn(id, "particle");
+
+    // also store track IDs if instructed
+    if (fStoreTrackID) {
+      ana_man->CreateNtupleIColumn(id, "trackid");
+      ana_man->CreateNtupleIColumn(id, "parent_trackid");
+    }
+    // store the floating points values
     CreateNtupleFOrDColumn(ana_man, id, "edep_in_keV", fStoreSinglePrecisionEnergy);
     ana_man->CreateNtupleDColumn(id, "time_in_ns");
     CreateNtupleFOrDColumn(ana_man, id, "xloc_in_m", fStoreSinglePrecisionPosition);
@@ -95,7 +103,6 @@ RMGGermaniumDetectorHitsCollection* RMGGermaniumOutputScheme::GetHitColl(const G
   return hit_coll;
 }
 
-// invoked in RMGEventAction::EndOfEventAction()
 bool RMGGermaniumOutputScheme::ShouldDiscardEvent(const G4Event* event) {
   // exit fast if no threshold is configured.
   if ((fEdepCutLow < 0 && fEdepCutHigh < 0) || fEdepCutDetectors.empty()) return false;
@@ -122,7 +129,6 @@ bool RMGGermaniumOutputScheme::ShouldDiscardEvent(const G4Event* event) {
   return false;
 }
 
-// invoked in RMGEventAction::EndOfEventAction()
 void RMGGermaniumOutputScheme::StoreEvent(const G4Event* event) {
   auto hit_coll = GetHitColl(event);
   if (!hit_coll) return;
@@ -146,11 +152,19 @@ void RMGGermaniumOutputScheme::StoreEvent(const G4Event* event) {
       auto ntupleid = rmg_man->GetNtupleID(hit->detector_uid);
 
       int col_id = 0;
+      // store the indices
       ana_man->FillNtupleIColumn(ntupleid, col_id++, event->GetEventID());
       if (!fNtuplePerDetector) {
         ana_man->FillNtupleIColumn(ntupleid, col_id++, hit->detector_uid);
       }
       ana_man->FillNtupleIColumn(ntupleid, col_id++, hit->particle_type);
+
+      // store track IDs if instructed
+      if (fStoreTrackID) {
+        ana_man->FillNtupleIColumn(ntupleid, col_id++, hit->track_id);
+        ana_man->FillNtupleIColumn(ntupleid, col_id++, hit->parent_track_id);
+      }
+
       FillNtupleFOrDColumn(ana_man, ntupleid, col_id++, hit->energy_deposition / u::keV,
           fStoreSinglePrecisionEnergy);
       ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_time / u::ns);
@@ -229,6 +243,10 @@ void RMGGermaniumOutputScheme::DefineCommands() {
   fMessenger->DeclareProperty("StoreSinglePrecisionEnergy", fStoreSinglePrecisionEnergy)
       .SetGuidance("Use float32 (instead of float64) for energy output.")
       .SetStates(G4State_Idle);
+
+  fMessenger->DeclareProperty("StoreTrackID", fStoreTrackID)
+      .SetGuidance("Store Track IDs for hits in the output file.")
+      .SetStates(G4State_Idle);
 }
 
 // vim: tabstop=2 shiftwidth=2 expandtab
diff --git a/tests/output/CMakeLists.txt b/tests/output/CMakeLists.txt
@@ -25,7 +25,7 @@ foreach(_file ${_aux})
 endforeach()
 
 # MAIN OUTPUT TESTS
-set(_macros ntuple-per-det.mac ntuple-per-det-vol.mac ntuple-flat.mac)
+set(_macros ntuple-per-det.mac ntuple-per-det-vol.mac ntuple-flat.mac ntuple-flat-with-tracks.mac)
 
 foreach(_mac ${_macros})
   add_test(NAME output/hdf5-${_mac} COMMAND ./run-test-hdf5.sh $<TARGET_FILE:remage-cli-cpp>

diff --git a/tests/output/macros/ntuple-flat-with-tracks.hdf5ls b/tests/output/macros/ntuple-flat-with-tracks.hdf5ls
@@ -0,0 +1,108 @@
+default_histograms
+header
+stp
+stp/germanium
+stp/germanium/columns
+stp/germanium/det_uid
+stp/germanium/det_uid/entries
+stp/germanium/det_uid/pages
+stp/germanium/dist_to_surf_in_m
+stp/germanium/dist_to_surf_in_m/entries
+stp/germanium/dist_to_surf_in_m/pages
+stp/germanium/edep_in_keV
+stp/germanium/edep_in_keV/entries
+stp/germanium/edep_in_keV/pages
+stp/germanium/entries
+stp/germanium/evtid
+stp/germanium/evtid/entries
+stp/germanium/evtid/pages
+stp/germanium/forms
+stp/germanium/names
+stp/germanium/parent_trackid
+stp/germanium/parent_trackid/entries
+stp/germanium/parent_trackid/pages
+stp/germanium/particle
+stp/germanium/particle/entries
+stp/germanium/particle/pages
+stp/germanium/time_in_ns
+stp/germanium/time_in_ns/entries
+stp/germanium/time_in_ns/pages
+stp/germanium/trackid
+stp/germanium/trackid/entries
+stp/germanium/trackid/pages
+stp/germanium/xloc_in_m
+stp/germanium/xloc_in_m/entries
+stp/germanium/xloc_in_m/pages
+stp/germanium/yloc_in_m
+stp/germanium/yloc_in_m/entries
+stp/germanium/yloc_in_m/pages
+stp/germanium/zloc_in_m
+stp/germanium/zloc_in_m/entries
+stp/germanium/zloc_in_m/pages
+stp/scintillator
+stp/scintillator/columns
+stp/scintillator/det_uid
+stp/scintillator/det_uid/entries
+stp/scintillator/det_uid/pages
+stp/scintillator/edep_in_keV
+stp/scintillator/edep_in_keV/entries
+stp/scintillator/edep_in_keV/pages
+stp/scintillator/entries
+stp/scintillator/evtid
+stp/scintillator/evtid/entries
+stp/scintillator/evtid/pages
+stp/scintillator/forms
+stp/scintillator/names
+stp/scintillator/particle
+stp/scintillator/particle/entries
+stp/scintillator/particle/pages
+stp/scintillator/time_in_ns
+stp/scintillator/time_in_ns/entries
+stp/scintillator/time_in_ns/pages
+stp/scintillator/v_post_in_m\ns
+stp/scintillator/v_post_in_m\ns/entries
+stp/scintillator/v_post_in_m\ns/pages
+stp/scintillator/v_pre_in_m\ns
+stp/scintillator/v_pre_in_m\ns/entries
+stp/scintillator/v_pre_in_m\ns/pages
+stp/scintillator/xloc_post_in_m
+stp/scintillator/xloc_post_in_m/entries
+stp/scintillator/xloc_post_in_m/pages
+stp/scintillator/xloc_pre_in_m
+stp/scintillator/xloc_pre_in_m/entries
+stp/scintillator/xloc_pre_in_m/pages
+stp/scintillator/yloc_post_in_m
+stp/scintillator/yloc_post_in_m/entries
+stp/scintillator/yloc_post_in_m/pages
+stp/scintillator/yloc_pre_in_m
+stp/scintillator/yloc_pre_in_m/entries
+stp/scintillator/yloc_pre_in_m/pages
+stp/scintillator/zloc_post_in_m
+stp/scintillator/zloc_post_in_m/entries
+stp/scintillator/zloc_post_in_m/pages
+stp/scintillator/zloc_pre_in_m
+stp/scintillator/zloc_pre_in_m/entries
+stp/scintillator/zloc_pre_in_m/pages
+stp/vertices
+stp/vertices/columns
+stp/vertices/entries
+stp/vertices/evtid
+stp/vertices/evtid/entries
+stp/vertices/evtid/pages
+stp/vertices/forms
+stp/vertices/n_part
+stp/vertices/n_part/entries
+stp/vertices/n_part/pages
+stp/vertices/names
+stp/vertices/time_in_ns
+stp/vertices/time_in_ns/entries
+stp/vertices/time_in_ns/pages
+stp/vertices/xloc_in_m
+stp/vertices/xloc_in_m/entries
+stp/vertices/xloc_in_m/pages
+stp/vertices/yloc_in_m
+stp/vertices/yloc_in_m/entries
+stp/vertices/yloc_in_m/pages
+stp/vertices/zloc_in_m
+stp/vertices/zloc_in_m/entries
+stp/vertices/zloc_in_m/pages
diff --git a/tests/output/macros/ntuple-flat-with-tracks.hdf5ls-lh5 b/tests/output/macros/ntuple-flat-with-tracks.hdf5ls-lh5
@@ -0,0 +1,34 @@
+stp : struct{germanium,scintillator,vertices}
+stp/germanium : table{evtid,det_uid,particle,trackid,parent_trackid,edep,time,xloc,yloc,zloc,dist_to_surf}
+stp/germanium/det_uid : array<1>{real}
+stp/germanium/dist_to_surf : array<1>{real}
+stp/germanium/edep : array<1>{real}
+stp/germanium/evtid : array<1>{real}
+stp/germanium/parent_trackid : array<1>{real}
+stp/germanium/particle : array<1>{real}
+stp/germanium/time : array<1>{real}
+stp/germanium/trackid : array<1>{real}
+stp/germanium/xloc : array<1>{real}
+stp/germanium/yloc : array<1>{real}
+stp/germanium/zloc : array<1>{real}
+stp/scintillator : table{evtid,det_uid,particle,edep,time,xloc_pre,yloc_pre,zloc_pre,xloc_post,yloc_post,zloc_post,v_pre,v_post}
+stp/scintillator/det_uid : array<1>{real}
+stp/scintillator/edep : array<1>{real}
+stp/scintillator/evtid : array<1>{real}
+stp/scintillator/particle : array<1>{real}
+stp/scintillator/time : array<1>{real}
+stp/scintillator/v_post : array<1>{real}
+stp/scintillator/v_pre : array<1>{real}
+stp/scintillator/xloc_post : array<1>{real}
+stp/scintillator/xloc_pre : array<1>{real}
+stp/scintillator/yloc_post : array<1>{real}
+stp/scintillator/yloc_pre : array<1>{real}
+stp/scintillator/zloc_post : array<1>{real}
+stp/scintillator/zloc_pre : array<1>{real}
+stp/vertices : table{evtid,time,xloc,yloc,zloc,n_part}
+stp/vertices/evtid : array<1>{real}
+stp/vertices/n_part : array<1>{real}
+stp/vertices/time : array<1>{real}
+stp/vertices/xloc : array<1>{real}
+stp/vertices/yloc : array<1>{real}
+stp/vertices/zloc : array<1>{real}
diff --git a/tests/output/macros/ntuple-flat-with-tracks.mac b/tests/output/macros/ntuple-flat-with-tracks.mac
@@ -0,0 +1,20 @@
+/RMG/Output/NtuplePerDetector false
+
+/RMG/Geometry/RegisterDetector Scintillator scint1 1
+/RMG/Geometry/RegisterDetector Scintillator scint2 2
+/RMG/Geometry/RegisterDetector Germanium det1 11
+/RMG/Geometry/RegisterDetector Germanium det2 12
+
+/control/execute macros/_init.mac
+/RMG/Output/Germanium/StoreTrackID true
+
+
+/RMG/Generator/Confine UnConfined
+
+/RMG/Generator/Select GPS
+/gps/position 0 0 0
+/gps/ang/type iso
+/gps/particle gamma
+/gps/energy 200 keV
+
+/run/beamOn 50