Adding handling of decayIndex and other variables used by GatePositroniumSource in the context of Singles and Coincidences

docs/data_output_management.rst

@@ -196,6 +196,67 @@ generation (Poisson law). Gate generates the time delay from the
 previous event, if it is out of the time slice it stops the event 
 processing for the current time slice and if needed it starts a new time slice.
 
+.. _analysis_output_modules-label:
+
+Analysis output modules
+-----------------------
+
+A sensitive detector records the energy deposit, the position and the process name of every hit, but not the Monte Carlo information describing the *history* of the detected photon: the event and source it belongs to, how many times it scattered in the phantom and in the crystal before being absorbed, and in which volume the last scattering took place.
+Those attributes are filled at the end of each event by an **analysis output module**, which then triggers the digitizer, so every Single and every Coincidence inherits them.
+
+Three modules implement that step and they are mutually exclusive: **exactly one has to be enabled** whenever Singles or Coincidences are written. With none of them enabled GATE stops the simulation with::
+
+   ***ERROR*** Digitizer Manager is not initialized properly. Please, enable analysis,
+   fastanalysis or multianalysis Output Modules to write down Singles or Coincidences.
+
+.. list-table:: Analysis output modules
+   :widths: 22 30 48
+   :header-rows: 1
+   :name: analysis_output_modules_table
+
+   * - Module
+     - Command
+     - Purpose
+   * - ``analysis``
+     - ``/gate/output/analysis/enable`` (enabled by default)
+     - Reference implementation. Fills every attribute, but assigns the interaction counters only to the two annihilation gammas of an event.
+   * - ``fastanalysis``
+     - ``/gate/output/fastanalysis/enable``
+     - Speed-oriented variant. Fills only ``runID``, ``eventID`` and ``sourceID``; every other attribute is set to ``-1`` or ``"NULL"`` to mark it as not computed.
+   * - ``multianalysis``
+     - ``/gate/output/multianalysis/enable``
+     - Multi-photon variant, see below. Fills the same attributes as ``analysis`` plus ``nInteractions``, and does so for every emitted gamma.
+
+Remember to disable the module you are replacing, since ``analysis`` is enabled by default::
+
+   /gate/output/analysis/disable
+   /gate/output/fastanalysis/disable
+   /gate/output/multianalysis/enable
+
+The multi-photon analysis
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``GateMultiPhotonAnalysis`` targets sources that emit more than two gammas per decay - ortho-positronium, which decays into three gammas, and de-excitation (prompt) gammas emitted alongside the annihilation pair. ``GateAnalysis`` resolves exactly two annihilation gammas per event, so hits belonging to a third gamma or to a prompt gamma keep zeroed counters; ``multianalysis`` builds the set of reference photons from the event itself and follows each of them separately.
+
+For decays into two gammas both modules agree: ``eventID``, ``trackID``, ``parentID``, ``processName``, ``nPhantomCompton``, ``nCrystalCompton``, ``nPhantomRayleigh``, ``nCrystalRayleigh``, ``comptVolName`` and ``RayleighVolName`` are filled identically.
+
+Three differences are deliberate:
+
+* ``nInteractions`` is filled **only** by ``multianalysis``. It counts the Compton and Rayleigh scatterings along the path of the photon, the current hit included, and grows along the whole history of that photon. On the ``analysis`` and ``fastanalysis`` paths the branch stays at ``-1``, which means "not computed".
+* ``photonID`` is always ``0``. In ``GateAnalysis`` the field distinguishes the first from the second annihilation gamma, which has no meaning once an event may contain three gammas or a prompt gamma; ``trackID`` identifies the photon instead.
+* septal penetration is configured through the ``analysis`` module even when it is disabled. ``/gate/output/analysis/setSeptalVolumeName`` and ``/gate/output/analysis/recordSeptalPenetration`` keep working, and ``multianalysis`` reads that configuration to fill ``septalNb`` exactly like ``GateAnalysis`` does.
+
+A typical macro looks like this::
+
+   /gate/output/analysis/disable
+   /gate/output/fastanalysis/disable
+   /gate/output/multianalysis/enable
+
+   /gate/output/root/enable
+   /gate/output/root/setFileName data
+   /gate/output/root/setRootHitFlag     1
+   /gate/output/root/setRootSinglesFlag 1
+
 .. _root_output-label:
 
 Root output
