diff --git a/CMakeLists.txt b/CMakeLists.txt index 1abf44bef..7677ab14e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,9 @@ CMAKE_POLICY(SET CMP0012 NEW) # if() recognizes numbers and boolean constants. CMAKE_POLICY(SET CMP0025 NEW) # Compiler id for Apple Clang is now AppleClang. CMAKE_POLICY(SET CMP0057 NEW) # Support new if() IN_LIST operator. CMAKE_POLICY(SET CMP0072 NEW) # FindOpenGL prefers GLVND by default when available. +#========================================================= +#set(CMAKE_CXX_COMPILER "/usr/bin/x86_64-linux-gnu-g++-13") +#set(CMAKE_EXPORT_COMPILE_COMMANDS ON) #========================================================= PROJECT(Gate) @@ -20,6 +23,37 @@ IF(BUILD_TESTING) ENABLE_TESTING() SET(BUILDNAME "${BUILDNAME}" CACHE STRING "Name of build on the dashboard") MARK_AS_ADVANCED(BUILDNAME) + #Stupid solution cause it links to everything but for now on + ADD_EXECUTABLE(test_dummy tests/test_dummy.cpp) + TARGET_LINK_LIBRARIES(test_dummy GateLib) + + ADD_TEST(NAME test_dummy COMMAND test_dummy) + + ADD_EXECUTABLE(test_GatePositroniumDecayModel tests/test_GatePositroniumDecayModel.cpp) + TARGET_LINK_LIBRARIES(test_GatePositroniumDecayModel GateLib) + + ADD_TEST(NAME test_GatePositroniumDecayModel COMMAND test_GatePositroniumDecayModel) + + ADD_EXECUTABLE(test_GatePositronium tests/test_GatePositronium.cpp) + TARGET_LINK_LIBRARIES(test_GatePositronium GateLib) + + ADD_TEST(NAME test_GatePositronium COMMAND test_GatePositronium) + + ADD_EXECUTABLE(test_GatePositroniumHelper tests/test_GatePositroniumHelper.cpp) + TARGET_LINK_LIBRARIES(test_GatePositroniumHelper GateLib) + + ADD_TEST(NAME test_GatePositroniumHelper COMMAND test_GatePositroniumHelper) + + ADD_EXECUTABLE(test_GatePositroniumDecyParamsGenerator tests/test_GatePositroniumDecayParamsGenerator.cpp) + TARGET_LINK_LIBRARIES(test_GatePositroniumDecyParamsGenerator GateLib) + + ADD_TEST(NAME test_GatePositroniumDecyParamsGenerator COMMAND test_GatePositroniumDecayParamsGenerator) + + ADD_EXECUTABLE(test_GatePositroniumSourceMessenger tests/test_GatePositroniumSourceMessenger.cpp) + TARGET_LINK_LIBRARIES(test_GatePositroniumSourceMessenger GateLib) + + ADD_TEST(NAME test_GatePositroniumSourceMessenger COMMAND test_GatePositroniumSourceMessenger) + ENDIF(BUILD_TESTING) #========================================================= @@ -226,6 +260,7 @@ FILE(GLOB sources ${PROJECT_SOURCE_DIR}/source/geometry/src/*.cc ${PROJECT_SOURCE_DIR}/source/digits_hits/src/*.cc ${PROJECT_SOURCE_DIR}/source/physics/src/*.cc + ${PROJECT_SOURCE_DIR}/source/physics/src/legacy/*.cc ${PROJECT_SOURCE_DIR}/source/general/src/*.cc ${PROJECT_SOURCE_DIR}/source/externals/clhep/src/CLHEP/Matrix/*.cc ${PROJECT_SOURCE_DIR}/source/externals/clhep/src/CLHEP/RandomObjects/*.cc) @@ -235,6 +270,7 @@ FILE(GLOB headers ${PROJECT_SOURCE_DIR}/source/arf/include/*.hh ${PROJECT_SOURCE_DIR}/source/geometry/include/*.hh ${PROJECT_SOURCE_DIR}/source/physics/include/*.hh + ${PROJECT_SOURCE_DIR}/source/physics/include/legacy/*.hh ${PROJECT_SOURCE_DIR}/source/digits_hits/include/*.hh ${PROJECT_SOURCE_DIR}/source/general/include/*.hh ${PROJECT_SOURCE_DIR}/source/externals/clhep/include/*.hh) diff --git a/source/digits_hits/include/GateHit.hh b/source/digits_hits/include/GateHit.hh index 658cc4f1e..6bbcaacc7 100644 --- a/source/digits_hits/include/GateHit.hh +++ b/source/digits_hits/include/GateHit.hh @@ -78,6 +78,7 @@ public: G4int m_nCrystalCompton; // # of compton processes in the crystal occurred to the photon G4int m_nPhantomRayleigh; // # of Rayleigh processes in the phantom occurred to the photon G4int m_nCrystalRayleigh; // # of Rayleigh processes in the crystal occurred to the photon + G4int m_nInteractions; // # of non-Transportation interactions in phantom + crystal; -1 when not computed G4String m_comptonVolumeName; // name of the volume of the last (if any) compton scattering G4String m_RayleighVolumeName; // name of the volume of the last (if any) Rayleigh scattering G4int m_primaryID; // primary that caused the hit @@ -108,6 +109,7 @@ public: G4int m_sourceType = 0;//0 means 'by default not known'; sourceType says what type of positronium (Ps) is used: pPs, oPs G4int m_decayType = 0;//0 means 'by default not known'; decayType says what type of Ps decay is used: standard (without prompt gamma), deexcitation (with prompt gamma) G4int m_gammaType = 0;//0 means 'by default not known'; gammaType says what type of gamma is emitted: annihilation, prompt, other + G4int m_decayIndex = -1;//decay channel index; -1 means 'by default not known' public: inline void SetEdep(G4double de) { m_edep = de; } @@ -170,6 +172,9 @@ public: inline void SetNCrystalRayleigh(G4int j) { m_nCrystalRayleigh = j; } inline G4int GetNCrystalRayleigh() const { return m_nCrystalRayleigh; } + inline void SetNInteractions(G4int j) { m_nInteractions = j; } + inline G4int GetNInteractions() const { return m_nInteractions; } + inline void SetComptonVolumeName(G4String name) { m_comptonVolumeName = name; } inline G4String GetComptonVolumeName() const { return m_comptonVolumeName; } @@ -256,6 +261,9 @@ public: inline void SetGammaType(G4int value){ m_gammaType = value; } inline G4int GetGammaType() const { return m_gammaType; } + + inline void SetDecayIndex(G4int value){ m_decayIndex = value; } + inline G4int GetDecayIndex() const { return m_decayIndex; } }; diff --git a/source/digits_hits/include/GateMultiPhotonAnalysis.hh b/source/digits_hits/include/GateMultiPhotonAnalysis.hh new file mode 100644 index 000000000..cb8ef7ab1 --- /dev/null +++ b/source/digits_hits/include/GateMultiPhotonAnalysis.hh @@ -0,0 +1,166 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +#ifndef GateMultiPhotonAnalysis_h +#define GateMultiPhotonAnalysis_h + +#include "GateMultiPhotonAnalysisHelpers.hh" +#include "GateVOutputModule.hh" + +#include + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About: Multi-photon analysis output-module interface with event processing callbacks and trajectory-missing policies. + **/ + +class GateVVolume; +class GateMultiPhotonTrajectoryNavigator; +class GateMultiPhotonAnalysisMessenger; + +/** + * @brief Multi-photon variant of event-end hit analysis. + * + * The class mirrors the legacy GateAnalysis role while removing the hardcoded + * two-photon mapping assumption. It aggregates per-photon interaction counters + * and assigns resolved values to crystal hits. + * + * Notes: + * - Iteration 1 supports TrackingMode::kBoth only. + * - Tracker/detector split mode handling is explicitly deferred. + */ +class GateMultiPhotonAnalysis : public GateVOutputModule { + public: + /** @brief Policy for handling missing trajectory containers. */ + enum MissingTrajectoryPolicy { + kStrict, + kResilient + }; + + /** + * @brief Constructs the output module. + * + * Args: + * name: Output module name. + * outputMgr: Output manager owner. + * digiMode: Runtime/offline mode. + */ + GateMultiPhotonAnalysis(const G4String &name, GateOutputMgr *outputMgr, DigiMode digiMode); + + virtual ~GateMultiPhotonAnalysis(); + + /** + * @brief Returns module file name placeholder. + * + * Returns: + * Dummy name since this module does not own a file. + */ + const G4String &GiveNameOfFile(); + + /** @brief Acquisition-begin callback. */ + void RecordBeginOfAcquisition(); + /** @brief Acquisition-end callback. */ + void RecordEndOfAcquisition(); + /** @brief Run-begin callback. */ + void RecordBeginOfRun(const G4Run *); + /** @brief Run-end callback. */ + void RecordEndOfRun(const G4Run *); + /** @brief Event-begin callback. */ + void RecordBeginOfEvent(const G4Event *); + + /** + * @brief Performs event-end multi-photon aggregation and hit assignment. + * + * Args: + * event: Event carrying trajectory and hit collections. + */ + void RecordEndOfEvent(const G4Event *event); + + /** + * @brief Step callback (unused by this module). + * + * Args: + * v: Current volume. + * step: Current step. + */ + void RecordStepWithVolume(const GateVVolume *v, const G4Step *step); + + /** @brief Voxel callback (not used). */ + void RecordVoxels(GateVGeometryVoxelStore *) {} + + /** + * @brief Sets module and navigator verbosity. + * + * Args: + * val: Verbose level. + */ + virtual void SetVerboseLevel(G4int val); + + /** + * @brief Sets policy for missing trajectory container handling. + * + * Args: + * policy: Desired behavior policy. + */ + void SetMissingTrajectoryPolicy(MissingTrajectoryPolicy policy); + + /** + * @brief Parses and applies missing trajectory policy from text. + * + * Args: + * policyName: Policy name (`strict` or `resilient`). + */ + void SetMissingTrajectoryPolicyFromString(const G4String &policyName); + + /** + * @brief Returns current missing trajectory policy. + * + * Returns: + * Active policy value. + */ + MissingTrajectoryPolicy GetMissingTrajectoryPolicy() const; + + /** + * @brief Returns current missing trajectory policy as string. + * + * Returns: + * `strict` or `resilient`. + */ + G4String GetMissingTrajectoryPolicyName() const; + + private: + /** + * @brief Checks whether current tracking mode is supported. + * + * Args: + * tracking_mode_code: Integer value of TrackingMode. + * + * Returns: + * True when event processing should continue for the current mode. + */ + bool IsTrackingModeSupported(int tracking_mode_code) const; + + /** + * @brief Returns policy value for legacy photonID field. + * + * Returns: + * Always 0 in multiphoton mode (explicit policy for iteration 1). + */ + int ResolveLegacyPhotonIDPolicy() const; + + GateMultiPhotonTrajectoryNavigator *m_trajectoryNavigator; + GateMultiPhotonAnalysisMessenger *m_messenger; + G4String m_noFileName; + MissingTrajectoryPolicy m_missingTrajectoryPolicy; + G4int m_missingTrajectoryEventCount; + G4int m_missingTrajectoryWithHitsCount; +}; + +#endif diff --git a/source/digits_hits/include/GateMultiPhotonAnalysisHelpers.hh b/source/digits_hits/include/GateMultiPhotonAnalysisHelpers.hh new file mode 100644 index 000000000..ef10a140a --- /dev/null +++ b/source/digits_hits/include/GateMultiPhotonAnalysisHelpers.hh @@ -0,0 +1,117 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#ifndef GateMultiPhotonAnalysisHelpers_h +#define GateMultiPhotonAnalysisHelpers_h + +#include +#include + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About namespace: Helper functions for multi-photon analysis, including a structure to hold gamma interaction statistics and functions to accumulate statistics based on interaction process names. + **/ +namespace MultiPhotonAnalysisHelpers { + + /** @brief Aggregated interaction counters for a single tracked gamma. */ + struct GammaStatistics { + int phantomCompton = 0; + int phantomRayleigh = 0; + int crystalCompton = 0; + int crystalRayleigh = 0; + int phantomInteractions = 0; + int crystalInteractions = 0; + std::string comptonVolumeName = "NULL"; + std::string rayleighVolumeName = "NULL"; + }; + + /** @brief Interaction process classes used by multi-photon analysis. */ + enum class InteractionProcess { + Unknown, + Compton, + Rayleigh, + Transportation + }; + + /** + * @brief Classifies process name into a compact process category. + * + * Args: + * processName: Geant4 process name. + * + * Returns: + * Recognized interaction category or Unknown. + */ + constexpr InteractionProcess GetInteractionProcess(const std::string_view processName) { + if (processName.size() >= 4) { + if (processName.substr(0, 4) == "Tran") { + return InteractionProcess::Transportation; + } + if (processName.substr(0, 4) == "Comp" || processName.substr(0, 4) == "comp") { + return InteractionProcess::Compton; + } + if (processName.substr(0, 4) == "Rayl" || processName.substr(0, 4) == "rayl") { + return InteractionProcess::Rayleigh; + } + } + return InteractionProcess::Unknown; + } + + /** + * @brief Accumulates crystal-side interaction statistics. + * + * Args: + * processName: Geant4 process name at crystal hit. + * gammaStatistics: Statistics object to update. + */ + inline void AccumulateCrystal(const std::string_view processName, GammaStatistics& gammaStatistics) { + const InteractionProcess process = GetInteractionProcess(processName); + if (process != InteractionProcess::Transportation) { + gammaStatistics.crystalInteractions++; + switch (process) { + case InteractionProcess::Compton: + gammaStatistics.crystalCompton++; + break; + case InteractionProcess::Rayleigh: + gammaStatistics.crystalRayleigh++; + break; + default: + break; + } + } + } + + /** + * @brief Accumulates phantom-side interaction statistics. + * + * Args: + * processName: Geant4 process name at phantom hit. + * gammaStatistics: Statistics object to update. + */ + inline void AccumulatePhantom(const std::string_view processName, GammaStatistics& gammaStatistics) { + const InteractionProcess process = GetInteractionProcess(processName); + if (process != InteractionProcess::Transportation) { + gammaStatistics.phantomInteractions++; + switch (process) { + case InteractionProcess::Compton: + gammaStatistics.phantomCompton++; + break; + case InteractionProcess::Rayleigh: + gammaStatistics.phantomRayleigh++; + break; + default: + break; + } + } + } + +}// namespace MultiPhotonAnalysisHelpers + + +#endif // GateMultiPhotonAnalysisHelpers_h diff --git a/source/digits_hits/include/GateMultiPhotonAnalysisMessenger.hh b/source/digits_hits/include/GateMultiPhotonAnalysisMessenger.hh new file mode 100644 index 000000000..8d3734b38 --- /dev/null +++ b/source/digits_hits/include/GateMultiPhotonAnalysisMessenger.hh @@ -0,0 +1,61 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +#ifndef GateMultiPhotonAnalysisMessenger_h +#define GateMultiPhotonAnalysisMessenger_h 1 + +#include "GateOutputModuleMessenger.hh" + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About: UI messenger for configuring and controlling the multi-photon analysis output module. + **/ + +class G4UIcmdWithAString; + +class GateMultiPhotonAnalysis; + +/** + * @brief Messenger for GateMultiPhotonAnalysis output module. + * + * The messenger inherits standard output-module commands (`enable`, `disable`, + * `verbose`, `describe`) and adds `setMissingTrajectoryPolicy` for configuring + * how missing trajectory containers are handled. + */ +class GateMultiPhotonAnalysisMessenger : public GateOutputModuleMessenger { + public: + /** + * @brief Creates a messenger bound to GateMultiPhotonAnalysis. + * + * Args: + * gateMultiPhotonAnalysis: Owner module. + */ + explicit GateMultiPhotonAnalysisMessenger(GateMultiPhotonAnalysis *gateMultiPhotonAnalysis); + + ~GateMultiPhotonAnalysisMessenger(); + + /** + * @brief Handles UI command updates. + * + * Args: + * command: UI command. + * newValue: Command value. + */ + virtual void SetNewValue(G4UIcommand *command, G4String newValue); + + protected: + /** @brief Bound multi-photon analysis module instance. */ + GateMultiPhotonAnalysis *m_gateMultiPhotonAnalysis; + + /** @brief UI command for setting missing trajectory policy. */ + G4UIcmdWithAString *m_setMissingTrajectoryPolicyCmd; +}; + +#endif diff --git a/source/digits_hits/include/GateMultiPhotonTrajectoryNavigator.hh b/source/digits_hits/include/GateMultiPhotonTrajectoryNavigator.hh new file mode 100644 index 000000000..bf2d71a41 --- /dev/null +++ b/source/digits_hits/include/GateMultiPhotonTrajectoryNavigator.hh @@ -0,0 +1,139 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +#ifndef GateMultiPhotonTrajectoryNavigator_h +#define GateMultiPhotonTrajectoryNavigator_h + +#include "G4ThreeVector.hh" + +#include +#include +#include + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About: Event-level trajectory index and ancestry resolution API for multi-photon analysis. + **/ + +class G4TrajectoryContainer; + +/** + * @brief Resolves multi-photon trajectory ancestry for one event. + * + * This class builds per-event lookup tables that map track IDs to parent IDs + * and provides cached ancestry queries for ancestor photon and primary track + * resolution. + */ +class GateMultiPhotonTrajectoryNavigator { + public: + /** + * @brief Constructs an empty event-scoped navigator. + */ + GateMultiPhotonTrajectoryNavigator(); + + virtual ~GateMultiPhotonTrajectoryNavigator(); + + /** + * @brief Sets trajectory container for current event and resets caches. + * + * Args: + * trajectoryContainer: Event trajectory container. + */ + void SetTrajectoryContainer(G4TrajectoryContainer *trajectoryContainer); + + /** + * @brief Builds all lookup indices for current trajectory container. + */ + void BuildIndex(); + + /** + * @brief Returns source position inferred from source track. + * + * Returns: + * Source position when available, or default vector otherwise. + */ + G4ThreeVector FindSourcePosition() const; + + /** + * @brief Returns list of reference photon track IDs for current event. + * + * Returns: + * Reference photon track identifiers. + */ + std::vector FindReferencePhotonTrackIDs() const; + + /** + * @brief Resolves ancestor reference photon for a track. + * + * Args: + * trackID: Track identifier to resolve. + * + * Returns: + * Ancestor reference photon track ID, or 0 if unresolved. + */ + int FindAncestorPhotonTrackID(int trackID) const; + + /** + * @brief Resolves primary track for a track. + * + * Args: + * trackID: Track identifier to resolve. + * + * Returns: + * Primary track ID, or 0 if unresolved. + */ + int FindPrimaryTrackID(int trackID) const; + + /** + * @brief Sets verbose level used for diagnostic warnings. + * + * Args: + * v: Verbose level. + */ + void SetVerboseLevel(int v); + + private: + /** + * @brief Clears all internal caches and indices. + */ + void Reset(); + + /** + * @brief Builds reference photon set according to legacy-compatible policy. + */ + void BuildReferencePhotonSet(); + + /** @brief PDG code for positron. */ + const int kPositronPDG = -11; + + /** @brief PDG code for photon. */ + const int kPhotonPDG = 22; + + G4TrajectoryContainer *m_tc; + + std::unordered_map m_parentByTrack; + std::unordered_map m_pdgByTrack; + std::unordered_map m_chargeByTrack; + std::unordered_set m_referencePhotonTrackIDs; + + mutable std::unordered_map m_ancestorPhotonCache; + mutable std::unordered_map m_primaryCache; + + int m_positronTrackID; + int m_ionID; + int m_verboseLevel; + + bool m_hasSourcePosition; + double m_sourceX; + double m_sourceY; + double m_sourceZ; +}; + +#endif diff --git a/source/digits_hits/include/GateMultiPhotonTrajectoryNavigatorHelpers.hh b/source/digits_hits/include/GateMultiPhotonTrajectoryNavigatorHelpers.hh new file mode 100644 index 000000000..41fe708b2 --- /dev/null +++ b/source/digits_hits/include/GateMultiPhotonTrajectoryNavigatorHelpers.hh @@ -0,0 +1,212 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +#ifndef GateMultiPhotonTrajectoryNavigatorHelpers_h +#define GateMultiPhotonTrajectoryNavigatorHelpers_h + +#include +#include +#include + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About: Pure helper algorithms for resolving ancestor and primary tracks from parent-track maps. + **/ + +namespace MultiphotonTrajectoryResolver { + +/** + * @brief Resolves ancestor and primary track identifiers from parent links. + * + * This helper is intentionally independent from Geant4 trajectory classes so + * that algorithmic behavior can be unit-tested with synthetic graphs. + */ + /** + * @brief Marks all visited tracks as unresolved in cache. + * + * Args: + * visited: Traversed track IDs. + * cache: Optional cache to update. + */ + inline void clean_cache(const std::vector &visited, std::unordered_map *cache) { + if (cache) { + for (const int v : visited) { + (*cache)[v] = 0; + } + } + } + + /** + * @brief Tries to resolve current track using cache and propagates result. + * + * Args: + * current: Current track ID being resolved. + * visited: Traversed track IDs. + * cache: Optional memoization cache. + * resolved: Output resolved track ID. + * + * Returns: + * True when cached resolution was found. + */ + inline bool try_get_cached_resolution( + int current, + const std::vector &visited, + std::unordered_map *cache, + int *resolved) { + if (!cache || !resolved) { + return false; + } + + const auto it_cache = cache->find(current); + if (it_cache == cache->end()) { + return false; + } + + *resolved = it_cache->second; + for (const int v : visited) { + (*cache)[v] = *resolved; + } + return true; + } + + /** + * @brief Writes resolved track ID for current and visited nodes. + * + * Args: + * current: Current track ID. + * visited: Traversed track IDs. + * cache: Optional memoization cache. + * resolved: Resolved track ID to store. + */ + inline void cache_resolution_for_current_and_visited( + int current, + const std::vector &visited, + std::unordered_map *cache, + int resolved) { + if (!cache) { + return; + } + + (*cache)[current] = resolved; + for (const int v : visited) { + (*cache)[v] = resolved; + } + } + + /** + * @brief Resolves the nearest ancestor belonging to a reference photon set. + * + * Args: + * track_id: Start track identifier to resolve. + * parent_by_track: Mapping track -> parent track. + * reference_photon_track_ids: Set of accepted ancestor photon tracks. + * cache: Optional memoization cache for resolved ancestors. + * max_steps: Optional traversal guard; <=0 enables auto-guard. + * + * Returns: + * Matching ancestor photon track ID, or 0 when unresolved. + */ + inline int ResolveAncestorPhotonTrackID( + int track_id, + const std::unordered_map &parent_by_track, + const std::unordered_set &reference_photon_track_ids, + std::unordered_map *cache, + int max_steps = -1) { + if (track_id <= 0) { + return 0; + } + + std::vector visited; + int current = track_id; + int resolved = 0; + const int guard = max_steps > 0 ? max_steps : static_cast(parent_by_track.size()) + 2; + + for (int step = 0; step < guard; ++step) { + if (try_get_cached_resolution(current, visited, cache, &resolved)) { + return resolved; + } + + if (reference_photon_track_ids.find(current) != reference_photon_track_ids.end()) { + cache_resolution_for_current_and_visited(current, visited, cache, current); + return current; + } + + visited.push_back(current); + + const auto it_parent = parent_by_track.find(current); + if (it_parent == parent_by_track.end()) { + break; + } + + current = it_parent->second; + if (current == 0) { + break; + } + } + + clean_cache(visited, cache); + + return 0; + } + + /** + * @brief Resolves the primary track ID for a given track. + * + * Args: + * track_id: Start track identifier to resolve. + * parent_by_track: Mapping track -> parent track. + * cache: Optional memoization cache for primary track IDs. + * max_steps: Optional traversal guard; <=0 enables auto-guard. + * + * Returns: + * Primary track ID, or 0 when unresolved. + */ + inline int ResolvePrimaryTrackID( + int track_id, + const std::unordered_map &parent_by_track, + std::unordered_map *cache, + int max_steps = -1) { + if (track_id <= 0) { + return 0; + } + + std::vector visited; + int current = track_id; + int resolved = 0; + const int guard = max_steps > 0 ? max_steps : static_cast(parent_by_track.size()) + 2; + + for (int step = 0; step < guard; ++step) { + if (try_get_cached_resolution(current, visited, cache, &resolved)) { + return resolved; + } + + visited.push_back(current); + + const auto it_parent = parent_by_track.find(current); + if (it_parent == parent_by_track.end()) { + break; + } + + const int parent = it_parent->second; + if (parent == 0) { + cache_resolution_for_current_and_visited(current, visited, cache, current); + return current; + } + + current = parent; + } + + clean_cache(visited, cache); + + return 0; + } +} // namespace MultiphotonTrajectoryResolver + +#endif diff --git a/source/digits_hits/include/GateOutputMgr.hh b/source/digits_hits/include/GateOutputMgr.hh index 854313bd2..a7e2b7a12 100644 --- a/source/digits_hits/include/GateOutputMgr.hh +++ b/source/digits_hits/include/GateOutputMgr.hh @@ -19,6 +19,16 @@ #include "GateDigi.hh" #include "GateCoincidenceDigi.hh" +/** + * @brief Evaluates if at least one analysis-capable output module is enabled. + */ +inline bool IsAnyAnalysisModuleEnabled( + bool analysis_enabled, + bool fastanalysis_enabled, + bool multianalysis_enabled) { + return analysis_enabled || fastanalysis_enabled || multianalysis_enabled; +} + class G4Run; class G4Step; class G4Event; diff --git a/source/digits_hits/include/GateToTree.hh b/source/digits_hits/include/GateToTree.hh index 46de07d5c..f04aa78c1 100644 --- a/source/digits_hits/include/GateToTree.hh +++ b/source/digits_hits/include/GateToTree.hh @@ -201,6 +201,7 @@ private: G4int m_sourceType = 0; G4int m_decayType = 0; G4int m_gammaType = 0; + G4int m_decayIndex = -1; G4float m_sourceEnergy; G4int m_sourcePDG; diff --git a/source/digits_hits/src/GateCrystalSD.cc b/source/digits_hits/src/GateCrystalSD.cc index 0e753befa..63d42e943 100644 --- a/source/digits_hits/src/GateCrystalSD.cc +++ b/source/digits_hits/src/GateCrystalSD.cc @@ -144,6 +144,7 @@ G4bool GateCrystalSD::ProcessHits(G4Step*aStep, G4TouchableHistory*) G4int source_type = static_cast(GateEmittedGammaInformation::SourceKind::NotDefined); G4int decay_type = static_cast(GateEmittedGammaInformation::DecayModel::None); G4int gamma_type = static_cast(GateEmittedGammaInformation::GammaKind::Unknown); + G4int decay_index = -1; G4PrimaryParticle* primary_particle = aTrack->GetDynamicParticle()->GetPrimaryParticle(); if( primary_particle != nullptr ) { @@ -153,6 +154,7 @@ G4bool GateCrystalSD::ProcessHits(G4Step*aStep, G4TouchableHistory*) source_type = static_cast( info->GetSourceKind() ); decay_type = static_cast( info->GetDecayModel() ); gamma_type = static_cast( info->GetGammaKind() ); + decay_index = info->GetDecayIndex(); } } @@ -289,6 +291,7 @@ G4bool GateCrystalSD::ProcessHits(G4Step*aStep, G4TouchableHistory*) aHit->SetSourceType( source_type ); aHit->SetDecayType( decay_type ); aHit->SetGammaType( gamma_type ); + aHit->SetDecayIndex( decay_index ); // OK GND 2023 for CC diff --git a/source/digits_hits/src/GateHit.cc b/source/digits_hits/src/GateHit.cc index b152bbf1b..99d307d83 100644 --- a/source/digits_hits/src/GateHit.cc +++ b/source/digits_hits/src/GateHit.cc @@ -30,6 +30,7 @@ GateHit::GateHit() m_systemID(-1), m_sourceEnergy(-1), m_sourcePDG(0), + m_nInteractions(-1), m_nCrystalConv(0) {;} //--------------------------------------------------------------------- diff --git a/source/digits_hits/src/GateMultiPhotonAnalysis.cc b/source/digits_hits/src/GateMultiPhotonAnalysis.cc new file mode 100644 index 000000000..36417e26c --- /dev/null +++ b/source/digits_hits/src/GateMultiPhotonAnalysis.cc @@ -0,0 +1,434 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +#include "GateMultiPhotonAnalysis.hh" + +#include "GateActions.hh" +#include "GateDigitizerMgr.hh" +#include "GateHit.hh" +#include "GateMultiPhotonAnalysisMessenger.hh" +#include "GateMultiPhotonTrajectoryNavigator.hh" +#include "GateOutputMgr.hh" +#include "GatePhantomHit.hh" +#include "GateRunManager.hh" +#include "GateSourceMgr.hh" + +#include "G4Event.hh" +#include "G4HCofThisEvent.hh" +#include "G4Run.hh" +#include "G4TrajectoryContainer.hh" + +#include +#include +#include +#include +#include + +namespace { + +struct TimelineEntry { + double time = 0.0; + std::size_t sequence = 0; + bool is_phantom = false; + int track_id = 0; + int ancestor_photon = 0; + GatePhantomHit *phantom_hit = 0; + GateHit *crystal_hit = 0; +}; + +struct EventContext { + G4int event_id = 0; + G4int run_id = 0; + G4int source_id = -1; + G4ThreeVector source_vertex = G4ThreeVector(-1, -1, -1); +}; + +bool EventHasProcessableHits(const std::vector &CHC_vector, + GatePhantomHitsCollection *PHC) { + if (PHC && PHC->entries() > 0) { + return true; + } + + for (std::size_t i = 0; i < CHC_vector.size(); ++i) { + GateHitsCollection *CHC = CHC_vector[i]; + if (CHC && CHC->entries() > 0) { + return true; + } + } + + return false; +} + +std::vector BuildBasePhantomTimeline( + GatePhantomHitsCollection *PHC, + GateMultiPhotonTrajectoryNavigator *trajectoryNavigator) { + std::vector basePhantomTimeline; + if (!PHC) { + return basePhantomTimeline; + } + + std::size_t sequence = 0; + const G4int NpHits = PHC->entries(); + basePhantomTimeline.reserve(static_cast(NpHits)); + + for (G4int iPHit = 0; iPHit < NpHits; ++iPHit) { + GatePhantomHit *phantomHit = (*PHC)[iPHit]; + if (!phantomHit) { + continue; + } + + const G4int trackID = phantomHit->GetTrackID(); + const int ancestorPhoton = trajectoryNavigator->FindAncestorPhotonTrackID(trackID); + if (ancestorPhoton == 0) { + continue; + } + + TimelineEntry entry; + entry.time = phantomHit->GetTime(); + entry.sequence = sequence++; + entry.is_phantom = true; + entry.track_id = trackID; + entry.ancestor_photon = ancestorPhoton; + entry.phantom_hit = phantomHit; + basePhantomTimeline.push_back(entry); + } + + return basePhantomTimeline; +} + +EventContext BuildEventContext( + const G4Event *event, + GateRunManager *runManager, + GateMultiPhotonTrajectoryNavigator *trajectoryNavigator) { + EventContext context; + context.event_id = event->GetEventID(); + context.run_id = runManager->GetCurrentRun()->GetRunID(); + + GateSourceMgr *sourceMgr = GateSourceMgr::GetInstance(); + const std::vector &sourcesForThisEvent = sourceMgr->GetSourcesForThisEvent(); + if (!sourcesForThisEvent.empty()) { + context.source_id = sourcesForThisEvent[0]->GetSourceID(); + context.source_vertex = trajectoryNavigator->FindSourcePosition(); + } + + return context; +} + +std::vector BuildTimelineForCrystalCollection( + const std::vector &basePhantomTimeline, + GateHitsCollection *CHC, + GateMultiPhotonTrajectoryNavigator *trajectoryNavigator) { + std::vector timeline = basePhantomTimeline; + std::size_t sequence = timeline.size(); + + const G4int NbHits = CHC->entries(); + timeline.reserve(basePhantomTimeline.size() + static_cast(NbHits)); + for (G4int iHit = 0; iHit < NbHits; ++iHit) { + GateHit *crystalHit = (*CHC)[iHit]; + if (!crystalHit) { + continue; + } + + const G4int trackID = crystalHit->GetTrackID(); + const int ancestorPhoton = trajectoryNavigator->FindAncestorPhotonTrackID(trackID); + if (ancestorPhoton == 0) { + continue; + } + + TimelineEntry entry; + entry.time = crystalHit->GetTime(); + entry.sequence = sequence++; + entry.is_phantom = false; + entry.track_id = trackID; + entry.ancestor_photon = ancestorPhoton; + entry.crystal_hit = crystalHit; + timeline.push_back(entry); + } + + std::sort(timeline.begin(), timeline.end(), [](const TimelineEntry &a, const TimelineEntry &b) { + if (a.time < b.time) { + return true; + } + if (a.time > b.time) { + return false; + } + return a.sequence < b.sequence; + }); + + return timeline; +} + +void ProcessTimeline( + std::vector *timeline, + GateMultiPhotonTrajectoryNavigator *trajectoryNavigator, + const EventContext &context, + int legacyPhotonIDPolicy) { + std::unordered_map runningStatsByPhotonTrackId; + runningStatsByPhotonTrackId.reserve(timeline->size()); + + for (std::size_t idx = 0; idx < timeline->size(); ++idx) { + TimelineEntry &entry = (*timeline)[idx]; + MultiPhotonAnalysisHelpers::GammaStatistics &runningStats = runningStatsByPhotonTrackId[entry.ancestor_photon]; + + if (!entry.is_phantom && entry.crystal_hit && entry.crystal_hit->GoodForAnalysis()) { + GateHit *hit = entry.crystal_hit; + const int primaryID = trajectoryNavigator->FindPrimaryTrackID(entry.track_id); + const int nInteractions = runningStats.phantomInteractions + runningStats.crystalInteractions; + + hit->SetSourceID(context.source_id); + hit->SetSourcePosition(context.source_vertex); + hit->SetNPhantomCompton(runningStats.phantomCompton); + hit->SetNPhantomRayleigh(runningStats.phantomRayleigh); + if (runningStats.comptonVolumeName != "NULL") { + hit->SetComptonVolumeName(runningStats.comptonVolumeName.c_str()); + } + if (runningStats.rayleighVolumeName != "NULL") { + hit->SetRayleighVolumeName(runningStats.rayleighVolumeName.c_str()); + } + hit->SetPhotonID(legacyPhotonIDPolicy); + hit->SetPrimaryID(primaryID); + hit->SetEventID(context.event_id); + hit->SetRunID(context.run_id); + hit->SetNCrystalCompton(runningStats.crystalCompton); + hit->SetNCrystalRayleigh(runningStats.crystalRayleigh); + // nInteractions is intentionally filled only in the multiphoton analysis path. + hit->SetNInteractions(nInteractions); + } + + if (entry.is_phantom) { + if (!entry.phantom_hit) { + continue; + } + + MultiPhotonAnalysisHelpers::AccumulatePhantom(entry.phantom_hit->GetProcess(), runningStats); + continue; + } + + if (entry.crystal_hit) { + MultiPhotonAnalysisHelpers::AccumulateCrystal(entry.crystal_hit->GetProcess(), runningStats); + } + } +} + +void RunDigitizersIfNeeded() { + GateDigitizerMgr *digitizerMgr = GateDigitizerMgr::GetInstance(); + if (!digitizerMgr->m_alreadyRun) { + if (digitizerMgr->m_recordSingles || digitizerMgr->m_recordCoincidences) { + digitizerMgr->RunDigitizers(); + digitizerMgr->RunCoincidenceSorters(); + digitizerMgr->RunCoincidenceDigitizers(); + } + } +} + +} // namespace + +GateMultiPhotonAnalysis::GateMultiPhotonAnalysis(const G4String &name, GateOutputMgr *outputMgr, DigiMode digiMode) + : GateVOutputModule(name, outputMgr, digiMode), + m_trajectoryNavigator(new GateMultiPhotonTrajectoryNavigator()), + m_messenger(new GateMultiPhotonAnalysisMessenger(this)), + m_missingTrajectoryPolicy(kResilient), + m_missingTrajectoryEventCount(0), + m_missingTrajectoryWithHitsCount(0) { + m_isEnabled = false; + SetVerboseLevel(0); +} + +GateMultiPhotonAnalysis::~GateMultiPhotonAnalysis() { + delete m_messenger; + m_messenger = 0; + delete m_trajectoryNavigator; + m_trajectoryNavigator = 0; +} + +const G4String &GateMultiPhotonAnalysis::GiveNameOfFile() { + m_noFileName = " "; + return m_noFileName; +} + +void GateMultiPhotonAnalysis::RecordBeginOfAcquisition() { + if (nVerboseLevel > 2) { + G4cout << "GateMultiPhotonAnalysis::RecordBeginOfAcquisition" << G4endl; + } +} + +void GateMultiPhotonAnalysis::RecordEndOfAcquisition() { + if (nVerboseLevel > 2) { + G4cout << "GateMultiPhotonAnalysis::RecordEndOfAcquisition" << G4endl; + } +} + +void GateMultiPhotonAnalysis::RecordBeginOfRun(const G4Run *) { + if (nVerboseLevel > 2) { + G4cout << "GateMultiPhotonAnalysis::RecordBeginOfRun" << G4endl; + } +} + +void GateMultiPhotonAnalysis::RecordEndOfRun(const G4Run *) { + if (nVerboseLevel > 2) { + G4cout << "GateMultiPhotonAnalysis::RecordEndOfRun" << G4endl; + } + + if (m_missingTrajectoryEventCount > 0) { + G4cout + << "[GateMultiPhotonAnalysis] Missing trajectory container in " + << m_missingTrajectoryEventCount + << " event(s), including " + << m_missingTrajectoryWithHitsCount + << " event(s) with processable hits. Policy=" + << GetMissingTrajectoryPolicyName() + << G4endl; + } +} + +void GateMultiPhotonAnalysis::RecordBeginOfEvent(const G4Event *) { + if (nVerboseLevel > 2) { + G4cout << "GateMultiPhotonAnalysis::RecordBeginOfEvent" << G4endl; + } +} + +void GateMultiPhotonAnalysis::RecordEndOfEvent(const G4Event *event) { + if (!event) { + return; + } + + GateRunManager *runManager = GateRunManager::GetRunManager(); + GateSteppingAction *steppingAction = (GateSteppingAction *)(runManager->GetUserSteppingAction()); + TrackingMode mode = steppingAction->GetMode(); + const int tracking_mode_code = static_cast(mode); + if (!IsTrackingModeSupported(tracking_mode_code)) { + G4Exception( + "GateMultiPhotonAnalysis::RecordEndOfEvent", + "GateMultiPhotonAnalysisUnsupportedTrackingMode", + FatalException, + "GateMultiPhotonAnalysis cannot process the current tracking mode. " + "Only TrackingMode::kBoth is supported, and continuing would skip the " + "remaining end-of-event processing, including digitizer output."); + return; + } + + std::vector CHC_vector = GetOutputMgr()->GetHitCollections(); + GatePhantomHitsCollection *PHC = GetOutputMgr()->GetPhantomHitCollection(); + + G4TrajectoryContainer *trajectoryContainer = event->GetTrajectoryContainer(); + if (!trajectoryContainer) { + ++m_missingTrajectoryEventCount; + + const bool hasProcessableHits = EventHasProcessableHits(CHC_vector, PHC); + if (hasProcessableHits) { + ++m_missingTrajectoryWithHitsCount; + } + + G4int runID = -1; + if (runManager->GetCurrentRun()) { + runID = runManager->GetCurrentRun()->GetRunID(); + } + const G4int eventID = event->GetEventID(); + + G4String details = "GateMultiPhotonAnalysis missing trajectory container for run=" + + std::to_string(runID) + + ", event=" + + std::to_string(eventID) + + ". hasProcessableHits=" + + (hasProcessableHits ? "true" : "false") + + ". Policy=" + + GetMissingTrajectoryPolicyName() + + "."; + + if (m_missingTrajectoryPolicy == kStrict) { + G4Exception( + "GateMultiPhotonAnalysis::RecordEndOfEvent", + "GateMultiPhotonAnalysisMissingTrajectoryContainer", + FatalException, + details.c_str()); + return; + } + + G4Exception( + "GateMultiPhotonAnalysis::RecordEndOfEvent", + "GateMultiPhotonAnalysisMissingTrajectoryContainer", + JustWarning, + details.c_str()); + return; + } + + m_trajectoryNavigator->SetTrajectoryContainer(trajectoryContainer); + m_trajectoryNavigator->BuildIndex(); + std::vector basePhantomTimeline = BuildBasePhantomTimeline(PHC, m_trajectoryNavigator); + const EventContext context = BuildEventContext(event, runManager, m_trajectoryNavigator); + const int legacyPhotonIDPolicy = ResolveLegacyPhotonIDPolicy(); + + for (size_t i = 0; i < CHC_vector.size(); ++i) { + GateHitsCollection *CHC = CHC_vector[i]; + if (!CHC) { + continue; + } + + std::vector timeline = BuildTimelineForCrystalCollection( + basePhantomTimeline, + CHC, + m_trajectoryNavigator); + ProcessTimeline(&timeline, m_trajectoryNavigator, context, legacyPhotonIDPolicy); + } + + RunDigitizersIfNeeded(); +} + +void GateMultiPhotonAnalysis::RecordStepWithVolume(const GateVVolume *, const G4Step *) { + if (nVerboseLevel > 2) { + G4cout << "GateMultiPhotonAnalysis::RecordStepWithVolume" << G4endl; + } +} + +void GateMultiPhotonAnalysis::SetVerboseLevel(G4int val) { + nVerboseLevel = val; + if (m_trajectoryNavigator) { + m_trajectoryNavigator->SetVerboseLevel(val); + } +} + +bool GateMultiPhotonAnalysis::IsTrackingModeSupported(int tracking_mode_code) const { + return tracking_mode_code == static_cast(TrackingMode::kBoth); +} + +int GateMultiPhotonAnalysis::ResolveLegacyPhotonIDPolicy() const { return 0; } + +void GateMultiPhotonAnalysis::SetMissingTrajectoryPolicy(MissingTrajectoryPolicy policy) { + m_missingTrajectoryPolicy = policy; +} + +void GateMultiPhotonAnalysis::SetMissingTrajectoryPolicyFromString(const G4String &policyName) { + if (policyName == "strict") { + m_missingTrajectoryPolicy = kStrict; + return; + } + if (policyName == "resilient") { + m_missingTrajectoryPolicy = kResilient; + return; + } + + G4String message = "Unsupported missing trajectory policy: " + policyName + + ". Expected one of: strict resilient."; + G4Exception( + "GateMultiPhotonAnalysis::SetMissingTrajectoryPolicyFromString", + "GateMultiPhotonAnalysisInvalidPolicy", + JustWarning, + message.c_str()); +} + +GateMultiPhotonAnalysis::MissingTrajectoryPolicy GateMultiPhotonAnalysis::GetMissingTrajectoryPolicy() const { + return m_missingTrajectoryPolicy; +} + +G4String GateMultiPhotonAnalysis::GetMissingTrajectoryPolicyName() const { + if (m_missingTrajectoryPolicy == kStrict) { + return "strict"; + } + return "resilient"; +} + diff --git a/source/digits_hits/src/GateMultiPhotonAnalysisMessenger.cc b/source/digits_hits/src/GateMultiPhotonAnalysisMessenger.cc new file mode 100644 index 000000000..21eaf3dd1 --- /dev/null +++ b/source/digits_hits/src/GateMultiPhotonAnalysisMessenger.cc @@ -0,0 +1,39 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +#include "GateMultiPhotonAnalysisMessenger.hh" + +#include "GateMultiPhotonAnalysis.hh" + +#include "G4UIcmdWithAString.hh" + +GateMultiPhotonAnalysisMessenger::GateMultiPhotonAnalysisMessenger( + GateMultiPhotonAnalysis *gateMultiPhotonAnalysis) + : GateOutputModuleMessenger(gateMultiPhotonAnalysis), + m_gateMultiPhotonAnalysis(gateMultiPhotonAnalysis), + m_setMissingTrajectoryPolicyCmd(0) { + G4String cmdName = GetDirectoryName() + "setMissingTrajectoryPolicy"; + m_setMissingTrajectoryPolicyCmd = new G4UIcmdWithAString(cmdName, this); + m_setMissingTrajectoryPolicyCmd->SetGuidance("Set policy for missing trajectory container handling."); + m_setMissingTrajectoryPolicyCmd->SetGuidance("Options: strict resilient"); + m_setMissingTrajectoryPolicyCmd->SetCandidates("strict resilient"); + m_setMissingTrajectoryPolicyCmd->SetParameterName("policy", false); +} + +GateMultiPhotonAnalysisMessenger::~GateMultiPhotonAnalysisMessenger() { + delete m_setMissingTrajectoryPolicyCmd; +} + +void GateMultiPhotonAnalysisMessenger::SetNewValue(G4UIcommand *command, G4String newValue) { + if (command == m_setMissingTrajectoryPolicyCmd) { + m_gateMultiPhotonAnalysis->SetMissingTrajectoryPolicyFromString(newValue); + return; + } + + GateOutputModuleMessenger::SetNewValue(command, newValue); +} diff --git a/source/digits_hits/src/GateMultiPhotonTrajectoryNavigator.cc b/source/digits_hits/src/GateMultiPhotonTrajectoryNavigator.cc new file mode 100644 index 000000000..70cdefb58 --- /dev/null +++ b/source/digits_hits/src/GateMultiPhotonTrajectoryNavigator.cc @@ -0,0 +1,180 @@ +/*---------------------- + Copyright (C): OpenGATE Collaboration + +This software is distributed under the terms +of the GNU Lesser General Public Licence (LGPL) +See LICENSE.md for further details +----------------------*/ + +#include "GateMultiPhotonTrajectoryNavigator.hh" + +#include "GateMultiPhotonTrajectoryNavigatorHelpers.hh" + +#include "G4ThreeVector.hh" +#include "G4Trajectory.hh" +#include "G4TrajectoryContainer.hh" + +#include "G4ios.hh" + +GateMultiPhotonTrajectoryNavigator::GateMultiPhotonTrajectoryNavigator() + : m_tc(0), + m_positronTrackID(0), + m_ionID(0), + m_verboseLevel(0), + m_hasSourcePosition(false), + m_sourceX(0.0), + m_sourceY(0.0), + m_sourceZ(0.0) {} + +GateMultiPhotonTrajectoryNavigator::~GateMultiPhotonTrajectoryNavigator() {} + +void GateMultiPhotonTrajectoryNavigator::SetTrajectoryContainer(G4TrajectoryContainer *trajectoryContainer) { + m_tc = trajectoryContainer; + Reset(); +} + +void GateMultiPhotonTrajectoryNavigator::BuildIndex() { + Reset(); + + if (!m_tc) { + if (m_verboseLevel > 0) { + G4cout << "GateMultiPhotonTrajectoryNavigator::BuildIndex: WARNING null trajectory container" << G4endl; + } + return; + } + + const int n_trajectories = m_tc->entries(); + const std::size_t expected_track_count = static_cast(n_trajectories); + m_parentByTrack.reserve(expected_track_count); + m_pdgByTrack.reserve(expected_track_count); + m_chargeByTrack.reserve(expected_track_count); + m_referencePhotonTrackIDs.reserve(expected_track_count); + m_ancestorPhotonCache.reserve(expected_track_count); + m_primaryCache.reserve(expected_track_count); + + for (int i = 0; i < n_trajectories; ++i) { + G4Trajectory *trj = (G4Trajectory *)((*m_tc)[i]); + if (!trj) { + continue; + } + + const int track_id = trj->GetTrackID(); + m_parentByTrack[track_id] = trj->GetParentID(); + m_pdgByTrack[track_id] = trj->GetPDGEncoding(); + m_chargeByTrack[track_id] = trj->GetCharge(); + + if (track_id == 1 && trj->GetPointEntries() > 0) { + const G4ThreeVector p = ((G4TrajectoryPoint *)(trj->GetPoint(0)))->GetPosition(); + m_sourceX = p.x(); + m_sourceY = p.y(); + m_sourceZ = p.z(); + m_hasSourcePosition = true; + } + } + + for (std::unordered_map::const_iterator it = m_chargeByTrack.begin(); it != m_chargeByTrack.end(); ++it) { + if (it->second > 2.0) { + m_ionID = 1; + break; + } + } + + for (std::unordered_map::const_iterator it = m_pdgByTrack.begin(); it != m_pdgByTrack.end(); ++it) { + const int track_id = it->first; + const int pdg = it->second; + const std::unordered_map::const_iterator parent_it = m_parentByTrack.find(track_id); + if (parent_it == m_parentByTrack.end()) { + continue; + } + if (parent_it->second == m_ionID && pdg == kPositronPDG) { + m_positronTrackID = track_id; + break; + } + } + + BuildReferencePhotonSet(); +} + +G4ThreeVector GateMultiPhotonTrajectoryNavigator::FindSourcePosition() const { + if (!m_hasSourcePosition) { + return G4ThreeVector(); + } + return G4ThreeVector(m_sourceX, m_sourceY, m_sourceZ); +} + +std::vector GateMultiPhotonTrajectoryNavigator::FindReferencePhotonTrackIDs() const { + std::vector out; + out.reserve(m_referencePhotonTrackIDs.size()); + for (std::unordered_set::const_iterator it = m_referencePhotonTrackIDs.begin(); it != m_referencePhotonTrackIDs.end(); ++it) { + out.push_back(*it); + } + return out; +} + +int GateMultiPhotonTrajectoryNavigator::FindAncestorPhotonTrackID(int trackID) const { + return MultiphotonTrajectoryResolver::ResolveAncestorPhotonTrackID( + trackID, + m_parentByTrack, + m_referencePhotonTrackIDs, + &m_ancestorPhotonCache); +} + +int GateMultiPhotonTrajectoryNavigator::FindPrimaryTrackID(int trackID) const { + return MultiphotonTrajectoryResolver::ResolvePrimaryTrackID(trackID, m_parentByTrack, &m_primaryCache); +} + +void GateMultiPhotonTrajectoryNavigator::SetVerboseLevel(int v) { m_verboseLevel = v; } + +void GateMultiPhotonTrajectoryNavigator::Reset() { + m_parentByTrack.clear(); + m_pdgByTrack.clear(); + m_chargeByTrack.clear(); + m_referencePhotonTrackIDs.clear(); + m_ancestorPhotonCache.clear(); + m_primaryCache.clear(); + m_positronTrackID = 0; + m_ionID = 0; + m_hasSourcePosition = false; + m_sourceX = 0.0; + m_sourceY = 0.0; + m_sourceZ = 0.0; +} + +void GateMultiPhotonTrajectoryNavigator::BuildReferencePhotonSet() { + for (std::unordered_map::const_iterator it = m_pdgByTrack.begin(); it != m_pdgByTrack.end(); ++it) { + const int track_id = it->first; + const int pdg = it->second; + if (pdg != kPhotonPDG) { + continue; + } + + const std::unordered_map::const_iterator parent_it = m_parentByTrack.find(track_id); + if (parent_it == m_parentByTrack.end()) { + continue; + } + + const int parent_id = parent_it->second; + + if (m_positronTrackID != 0) { + if (parent_id == m_positronTrackID) { + m_referencePhotonTrackIDs.insert(track_id); + } + } else { + bool is_reference_photon = false; + if (parent_id == 0) { + is_reference_photon = true; + } else if (m_ionID != 0 && parent_id == m_ionID) { + is_reference_photon = true; + } else { + const std::unordered_map::const_iterator parent_pdg_it = m_pdgByTrack.find(parent_id); + if (parent_pdg_it != m_pdgByTrack.end() && parent_pdg_it->second == kPositronPDG) { + is_reference_photon = true; + } + } + + if (is_reference_photon) { + m_referencePhotonTrackIDs.insert(track_id); + } + } + } +} diff --git a/source/digits_hits/src/GateOutputMgr.cc b/source/digits_hits/src/GateOutputMgr.cc index e587b3f32..bf1c4b3f4 100644 --- a/source/digits_hits/src/GateOutputMgr.cc +++ b/source/digits_hits/src/GateOutputMgr.cc @@ -13,6 +13,7 @@ #include "GateOutputMgrMessenger.hh" #include "GateConfiguration.h" #include "GateAnalysis.hh" +#include "GateMultiPhotonAnalysis.hh" #ifdef GATE_USE_OPTICAL #include "GateFastAnalysis.hh" #endif @@ -93,6 +94,9 @@ GateOutputMgr::GateOutputMgr(const G4String name) GateAnalysis* gateAnalysis = new GateAnalysis("analysis", this,m_digiMode); AddOutputModule((GateVOutputModule*)gateAnalysis); + GateMultiPhotonAnalysis* gateMultiPhotonAnalysis = new GateMultiPhotonAnalysis("multianalysis", this, m_digiMode); + AddOutputModule((GateVOutputModule*)gateMultiPhotonAnalysis); + } @@ -262,11 +266,23 @@ void GateOutputMgr::RecordBeginOfAcquisition() G4cout << "GateOutputMgr::RecordBeginOfAcquisition\n"; //OK GND GateDigitizerMgr* digitizerMgr=GateDigitizerMgr::GetInstance(); - if((digitizerMgr->m_recordSingles|| digitizerMgr->m_recordCoincidences) - && !this->FindOutputModule("analysis")->IsEnabled() - && !this->FindOutputModule("fastanalysis")->IsEnabled()) + bool analysisEnabled = false; + bool fastanalysisEnabled = false; + bool multianalysisEnabled = false; + for (size_t iMod=0; iModGetName(); + if (moduleName == "analysis") { + analysisEnabled = m_outputModules[iMod]->IsEnabled(); + } else if (moduleName == "fastanalysis") { + fastanalysisEnabled = m_outputModules[iMod]->IsEnabled(); + } else if (moduleName == "multianalysis") { + multianalysisEnabled = m_outputModules[iMod]->IsEnabled(); + } + } + if((digitizerMgr->m_recordSingles|| digitizerMgr->m_recordCoincidences) + && !IsAnyAnalysisModuleEnabled(analysisEnabled, fastanalysisEnabled, multianalysisEnabled)) { - GateError("***ERROR*** Digitizer Manager is not initialized properly. Please, enable analysis or fastanalysis Output Modules to write down Singles or Coincidences.\n Use, /gate/output/analysis/enable or /gate/output/fastanalysis/enable.\n"); + GateError("***ERROR*** Digitizer Manager is not initialized properly. Please, enable analysis, fastanalysis or multianalysis Output Modules to write down Singles or Coincidences.\n Use, /gate/output/analysis/enable or /gate/output/fastanalysis/enable or /gate/output/multianalysis/enable.\n"); } diff --git a/source/digits_hits/src/GateToTree.cc b/source/digits_hits/src/GateToTree.cc index 0fa5671ef..51ae649d6 100644 --- a/source/digits_hits/src/GateToTree.cc +++ b/source/digits_hits/src/GateToTree.cc @@ -103,6 +103,7 @@ GateToTree::GateToTree(const G4String &name, GateOutputMgr *outputMgr, DigiMode m_hitsParams_to_write.emplace("sourceType", SaveDataParam()); m_hitsParams_to_write.emplace("decayType", SaveDataParam()); m_hitsParams_to_write.emplace("gammaType", SaveDataParam()); + m_hitsParams_to_write.emplace("decayIndex", SaveDataParam()); //CC m_hitsParams_to_write.emplace("sourceEnergy",SaveDataParam()); m_hitsParams_to_write.emplace("sourcePDG",SaveDataParam()); @@ -360,6 +361,9 @@ void GateToTree::RecordBeginOfAcquisition() { if (m_hitsParams_to_write.at("gammaType").toSave()) m_manager_hits.write_variable("gammaType", &m_gammaType); + if (m_hitsParams_to_write.at("decayIndex").toSave()) + m_manager_hits.write_variable("decayIndex", &m_decayIndex); + m_manager_hits.write_header(); @@ -517,6 +521,9 @@ void GateToTree::RecordBeginOfAcquisition() { if (m_hitsParams_to_write.at("gammaType").toSave()) mm.write_variable("gammaType", &m_gammaType); + if (m_hitsParams_to_write.at("decayIndex").toSave()) + mm.write_variable("decayIndex", &m_decayIndex); + if(m_cc_enabled) { if (m_hitsParams_to_write.at("sourceEnergy").toSave()) @@ -980,6 +987,7 @@ void GateToTree::RecordEndOfEvent(const G4Event *event) { m_sourceType = hit->GetSourceType(); m_decayType = hit->GetDecayType(); m_gammaType = hit->GetGammaType(); + m_decayIndex = hit->GetDecayIndex(); if(m_cc_enabled) { diff --git a/source/general/include/GateRootDefs.hh b/source/general/include/GateRootDefs.hh index cc199af8f..349b85548 100644 --- a/source/general/include/GateRootDefs.hh +++ b/source/general/include/GateRootDefs.hh @@ -222,6 +222,7 @@ class GateRootHitBuffer Int_t nCrystalCompton; //!< Number of Compton interactions in the crystam Int_t nPhantomRayleigh; //!< Number of Rayleigh interactions in the phantom Int_t nCrystalRayleigh; //!< Number of Rayleigh interactions in the crystam + Int_t nInteractions; //!< Number of non-Transportation interactions in phantom + crystal Int_t primaryID; //!< Primary ID Float_t sourcePosX,sourcePosY,sourcePosZ; //!< Global decay position (in millimeters) Int_t sourceID; //!< Source ID @@ -237,8 +238,7 @@ class GateRootHitBuffer Int_t sourceType = 0; //Type of gamma source (check ExtendedVSource) Int_t decayType = 0; //Type of positronium decay (check ExtendedVSource) Int_t gammaType = 0; //Gamma type - single, annhilation, prompt (check ExtendedVSource) - - + Int_t decayIndex = -1; //Decay channel index //OK GND for CC G4bool m_CCflag; Float_t sourceEnergy; diff --git a/source/general/src/GateRootDefs.cc b/source/general/src/GateRootDefs.cc index b5840963b..aa8028855 100644 --- a/source/general/src/GateRootDefs.cc +++ b/source/general/src/GateRootDefs.cc @@ -133,12 +133,14 @@ void GateRootHitBuffer::Clear() nCrystalCompton = -1; nPhantomRayleigh = -1; nCrystalRayleigh = -1; + nInteractions = -1; primaryID = -1; axialPos = 0.; rotationAngle = 0.; sourceType = 0; decayType = 0; gammaType = 0; + decayIndex = -1; strcpy (comptonVolumeName," "); strcpy (RayleighVolumeName," "); } @@ -207,6 +209,7 @@ void GateRootHitBuffer::Fill(GateHit* aHit) nCrystalCompton = aHit->GetNCrystalCompton(); nPhantomRayleigh = aHit->GetNPhantomRayleigh(); nCrystalRayleigh = aHit->GetNCrystalRayleigh(); + nInteractions = aHit->GetNInteractions(); primaryID = aHit->GetPrimaryID(); momDirX = aHit->GetMomentumDir().x(); momDirY = aHit->GetMomentumDir().y(); @@ -216,6 +219,7 @@ void GateRootHitBuffer::Fill(GateHit* aHit) sourceType = aHit->GetSourceType(); decayType = aHit->GetDecayType(); gammaType = aHit->GetGammaType(); + decayIndex = aHit->GetDecayIndex(); } else { @@ -271,6 +275,7 @@ GateHit* GateRootHitBuffer::CreateHit() aHit->SetNCrystalCompton( nCrystalCompton ); aHit->SetNPhantomRayleigh( nPhantomRayleigh); aHit->SetNCrystalRayleigh( nCrystalRayleigh ); + aHit->SetNInteractions( nInteractions ); aHit->SetComptonVolumeName( comptonVolumeName ); aHit->SetRayleighVolumeName( RayleighVolumeName ); aHit->SetPrimaryID( primaryID ); @@ -284,6 +289,7 @@ GateHit* GateRootHitBuffer::CreateHit() aHit->SetSourceType(sourceType); aHit->SetDecayType(decayType); aHit->SetGammaType(gammaType); + aHit->SetDecayIndex(decayIndex); aHit->SetSourceEnergy(GetSourceEnergy()); aHit->SetSourcePDG(GetSourcePDG()); @@ -341,6 +347,7 @@ void GateHitTree::Init(GateRootHitBuffer& buffer) Branch("nCrystalCompton",&buffer.nCrystalCompton,"nCrystalCompton/I"); Branch("nPhantomRayleigh",&buffer.nPhantomRayleigh,"nPhantomRayleigh/I"); Branch("nCrystalRayleigh",&buffer.nCrystalRayleigh,"nCrystalRayleigh/I"); + Branch("nInteractions",&buffer.nInteractions,"nInteractions/I"); Branch("primaryID", &buffer.primaryID,"primaryID/I"); Branch("axialPos", &buffer.axialPos,"axialPos/F"); @@ -354,6 +361,7 @@ void GateHitTree::Init(GateRootHitBuffer& buffer) Branch("sourceType", &buffer.sourceType,"sourceType/I"); Branch("decayType", &buffer.decayType,"decayType/I"); Branch("gammaType", &buffer.gammaType,"gammaType/I"); + Branch("decayIndex", &buffer.decayIndex,"decayIndex/I"); } else { @@ -376,6 +384,7 @@ void GateHitTree::SetBranchAddresses(TTree* hitTree,GateRootHitBuffer& buffer) hitTree->SetBranchAddress("PDGEncoding",&buffer.PDGEncoding); hitTree->SetBranchAddress("trackID",&buffer.trackID); hitTree->SetBranchAddress("parentID",&buffer.parentID); + hitTree->SetBranchAddress("trackLocalTime",&buffer.trackLocalTime); hitTree->SetBranchAddress("time",&buffer.time); hitTree->SetBranchAddress("edep",&buffer.edep); hitTree->SetBranchAddress("stepLength",&buffer.stepLength); @@ -413,6 +422,11 @@ void GateHitTree::SetBranchAddresses(TTree* hitTree,GateRootHitBuffer& buffer) hitTree->SetBranchAddress("nCrystalCompton",&buffer.nCrystalCompton); hitTree->SetBranchAddress("nPhantomRayleigh",&buffer.nPhantomRayleigh); hitTree->SetBranchAddress("nCrystalRayleigh",&buffer.nCrystalRayleigh); + if (hitTree->GetBranch("nInteractions")) { + hitTree->SetBranchAddress("nInteractions",&buffer.nInteractions); + } else { + buffer.nInteractions = -1; + } hitTree->SetBranchAddress("primaryID",&buffer.primaryID); hitTree->SetBranchAddress("axialPos",&buffer.axialPos); @@ -424,6 +438,7 @@ void GateHitTree::SetBranchAddresses(TTree* hitTree,GateRootHitBuffer& buffer) hitTree->SetBranchAddress("sourceType",&buffer.sourceType); hitTree->SetBranchAddress("decayType",&buffer.decayType); hitTree->SetBranchAddress("gammaType",&buffer.gammaType); + hitTree->SetBranchAddress("decayIndex",&buffer.decayIndex); } else diff --git a/source/physics/include/GateEmittedGammaInformation.hh b/source/physics/include/GateEmittedGammaInformation.hh index 157b1d686..d731e9c70 100644 --- a/source/physics/include/GateEmittedGammaInformation.hh +++ b/source/physics/include/GateEmittedGammaInformation.hh @@ -30,7 +30,8 @@ class GateEmittedGammaInformation : public G4VUserPrimaryParticleInformation NotDefined = 0, // by default SingleGammaEmitter = 1, // just emitted single gamma with specfic energy ( by class GateGammaEmissionModel ) ParaPositronium = 2, // 2 gammas ( plus prompt if it is required ) from pPs decay ( by GatePositroniumDecayModel ) - OrthoPositronium = 3 // 3 gammas ( plus prompt if it is required ) from oPs decay ( by GatePositroniumDecayModel ) + OrthoPositronium = 3, // 3 gammas ( plus prompt if it is required ) from oPs decay ( by GatePositroniumDecayModel ) + DirectAnnihilation = 4 // 2 gammas from direct annihilation of positron without formation of positronium }; /** This enum specifies model of source decay ( if it is present ) @@ -61,6 +62,13 @@ class GateEmittedGammaInformation : public G4VUserPrimaryParticleInformation void SetGammaKind( GammaKind gamma_kind ); GammaKind GetGammaKind() const; + /** Set decay channel index + **/ + void SetDecayIndex( G4int decay_index ); + /** Get decay channel index + **/ + G4int GetDecayIndex() const; + /** Set polarization of gamma at the moment when it was emitted **/ void SetInitialPolarization( const G4ThreeVector& polarization ); @@ -80,6 +88,7 @@ class GateEmittedGammaInformation : public G4VUserPrimaryParticleInformation SourceKind fSourceKind = SourceKind::NotDefined; DecayModel fDecayModel = DecayModel::None; GammaKind fGammaKind = GammaKind::Unknown; + G4int fDecayIndex = -1; G4ThreeVector fInitialPolarization = G4ThreeVector( 0.0, 0.0, 0.0 ); G4double fTimeShift = 0.0;//[ns] }; diff --git a/source/physics/include/GateGammaEmissionModel.hh b/source/physics/include/GateGammaEmissionModel.hh index c747cc4d4..72c52f94b 100644 --- a/source/physics/include/GateGammaEmissionModel.hh +++ b/source/physics/include/GateGammaEmissionModel.hh @@ -8,15 +8,12 @@ #define GateGammaEmissionModel_hh #include "G4PrimaryParticle.hh" -#include -#include #include "GateEmittedGammaInformation.hh" #include "G4Event.hh" #include "G4SystemOfUnits.hh" /** Author: Mateusz Bała * Email: bala.mateusz@gmail.com - * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) * About class: Basic class for other model of gammas emission. Provides basic tools for calculation and can generate single gamma. **/ class GateGammaEmissionModel @@ -40,13 +37,6 @@ class GateGammaEmissionModel **/ G4double GetEmissionEnergy() const; - /** Set seed for generators from "Randomize.hh" - **/ - void SetSeed( G4long seed ); - /** Get seed for generators from "Randomize.hh" - **/ - G4long GetSeed() const; - /** Generate single vertex with single gamma **/ virtual G4int GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position); diff --git a/source/physics/include/GatePositronium.hh b/source/physics/include/GatePositronium.hh new file mode 100644 index 000000000..325a3addb --- /dev/null +++ b/source/physics/include/GatePositronium.hh @@ -0,0 +1,42 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About class: Represents a single positronium species (para-Ps or ortho-Ps), wrapping the Geant4 particle decay table to provide access to lifetime, annihilation gamma count, and decay products. + **/ + +#ifndef GatePositronium_hh +#define GatePositronium_hh + +#include "G4DecayProducts.hh" +#include "GatePositroniumDecayChannel.hh" + +class GatePositronium { +public: + GatePositronium(const G4String& name, G4double life_time); + ~GatePositronium() = default; + + GatePositronium(const GatePositronium&) = delete; + GatePositronium& operator=(const GatePositronium&) = delete; + GatePositronium(GatePositronium&&) = default; + GatePositronium& operator=(GatePositronium&&) = default; + + + G4double GetLifeTime() const; + const G4String& GetName() const; + G4int GetAnnihilationGammasNumber() const; + G4DecayProducts *GetDecayProducts() const; + +private: + G4String fName; + G4double fLifeTime = 0.0; //[ns] + GatePositroniumDecayChannel *pDecayChannel = nullptr; // Todo check who owns pDecayChannel? +}; +#endif diff --git a/source/physics/include/GatePositroniumConstants.hh b/source/physics/include/GatePositroniumConstants.hh new file mode 100644 index 000000000..c8071d50e --- /dev/null +++ b/source/physics/include/GatePositroniumConstants.hh @@ -0,0 +1,27 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About class: Namespace of physical constants for positronium physics: para-Ps lifetime (0.1244 ns), ortho-Ps mean lifetime (142 ns), para-to-ortho fraction (1/3), and hyperfine coefficient (372). + **/ + +#ifndef GatePositroniumConstants_hh +#define GatePositroniumConstants_hh + +namespace GatePositroniumConstants +{ + constexpr double kParaPsLifetime_ns = 0.1244; + constexpr double kParaToOrthoPsFraction = 1.0/3.0; + constexpr double kOrthoPsMeanLifetime_ns = 142.; + constexpr double kHyperfineCoefficient = 372; +} + +#endif + diff --git a/source/physics/include/GatePositroniumDecayChannel.hh b/source/physics/include/GatePositroniumDecayChannel.hh index 290970cf3..3af067371 100644 --- a/source/physics/include/GatePositroniumDecayChannel.hh +++ b/source/physics/include/GatePositroniumDecayChannel.hh @@ -7,16 +7,14 @@ #ifndef GatePositroniumDecayChannel_hh #define GatePositroniumDecayChannel_hh -#include "globals.hh" #include "G4GeneralPhaseSpaceDecay.hh" #include "G4PhysicalConstants.hh" -#include "G4SystemOfUnits.hh" /** Author: Mateusz Bała * Email: bala.mateusz@gmail.com - * Theorem author for oPs decay: Daria Kamińska ( Eur. Phys. J. C (2016) 76:445 ) - * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) - * About class: Implements decay of positronium ( pPs and oPs ). Provides support for polarization. + * Original author of the oPs decay model: Daria Kamińska et al. ( Eur. Phys. J. C (2016) 76:445 ) + * Refactored by: Wojciech Krzemien + * About class: Implements pPs and oPs positronium decays. Provides support for polarization. **/ class GatePositroniumDecayChannel : public G4GeneralPhaseSpaceDecay { @@ -25,11 +23,11 @@ class GatePositroniumDecayChannel : public G4GeneralPhaseSpaceDecay //Describes for which positronium we need decay enum PositroniumKind { NotDefined, ParaPositronium, OrthoPositronium }; - GatePositroniumDecayChannel( const G4String& theParentName, G4double theBR); - virtual ~GatePositroniumDecayChannel(); + GatePositroniumDecayChannel( const G4String& parentName, G4double BR); + virtual ~GatePositroniumDecayChannel() = default; /** Return gammas from positronium decay **/ - virtual G4DecayProducts* DecayIt(G4double) override; + virtual G4DecayProducts* DecayIt(G4double m=0.0) override; protected: /** Return gammas from para-positronium decay @@ -41,28 +39,31 @@ class GatePositroniumDecayChannel : public G4GeneralPhaseSpaceDecay /** Calculate cross section Mij matrix element * Based on "Quantum electrodynamics" V. B. BERESTETSKY. * Chapter: 89. Annihilation of positronium - * Exquantation: 89.14 + * Equation: 89.14 **/ - G4double GetOrthoPsM( const G4double w1, const G4double w2, const G4double w3 ) const; + G4double GetOrthoPsM(G4double w1, G4double w2, G4double w3) const; /** Calculate polarization orthogonal to momentum direction **/ - G4ThreeVector GetPolarization( const G4ThreeVector& momentum ) const; + G4ThreeVector GetPolarization(const G4ThreeVector& momentum) const; /** Generate perpendiculator vector ( to calculate orthogonal polarization ) **/ G4ThreeVector GetPerpendicularVector(const G4ThreeVector& v) const; protected: - //Decay constants - const G4String kParaPositroniumName = "pPs"; - const G4String kOrthoPositroniumName = "oPs"; - const G4String kDaughterName = "gamma"; - const G4int kParaPositroniumAnnihilationGammasNumber = 2; - const G4int kOrthoPositroniumAnnihilationGammasNumber = 3; - const G4double kPositroniumMass = 2.0 * electron_mass_c2; + static inline const G4String kParaPositroniumName = "pPs"; + static inline const G4String kOrthoPositroniumName = "oPs"; + static inline const G4String kDaughterName = "gamma"; + + //Decay paramters + static constexpr G4int kParaPositroniumAnnihilationGammasNumber = 2; + static constexpr G4int kOrthoPositroniumAnnihilationGammasNumber = 3; + + static constexpr G4double kPositroniumMass = 2.0 * electron_mass_c2; + ///This is maximal number which can be calculated by function GetOrthoPsM() - determined based on 10^7 iterations + static constexpr G4double kOrthoPsMMax = 7.65928; + static constexpr G4double kElectronMass = electron_mass_c2; //[MeV] + PositroniumKind fPositroniumKind = PositroniumKind::NotDefined; - ///This is maxiaml number which can be calculated by function GetOrthoPsM() - based on 10^7 iterations - const G4double kOrthoPsMMax = 7.65928; - const G4double kElectronMass = electron_mass_c2; //[MeV] }; #endif diff --git a/source/physics/include/GatePositroniumDecayModel.hh b/source/physics/include/GatePositroniumDecayModel.hh index b64c6b97c..df704e081 100644 --- a/source/physics/include/GatePositroniumDecayModel.hh +++ b/source/physics/include/GatePositroniumDecayModel.hh @@ -4,131 +4,46 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#ifndef PositroniumDecayModel_hh -#define PositroniumDecayModel_hh +#ifndef GatePositroniumDecayModel_hh +#define GatePositroniumDecayModel_hh #include -#include "G4DecayTable.hh" -#include "G4ParticleTable.hh" -#include "G4PrimaryParticle.hh" -#include "G4GeneralPhaseSpaceDecay.hh" -#include "G4DecayTable.hh" -#include "G4ParticleDefinition.hh" -#include "GateEmittedGammaInformation.hh" -#include "GateGammaEmissionModel.hh" -#include "G4PhysicalConstants.hh" -#include "G4SystemOfUnits.hh" - -#include "G4VDecayChannel.hh" -#include "G4ParticleDefinition.hh" +#include "G4PrimaryParticle.hh" #include "G4PrimaryVertex.hh" -/** Author: Mateusz Bała - * Email: bala.mateusz@gmail.com - * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) - * About class: Class generate gammas from positronium decay with including deexcitation gamma ( prompt gamma ) if is required. - **/ -class GatePositroniumDecayModel : public GateGammaEmissionModel -{ - public: - - /** About class: representation of positronium for main class. Provides access to positronium decay channel. - **/ - class Positronium - { - public: - Positronium( G4String name, G4double life_time, G4int annihilation_gammas_number ); - void SetLifeTime( const G4double& life_time ); - G4double GetLifeTime() const; - G4String GetName() const; - G4int GetAnnihilationGammasNumber() const; - G4DecayProducts* GetDecayProducts(); - private: - G4String fName = ""; - G4double fLifeTime = 0.0;//[ns] - G4int fAnnihilationGammasNumber = 0; - G4VDecayChannel* pDecayChannel = nullptr; - }; - - /** Positronium kind tells which positronium we will use. - * Depends on used Ps gammas number and time will be different. - **/ - enum PositroniumKind { pPs, oPs }; - /** Decay model descibes decay of positronium. - * For example for PositroniumKind::pPs we have decays: - * 1) Standard: pPs -> 2 gamma - * 2) WithPrompt: pPs* -> 2 gamma + prompt_gamma - * Only prompt gamma has nod modified time, other gammas always has modified time. - **/ - enum DecayModel { Standard, WithPrompt }; - - GatePositroniumDecayModel(); - virtual ~GatePositroniumDecayModel(); - - void SetPositroniumKind( PositroniumKind positronium_kind ); - PositroniumKind GetPositroniumKind() const; - - void SetDecayModel( DecayModel decay_model ); - DecayModel GetDecayModel() const; - - void SetPostroniumLifetime( G4String positronium_name, G4double life_time ); //[ns] +#include "GateGammaEmissionModel.hh" +#include "GatePositroniumDecayModelParams.hh" +#include "GatePositronium.hh" - void SetPromptGammaEnergy( G4double prompt_energy ); //[keV] - G4double GetPromptGammaEnergy() const; //[keV] +class GatePositroniumDecayModel:public GateGammaEmissionModel +{ + public: + static int getPositroniumDecayIndex(const std::vector& fractions); + static G4ThreeVector AddPositronRangeShift(const G4ThreeVector& original_position, G4double mean_positron_range); - /** Set probability of gammas emission from para-positronium decay ( and ortho-positronium too ) - * @param: fraction - number in range from 0.0 to 1.0 ( 0.0 - generate only gammas from oPs decay, 1.0 - generate only gammas from pPs decay ) - **/ - void SetParaPositroniumFraction( G4double fraction ); + explicit GatePositroniumDecayModel(const PositroniumDecayModelParams& modelParams); - /** Generate vertices and fill them with gammas. Vertices number depends on used DecayModel: - * - one ( if DecayModel::Standard ) - * - two ( if DecayModel::WithPrompt ) - **/ + protected: virtual G4int GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position) override; - - protected: - /** Provides additional information for user about gamma - **/ - virtual GateEmittedGammaInformation* GetPrimaryParticleInformation( const G4PrimaryParticle* pp, const GateEmittedGammaInformation::GammaKind& gamma_kind ) const override; - /** Depends on used model and setted fractions it chooses positronium which decay will be used to generate gammas - **/ - void PreparePositroniumParametrization(); - /** Generate vertex for deexcitation gamma ( prompt gamma ) - position and time is the same as generted by source - **/ - G4PrimaryVertex* GetPrimaryVertexFromDeexcitation(const G4double& particle_time, const G4ThreeVector& particle_position ); - /** Generate vertex for annihilation gammas - position is the same as generted by source, but time is shifted by positronium lifetime ( T0 + f(lifetime)) - **/ - G4PrimaryVertex* GetPrimaryVertexFromPositroniumAnnihilation( const G4double& particle_time, const G4ThreeVector& particle_position ); - /** Generate deexcitation ( prompt ) gamma - **/ - G4PrimaryParticle* GetGammaFromDeexcitation(); - /** Generate annihilation gammas - **/ - std::vector GetGammasFromPositroniumAnnihilation(); - - protected: - //Positronium model - for para-positronium - Positronium fParaPs = Positronium( "pPs", 0.1244 * ns, 2 ); - //Positronium model - for ortho-positronium - Positronium fOrthoPs = Positronium( "oPs", 138.6 * ns, 3 ); - //Positronium model - for current event - Positronium* pInfoPs = nullptr; - //Default deexcitation gamma energy - if user didn't set prompt gamma energy this value will be used - const G4double kSodium22DeexcitationGammaEnergy = 1.022 * MeV;//[MeV] - - //Which positronium use for gammas generation - PositroniumKind fPositroniumKind = PositroniumKind::pPs; - //Which decay model use for gammas generation - DecayModel fDecayModel = DecayModel::Standard; - //Dexcitation gamma energy - G4double fPromptGammaEnergy = kSodium22DeexcitationGammaEnergy;//[MeV] - //Propability of emiiting gammas from para-positronium ( number in range from 0.0 to 1.0 ) - G4double fParaPositroniumFraction = 1.0; - //It is required to generate mixed positronium decays ( pPs and oPs witch propability controled by varaible fParaPositroniumFraction ) - G4bool fUsePositroniumFractions = false; + G4PrimaryVertex* GetPrimaryVertexFromDeexcitation(G4double particle_time, const G4ThreeVector& particle_position, int decayIndex); + G4PrimaryVertex* GetPrimaryVertexFromPositroniumAnnihilation(G4double particle_time, const G4ThreeVector &particle_position, int decayIndex); + G4PrimaryParticle* GetGammaFromDeexcitation(int decayIndex); + std::vector GetGammasFromPositroniumAnnihilation(int decayIndex); + /* The decay model is determined by presence of prompt gamma - if it is present (probability > 0) then it is Deexcitation model, + * otherwise it is Standard model. + */ + GateEmittedGammaInformation::DecayModel GetDecayModel(const int decayIndex) const; + /* + * The source kind is determined by type of positron-electron interaction + * ( paraPs, orthoPs or direct annihilation without positronium formation ) + */ + GateEmittedGammaInformation::SourceKind GetSourceKind(int decayIndex) const; + + private: + PositroniumDecayModelParams fModelParams; + std::vector fPositroniumDecayChannel; }; #endif diff --git a/source/physics/include/GatePositroniumDecayModelParams.hh b/source/physics/include/GatePositroniumDecayModelParams.hh new file mode 100644 index 000000000..717ec8e00 --- /dev/null +++ b/source/physics/include/GatePositroniumDecayModelParams.hh @@ -0,0 +1,37 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About class: Data structure holding all per-component parameters for the positronium decay model: fractions, lifetimes, decay kinds (2γ/3γ), positron-electron interaction types, electron capture probabilities, prompt gamma parameters, and positron range settings. + **/ + +#ifndef GatePositroniumDecayModelParams_hh +#define GatePositroniumDecayModelParams_hh + +#include + +enum PositroniumDecayKind {k2Gamma, k3Gamma}; + +enum PositronElectronInteraction {kParaPs, kDirect, kOrthoPs}; + +struct PositroniumDecayModelParams +{ + std::vector fFractions; + std::vector fLifetimes; + std::vector fElectronCaptureProbabilities; + std::vector fPromptGammaProbabilities; + std::vector fPromptGammaEnergy; + std::vector fDecayKind; + std::vector fPositronInteractions; + std::vector fMeanPositronRangeEnabled; + std::vector fMeanPositronRange; +}; + +#endif diff --git a/source/physics/include/GatePositroniumDecayParamsGenerator.hh b/source/physics/include/GatePositroniumDecayParamsGenerator.hh new file mode 100644 index 000000000..77d92b2d5 --- /dev/null +++ b/source/physics/include/GatePositroniumDecayParamsGenerator.hh @@ -0,0 +1,61 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About class: Builder that assembles and validates a PositroniumDecayModelParams from individually set parameters; supports three preset modes (para-Ps only, ortho-Ps only, and mixed Ps with user-defined components including electron capture). + **/ + +#ifndef GatePositroniumDecayParamsGenerator_hh +#define GatePositroniumDecayParamsGenerator_hh + +#include +#include + +#include "GatePositroniumDecayModelParams.hh" + +/*! class GatePositroniumDecayParamsGenerator + * brief Brief class description + * + * Detailed description + */ + + +class GatePositroniumDecayParamsGenerator +{ +public: + + GatePositroniumDecayParamsGenerator()=default; + virtual ~GatePositroniumDecayParamsGenerator()=default; + + void SetElectronCaptureProbabilities(const std::vector& electronCaptureProb); + void SetPromptGammaProbabilities(const std::vector& promptGammaProb); + void SetPromptGammaEnergies(const std::vector& energies); + void SetPositroniumLifetimes(const std::vector& fPositroniumLifetimes); + void SetDecayKinds(const std::vector& decayKinds); + void SetPositronInteractions(const std::vector& positronInteractions); + void SetPositroniumFraction(const std::vector& positroniumFractions); + void SetMeanPositronRange(const std::vector& meanPositronRange); + + PositroniumDecayModelParams generatePositroniumDecayParams() const; + void validatePositroniumDecayParams(const PositroniumDecayModelParams& params) const; + + +private: + std::optional> fPositroniumFractions; + std::optional> fPositroniumLifetimes; + std::optional> fPromptGammaProbabilities; + std::optional> fPromptGammaEnergies; + std::optional> fElectronCaptureProbabilities; + std::optional> fDecayKinds; + std::optional> fPositronInteractions; + std::optional> fMeanPositronRangeEnabled; + std::optional> fMeanPositronRange; +}; +#endif diff --git a/source/physics/include/GatePositroniumHelper.hh b/source/physics/include/GatePositroniumHelper.hh new file mode 100644 index 000000000..658dd179c --- /dev/null +++ b/source/physics/include/GatePositroniumHelper.hh @@ -0,0 +1,34 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About class: Utility class computing per-component decay fractions from lifetimes and positron-electron interaction types (ortho-Ps, direct annihilation, para-Ps), using the hyperfine coefficient and mean ortho-Ps lifetime; also normalizes fraction vectors. + **/ + +#ifndef GatePositroniumHelper_h +#define GatePositroniumHelper_h 1 + +#include "GatePositroniumDecayModelParams.hh" + +class GatePositroniumHelper +{ +public: + GatePositroniumHelper() {}; + virtual ~GatePositroniumHelper() {}; + + struct PositroniumDecayModelParams CalculateFractionsFromLifetimes(struct PositroniumDecayModelParams params); + float CalcPPsFractionFromOPs(std::vector fractions, std::vector decays); + std::pair CalcFractionsFromLifetime(float intensity, float lifetime, PositronElectronInteraction inter); + float CalcFractionFromOPsLifetime(float intensity, float lifetime, PositroniumDecayKind decay); + float CalcFractionFromDirectLifetime(float intensity, PositroniumDecayKind decay); + std::vector NormalizeFractions (std::vector fractions); +}; + +#endif diff --git a/source/physics/include/GatePositroniumSource.hh b/source/physics/include/GatePositroniumSource.hh new file mode 100644 index 000000000..e4e76ccf9 --- /dev/null +++ b/source/physics/include/GatePositroniumSource.hh @@ -0,0 +1,40 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About class: Geant4 primary event source for positronium annihilation; lazily initializes a GatePositroniumDecayModel from messenger-supplied parameters on the first event and delegates primary vertex generation to the model. + **/ + +#ifndef GatePositroniumSource_hh +#define GatePositroniumSource_hh + +#include + +#include "GateVSource.hh" +#include "GatePositroniumSourceMessenger.hh" +#include "GateGammaEmissionModel.hh" + +class GatePositroniumSource : public GateVSource +{ +public: + explicit GatePositroniumSource(const G4String& name); + virtual ~GatePositroniumSource() = default; + + virtual G4int GeneratePrimaries( G4Event* event ) override; + + protected: + void PrepareModel(); + + protected: + std::unique_ptr pModel; + std::unique_ptr pMessenger; +}; + +#endif diff --git a/source/physics/include/GatePositroniumSourceMessenger.hh b/source/physics/include/GatePositroniumSourceMessenger.hh new file mode 100644 index 000000000..3c3e77b1b --- /dev/null +++ b/source/physics/include/GatePositroniumSourceMessenger.hh @@ -0,0 +1,66 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#ifndef GatePositroniumSourceMessenger_hh +#define GatePositroniumSourceMessenger_hh + +#include + +#include "GateVSourceMessenger.hh" +#include "G4UIcmdWithAnInteger.hh" +#include "G4UIcmdWithADoubleAndUnit.hh" +#include "G4UIcmdWithABool.hh" +#include "G4UIcmdWith3Vector.hh" +#include "G4UIcmdWithAString.hh" +#include "G4UIcmdWith3VectorAndUnit.hh" + +#include "GatePositroniumDecayParamsGenerator.hh" + +class GatePositroniumSource; + +/** + * About class: Messenger for GatePositroniumSource class + **/ +class GatePositroniumSourceMessenger: public GateVSourceMessenger +{ + public: + explicit GatePositroniumSourceMessenger(GatePositroniumSource *source); + ~GatePositroniumSourceMessenger()=default; + + void SetNewValue(G4UIcommand *command, G4String newValue) override; + + PositroniumDecayModelParams generatePositroniumDecayParams() const; + + protected: + void InitCommands(); + + G4UIcmdWithABool* GetBoolCmd(const G4String& cmd_name, const G4String& cmd_guidance); + G4UIcmdWithADoubleAndUnit* GetDoubleCmdWithUnit(const G4String& cmd_name, const G4String& cmd_guidance, const G4String& default_unit, const G4String& unit_candidates); + G4UIcmdWith3Vector* GetVectorCmd(const G4String& cmd_name, const G4String& cmd_guidance); + G4UIcmdWithAnInteger* GetIntCmd(const G4String& cmd_name, const G4String& cmd_guidance); + G4UIcmdWithAString* GetStringCmd(const G4String& cmd_name, const G4String& cmd_guidance); + G4UIcmdWith3VectorAndUnit* GetVectorCmdWithUnit(const G4String& cmd_name, const G4String& cmd_guidance, const G4String& default_unit, const G4String& unit_candidates); + + protected: + GatePositroniumSource* pSource = nullptr; + + std::unique_ptr upCmdSetPositroniumFractions; + std::unique_ptr upCmdSetPositroniumLifetimes; + std::unique_ptr upCmdSetPromptPhotonProbabilites; + std::unique_ptr upCmdSetPromptPhotonEnergies; + std::unique_ptr upCmdSetElectronCaptureProbabilities; + std::unique_ptr upCmdSetDecayKinds; + std::unique_ptr upCmdSetPositronInteractions; + + std::unique_ptr upCmdSetMeanPositronRange; + + std::vector parseListOfParamsWithUnit(const G4String& input) const; + + GatePositroniumDecayParamsGenerator fParamGenerator; + +}; + +#endif diff --git a/source/physics/include/legacy/GateEmittedGammaInformation.hh b/source/physics/include/legacy/GateEmittedGammaInformation.hh new file mode 100644 index 000000000..aad9a10e3 --- /dev/null +++ b/source/physics/include/legacy/GateEmittedGammaInformation.hh @@ -0,0 +1,90 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#ifndef Legacy_GateEmittedGammaInformation_hh +#define Legacy_GateEmittedGammaInformation_hh + +#include "G4VUserPrimaryParticleInformation.hh" +#include "G4ThreeVector.hh" + +namespace GateLegacy{ + +/** Author: Mateusz Bała + * Email: bala.mateusz@gmail.com + * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) + * About class: + **/ +class GateEmittedGammaInformation : public G4VUserPrimaryParticleInformation +{ + public: + GateEmittedGammaInformation(); + virtual ~GateEmittedGammaInformation(); + + /** This enum specifies what is source of gamma. + Each class which inherits from GateEmissionModel class should add + here its own kinds. + **/ + enum SourceKind + { + NotDefined = 0, // by default + SingleGammaEmitter = 1, // just emitted single gamma with specfic energy ( by class GateGammaEmissionModel ) + ParaPositronium = 2, // 2 gammas ( plus prompt if it is required ) from pPs decay ( by GatePositroniumDecayModel ) + OrthoPositronium = 3 // 3 gammas ( plus prompt if it is required ) from oPs decay ( by GatePositroniumDecayModel ) + }; + + /** This enum specifies model of source decay ( if it is present ) + **/ + enum DecayModel + { + None = 0, // Not specify - for example when there is not decay ( e.g. class GateGammaEmissionModel ) + Standard = 1, // Standard decay of source ( e.g. pPs -> 2 gamma ) + Deexcitation = 2 // If only deexcitation is present or when is part of standard decay changel ( Na22* -> Na22 + prompt gamma, Na22 -> Ne22 + e+, pPs-> 2 gamma <==> pPs* --> 2 gamma + prompt ) + }; + + /** This enum specifies of gamma emitted by source + **/ + enum GammaKind + { + Unknown = 0, // by default + Single = 1, // gamma is emitted alone ( by class GateGammaEmissionModel ) + Annihilation = 2, // gamma is a product of annihilation + Prompt = 3 // gamma is emitted from deexcitation process + }; + + void SetSourceKind( SourceKind source_kind ); + SourceKind GetSourceKind() const; + + void SetDecayModel( DecayModel decay_model ); + DecayModel GetDecayModel() const; + + void SetGammaKind( GammaKind gamma_kind ); + GammaKind GetGammaKind() const; + + /** Set polarization of gamma at the moment when it was emitted + **/ + void SetInitialPolarization( const G4ThreeVector& polarization ); + /** Get polarization of gamma at the moment when it was emitted + **/ + G4ThreeVector GetInitialPolarization() const; + /** Set time shift - caused by positronium lifetime + **/ + void SetTimeShift( const G4double& time_shift ); + /** Get time shift - caused by positronium lifetime + **/ + G4double GetTimeShift() const; + + virtual void Print() const; + + protected: + SourceKind fSourceKind = SourceKind::NotDefined; + DecayModel fDecayModel = DecayModel::None; + GammaKind fGammaKind = GammaKind::Unknown; + G4ThreeVector fInitialPolarization = G4ThreeVector( 0.0, 0.0, 0.0 ); + G4double fTimeShift = 0.0;//[ns] +}; +} + +#endif diff --git a/source/physics/include/GateExtendedVSource.hh b/source/physics/include/legacy/GateExtendedVSource.hh similarity index 95% rename from source/physics/include/GateExtendedVSource.hh rename to source/physics/include/legacy/GateExtendedVSource.hh index be4375824..6b2e8c9ad 100644 --- a/source/physics/include/GateExtendedVSource.hh +++ b/source/physics/include/legacy/GateExtendedVSource.hh @@ -4,12 +4,14 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#ifndef GateExtendedVSource_hh -#define GateExtendedVSource_hh +#ifndef Legacy_GateExtendedVSource_hh +#define Legacy_GateExtendedVSource_hh #include "GateVSource.hh" -#include "GateExtendedVSourceMessenger.hh" -#include "GateGammaEmissionModel.hh" +#include "legacy/GateExtendedVSourceMessenger.hh" +#include "legacy/GateGammaEmissionModel.hh" + +namespace GateLegacy{ /** About class: this is helper class to control if setting is in use **/ @@ -116,5 +118,6 @@ public: G4bool fBehaveLikeVSource = false; }; +} #endif diff --git a/source/physics/include/GateExtendedVSourceMessenger.hh b/source/physics/include/legacy/GateExtendedVSourceMessenger.hh similarity index 95% rename from source/physics/include/GateExtendedVSourceMessenger.hh rename to source/physics/include/legacy/GateExtendedVSourceMessenger.hh index e3cd9ffad..b5adc57e9 100644 --- a/source/physics/include/GateExtendedVSourceMessenger.hh +++ b/source/physics/include/legacy/GateExtendedVSourceMessenger.hh @@ -4,8 +4,8 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#ifndef GateExtendedVSourceMessenger_hh -#define GateExtendedVSourceMessenger_hh +#ifndef Legacy_GateExtendedVSourceMessenger_hh +#define Legacy_GateExtendedVSourceMessenger_hh #include "GateVSourceMessenger.hh" #include "G4UImessenger.hh" @@ -17,6 +17,9 @@ #include "G4UIcmdWith3VectorAndUnit.hh" #include + +namespace GateLegacy +{ class GateExtendedVSource; /** Author: Mateusz Bała @@ -55,5 +58,6 @@ class GateExtendedVSourceMessenger: public GateVSourceMessenger std::unique_ptr upCmdSetLifetime; }; +} #endif diff --git a/source/physics/include/legacy/GateGammaEmissionModel.hh b/source/physics/include/legacy/GateGammaEmissionModel.hh new file mode 100644 index 000000000..a2c4eee90 --- /dev/null +++ b/source/physics/include/legacy/GateGammaEmissionModel.hh @@ -0,0 +1,90 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#ifndef Legacy_GateGammaEmissionModel_hh +#define Legacy_GateGammaEmissionModel_hh + +#include "G4PrimaryParticle.hh" +#include +#include +#include "legacy/GateEmittedGammaInformation.hh" +#include "G4Event.hh" +#include "G4SystemOfUnits.hh" + +namespace GateLegacy{ + +/** Author: Mateusz Bała + * Email: bala.mateusz@gmail.com + * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) + * About class: Basic class for other model of gammas emission. Provides basic tools for calculation and can generate single gamma. + **/ +class GateGammaEmissionModel +{ + public: + + GateGammaEmissionModel(); + virtual ~GateGammaEmissionModel(); + + /** Force gamma emissions in only one direction. When this method is called fUseFixedEmissionDirection flag is switched to TRUE. + **/ + void SetFixedEmissionDirection( const G4ThreeVector& momentum_direction ); + /** Force gamma emissions in only one direction. When this method is fUseFixedEmissionDirection flag is switched to value=enable. + **/ + void SetEnableFixedEmissionDirection( const G4bool enable ); + + /** Set single gamma kinematic energy. + **/ + void SetEmissionEnergy( const G4double& energy ); + /** Get single gamma kinematic energy. + **/ + G4double GetEmissionEnergy() const; + + /** Set seed for generators from "Randomize.hh" + **/ + void SetSeed( G4long seed ); + /** Get seed for generators from "Randomize.hh" + **/ + G4long GetSeed() const; + + /** Generate single vertex with single gamma + **/ + virtual G4int GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position); + + protected: + /** Provides additional information for user about gamma + **/ + virtual GateEmittedGammaInformation* GetPrimaryParticleInformation( const G4PrimaryParticle* pp, const GateEmittedGammaInformation::GammaKind& gamma_kind ) const; + /** Generate random unit vector from uniform sphere distribution + **/ + G4ThreeVector GetUniformOnSphere() const; + /** Generate gamma polarization + **/ + G4ThreeVector GetPolarization( const G4ThreeVector& momentum ) const; + /** Generate perpendicular vector ( required for polarization ) + **/ + G4ThreeVector GetPerpendicularVector(const G4ThreeVector& v) const; + /** Set model name ( class name ) + **/ + void SetModelName( const G4String model_name ); + /** Report error + **/ + void NoticeError( G4String method_name, G4String exception_description ) const; + /** Generate single gamma with desired energy + **/ + G4PrimaryParticle* GetSingleGamma( const G4double& energy ) const; + + protected: + G4ThreeVector fFixedEmissionDirection = G4ThreeVector( 0.0, 0.0, 0.0 ); + G4bool fUseFixedEmissionDirection = false; + G4double fEmissionEnergy = 0.511 * MeV;//[MeV] + G4String fModelName = "GateGammaEmissionModel"; + //Gamma definiotion - the same is used in each one model + G4ParticleDefinition* pGammaDefinition = nullptr; + +}; +} + +#endif diff --git a/source/physics/include/GateOrthoPositronium.hh b/source/physics/include/legacy/GateOrthoPositronium.hh similarity index 88% rename from source/physics/include/GateOrthoPositronium.hh rename to source/physics/include/legacy/GateOrthoPositronium.hh index 73b224c0c..47a742b96 100644 --- a/source/physics/include/GateOrthoPositronium.hh +++ b/source/physics/include/legacy/GateOrthoPositronium.hh @@ -4,11 +4,13 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#ifndef GateOrthoPositronium_hh -#define GateOrthoPositronium_hh +#ifndef Legacy_GateOrthoPositronium_hh +#define Legacy_GateOrthoPositronium_hh #include "G4ParticleDefinition.hh" +namespace GateLegacy{ + /** Author: Mateusz Bała * Email: bala.mateusz@gmail.com * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) @@ -25,5 +27,6 @@ class GateOrthoPositronium : public G4ParticleDefinition static GateOrthoPositronium* OrthoPositroniumDefinition(); static GateOrthoPositronium* OrthoPositronium(); }; +} #endif diff --git a/source/physics/include/GateParaPositronium.hh b/source/physics/include/legacy/GateParaPositronium.hh similarity index 88% rename from source/physics/include/GateParaPositronium.hh rename to source/physics/include/legacy/GateParaPositronium.hh index 56e14c771..1efddd5fd 100644 --- a/source/physics/include/GateParaPositronium.hh +++ b/source/physics/include/legacy/GateParaPositronium.hh @@ -4,11 +4,13 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#ifndef GateParaPositronium_hh -#define GateParaPositronium_hh +#ifndef Legacy_GateParaPositronium_hh +#define Legacy_GateParaPositronium_hh #include "G4ParticleDefinition.hh" +namespace GateLegacy{ + /** Author: Mateusz Bała * Email: bala.mateusz@gmail.com * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) @@ -25,5 +27,6 @@ class GateParaPositronium : public G4ParticleDefinition static GateParaPositronium* ParaPositroniumDefinition(); static GateParaPositronium* ParaPositronium(); }; +} #endif diff --git a/source/physics/include/legacy/GatePositroniumDecayChannel.hh b/source/physics/include/legacy/GatePositroniumDecayChannel.hh new file mode 100644 index 000000000..859a22d52 --- /dev/null +++ b/source/physics/include/legacy/GatePositroniumDecayChannel.hh @@ -0,0 +1,70 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#ifndef Legacy_GatePositroniumDecayChannel_hh +#define Legacy_GatePositroniumDecayChannel_hh + +#include "globals.hh" +#include "G4GeneralPhaseSpaceDecay.hh" +#include "G4PhysicalConstants.hh" +#include "G4SystemOfUnits.hh" + +namespace GateLegacy{ + +/** Author: Mateusz Bała + * Email: bala.mateusz@gmail.com + * Theorem author for oPs decay: Daria Kamińska ( Eur. Phys. J. C (2016) 76:445 ) + * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) + * About class: Implements decay of positronium ( pPs and oPs ). Provides support for polarization. + **/ +class GatePositroniumDecayChannel : public G4GeneralPhaseSpaceDecay +{ + public: + + //Describes for which positronium we need decay + enum PositroniumKind { NotDefined, ParaPositronium, OrthoPositronium }; + + GatePositroniumDecayChannel( const G4String& theParentName, G4double theBR); + virtual ~GatePositroniumDecayChannel(); + /** Return gammas from positronium decay + **/ + virtual G4DecayProducts* DecayIt(G4double) override; + + protected: + /** Return gammas from para-positronium decay + **/ + G4DecayProducts* DecayParaPositronium(); + /** Return gammas from ortho-positronium decay + **/ + G4DecayProducts* DecayOrthoPositronium(); + /** Calculate cross section Mij matrix element + * Based on "Quantum electrodynamics" V. B. BERESTETSKY. + * Chapter: 89. Annihilation of positronium + * Exquantation: 89.14 + **/ + G4double GetOrthoPsM( const G4double w1, const G4double w2, const G4double w3 ) const; + /** Calculate polarization orthogonal to momentum direction + **/ + G4ThreeVector GetPolarization( const G4ThreeVector& momentum ) const; + /** Generate perpendiculator vector ( to calculate orthogonal polarization ) + **/ + G4ThreeVector GetPerpendicularVector(const G4ThreeVector& v) const; + + protected: + //Decay constants + const G4String kParaPositroniumName = "pPs"; + const G4String kOrthoPositroniumName = "oPs"; + const G4String kDaughterName = "gamma"; + const G4int kParaPositroniumAnnihilationGammasNumber = 2; + const G4int kOrthoPositroniumAnnihilationGammasNumber = 3; + const G4double kPositroniumMass = 2.0 * electron_mass_c2; + PositroniumKind fPositroniumKind = PositroniumKind::NotDefined; + ///This is maxiaml number which can be calculated by function GetOrthoPsM() - based on 10^7 iterations + const G4double kOrthoPsMMax = 7.65928; + const G4double kElectronMass = electron_mass_c2; //[MeV] +}; +} +#endif diff --git a/source/physics/include/legacy/GatePositroniumDecayModel.hh b/source/physics/include/legacy/GatePositroniumDecayModel.hh new file mode 100644 index 000000000..df7db0a37 --- /dev/null +++ b/source/physics/include/legacy/GatePositroniumDecayModel.hh @@ -0,0 +1,136 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#ifndef Legacy_PositroniumDecayModel_hh +#define Legacy_PositroniumDecayModel_hh + +#include +#include "G4DecayTable.hh" +#include "G4ParticleTable.hh" +#include "G4PrimaryParticle.hh" +#include "G4GeneralPhaseSpaceDecay.hh" +#include "G4DecayTable.hh" +#include "G4ParticleDefinition.hh" +#include "legacy/GateEmittedGammaInformation.hh" +#include "legacy/GateGammaEmissionModel.hh" + +#include "G4PhysicalConstants.hh" +#include "G4SystemOfUnits.hh" + +#include "G4VDecayChannel.hh" +#include "G4ParticleDefinition.hh" +#include "G4PrimaryVertex.hh" + +namespace GateLegacy { + +/** Author: Mateusz Bała + * Email: bala.mateusz@gmail.com + * Organization: J-PET (http://koza.if.uj.edu.pl/pet/) + * About class: Class generate gammas from positronium decay with including deexcitation gamma ( prompt gamma ) if is required. + **/ +class GatePositroniumDecayModel : public GateGammaEmissionModel +{ + public: + + /** About class: representation of positronium for main class. Provides access to positronium decay channel. + **/ + class Positronium + { + public: + Positronium( G4String name, G4double life_time, G4int annihilation_gammas_number ); + void SetLifeTime( const G4double& life_time ); + G4double GetLifeTime() const; + G4String GetName() const; + G4int GetAnnihilationGammasNumber() const; + G4DecayProducts* GetDecayProducts(); + private: + G4String fName = ""; + G4double fLifeTime = 0.0;//[ns] + G4int fAnnihilationGammasNumber = 0; + G4VDecayChannel* pDecayChannel = nullptr; + }; + + /** Positronium kind tells which positronium we will use. + * Depends on used Ps gammas number and time will be different. + **/ + enum PositroniumKind { pPs, oPs }; + /** Decay model descibes decay of positronium. + * For example for PositroniumKind::pPs we have decays: + * 1) Standard: pPs -> 2 gamma + * 2) WithPrompt: pPs* -> 2 gamma + prompt_gamma + * Only prompt gamma has nod modified time, other gammas always has modified time. + **/ + enum DecayModel { Standard, WithPrompt }; + + GatePositroniumDecayModel(); + virtual ~GatePositroniumDecayModel(); + + void SetPositroniumKind( PositroniumKind positronium_kind ); + PositroniumKind GetPositroniumKind() const; + + void SetDecayModel( DecayModel decay_model ); + DecayModel GetDecayModel() const; + + void SetPostroniumLifetime( G4String positronium_name, G4double life_time ); //[ns] + + void SetPromptGammaEnergy( G4double prompt_energy ); //[keV] + G4double GetPromptGammaEnergy() const; //[keV] + + /** Set probability of gammas emission from para-positronium decay ( and ortho-positronium too ) + * @param: fraction - number in range from 0.0 to 1.0 ( 0.0 - generate only gammas from oPs decay, 1.0 - generate only gammas from pPs decay ) + **/ + void SetParaPositroniumFraction( G4double fraction ); + + /** Generate vertices and fill them with gammas. Vertices number depends on used DecayModel: + * - one ( if DecayModel::Standard ) + * - two ( if DecayModel::WithPrompt ) + **/ + virtual G4int GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position) override; + + protected: + /** Provides additional information for user about gamma + **/ + virtual GateEmittedGammaInformation* GetPrimaryParticleInformation( const G4PrimaryParticle* pp, const GateEmittedGammaInformation::GammaKind& gamma_kind ) const override; + /** Depends on used model and setted fractions it chooses positronium which decay will be used to generate gammas + **/ + void PreparePositroniumParametrization(); + /** Generate vertex for deexcitation gamma ( prompt gamma ) - position and time is the same as generted by source + **/ + G4PrimaryVertex* GetPrimaryVertexFromDeexcitation(const G4double& particle_time, const G4ThreeVector& particle_position ); + /** Generate vertex for annihilation gammas - position is the same as generted by source, but time is shifted by positronium lifetime ( T0 + f(lifetime)) + **/ + G4PrimaryVertex* GetPrimaryVertexFromPositroniumAnnihilation( const G4double& particle_time, const G4ThreeVector& particle_position ); + /** Generate deexcitation ( prompt ) gamma + **/ + G4PrimaryParticle* GetGammaFromDeexcitation(); + /** Generate annihilation gammas + **/ + std::vector GetGammasFromPositroniumAnnihilation(); + + protected: + //Positronium model - for para-positronium + Positronium fParaPs = Positronium( "pPs", 0.1244 * ns, 2 ); + //Positronium model - for ortho-positronium + Positronium fOrthoPs = Positronium( "oPs", 138.6 * ns, 3 ); + //Positronium model - for current event + Positronium* pInfoPs = nullptr; + //Default deexcitation gamma energy - if user didn't set prompt gamma energy this value will be used + const G4double kSodium22DeexcitationGammaEnergy = 1.022 * MeV;//[MeV] + + //Which positronium use for gammas generation + PositroniumKind fPositroniumKind = PositroniumKind::pPs; + //Which decay model use for gammas generation + DecayModel fDecayModel = DecayModel::Standard; + //Dexcitation gamma energy + G4double fPromptGammaEnergy = kSodium22DeexcitationGammaEnergy;//[MeV] + //Propability of emiiting gammas from para-positronium ( number in range from 0.0 to 1.0 ) + G4double fParaPositroniumFraction = 1.0; + //It is required to generate mixed positronium decays ( pPs and oPs witch propability controled by varaible fParaPositroniumFraction ) + G4bool fUsePositroniumFractions = false; + +}; +} +#endif diff --git a/source/physics/src/GateEmittedGammaInformation.cc b/source/physics/src/GateEmittedGammaInformation.cc index 651dd2b37..40f37c963 100644 --- a/source/physics/src/GateEmittedGammaInformation.cc +++ b/source/physics/src/GateEmittedGammaInformation.cc @@ -22,6 +22,10 @@ void GateEmittedGammaInformation::SetGammaKind( GateEmittedGammaInformation::Gam GateEmittedGammaInformation::GammaKind GateEmittedGammaInformation::GetGammaKind() const { return fGammaKind; } +void GateEmittedGammaInformation::SetDecayIndex( G4int decay_index ) { fDecayIndex = decay_index; } + +G4int GateEmittedGammaInformation::GetDecayIndex() const { return fDecayIndex; } + void GateEmittedGammaInformation::SetInitialPolarization( const G4ThreeVector& polarization ) { fInitialPolarization = polarization; } G4ThreeVector GateEmittedGammaInformation::GetInitialPolarization() const { return fInitialPolarization; } diff --git a/source/physics/src/GateGammaEmissionModel.cc b/source/physics/src/GateGammaEmissionModel.cc index 88a7bd4e2..2e27a6e15 100644 --- a/source/physics/src/GateGammaEmissionModel.cc +++ b/source/physics/src/GateGammaEmissionModel.cc @@ -70,14 +70,6 @@ void GateGammaEmissionModel::SetEmissionEnergy( const G4double& energy ) G4double GateGammaEmissionModel::GetEmissionEnergy() const { return fEmissionEnergy; } -void GateGammaEmissionModel::SetSeed( G4long seed ) -{ - if ( seed < 0 ) { NoticeError( G4String( __FUNCTION__ ), "seed should be positive value." ); } - G4Random::setTheSeed( seed ); -} - -G4long GateGammaEmissionModel::GetSeed() const { return G4Random::getTheSeed (); } - G4ThreeVector GateGammaEmissionModel::GetUniformOnSphere() const { //Based on TRandom::Sphere diff --git a/source/physics/src/GatePhysicsList.cc b/source/physics/src/GatePhysicsList.cc index 8c89e0943..e3f4e76c5 100644 --- a/source/physics/src/GatePhysicsList.cc +++ b/source/physics/src/GatePhysicsList.cc @@ -72,8 +72,8 @@ #include "G4OpticalPhoton.hh" #include "G4OpticalPhysics.hh" -#include "GateParaPositronium.hh" -#include "GateOrthoPositronium.hh" +#include "legacy/GateParaPositronium.hh" +#include "legacy/GateOrthoPositronium.hh" //----------------------------------------------------------------------------------------- @@ -498,10 +498,10 @@ void GatePhysicsList::ConstructParticle() emDNAActivator->ConstructParticle(); //Construct positroniums - GateParaPositronium::ParaPositroniumDefinition(); - GateOrthoPositronium::OrthoPositroniumDefinition(); G4QuasiOpticalPhoton::QuasiOpticalPhotonDefinition(); + GateLegacy::GateParaPositronium::ParaPositroniumDefinition(); + GateLegacy::GateOrthoPositronium::OrthoPositroniumDefinition(); } //----------------------------------------------------------------------------------------- diff --git a/source/physics/src/GatePositronium.cc b/source/physics/src/GatePositronium.cc new file mode 100644 index 000000000..1f977cd16 --- /dev/null +++ b/source/physics/src/GatePositronium.cc @@ -0,0 +1,28 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include "G4DecayProducts.hh" +#include "G4ParticleTable.hh" +#include "G4DecayTable.hh" +#include "G4ParticleDefinition.hh" + +#include "GatePositronium.hh" + +GatePositronium::GatePositronium(const G4String& name, G4double life_time): fName(name), fLifeTime(life_time) +{ + G4ParticleDefinition *positronium_def = G4ParticleTable::GetParticleTable()->FindParticle(name); + G4DecayTable *positronium_decay_table = positronium_def->GetDecayTable(); + pDecayChannel = static_cast(positronium_decay_table->GetDecayChannel(0)); +} + +G4double GatePositronium::GetLifeTime() const { return fLifeTime; } + +const G4String& GatePositronium::GetName() const { return fName; } + +G4int GatePositronium::GetAnnihilationGammasNumber() const { return pDecayChannel->GetNumberOfDaughters(); } + +G4DecayProducts* GatePositronium::GetDecayProducts() const { return pDecayChannel->DecayIt(); } diff --git a/source/physics/src/GatePositroniumDecayChannel.cc b/source/physics/src/GatePositroniumDecayChannel.cc index 70b5bc282..29eb8f0af 100644 --- a/source/physics/src/GatePositroniumDecayChannel.cc +++ b/source/physics/src/GatePositroniumDecayChannel.cc @@ -10,55 +10,51 @@ #include "Randomize.hh" #include "G4LorentzVector.hh" -GatePositroniumDecayChannel::GatePositroniumDecayChannel(const G4String& theParentName, G4double theBR) +GatePositroniumDecayChannel::GatePositroniumDecayChannel(const G4String& parentName, G4double BR) { - G4int daughters_number = 0; - - if ( theParentName == kParaPositroniumName ) - { - daughters_number = kParaPositroniumAnnihilationGammasNumber; - fPositroniumKind = GatePositroniumDecayChannel::PositroniumKind::ParaPositronium; - } - else if ( theParentName == kOrthoPositroniumName ) - { - daughters_number = kOrthoPositroniumAnnihilationGammasNumber; - fPositroniumKind = GatePositroniumDecayChannel::PositroniumKind::OrthoPositronium; - } - else - { - #ifdef G4VERBOSE - if (GetVerboseLevel()>0) - { - G4cout << "GatePositroniumDecayChannel:: constructor :"; - G4cout << " parent particle is not positronium (pPs,oPs) but "; - G4cout << theParentName << G4endl; + G4int daughters_number = 0; + + if (parentName == kParaPositroniumName) { + daughters_number = kParaPositroniumAnnihilationGammasNumber; + fPositroniumKind = + GatePositroniumDecayChannel::PositroniumKind::ParaPositronium; + } else if (parentName == kOrthoPositroniumName) { + daughters_number = kOrthoPositroniumAnnihilationGammasNumber; + fPositroniumKind = + GatePositroniumDecayChannel::PositroniumKind::OrthoPositronium; + } else { +#ifdef G4VERBOSE + if (GetVerboseLevel() > 0) { + G4cout << "GatePositroniumDecayChannel:: constructor :"; + G4cout << " parent particle is not positronium (pPs,oPs) but "; + G4cout << parentName << G4endl; + } +#endif } - #endif - } - SetParentMass( kPositroniumMass ); - SetBR( theBR ); - SetParent( theParentName ); - SetNumberOfDaughters( daughters_number ); - for ( G4int daughter_index = 0; daughter_index < daughters_number; ++daughter_index ) { SetDaughter( daughter_index, kDaughterName ); } + SetParentMass(kPositroniumMass); + SetBR(BR); + SetParent(parentName); + SetNumberOfDaughters(daughters_number); + for (G4int daughter_index = 0; daughter_index < daughters_number; + ++daughter_index) { + SetDaughter(daughter_index, kDaughterName); + } } -GatePositroniumDecayChannel::~GatePositroniumDecayChannel() {} - G4DecayProducts* GatePositroniumDecayChannel::DecayIt(G4double) { - switch ( fPositroniumKind ) - { - case GatePositroniumDecayChannel::PositroniumKind::ParaPositronium: - return DecayParaPositronium(); - case GatePositroniumDecayChannel::PositroniumKind::OrthoPositronium: - return DecayOrthoPositronium(); - default: - G4cout << "GatePositroniumDecayChannel::DecayIt "; - G4cout << *parent_name << " can not decay " << G4endl; - DumpInfo(); - return nullptr; - }; + switch (fPositroniumKind) { + case GatePositroniumDecayChannel::PositroniumKind::ParaPositronium: + return DecayParaPositronium(); + case GatePositroniumDecayChannel::PositroniumKind::OrthoPositronium: + return DecayOrthoPositronium(); + default: + G4cout << "GatePositroniumDecayChannel::DecayIt "; + G4cout << *parent_name << " can not decay " << G4endl; + DumpInfo(); + return nullptr; + } } G4DecayProducts* GatePositroniumDecayChannel::DecayParaPositronium() @@ -89,7 +85,10 @@ G4DecayProducts* GatePositroniumDecayChannel::DecayOrthoPositronium() do { - if ( decay_products != nullptr ) { delete decay_products; } + if (decay_products) { + delete decay_products; + decay_products = nullptr; + } decay_products = G4GeneralPhaseSpaceDecay::DecayIt(); @@ -104,39 +103,35 @@ G4DecayProducts* GatePositroniumDecayChannel::DecayOrthoPositronium() weight = GetOrthoPsM( lv_gamma_1.e(), lv_gamma_2.e(), lv_gamma_3.e() ); random_weight = kOrthoPsMMax * G4UniformRand(); } - while( random_weight > weight ); + while (random_weight > weight); ///Polarization - G4ThreeVector polarization_gamma_1 = GetPolarization( gamma_1->GetMomentumDirection() ); - G4ThreeVector polarization_gamma_2 = GetPolarization( gamma_2->GetMomentumDirection() ); - G4ThreeVector polarization_gamma_3 = GetPolarization( gamma_3->GetMomentumDirection() ); - - gamma_1->SetPolarization( polarization_gamma_1.x(), polarization_gamma_1.y(), polarization_gamma_1.z() ); - gamma_2->SetPolarization( polarization_gamma_2.x(), polarization_gamma_2.y(), polarization_gamma_2.z() ); - gamma_3->SetPolarization( polarization_gamma_3.x(), polarization_gamma_3.y(), polarization_gamma_3.z() ); - + for (size_t i = 0; i < decay_products->entries(); i++) + { + auto* gamma = (*decay_products)[i]; + auto polarization = GetPolarization(gamma->GetMomentumDirection()); + gamma->SetPolarization(polarization.x(), polarization.y(), polarization.z()); + } return decay_products; } -G4double GatePositroniumDecayChannel::GetOrthoPsM( const G4double w1, const G4double w2, const G4double w3 ) const +G4double GatePositroniumDecayChannel::GetOrthoPsM(G4double w1, G4double w2, G4double w3) const { - return pow( ( kElectronMass - w1 ) / ( w2 * w3 ), 2 ) + pow( ( kElectronMass - w2 ) / ( w1 * w3 ), 2 ) + pow( ( kElectronMass - w3 ) / ( w1 * w2 ), 2 ); + return pow((kElectronMass - w1) / (w2 * w3), 2) + + pow((kElectronMass - w2) / (w1 * w3), 2) + + pow((kElectronMass - w3) / (w1 * w2), 2); } -G4ThreeVector GatePositroniumDecayChannel::GetPolarization( const G4ThreeVector& momentum ) const +G4ThreeVector GatePositroniumDecayChannel::GetPolarization(const G4ThreeVector& momentum) const { - G4ThreeVector polarization(0.0,0.0,0.0); - - G4ThreeVector a0,b0,d0; - d0 = momentum.unit(); - a0 = GetPerpendicularVector( d0 ).unit(); - b0 = d0.cross( a0 ).unit(); - G4double angle_radians = G4UniformRand() * M_PI; - polarization = std::cos( angle_radians ) * a0 + std::sin( angle_radians ) * b0; - polarization.unit(); - return polarization; + auto d0 = momentum.unit(); + auto a0 = GetPerpendicularVector(d0).unit(); + auto b0 = d0.cross(a0).unit(); + G4double angle_radians = G4UniformRand() * M_PI; + auto polarization = std::cos(angle_radians) * a0 + std::sin(angle_radians) * b0; + return polarization.unit(); } - + G4ThreeVector GatePositroniumDecayChannel::GetPerpendicularVector(const G4ThreeVector& v) const { G4double dx = v.x(); diff --git a/source/physics/src/GatePositroniumDecayModel.cc b/source/physics/src/GatePositroniumDecayModel.cc index a30e64a4e..410a91b61 100644 --- a/source/physics/src/GatePositroniumDecayModel.cc +++ b/source/physics/src/GatePositroniumDecayModel.cc @@ -4,168 +4,144 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#include "GatePositroniumDecayModel.hh" -#include "Randomize.hh" #include #include -#include +#include + +#include "Randomize.hh" #include "G4DecayProducts.hh" #include "G4LorentzVector.hh" -#include "G4ParticleTable.hh" - -GatePositroniumDecayModel::Positronium::Positronium( G4String name, G4double life_time, G4int annihilation_gammas_number ) : fName( name ), fLifeTime( life_time ), fAnnihilationGammasNumber( annihilation_gammas_number ) -{ - G4ParticleDefinition* positronium_def = G4ParticleTable::GetParticleTable()->FindParticle( name ); - G4DecayTable* positronium_decay_table = positronium_def->GetDecayTable(); - pDecayChannel = positronium_decay_table->GetDecayChannel(0); -} -void GatePositroniumDecayModel::Positronium::SetLifeTime( const G4double& life_time ) { fLifeTime = life_time; } - -G4double GatePositroniumDecayModel::Positronium::GetLifeTime() const { return fLifeTime; } - -G4String GatePositroniumDecayModel::Positronium::GetName() const { return fName; } - -G4int GatePositroniumDecayModel::Positronium::GetAnnihilationGammasNumber() const { return fAnnihilationGammasNumber; } - -G4DecayProducts* GatePositroniumDecayModel::Positronium::GetDecayProducts() { return pDecayChannel->DecayIt(); } +#include "GatePositroniumDecayModel.hh" -GatePositroniumDecayModel::GatePositroniumDecayModel() -{ - SetModelName( "GatePositroniumDecayModel" ); +int GatePositroniumDecayModel::getPositroniumDecayIndex(const std::vector& fractions) { + auto r = G4UniformRand(); + float curr_frac_cumulative = 0.0; + for (int i = 0; i < fractions.size(); ++i) { + curr_frac_cumulative = curr_frac_cumulative + fractions[i]; + if(r<= curr_frac_cumulative) return i; + } + return static_cast(fractions.size()) - 1; } -GatePositroniumDecayModel::~GatePositroniumDecayModel() {} - -void GatePositroniumDecayModel::SetPositroniumKind( GatePositroniumDecayModel::PositroniumKind positronium_kind ) { fPositroniumKind = positronium_kind; } - -GatePositroniumDecayModel::PositroniumKind GatePositroniumDecayModel::GetPositroniumKind() const { return fPositroniumKind; } - -void GatePositroniumDecayModel::SetDecayModel( GatePositroniumDecayModel::DecayModel decay_model ) { fDecayModel = decay_model; } - -GatePositroniumDecayModel::DecayModel GatePositroniumDecayModel::GetDecayModel() const { return fDecayModel; } - -void GatePositroniumDecayModel::SetPostroniumLifetime( G4String positronium_name, G4double life_time ) -{ - if ( !( life_time > 0.0 ) ) { NoticeError( G4String( __FUNCTION__ ), "positronium life-time should be positive value." ); } - - if ( positronium_name == fParaPs.GetName() ) { fParaPs.SetLifeTime( life_time ); } - else if ( positronium_name == fOrthoPs.GetName() ) { fOrthoPs.SetLifeTime( life_time ); } - else { NoticeError( G4String( __FUNCTION__ ), "Unknown positronium name." ); } +G4ThreeVector GatePositroniumDecayModel::AddPositronRangeShift(const G4ThreeVector& original_position, G4double mean_positron_range) +{ + // r = sqrt(x**2+y**2+z**2) + // = sigma * sqrt(8/Pi) // matching mean for 3-D Gaussian + const G4double sqrt8_over_pi = std::sqrt(8.0/CLHEP::pi); + G4double sigma = mean_positron_range/sqrt8_over_pi ; + G4ThreeVector shift(G4RandGauss::shoot(0., sigma), + G4RandGauss::shoot(0., sigma), + G4RandGauss::shoot(0., sigma)); + return original_position + shift; } -void GatePositroniumDecayModel::SetPromptGammaEnergy( G4double prompt_energy ) -{ - if ( !( prompt_energy > 0.0 ) ) { NoticeError( G4String( __FUNCTION__ ), "prompt gamma energy should be positive value." ); } - fPromptGammaEnergy = prompt_energy; +GatePositroniumDecayModel::GatePositroniumDecayModel(const PositroniumDecayModelParams& modelParams):fModelParams(modelParams) +{ + auto num_of_decay_channels = fModelParams.fDecayKind.size(); + for (int i = 0; i < num_of_decay_channels; i++) { + if (fModelParams.fDecayKind[i] == PositroniumDecayKind::k2Gamma) + { + fPositroniumDecayChannel.push_back(std::move(GatePositronium("pPs", fModelParams.fLifetimes[i]* ns))); + } else { + fPositroniumDecayChannel.push_back(std::move(GatePositronium("oPs", fModelParams.fLifetimes[i]* ns))); + } + } } -G4double GatePositroniumDecayModel::GetPromptGammaEnergy() const { return fPromptGammaEnergy; } -void GatePositroniumDecayModel::SetParaPositroniumFraction( G4double fraction ) +G4PrimaryVertex* GatePositroniumDecayModel::GetPrimaryVertexFromPositroniumAnnihilation(G4double particle_time, const G4ThreeVector& particle_position, int decayIndex) { - fUsePositroniumFractions = true; - fParaPositroniumFraction = fraction; -} + bool is_positron_range_enabled = fModelParams.fMeanPositronRangeEnabled[decayIndex]; -void GatePositroniumDecayModel::PreparePositroniumParametrization() -{ - if ( pInfoPs != nullptr ) + G4double shifted_particle_time = particle_time + G4RandExponential::shoot(fModelParams.fLifetimes[decayIndex]); + + auto shifted_particle_position = particle_position; + if (is_positron_range_enabled) { - if ( !fUsePositroniumFractions ) { return; } - - //Let's draw a positronium decay for current event - - if ( fParaPositroniumFraction >= G4UniformRand() ) { fPositroniumKind = PositroniumKind::pPs; } - else { fPositroniumKind = PositroniumKind::oPs; } + shifted_particle_position = AddPositronRangeShift(particle_position, fModelParams.fMeanPositronRange[decayIndex]); } - switch ( fPositroniumKind ) - { - case PositroniumKind::pPs: - pInfoPs = &fParaPs; - break; - case PositroniumKind::oPs: - pInfoPs = &fOrthoPs; - break; - default: - NoticeError( G4String( __FUNCTION__ ), "improper chosen positronium kind." ); - break; - }; + G4PrimaryVertex* vertex = new G4PrimaryVertex( shifted_particle_position, shifted_particle_time ); + std::vector gammas = GetGammasFromPositroniumAnnihilation(decayIndex); + std::for_each( gammas.begin(), gammas.end(), [&]( G4PrimaryParticle* gamma ) { vertex->SetPrimary( gamma ); } ); + return vertex; } -GateEmittedGammaInformation* GatePositroniumDecayModel::GetPrimaryParticleInformation( const G4PrimaryParticle* pp, const GateEmittedGammaInformation::GammaKind& gamma_kind ) const +G4int GatePositroniumDecayModel::GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position) { - GateEmittedGammaInformation* egi = new GateEmittedGammaInformation(); - - GateEmittedGammaInformation::SourceKind source_kind = GateEmittedGammaInformation::SourceKind::ParaPositronium; - GateEmittedGammaInformation::DecayModel decay_model = GateEmittedGammaInformation::DecayModel::Standard; - - if ( fPositroniumKind == PositroniumKind::oPs ) { source_kind = GateEmittedGammaInformation::SourceKind::OrthoPositronium; } - if ( fDecayModel == DecayModel::WithPrompt ) { decay_model = GateEmittedGammaInformation::DecayModel::Deexcitation; } - - egi->SetSourceKind( source_kind ); - egi->SetDecayModel( decay_model ); - egi->SetGammaKind( gamma_kind ); - egi->SetInitialPolarization( pp->GetPolarization() ); - - if ( gamma_kind == GateEmittedGammaInformation::GammaKind::Annihilation ){ egi->SetTimeShift( pInfoPs->GetLifeTime() ); } - - return egi; -} - -G4int GatePositroniumDecayModel::GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position ) + G4int number_of_vertices = 0; + while (number_of_vertices <=0) { + auto decayIndex = GatePositroniumDecayModel::getPositroniumDecayIndex(fModelParams.fFractions); + auto no_electron_capture_prob = 1- fModelParams.fElectronCaptureProbabilities[decayIndex]; + assert(no_electron_capture_prob>=0); + + if(G4UniformRand() <= fModelParams.fPromptGammaProbabilities[decayIndex]) + { + ++number_of_vertices; + event->AddPrimaryVertex(GetPrimaryVertexFromDeexcitation(particle_time, particle_position, decayIndex)); + } + + if(G4UniformRand() <= no_electron_capture_prob) + { + ++number_of_vertices; + event->AddPrimaryVertex(GetPrimaryVertexFromPositroniumAnnihilation(particle_time, particle_position, decayIndex)); + } + } + return number_of_vertices; +} + +GateEmittedGammaInformation::DecayModel GatePositroniumDecayModel::GetDecayModel(const int decayIndex) const { - PreparePositroniumParametrization(); - G4int vertexes_number = 1; - - if ( fDecayModel == DecayModel::WithPrompt ) - { - ++vertexes_number; - event->AddPrimaryVertex( GetPrimaryVertexFromDeexcitation(particle_time, particle_position) ); - } - - event->AddPrimaryVertex( GetPrimaryVertexFromPositroniumAnnihilation(particle_time, particle_position) ); - - //Do testu - /*G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position, particle_time); - std::vector gammas_ps = GetGammasFromPositroniumAnnihilation(); - vertex->SetPrimary( GetGammaFromDeexcitation() ); - std::for_each( gammas_ps.begin(), gammas_ps.end(), [&]( G4PrimaryParticle* gamma ) { vertex->SetPrimary( gamma ); } );*/ - - return vertexes_number; + if (fModelParams.fPromptGammaProbabilities[decayIndex] > 0) { + return GateEmittedGammaInformation::DecayModel::Deexcitation; + } + return GateEmittedGammaInformation::DecayModel::Standard; } -G4PrimaryVertex* GatePositroniumDecayModel::GetPrimaryVertexFromDeexcitation(const G4double& particle_time, const G4ThreeVector& particle_position ) +GateEmittedGammaInformation::SourceKind GatePositroniumDecayModel::GetSourceKind(int decayIndex) const { - G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position, particle_time); - vertex->SetPrimary( GetGammaFromDeexcitation() ); - return vertex; + const PositronElectronInteraction interaction = fModelParams.fPositronInteractions[decayIndex]; + + switch (interaction) + { + case PositronElectronInteraction::kParaPs: + return GateEmittedGammaInformation::SourceKind::ParaPositronium; + case PositronElectronInteraction::kOrthoPs: + return GateEmittedGammaInformation::SourceKind::OrthoPositronium; + case PositronElectronInteraction::kDirect: + return GateEmittedGammaInformation::SourceKind::DirectAnnihilation; + default: + break; + } + return GateEmittedGammaInformation::SourceKind::NotDefined; } -G4PrimaryVertex* GatePositroniumDecayModel::GetPrimaryVertexFromPositroniumAnnihilation( const G4double& particle_time, const G4ThreeVector& particle_position ) +G4PrimaryParticle* GatePositroniumDecayModel::GetGammaFromDeexcitation(int decayIndex) { - G4double shifted_particle_time = particle_time + G4RandExponential::shoot( pInfoPs->GetLifeTime() ); - - G4PrimaryVertex* vertex = new G4PrimaryVertex( particle_position, shifted_particle_time ); - std::vector gammas = GetGammasFromPositroniumAnnihilation(); - std::for_each( gammas.begin(), gammas.end(), [&]( G4PrimaryParticle* gamma ) { vertex->SetPrimary( gamma ); } ); - return vertex; + G4PrimaryParticle* gamma = GetSingleGamma(fModelParams.fPromptGammaEnergy[decayIndex]); + GateEmittedGammaInformation* info = GetPrimaryParticleInformation( gamma, GateEmittedGammaInformation::GammaKind::Prompt ); + info->SetDecayIndex( decayIndex ); + info->SetDecayModel( GetDecayModel(decayIndex) ); + info->SetSourceKind( GetSourceKind(decayIndex) ); + gamma->SetUserInformation( info ); + return gamma; } -G4PrimaryParticle* GatePositroniumDecayModel::GetGammaFromDeexcitation() +G4PrimaryVertex* GatePositroniumDecayModel::GetPrimaryVertexFromDeexcitation(G4double particle_time, const G4ThreeVector& particle_position, int decayIndex) { - G4PrimaryParticle* gamma = GetSingleGamma( fPromptGammaEnergy ); - gamma->SetUserInformation( GetPrimaryParticleInformation( gamma, GateEmittedGammaInformation::GammaKind::Prompt ) ); - return gamma; + G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position, particle_time); + vertex->SetPrimary(GetGammaFromDeexcitation(decayIndex)); + return vertex; } -std::vector GatePositroniumDecayModel::GetGammasFromPositroniumAnnihilation() +std::vector GatePositroniumDecayModel::GetGammasFromPositroniumAnnihilation(int decayIndex) { - std::vector gammas( pInfoPs->GetAnnihilationGammasNumber() ) ; + int annihilation_gammas_number = fPositroniumDecayChannel[decayIndex].GetAnnihilationGammasNumber(); + std::vector gammas(annihilation_gammas_number); - G4DecayProducts* decay_products = pInfoPs->GetDecayProducts(); - for ( G4int i = 0; i < pInfoPs->GetAnnihilationGammasNumber(); ++i ) + G4DecayProducts* decay_products = fPositroniumDecayChannel[decayIndex].GetDecayProducts(); + for ( G4int i = 0; i < annihilation_gammas_number; ++i ) { G4PrimaryParticle* gamma = new G4PrimaryParticle( pGammaDefinition ); @@ -173,7 +149,11 @@ std::vector GatePositroniumDecayModel::GetGammasFromPositron G4LorentzVector lv = dynamic_gamma->Get4Momentum(); gamma->Set4Momentum( lv.px(), lv.py(), lv.pz(), lv.e() ); gamma->SetPolarization( dynamic_gamma->GetPolarization() ); - gamma->SetUserInformation( GetPrimaryParticleInformation( gamma, GateEmittedGammaInformation::GammaKind::Annihilation ) ); + GateEmittedGammaInformation* info = GetPrimaryParticleInformation( gamma, GateEmittedGammaInformation::GammaKind::Annihilation ); + info->SetDecayIndex( decayIndex ); + info->SetDecayModel( GetDecayModel(decayIndex) ); + info->SetSourceKind( GetSourceKind(decayIndex) ); + gamma->SetUserInformation( info ); gammas[i] = gamma; } delete decay_products; @@ -181,4 +161,3 @@ std::vector GatePositroniumDecayModel::GetGammasFromPositron return gammas; } - diff --git a/source/physics/src/GatePositroniumDecayParamsGenerator.cc b/source/physics/src/GatePositroniumDecayParamsGenerator.cc new file mode 100644 index 000000000..26c042d72 --- /dev/null +++ b/source/physics/src/GatePositroniumDecayParamsGenerator.cc @@ -0,0 +1,135 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include + +#include "GateMessageManager.hh" +#include "GatePositroniumHelper.hh" +#include "GatePositroniumDecayParamsGenerator.hh" + +void GatePositroniumDecayParamsGenerator::SetElectronCaptureProbabilities(const std::vector& electronCaptureProb) +{ + fElectronCaptureProbabilities = electronCaptureProb; +} + +void GatePositroniumDecayParamsGenerator::SetPromptGammaProbabilities(const std::vector& promptGammaProb) +{ + fPromptGammaProbabilities = promptGammaProb; +} +void GatePositroniumDecayParamsGenerator::SetDecayKinds(const std::vector& decayKinds) +{ + fDecayKinds = decayKinds; +} +void GatePositroniumDecayParamsGenerator::SetPositronInteractions(const std::vector& positronInteractions) +{ + fPositronInteractions = positronInteractions; +} +void GatePositroniumDecayParamsGenerator::SetPositroniumLifetimes(const std::vector& positroniumLifetimes) +{ + fPositroniumLifetimes = positroniumLifetimes; +} +void GatePositroniumDecayParamsGenerator::SetPromptGammaEnergies(const std::vector& energies) +{ + fPromptGammaEnergies = energies; +} +void GatePositroniumDecayParamsGenerator::SetPositroniumFraction(const std::vector& positroniumFractions) +{ + fPositroniumFractions = positroniumFractions; +} + +void GatePositroniumDecayParamsGenerator::SetMeanPositronRange(const std::vector& meanPositironRange) +{ + fMeanPositronRange = meanPositironRange; +} + + +PositroniumDecayModelParams GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams() const +{ + PositroniumDecayModelParams params; + GatePositroniumHelper positronHelper; + if(fPositroniumFractions.has_value()) { + params.fFractions=positronHelper.NormalizeFractions(fPositroniumFractions.value()); + } else { + GateError("GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: Positronium decay fractions are not set"); + } + + if(fPositroniumLifetimes.has_value()) { + params.fLifetimes=fPositroniumLifetimes.value(); + } else { + GateError("GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: Positronium lifetimes are not set"); + } + + if(fPositronInteractions.has_value()) { + params.fPositronInteractions=fPositronInteractions.value(); + } else if (!fDecayKinds.has_value()) { + GateError("GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: Positronium interactions are not set and decay kinds are also empty"); + } + + if(fPromptGammaProbabilities.has_value()) { + params.fPromptGammaProbabilities=fPromptGammaProbabilities.value(); + } else { + GateError("GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: Positronium prompt gamma probabilities are not set"); + } + + if(fPromptGammaEnergies.has_value()) { + params.fPromptGammaEnergy=fPromptGammaEnergies.value(); + } else { + GateError("GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: Prompt gamma energies are not set"); + } + + auto nElements = fPositroniumFractions.value().size(); + if(fElectronCaptureProbabilities.has_value()) { + params.fElectronCaptureProbabilities=fElectronCaptureProbabilities.value(); + } else { + params.fElectronCaptureProbabilities.assign(nElements, 0.0); + } + + if(fDecayKinds.has_value()) { + params.fDecayKind=fDecayKinds.value(); + } else if (!params.fPositronInteractions.empty() && !params.fLifetimes.empty() && !params.fFractions.empty()) { + params = positronHelper.CalculateFractionsFromLifetimes(params); + } else { + GateError("GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: Positronium decay kinds are not set and one or more sets that can calculate them (positronInteractions, lifetimes or fractions) is/are empty"); + } + + if(fMeanPositronRange.has_value()) { + params.fMeanPositronRangeEnabled.assign(nElements, true); + params.fMeanPositronRange=fMeanPositronRange.value(); + } else { + params.fMeanPositronRangeEnabled.assign(nElements, false); + params.fMeanPositronRange.assign(nElements, 0.0); + } + +validatePositroniumDecayParams(params); + +return params; +} + +void GatePositroniumDecayParamsGenerator::validatePositroniumDecayParams(const PositroniumDecayModelParams& params) const +{ + auto ref_param_number = params.fDecayKind.size(); + bool size_mismatch = (ref_param_number != params.fFractions.size()) || + (ref_param_number != params.fPromptGammaProbabilities.size()) || + (ref_param_number != params.fLifetimes.size()) || + (ref_param_number != params.fPromptGammaEnergy.size()) || + (ref_param_number != params.fElectronCaptureProbabilities.size()); + if (size_mismatch) { + std::cout << ref_param_number << " " << params.fFractions.size() << " " << params.fPromptGammaProbabilities.size() << " " << params.fLifetimes.size() << " " << params.fPromptGammaEnergy.size() << params.fElectronCaptureProbabilities.size() << std::endl; + GateError( + "GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: " + "number of provided parameters in Fractions, PromptGamma, Lifetimes, " + "Gamma Energies, Electron Capture Probabilites are not the same"); + } + constexpr double kProb_of_no_particle_limit = 0.9; + for (int i = 0; i < ref_param_number; i++) { + auto electron_capture_prob = params.fElectronCaptureProbabilities[i]; + auto no_prompt_prob = 1 -params.fPromptGammaProbabilities[i]; + if (electron_capture_prob *no_prompt_prob >= kProb_of_no_particle_limit) { + GateWarning("The probability of 0 prompt emission times probability of electron capture (no anihillation) is larger than 90%. Please check channel definitions! Otherwise simulations can take a lot of time"); + } + } +} diff --git a/source/physics/src/GatePositroniumHelper.cc b/source/physics/src/GatePositroniumHelper.cc new file mode 100644 index 000000000..8685d0a87 --- /dev/null +++ b/source/physics/src/GatePositroniumHelper.cc @@ -0,0 +1,131 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include +#include +#include + +#include "GateMessageManager.hh" +#include "GatePositroniumHelper.hh" +#include "GatePositroniumConstants.hh" + +using namespace GatePositroniumConstants; + +PositroniumDecayModelParams GatePositroniumHelper::CalculateFractionsFromLifetimes(PositroniumDecayModelParams params) +{ + PositroniumDecayModelParams paramsOut; + bool pPsExist = false; + for (unsigned i=0; i intens = CalcFractionsFromLifetime(params.fFractions.at(i), params.fLifetimes.at(i), params.fPositronInteractions.at(i)); + paramsOut.fFractions.push_back(intens.first); //2G intens + paramsOut.fFractions.push_back(intens.second); //3G intens + paramsOut.fLifetimes.push_back(params.fLifetimes.at(i)); + paramsOut.fLifetimes.push_back(params.fLifetimes.at(i)); + paramsOut.fPromptGammaProbabilities.push_back(params.fPromptGammaProbabilities.at(i)); + paramsOut.fPromptGammaProbabilities.push_back(params.fPromptGammaProbabilities.at(i)); + paramsOut.fPromptGammaEnergy.push_back(params.fPromptGammaEnergy.at(i)); + paramsOut.fPromptGammaEnergy.push_back(params.fPromptGammaEnergy.at(i)); + paramsOut.fDecayKind.push_back(PositroniumDecayKind::k2Gamma); + paramsOut.fDecayKind.push_back(PositroniumDecayKind::k3Gamma); + paramsOut.fPositronInteractions.push_back(params.fPositronInteractions.at(i)); + paramsOut.fPositronInteractions.push_back(params.fPositronInteractions.at(i)); + } + } + if (!paramsOut.fPositronInteractions.size()) { + GateError("GatePositroniumHelper::CalculateFractionsFromLifetimes: Could not calculate fraction from lifetimes"); + return params; + } + +// setting pPs + float pPsIntens = CalcPPsFractionFromOPs(paramsOut.fFractions, paramsOut.fPositronInteractions); + if (pPsIntens > 0) { + paramsOut.fFractions.push_back(pPsIntens); + paramsOut.fLifetimes.push_back(kParaPsLifetime_ns); + paramsOut.fPromptGammaProbabilities.push_back(paramsOut.fPromptGammaProbabilities.at(0)); + paramsOut.fPromptGammaEnergy.push_back(paramsOut.fPromptGammaEnergy.at(0)); + paramsOut.fDecayKind.push_back(PositroniumDecayKind::k2Gamma); + paramsOut.fPositronInteractions.push_back(PositronElectronInteraction::kParaPs); + } + + paramsOut.fFractions = NormalizeFractions(paramsOut.fFractions); + + return paramsOut; +} + +float GatePositroniumHelper::CalcPPsFractionFromOPs(std::vector fractions, std::vector decays) { + float sum = 0.; + auto itFrac = fractions.begin(); + auto itDec = decays.begin(); + while (itFrac != fractions.end() && itDec != decays.end()) { + if (*itDec == PositronElectronInteraction::kOrthoPs) + sum += *itFrac; + ++itFrac; + ++itDec; + } + return sum*kParaToOrthoPsFraction ; +} + +std::pair GatePositroniumHelper::CalcFractionsFromLifetime(float intensity, float lifetime, PositronElectronInteraction inter) { + float intens2G = 1., intens3G = 0.; + switch (inter) { + case PositronElectronInteraction::kDirect: + intens2G = intensity*(kHyperfineCoefficient - 1.)/kHyperfineCoefficient; + intens3G = intensity/kHyperfineCoefficient; + break; + case PositronElectronInteraction::kOrthoPs: + intens2G = intensity*(kOrthoPsMeanLifetime_ns - lifetime)/kOrthoPsMeanLifetime_ns; + intens3G = intensity*lifetime/kOrthoPsMeanLifetime_ns; + break; + } + return std::make_pair(intens2G, intens3G); +} + +float GatePositroniumHelper::CalcFractionFromOPsLifetime(float intensity, float lifetime, PositroniumDecayKind decay) { + float nominator = 0.; + if (intensity > 0 && lifetime > 0) { + switch (decay) { + case PositroniumDecayKind::k2Gamma: + nominator = kOrthoPsMeanLifetime_ns - lifetime; + break; + case PositroniumDecayKind::k3Gamma: + nominator = lifetime; + break; +// If there will be more decays one needs to modify it + } + return intensity*nominator/kOrthoPsMeanLifetime_ns ; + } else + return nominator; +} + +float GatePositroniumHelper::CalcFractionFromDirectLifetime(float intensity, PositroniumDecayKind decay) { + float nominator = 0.; + if (intensity > 0) { + switch (decay) { + case PositroniumDecayKind::k2Gamma: + nominator = kHyperfineCoefficient - 1.; + break; + case PositroniumDecayKind::k3Gamma: + nominator = 1.; + break; + } + return intensity*nominator/kHyperfineCoefficient ; + } else + return nominator; +} + +std::vector GatePositroniumHelper::NormalizeFractions(std::vector fractions) { + std::vector normalizedVector = fractions; + double sum = std::accumulate(fractions.begin(), fractions.end(), 0.0); + if (sum > 0) + std::transform(fractions.begin(), fractions.end(), normalizedVector.begin(), [sum](float element){ return element/sum; }); + else { +//ThrowException + normalizedVector = fractions; + } + return normalizedVector; +} diff --git a/source/physics/src/GatePositroniumSource.cc b/source/physics/src/GatePositroniumSource.cc new file mode 100644 index 000000000..317ecf17e --- /dev/null +++ b/source/physics/src/GatePositroniumSource.cc @@ -0,0 +1,41 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#include + +#include "G4Event.hh" + +#include "GatePositroniumSource.hh" +#include "GatePositroniumDecayModel.hh" + + +GatePositroniumSource::GatePositroniumSource(const G4String &name) + : GateVSource(name), + pMessenger(std::make_unique(this)) +{ +} + +void GatePositroniumSource::PrepareModel() +{ + if (GetType() != "Ps") { + GateError("GatePositroniumSource::PrepareModel - model type is not Positronium. Current type: " + GetType()); + } + + auto params = pMessenger->generatePositroniumDecayParams(); + pModel = std::make_unique(params); +} + +G4int GatePositroniumSource::GeneratePrimaries(G4Event* event) +{ + + if (!pModel) { PrepareModel(); } + + G4double particle_time = GetTime(); + G4ThreeVector particle_position = GetPosDist()->GenerateOne(); + ChangeParticlePositionRelativeToAttachedVolume(particle_position); + return pModel->GeneratePrimaryVertices(event, particle_time, particle_position); +} + diff --git a/source/physics/src/GatePositroniumSourceMessenger.cc b/source/physics/src/GatePositroniumSourceMessenger.cc new file mode 100644 index 000000000..0e3abb8bd --- /dev/null +++ b/source/physics/src/GatePositroniumSourceMessenger.cc @@ -0,0 +1,195 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include + +#include "GatePositroniumSourceMessenger.hh" +#include "GatePositroniumSource.hh" + +GatePositroniumSourceMessenger::GatePositroniumSourceMessenger(GatePositroniumSource *source): GateVSourceMessenger(source), pSource(source) +{ + InitCommands(); +} + +G4UIcmdWithABool* GatePositroniumSourceMessenger::GetBoolCmd(const G4String& cmd_name, const G4String& cmd_guidance ) +{ + G4String cmd_path = GetDirectoryName() + cmd_name; + G4UIcmdWithABool* cmd = new G4UIcmdWithABool( cmd_path, this ); + cmd->SetGuidance( cmd_guidance ); + cmd->SetParameterName( cmd_name, false ); + return cmd; +} + +G4UIcmdWithADoubleAndUnit* GatePositroniumSourceMessenger::GetDoubleCmdWithUnit( const G4String& cmd_name, const G4String& cmd_guidance, const G4String& default_unit, const G4String& unit_candidates ) +{ + G4String cmd_path = GetDirectoryName() + cmd_name; + G4UIcmdWithADoubleAndUnit* cmd = new G4UIcmdWithADoubleAndUnit( cmd_path , this ); + cmd->SetGuidance( cmd_guidance ); + cmd->SetParameterName( cmd_name, false ); + cmd->SetDefaultUnit( default_unit.c_str() ); + cmd->SetUnitCandidates( unit_candidates.c_str() ); + return cmd; +} + +G4UIcmdWith3Vector* GatePositroniumSourceMessenger::GetVectorCmd( const G4String& cmd_name, const G4String& cmd_guidance ) +{ + G4String cmd_path = GetDirectoryName() + cmd_name; + G4UIcmdWith3Vector* cmd = new G4UIcmdWith3Vector( cmd_path, this ); + cmd->SetGuidance( cmd_guidance ); + cmd->SetParameterName( G4String( cmd_name + "_x" ), G4String( cmd_name + "_y" ), G4String( cmd_name + "_z" ), false ); + return cmd; +} + +G4UIcmdWithAnInteger* GatePositroniumSourceMessenger::GetIntCmd( const G4String& cmd_name, const G4String& cmd_guidance ) +{ + G4String cmd_path = GetDirectoryName() + cmd_name; + G4UIcmdWithAnInteger* cmd = new G4UIcmdWithAnInteger( cmd_path, this ); + cmd->SetGuidance( cmd_guidance ); + cmd->SetParameterName( cmd_name, false ); + return cmd; +} + +G4UIcmdWithAString* GatePositroniumSourceMessenger::GetStringCmd(const G4String& cmd_name, const G4String& cmd_guidance ) +{ + G4String cmd_path = GetDirectoryName() + cmd_name; + G4UIcmdWithAString* cmd = new G4UIcmdWithAString( cmd_path, this ); + cmd->SetGuidance( cmd_guidance ); + cmd->SetParameterName( cmd_name, false ); + return cmd; +} + +G4UIcmdWith3VectorAndUnit* GatePositroniumSourceMessenger::GetVectorCmdWithUnit( const G4String& cmd_name, const G4String& cmd_guidance, const G4String& default_unit, const G4String& unit_candidates ) +{ + G4String cmd_path = GetDirectoryName() + cmd_name; + G4UIcmdWith3VectorAndUnit* cmd = new G4UIcmdWith3VectorAndUnit( cmd_path , this ); + cmd->SetGuidance( cmd_guidance ); + cmd->SetParameterName( G4String( cmd_name + "_x" ), G4String( cmd_name + "_y" ), G4String( cmd_name + "_z" ), false ); + cmd->SetDefaultUnit( default_unit.c_str() ); + cmd->SetUnitCandidates( unit_candidates.c_str() ); + return cmd; +} + +void GatePositroniumSourceMessenger::InitCommands() +{ + upCmdSetPositroniumFractions.reset(GetStringCmd( "setPositroniumFractions", "\"f1, f2, f3 .., fn\" - where fi in [0.0, 1.0] and sum of all fi ==1" ) ); + upCmdSetPositroniumLifetimes.reset(GetStringCmd( "setPositroniumLifetimes", "\"t1, t2, t3 .., tn t_unit\" - where ti corresponds to lifetime constants and t_unit is one of the Geant4 time units e.g. ns" ) ); + upCmdSetDecayKinds.reset(GetStringCmd( "setDecayKinds", "\"k1, k2, k3 .., kn\" - where ki is k2Gamma or k3Gamma" ) ); + upCmdSetPositronInteractions.reset(GetStringCmd( "setPositronInteractions", "\"k1, k2, k3 .., kn\" - where ki is kParaPs, kDirect or kOrthoPs, of a given element in vector of components. Used to properly recalculate intensities of the components from the theory" ) ); + upCmdSetPromptPhotonProbabilites.reset(GetStringCmd( "setPromptPhotonProbabilites", "\"f1, f2, f3 .., fn\" - where fi in [0.0, 1.0] " ) ); + upCmdSetPromptPhotonEnergies.reset(GetStringCmd( "setPromptPhotonEnergies", "\"e1, e2, e3 .., en e_unit\" - where ei are energies and e_unit is one of the Geant4 energy units e.g. MeV" ) ); + upCmdSetMeanPositronRange.reset(GetStringCmd( "setMeanPositronRange", "\"r1, r2, r3 .., rn r_unit\" - where ri are mean positron range and r_unit is one of the Geant 4 distance units e.g. mm " ) ); + upCmdSetElectronCaptureProbabilities.reset(GetStringCmd( "setElectronCaptureProbabilities", "\"f1, f2, f3 .., fn\" - where fi in [0.0, 1.0] " ) ); +} + +std::vector GatePositroniumSourceMessenger::parseListOfParamsWithUnit(const G4String& input) const +{ + std::vector values; + std::stringstream ss(input); + std::vector tokens; + std::string token; + + while (ss >> token) { + tokens.push_back(token); + } + + if (tokens.size() < 2) { + GateError("parseListOfParamsWithUnit(): Need at least one value and one unit."); + } + + G4String unitStr = tokens.back(); + G4double unitValue = G4UnitDefinition::GetValueOf(unitStr); + + for (unsigned i=0; i fractions; + std::stringstream ss(new_value); + G4double num; + while (ss >> num) { + fractions.push_back(num); + } + fParamGenerator.SetPositroniumFraction(fractions); + } else if (command == upCmdSetPositroniumLifetimes.get()) { + auto lifetimes = parseListOfParamsWithUnit(new_value); + fParamGenerator.SetPositroniumLifetimes(lifetimes); + } else if (command == upCmdSetPromptPhotonProbabilites.get()) { + std::vector promptPhotonProb; + std::stringstream ss(new_value); + float prob; + while (ss >> prob) { + prob = prob < 0 ? 0 : prob; + prob = prob > 1 ? 1 : prob; + promptPhotonProb.push_back(prob); + } + fParamGenerator.SetPromptGammaProbabilities(promptPhotonProb); + } else if (command == upCmdSetPromptPhotonEnergies.get()) { + auto promptPhotonEnergies = parseListOfParamsWithUnit(new_value); + fParamGenerator.SetPromptGammaEnergies(promptPhotonEnergies); + } else if (command == upCmdSetMeanPositronRange.get()) { + auto positronRange = parseListOfParamsWithUnit(new_value); + fParamGenerator.SetMeanPositronRange(positronRange); + } else if (command == upCmdSetElectronCaptureProbabilities.get()) { + std::vector electronCaptureProb; + std::stringstream ss(new_value); + float prob; + while (ss >> prob) { + prob = prob < 0 ? 0 : prob; + prob = prob > 1 ? 1 : prob; + electronCaptureProb.push_back(prob); + } + fParamGenerator.SetElectronCaptureProbabilities(electronCaptureProb); + } else if (command == upCmdSetDecayKinds.get()) { + std::vector decayKinds; + std::stringstream ss(new_value); + std::string kind; + while (ss >> kind) { + if (kind == "k2Gamma") { + decayKinds.push_back(k2Gamma); + } else { + if (kind == "k3Gamma") { + decayKinds.push_back(k3Gamma); + } else { + GateError("GatePositroniumSourceMessenger::SetNewValue: unknown " + "decay kind read from macro."); + } + } + } + fParamGenerator.SetDecayKinds(decayKinds); + } else if (command == upCmdSetPositronInteractions.get()) { + std::vector positronInteractions; + std::stringstream ss(new_value); + std::string inter; + while (ss >> inter) { + if (inter == "kParaPs") { + positronInteractions.push_back(PositronElectronInteraction::kParaPs); + } else if (inter == "kDirect") { + positronInteractions.push_back(PositronElectronInteraction::kDirect); + } else if (inter == "kOrthoPs") { + positronInteractions.push_back(PositronElectronInteraction::kOrthoPs); + } else { + + GateError("GatePositroniumSourceMessenger::SetNewValue: unknown " + "interaction type read from macro."); + } + } + fParamGenerator.SetPositronInteractions(positronInteractions); + } else { + GateVSourceMessenger::SetNewValue(command, new_value); + } +} + +PositroniumDecayModelParams GatePositroniumSourceMessenger::generatePositroniumDecayParams() const +{ + return fParamGenerator.generatePositroniumDecayParams(); +} diff --git a/source/physics/src/GateSourceMgr.cc b/source/physics/src/GateSourceMgr.cc index db8f192ce..b71f162f2 100644 --- a/source/physics/src/GateSourceMgr.cc +++ b/source/physics/src/GateSourceMgr.cc @@ -29,7 +29,8 @@ #include "G4RunManager.hh" #include "GateSourceOfPromptGamma.hh" #include "GateSourcePhaseSpace.hh" -#include "GateExtendedVSource.hh" +#include "legacy/GateExtendedVSource.hh" +#include "GatePositroniumSource.hh" //---------------------------------------------------------------------------------------- GateSourceMgr* GateSourceMgr::mInstance = 0; @@ -234,7 +235,11 @@ G4int GateSourceMgr::AddSource( std::vector sourceVec ) source->SetIfSourceVoxelized(false); // added by I. Martinez-Rovira (immamartinez@gmail.com) } else if (sourceGeomType == "Extended"){ - source = new GateExtendedVSource( sourceName ); + source = new GateLegacy::GateExtendedVSource( sourceName ); + source->SetSourceID( m_sourceProgressiveNumber ); + } + else if (sourceGeomType == "PositroniumSource"){ + source = new GatePositroniumSource( sourceName ); source->SetSourceID( m_sourceProgressiveNumber ); } else { diff --git a/source/physics/src/legacy/GateEmittedGammaInformation.cc b/source/physics/src/legacy/GateEmittedGammaInformation.cc new file mode 100644 index 000000000..912984197 --- /dev/null +++ b/source/physics/src/legacy/GateEmittedGammaInformation.cc @@ -0,0 +1,37 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#include "legacy/GateEmittedGammaInformation.hh" + +namespace GateLegacy{ + +GateEmittedGammaInformation::GateEmittedGammaInformation() {} + +GateEmittedGammaInformation::~GateEmittedGammaInformation() {} + +void GateEmittedGammaInformation::SetSourceKind( GateEmittedGammaInformation::SourceKind source_kind ) { fSourceKind = source_kind; } + +GateEmittedGammaInformation::SourceKind GateEmittedGammaInformation::GetSourceKind() const { return fSourceKind; } + +void GateEmittedGammaInformation::SetDecayModel( GateEmittedGammaInformation::DecayModel decay_model ) { fDecayModel = decay_model; } + +GateEmittedGammaInformation::DecayModel GateEmittedGammaInformation::GetDecayModel() const { return fDecayModel; } + +void GateEmittedGammaInformation::SetGammaKind( GateEmittedGammaInformation::GammaKind gamma_kind ) { fGammaKind = gamma_kind; } + +GateEmittedGammaInformation::GammaKind GateEmittedGammaInformation::GetGammaKind() const { return fGammaKind; } + +void GateEmittedGammaInformation::SetInitialPolarization( const G4ThreeVector& polarization ) { fInitialPolarization = polarization; } + +G4ThreeVector GateEmittedGammaInformation::GetInitialPolarization() const { return fInitialPolarization; } + +void GateEmittedGammaInformation::SetTimeShift( const G4double& time_shift ) { fTimeShift = time_shift; } + +G4double GateEmittedGammaInformation::GetTimeShift() const { return fTimeShift; } + +void GateEmittedGammaInformation::Print() const { G4cout << "GateEmittedGammaInformation" << G4endl; } + +} diff --git a/source/physics/src/GateExtendedVSource.cc b/source/physics/src/legacy/GateExtendedVSource.cc similarity index 97% rename from source/physics/src/GateExtendedVSource.cc rename to source/physics/src/legacy/GateExtendedVSource.cc index cb80afd2f..ddf840912 100644 --- a/source/physics/src/GateExtendedVSource.cc +++ b/source/physics/src/legacy/GateExtendedVSource.cc @@ -4,11 +4,13 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#include "GateExtendedVSource.hh" +#include "legacy/GateExtendedVSource.hh" #include -#include "GatePositroniumDecayModel.hh" +#include "legacy/GatePositroniumDecayModel.hh" #include "G4Event.hh" +namespace GateLegacy{ + GateExtendedVSource::GateExtendedVSource( G4String name ) : GateVSource( name ) { pMessenger = new GateExtendedVSourceMessenger( this ); @@ -92,7 +94,6 @@ void GateExtendedVSource::PrepareModel() if ( fFixedEmissionDirection.IsSetted() ) { pModel->SetFixedEmissionDirection( fFixedEmissionDirection.Get() ); } if ( fEnableFixedEmissionDirection.IsSetted() ) { pModel->SetEnableFixedEmissionDirection( fEnableFixedEmissionDirection.Get() ); } if ( fEmissionEnergy.IsSetted() ) { pModel->SetEmissionEnergy( fEmissionEnergy.Get() ); } - if ( fSeed.IsSetted() ) { pModel->SetSeed( fSeed.Get() ); } } @@ -106,4 +107,5 @@ G4int GateExtendedVSource::GeneratePrimaries( G4Event* event ) ChangeParticlePositionRelativeToAttachedVolume( particle_position ); return pModel->GeneratePrimaryVertices( event, particle_time, particle_position); } +} diff --git a/source/physics/src/GateExtendedVSourceMessenger.cc b/source/physics/src/legacy/GateExtendedVSourceMessenger.cc similarity index 98% rename from source/physics/src/GateExtendedVSourceMessenger.cc rename to source/physics/src/legacy/GateExtendedVSourceMessenger.cc index ff9c72f5d..00821e10e 100644 --- a/source/physics/src/GateExtendedVSourceMessenger.cc +++ b/source/physics/src/legacy/GateExtendedVSourceMessenger.cc @@ -5,10 +5,12 @@ See LICENSE.md for further details ----------------------*/ -#include "GateExtendedVSourceMessenger.hh" -#include "GateExtendedVSource.hh" +#include "legacy/GateExtendedVSourceMessenger.hh" +#include "legacy/GateExtendedVSource.hh" #include +namespace GateLegacy +{ GateExtendedVSourceMessenger::GateExtendedVSourceMessenger( GateExtendedVSource* source ) : GateVSourceMessenger( source ) { pSource = source; @@ -141,6 +143,6 @@ void GateExtendedVSourceMessenger::SetNewValue( G4UIcommand* command, G4String n } } - +} diff --git a/source/physics/src/legacy/GateGammaEmissionModel.cc b/source/physics/src/legacy/GateGammaEmissionModel.cc new file mode 100644 index 000000000..45f4845f5 --- /dev/null +++ b/source/physics/src/legacy/GateGammaEmissionModel.cc @@ -0,0 +1,146 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#include "legacy/GateGammaEmissionModel.hh" +#include +#include "G4LorentzVector.hh" +#include "Randomize.hh" +#include +#include "G4PrimaryVertex.hh" +#include "G4PrimaryParticle.hh" +#include "G4ParticleTable.hh" + +namespace GateLegacy{ + +GateGammaEmissionModel::GateGammaEmissionModel() { pGammaDefinition = G4ParticleTable::GetParticleTable()->FindParticle( "gamma" ); } + +GateGammaEmissionModel::~GateGammaEmissionModel() {} + +G4int GateGammaEmissionModel::GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position) +{ + G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position, particle_time); + G4PrimaryParticle* gamma = GetSingleGamma( fEmissionEnergy ); + gamma->SetUserInformation( GetPrimaryParticleInformation( gamma, GateEmittedGammaInformation::GammaKind::Single ) ); + vertex->SetPrimary( gamma ); + event->AddPrimaryVertex( vertex ); + + return 1; +} + +G4PrimaryParticle* GateGammaEmissionModel::GetSingleGamma( const G4double& energy ) const +{ + G4PrimaryParticle* gamma = new G4PrimaryParticle( pGammaDefinition ); + + G4ThreeVector momentum_direction; + + if ( fUseFixedEmissionDirection ) { momentum_direction = fFixedEmissionDirection; } + else { momentum_direction = GetUniformOnSphere(); } + + G4LorentzVector lv_gamma( momentum_direction.x(), momentum_direction.y(), momentum_direction.z(), 1.0 ); + lv_gamma *= energy; + //Update gamma + ///Momentum + gamma->Set4Momentum( lv_gamma.px(), lv_gamma.py(), lv_gamma.pz(), lv_gamma.e() ); + ///Polarization + gamma->SetPolarization( GetPolarization( gamma->GetMomentumDirection() ) ); + + return gamma; +} + +void GateGammaEmissionModel::SetFixedEmissionDirection( const G4ThreeVector& momentum_direction ) +{ + if ( momentum_direction.mag() == 0 ) + { + NoticeError( G4String( __FUNCTION__ ), "fixed momemntum direction vector should be non zero vector" ); + } + fFixedEmissionDirection = momentum_direction; + fFixedEmissionDirection = fFixedEmissionDirection.unit(); + fUseFixedEmissionDirection = true; +} + +void GateGammaEmissionModel::SetEnableFixedEmissionDirection( const G4bool enable ) { fUseFixedEmissionDirection = enable; } + +void GateGammaEmissionModel::SetEmissionEnergy( const G4double& energy ) +{ + if ( !( energy > 0.0 ) ) { NoticeError( G4String( __FUNCTION__ ), "emission energy should be positive value." ); } + fEmissionEnergy = energy; +} + +G4double GateGammaEmissionModel::GetEmissionEnergy() const { return fEmissionEnergy; } + +void GateGammaEmissionModel::SetSeed( G4long seed ) +{ + if ( seed < 0 ) { NoticeError( G4String( __FUNCTION__ ), "seed should be positive value." ); } + G4Random::setTheSeed( seed ); +} + +G4long GateGammaEmissionModel::GetSeed() const { return G4Random::getTheSeed (); } + +G4ThreeVector GateGammaEmissionModel::GetUniformOnSphere() const +{ + //Based on TRandom::Sphere + G4double a = 0,b = 0, r2 = 1; + while ( r2 > 0.25 ) + { + a = G4UniformRand() - 0.5; + b = G4UniformRand() - 0.5; + r2 = a*a + b*b; + } + + G4double scale = 8.0 * sqrt(0.25 - r2); + return G4ThreeVector( a * scale, b * scale, -1. + 8.0 * r2 ); +} + +G4ThreeVector GateGammaEmissionModel::GetPolarization( const G4ThreeVector& momentum ) const +{ + G4ThreeVector polarization(0.0,0.0,0.0); + + G4ThreeVector a0,b0,d0; + d0 = momentum.unit(); + a0 = GetPerpendicularVector( d0 ).unit(); + b0 = d0.cross( a0 ).unit(); + G4double angle_radians = G4UniformRand() * M_PI; + polarization = std::cos( angle_radians ) * a0 + std::sin( angle_radians ) * b0; + polarization.unit(); + return polarization; +} + +G4ThreeVector GateGammaEmissionModel::GetPerpendicularVector(const G4ThreeVector& v) const +{ + G4double dx = v.x(); + G4double dy = v.y(); + G4double dz = v.z(); + + G4double x = dx < 0.0 ? -dx : dx; + G4double y = dy < 0.0 ? -dy : dy; + G4double z = dz < 0.0 ? -dz : dz; + + if (x < y) { return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy); } + else { return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0); } +} + +void GateGammaEmissionModel::SetModelName( const G4String model_name ) +{ + if ( model_name.size() == 0 ) { NoticeError( G4String( __FUNCTION__ ), "not provided model name." ); } + fModelName = model_name; +} + +void GateGammaEmissionModel::NoticeError( G4String method_name, G4String exception_description ) const +{ + G4String error_msg = fModelName + "::" + method_name + " : " + exception_description; + GateError( error_msg ); +} + +GateEmittedGammaInformation* GateGammaEmissionModel::GetPrimaryParticleInformation( const G4PrimaryParticle* pp, const GateEmittedGammaInformation::GammaKind& gamma_kind ) const +{ + GateEmittedGammaInformation* egi = new GateEmittedGammaInformation(); + egi->SetSourceKind( GateEmittedGammaInformation::SourceKind::SingleGammaEmitter ); + egi->SetDecayModel( GateEmittedGammaInformation::DecayModel::None ); + egi->SetGammaKind( gamma_kind ); + egi->SetInitialPolarization( pp->GetPolarization() ); + return egi; +} +} diff --git a/source/physics/src/GateOrthoPositronium.cc b/source/physics/src/legacy/GateOrthoPositronium.cc similarity index 94% rename from source/physics/src/GateOrthoPositronium.cc rename to source/physics/src/legacy/GateOrthoPositronium.cc index 9c3fcf9ff..6fbb33619 100644 --- a/source/physics/src/GateOrthoPositronium.cc +++ b/source/physics/src/legacy/GateOrthoPositronium.cc @@ -4,13 +4,14 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#include "GateOrthoPositronium.hh" +#include "legacy/GateOrthoPositronium.hh" #include "G4PhysicalConstants.hh" #include "G4SystemOfUnits.hh" #include "G4ParticleTable.hh" #include "G4DecayTable.hh" -#include "GatePositroniumDecayChannel.hh" +#include "legacy/GatePositroniumDecayChannel.hh" +namespace GateLegacy{ GateOrthoPositronium* GateOrthoPositronium::theInstance = 0; @@ -62,3 +63,4 @@ GateOrthoPositronium* GateOrthoPositronium::Definition() GateOrthoPositronium* GateOrthoPositronium::OrthoPositroniumDefinition() { return Definition(); } GateOrthoPositronium* GateOrthoPositronium::OrthoPositronium() { return Definition(); } +} diff --git a/source/physics/src/GateParaPositronium.cc b/source/physics/src/legacy/GateParaPositronium.cc similarity index 94% rename from source/physics/src/GateParaPositronium.cc rename to source/physics/src/legacy/GateParaPositronium.cc index 1f304c4e6..dc58bc895 100644 --- a/source/physics/src/GateParaPositronium.cc +++ b/source/physics/src/legacy/GateParaPositronium.cc @@ -4,13 +4,14 @@ of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details ----------------------*/ -#include "GateParaPositronium.hh" +#include "legacy/GateParaPositronium.hh" #include "G4PhysicalConstants.hh" #include "G4SystemOfUnits.hh" #include "G4ParticleTable.hh" #include "G4DecayTable.hh" -#include "GatePositroniumDecayChannel.hh" +#include "legacy/GatePositroniumDecayChannel.hh" +namespace GateLegacy{ GateParaPositronium* GateParaPositronium::theInstance = 0; @@ -62,3 +63,4 @@ GateParaPositronium* GateParaPositronium::Definition() GateParaPositronium* GateParaPositronium::ParaPositroniumDefinition() { return Definition(); } GateParaPositronium* GateParaPositronium::ParaPositronium() { return Definition(); } +} diff --git a/source/physics/src/legacy/GatePositroniumDecayChannel.cc b/source/physics/src/legacy/GatePositroniumDecayChannel.cc new file mode 100644 index 000000000..04ef926bd --- /dev/null +++ b/source/physics/src/legacy/GatePositroniumDecayChannel.cc @@ -0,0 +1,160 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#include "legacy/GatePositroniumDecayChannel.hh" +#include "G4DynamicParticle.hh" +#include "G4DecayProducts.hh" +#include "Randomize.hh" +#include "G4LorentzVector.hh" + +namespace GateLegacy{ + +GatePositroniumDecayChannel::GatePositroniumDecayChannel(const G4String& theParentName, G4double theBR) +{ + G4int daughters_number = 0; + + if ( theParentName == kParaPositroniumName ) + { + daughters_number = kParaPositroniumAnnihilationGammasNumber; + fPositroniumKind = GatePositroniumDecayChannel::PositroniumKind::ParaPositronium; + } + else if ( theParentName == kOrthoPositroniumName ) + { + daughters_number = kOrthoPositroniumAnnihilationGammasNumber; + fPositroniumKind = GatePositroniumDecayChannel::PositroniumKind::OrthoPositronium; + } + else + { + #ifdef G4VERBOSE + if (GetVerboseLevel()>0) + { + G4cout << "GatePositroniumDecayChannel:: constructor :"; + G4cout << " parent particle is not positronium (pPs,oPs) but "; + G4cout << theParentName << G4endl; + } + #endif + } + + SetParentMass( kPositroniumMass ); + SetBR( theBR ); + SetParent( theParentName ); + SetNumberOfDaughters( daughters_number ); + for ( G4int daughter_index = 0; daughter_index < daughters_number; ++daughter_index ) { SetDaughter( daughter_index, kDaughterName ); } +} + +GatePositroniumDecayChannel::~GatePositroniumDecayChannel() {} + +G4DecayProducts* GatePositroniumDecayChannel::DecayIt(G4double) +{ + switch ( fPositroniumKind ) + { + case GatePositroniumDecayChannel::PositroniumKind::ParaPositronium: + return DecayParaPositronium(); + case GatePositroniumDecayChannel::PositroniumKind::OrthoPositronium: + return DecayOrthoPositronium(); + default: + G4cout << "GatePositroniumDecayChannel::DecayIt "; + G4cout << *parent_name << " can not decay " << G4endl; + DumpInfo(); + return nullptr; + }; +} + +G4DecayProducts* GatePositroniumDecayChannel::DecayParaPositronium() +{ + G4DecayProducts* decay_products = G4GeneralPhaseSpaceDecay::DecayIt(); + + G4DynamicParticle* gamma_1 = (*decay_products)[0]; + G4DynamicParticle* gamma_2 = (*decay_products)[1]; + + ///Polarization + G4ThreeVector polarization_gamma_1 = GetPolarization( gamma_1->GetMomentumDirection() ); + G4ThreeVector polarization_gamma_2 = gamma_1->GetMomentumDirection().cross( polarization_gamma_1 ); + + gamma_1->SetPolarization( polarization_gamma_1.x(), polarization_gamma_1.y(), polarization_gamma_1.z() ); + gamma_2->SetPolarization( polarization_gamma_2.x(), polarization_gamma_2.y(), polarization_gamma_2.z() ); + + return decay_products; +} + +G4DecayProducts* GatePositroniumDecayChannel::DecayOrthoPositronium() +{ + G4LorentzVector lv_gamma_1, lv_gamma_2, lv_gamma_3; + G4double weight = 0.0, random_weight = 0.0; + G4DecayProducts* decay_products = nullptr; + G4DynamicParticle* gamma_1 = nullptr; + G4DynamicParticle* gamma_2 = nullptr; + G4DynamicParticle* gamma_3 = nullptr; + + do + { + if ( decay_products != nullptr ) { delete decay_products; } + + decay_products = G4GeneralPhaseSpaceDecay::DecayIt(); + + gamma_1 = (*decay_products)[0]; + gamma_2 = (*decay_products)[1]; + gamma_3 = (*decay_products)[2]; + + lv_gamma_1 = gamma_1->Get4Momentum(); + lv_gamma_2 = gamma_2->Get4Momentum(); + lv_gamma_3 = gamma_3->Get4Momentum(); + + weight = GetOrthoPsM( lv_gamma_1.e(), lv_gamma_2.e(), lv_gamma_3.e() ); + random_weight = kOrthoPsMMax * G4UniformRand(); + } + while( random_weight > weight ); + + ///Polarization + G4ThreeVector polarization_gamma_1 = GetPolarization( gamma_1->GetMomentumDirection() ); + G4ThreeVector polarization_gamma_2 = GetPolarization( gamma_2->GetMomentumDirection() ); + G4ThreeVector polarization_gamma_3 = GetPolarization( gamma_3->GetMomentumDirection() ); + + gamma_1->SetPolarization( polarization_gamma_1.x(), polarization_gamma_1.y(), polarization_gamma_1.z() ); + gamma_2->SetPolarization( polarization_gamma_2.x(), polarization_gamma_2.y(), polarization_gamma_2.z() ); + gamma_3->SetPolarization( polarization_gamma_3.x(), polarization_gamma_3.y(), polarization_gamma_3.z() ); + + return decay_products; +} + +G4double GatePositroniumDecayChannel::GetOrthoPsM( const G4double w1, const G4double w2, const G4double w3 ) const +{ + return pow( ( kElectronMass - w1 ) / ( w2 * w3 ), 2 ) + pow( ( kElectronMass - w2 ) / ( w1 * w3 ), 2 ) + pow( ( kElectronMass - w3 ) / ( w1 * w2 ), 2 ); +} + +G4ThreeVector GatePositroniumDecayChannel::GetPolarization( const G4ThreeVector& momentum ) const +{ + G4ThreeVector polarization(0.0,0.0,0.0); + + G4ThreeVector a0,b0,d0; + d0 = momentum.unit(); + a0 = GetPerpendicularVector( d0 ).unit(); + b0 = d0.cross( a0 ).unit(); + G4double angle_radians = G4UniformRand() * M_PI; + polarization = std::cos( angle_radians ) * a0 + std::sin( angle_radians ) * b0; + polarization.unit(); + return polarization; +} + +G4ThreeVector GatePositroniumDecayChannel::GetPerpendicularVector(const G4ThreeVector& v) const +{ + G4double dx = v.x(); + G4double dy = v.y(); + G4double dz = v.z(); + + G4double x = dx < 0.0 ? -dx : dx; + G4double y = dy < 0.0 ? -dy : dy; + G4double z = dz < 0.0 ? -dz : dz; + + if (x < y) { return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy); } + else { return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0); } +} +} + + + + + diff --git a/source/physics/src/legacy/GatePositroniumDecayModel.cc b/source/physics/src/legacy/GatePositroniumDecayModel.cc new file mode 100644 index 000000000..c6f62b056 --- /dev/null +++ b/source/physics/src/legacy/GatePositroniumDecayModel.cc @@ -0,0 +1,187 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ +#include "legacy/GatePositroniumDecayModel.hh" +#include "Randomize.hh" +#include +#include +#include +#include "G4DecayProducts.hh" +#include "G4LorentzVector.hh" +#include "G4ParticleTable.hh" + +namespace GateLegacy +{ + +GatePositroniumDecayModel::Positronium::Positronium( G4String name, G4double life_time, G4int annihilation_gammas_number ) : fName( name ), fLifeTime( life_time ), fAnnihilationGammasNumber( annihilation_gammas_number ) +{ + G4ParticleDefinition* positronium_def = G4ParticleTable::GetParticleTable()->FindParticle( name ); + G4DecayTable* positronium_decay_table = positronium_def->GetDecayTable(); + pDecayChannel = positronium_decay_table->GetDecayChannel(0); +} + +void GatePositroniumDecayModel::Positronium::SetLifeTime( const G4double& life_time ) { fLifeTime = life_time; } + +G4double GatePositroniumDecayModel::Positronium::GetLifeTime() const { return fLifeTime; } + +G4String GatePositroniumDecayModel::Positronium::GetName() const { return fName; } + +G4int GatePositroniumDecayModel::Positronium::GetAnnihilationGammasNumber() const { return fAnnihilationGammasNumber; } + +G4DecayProducts* GatePositroniumDecayModel::Positronium::GetDecayProducts() { return pDecayChannel->DecayIt(); } + +GatePositroniumDecayModel::GatePositroniumDecayModel() +{ + SetModelName( "GatePositroniumDecayModel" ); +} + +GatePositroniumDecayModel::~GatePositroniumDecayModel() {} + +void GatePositroniumDecayModel::SetPositroniumKind( GatePositroniumDecayModel::PositroniumKind positronium_kind ) { fPositroniumKind = positronium_kind; } + +GatePositroniumDecayModel::PositroniumKind GatePositroniumDecayModel::GetPositroniumKind() const { return fPositroniumKind; } + +void GatePositroniumDecayModel::SetDecayModel( GatePositroniumDecayModel::DecayModel decay_model ) { fDecayModel = decay_model; } + +GatePositroniumDecayModel::DecayModel GatePositroniumDecayModel::GetDecayModel() const { return fDecayModel; } + +void GatePositroniumDecayModel::SetPostroniumLifetime( G4String positronium_name, G4double life_time ) +{ + if ( !( life_time > 0.0 ) ) { NoticeError( G4String( __FUNCTION__ ), "positronium life-time should be positive value." ); } + + if ( positronium_name == fParaPs.GetName() ) { fParaPs.SetLifeTime( life_time ); } + else if ( positronium_name == fOrthoPs.GetName() ) { fOrthoPs.SetLifeTime( life_time ); } + else { NoticeError( G4String( __FUNCTION__ ), "Unknown positronium name." ); } +} + +void GatePositroniumDecayModel::SetPromptGammaEnergy( G4double prompt_energy ) +{ + if ( !( prompt_energy > 0.0 ) ) { NoticeError( G4String( __FUNCTION__ ), "prompt gamma energy should be positive value." ); } + fPromptGammaEnergy = prompt_energy; +} + +G4double GatePositroniumDecayModel::GetPromptGammaEnergy() const { return fPromptGammaEnergy; } + +void GatePositroniumDecayModel::SetParaPositroniumFraction( G4double fraction ) +{ + fUsePositroniumFractions = true; + fParaPositroniumFraction = fraction; +} + +void GatePositroniumDecayModel::PreparePositroniumParametrization() +{ + if ( pInfoPs != nullptr ) + { + if ( !fUsePositroniumFractions ) { return; } + + //Let's draw a positronium decay for current event + + if ( fParaPositroniumFraction >= G4UniformRand() ) { fPositroniumKind = PositroniumKind::pPs; } + else { fPositroniumKind = PositroniumKind::oPs; } + } + + switch ( fPositroniumKind ) + { + case PositroniumKind::pPs: + pInfoPs = &fParaPs; + break; + case PositroniumKind::oPs: + pInfoPs = &fOrthoPs; + break; + default: + NoticeError( G4String( __FUNCTION__ ), "improper chosen positronium kind." ); + break; + }; +} + +GateEmittedGammaInformation* GatePositroniumDecayModel::GetPrimaryParticleInformation( const G4PrimaryParticle* pp, const GateEmittedGammaInformation::GammaKind& gamma_kind ) const +{ + GateEmittedGammaInformation* egi = new GateEmittedGammaInformation(); + + GateEmittedGammaInformation::SourceKind source_kind = GateEmittedGammaInformation::SourceKind::ParaPositronium; + GateEmittedGammaInformation::DecayModel decay_model = GateEmittedGammaInformation::DecayModel::Standard; + + if ( fPositroniumKind == PositroniumKind::oPs ) { source_kind = GateEmittedGammaInformation::SourceKind::OrthoPositronium; } + if ( fDecayModel == DecayModel::WithPrompt ) { decay_model = GateEmittedGammaInformation::DecayModel::Deexcitation; } + + egi->SetSourceKind( source_kind ); + egi->SetDecayModel( decay_model ); + egi->SetGammaKind( gamma_kind ); + egi->SetInitialPolarization( pp->GetPolarization() ); + + if ( gamma_kind == GateEmittedGammaInformation::GammaKind::Annihilation ){ egi->SetTimeShift( pInfoPs->GetLifeTime() ); } + + return egi; +} + +G4int GatePositroniumDecayModel::GeneratePrimaryVertices(G4Event* event, G4double& particle_time, G4ThreeVector& particle_position ) +{ + PreparePositroniumParametrization(); + G4int vertexes_number = 1; + + if ( fDecayModel == DecayModel::WithPrompt ) + { + ++vertexes_number; + event->AddPrimaryVertex( GetPrimaryVertexFromDeexcitation(particle_time, particle_position) ); + } + + event->AddPrimaryVertex( GetPrimaryVertexFromPositroniumAnnihilation(particle_time, particle_position) ); + + //Do testu + /*G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position, particle_time); + std::vector gammas_ps = GetGammasFromPositroniumAnnihilation(); + vertex->SetPrimary( GetGammaFromDeexcitation() ); + std::for_each( gammas_ps.begin(), gammas_ps.end(), [&]( G4PrimaryParticle* gamma ) { vertex->SetPrimary( gamma ); } );*/ + + return vertexes_number; +} + +G4PrimaryVertex* GatePositroniumDecayModel::GetPrimaryVertexFromDeexcitation(const G4double& particle_time, const G4ThreeVector& particle_position ) +{ + G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position, particle_time); + vertex->SetPrimary( GetGammaFromDeexcitation() ); + return vertex; +} + +G4PrimaryVertex* GatePositroniumDecayModel::GetPrimaryVertexFromPositroniumAnnihilation( const G4double& particle_time, const G4ThreeVector& particle_position ) +{ + G4double shifted_particle_time = particle_time + G4RandExponential::shoot( pInfoPs->GetLifeTime() ); + + G4PrimaryVertex* vertex = new G4PrimaryVertex( particle_position, shifted_particle_time ); + std::vector gammas = GetGammasFromPositroniumAnnihilation(); + std::for_each( gammas.begin(), gammas.end(), [&]( G4PrimaryParticle* gamma ) { vertex->SetPrimary( gamma ); } ); + return vertex; +} + +G4PrimaryParticle* GatePositroniumDecayModel::GetGammaFromDeexcitation() +{ + G4PrimaryParticle* gamma = GetSingleGamma( fPromptGammaEnergy ); + gamma->SetUserInformation( GetPrimaryParticleInformation( gamma, GateEmittedGammaInformation::GammaKind::Prompt ) ); + return gamma; +} + +std::vector GatePositroniumDecayModel::GetGammasFromPositroniumAnnihilation() +{ + std::vector gammas( pInfoPs->GetAnnihilationGammasNumber() ) ; + + G4DecayProducts* decay_products = pInfoPs->GetDecayProducts(); + for ( G4int i = 0; i < pInfoPs->GetAnnihilationGammasNumber(); ++i ) + { + G4PrimaryParticle* gamma = new G4PrimaryParticle( pGammaDefinition ); + + G4DynamicParticle* dynamic_gamma = (*decay_products)[i]; + G4LorentzVector lv = dynamic_gamma->Get4Momentum(); + gamma->Set4Momentum( lv.px(), lv.py(), lv.pz(), lv.e() ); + gamma->SetPolarization( dynamic_gamma->GetPolarization() ); + gamma->SetUserInformation( GetPrimaryParticleInformation( gamma, GateEmittedGammaInformation::GammaKind::Annihilation ) ); + gammas[i] = gamma; + } + delete decay_products; + + return gammas; +} +} + diff --git a/tests/TestingTools.h b/tests/TestingTools.h new file mode 100644 index 000000000..652f77e5e --- /dev/null +++ b/tests/TestingTools.h @@ -0,0 +1,29 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +/** Authors: Wojciech Krzemień, Mateusz Bała and Kamil Dulski + * Emails: wojciech.krzemien@ncbj.gov.pl, mateusz.bala@ncbj.gov.pl and kamil.dulski@gmail.com + * Organization: National Centre For Nuclear Research (NCBJ, https://ncbj.gov.pl), Poland + * Developed within the IMPET project: https://pet.ncbj.gov.pl/ + * About class: Testing utility header providing the CHECK macro for assertion-style test failure reporting with function name and line number; used by all positronium unit tests. + **/ + +#ifndef TestingTools_h +#define TestingTools_h + + +// ---- Test helper ---- +#define CHECK(cond, msg) \ + do { \ + if (!(cond)) { \ + std::cerr << "Test failure: " << msg \ + << " (in " << __FUNCTION__ << ", line " << __LINE__ << ")" << "\n"; \ + return false; \ + } \ + } while(0) + +#endif diff --git a/tests/test_GateExtendedVSource.cpp b/tests/test_GateExtendedVSource.cpp new file mode 100644 index 000000000..fcbb755c5 --- /dev/null +++ b/tests/test_GateExtendedVSource.cpp @@ -0,0 +1,63 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include +#include + +#include + +#include + +#include "GateRunManager.hh" +#include "GatePhysicsList.hh" +#include "GateDetectorConstruction.hh" + +void initializeGateRunManager(GateRunManager* runManager) +{ + // Set the DetectorConstruction + GateDetectorConstruction* gateDC = new GateDetectorConstruction(); + runManager->SetUserInitialization( gateDC ); + // Set the PhysicsList + runManager->SetUserInitialization( GatePhysicsList::GetInstance() ); + // Initialize G4 kernel + runManager->InitializeAll(); +} + +bool run_tests() +{ + GateExtendedVSource source("source1"); + if(source.GetName() !="source1") { + return false; + } + source.SetType("Ps"); + if(source.GetType() !="Ps") { + return false; + } + + G4Event* pEvent = nullptr; + source.GeneratePrimaries(pEvent); + return true; +} + + +int main() +{ + // Initialization + std::unique_ptr runManager(new GateRunManager); + initializeGateRunManager(runManager.get()); + + bool res = true; + res = res & run_tests(); + + if (res) { + std::cout << "All tests have passed" << std::endl; + return 0; + } else { + std::cerr << "Some tests failed" << std::endl; + return -1; + } +} diff --git a/tests/test_GatePositronium.cpp b/tests/test_GatePositronium.cpp new file mode 100644 index 000000000..a20bac170 --- /dev/null +++ b/tests/test_GatePositronium.cpp @@ -0,0 +1,164 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include +#include + +#include + +#include + +#include "GateRunManager.hh" +#include "GatePhysicsList.hh" +#include "GateDetectorConstruction.hh" + +#include "TestingTools.h" + + +void initializeGateRunManager(GateRunManager* runManager) +{ + // Set the DetectorConstruction + GateDetectorConstruction* gateDC = new GateDetectorConstruction(); + runManager->SetUserInitialization( gateDC ); + // Set the PhysicsList + runManager->SetUserInitialization( GatePhysicsList::GetInstance() ); + // Initialize G4 kernel + runManager->InitializeAll(); +} + +bool test_pPs_properties() +{ + std::cout << "test_pPs_properties\n"; + GatePositronium pPs("pPs", 0.1 * ns); + + CHECK(pPs.GetName() == "pPs", "Name mismatch"); + CHECK(std::abs(pPs.GetLifeTime() - 0.1 * ns) < 1e-12, "Lifetime mismatch"); + CHECK(pPs.GetAnnihilationGammasNumber() == 2, "pPs should decay to 2 gammas"); + + return true; +} + +bool test_oPs_properties() +{ + std::cout << "test_oPs_properties\n"; + GatePositronium oPs("oPs", 142.0 * ns); + + CHECK(oPs.GetName() == "oPs", "Name mismatch"); + CHECK(oPs.GetAnnihilationGammasNumber() == 3, "oPs should decay to 3 gammas"); + + return true; +} + +bool test_move_semantics() +{ + std::cout << "test_move_semantics\n"; + + GatePositronium oPs("oPs", 1000); + GatePositronium moved = std::move(oPs); + + CHECK(moved.GetName() == "oPs", "move: name lost"); + CHECK(moved.GetAnnihilationGammasNumber() == 3, "move: wrong gamma count"); + + return true; +} + + +/// todo: add test to check what happens if pPs with 3 ? +bool test_in_vector() +{ + std::vector vect; + vect.push_back(std::move(GatePositronium("oPs", 1000))); + vect.push_back(std::move(GatePositronium("pPs", 0.1))); + vect.push_back(std::move(GatePositronium("oPs", 2000))); + vect.push_back(std::move(GatePositronium("oPs", 5000))); + if(vect[0].GetName()!="oPs") { + return false; + } + if(vect[1].GetName()!="pPs") { + return false; + } + if(vect[2].GetName()!="oPs") { + return false; + } + if(vect[3].GetName()!="oPs") { + return false; + } + + if(vect[0].GetLifeTime()!= 1000) { + return false; + } + if(vect[1].GetLifeTime()!= 0.1) { + return false; + } + if(vect[2].GetLifeTime()!=2000) { + return false; + } + if(vect[3].GetLifeTime()!=5000) { + return false; + } + + if(vect[0].GetAnnihilationGammasNumber()!= 3) { + return false; + } + if(vect[1].GetAnnihilationGammasNumber()!= 2) { + return false; + } + if(vect[2].GetAnnihilationGammasNumber()!= 3) { + return false; + } + if(vect[3].GetAnnihilationGammasNumber()!= 3) { + return false; + } + + return true; +} + +bool test_decay_products_pPs() +{ + std::cout << "test_decay_products_pPs\n"; + GatePositronium pPs("pPs", 0.1 * ns); + + auto* products = pPs.GetDecayProducts(); + CHECK(products != nullptr, "DecayProducts should not be null"); + CHECK(products->entries() == 2, "pPs should produce 2 daughters"); + + return true; +} + +bool test_decay_products_oPs() +{ + std::cout << "test_decay_products_oPs\n"; + GatePositronium pPs("oPs", 0.1 * ns); + + auto* products = pPs.GetDecayProducts(); + CHECK(products != nullptr, "DecayProducts should not be null"); + CHECK(products->entries() == 3, "oPs should produce 3 daughters"); + return true; +} + +int main() +{ + // Initialization + std::unique_ptr runManager(new GateRunManager); + initializeGateRunManager(runManager.get()); + + bool ok = true; + ok = ok & test_pPs_properties(); + ok = ok & test_oPs_properties(); + ok = ok & test_move_semantics(); + ok = ok & test_in_vector(); + ok = ok & test_decay_products_pPs(); + ok = ok & test_decay_products_oPs(); + + if (ok) { + std::cout << "All tests have passed" << std::endl; + return 0; + } else { + std::cerr << "Some tests failed" << std::endl; + return -1; + } +} diff --git a/tests/test_GatePositroniumDecayModel.cpp b/tests/test_GatePositroniumDecayModel.cpp new file mode 100644 index 000000000..76be4cda5 --- /dev/null +++ b/tests/test_GatePositroniumDecayModel.cpp @@ -0,0 +1,169 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include +#include +#include +#include +#include +#include + +#include + +#include "GateRunManager.hh" +#include "GatePhysicsList.hh" +#include "GateDetectorConstruction.hh" + +#include "TestingTools.h" + +void initializeGateRunManager(GateRunManager* runManager) +{ + // Set the DetectorConstruction + GateDetectorConstruction* gateDC = new GateDetectorConstruction(); + runManager->SetUserInitialization( gateDC ); + // Set the PhysicsList + runManager->SetUserInitialization( GatePhysicsList::GetInstance() ); + // Initialize G4 kernel + runManager->InitializeAll(); +} + +bool run_tests2() +{ + std::unique_ptr runManager(new GateRunManager); + initializeGateRunManager(runManager.get()); + + PositroniumDecayModelParams params; + params.fFractions={0.3,0.7}; + params.fLifetimes={0.1244 ,138.6}; + params.fDecayKind={PositroniumDecayKind::k2Gamma, PositroniumDecayKind::k3Gamma}; + GatePositroniumDecayModel model(params); + + return true; +} + +/// It should be separated in several subtests +bool run_tests() +{ + bool res = true; + PositroniumDecayModelParams params; + GatePositroniumDecayModel model(params); + int index = GatePositroniumDecayModel::getPositroniumDecayIndex({1}); + if (index != 0) { + res = false; + std::cerr << "getPositroniumDecayDecayIndex({1})!=0" << std::endl; + } + index = GatePositroniumDecayModel::getPositroniumDecayIndex({}); + if (index != -1) { + res = false; + std::cerr << "getPositroniumDecayDecayIndex({})!=-1" << std::endl; + } + + float epsilon = 1e-3; + int num_of_trials = 1e6; + std::vector fractions={0.5,0.4,0.1}; + std::vector indices = {0,0,0}; + for (int i = 0; i < num_of_trials; i++) { + index = GatePositroniumDecayModel::getPositroniumDecayIndex(fractions); + indices[index] = indices[index] +1; + } + std::vector estimated_fractions = {0.,0., 0.}; + for (int i = 0; i < estimated_fractions.size(); i++) { + estimated_fractions[i] = float(indices[i])/num_of_trials; + } + + for (int i = 0; i < estimated_fractions.size(); i++) { + if(std::abs(fractions[i] - estimated_fractions[i])>epsilon) + { + res = false; + std::cerr << "assumed fractions and estimated fractions differ:" << std::abs(fractions[i] - estimated_fractions[i]) << std::endl; + } + } + return res; +} + +bool run_tests_positron_range() { + bool res = true; + G4ThreeVector vec(0, 0, 0); + G4double expectedRange = 0.1 * mm; + double epsilon = 3e-3 * mm; + G4int ntrials = 1000000; + std::vector shiftedPos; + std::vector shiftedRadii(ntrials); + for (int i = 0; i < ntrials; i++) { + shiftedPos.push_back( + GatePositroniumDecayModel::AddPositronRangeShift(vec, expectedRange)); + } + std::transform(shiftedPos.begin(), shiftedPos.end(), shiftedRadii.begin(), + [](const G4ThreeVector &vec) { return vec.mag(); }); + auto radiiSum = + std::accumulate(shiftedRadii.begin(), shiftedRadii.end(), 0.0); + auto radiusMean = radiiSum / ntrials; + if (std::abs(radiusMean - expectedRange) > epsilon) { + res = false; + std::cerr << "assumed range and estimated range differ, expected:" + << expectedRange << ", determined:" << radiusMean << std::endl; + } + return res; +} + +// Bug #2: getPositroniumDecayIndex returns -1 when fractions do not sum to 1.0. +// The caller (GeneratePrimaryVertices) uses the return value directly as a vector +// index on fElectronCaptureProbabilities, fMeanPositronRangeEnabled, +// fMeanPositronRange, and fPositroniumDecayChannel (lines 55, 76, 96, 112 of +// GatePositroniumDecayModel.cc), so -1 causes undefined behaviour. + +bool test_decay_index_in_bounds_for_sub_unity_fractions() +{ + // Fractions {0.3, 0.3} sum to 0.6: roughly 40 % of calls return -1. + // The test checks that the return value is always a valid index; + // it fails reliably under the current implementation. + const std::vector fractions = {0.3f, 0.3f}; + for (int i = 0; i < 10000; ++i) { + int index = GatePositroniumDecayModel::getPositroniumDecayIndex(fractions); + CHECK(index >= 0 && index < static_cast(fractions.size()), + "getPositroniumDecayIndex returned " + std::to_string(index) + + " for a vector of size " + std::to_string(fractions.size()) + + " — out-of-bounds index causes UB in GeneratePrimaryVertices (Bug #2)"); + } + return true; +} + +bool test_decay_index_float_fractions_sum_below_unity() +{ + // Even fractions that appear normalized can fail: in IEEE 754 single precision + // (1.0f/3.0f)*3 = 1 - 2^-24 < 1.0, leaving a gap that G4UniformRand (a double) + // can land in, returning -1 from getPositroniumDecayIndex. + // This test confirms the precondition: the cumulative float sum of three equal + // thirds is strictly below 1.0, so the gap exists. + const std::vector fractions = {1.0f/3.0f, 1.0f/3.0f, 1.0f/3.0f}; + float cumulative = 0.0f; + for (float f : fractions) cumulative += f; + CHECK(cumulative == 1.0f, + "Float fractions do not sum to exactly 1.0 (actual sum = " + + std::to_string(cumulative) + + "). Any G4UniformRand() value in (sum, 1.0) causes getPositroniumDecayIndex" + " to return -1 instead of the last valid index (Bug #2)"); + return true; +} + +int main() +{ + bool res = true; + res = res & run_tests(); + res = res & run_tests2(); + res = res & run_tests_positron_range(); + res = res & test_decay_index_in_bounds_for_sub_unity_fractions(); + res = res & test_decay_index_float_fractions_sum_below_unity(); + + if (res) { + std::cout << "All tests have passed" << std::endl; + return 0; + } else { + std::cerr << "Some tests failed" << std::endl; + return -1; + } +} diff --git a/tests/test_GatePositroniumDecayParamsGenerator.cpp b/tests/test_GatePositroniumDecayParamsGenerator.cpp new file mode 100644 index 000000000..f6c5f9e95 --- /dev/null +++ b/tests/test_GatePositroniumDecayParamsGenerator.cpp @@ -0,0 +1,120 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include +#include +#include + +#include "GatePositroniumDecayParamsGenerator.hh" +#include "GatePositroniumDecayModelParams.hh" +#include "GateRunManager.hh" +#include "GatePhysicsList.hh" +#include "GateDetectorConstruction.hh" + +#include "TestingTools.h" + +void initializeGateRunManager(GateRunManager* runManager) +{ + GateDetectorConstruction* gateDC = new GateDetectorConstruction(); + runManager->SetUserInitialization(gateDC); + runManager->SetUserInitialization(GatePhysicsList::GetInstance()); + runManager->InitializeAll(); +} + + + +bool test_positronium_custom() +{ + std::cout << "test_positronium_custom" << std::endl; + GatePositroniumDecayParamsGenerator gen; + + gen.SetPositroniumFraction({0.3f, 0.7f}); + gen.SetPositroniumLifetimes({0.12f, 140.0f}); + gen.SetDecayKinds({k2Gamma, k3Gamma}); + gen.SetPromptGammaProbabilities({0, 1.0}); + gen.SetPromptGammaEnergies({0.0f, 1.2f}); + gen.SetElectronCaptureProbabilities({0.0, 0.5}); + + auto p = gen.generatePositroniumDecayParams(); + + CHECK(p.fFractions[0] == 0.3f, "PS fraction mismatch"); + CHECK(p.fLifetimes[0] == 0.12f, "PS lifetime mismatch"); + CHECK(p.fDecayKind[0] == k2Gamma, "PS decay kind mismatch"); + CHECK(p.fPromptGammaProbabilities[0] == 0.0, "PS prompt probability mismatch"); + CHECK(p.fPromptGammaEnergy[0] == 0.0f, "PS prompt energy mismatch"); + CHECK(p.fElectronCaptureProbabilities[0] == 0.0f, "PS electron capture mismatch"); + + CHECK(p.fFractions[1] == 0.7f, "PS fraction mismatch"); + CHECK(p.fLifetimes[1] == 140.0f, "PS lifetime mismatch"); + CHECK(p.fDecayKind[1] == k3Gamma, "PS decay kind mismatch"); + CHECK(p.fPromptGammaProbabilities[1] == 1.0, "PS prompt probability mismatch"); + CHECK(p.fPromptGammaEnergy[1] == 1.2f, "PS prompt energy mismatch"); + CHECK(p.fElectronCaptureProbabilities[1] == 0.5f, "PS electron capture mismatch"); + return true; +} + +bool test_missing_params_should_fail() +// This test cannot run — GateError calls exit(-1) +{ + GatePositroniumDecayParamsGenerator gen; + std::cout << "test_missing_params_should_fail" << std::endl; + std::cout << "but we cannot run it because GateError() calls exit(-1)" << std::endl; + //bool caught = false; + //try { + //auto p = gen.generatePositroniumDecayParams(GatePositroniumDecayParamsGenerator::kPositronium); + //} catch (...) { + //std::cout << "we caught the exception" << std::endl; + //caught = true; + //} + //CHECK(caught, "Missing parameters did NOT throw"); + return true; +} + + +bool test_vector_size_mismatch() +// This test cannot run — GateError calls exit(-1) +{ + std::cout << "test_vector_size_mismatch" << std::endl; + std::cout << "but we cannot run it because GateError() calls exit(-1)" << std::endl; + + //GatePositroniumDecayParamsGenerator gen; + //gen.SetPositroniumFraction({0.5f, 0.5f}); + //gen.SetPositroniumLifetimes({0.12f}); // mismatch! + //gen.SetDecayKinds({k2Gamma, k3Gamma}); + //gen.SetPromptGammaProbabilities({0, 0}); + //gen.SetPromptGammaEnergies({0.0f, 0.0f}); + + //bool caught = false; + //try { + //auto p = gen.generatePositroniumDecayParams(); + //} catch (...) { + //std::cout << "we caught the exception" << std::endl; + //caught = true; + //} + //CHECK(caught, "Mismatch in vector sizes did NOT throw"); + return true; +} + +int main() +{ + std::unique_ptr runManager(new GateRunManager); + initializeGateRunManager(runManager.get()); + + bool res = true; + res &= test_positronium_custom(); + res &= test_missing_params_should_fail(); + res &= test_vector_size_mismatch(); + + if (res) { + std::cout << " All PositroniumDecayParamsGenerator tests passed\n"; + return 0; + } else { + std::cerr << " Some PositroniumDecayParamsGenerator tests failed\n"; + return -1; + } +} + diff --git a/tests/test_GatePositroniumHelper.cpp b/tests/test_GatePositroniumHelper.cpp new file mode 100644 index 000000000..4cc470afe --- /dev/null +++ b/tests/test_GatePositroniumHelper.cpp @@ -0,0 +1,89 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include +#include + +#include + +#include +#include +#include + +#include "GateRunManager.hh" +#include "GatePhysicsList.hh" +#include "GateDetectorConstruction.hh" + +void initializeGateRunManager(GateRunManager* runManager) +{ + // Set the DetectorConstruction + GateDetectorConstruction* gateDC = new GateDetectorConstruction(); + runManager->SetUserInitialization( gateDC ); + // Set the PhysicsList + runManager->SetUserInitialization( GatePhysicsList::GetInstance() ); + // Initialize G4 kernel + runManager->InitializeAll(); +} + +bool run_tests() +{ + bool res = true; + std::unique_ptr runManager(new GateRunManager); + initializeGateRunManager(runManager.get()); + + GatePositroniumHelper helper; + + std::vector fractions = {2, 12, 4, 2}; + std::vector lifetimes = {0.125, 0.4, 2.0, 50.0}; + std::vector decays = {kParaPs, kDirect, kOrthoPs, kOrthoPs}; + std::vector normFractions = helper.NormalizeFractions(fractions); + std::vector goodFractions = {0.1, 0.6, 0.2, 0.1}; + for (int i=0; i +#include + +#include "GatePositroniumSourceMessenger.hh" +#include "GatePositroniumSource.hh" +#include "GateRunManager.hh" +#include "GateDetectorConstruction.hh" +#include "GatePhysicsList.hh" + +#include "TestingTools.h" + +class DummyPositroniumSource : public GatePositroniumSource { +public: + DummyPositroniumSource() : GatePositroniumSource("Ps") {} +}; + +class TestablePositroniumMessenger : public GatePositroniumSourceMessenger { +public: + using GatePositroniumSourceMessenger::GatePositroniumSourceMessenger; + + G4UIcommand* CmdFractions() { return upCmdSetPositroniumFractions.get(); } + G4UIcommand* CmdLifetimes() { return upCmdSetPositroniumLifetimes.get(); } + G4UIcommand* CmdDecayKinds() { return upCmdSetDecayKinds.get(); } + G4UIcommand* CmdPromptProb() { return upCmdSetPromptPhotonProbabilites.get(); } + G4UIcommand* CmdElectronCaptureProb() { return upCmdSetElectronCaptureProbabilities.get(); } + G4UIcommand* CmdPromptEnergy() { return upCmdSetPromptPhotonEnergies.get(); } + G4UIcommand* CmdInteractions() { return upCmdSetPositronInteractions.get(); } + G4UIcommand* CmdMeanPositronRange() { return upCmdSetMeanPositronRange.get(); } +}; + +void initializeGateRunManager(GateRunManager* runManager) +{ + GateDetectorConstruction* gateDC = new GateDetectorConstruction(); + runManager->SetUserInitialization(gateDC); + runManager->SetUserInitialization(GatePhysicsList::GetInstance()); + runManager->InitializeAll(); +} + +bool test_full_custom() +{ + std::cout << "test_full_custom\n"; + + DummyPositroniumSource src; + TestablePositroniumMessenger msg(&src); + + msg.SetNewValue(msg.CmdFractions(), "0.3 0.7"); + msg.SetNewValue(msg.CmdLifetimes(), "0.1 100 ns"); + msg.SetNewValue(msg.CmdDecayKinds(), "k2Gamma k3Gamma"); + msg.SetNewValue(msg.CmdPromptProb(), "0.2 0.8"); + msg.SetNewValue(msg.CmdPromptEnergy(), "0.5 1.0 MeV"); + msg.SetNewValue(msg.CmdElectronCaptureProb(), "0.3 0.0"); + + auto p = msg.generatePositroniumDecayParams(); + + CHECK(p.fMeanPositronRangeEnabled.size() == 2, "wrong size of fMeanPositronRangeEnabled"); + CHECK(p.fMeanPositronRange.size() == 2, "wrong size of fMeanPositronRange"); + + CHECK(p.fMeanPositronRangeEnabled[0] == false, "positron range enable wrong"); + CHECK(p.fMeanPositronRange[0] == 0.0, "positron range enable wrong"); + + CHECK(p.fFractions.size() == 2, "wrong size"); + CHECK(p.fFractions[1] == 0.7f, "fraction wrong"); + CHECK(p.fDecayKind[1] == k3Gamma, "decay kind wrong"); + CHECK(p.fPromptGammaEnergy[1] == 1.0f, "energy wrong"); + CHECK(p.fMeanPositronRangeEnabled[1] == false, "positron range enable wrong"); + CHECK(p.fMeanPositronRange[1] == 0.0, "positron range enable wrong"); + + return true; +} + + +bool test_unit_handling_in_commands() +{ + std::cout << "test_unit_handling_in_commands\n"; + + DummyPositroniumSource src; + TestablePositroniumMessenger msg(&src); + + msg.SetNewValue(msg.CmdFractions(), "0.3 0.7"); + msg.SetNewValue(msg.CmdLifetimes(), "0.1 100 ns"); + msg.SetNewValue(msg.CmdDecayKinds(), "k2Gamma k3Gamma"); + msg.SetNewValue(msg.CmdPromptProb(), "0.2 0.8"); + msg.SetNewValue(msg.CmdPromptEnergy(), "0.5 1.0 keV"); + msg.SetNewValue(msg.CmdMeanPositronRange(), "0.1 0.0 cm"); + msg.SetNewValue(msg.CmdElectronCaptureProb(), "0.0 0.0"); + + auto p = msg.generatePositroniumDecayParams(); + + CHECK(p.fMeanPositronRangeEnabled[0] == true, "positron range enabled wrong"); + CHECK(p.fMeanPositronRange[0] == 0.1 * 10, "positron range wrong"); // Cause mm is the default unit + + CHECK(p.fFractions.size() == 2, "wrong size"); + CHECK(p.fFractions[1] == 0.7f, "fraction wrong"); + CHECK(p.fDecayKind[1] == k3Gamma, "decay kind wrong"); + CHECK(p.fPromptGammaEnergy[1] == 1.0f /1000, "energy wrong"); // Cause we used keV and MeV is the default unit + CHECK(p.fMeanPositronRangeEnabled[1] == true, "positron range enabled wrong"); + CHECK(p.fMeanPositronRange[1] == 0.0 *10 , "positron range wrong"); + + return true; +} +// ========================================================== +// main() +// ========================================================== +int main() +{ + std::unique_ptr runManager(new GateRunManager); + initializeGateRunManager(runManager.get()); + + bool res = true; + + res &= test_full_custom(); + res &= test_unit_handling_in_commands(); + + if (res) { + std::cout << " All GatePositroniumSourceMessenger tests passed\n"; + return 0; + } else { + std::cerr << " Some GatePositroniumSourceMessenger tests failed\n"; + return -1; + } +} + + diff --git a/tests/test_dummy.cpp b/tests/test_dummy.cpp new file mode 100644 index 000000000..d8c58cf31 --- /dev/null +++ b/tests/test_dummy.cpp @@ -0,0 +1,29 @@ +/** ---------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + ----------------------*/ + +#include +#include +bool run_tests() +{ + if (1==0) { + return true; + } + return false; +} + +int main() +{ + bool res = true; + res = res & run_tests(); + if (res) { + std::cout << "All tests have passed" << std::endl; + return 0; + } else { + std::cerr << "Some tests failed" << std::endl; + return -1; + } +}