modular stacking action with RMGVOutputScheme

include/RMGGermaniumOutputScheme.hh

@@ -16,6 +16,7 @@
 #ifndef _RMG_GERMANIUM_OUTPUT_SCHEME_HH_
 #define _RMG_GERMANIUM_OUTPUT_SCHEME_HH_
 
+#include <optional>
 #include <set>
 
 #include "G4AnalysisManager.hh"
@@ -34,15 +35,13 @@ class RMGGermaniumOutputScheme : public RMGVOutputScheme {
     void AssignOutputNames(G4AnalysisManager* ana_man) override;
     void StoreEvent(const G4Event*) override;
     bool ShouldDiscardEvent(const G4Event*) override;
+    std::optional<bool> StackingActionNewStage(int) override;
+    std::optional<G4ClassificationOfNewTrack> StackingActionClassify(const G4Track*, int) override;
 
     inline void SetEdepCutLow(double threshold) { fEdepCutLow = threshold; }
     inline void SetEdepCutHigh(double threshold) { fEdepCutHigh = threshold; }
     inline void AddEdepCutDetector(int det_uid) { fEdepCutDetectors.insert(det_uid); }
 
-    [[nodiscard]] inline bool GetDiscardPhotonsIfNoGermaniumEdep() const {
-      return fDiscardPhotonsIfNoGermaniumEdep;
-    }
-
   private:
 
     RMGGermaniumDetectorHitsCollection* GetHitColl(const G4Event*);

include/RMGRunAction.hh

@@ -52,9 +52,10 @@ class RMGRunAction : public G4UserRunAction {
 
     inline std::shared_ptr<RMGVOutputScheme> GetOutputDataFields(RMGHardware::DetectorType d_type) {
       auto it = fOutputDataFields.find(d_type);
-      if (it != fOutputDataFields.end()) return (*it).second;
+      if (it != fOutputDataFields.end()) return it->second;
       return nullptr;
     }
+    [[nodiscard]] inline const auto& GetAllOutputDataFields() { return fOutputDataFields; }
     inline void ClearOutputDataFields() {
       for (auto& el : fOutputDataFields) el.second->ClearBeforeEvent();
     }

include/RMGStackingAction.hh

@@ -38,6 +38,7 @@ class RMGStackingAction : public G4UserStackingAction {
 
   private:
 
+    int fStage = 0;
     RMGRunAction* fRunAction = nullptr;
 };
 

include/RMGVOutputScheme.hh

@@ -16,9 +16,12 @@
 #ifndef _RMG_V_OUTPUT_SCHEME_HH_
 #define _RMG_V_OUTPUT_SCHEME_HH_
 
+#include <optional>
 #include <string>
 
 #include "G4AnalysisManager.hh"
+#include "G4Track.hh"
+#include "G4UserStackingAction.hh"
 
 #include "fmt/format.h"
 
@@ -36,6 +39,12 @@ class RMGVOutputScheme {
     virtual inline void ClearBeforeEvent() {}
     virtual inline bool ShouldDiscardEvent(const G4Event*) { return false; }
     virtual inline void StoreEvent(const G4Event*) {}
+    virtual inline void StackingAction() {}
+    virtual inline std::optional<G4ClassificationOfNewTrack> StackingActionClassify(const G4Track*,
+        const int) {
+      return std::nullopt;
+    }
+    virtual inline std::optional<bool> StackingActionNewStage(const int) { return std::nullopt; }
 
     inline std::string GetNtupleName(int det_uid) { return fmt::format("det{:03}", det_uid); }
 };

src/RMGGermaniumOutputScheme.cc

@@ -20,6 +20,7 @@
 #include "G4AnalysisManager.hh"
 #include "G4Event.hh"
 #include "G4HCtable.hh"
+#include "G4OpticalPhoton.hh"
 #include "G4SDManager.hh"
 
 #include "RMGGermaniumDetector.hh"
@@ -131,20 +132,38 @@ void RMGGermaniumOutputScheme::StoreEvent(const G4Event* event) {
       const auto ana_man = G4AnalysisManager::Instance();
       auto ntupleid = rmg_man->GetNtupleID(hit->detector_uid);
 
-      ana_man->FillNtupleIColumn(ntupleid, 0, event->GetEventID());
-      ana_man->FillNtupleIColumn(ntupleid, 1, hit->particle_type);
-      ana_man->FillNtupleDColumn(ntupleid, 2, hit->energy_deposition / u::keV);
-      ana_man->FillNtupleDColumn(ntupleid, 3, hit->global_time / u::ns);
-      ana_man->FillNtupleDColumn(ntupleid, 4, hit->global_position.getX() / u::m);
-      ana_man->FillNtupleDColumn(ntupleid, 5, hit->global_position.getY() / u::m);
-      ana_man->FillNtupleDColumn(ntupleid, 6, hit->global_position.getZ() / u::m);
+      int col_id = 0;
+      ana_man->FillNtupleIColumn(ntupleid, col_id++, event->GetEventID());
+      ana_man->FillNtupleIColumn(ntupleid, col_id++, hit->particle_type);
+      ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->energy_deposition / u::keV);
+      ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_time / u::ns);
+      ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position.getX() / u::m);
+      ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position.getY() / u::m);
+      ana_man->FillNtupleDColumn(ntupleid, col_id++, hit->global_position.getZ() / u::m);
 
       // NOTE: must be called here for hit-oriented output
       ana_man->AddNtupleRow(ntupleid);
     }
   }
 }
 