@@ -232,6 +293,46 @@ If you want to disable the whole ROOT output, just do not call it, or use the fo
    /gate/output/root/disable
 
 
+.. _decay_branches_root_output-label:
+
+Branches describing the decay and the interaction counters
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``PositroniumSource`` attaches four integer fields to every emitted gamma: ``sourceType``, ``decayType``, ``gammaType`` and ``decayIndex``. Their values are listed in :ref:`positronium_source_branches-label`. A fifth field, ``nInteractions``, counts the Compton and Rayleigh scatterings along the path of the photon and is filled by the ``multianalysis`` output module (``GateMultiPhotonAnalysis``) only.
+
+All five fields describe the decay and the photon rather than the individual energy deposit, so they survive digitization and are written to the three ROOT trees:
+
+.. list-table:: Branches propagated to the ROOT trees
+   :widths: 20 40 40
+   :header-rows: 1
+   :name: decay_branches_root_trees
+
+   * - Tree
+     - Branch names
+     - Source of the value
+   * - ``Hits``
+     - ``sourceType``, ``decayType``, ``gammaType``, ``decayIndex``, ``nInteractions``
+     - the hit itself
+   * - ``Singles``
+     - same names
+     - copied from the hit, then merged by the digitizer
+   * - ``Coincidences``
+     - same names with the ``1`` and ``2`` suffix, one per arm
+     - the two digis forming the pair
+
+Merging rules follow the convention already used for ``sourceEnergy`` and ``sourcePDG``:
+
+* ``sourceType`` and ``decayIndex`` describe the decay, and all gammas of one event come from the same decay, so merging cannot produce a conflict;
+* ``decayType`` and ``gammaType`` describe a single gamma. When hits of gammas of different kinds are merged in one crystal - an annihilation gamma and a prompt gamma, for instance - the value falls back to the "not known" value of the enum, that is ``0``;
+* ``decayIndex`` falls back to ``-1`` for the same reason;
+* ``nInteractions`` is merged with the maximum, like ``nPhantomCompton``.
+
+Files written by older versions of GATE do not contain these branches in ``Singles`` and ``Coincidences``; reading such a file keeps the buffer at its "not known" value instead of failing.
+
+.. note::
+
+   The unified tree output (:ref:`new_unified_tree_output_general_set-label`) writes ``sourceType``, ``decayType``, ``gammaType`` and ``decayIndex`` for hits only, and does not write ``nInteractions`` at all. Use the ROOT output described here if you need those fields at the Singles or Coincidences level.
+
 Using TBrowser To Browse ROOT Objects
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

docs/source_and_particle_management.rst

@@ -822,7 +822,7 @@ Required commands
 ``/gate/source/NAME/setPositroniumLifetimes t1 t2 ... tn unit``
    Lifetime of each component. The last token is a Geant4 time unit such as ``ns``.
 
-``/gate/source/NAME/setPromptPhotonProbabilites p1 p2 ... pn``
+``/gate/source/NAME/setPromptPhotonProbabilities p1 p2 ... pn``
    Probability of prompt-gamma emission for each component. Values are restricted to the ``[0, 1]`` range by the messenger.
 
 ``/gate/source/NAME/setPromptPhotonEnergies e1 e2 ... en unit``
@@ -834,11 +834,11 @@ One of the following commands must also be provided:
    Explicit annihilation channel selection for each component. Allowed values are ``k2Gamma`` and ``k3Gamma``.
 
 ``/gate/source/NAME/setPositronInteractions k1 k2 ... kn``
