From da664fbd0b061253b582d086a47447b168998244 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Mon, 6 May 2024 14:08:50 +0200 Subject: [PATCH] modular stacking action with RMGVOutputScheme --- include/RMGGermaniumOutputScheme.hh | 7 ++--- include/RMGRunAction.hh | 3 +- include/RMGStackingAction.hh | 1 + include/RMGVOutputScheme.hh | 9 ++++++ src/RMGGermaniumOutputScheme.cc | 33 +++++++++++++++----- src/RMGStackingAction.cc | 47 ++++++++++++++++++++--------- 6 files changed, 73 insertions(+), 27 deletions(-) diff --git a/include/RMGGermaniumOutputScheme.hh b/include/RMGGermaniumOutputScheme.hh index ed60b19..12eacf1 100644 --- a/include/RMGGermaniumOutputScheme.hh +++ b/include/RMGGermaniumOutputScheme.hh @@ -16,6 +16,7 @@ #ifndef _RMG_GERMANIUM_OUTPUT_SCHEME_HH_ #define _RMG_GERMANIUM_OUTPUT_SCHEME_HH_ +#include #include #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 StackingActionNewStage(int) override; + std::optional 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*); diff --git a/include/RMGRunAction.hh b/include/RMGRunAction.hh index 0070f59..3505e2e 100644 --- a/include/RMGRunAction.hh +++ b/include/RMGRunAction.hh @@ -52,9 +52,10 @@ class RMGRunAction : public G4UserRunAction { inline std::shared_ptr 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(); } diff --git a/include/RMGStackingAction.hh b/include/RMGStackingAction.hh index 719879f..e55fec7 100644 --- a/include/RMGStackingAction.hh +++ b/include/RMGStackingAction.hh @@ -38,6 +38,7 @@ class RMGStackingAction : public G4UserStackingAction { private: + int fStage = 0; RMGRunAction* fRunAction = nullptr; }; diff --git a/include/RMGVOutputScheme.hh b/include/RMGVOutputScheme.hh index 3693852..28233d9 100644 --- a/include/RMGVOutputScheme.hh +++ b/include/RMGVOutputScheme.hh @@ -16,9 +16,12 @@ #ifndef _RMG_V_OUTPUT_SCHEME_HH_ #define _RMG_V_OUTPUT_SCHEME_HH_ +#include #include #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 StackingActionClassify(const G4Track*, + const int) { + return std::nullopt; + } + virtual inline std::optional StackingActionNewStage(const int) { return std::nullopt; } inline std::string GetNtupleName(int det_uid) { return fmt::format("det{:03}", det_uid); } }; diff --git a/src/RMGGermaniumOutputScheme.cc b/src/RMGGermaniumOutputScheme.cc index 7dc2590..9091add 100644 --- a/src/RMGGermaniumOutputScheme.cc +++ b/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,13 +132,14 @@ 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); @@ -145,6 +147,23 @@ void RMGGermaniumOutputScheme::StoreEvent(const G4Event* event) { } } +std::optional 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 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(this, "/RMG/Output/Germanium/", diff --git a/src/RMGStackingAction.cc b/src/RMGStackingAction.cc index 1b9504f..dbe9a6a 100644 --- a/src/RMGStackingAction.cc +++ b/src/RMGStackingAction.cc @@ -15,42 +15,59 @@ #include "RMGStackingAction.hh" +#include + #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 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 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(*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