[WIP] Multi-Channel Positronium Decay Source Model

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)

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; }
 };
 
 

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 <vector>
+
+/** 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

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 <string>
+#include <string_view>
+
+/** 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