-   Interaction type assigned to each component. Allowed values are ``kParaPs``, ``kDirect`` and ``kOrthoPs``. This command is used by the model to derive the ``sourceType`` metadata stored in the **Hits** tree. It can also be used instead of ``setDecayKinds`` to derive the 2-gamma and 3-gamma contributions from the theory.
+   Interaction type assigned to each component. Allowed values are ``kParaPs``, ``kDirect`` and ``kOrthoPs``. This command is used by the model to derive the ``sourceType`` metadata stored in the **Hits**, **Singles** and **Coincidences** trees. It can also be used instead of ``setDecayKinds`` to derive the 2-gamma and 3-gamma contributions from the theory.
 
 .. note::
 
-   The current implementation can build the decay model from ``setDecayKinds`` alone, but the ``sourceType`` metadata in the **Hits** tree is derived from ``setPositronInteractions``. In practice, if you want a fully consistent ``PositroniumSource`` configuration and meaningful ``sourceType`` values, define ``setPositronInteractions`` as well.
+   The current implementation can build the decay model from ``setDecayKinds`` alone, but the ``sourceType`` metadata written to the output trees is derived from ``setPositronInteractions``. In practice, if you want a fully consistent ``PositroniumSource`` configuration and meaningful ``sourceType`` values, define ``setPositronInteractions`` as well.
 
 Optional commands
 ^^^^^^^^^^^^^^^^^
@@ -855,7 +855,7 @@ Optional commands
 
 .. note::
 
-   The current implementation requires both ``setPromptPhotonProbabilites`` and ``setPromptPhotonEnergies`` even when the prompt-gamma probability is zero for all components. A common convention is to set the corresponding energies to ``0``.
+   The current implementation requires both ``setPromptPhotonProbabilities`` and ``setPromptPhotonEnergies`` even when the prompt-gamma probability is zero for all components. A common convention is to set the corresponding energies to ``0``.
 
 Event timing and generated vertices
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -869,8 +869,10 @@ If ``setMeanPositronRange`` is used, the annihilation vertex position is additio
 
 If the electron-capture probability is very high and the prompt-gamma probability is very low, event generation can become inefficient because the source retries until at least one primary vertex is created. The code emits a warning when the probability of producing neither a prompt gamma nor an annihilation exceeds 90% for a component.
 
-Dedicated Hits branches
-~~~~~~~~~~~~~~~~~~~~~~~
+.. _positronium_source_branches-label:
+
+Dedicated output branches
+~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Hits produced by ``PositroniumSource`` carry four additional integer branches in the **Hits** tree:
 
@@ -879,6 +881,8 @@ Hits produced by ``PositroniumSource`` carry four additional integer branches in
 * **decayType** - decay model attached to the sampled component,
 * **decayIndex** - zero-based index of the sampled decay component.
 
+The same four branches are also written to the **Singles** and **Coincidences** trees, in the latter with the ``1`` and ``2`` suffix for the two arms of the pair; the merging rules applied by the digitizer are described in :ref:`decay_branches_root_output-label`.
+
 The ``decayIndex`` value is copied to all gammas emitted from the same sampled component. In the explicit ``setDecayKinds`` workflow, the index follows the order of components in the user-defined vectors. When the source is configured through ``setPositronInteractions``, the helper may internally expand one user component into several annihilation channels, so ``decayIndex`` refers to the final internal component list used by the model.
 
 .. table:: Description of values in gammaType branch for PositroniumSource
@@ -958,7 +962,7 @@ The example below defines a two-component source with explicit 2-gamma and 3-gam
    /gate/source/psSource/setPositroniumLifetimes 0.125 142.0 ns
    /gate/source/psSource/setDecayKinds k2Gamma k3Gamma
    /gate/source/psSource/setPositronInteractions kParaPs kOrthoPs
