diff --git a/docs/data_output_management.rst b/docs/data_output_management.rst index 9f9202986..16fb768f7 100644 --- a/docs/data_output_management.rst +++ b/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 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/source_and_particle_management.rst b/docs/source_and_particle_management.rst index 35b2b45e1..747007df6 100644 --- a/docs/source_and_particle_management.rst +++ b/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 diff --git a/source/digits_hits/include/GateCrystalSD.hh b/source/digits_hits/include/GateCrystalSD.hh index 0354bde18..b6c4f7091 100644 --- a/source/digits_hits/include/GateCrystalSD.hh +++ b/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 diff --git a/source/digits_hits/include/GateDigi.hh b/source/digits_hits/include/GateDigi.hh index 95908f6df..472dc1715 100644 --- a/source/digits_hits/include/GateDigi.hh +++ b/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 diff --git a/source/digits_hits/include/GateMultiPhotonAnalysisHelpers.hh b/source/digits_hits/include/GateMultiPhotonAnalysisHelpers.hh index ef10a140a..c40d4fd83 100644 --- a/source/digits_hits/include/GateMultiPhotonAnalysisHelpers.hh +++ b/source/digits_hits/include/GateMultiPhotonAnalysisHelpers.hh @@ -15,99 +15,137 @@ * 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. + * About namespace: Helper functions for multi-photon analysis, including structures holding + * gamma interaction statistics and functions accumulating them from interaction process names. **/ namespace MultiPhotonAnalysisHelpers { - /** @brief Aggregated interaction counters for a single tracked gamma. */ - struct GammaStatistics { - int phantomCompton = 0; - int phantomRayleigh = 0; + /** @brief Default name stored when no scattering volume has been recorded yet. */ + inline const std::string kNoVolumeName = "NULL"; + + /** + * @brief Phantom scattering counters of a single gamma, summed over the whole event. + * + * GateAnalysis counts phantom interactions in a loop that runs before crystal hits are + * processed, so every crystal hit of a given gamma receives the same event total. These + * counters reproduce that behaviour. + */ + struct PhantomStatistics { + int compton = 0; + int rayleigh = 0; + std::string comptonVolumeName = kNoVolumeName; + std::string rayleighVolumeName = kNoVolumeName; + }; + + /** + * @brief Running counters of a single gamma, accumulated along the event time axis. + * + * Crystal counters reproduce GateAnalysis: they grow hit by hit and include the current + * hit. `scatters` feeds nInteractions and counts Compton and Rayleigh interactions in both + * the phantom and the crystal. + */ + struct RunningStatistics { int crystalCompton = 0; int crystalRayleigh = 0; - int phantomInteractions = 0; - int crystalInteractions = 0; - std::string comptonVolumeName = "NULL"; - std::string rayleighVolumeName = "NULL"; + int scatters = 0; }; /** @brief Interaction process classes used by multi-photon analysis. */ enum class InteractionProcess { - Unknown, + Other, Compton, - Rayleigh, - Transportation + Rayleigh }; /** * @brief Classifies process name into a compact process category. * + * Matching is done by substring, exactly like GateAnalysis (`find("ompt")` and + * `find("Rayl")`), so names such as `compt`, `Compton` or `LowEnCompt` are all recognized. + * * Args: * processName: Geant4 process name. * * Returns: - * Recognized interaction category or Unknown. + * Recognized interaction category or Other. */ 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; - } + if (processName.find("ompt") != std::string_view::npos) { + return InteractionProcess::Compton; } - return InteractionProcess::Unknown; + if (processName.find("Rayl") != std::string_view::npos) { + return InteractionProcess::Rayleigh; + } + return InteractionProcess::Other; } /** - * @brief Accumulates crystal-side interaction statistics. + * @brief Adds a phantom interaction to the event totals of one gamma. * * Args: - * processName: Geant4 process name at crystal hit. - * gammaStatistics: Statistics object to update. + * processName: Geant4 process name at the phantom hit. + * volumeName: Name of the volume the phantom hit occurred in. + * phantomStatistics: 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; - } + inline void AccumulatePhantomTotals(const std::string_view processName, + const std::string& volumeName, + PhantomStatistics& phantomStatistics) { + switch (GetInteractionProcess(processName)) { + case InteractionProcess::Compton: + phantomStatistics.compton++; + phantomStatistics.comptonVolumeName = volumeName; + break; + case InteractionProcess::Rayleigh: + phantomStatistics.rayleigh++; + phantomStatistics.rayleighVolumeName = volumeName; + break; + default: + break; } } /** - * @brief Accumulates phantom-side interaction statistics. + * @brief Accumulates a crystal interaction into the running counters. + * + * Args: + * processName: Geant4 process name at the crystal hit. + * runningStatistics: Statistics object to update. + */ + inline void AccumulateCrystal(const std::string_view processName, + RunningStatistics& runningStatistics) { + switch (GetInteractionProcess(processName)) { + case InteractionProcess::Compton: + runningStatistics.crystalCompton++; + runningStatistics.scatters++; + break; + case InteractionProcess::Rayleigh: + runningStatistics.crystalRayleigh++; + runningStatistics.scatters++; + break; + default: + break; + } + } + + /** + * @brief Accumulates a phantom interaction into the running counters. + * + * Only the nInteractions counter is affected - phantom Compton and Rayleigh columns are + * event totals and are computed separately (see PhantomStatistics). * * Args: - * processName: Geant4 process name at phantom hit. - * gammaStatistics: Statistics object to update. + * processName: Geant4 process name at the phantom hit. + * runningStatistics: 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; - } + inline void AccumulatePhantom(const std::string_view processName, + RunningStatistics& runningStatistics) { + switch (GetInteractionProcess(processName)) { + case InteractionProcess::Compton: + case InteractionProcess::Rayleigh: + runningStatistics.scatters++; + break; + default: + break; } } diff --git a/source/digits_hits/include/GateMultiPhotonTrajectoryNavigator.hh b/source/digits_hits/include/GateMultiPhotonTrajectoryNavigator.hh index bf2d71a41..b244cb987 100644 --- a/source/digits_hits/include/GateMultiPhotonTrajectoryNavigator.hh +++ b/source/digits_hits/include/GateMultiPhotonTrajectoryNavigator.hh @@ -11,6 +11,7 @@ See LICENSE.md for further details #include "G4ThreeVector.hh" +#include #include #include #include @@ -33,6 +34,16 @@ class G4TrajectoryContainer; */ class GateMultiPhotonTrajectoryNavigator { public: + /** + * @brief Maximum number of reference photons kept for a single event. + * + * Four is the largest number GatePositroniumSource can emit from one decay: an ortho + * positronium decay gives three annihilation gammas and the source may additionally emit + * a prompt (de-excitation) gamma. Reference photons above this limit are dropped with + * a warning. If a more complex process is ever simulated, this limit has to be raised. + */ + static constexpr std::size_t kMaxReferencePhotons = 4; + /** * @brief Constructs an empty event-scoped navigator. */ @@ -69,6 +80,17 @@ class GateMultiPhotonTrajectoryNavigator { */ std::vector FindReferencePhotonTrackIDs() const; + /** + * @brief Tells whether a track is itself a reference photon. + * + * Args: + * trackID: Track identifier to check. + * + * Returns: + * True when the track belongs to the reference photon set of the current event. + */ + bool IsReferencePhoton(int trackID) const; + /** * @brief Resolves ancestor reference photon for a track. * diff --git a/source/digits_hits/include/GateVDigitizerModule.hh b/source/digits_hits/include/GateVDigitizerModule.hh index 46c652179..8e19c092f 100755 --- a/source/digits_hits/include/GateVDigitizerModule.hh +++ b/source/digits_hits/include/GateVDigitizerModule.hh @@ -40,6 +40,11 @@ public: GateDigi* CentroidMerge(GateDigi* right, GateDigi* output ); GateDigi* MergePositionEnergyWin(GateDigi *right, GateDigi *output); + //! Merge the fields describing the decay, the kind of gamma and the interaction counter. + //! Shared by every adder so that all of them follow the same rule; see the definition + //! for the rule itself. + static void MergeEmittedGammaInformation(GateDigi *right, GateDigi *output); + //! Method overloading GateClockDependent::Describe() diff --git a/source/digits_hits/src/GateAdderComptPhotIdeal.cc b/source/digits_hits/src/GateAdderComptPhotIdeal.cc index 4c239f3b3..8bd9017fc 100755 --- a/source/digits_hits/src/GateAdderComptPhotIdeal.cc +++ b/source/digits_hits/src/GateAdderComptPhotIdeal.cc @@ -329,6 +329,8 @@ GateDigi* GateAdderComptPhotIdeal::CentroidMergeComptPhotIdeal(GateDigi *right, output->m_nSeptal = right->m_nSeptal; } + MergeEmittedGammaInformation(right, output); + // VolumeID: should be identical for both pulses, we do nothing // m_scannerPos: identical for both pulses, nothing to do // m_scannerRotAngle: identical for both pulses, nothing to do diff --git a/source/digits_hits/src/GateAdderCompton.cc b/source/digits_hits/src/GateAdderCompton.cc index 075294d04..5db610cc3 100755 --- a/source/digits_hits/src/GateAdderCompton.cc +++ b/source/digits_hits/src/GateAdderCompton.cc @@ -226,6 +226,8 @@ GateDigi* GateAdderCompton::CentroidMergeCompton(GateDigi *right, GateDigi *outp output->m_nSeptal = right->m_nSeptal; } + MergeEmittedGammaInformation(right, output); + // VolumeID: should be identical for both pulses, we do nothing // m_scannerPos: identical for both pulses, nothing to do // m_scannerRotAngle: identical for both pulses, nothing to do diff --git a/source/digits_hits/src/GateCoincidenceSorter.cc b/source/digits_hits/src/GateCoincidenceSorter.cc index 4a52c0e0b..5a38bfc1e 100644 --- a/source/digits_hits/src/GateCoincidenceSorter.cc +++ b/source/digits_hits/src/GateCoincidenceSorter.cc @@ -514,10 +514,19 @@ bool GateCoincidenceSorter::ProcessCompletedCoincidenceWindow(GateCoincidenceDig return; } */ - if ( (m_multiplesPolicy==kTakeWinnerIfOnlyOneGood) && (nGoods==1)) + if ( (m_multiplesPolicy==kTakeWinnerIfOnlyOneGood) && (nGoods==1)) { - m_OutputCoincidenceDigiCollection->insert(coincidence); - return true; + // Store the one good pair, not the multiple it was found in. A GateCoincidenceDigi is a + // vector of digis, but the ROOT writer only ever reads GetDigi(0) and GetDigi(1), so an + // unsplit multiple would put the first two digis of the window on the output - a pair + // IsForbiddenCoincidence had just rejected - and drop the remaining singles silently. + // Same loop as in kKillAllIfMultipleGoods below, which already does this correctly. + for(i=0; i<(coincidence->IsDelayed()?1:(nDigis-1)); i++) + for(j=i+1; jat(i),coincidence->at(j))) + m_OutputCoincidenceDigiCollection->insert(CreateSubDigi(coincidence, i, j)); + + return false; } if( (m_multiplesPolicy==kTakeWinnerIfOnlyOneGood) ) { diff --git a/source/digits_hits/src/GateCrystalSD.cc b/source/digits_hits/src/GateCrystalSD.cc index 415dc1bf7..10a7dd462 100644 --- a/source/digits_hits/src/GateCrystalSD.cc +++ b/source/digits_hits/src/GateCrystalSD.cc @@ -362,20 +362,35 @@ G4bool GateCrystalSD::ProcessHits(G4Step*aStep, G4TouchableHistory*) //virtual void EndOfEvent (G4HCofThisEvent *) +//------------------------------------------------------------------------------ +//! Asked before the SD is built, so that a refused attachment leaves nothing behind. +//! The message keeps its original wording and prefix: user macros and logs may match on it. +G4bool GateCrystalSD::CanAttachToCreator(GateVVolume* aCreator) +{ + if (GateSystemListManager::GetInstance()->FindSystemOfCreator(aCreator->GetCreator())) { + return true; + } + + G4cout << Gateendl << Gateendl << "[GateCrystalSD::PrepareCreatorAttachment]:\n" + << "Volume '" << aCreator->GetObjectName() << "' does not belong to any system.\n" + << "Your volume must belong to a system to be used with a crystalSD.\n" + << "Attachment request ignored --> you won't have any hit output from this volume!!!\n" << Gateendl; + return false; +} +//------------------------------------------------------------------------------ + + //------------------------------------------------------------------------------ //! Next method underwent an important modification to be compatible with the multi-system approach G4int GateCrystalSD::PrepareCreatorAttachment(GateVVolume* aCreator) { - GateVSystem* creatorSystem = GateSystemListManager::GetInstance()->FindSystemOfCreator(aCreator->GetCreator()); - if (!creatorSystem) { - G4cout << Gateendl << Gateendl << "[GateCrystalSD::PrepareCreatorAttachment]:\n" - << "Volume '" << aCreator->GetObjectName() << "' does not belong to any system.\n" - << "Your volume must belong to a system to be used with a crystalSD.\n" - << "Attachment request ignored --> you won't have any hit output from this volume!!!\n" << Gateendl; + if (!CanAttachToCreator(aCreator)) { return -1; } + GateVSystem* creatorSystem = GateSystemListManager::GetInstance()->FindSystemOfCreator(aCreator->GetCreator()); + if (!m_systemList) { m_systemList = new GateSystemList; diff --git a/source/digits_hits/src/GateDigi.cc b/source/digits_hits/src/GateDigi.cc index 78a09ad5d..3855fde5f 100644 --- a/source/digits_hits/src/GateDigi.cc +++ b/source/digits_hits/src/GateDigi.cc @@ -39,6 +39,11 @@ GateDigi::GateDigi(const void* itsMother): #ifdef GATE_USE_OPTICAL m_optical(false), #endif + m_sourceType(0), + m_decayType(0), + m_gammaType(0), + m_decayIndex(-1), + m_nInteractions(-1), m_energyError(0.0), m_globalPosError(0.0), m_localPosError(0.0), diff --git a/source/digits_hits/src/GateDigitizerInitializationModule.cc b/source/digits_hits/src/GateDigitizerInitializationModule.cc index 9815536d9..6b2ef2b81 100644 --- a/source/digits_hits/src/GateDigitizerInitializationModule.cc +++ b/source/digits_hits/src/GateDigitizerInitializationModule.cc @@ -133,16 +133,25 @@ void GateDigitizerInitializationModule::Digitize() Digi->SetSourcePDG((*inHC)[i]->GetSourcePDG()); Digi->SetNCrystalConv( (*inHC)[i]->GetNCrystalConv() ); + // Information from ExtendedVSource plus the interaction counter: they describe + // the decay and the photon, not the deposit, so they survive digitization + Digi->SetSourceType( (*inHC)[i]->GetSourceType() ); + Digi->SetDecayType( (*inHC)[i]->GetDecayType() ); + Digi->SetGammaType( (*inHC)[i]->GetGammaType() ); + Digi->SetDecayIndex( (*inHC)[i]->GetDecayIndex() ); + Digi->SetNInteractions( (*inHC)[i]->GetNInteractions() ); + //------------------------------------------------- + // An empty volume name only means that the output module in use did not + // fill it, so it is replaced by the "NULL" placeholder. It tells nothing + // about the source of the event, hence sourceID is left untouched. if ((*inHC)[i]->GetComptonVolumeName().empty()) { Digi->SetComptonVolumeName( "NULL" ); - Digi->SetSourceID( -1 ); } if ((*inHC)[i]->GetRayleighVolumeName().empty()) { Digi->SetRayleighVolumeName( "NULL" ); - Digi->SetSourceID( -1 ); } /* // if (nVerboseLevel>1) diff --git a/source/digits_hits/src/GateHit.cc b/source/digits_hits/src/GateHit.cc index 99d307d83..47183854f 100644 --- a/source/digits_hits/src/GateHit.cc +++ b/source/digits_hits/src/GateHit.cc @@ -27,6 +27,7 @@ GateHit::GateHit() m_PDGEncoding(0), m_trackID(0), m_parentID(0), + m_sourceID(-1), m_systemID(-1), m_sourceEnergy(-1), m_sourcePDG(0), diff --git a/source/digits_hits/src/GateHitConvertor.cc b/source/digits_hits/src/GateHitConvertor.cc index 818a99ef2..f76330d9e 100644 --- a/source/digits_hits/src/GateHitConvertor.cc +++ b/source/digits_hits/src/GateHitConvertor.cc @@ -179,14 +179,15 @@ void GateHitConvertor::ProcessOneHit(const GateHit* hit,GatePulseList* pulseList + // An empty volume name only means that the output module in use did not fill it, + // so it is replaced by the "NULL" placeholder. It tells nothing about the source + // of the event, hence sourceID is left untouched. if (hit->GetComptonVolumeName().empty()) { pulse->SetComptonVolumeName( "NULL" ); - pulse->SetSourceID( -1 ); } if (hit->GetRayleighVolumeName().empty()) { pulse->SetRayleighVolumeName( "NULL" ); - pulse->SetSourceID( -1 ); } if (nVerboseLevel>1) diff --git a/source/digits_hits/src/GateMultiPhotonAnalysis.cc b/source/digits_hits/src/GateMultiPhotonAnalysis.cc index 36417e26c..0e3d25f13 100644 --- a/source/digits_hits/src/GateMultiPhotonAnalysis.cc +++ b/source/digits_hits/src/GateMultiPhotonAnalysis.cc @@ -9,6 +9,7 @@ See LICENSE.md for further details #include "GateMultiPhotonAnalysis.hh" #include "GateActions.hh" +#include "GateAnalysis.hh" #include "GateDigitizerMgr.hh" #include "GateHit.hh" #include "GateMultiPhotonAnalysisMessenger.hh" @@ -20,8 +21,10 @@ See LICENSE.md for further details #include "G4Event.hh" #include "G4HCofThisEvent.hh" +#include "G4Navigator.hh" #include "G4Run.hh" #include "G4TrajectoryContainer.hh" +#include "G4TransportationManager.hh" #include #include @@ -36,7 +39,6 @@ struct TimelineEntry { 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; }; @@ -64,9 +66,7 @@ bool EventHasProcessableHits(const std::vector &CHC_vector return false; } -std::vector BuildBasePhantomTimeline( - GatePhantomHitsCollection *PHC, - GateMultiPhotonTrajectoryNavigator *trajectoryNavigator) { +std::vector BuildBasePhantomTimeline(GatePhantomHitsCollection *PHC) { std::vector basePhantomTimeline; if (!PHC) { return basePhantomTimeline; @@ -82,18 +82,11 @@ std::vector BuildBasePhantomTimeline( 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.track_id = phantomHit->GetTrackID(); entry.phantom_hit = phantomHit; basePhantomTimeline.push_back(entry); } @@ -101,6 +94,79 @@ std::vector BuildBasePhantomTimeline( return basePhantomTimeline; } +std::string LocateVolumeName(const G4ThreeVector &position) { + G4Navigator *navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4ThreeVector direction(0., 0., 0.); + G4VPhysicalVolume *volume = navigator->LocateGlobalPointAndSetup(position, &direction, false); + if (!volume) { + return MultiPhotonAnalysisHelpers::kNoVolumeName; + } + return volume->GetName(); +} + +int CountSeptalHits(GatePhantomHitsCollection *PHC) { + // Septal penetration is configured on the GateAnalysis module (setSeptalVolumeName and + // recordSeptalPenetration), so the settings are read from there - the same way GateRootDefs + // does it when deciding whether to create the septalNb branch. + // The counting rule is copied from GateAnalysis: every phantom hit recorded in the septal + // volume increments the counter, with no per-photon attribution, and the resulting event + // total is written to all crystal hits. + GateAnalysis *analysis = + dynamic_cast(GateOutputMgr::GetInstance()->GetModule("analysis")); + if (!analysis || !analysis->GetRecordSeptalFlag() || !PHC) { + return 0; + } + + const G4String &septalVolumeName = analysis->GetSeptalPhysVolumeName(); + + int septalNb = 0; + const G4int NpHits = PHC->entries(); + for (G4int iPHit = 0; iPHit < NpHits; ++iPHit) { + GatePhantomHit *phantomHit = (*PHC)[iPHit]; + if (phantomHit && phantomHit->GetPhysVolName() == septalVolumeName) { + ++septalNb; + } + } + + return septalNb; +} + +std::unordered_map ComputePhantomTotals( + GatePhantomHitsCollection *PHC, + GateMultiPhotonTrajectoryNavigator *trajectoryNavigator) { + // Phantom counters are event totals, exactly like in GateAnalysis: they are computed in + // a pass that precedes the crystal hits, so every crystal hit of a gamma gets the same value. + // Only hits of the reference photon itself are counted (attribution by exact track ID). + std::unordered_map totals; + if (!PHC) { + return totals; + } + + const G4int NpHits = PHC->entries(); + for (G4int iPHit = 0; iPHit < NpHits; ++iPHit) { + GatePhantomHit *phantomHit = (*PHC)[iPHit]; + if (!phantomHit) { + continue; + } + + const G4int trackID = phantomHit->GetTrackID(); + if (!trajectoryNavigator->IsReferencePhoton(trackID)) { + continue; + } + + const MultiPhotonAnalysisHelpers::InteractionProcess process = + MultiPhotonAnalysisHelpers::GetInteractionProcess(phantomHit->GetProcess()); + if (process == MultiPhotonAnalysisHelpers::InteractionProcess::Other) { + continue; + } + + MultiPhotonAnalysisHelpers::AccumulatePhantomTotals( + phantomHit->GetProcess(), LocateVolumeName(phantomHit->GetPos()), totals[trackID]); + } + + return totals; +} + EventContext BuildEventContext( const G4Event *event, GateRunManager *runManager, @@ -121,8 +187,7 @@ EventContext BuildEventContext( std::vector BuildTimelineForCrystalCollection( const std::vector &basePhantomTimeline, - GateHitsCollection *CHC, - GateMultiPhotonTrajectoryNavigator *trajectoryNavigator) { + GateHitsCollection *CHC) { std::vector timeline = basePhantomTimeline; std::size_t sequence = timeline.size(); @@ -134,18 +199,13 @@ std::vector BuildTimelineForCrystalCollection( continue; } - const G4int trackID = crystalHit->GetTrackID(); - const int ancestorPhoton = trajectoryNavigator->FindAncestorPhotonTrackID(trackID); - if (ancestorPhoton == 0) { - continue; - } - + // Hits whose gamma cannot be resolved are kept as well: they receive zeroed counters, + // but eventID, runID and the remaining attributes are filled in like in GateAnalysis. TimelineEntry entry; entry.time = crystalHit->GetTime(); entry.sequence = sequence++; entry.is_phantom = false; - entry.track_id = trackID; - entry.ancestor_photon = ancestorPhoton; + entry.track_id = crystalHit->GetTrackID(); entry.crystal_hit = crystalHit; timeline.push_back(entry); } @@ -167,51 +227,75 @@ void ProcessTimeline( std::vector *timeline, GateMultiPhotonTrajectoryNavigator *trajectoryNavigator, const EventContext &context, + const std::unordered_map &phantomTotals, + int septalNb, int legacyPhotonIDPolicy) { - std::unordered_map runningStatsByPhotonTrackId; + const MultiPhotonAnalysisHelpers::PhantomStatistics kNoPhantomStatistics; + + 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; + // Counters are incremented only by hits of the reference photon itself (attribution by + // exact track ID, like in GateAnalysis) and BEFORE the hit is written, so that the value + // stored in the hit includes the current interaction. + if (trajectoryNavigator->IsReferencePhoton(entry.track_id)) { + MultiPhotonAnalysisHelpers::RunningStatistics &runningStats = + runningStatsByPhotonTrackId[entry.track_id]; + if (entry.is_phantom) { + if (entry.phantom_hit) { + MultiPhotonAnalysisHelpers::AccumulatePhantom(entry.phantom_hit->GetProcess(), runningStats); + } + } else if (entry.crystal_hit) { + MultiPhotonAnalysisHelpers::AccumulateCrystal(entry.crystal_hit->GetProcess(), runningStats); } + } - MultiPhotonAnalysisHelpers::AccumulatePhantom(entry.phantom_hit->GetProcess(), runningStats); + if (entry.is_phantom || !entry.crystal_hit || !entry.crystal_hit->GoodForAnalysis()) { continue; } - if (entry.crystal_hit) { - MultiPhotonAnalysisHelpers::AccumulateCrystal(entry.crystal_hit->GetProcess(), runningStats); + GateHit *hit = entry.crystal_hit; + + // Which gamma the hit belongs to is resolved through the ancestry chain, so hits of + // secondary tracks receive the counters of their parent gamma without incrementing them. + const int ancestorPhoton = trajectoryNavigator->FindAncestorPhotonTrackID(entry.track_id); + + const MultiPhotonAnalysisHelpers::PhantomStatistics *phantomStats = &kNoPhantomStatistics; + MultiPhotonAnalysisHelpers::RunningStatistics runningStats; + if (ancestorPhoton != 0) { + const std::unordered_map::const_iterator + phantom_it = phantomTotals.find(ancestorPhoton); + if (phantom_it != phantomTotals.end()) { + phantomStats = &(phantom_it->second); + } + + const std::unordered_map::const_iterator + running_it = runningStatsByPhotonTrackId.find(ancestorPhoton); + if (running_it != runningStatsByPhotonTrackId.end()) { + runningStats = running_it->second; + } } + + hit->SetSourceID(context.source_id); + hit->SetSourcePosition(context.source_vertex); + hit->SetNPhantomCompton(phantomStats->compton); + hit->SetNPhantomRayleigh(phantomStats->rayleigh); + // Volume names are always written, also when no scattering happened - an empty string would + // make the digitizer overwrite sourceID with -1 in the Singles and Coincidences trees. + hit->SetComptonVolumeName(phantomStats->comptonVolumeName.c_str()); + hit->SetRayleighVolumeName(phantomStats->rayleighVolumeName.c_str()); + hit->SetPhotonID(legacyPhotonIDPolicy); + hit->SetPrimaryID(trajectoryNavigator->FindPrimaryTrackID(entry.track_id)); + hit->SetEventID(context.event_id); + hit->SetRunID(context.run_id); + hit->SetNCrystalCompton(runningStats.crystalCompton); + hit->SetNCrystalRayleigh(runningStats.crystalRayleigh); + hit->SetNSeptal(septalNb); + // nInteractions is intentionally filled only in the multiphoton analysis path. + hit->SetNInteractions(runningStats.scatters); } } @@ -226,6 +310,20 @@ void RunDigitizersIfNeeded() { } } +//! Runs the digitizer chain when the enclosing scope is left, whichever path is taken. +//! GateAnalysis calls the digitizer outside of its "no trajectory container" branch, so +//! whether Singles and Coincidences are produced for an event does not depend on the +//! analysis being able to process that event. The guard gives the same guarantee here, +//! where the event is abandoned in more than one place. +class ScopedDigitizerRunner { + public: + ScopedDigitizerRunner() = default; + ~ScopedDigitizerRunner() { RunDigitizersIfNeeded(); } + + ScopedDigitizerRunner(const ScopedDigitizerRunner &) = delete; + ScopedDigitizerRunner &operator=(const ScopedDigitizerRunner &) = delete; +}; + } // namespace GateMultiPhotonAnalysis::GateMultiPhotonAnalysis(const G4String &name, GateOutputMgr *outputMgr, DigiMode digiMode) @@ -312,6 +410,18 @@ void GateMultiPhotonAnalysis::RecordEndOfEvent(const G4Event *event) { return; } + // From here on the event is digitized whatever happens to the analysis itself. + ScopedDigitizerRunner digitizerRunner; + + // An event without any primary vertex carries no track, hence no hit and no trajectory + // container. GateSourceMgr stops generating vertices once the time limit of the run is + // reached, so the last event of every run looks exactly like this. There is nothing to + // analyse and nothing anomalous about it, so it is skipped quietly instead of being + // reported as a missing trajectory container. + if (event->GetNumberOfPrimaryVertex() == 0) { + return; + } + std::vector CHC_vector = GetOutputMgr()->GetHitCollections(); GatePhantomHitsCollection *PHC = GetOutputMgr()->GetPhantomHitCollection(); @@ -359,7 +469,10 @@ void GateMultiPhotonAnalysis::RecordEndOfEvent(const G4Event *event) { m_trajectoryNavigator->SetTrajectoryContainer(trajectoryContainer); m_trajectoryNavigator->BuildIndex(); - std::vector basePhantomTimeline = BuildBasePhantomTimeline(PHC, m_trajectoryNavigator); + std::vector basePhantomTimeline = BuildBasePhantomTimeline(PHC); + const std::unordered_map phantomTotals = + ComputePhantomTotals(PHC, m_trajectoryNavigator); + const int septalNb = CountSeptalHits(PHC); const EventContext context = BuildEventContext(event, runManager, m_trajectoryNavigator); const int legacyPhotonIDPolicy = ResolveLegacyPhotonIDPolicy(); @@ -369,14 +482,10 @@ void GateMultiPhotonAnalysis::RecordEndOfEvent(const G4Event *event) { continue; } - std::vector timeline = BuildTimelineForCrystalCollection( - basePhantomTimeline, - CHC, - m_trajectoryNavigator); - ProcessTimeline(&timeline, m_trajectoryNavigator, context, legacyPhotonIDPolicy); + std::vector timeline = BuildTimelineForCrystalCollection(basePhantomTimeline, CHC); + ProcessTimeline(&timeline, m_trajectoryNavigator, context, phantomTotals, septalNb, + legacyPhotonIDPolicy); } - - RunDigitizersIfNeeded(); } void GateMultiPhotonAnalysis::RecordStepWithVolume(const GateVVolume *, const G4Step *) { @@ -396,7 +505,14 @@ bool GateMultiPhotonAnalysis::IsTrackingModeSupported(int tracking_mode_code) co return tracking_mode_code == static_cast(TrackingMode::kBoth); } -int GateMultiPhotonAnalysis::ResolveLegacyPhotonIDPolicy() const { return 0; } +int GateMultiPhotonAnalysis::ResolveLegacyPhotonIDPolicy() const { + // photonID is a GateAnalysis-specific field: it indexes the two gammas of a back-to-back + // annihilation (1 or 2, 0 when the hit comes from neither). There is no meaningful + // generalization for an arbitrary number of reference photons, and the field is not + // propagated to the Singles or Coincidences trees, so multi-photon analysis always writes 0. + // The value therefore differs from GateAnalysis in the Hits tree - this is intentional. + return 0; +} void GateMultiPhotonAnalysis::SetMissingTrajectoryPolicy(MissingTrajectoryPolicy policy) { m_missingTrajectoryPolicy = policy; diff --git a/source/digits_hits/src/GateMultiPhotonTrajectoryNavigator.cc b/source/digits_hits/src/GateMultiPhotonTrajectoryNavigator.cc index 70cdefb58..64309e2b7 100644 --- a/source/digits_hits/src/GateMultiPhotonTrajectoryNavigator.cc +++ b/source/digits_hits/src/GateMultiPhotonTrajectoryNavigator.cc @@ -16,6 +16,9 @@ See LICENSE.md for further details #include "G4ios.hh" +#include +#include + GateMultiPhotonTrajectoryNavigator::GateMultiPhotonTrajectoryNavigator() : m_tc(0), m_positronTrackID(0), @@ -111,6 +114,10 @@ std::vector GateMultiPhotonTrajectoryNavigator::FindReferencePhotonTrackIDs return out; } +bool GateMultiPhotonTrajectoryNavigator::IsReferencePhoton(int trackID) const { + return m_referencePhotonTrackIDs.find(trackID) != m_referencePhotonTrackIDs.end(); +} + int GateMultiPhotonTrajectoryNavigator::FindAncestorPhotonTrackID(int trackID) const { return MultiphotonTrajectoryResolver::ResolveAncestorPhotonTrackID( trackID, @@ -141,6 +148,12 @@ void GateMultiPhotonTrajectoryNavigator::Reset() { } void GateMultiPhotonTrajectoryNavigator::BuildReferencePhotonSet() { + // Reference photons are gammas that either start the event (primary particles - this covers + // GatePositroniumSource, back-to-back and single gamma sources, including the prompt gamma) + // or are produced by the positron (ion sources with radioactive decay, e+ sources). + // No common-vertex criterion is applied: the prompt gamma and the annihilation gammas are + // emitted from two different vertices whenever the positron range is enabled. + std::vector candidates; 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; @@ -154,27 +167,22 @@ void GateMultiPhotonTrajectoryNavigator::BuildReferencePhotonSet() { } 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); - } + const bool is_primary = (parent_id == 0); + const bool comes_from_positron = (m_positronTrackID != 0) && (parent_id == m_positronTrackID); + if (is_primary || comes_from_positron) { + candidates.push_back(track_id); } } + + // Deterministic order, so that the cut below always keeps the same photons. + std::sort(candidates.begin(), candidates.end()); + + if (candidates.size() > kMaxReferencePhotons) { + G4cout << "GateMultiPhotonTrajectoryNavigator::BuildReferencePhotonSet: WARNING found " + << candidates.size() << " reference photons, keeping the first " << kMaxReferencePhotons + << " (by track ID). Raise kMaxReferencePhotons if this process is expected." << G4endl; + candidates.resize(kMaxReferencePhotons); + } + + m_referencePhotonTrackIDs.insert(candidates.begin(), candidates.end()); } diff --git a/source/digits_hits/src/GateToRoot.cc b/source/digits_hits/src/GateToRoot.cc index 24c19de81..a3fa6a925 100644 --- a/source/digits_hits/src/GateToRoot.cc +++ b/source/digits_hits/src/GateToRoot.cc @@ -628,9 +628,17 @@ void GateToRoot::RecordEndOfAcquisition() { // m_hfile = m_treeHit->GetCurrentFile(); - // Get file only for the 1st tree as it is the same file for all trees - m_hfile = m_treesHit[0]->GetCurrentFile(); + // Get file only for the 1st tree as it is the same file for all trees. + // The vector may be empty here: RecordEndOfRun() drops the trees (it saves the + // current file beforehand) and a setup without any sensitive detector never fills + // it. In both cases m_hfile already points to the file to be written. + if (!m_treesHit.empty() && m_treesHit[0]) + m_hfile = m_treesHit[0]->GetCurrentFile(); + if (!m_hfile) { + G4cerr << "GateToRoot::RecordEndOfAcquisition(): no ROOT file to write !\n"; + return; + } if (nVerboseLevel > 0) G4cout << "GateToRoot: ROOT: files writing...\n"; @@ -686,6 +694,12 @@ void GateToRoot::RecordEndOfRun(const G4Run *) { nbPrimaries -= 1.; // Number of primaries increase too much at each end of run ! + // ROOT may have switched to another file while the run was written (files are split + // above 1.9 GBytes), so remember where the trees actually live before dropping them - + // RecordEndOfAcquisition() has no tree left to ask. + if (!m_treesHit.empty() && m_treesHit[0]) + m_hfile = m_treesHit[0]->GetCurrentFile(); + m_hitBuffers.clear(); m_treesHit.clear(); } diff --git a/source/digits_hits/src/GateVDigitizerModule.cc b/source/digits_hits/src/GateVDigitizerModule.cc index 89839b519..e1123ec68 100755 --- a/source/digits_hits/src/GateVDigitizerModule.cc +++ b/source/digits_hits/src/GateVDigitizerModule.cc @@ -219,6 +219,26 @@ void GateVDigitizerModule::InputCollectionID() } +void GateVDigitizerModule::MergeEmittedGammaInformation(GateDigi *right, GateDigi *output) +{ + // Decay level (sourceType, decayIndex) and gamma level (decayType, gammaType) fields. + // Two digis of the same event always come from the same decay, so the decay level + // fields cannot conflict; the gamma level ones can, when hits of different gammas are + // merged. Both are handled the same way as m_sourceEnergy / m_sourcePDG above: keep the + // value when it agrees, fall back to the "not known" value of the enum when it does not. + if (output->m_sourceType != right->m_sourceType) output->m_sourceType = 0; + if (output->m_decayType != right->m_decayType) output->m_decayType = 0; + if (output->m_gammaType != right->m_gammaType) output->m_gammaType = 0; + if (output->m_decayIndex != right->m_decayIndex) output->m_decayIndex = -1; + + // # of Compton and Rayleigh interactions: store the max nb, like the counters above + if ( right->m_nInteractions > output->m_nInteractions ) + { + output->m_nInteractions = right->m_nInteractions; + } +} + + GateDigi* GateVDigitizerModule::CentroidMerge(GateDigi* right, GateDigi* output ) { @@ -283,6 +303,8 @@ GateDigi* GateVDigitizerModule::CentroidMerge(GateDigi* right, GateDigi* output output->m_nSeptal = right->m_nSeptal; } + MergeEmittedGammaInformation(right, output); + // VolumeID: should be identical for both pulses, we do nothing // m_scannerPos: identical for both pulses, nothing to do // m_scannerRotAngle: identical for both pulses, nothing to do @@ -351,6 +373,8 @@ GateDigi* GateVDigitizerModule::MergePositionEnergyWin(GateDigi *right, GateDigi output->m_nSeptal = right->m_nSeptal; } + MergeEmittedGammaInformation(right, output); + // VolumeID: should be identical for both pulses, we do nothing // m_scannerPos: identical for both pulses, nothing to do // m_scannerRotAngle: identical for both pulses, nothing to do diff --git a/source/general/include/GateRootDefs.hh b/source/general/include/GateRootDefs.hh index 16d7d42ca..ad2f5bd5f 100644 --- a/source/general/include/GateRootDefs.hh +++ b/source/general/include/GateRootDefs.hh @@ -235,9 +235,9 @@ class GateRootHitBuffer Char_t RayleighVolumeName[40]; //!< Name of the last phantom-volume generating a Rayleigh Int_t volumeID[ROOT_VOLUMEIDSIZE]; //!< Volume ID Int_t septalNb; //!< HDS : septal penetration - 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 sourceType = 0; //Type of gamma source (check PositroniumSource) + Int_t decayType = 0; //Type of positronium decay (check PositroniumSource) + Int_t gammaType = 0; //Gamma type - single, annhilation, prompt (check PositroniumSource) Int_t decayIndex = -1; //Decay channel index //OK GND for CC G4bool m_CCflag; @@ -343,6 +343,14 @@ class GateRootSingleBuffer Float_t energyIni; Int_t septalNb; //!< HDS : septal penetration + + // Information from PositroniumSource, propagated from the hit through GateDigi + Int_t sourceType; //!< Type of gamma source (check PositroniumSource) + Int_t decayType; //!< Type of positronium decay (check PositroniumSource) + Int_t gammaType; //!< Gamma type - single, annihilation, prompt (check PositroniumSource) + Int_t decayIndex; //!< Decay channel index + Int_t nInteractions; //!< Number of Compton and Rayleigh interactions; -1 when not computed + // Spatial resolution standard deviations (from 2D FWHM distribution) in mm Float_t spatialRes2DStdDevX; Float_t spatialRes2DStdDevY; @@ -429,6 +437,11 @@ class GateRootCoincBuffer Int_t RayleighCrystal1; Char_t comptonVolumeName1[40]; Char_t RayleighVolumeName1[40]; + Int_t sourceType1; //!< Type of gamma source (check PositroniumSource) + Int_t decayType1; //!< Type of positronium decay (check PositroniumSource) + Int_t gammaType1; //!< Gamma type - single, annihilation, prompt (check PositroniumSource) + Int_t decayIndex1; //!< Decay channel index + Int_t nInteractions1; //!< Number of Compton and Rayleigh interactions; -1 when not computed Int_t eventID2; Int_t sourceID2; @@ -447,6 +460,11 @@ class GateRootCoincBuffer Int_t RayleighCrystal2; Char_t comptonVolumeName2[40]; Char_t RayleighVolumeName2[40]; + Int_t sourceType2; //!< Type of gamma source (check PositroniumSource) + Int_t decayType2; //!< Type of positronium decay (check PositroniumSource) + Int_t gammaType2; //!< Gamma type - single, annihilation, prompt (check PositroniumSource) + Int_t decayIndex2; //!< Decay channel index + Int_t nInteractions2; //!< Number of Compton and Rayleigh interactions; -1 when not computed // Spatial resolution standard deviations (averaged from both digis) in mm Float_t spatialRes2DStdDevX; diff --git a/source/general/src/GateRootDefs.cc b/source/general/src/GateRootDefs.cc index 736c97225..4109d5d4f 100644 --- a/source/general/src/GateRootDefs.cc +++ b/source/general/src/GateRootDefs.cc @@ -478,6 +478,13 @@ void GateRootSingleBuffer::Clear() // HDS : septal septalNb = 0; + // Information from PositroniumSource: zero and -1 mean 'not known', as in the hit buffer + sourceType = 0; + decayType = 0; + gammaType = 0; + decayIndex = -1; + nInteractions = -1; + // initialize spatial resolution stddev fields if(GetSpatialRes2DStdDevFlag()) { @@ -551,6 +558,13 @@ void GateRootSingleBuffer::Fill(GateDigi* aDigi) // HDS : septal penetration septalNb = aDigi->GetNSeptal(); + // Information from PositroniumSource plus the interaction counter + sourceType = aDigi->GetSourceType(); + decayType = aDigi->GetDecayType(); + gammaType = aDigi->GetGammaType(); + decayIndex = aDigi->GetDecayIndex(); + nInteractions = aDigi->GetNInteractions(); + if (!GetCCFlag()) { @@ -616,6 +630,19 @@ void GateSingleTree::SetBranchAddresses(TTree* singleTree, if (GateRootDefs::GetRecordSeptalFlag()) singleTree->SetBranchAddress("septalNb", &buffer.septalNb); + // Branches added later than the rest, so files written by older versions do not have + // them - reading such a file leaves the buffer at its 'not known' value + if (singleTree->GetBranch("sourceType")) + singleTree->SetBranchAddress("sourceType", &buffer.sourceType); + if (singleTree->GetBranch("decayType")) + singleTree->SetBranchAddress("decayType", &buffer.decayType); + if (singleTree->GetBranch("gammaType")) + singleTree->SetBranchAddress("gammaType", &buffer.gammaType); + if (singleTree->GetBranch("decayIndex")) + singleTree->SetBranchAddress("decayIndex", &buffer.decayIndex); + if (singleTree->GetBranch("nInteractions")) + singleTree->SetBranchAddress("nInteractions", &buffer.nInteractions); + if (!buffer.GetCCFlag()) { singleTree->SetBranchAddress("comptonPhantom", &buffer.comptonPhantom); @@ -691,6 +718,13 @@ GateDigi* GateRootSingleBuffer::CreateDigi() // Septal penetration aDigi->SetNSeptal(septalNb); + // Information from PositroniumSource plus the interaction counter + aDigi->SetSourceType(sourceType); + aDigi->SetDecayType(decayType); + aDigi->SetGammaType(gammaType); + aDigi->SetDecayIndex(decayIndex); + aDigi->SetNInteractions(nInteractions); + if (!GetCCFlag()) { aDigi->SetNPhantomCompton(comptonPhantom); @@ -770,6 +804,15 @@ void GateSingleTree::Init(GateRootSingleBuffer& buffer) // HDS : record septal penetration if (GateRootDefs::GetRecordSeptalFlag()) Branch("septalNb", &buffer.septalNb,"septalNb/I"); + // Information from PositroniumSource plus the interaction counter, propagated from the hit. + // The same branches exist in the Hits tree; nInteractions is filled by + // GateMultiPhotonAnalysis only and stays at -1 on the classical GateAnalysis path. + Branch("sourceType", &buffer.sourceType,"sourceType/I"); + Branch("decayType", &buffer.decayType,"decayType/I"); + Branch("gammaType", &buffer.gammaType,"gammaType/I"); + Branch("decayIndex", &buffer.decayIndex,"decayIndex/I"); + Branch("nInteractions", &buffer.nInteractions,"nInteractions/I"); + if (!buffer.GetCCFlag()) { if ( GateDigi::GetSingleASCIIMask(12) ) @@ -843,6 +886,11 @@ void GateRootCoincBuffer::Clear() RayleighCrystal1 = -1; strcpy (comptonVolumeName1," "); strcpy (RayleighVolumeName1," "); + sourceType1 = 0; + decayType1 = 0; + gammaType1 = 0; + decayIndex1 = -1; + nInteractions1 = -1; eventID2 = -1; sourceID2 = -1; @@ -863,7 +911,12 @@ void GateRootCoincBuffer::Clear() RayleighCrystal2 = -1; strcpy (comptonVolumeName2," "); strcpy (RayleighVolumeName2," "); - + sourceType2 = 0; + decayType2 = 0; + gammaType2 = 0; + decayIndex2 = -1; + nInteractions2 = -1; + if(GetSpatialRes2DStdDevFlag()) { // initialize spatial resolution stddev fields for coincidences @@ -903,6 +956,11 @@ void GateRootCoincBuffer::Fill(GateCoincidenceDigi* aDigi) strcpy (comptonVolumeName1,((aDigi->GetDigi(0))->GetComptonVolumeName()).c_str()); strcpy (RayleighVolumeName1,((aDigi->GetDigi(0))->GetRayleighVolumeName()).c_str()); + sourceType1 = (aDigi->GetDigi(0))->GetSourceType(); + decayType1 = (aDigi->GetDigi(0))->GetDecayType(); + gammaType1 = (aDigi->GetDigi(0))->GetGammaType(); + decayIndex1 = (aDigi->GetDigi(0))->GetDecayIndex(); + nInteractions1 = (aDigi->GetDigi(0))->GetNInteractions(); eventID2 = (aDigi->GetDigi(1))->GetEventID(); sourceID2 = (aDigi->GetDigi(1))->GetSourceID(); @@ -924,6 +982,11 @@ void GateRootCoincBuffer::Fill(GateCoincidenceDigi* aDigi) strcpy (comptonVolumeName2,((aDigi->GetDigi(1))->GetComptonVolumeName()).c_str()); strcpy (RayleighVolumeName2,((aDigi->GetDigi(1))->GetRayleighVolumeName()).c_str()); + sourceType2 = (aDigi->GetDigi(1))->GetSourceType(); + decayType2 = (aDigi->GetDigi(1))->GetDecayType(); + gammaType2 = (aDigi->GetDigi(1))->GetGammaType(); + decayIndex2 = (aDigi->GetDigi(1))->GetDecayIndex(); + nInteractions2 = (aDigi->GetDigi(1))->GetNInteractions(); if(GetSpatialRes2DStdDevFlag()) @@ -1086,6 +1149,20 @@ void GateCoincTree::Init(GateRootCoincBuffer& buffer) Branch("RayleighVolName1", (void *)buffer.RayleighVolumeName1,"RayleighVolName1/C"); if ( GateCoincidenceDigi::GetCoincidenceASCIIMask(21) ) Branch("RayleighVolName2", (void *)buffer.RayleighVolumeName2,"RayleighVolName2/C"); + + // Information from PositroniumSource plus the interaction counter, propagated from the hit + // through the digi of each arm; same branches as in the Singles tree, with the 1/2 suffix + Branch("sourceType1", &buffer.sourceType1,"sourceType1/I"); + Branch("decayType1", &buffer.decayType1,"decayType1/I"); + Branch("gammaType1", &buffer.gammaType1,"gammaType1/I"); + Branch("decayIndex1", &buffer.decayIndex1,"decayIndex1/I"); + Branch("nInteractions1", &buffer.nInteractions1,"nInteractions1/I"); + + Branch("sourceType2", &buffer.sourceType2,"sourceType2/I"); + Branch("decayType2", &buffer.decayType2,"decayType2/I"); + Branch("gammaType2", &buffer.gammaType2,"gammaType2/I"); + Branch("decayIndex2", &buffer.decayIndex2,"decayIndex2/I"); + Branch("nInteractions2", &buffer.nInteractions2,"nInteractions2/I"); } diff --git a/source/geometry/src/GateVVolume.cc b/source/geometry/src/GateVVolume.cc index 552117db8..468440219 100644 --- a/source/geometry/src/GateVVolume.cc +++ b/source/geometry/src/GateVVolume.cc @@ -385,6 +385,17 @@ void GateVVolume::AttachCrystalSD() { // GateCrystalSD *crystalSD = GateDetectorConstruction::GetGateDetectorConstruction()->GetCrystalSD(); + // Check whether this attachment is allowed or forbidden *before* building the SD: the + // GateCrystalSD constructor registers the detector in GateDigitizerMgr and nothing takes + // it back, so a refused attachment used to leave an empty hit collection behind. That + // collection renames the ROOT trees from Hits to Hits_ and changes the default + // name of the singles collection, which in turn breaks macros using a CoincidenceSorter. + if (!GateCrystalSD::CanAttachToCreator(this)) { + G4cout << "[GateVVolume::AttachCrystalSD]:\n" + << "Can not attach crystalSD!\n"; + return; + } + // OK GND 2022 //---- G4SDManager* SDman = G4SDManager::GetSDMpointer(); @@ -394,7 +405,7 @@ void GateVVolume::AttachCrystalSD() { SDman->AddNewDetector(crystalSD); //---- - // Check whether this attachement is allowed or forbidden + // Registers the system of this volume in the SD; cannot fail, the check above passed if (crystalSD->PrepareCreatorAttachment(this)) { G4cout << "[GateVVolume::AttachCrystalSD]:\n" << "Can not attach crystalSD!\n"; diff --git a/source/physics/include/GatePositroniumSourceMessenger.hh b/source/physics/include/GatePositroniumSourceMessenger.hh index 3c3e77b1b..af7d07cc1 100644 --- a/source/physics/include/GatePositroniumSourceMessenger.hh +++ b/source/physics/include/GatePositroniumSourceMessenger.hh @@ -49,7 +49,7 @@ class GatePositroniumSourceMessenger: public GateVSourceMessenger std::unique_ptr upCmdSetPositroniumFractions; std::unique_ptr upCmdSetPositroniumLifetimes; - std::unique_ptr upCmdSetPromptPhotonProbabilites; + std::unique_ptr upCmdSetPromptPhotonProbabilities; std::unique_ptr upCmdSetPromptPhotonEnergies; std::unique_ptr upCmdSetElectronCaptureProbabilities; std::unique_ptr upCmdSetDecayKinds; diff --git a/source/physics/src/GatePositroniumDecayParamsGenerator.cc b/source/physics/src/GatePositroniumDecayParamsGenerator.cc index b30471506..1c08bf5ef 100644 --- a/source/physics/src/GatePositroniumDecayParamsGenerator.cc +++ b/source/physics/src/GatePositroniumDecayParamsGenerator.cc @@ -122,7 +122,7 @@ void GatePositroniumDecayParamsGenerator::validatePositroniumDecayParams(const P GateError( "GatePositroniumDecayParamsGenerator::generatePositroniumDecayParams: " "number of provided parameters in Fractions, PromptGamma, Lifetimes, " - "Gamma Energies, Electron Capture Probabilites are not the same"); + "Gamma Energies, Electron Capture Probabilities are not the same"); } if (params.fPositronInteractions.empty()) GateWarning("PositroniumSource: setPositronInteractions was not called. " diff --git a/source/physics/src/GatePositroniumSourceMessenger.cc b/source/physics/src/GatePositroniumSourceMessenger.cc index 0e3abb8bd..f30a4e1d4 100644 --- a/source/physics/src/GatePositroniumSourceMessenger.cc +++ b/source/physics/src/GatePositroniumSourceMessenger.cc @@ -79,7 +79,7 @@ void GatePositroniumSourceMessenger::InitCommands() 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] " ) ); + upCmdSetPromptPhotonProbabilities.reset(GetStringCmd( "setPromptPhotonProbabilities", "\"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] " ) ); @@ -123,7 +123,7 @@ void GatePositroniumSourceMessenger::SetNewValue(G4UIcommand *command, G4String } else if (command == upCmdSetPositroniumLifetimes.get()) { auto lifetimes = parseListOfParamsWithUnit(new_value); fParamGenerator.SetPositroniumLifetimes(lifetimes); - } else if (command == upCmdSetPromptPhotonProbabilites.get()) { + } else if (command == upCmdSetPromptPhotonProbabilities.get()) { std::vector promptPhotonProb; std::stringstream ss(new_value); float prob;