From d91a08ad718703c835a2310ffb906bb204527e13 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Tue, 22 Oct 2024 16:09:36 +0200 Subject: [PATCH 1/2] Add option to reset global time after initial decay --- docs/rmg-commands.md | 10 ++++++++ include/RMGTrackingAction.hh | 9 +++++++ src/RMGTrackingAction.cc | 48 ++++++++++++++++++++++++++++++++++-- 3 files changed, 65 insertions(+), 2 deletions(-) diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index 268712f..896945c 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -502,6 +502,7 @@ Commands for controlling physics processes ### Commands * `DaughterNucleusMaxLifetime` – Determines which unstable daughter nuclei will be killed, if they are at rest, depending on their lifetime. +* `ResetInitialDecayTime` – If the initial step is a radioactive decay, reset the global time of all its secondary tracks to 0. ### `/RMG/Processes/Stepping/DaughterNucleusMaxLifetime` @@ -521,6 +522,15 @@ Set to -1 to disable this feature. * **Default value** – `us` * **Candidates** – `s ms us ns ps min h d y second millisecond microsecond nanosecond picosecond minute hour day year` +### `/RMG/Processes/Stepping/ResetInitialDecayTime` + +If the initial step is a radioactive decay, reset the global time of all its secondary tracks to 0. + +* **Parameter** – `value` + * **Parameter type** – `b` + * **Omittable** – `False` + * **Default value** – `false` + ## `/RMG/Geometry/` Commands for controlling geometry definitions diff --git a/include/RMGTrackingAction.hh b/include/RMGTrackingAction.hh index f4b6b9a..9f587f4 100644 --- a/include/RMGTrackingAction.hh +++ b/include/RMGTrackingAction.hh @@ -16,6 +16,9 @@ #ifndef _RMG_TRACKING_ACTION_HH_ #define _RMG_TRACKING_ACTION_HH_ +#include + +#include "G4GenericMessenger.hh" #include "G4UserTrackingAction.hh" class RMGRunAction; @@ -40,6 +43,12 @@ class RMGTrackingAction : public G4UserTrackingAction { private: RMGRunAction* fRunAction = nullptr; + bool fResetInitialDecayTime = false; + + void ResetInitialDecayTime(const G4Track*); + + std::unique_ptr fMessenger; + void DefineCommands(); }; #endif diff --git a/src/RMGTrackingAction.cc b/src/RMGTrackingAction.cc index cab6ff7..950ab73 100644 --- a/src/RMGTrackingAction.cc +++ b/src/RMGTrackingAction.cc @@ -15,17 +15,61 @@ #include "RMGTrackingAction.hh" +#include "G4Event.hh" +#include "G4EventManager.hh" +#include "G4RadioactiveDecay.hh" #include "G4Track.hh" +#include "G4TrackingManager.hh" +#include "RMGLog.hh" #include "RMGRunAction.hh" -RMGTrackingAction::RMGTrackingAction(RMGRunAction* run_action) : fRunAction(run_action) {} +RMGTrackingAction::RMGTrackingAction(RMGRunAction* run_action) : fRunAction(run_action) { + + this->DefineCommands(); +} void RMGTrackingAction::PreUserTrackingAction(const G4Track* aTrack) { for (auto& el : fRunAction->GetAllOutputDataFields()) { el->TrackingActionPre(aTrack); } } -void RMGTrackingAction::PostUserTrackingAction(const G4Track* /*aTrack*/) {} +void RMGTrackingAction::PostUserTrackingAction(const G4Track* aTrack) { + + if (fResetInitialDecayTime) ResetInitialDecayTime(aTrack); +} + +void RMGTrackingAction::ResetInitialDecayTime(const G4Track* aTrack) { + + // only nuclei in the first step are eligible to be reset. + if (aTrack->GetTrackID() != 1 || aTrack->GetParentID() != 0) return; + if (aTrack->GetDefinition()->GetParticleType() != "nucleus") return; + + // only reset the time if the last process is a radioactive decay. + auto creator_process = aTrack->GetStep()->GetPostStepPoint()->GetProcessDefinedStep(); + if (!dynamic_cast(creator_process)) return; + + const auto secondaries = fpTrackingManager->GimmeSecondaries(); + auto secondaries_in_current_step = aTrack->GetStep()->GetNumberOfSecondariesInCurrentStep(); + // if we have more secondaries than from the final decay, the earlier ones will have larger times. + if (secondaries_in_current_step != secondaries->size()) { + RMGLog::Out(RMGLog::warning, "inconsistent (non-monotonous) timing in event ", + G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID()); + } + + for (auto sec : *secondaries) { sec->SetGlobalTime(0.); } +} + +void RMGTrackingAction::DefineCommands() { + + fMessenger = std::make_unique(this, "/RMG/Processes/Stepping/", + "Commands for controlling physics processes"); + + fMessenger->DeclareProperty("ResetInitialDecayTime", fResetInitialDecayTime) + .SetGuidance("If the initial step is a radioactive decay, reset the global time of all its " + "secondary tracks to 0.") + .SetDefaultValue("false") + .SetStates(G4State_PreInit); +} // vim: tabstop=2 shiftwidth=2 expandtab From 3e4a9025754ef33ab49ff12fd73ceffa218d1480 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Wed, 13 Nov 2024 12:40:42 +0100 Subject: [PATCH 2/2] add warning for large global times --- include/RMGTrackingAction.hh | 3 ++- src/RMGTrackingAction.cc | 33 ++++++++++++++++++++++++++++----- 2 files changed, 30 insertions(+), 6 deletions(-) diff --git a/include/RMGTrackingAction.hh b/include/RMGTrackingAction.hh index 9f587f4..894b561 100644 --- a/include/RMGTrackingAction.hh +++ b/include/RMGTrackingAction.hh @@ -44,8 +44,9 @@ class RMGTrackingAction : public G4UserTrackingAction { RMGRunAction* fRunAction = nullptr; bool fResetInitialDecayTime = false; + bool fHadLongTimeWarning = false; - void ResetInitialDecayTime(const G4Track*); + bool ResetInitialDecayTime(const G4Track*); std::unique_ptr fMessenger; void DefineCommands(); diff --git a/src/RMGTrackingAction.cc b/src/RMGTrackingAction.cc index 950ab73..230e25a 100644 --- a/src/RMGTrackingAction.cc +++ b/src/RMGTrackingAction.cc @@ -24,6 +24,12 @@ #include "RMGLog.hh" #include "RMGRunAction.hh" +namespace { + template constexpr double const_pow(T base, T exp) { + return exp == 0 ? 1 : base * const_pow(base, exp - 1); + } +} // namespace + RMGTrackingAction::RMGTrackingAction(RMGRunAction* run_action) : fRunAction(run_action) { this->DefineCommands(); @@ -36,18 +42,33 @@ void RMGTrackingAction::PreUserTrackingAction(const G4Track* aTrack) { void RMGTrackingAction::PostUserTrackingAction(const G4Track* aTrack) { - if (fResetInitialDecayTime) ResetInitialDecayTime(aTrack); + bool check_global_time = true; + if (fResetInitialDecayTime) { check_global_time = !ResetInitialDecayTime(aTrack); } + + // this is just a good "guess" that might not hold true in all cases, i.e. some us values + // might still not be unique below this. + constexpr double max_representable_time_with_us_prec = + const_pow(2, std::numeric_limits::digits) * CLHEP::us; + + if (check_global_time && !fHadLongTimeWarning) { + if (aTrack->GetGlobalTime() > max_representable_time_with_us_prec) { + RMGLog::Out(RMGLog::warning, "encountered long global time (> ", + max_representable_time_with_us_prec / CLHEP::year, + " yr). Global time precision might be worse than 1 us."); + fHadLongTimeWarning = true; + } + } } -void RMGTrackingAction::ResetInitialDecayTime(const G4Track* aTrack) { +bool RMGTrackingAction::ResetInitialDecayTime(const G4Track* aTrack) { // only nuclei in the first step are eligible to be reset. - if (aTrack->GetTrackID() != 1 || aTrack->GetParentID() != 0) return; - if (aTrack->GetDefinition()->GetParticleType() != "nucleus") return; + if (aTrack->GetTrackID() != 1 || aTrack->GetParentID() != 0) return false; + if (aTrack->GetDefinition()->GetParticleType() != "nucleus") return false; // only reset the time if the last process is a radioactive decay. auto creator_process = aTrack->GetStep()->GetPostStepPoint()->GetProcessDefinedStep(); - if (!dynamic_cast(creator_process)) return; + if (!dynamic_cast(creator_process)) return false; const auto secondaries = fpTrackingManager->GimmeSecondaries(); auto secondaries_in_current_step = aTrack->GetStep()->GetNumberOfSecondariesInCurrentStep(); @@ -58,6 +79,8 @@ void RMGTrackingAction::ResetInitialDecayTime(const G4Track* aTrack) { } for (auto sec : *secondaries) { sec->SetGlobalTime(0.); } + + return true; } void RMGTrackingAction::DefineCommands() {