-   /gate/source/psSource/setPromptPhotonProbabilites 1.0 1.0
+   /gate/source/psSource/setPromptPhotonProbabilities 1.0 1.0
    /gate/source/psSource/setPromptPhotonEnergies 1.274 1.274 MeV
    /gate/source/psSource/setElectronCaptureProbabilities 0.0 0.0
    /gate/source/psSource/setMeanPositronRange 0.2 0.2 mm

source/digits_hits/include/GateCrystalSD.hh

@@ -88,6 +88,12 @@ class GateCrystalSD : public G4VSensitiveDetector
 
       G4int PrepareCreatorAttachment(GateVVolume* aCreator);
 
+      //! Tells whether a crystal SD may be attached to this volume, i.e. whether the volume
+      //! belongs to a system, and explains on G4cout when it may not. Has to be asked
+      //! *before* the SD is built: the constructor registers the SD in GateDigitizerMgr and
+      //! nothing takes it back, so a refused attachment would still leave a collection behind.
+      static G4bool CanAttachToCreator(GateVVolume* aCreator);
+
   protected:
      GateVSystem* m_system;                           //! System to which the SD is attached //mhadi_obso obsollete, because we use the multi-system approach
      GateSystemList* m_systemList = nullptr;          //! System list instead of one system

source/digits_hits/include/GateDigi.hh

@@ -142,6 +142,29 @@ public:
       inline G4int GetNSeptal() const { return m_nSeptal; }
       inline void SetNSeptal(G4int septalNb) { m_nSeptal = septalNb; }
 
+      //! \name Information coming from PositroniumSource, copied from the hit
+      //! The enums are defined in GateEmittedGammaInformation.hh. Merging keeps the value
+      //! when the merged digis agree and falls back to the "not known" value when they do
+      //! not, exactly like m_sourceEnergy / m_sourcePDG.
+      //@{
+      inline void  SetSourceType(G4int sourceType) { m_sourceType = sourceType; }
+      inline G4int GetSourceType() const { return m_sourceType; }
+
+      inline void  SetDecayType(G4int decayType) { m_decayType = decayType; }
+      inline G4int GetDecayType() const { return m_decayType; }
+
+      inline void  SetGammaType(G4int gammaType) { m_gammaType = gammaType; }
+      inline G4int GetGammaType() const { return m_gammaType; }
+
+      inline void  SetDecayIndex(G4int decayIndex) { m_decayIndex = decayIndex; }
+      inline G4int GetDecayIndex() const { return m_decayIndex; }
+      //@}
+
+      //! Number of Compton and Rayleigh interactions along the photon path; filled by
+      //! GateMultiPhotonAnalysis only, so it stays at -1 on the classical GateAnalysis path
+      inline void  SetNInteractions(G4int nInteractions) { m_nInteractions = nInteractions; }
+      inline G4int GetNInteractions() const { return m_nInteractions; }
+
 
 
 
@@ -239,6 +262,13 @@ public:
   #endif
   G4int m_nSeptal;				  //!< HDS : record septal penetration
 
+  // Information from PositroniumSource, copied from the hit (see GateEmittedGammaInformation.hh)
+  G4int m_sourceType;             //!< kind of gamma source; 0 means 'not known'
+  G4int m_decayType;              //!< kind of positronium decay; 0 means 'not known'
+  G4int m_gammaType;              //!< kind of gamma: single, annihilation, prompt; 0 means 'not known'
+  G4int m_decayIndex;             //!< decay channel index; -1 means 'not known'
+  G4int m_nInteractions;          //!< # of Compton and Rayleigh interactions; -1 when not computed
+
   // AE : Added for IdealComptonPhot adder which take into account several Comptons in the same volume
   //These variables no sense for a general pulse but I need them to  process idealy the hits. or create another structure
   G4String m_Postprocess;         // PostStep process