+std::optional<G4ClassificationOfNewTrack> RMGGermaniumOutputScheme::StackingActionClassify(const G4Track* aTrack,
+    int stage) {
+  if (stage != 0) return std::nullopt;
+  // defer tracking of optical photons.
+  if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) return fWaiting;
+  return fUrgent;
+}
+
+std::optional<bool> RMGGermaniumOutputScheme::StackingActionNewStage(const int stage) {
+  if (!fDiscardPhotonsIfNoGermaniumEdep || stage != 0) return true;
+
+  auto run_man = RMGManager::Instance()->GetG4RunManager();
+  const auto event = run_man->GetCurrentEvent();
+  // discard all waiting events, as there was no energy deposition in Germanium.
+  return !ShouldDiscardEvent(event);
+}
+
 void RMGGermaniumOutputScheme::DefineCommands() {
 
   fMessenger = std::make_unique<G4GenericMessenger>(this, "/RMG/Output/Germanium/",

src/RMGStackingAction.cc

@@ -15,42 +15,59 @@
 
 #include "RMGStackingAction.hh"
 
+#include <optional>
+
 #include "G4EventManager.hh"
-#include "G4OpticalPhoton.hh"
 #include "G4Track.hh"
 
 #include "RMGEventAction.hh"
-#include "RMGGermaniumOutputScheme.hh"
 #include "RMGHardware.hh"
+#include "RMGLog.hh"
 #include "RMGManager.hh"
-#include "RMGOpticalOutputScheme.hh"
 #include "RMGRunAction.hh"
 
 RMGStackingAction::RMGStackingAction(RMGRunAction* runaction) : fRunAction(runaction) {}
 
 G4ClassificationOfNewTrack RMGStackingAction::ClassifyNewTrack(const G4Track* aTrack) {
-  // defer tracking of optical photons.
-  if (aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) return fWaiting;
+  std::optional<G4ClassificationOfNewTrack> new_status = std::nullopt;
+  for (auto& el : fRunAction->GetAllOutputDataFields()) {
+    auto request_status = el.second->StackingActionClassify(aTrack, fStage);
+    if (!request_status.has_value()) continue; // this output scheme does not care.
+
+    if (!new_status.has_value() || new_status.value() == request_status.value()) {
+      new_status = request_status;
+    } else {
+      RMGLog::Out(RMGLog::error,
+          "Conflicting requests for new track classification in RMGStackingAction.");
+    }
+  }
+
+  if (new_status.has_value()) return new_status.value();
   return fUrgent;
 }
 
 void RMGStackingAction::NewStage() {
-  auto ge_output = fRunAction->GetOutputDataFields(RMGHardware::DetectorType::kGermanium);
-  if (ge_output == nullptr) return;
+  // we can have only one result from all output schemes, if we have
+  std::optional<bool> should_do_stage = std::nullopt;
+  for (auto& el : fRunAction->GetAllOutputDataFields()) {
+    auto request_stage = el.second->StackingActionNewStage(fStage);
+    if (!request_stage.has_value()) continue; // this output scheme does not care.
 
-  bool discard_photons =
-      (dynamic_cast<RMGGermaniumOutputScheme&>(*ge_output)).GetDiscardPhotonsIfNoGermaniumEdep();
-  if (!discard_photons) return;
+    if (!should_do_stage.has_value() || should_do_stage.value() == request_stage.value()) {
+      should_do_stage = request_stage;
+    } else {
+      RMGLog::Out(RMGLog::error,
+          "Conflicting requests for new stage termination in RMGStackingAction.");
+    }
+  }
 
-  auto run_man = RMGManager::Instance()->GetG4RunManager();
-  const auto event = run_man->GetCurrentEvent();
-  if (ge_output->ShouldDiscardEvent(event)) {
-    // discard all waiting events, as there was no energy deposition in Germanium.
+  if (should_do_stage.has_value() && !should_do_stage.value()) {
     auto stack_man = G4EventManager::GetEventManager()->GetStackManager();
     stack_man->clear();
   }
+  fStage++;
 }
 
-void RMGStackingAction::PrepareNewEvent() {}
+void RMGStackingAction::PrepareNewEvent() { fStage = 0; }
 
 // vim: tabstop=2 shiftwidth=2 expandtab