From 088fdd5a9595adfcef2798bb374678881601920f Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Tue, 4 Mar 2025 19:12:06 +0000 Subject: [PATCH 01/26] Add actor that passes data unchanged (placeholder for implementing a pile-up actor) --- core/opengate_core/opengate_core.cpp | 3 ++ .../digitizer/GateDigitizerPileupActor.cpp | 43 +++++++++++++++ .../digitizer/GateDigitizerPileupActor.h | 48 +++++++++++++++++ .../digitizer/pyGateDigitizerPileupActor.cpp | 21 ++++++++ opengate/actors/digitizers.py | 54 +++++++++++++++++++ opengate/managers.py | 2 + 6 files changed, 171 insertions(+) create mode 100644 core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp create mode 100644 core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h create mode 100644 core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index bd0294893..3ff7b2064 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -364,6 +364,8 @@ void init_GateHitsCollectionActor(py::module &); void init_GateHitsAdderActor(py::module &); +void init_GateDigitizerPileupActor(py::module &); + void init_GateDigitizerReadoutActor(py::module &m); void init_GateDigitizerBlurringActor(py::module &m); @@ -621,6 +623,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateDigiAttributeManager(m); init_GateVDigiAttribute(m); init_GateHitsAdderActor(m); + init_GateDigitizerPileupActor(m); init_GateDigitizerReadoutActor(m); init_GateDigitizerBlurringActor(m); init_GateDigitizerEfficiencyActor(m); diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp new file mode 100644 index 000000000..5d7d31cf7 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -0,0 +1,43 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include "GateDigitizerPileupActor.h" +#include "../GateHelpersDict.h" + +GateDigitizerPileupActor::GateDigitizerPileupActor(py::dict &user_info) + : GateVDigitizerWithOutputActor(user_info, false) { + + // Actions + fActions.insert("EndOfEventAction"); +} + +GateDigitizerPileupActor::~GateDigitizerPileupActor() = default; + +void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { + GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); +} + +void GateDigitizerPileupActor::DigitInitialize( + const std::vector &attributes_not_in_filler) { + auto a = attributes_not_in_filler; + GateVDigitizerWithOutputActor::DigitInitialize(a); +} + +void GateDigitizerPileupActor::BeginOfRunAction(const G4Run *run) { + GateVDigitizerWithOutputActor::BeginOfRunAction(run); +} + +void GateDigitizerPileupActor::EndOfEventAction(const G4Event * /*unused*/) { + auto &lr = fThreadLocalVDigitizerData.Get(); + auto &iter = lr.fInputIter; + iter.GoToBegin(); + while (!iter.IsAtEnd()) { + auto &i = lr.fInputIter.fIndex; + lr.fDigiAttributeFiller->Fill(i); + iter++; + } +} diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h new file mode 100644 index 000000000..8c0989872 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -0,0 +1,48 @@ +/* -------------------------------------------------- + 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 GateDigitizerPileupActor_h +#define GateDigitizerPileupActor_h + +#include "GateVDigitizerWithOutputActor.h" +#include +#include +#include + +namespace py = pybind11; + +/* + * Digitizer module for pile-up. + */ + +class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { + +public: + // Constructor + explicit GateDigitizerPileupActor(py::dict &user_info); + + // Destructor + ~GateDigitizerPileupActor() override; + + void InitializeUserInfo(py::dict &user_info) override; + + // Called every time a Run starts (all threads) + void BeginOfRunAction(const G4Run *run) override; + + // Called every time an Event ends (all threads) + void EndOfEventAction(const G4Event *event) override; + +protected: + void DigitInitialize( + const std::vector &attributes_not_in_filler) override; + + // During computation (thread local) + struct threadLocalT {}; + G4Cache fThreadLocalData; +}; + +#endif // GateDigitizerPileupActor_h diff --git a/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp new file mode 100644 index 000000000..e7c48ca00 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp @@ -0,0 +1,21 @@ +/* -------------------------------------------------- + Copyright (C): OpenGATE Collaboration + This software is distributed under the terms + of the GNU Lesser General Public Licence (LGPL) + See LICENSE.md for further details + -------------------------------------------------- */ + +#include +#include + +namespace py = pybind11; + +#include "GateDigitizerPileupActor.h" + +void init_GateDigitizerPileupActor(py::module &m) { + + py::class_, + GateVDigitizerWithOutputActor>(m, "GateDigitizerPileupActor") + .def(py::init()); +} diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index 3882f7191..44f51daa8 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -534,6 +534,59 @@ def EndSimulationAction(self): g4.GateDigitizerBlurringActor.EndSimulationAction(self) +class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor): + """ + TODO documentation + """ + + user_info_defaults = { + "attributes": ( + [], + { + "doc": "Attributes to be considered. ", + }, + ), + "input_digi_collection": ( + "Singles", + { + "doc": "Digi collection to be used as input. ", + }, + ), + "skip_attributes": ( + [], + { + "doc": "FIXME", + }, + ), + "clear_every": ( + 1e5, + { + "doc": "FIXME", + }, + ), + } + + def __init__(self, *args, **kwargs): + DigitizerBase.__init__(self, *args, **kwargs) + self.__initcpp__() + + def __initcpp__(self): + g4.GateDigitizerPileupActor.__init__(self, self.user_info) + self.AddActions({"StartSimulationAction", "EndSimulationAction"}) + + def initialize(self): + DigitizerBase.initialize(self) + self.InitializeUserInfo(self.user_info) + self.InitializeCpp() + + def StartSimulationAction(self): + DigitizerBase.StartSimulationAction(self) + g4.GateDigitizerPileupActor.StartSimulationAction(self) + + def EndSimulationAction(self): + g4.GateDigitizerPileupActor.EndSimulationAction(self) + + class DigitizerSpatialBlurringActor( DigitizerWithRootOutput, g4.GateDigitizerSpatialBlurringActor ): @@ -1276,6 +1329,7 @@ def volume_name(self, value): process_cls(DigitizerWithRootOutput) process_cls(DigitizerAdderActor) process_cls(DigitizerBlurringActor) +process_cls(DigitizerPileupActor) process_cls(DigitizerSpatialBlurringActor) process_cls(DigitizerEfficiencyActor) process_cls(DigitizerEnergyWindowsActor) diff --git a/opengate/managers.py b/opengate/managers.py index e8f656414..c89b26042 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -112,6 +112,7 @@ DigitizerProjectionActor, DigitizerEnergyWindowsActor, DigitizerHitsCollectionActor, + DigitizerPileupActor, PhaseSpaceActor, DigiAttributeProcessDefinedStepInVolumeActor, ) @@ -153,6 +154,7 @@ "DigitizerProjectionActor": DigitizerProjectionActor, "DigitizerEnergyWindowsActor": DigitizerEnergyWindowsActor, "DigitizerHitsCollectionActor": DigitizerHitsCollectionActor, + "DigitizerPileupActor": DigitizerPileupActor, "DigiAttributeProcessDefinedStepInVolumeActor": DigiAttributeProcessDefinedStepInVolumeActor, # biasing "BremsstrahlungSplittingActor": BremsstrahlungSplittingActor, From 0b3fe455b499c49812ec9df1b8af0474ed86b20d Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Tue, 4 Mar 2025 19:35:00 +0000 Subject: [PATCH 02/26] Add test for pile-up actor --- opengate/tests/src/test097_pileup.py | 129 +++++++++++++++++++++++++++ 1 file changed, 129 insertions(+) create mode 100644 opengate/tests/src/test097_pileup.py diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py new file mode 100644 index 000000000..097ccdc74 --- /dev/null +++ b/opengate/tests/src/test097_pileup.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from pathlib import Path +from opengate.tests import utility +import uproot + +green = [0, 1, 0, 1] +gray = [0.5, 0.5, 0.5, 1] +white = [1, 1, 1, 0.8] + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "gate_test089", "test089") + + sim = gate.Simulation() + + sim.visu = False + sim.visu_type = "vrml" + sim.random_seed = 1234 + sim.number_of_threads = 1 + + # Units + mm = gate.g4_units.mm + sec = gate.g4_units.s + keV = gate.g4_units.keV + Bq = gate.g4_units.Bq + gcm3 = gate.g4_units.g_cm3 + deg = gate.g4_units.deg + + # Folders + data_path = paths.data + output_path = paths.output + + # World + world = sim.world + world.size = [450 * mm, 450 * mm, 70 * mm] + world.material = "G4_AIR" + + sim.volume_manager.material_database.add_material_weights( + "LYSO", + ["Lu", "Y", "Si", "O"], + [0.31101534, 0.368765605, 0.083209699, 0.237009356], + 5.37 * gcm3, + ) + + # Ring volume + pet = sim.add_volume("Tubs", "pet") + pet.rmax = 200 * mm + pet.rmin = 127 * mm + pet.dz = 32 * mm + pet.color = gray + pet.material = "G4_AIR" + + # Block + block = sim.add_volume("Box", "block") + block.mother = pet.name + block.size = [60 * mm, 10 * mm, 10 * mm] + translations_ring, rotations_ring = gate.geometry.utility.get_circular_repetition( + 80, [160 * mm, 0.0 * mm, 0], start_angle_deg=180, axis=[0, 0, 1] + ) + block.translation = translations_ring + block.rotation = rotations_ring + block.material = "G4_AIR" + block.color = white + + # Crystal + crystal = sim.add_volume("Box", "crystal") + crystal.mother = block.name + crystal.size = [60 * mm, 10 * mm, 10 * mm] + crystal.material = "LYSO" + crystal.color = green + + source = sim.add_source("GenericSource", "b2b") + source.particle = "back_to_back" + source.activity = 200 * 1e6 * Bq + source.position.type = "sphere" + source.position.radius = 0.5 * 1e-15 * mm + source.direction.theta = [90 * deg, 90 * deg] + source.direction.phi = [0, 360 * deg] + + # Physics + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + + stats = sim.add_actor("SimulationStatisticsActor", "Stats") + + # Hits + hc = sim.add_actor("DigitizerHitsCollectionActor", "Hits") + hc.attached_to = crystal.name + hc.authorize_repeated_volumes = True + hc.root_output.write_to_disk = False + hc.attributes = [ + "EventID", + "PostPosition", + "TotalEnergyDeposit", + "PreStepUniqueVolumeID", + "GlobalTime", + ] + + # Singles + sc = sim.add_actor("DigitizerAdderActor", "Singles_before_pileup") + sc.attached_to = hc.attached_to + sc.authorize_repeated_volumes = True + sc.input_digi_collection = hc.name + sc.policy = "EnergyWinnerPosition" + sc.output_filename = output_path / "output_singles.root" + + # Pile-up + pu = sim.add_actor("DigitizerPileupActor", "Singles_after_pileup") + pu.attached_to = hc.attached_to + pu.authorize_repeated_volumes = True + pu.input_digi_collection = sc.name + pu.output_filename = sc.output_filename + + # Timing + sim.run_timing_intervals = [[0, 0.002 * sec]] + + sim.run() + + print(stats) + + with uproot.open(sc.output_filename) as root_file: + singles_tree = root_file["Singles_before_pileup"] + pileup_tree = root_file["Singles_after_pileup"] + + print(f"{int(singles_tree.num_entries)} singles before") + print(singles_tree.arrays(library="pd")) + print(f"{int(pileup_tree.num_entries)} singles after pileup") + print(pileup_tree.arrays(library="pd")) From 09ec690f7051a363519b24f29002e92c25410997 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Thu, 6 Nov 2025 19:58:41 +0000 Subject: [PATCH 03/26] Check number of singles after pile-up and check preservation of total deposited energy --- opengate/tests/src/test097_pileup.py | 79 ++++++++++++++++++++++++---- 1 file changed, 68 insertions(+), 11 deletions(-) diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py index 097ccdc74..1d7c2c4cf 100644 --- a/opengate/tests/src/test097_pileup.py +++ b/opengate/tests/src/test097_pileup.py @@ -2,7 +2,6 @@ # -*- coding: utf-8 -*- import opengate as gate -from pathlib import Path from opengate.tests import utility import uproot @@ -10,8 +9,53 @@ gray = [0.5, 0.5, 0.5, 1] white = [1, 1, 1, 0.8] + +def num_singles_after_pileup(singles_tree, time_window): + """ + This function calculates the number of singles remaining at the output of the pile-up actor, + for a given root file containing the singles that are fed into the pile-up actor. + The time window is specified in ns. + """ + df = singles_tree.arrays(library="pd") + grouped = df.groupby("PreStepUniqueVolumeID") + + total_num_pileup_singles = 0 + total_num_pileup_groups = 0 + for volume_id, group in grouped: + num_piledup_singles = 0 + num_pileup_groups = 0 + # For each volume, start with a sorted list of singles time values in ns + times = group["GlobalTime"].sort_values().values + current = 0 # index of a single that opens the current time window + next = 1 + while next < len(times): + # Increment next until it points to the single that opens the next time window + while next < len(times) and times[next] <= times[current] + time_window: + next += 1 + if next > current + 1: + # We have found a group of at least two singles in the same time window. + num_pileup_groups += 1 + num_piledup_singles += next - current + current = next + next = current + 1 + else: + # The time window opened by current contains only one event, proceed + current += 1 + next += 1 + total_num_pileup_singles += num_piledup_singles + total_num_pileup_groups += num_pileup_groups + + print(f"Number of singles before pile-up is {len(df)}") + print( + f"{total_num_pileup_singles} singles are piled up in {total_num_pileup_groups} groups" + ) + num_singles_remaining = len(df) - total_num_pileup_singles + total_num_pileup_groups + + return num_singles_remaining + + if __name__ == "__main__": - paths = utility.get_default_test_paths(__file__, "gate_test089", "test089") + paths = utility.get_default_test_paths(__file__, "gate_test097", "test097") sim = gate.Simulation() @@ -55,9 +99,9 @@ # Block block = sim.add_volume("Box", "block") block.mother = pet.name - block.size = [60 * mm, 10 * mm, 10 * mm] + block.size = [60 * mm, 20 * mm, 20 * mm] translations_ring, rotations_ring = gate.geometry.utility.get_circular_repetition( - 80, [160 * mm, 0.0 * mm, 0], start_angle_deg=180, axis=[0, 0, 1] + 40, [160 * mm, 0.0 * mm, 0], start_angle_deg=180, axis=[0, 0, 1] ) block.translation = translations_ring block.rotation = rotations_ring @@ -67,13 +111,13 @@ # Crystal crystal = sim.add_volume("Box", "crystal") crystal.mother = block.name - crystal.size = [60 * mm, 10 * mm, 10 * mm] + crystal.size = [60 * mm, 20 * mm, 20 * mm] crystal.material = "LYSO" crystal.color = green source = sim.add_source("GenericSource", "b2b") source.particle = "back_to_back" - source.activity = 200 * 1e6 * Bq + source.activity = 10 * 1e6 * Bq source.position.type = "sphere" source.position.radius = 0.5 * 1e-15 * mm source.direction.theta = [90 * deg, 90 * deg] @@ -111,9 +155,10 @@ pu.authorize_repeated_volumes = True pu.input_digi_collection = sc.name pu.output_filename = sc.output_filename + time_window = 1000.0 # Timing - sim.run_timing_intervals = [[0, 0.002 * sec]] + sim.run_timing_intervals = [[0, 0.001 * sec]] sim.run() @@ -123,7 +168,19 @@ singles_tree = root_file["Singles_before_pileup"] pileup_tree = root_file["Singles_after_pileup"] - print(f"{int(singles_tree.num_entries)} singles before") - print(singles_tree.arrays(library="pd")) - print(f"{int(pileup_tree.num_entries)} singles after pileup") - print(pileup_tree.arrays(library="pd")) + n_expected = num_singles_after_pileup(singles_tree, time_window) + e_expected = singles_tree.arrays(library="pd")["TotalEnergyDeposit"].sum() + + print(f"Expected number of singles at output of pile-up actor is {n_expected}") + print(f"Expected total deposited energy is {e_expected:.03f} Mev") + + n_actual = int(pileup_tree.num_entries) + e_actual = pileup_tree.arrays(library="pd")["TotalEnergyDeposit"].sum() + + print(f"Actual number of singles at output of pile-up actor is {n_actual}") + print(f"Actual total deposited energy is {e_actual:.03f} MeV") + + n_ok = n_actual == n_expected + e_ok = abs(e_actual - e_expected) / e_expected < 1e-6 + + utility.test_ok(n_ok and e_ok) From 1c38ea44122266c9cf513ee1d5d5b7635cee60c1 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 9 Nov 2025 10:30:57 +0000 Subject: [PATCH 04/26] Add Python implementation of pile-up for testing + check all attributes of singles after pile-up --- opengate/tests/src/test097_pileup.py | 72 +----------- opengate/tests/src/test097_pileup_helpers.py | 117 +++++++++++++++++++ 2 files changed, 122 insertions(+), 67 deletions(-) create mode 100644 opengate/tests/src/test097_pileup_helpers.py diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py index 1d7c2c4cf..685ce5c67 100644 --- a/opengate/tests/src/test097_pileup.py +++ b/opengate/tests/src/test097_pileup.py @@ -3,57 +3,12 @@ import opengate as gate from opengate.tests import utility -import uproot +from test097_pileup_helpers import check_gate_pileup green = [0, 1, 0, 1] gray = [0.5, 0.5, 0.5, 1] white = [1, 1, 1, 0.8] - -def num_singles_after_pileup(singles_tree, time_window): - """ - This function calculates the number of singles remaining at the output of the pile-up actor, - for a given root file containing the singles that are fed into the pile-up actor. - The time window is specified in ns. - """ - df = singles_tree.arrays(library="pd") - grouped = df.groupby("PreStepUniqueVolumeID") - - total_num_pileup_singles = 0 - total_num_pileup_groups = 0 - for volume_id, group in grouped: - num_piledup_singles = 0 - num_pileup_groups = 0 - # For each volume, start with a sorted list of singles time values in ns - times = group["GlobalTime"].sort_values().values - current = 0 # index of a single that opens the current time window - next = 1 - while next < len(times): - # Increment next until it points to the single that opens the next time window - while next < len(times) and times[next] <= times[current] + time_window: - next += 1 - if next > current + 1: - # We have found a group of at least two singles in the same time window. - num_pileup_groups += 1 - num_piledup_singles += next - current - current = next - next = current + 1 - else: - # The time window opened by current contains only one event, proceed - current += 1 - next += 1 - total_num_pileup_singles += num_piledup_singles - total_num_pileup_groups += num_pileup_groups - - print(f"Number of singles before pile-up is {len(df)}") - print( - f"{total_num_pileup_singles} singles are piled up in {total_num_pileup_groups} groups" - ) - num_singles_remaining = len(df) - total_num_pileup_singles + total_num_pileup_groups - - return num_singles_remaining - - if __name__ == "__main__": paths = utility.get_default_test_paths(__file__, "gate_test097", "test097") @@ -161,26 +116,9 @@ def num_singles_after_pileup(singles_tree, time_window): sim.run_timing_intervals = [[0, 0.001 * sec]] sim.run() - print(stats) - with uproot.open(sc.output_filename) as root_file: - singles_tree = root_file["Singles_before_pileup"] - pileup_tree = root_file["Singles_after_pileup"] - - n_expected = num_singles_after_pileup(singles_tree, time_window) - e_expected = singles_tree.arrays(library="pd")["TotalEnergyDeposit"].sum() - - print(f"Expected number of singles at output of pile-up actor is {n_expected}") - print(f"Expected total deposited energy is {e_expected:.03f} Mev") - - n_actual = int(pileup_tree.num_entries) - e_actual = pileup_tree.arrays(library="pd")["TotalEnergyDeposit"].sum() - - print(f"Actual number of singles at output of pile-up actor is {n_actual}") - print(f"Actual total deposited energy is {e_actual:.03f} MeV") - - n_ok = n_actual == n_expected - e_ok = abs(e_actual - e_expected) / e_expected < 1e-6 - - utility.test_ok(n_ok and e_ok) + all_match = check_gate_pileup( + sc.output_filename, "Singles_before_pileup", "Singles_after_pileup", time_window + ) + utility.test_ok(all_match) diff --git a/opengate/tests/src/test097_pileup_helpers.py b/opengate/tests/src/test097_pileup_helpers.py new file mode 100644 index 000000000..5e39d3451 --- /dev/null +++ b/opengate/tests/src/test097_pileup_helpers.py @@ -0,0 +1,117 @@ +import pandas as pd +import uproot + + +def pileup(singles_before_pileup: pd.DataFrame, time_window: float): + """ + This function simulates pile-up with the given time window in ns. + The singles_before_pileup are in a pandas DataFrame. + It returns a dict in which the keys are volume IDs, and the values are a pandas DataFrame + representing the singles for that volume ID after pile-up. + Each single after pile-up has a TotalEnergyDeposit equal to the sum of all singles in the same time window. + The other attributes of a single after pile-up are taken from the single in the time window that has the + highest TotalEnergyDeposit. + """ + df = singles_before_pileup # df for DataFrame + grouped = df.groupby("PreStepUniqueVolumeID") + + singles_after_pileup = {} + for volume_id, group in grouped: + # For each volume, start with a sorted list of singles time values in ns. + group = group.sort_values("GlobalTime") + times = group["GlobalTime"].values + + current = 0 # index of a single that opens the current time window + next = 1 + while next < len(times): + # Increment next until it points to the single that opens the next time window. + while next < len(times) and times[next] <= times[current] + time_window: + next += 1 + if next > current + 1: + # We have found a group of at least two singles in the same time window. + # Find the single with the highest TotalEnergyDeposit in the time window. + group_slice = group.iloc[current:next] + max_energy_idx = group_slice["TotalEnergyDeposit"].idxmax() + # Create a single with the attribute values from the max energy single, + # except for the TotalEnergyDeposit, take the sum over all singles in the time window. + pileup_row = group.loc[max_energy_idx].copy() + pileup_row["TotalEnergyDeposit"] = group_slice[ + "TotalEnergyDeposit" + ].sum() + # Add the combined single to the output. + singles_after_pileup.setdefault(volume_id, []).append(pileup_row) + # Update current and next for the next time window. + current = next + next = current + 1 + else: + # The time window opened by current contains only contains one event. + # Add the original single to the output unchanged. + singles_after_pileup.setdefault(volume_id, []).append( + group.iloc[current] + ) + # Update current and next for the next time window. + current += 1 + next += 1 + # If there is only one single left, add it to the output unchanged. + if current == len(times) - 1: + singles_after_pileup.setdefault(volume_id, []).append( + group.iloc[current] + ) + + return singles_after_pileup + + +def check_gate_pileup( + root_file_path: str, + name_before_pileup: str, + name_after_pileup: str, + time_window: float, +): + """ + This function checks that the singles generated by GateDigitizerPileupActor are identical + to the ones generated by the Python pile-up implementation in the pileup() function above. + """ + with uproot.open(root_file_path) as root_file: + singles_before_pileup = root_file[name_before_pileup].arrays(library="pd") + actual_singles_after_pileup = root_file[name_after_pileup].arrays(library="pd") + + expected_singles_after_pileup = pileup(singles_before_pileup, time_window) + num_expected_singles = sum( + len(entries) for entries in expected_singles_after_pileup.values() + ) + + print(f"Singles before pile-up: {len(singles_before_pileup)}") + print( + f"Expected singles after pile-up: {num_expected_singles} ({num_expected_singles / len(singles_before_pileup) * 100:.01f}%)" + ) + print(f"Actual singles after pile-up: {len(actual_singles_after_pileup)}") + + all_match = True + for volume_id in expected_singles_after_pileup.keys(): + + # Get the expected singles (Python) and actual singles (GateDigitizerPileupActor) for the current volume. + expected_singles = pd.DataFrame(expected_singles_after_pileup[volume_id]) + actual_singles = actual_singles_after_pileup[ + actual_singles_after_pileup["PreStepUniqueVolumeID"] == volume_id + ].reset_index(drop=True) + + # Compare the number of singles + if len(expected_singles) != len(actual_singles): + print( + f"Volume {volume_id}: Expected {len(expected_singles)} singles, got {len(actual_singles)}" + ) + all_match = False + continue + + # Compare all attributes, except the volume ID + for attr in expected_singles.columns: + if attr == "PreStepUniqueVolumeID": + continue + expected_values = expected_singles[attr].values + actual_values = actual_singles[attr].values + if not all(expected_values == actual_values): + print(f"Volume {volume_id}: Attribute {attr} does not match") + all_match = False + break + + return all_match From c668fddbb94054815908683637a9492b419b24c6 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 9 Nov 2025 10:31:24 +0000 Subject: [PATCH 05/26] Implement pile-up --- .../digitizer/GateDigitizerPileupActor.cpp | 190 +++++++++++++++++- .../digitizer/GateDigitizerPileupActor.h | 43 +++- .../digitizer/pyGateDigitizerPileupActor.cpp | 4 +- opengate/actors/digitizers.py | 14 ++ 4 files changed, 244 insertions(+), 7 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 5d7d31cf7..7e1be24c7 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -10,34 +10,214 @@ GateDigitizerPileupActor::GateDigitizerPileupActor(py::dict &user_info) : GateVDigitizerWithOutputActor(user_info, false) { - - // Actions fActions.insert("EndOfEventAction"); + fActions.insert("EndOfRunAction"); } GateDigitizerPileupActor::~GateDigitizerPileupActor() = default; void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); + + // Get time window parameter in ns. + fTimeWindow = 1000.0; // default value + if (py::len(user_info) > 0 && user_info.contains("time_window")) { + fTimeWindow = DictGetDouble(user_info, "time_window"); + } + fGroupVolumeDepth = -1; +} + +void GateDigitizerPileupActor::SetGroupVolumeDepth(const int depth) { + fGroupVolumeDepth = depth; } void GateDigitizerPileupActor::DigitInitialize( const std::vector &attributes_not_in_filler) { + // We handle all attributes explicitly by storing their values, + // because we need to keep singles across consecutive events. auto a = attributes_not_in_filler; + for (const auto &att_name : fInputDigiCollection->GetDigiAttributeNames()) { + a.push_back(att_name); + } + GateVDigitizerWithOutputActor::DigitInitialize(a); + + // Get output attribute pointers + fOutputEdepAttribute = + fOutputDigiCollection->GetDigiAttribute("TotalEnergyDeposit"); + fOutputGlobalTimeAttribute = + fOutputDigiCollection->GetDigiAttribute("GlobalTime"); + fOutputVolumeIDAttribute = + fOutputDigiCollection->GetDigiAttribute("PreStepUniqueVolumeID"); + + // Set up pointers to track specific attributes + auto &lr = fThreadLocalVDigitizerData.Get(); + auto &l = fThreadLocalData.Get(); + + lr.fInputIter = fInputDigiCollection->NewIterator(); + lr.fInputIter.TrackAttribute("TotalEnergyDeposit", &l.edep); + lr.fInputIter.TrackAttribute("GlobalTime", &l.time); + lr.fInputIter.TrackAttribute("PreStepUniqueVolumeID", &l.volID); } void GateDigitizerPileupActor::BeginOfRunAction(const G4Run *run) { GateVDigitizerWithOutputActor::BeginOfRunAction(run); } -void GateDigitizerPileupActor::EndOfEventAction(const G4Event * /*unused*/) { +void GateDigitizerPileupActor::StoreAttributeValues( + threadLocalT::PileupGroup &group, size_t index) { + // Clear previous values + group.stored_attributes.clear(); + + // Store values from all attributes in the input collection + for (auto *att : fInputDigiCollection->GetDigiAttributes()) { + auto name = att->GetDigiAttributeName(); + auto type = att->GetDigiAttributeType(); + + // Skip the attributes we handle explicitly + if (name == "TotalEnergyDeposit" || name == "GlobalTime" || + name == "PreStepUniqueVolumeID") { + continue; + } + + // Store value based on type + switch (type) { + case 'D': // double + group.stored_attributes[name] = att->GetDValues()[index]; + break; + case 'I': // int + group.stored_attributes[name] = att->GetIValues()[index]; + break; + case 'L': // int64_t + group.stored_attributes[name] = att->GetLValues()[index]; + break; + case 'S': // string + group.stored_attributes[name] = att->GetSValues()[index]; + break; + case '3': // G4ThreeVector + group.stored_attributes[name] = att->Get3Values()[index]; + break; + case 'U': // GateUniqueVolumeID::Pointer + group.stored_attributes[name] = att->GetUValues()[index]; + break; + } + } +} + +void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { + ProcessPileup(); + + // Create output singles from piled-up groups of input singles. + + auto &l = fThreadLocalData.Get(); + + for (auto &[volume_id, groups] : l.volume_groups) { + if (!groups.empty()) { + // Handle all groups, except the last one, because the last group may + // still get contributions from upcoming events. + for (auto it = groups.begin(); it != std::prev(groups.end()); ++it) { + FillAttributeValues(*it); + } + groups.erase(groups.begin(), std::prev(groups.end())); + } + } +} + +void GateDigitizerPileupActor::ProcessPileup() { + auto &lr = fThreadLocalVDigitizerData.Get(); + auto &l = fThreadLocalData.Get(); auto &iter = lr.fInputIter; + iter.GoToBegin(); while (!iter.IsAtEnd()) { - auto &i = lr.fInputIter.fIndex; - lr.fDigiAttributeFiller->Fill(i); + + const auto current_time = *l.time; + const auto current_edep = *l.edep; + const auto current_vol = + l.volID->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth); + const auto current_index = iter.fIndex; + + bool added_to_existing_group = false; + auto &groups = l.volume_groups[current_vol]; + if (groups.size() > 0) { + auto &group = groups.back(); + if (std::abs(current_time - group.first_time) <= fTimeWindow) { + // Accumulate deposited energy of all singles in the same time window. + group.total_edep += current_edep; + // Keep the attributes of the highest energy single. + if (current_edep > group.highest_edep) { + group.highest_edep = current_edep; + group.time = current_time; + // Store all other attribute values from this single. + StoreAttributeValues(group, current_index); + } + added_to_existing_group = true; + } + } + + // If not added to an existing group, create a new group. + if (!added_to_existing_group) { + typename threadLocalT::PileupGroup new_group; + new_group.total_edep = current_edep; + new_group.highest_edep = current_edep; + new_group.first_time = current_time; + new_group.time = current_time; + new_group.volume_id = *l.volID; + // Store all other attribute values from this single + StoreAttributeValues(new_group, current_index); + groups.push_back(new_group); + } + iter++; } } + +void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) { + + auto &l = fThreadLocalData.Get(); + + // Output any unfinished groups that are still present. + for (auto &[volume_id, groups] : l.volume_groups) { + if (!groups.empty()) { + for (auto &group : groups) { + FillAttributeValues(group); + } + groups.erase(groups.begin(), groups.end()); + } + } + + // Make sure everything is output into the root file. + fOutputDigiCollection->FillToRootIfNeeded(true); +} + +void GateDigitizerPileupActor::FillAttributeValues( + const threadLocalT::PileupGroup &group) { + + fOutputEdepAttribute->FillDValue(group.total_edep); + fOutputGlobalTimeAttribute->FillDValue(group.time); + fOutputVolumeIDAttribute->FillUValue(group.volume_id); + + // Fill all other stored attributes + for (const auto &[name, value] : group.stored_attributes) { + auto *att = fOutputDigiCollection->GetDigiAttribute(name); + std::visit( + [att](auto &&arg) { + using T = std::decay_t; + if constexpr (std::is_same_v) { + att->FillDValue(arg); + } else if constexpr (std::is_same_v) { + att->FillIValue(arg); + } else if constexpr (std::is_same_v) { + att->FillLValue(arg); + } else if constexpr (std::is_same_v) { + att->FillSValue(arg); + } else if constexpr (std::is_same_v) { + att->Fill3Value(arg); + } else if constexpr (std::is_same_v) { + att->FillUValue(arg); + } + }, + value); + } +} \ No newline at end of file diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index 8c0989872..63e530708 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -11,7 +11,9 @@ #include "GateVDigitizerWithOutputActor.h" #include #include +#include #include +#include namespace py = pybind11; @@ -36,13 +38,52 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { // Called every time an Event ends (all threads) void EndOfEventAction(const G4Event *event) override; + // Called every time a Run ends (all threads) + void EndOfRunAction(const G4Run *run) override; + + void SetGroupVolumeDepth(int depth); + protected: void DigitInitialize( const std::vector &attributes_not_in_filler) override; + // User parameters + int fGroupVolumeDepth; + double fTimeWindow; + + // Output attribute pointers + GateVDigiAttribute *fOutputEdepAttribute{}; + GateVDigiAttribute *fOutputGlobalTimeAttribute{}; + GateVDigiAttribute *fOutputVolumeIDAttribute{}; + + void ProcessPileup(); + // During computation (thread local) - struct threadLocalT {}; + struct threadLocalT { + double *edep; + double *time; + GateUniqueVolumeID::Pointer *volID; + + // Storage for piled up singles + struct PileupGroup { + double highest_edep; + double total_edep; + double first_time; + double time; + GateUniqueVolumeID::Pointer volume_id; + + using AttributeValue = + std::variant; + std::map stored_attributes; + }; + + std::map> volume_groups; + }; G4Cache fThreadLocalData; + + void StoreAttributeValues(threadLocalT::PileupGroup &group, size_t index); + void FillAttributeValues(const threadLocalT::PileupGroup &group); }; #endif // GateDigitizerPileupActor_h diff --git a/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp index e7c48ca00..6a8b7d3db 100644 --- a/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.cpp @@ -17,5 +17,7 @@ void init_GateDigitizerPileupActor(py::module &m) { py::class_, GateVDigitizerWithOutputActor>(m, "GateDigitizerPileupActor") - .def(py::init()); + .def(py::init()) + .def("SetGroupVolumeDepth", + &GateDigitizerPileupActor::SetGroupVolumeDepth); } diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index 44f51daa8..819d706d9 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -564,6 +564,12 @@ class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor) "doc": "FIXME", }, ), + "group_volume": ( + None, + { + "doc": "FIXME", + }, + ), } def __init__(self, *args, **kwargs): @@ -579,6 +585,14 @@ def initialize(self): self.InitializeUserInfo(self.user_info) self.InitializeCpp() + def set_group_by_depth(self): + depth = -1 + if self.user_info.group_volume is not None: + depth = self.simulation.volume_manager.get_volume( + self.user_info.group_volume + ).volume_depth_in_tree + self.SetGroupVolumeDepth(depth) + def StartSimulationAction(self): DigitizerBase.StartSimulationAction(self) g4.GateDigitizerPileupActor.StartSimulationAction(self) From 64edabcb066e4b55c2faac64f0a341da46413ff4 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 9 Nov 2025 11:45:42 +0000 Subject: [PATCH 06/26] Fix missing time_window parameter --- .../opengate_lib/digitizer/GateDigitizerPileupActor.cpp | 2 +- opengate/actors/digitizers.py | 6 ++++++ opengate/tests/src/test097_pileup.py | 7 +++++-- 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 7e1be24c7..8d76afbd0 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -20,7 +20,7 @@ void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); // Get time window parameter in ns. - fTimeWindow = 1000.0; // default value + fTimeWindow = 0.0; // default value, no pile-up if (py::len(user_info) > 0 && user_info.contains("time_window")) { fTimeWindow = DictGetDouble(user_info, "time_window"); } diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index 819d706d9..ce2227990 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -564,6 +564,12 @@ class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor) "doc": "FIXME", }, ), + "time_window": ( + None, + { + "doc": "FIXME", + }, + ), "group_volume": ( None, { diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py index 685ce5c67..b325e6c92 100644 --- a/opengate/tests/src/test097_pileup.py +++ b/opengate/tests/src/test097_pileup.py @@ -110,7 +110,7 @@ pu.authorize_repeated_volumes = True pu.input_digi_collection = sc.name pu.output_filename = sc.output_filename - time_window = 1000.0 + pu.time_window = 1000.0 # Timing sim.run_timing_intervals = [[0, 0.001 * sec]] @@ -119,6 +119,9 @@ print(stats) all_match = check_gate_pileup( - sc.output_filename, "Singles_before_pileup", "Singles_after_pileup", time_window + sc.output_filename, + "Singles_before_pileup", + "Singles_after_pileup", + pu.time_window, ) utility.test_ok(all_match) From de92812584f5a730e155e5e644dad46f84b4093b Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Wed, 3 Dec 2025 08:21:11 +0100 Subject: [PATCH 07/26] make FillAttributeValues const --- .../opengate_lib/digitizer/GateDigitizerPileupActor.cpp | 8 +++++--- .../opengate_lib/digitizer/GateDigitizerPileupActor.h | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 8d76afbd0..3410a68dd 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -33,7 +33,7 @@ void GateDigitizerPileupActor::SetGroupVolumeDepth(const int depth) { void GateDigitizerPileupActor::DigitInitialize( const std::vector &attributes_not_in_filler) { - // We handle all attributes explicitly by storing their values, + // We handle all attributes explicitly by storing their values // because we need to keep singles across consecutive events. auto a = attributes_not_in_filler; for (const auto &att_name : fInputDigiCollection->GetDigiAttributeNames()) { @@ -100,6 +100,9 @@ void GateDigitizerPileupActor::StoreAttributeValues( case 'U': // GateUniqueVolumeID::Pointer group.stored_attributes[name] = att->GetUValues()[index]; break; + default: + DDD(type); + Fatal("Error GateDigitizerPileupActor while storing attribute value"); } } } @@ -174,7 +177,6 @@ void GateDigitizerPileupActor::ProcessPileup() { } void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) { - auto &l = fThreadLocalData.Get(); // Output any unfinished groups that are still present. @@ -192,7 +194,7 @@ void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) { } void GateDigitizerPileupActor::FillAttributeValues( - const threadLocalT::PileupGroup &group) { + const threadLocalT::PileupGroup &group) const { fOutputEdepAttribute->FillDValue(group.total_edep); fOutputGlobalTimeAttribute->FillDValue(group.time); diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index 63e530708..48054e1b9 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -83,7 +83,7 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { G4Cache fThreadLocalData; void StoreAttributeValues(threadLocalT::PileupGroup &group, size_t index); - void FillAttributeValues(const threadLocalT::PileupGroup &group); + void FillAttributeValues(const threadLocalT::PileupGroup &group) const; }; #endif // GateDigitizerPileupActor_h From 51666af88e82422ad6b39a227ccd5c8877408e0e Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Wed, 3 Dec 2025 08:21:26 +0100 Subject: [PATCH 08/26] make test executable --- opengate/tests/src/test097_pileup.py | 0 opengate/tests/src/test097_pileup_helpers.py | 21 +++++++++++--------- 2 files changed, 12 insertions(+), 9 deletions(-) mode change 100644 => 100755 opengate/tests/src/test097_pileup.py diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py old mode 100644 new mode 100755 diff --git a/opengate/tests/src/test097_pileup_helpers.py b/opengate/tests/src/test097_pileup_helpers.py index 5e39d3451..ab5933cec 100644 --- a/opengate/tests/src/test097_pileup_helpers.py +++ b/opengate/tests/src/test097_pileup_helpers.py @@ -22,15 +22,18 @@ def pileup(singles_before_pileup: pd.DataFrame, time_window: float): times = group["GlobalTime"].values current = 0 # index of a single that opens the current time window - next = 1 - while next < len(times): + next_single = 1 + while next_single < len(times): # Increment next until it points to the single that opens the next time window. - while next < len(times) and times[next] <= times[current] + time_window: - next += 1 - if next > current + 1: + while ( + next_single < len(times) + and times[next_single] <= times[current] + time_window + ): + next_single += 1 + if next_single > current + 1: # We have found a group of at least two singles in the same time window. # Find the single with the highest TotalEnergyDeposit in the time window. - group_slice = group.iloc[current:next] + group_slice = group.iloc[current:next_single] max_energy_idx = group_slice["TotalEnergyDeposit"].idxmax() # Create a single with the attribute values from the max energy single, # except for the TotalEnergyDeposit, take the sum over all singles in the time window. @@ -41,8 +44,8 @@ def pileup(singles_before_pileup: pd.DataFrame, time_window: float): # Add the combined single to the output. singles_after_pileup.setdefault(volume_id, []).append(pileup_row) # Update current and next for the next time window. - current = next - next = current + 1 + current = next_single + next_single = current + 1 else: # The time window opened by current contains only contains one event. # Add the original single to the output unchanged. @@ -51,7 +54,7 @@ def pileup(singles_before_pileup: pd.DataFrame, time_window: float): ) # Update current and next for the next time window. current += 1 - next += 1 + next_single += 1 # If there is only one single left, add it to the output unchanged. if current == len(times) - 1: singles_after_pileup.setdefault(volume_id, []).append( From a2e50990586818c600ae188579db4703912e473a Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 21 Dec 2025 05:32:35 +0000 Subject: [PATCH 09/26] use GateDigiCollection as temporary storage for digis in PileupActor --- .../digitizer/GateDigitizerPileupActor.cpp | 256 +++++++----------- .../digitizer/GateDigitizerPileupActor.h | 55 ++-- 2 files changed, 126 insertions(+), 185 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 3410a68dd..730699cc9 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -7,6 +7,9 @@ #include "GateDigitizerPileupActor.h" #include "../GateHelpersDict.h" +#include "GateDigiCollectionManager.h" +#include "GateHelpersDigitizer.h" +#include GateDigitizerPileupActor::GateDigitizerPileupActor(py::dict &user_info) : GateVDigitizerWithOutputActor(user_info, false) { @@ -17,14 +20,15 @@ GateDigitizerPileupActor::GateDigitizerPileupActor(py::dict &user_info) GateDigitizerPileupActor::~GateDigitizerPileupActor() = default; void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { - GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); + GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); // Get time window parameter in ns. fTimeWindow = 0.0; // default value, no pile-up if (py::len(user_info) > 0 && user_info.contains("time_window")) { fTimeWindow = DictGetDouble(user_info, "time_window"); } fGroupVolumeDepth = -1; + fInputDigiCollectionName = DictGetStr(user_info, "input_digi_collection"); } void GateDigitizerPileupActor::SetGroupVolumeDepth(const int depth) { @@ -33,193 +37,137 @@ void GateDigitizerPileupActor::SetGroupVolumeDepth(const int depth) { void GateDigitizerPileupActor::DigitInitialize( const std::vector &attributes_not_in_filler) { - // We handle all attributes explicitly by storing their values - // because we need to keep singles across consecutive events. - auto a = attributes_not_in_filler; - for (const auto &att_name : fInputDigiCollection->GetDigiAttributeNames()) { - a.push_back(att_name); - } + auto a = attributes_not_in_filler; + a.push_back("TotalEnergyDeposit"); GateVDigitizerWithOutputActor::DigitInitialize(a); - // Get output attribute pointers + // Get output attribute pointer fOutputEdepAttribute = fOutputDigiCollection->GetDigiAttribute("TotalEnergyDeposit"); - fOutputGlobalTimeAttribute = - fOutputDigiCollection->GetDigiAttribute("GlobalTime"); - fOutputVolumeIDAttribute = - fOutputDigiCollection->GetDigiAttribute("PreStepUniqueVolumeID"); // Set up pointers to track specific attributes auto &lr = fThreadLocalVDigitizerData.Get(); auto &l = fThreadLocalData.Get(); - lr.fInputIter = fInputDigiCollection->NewIterator(); - lr.fInputIter.TrackAttribute("TotalEnergyDeposit", &l.edep); lr.fInputIter.TrackAttribute("GlobalTime", &l.time); lr.fInputIter.TrackAttribute("PreStepUniqueVolumeID", &l.volID); } -void GateDigitizerPileupActor::BeginOfRunAction(const G4Run *run) { - GateVDigitizerWithOutputActor::BeginOfRunAction(run); -} - -void GateDigitizerPileupActor::StoreAttributeValues( - threadLocalT::PileupGroup &group, size_t index) { - // Clear previous values - group.stored_attributes.clear(); +void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { + auto &lr = fThreadLocalVDigitizerData.Get(); + auto &l = fThreadLocalData.Get(); + auto &inputIter = lr.fInputIter; - // Store values from all attributes in the input collection - for (auto *att : fInputDigiCollection->GetDigiAttributes()) { - auto name = att->GetDigiAttributeName(); - auto type = att->GetDigiAttributeType(); + inputIter.GoToBegin(); + while (!inputIter.IsAtEnd()) { + // Look up or create the pile-up window object for the volume to which the + // current digi belongs. + auto &window = GetPileupWindowForVolume(l.volID); - // Skip the attributes we handle explicitly - if (name == "TotalEnergyDeposit" || name == "GlobalTime" || - name == "PreStepUniqueVolumeID") { - continue; + const auto current_time = *l.time; + if (window.digis->GetSize() == 0) { + // The window has no digis yet: make the window start at the time of the + // current digi. + window.startTime = current_time; + } else if (current_time - window.startTime > fTimeWindow) { + // The current digi is beyond the time window: process the digis that are + // currently in the window, then make the window start at the time of the + // current digi. + ProcessPileupWindow(window); + window.startTime = current_time; } - // Store value based on type - switch (type) { - case 'D': // double - group.stored_attributes[name] = att->GetDValues()[index]; - break; - case 'I': // int - group.stored_attributes[name] = att->GetIValues()[index]; - break; - case 'L': // int64_t - group.stored_attributes[name] = att->GetLValues()[index]; - break; - case 'S': // string - group.stored_attributes[name] = att->GetSValues()[index]; - break; - case '3': // G4ThreeVector - group.stored_attributes[name] = att->Get3Values()[index]; - break; - case 'U': // GateUniqueVolumeID::Pointer - group.stored_attributes[name] = att->GetUValues()[index]; - break; - default: - DDD(type); - Fatal("Error GateDigitizerPileupActor while storing attribute value"); - } + // Add the current digi to the window. + window.fillerIn->Fill(inputIter.fIndex); + + inputIter++; } } -void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { - ProcessPileup(); - - // Create output singles from piled-up groups of input singles. - +void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) { auto &l = fThreadLocalData.Get(); - - for (auto &[volume_id, groups] : l.volume_groups) { - if (!groups.empty()) { - // Handle all groups, except the last one, because the last group may - // still get contributions from upcoming events. - for (auto it = groups.begin(); it != std::prev(groups.end()); ++it) { - FillAttributeValues(*it); - } - groups.erase(groups.begin(), std::prev(groups.end())); - } + for (auto &[_vol_hash, window] : l.fVolumePileupWindows) { + ProcessPileupWindow(window); } + // Make sure everything is output into the root file. + fOutputDigiCollection->FillToRootIfNeeded(true); } -void GateDigitizerPileupActor::ProcessPileup() { +GateDigitizerPileupActor::PileupWindow & +GateDigitizerPileupActor::GetPileupWindowForVolume( + GateUniqueVolumeID::Pointer *volume) { + // This function looks up the PileupWindow object for the given volume. If it + // does not yet exist for the volume, it creates a PileupWindow. - auto &lr = fThreadLocalVDigitizerData.Get(); + const auto vol_hash = volume->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth); auto &l = fThreadLocalData.Get(); - auto &iter = lr.fInputIter; - - iter.GoToBegin(); - while (!iter.IsAtEnd()) { - - const auto current_time = *l.time; - const auto current_edep = *l.edep; - const auto current_vol = - l.volID->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth); - const auto current_index = iter.fIndex; - - bool added_to_existing_group = false; - auto &groups = l.volume_groups[current_vol]; - if (groups.size() > 0) { - auto &group = groups.back(); - if (std::abs(current_time - group.first_time) <= fTimeWindow) { - // Accumulate deposited energy of all singles in the same time window. - group.total_edep += current_edep; - // Keep the attributes of the highest energy single. - if (current_edep > group.highest_edep) { - group.highest_edep = current_edep; - group.time = current_time; - // Store all other attribute values from this single. - StoreAttributeValues(group, current_index); - } - added_to_existing_group = true; - } - } - // If not added to an existing group, create a new group. - if (!added_to_existing_group) { - typename threadLocalT::PileupGroup new_group; - new_group.total_edep = current_edep; - new_group.highest_edep = current_edep; - new_group.first_time = current_time; - new_group.time = current_time; - new_group.volume_id = *l.volID; - // Store all other attribute values from this single - StoreAttributeValues(new_group, current_index); - groups.push_back(new_group); - } - - iter++; + // Look up the window based + auto it = l.fVolumePileupWindows.find(vol_hash); + if (it != l.fVolumePileupWindows.end()) { + // Return a reference to the existing PileupWindow object for the volume. + return it->second; + } else { + // A PileupWindow object does not yet exist for this volume: create one. + PileupWindow window; + const auto vol_id = volume->get()->GetIdUpToDepth(fGroupVolumeDepth); + // Create a GateDigiCollection for this volume, as a temporary storage for + // digis that belong to the same time window (the name must be unique). + window.digis = GateDigiCollectionManager::GetInstance()->NewDigiCollection( + GetName() + "_" + vol_id); + window.digis->InitDigiAttributesFromCopy(fInputDigiCollection); + // Create an iterator to be used when digis will be combined into one digi, + // due to pile-up. + window.digiIter = window.digis->NewIterator(); + window.digiIter.TrackAttribute("TotalEnergyDeposit", &l.edep); + // Create a filler to copy all digi attributes from the input collection + // into the collection of the window. + window.fillerIn = std::make_unique( + fInputDigiCollection, window.digis, + window.digis->GetDigiAttributeNames()); + // Create a filler to copy digi attributes from the collection of the window + // to the output collection (used for the digis that will result from + // pile-up). + auto filler_out_attributes = window.digis->GetDigiAttributeNames(); + filler_out_attributes.erase("TotalEnergyDeposit"); + window.fillerOut = std::make_unique( + window.digis, fOutputDigiCollection, filler_out_attributes); + + // Store the PileupWindow in the map and return a reference. + l.fVolumePileupWindows[vol_hash] = std::move(window); + return l.fVolumePileupWindows[vol_hash]; } } -void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) { +void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { + // This function simulates pile-up by combining the digis in the given window + // into one digi. auto &l = fThreadLocalData.Get(); - // Output any unfinished groups that are still present. - for (auto &[volume_id, groups] : l.volume_groups) { - if (!groups.empty()) { - for (auto &group : groups) { - FillAttributeValues(group); - } - groups.erase(groups.begin(), groups.end()); + std::optional highest_edep{}; + double total_edep = 0.0; + size_t highest_edep_index = 0; + + // Iterate over all digis in the window from the beginning. + window.digiIter.Reset(); + while (!window.digiIter.IsAtEnd()) { + const auto current_edep = *l.edep; + // Remember the value and index of the highest deposited energy so far. + if (!highest_edep || current_edep > highest_edep) { + highest_edep = current_edep; + highest_edep_index = window.digiIter.fIndex; } + // Accumulate all deposited energy values. + total_edep += current_edep; + window.digiIter++; } - // Make sure everything is output into the root file. - fOutputDigiCollection->FillToRootIfNeeded(true); + // The resulting pile-up digi gets the total edep value as its edep value. + fOutputEdepAttribute->FillDValue(total_edep); + // All the other attribute values are taken from the digi with the highest + // edep. + window.fillerOut->Fill(highest_edep_index); + // Remove all processed digis from the window. + window.digis->Clear(); } - -void GateDigitizerPileupActor::FillAttributeValues( - const threadLocalT::PileupGroup &group) const { - - fOutputEdepAttribute->FillDValue(group.total_edep); - fOutputGlobalTimeAttribute->FillDValue(group.time); - fOutputVolumeIDAttribute->FillUValue(group.volume_id); - - // Fill all other stored attributes - for (const auto &[name, value] : group.stored_attributes) { - auto *att = fOutputDigiCollection->GetDigiAttribute(name); - std::visit( - [att](auto &&arg) { - using T = std::decay_t; - if constexpr (std::is_same_v) { - att->FillDValue(arg); - } else if constexpr (std::is_same_v) { - att->FillIValue(arg); - } else if constexpr (std::is_same_v) { - att->FillLValue(arg); - } else if constexpr (std::is_same_v) { - att->FillSValue(arg); - } else if constexpr (std::is_same_v) { - att->Fill3Value(arg); - } else if constexpr (std::is_same_v) { - att->FillUValue(arg); - } - }, - value); - } -} \ No newline at end of file diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index 48054e1b9..608fb4da8 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -12,8 +12,8 @@ #include #include #include +#include #include -#include namespace py = pybind11; @@ -24,17 +24,12 @@ namespace py = pybind11; class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { public: - // Constructor explicit GateDigitizerPileupActor(py::dict &user_info); - // Destructor ~GateDigitizerPileupActor() override; void InitializeUserInfo(py::dict &user_info) override; - // Called every time a Run starts (all threads) - void BeginOfRunAction(const G4Run *run) override; - // Called every time an Event ends (all threads) void EndOfEventAction(const G4Event *event) override; @@ -51,39 +46,37 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { int fGroupVolumeDepth; double fTimeWindow; - // Output attribute pointers + // Output attribute pointer GateVDigiAttribute *fOutputEdepAttribute{}; - GateVDigiAttribute *fOutputGlobalTimeAttribute{}; - GateVDigiAttribute *fOutputVolumeIDAttribute{}; - void ProcessPileup(); + // Struct for storing digis in one particular volume which belong to the same + // time window. + struct PileupWindow { + // Time at which the time window opens. + double startTime{}; + // Collection of digis in the same time window. + GateDigiCollection *digis; + // Iterator used to loop over the digis for simulating pile-up. + GateDigiCollectionIterator digiIter; + // Filler used to copy digi attributes from the input collection into the + // window. + std::unique_ptr fillerIn; + // Filler used to copy the pile-up digi from the window into the output + // collection. + std::unique_ptr fillerOut; + }; + + PileupWindow &GetPileupWindowForVolume(GateUniqueVolumeID::Pointer *volume); + void ProcessPileupWindow(PileupWindow &window); - // During computation (thread local) struct threadLocalT { - double *edep; - double *time; GateUniqueVolumeID::Pointer *volID; + double *time; + double *edep; - // Storage for piled up singles - struct PileupGroup { - double highest_edep; - double total_edep; - double first_time; - double time; - GateUniqueVolumeID::Pointer volume_id; - - using AttributeValue = - std::variant; - std::map stored_attributes; - }; - - std::map> volume_groups; + std::map fVolumePileupWindows; }; G4Cache fThreadLocalData; - - void StoreAttributeValues(threadLocalT::PileupGroup &group, size_t index); - void FillAttributeValues(const threadLocalT::PileupGroup &group) const; }; #endif // GateDigitizerPileupActor_h From e75b8341c99b7c41974f9661557b6cd082cdfd5c Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 21 Dec 2025 18:39:14 +0000 Subject: [PATCH 10/26] reduce number of calls to fThreadLocalData.Get() for performance --- .../digitizer/GateDigitizerPileupActor.cpp | 20 ++++++++++--------- .../digitizer/GateDigitizerPileupActor.h | 5 ++++- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 730699cc9..0e39506bc 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -63,7 +63,8 @@ void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { while (!inputIter.IsAtEnd()) { // Look up or create the pile-up window object for the volume to which the // current digi belongs. - auto &window = GetPileupWindowForVolume(l.volID); + auto &window = + GetPileupWindowForCurrentVolume(l.volID, l.fVolumePileupWindows); const auto current_time = *l.time; if (window.digis->GetSize() == 0) { @@ -95,17 +96,17 @@ void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) { } GateDigitizerPileupActor::PileupWindow & -GateDigitizerPileupActor::GetPileupWindowForVolume( - GateUniqueVolumeID::Pointer *volume) { +GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( + GateUniqueVolumeID::Pointer *volume, + std::map &windows) { // This function looks up the PileupWindow object for the given volume. If it // does not yet exist for the volume, it creates a PileupWindow. const auto vol_hash = volume->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth); - auto &l = fThreadLocalData.Get(); // Look up the window based - auto it = l.fVolumePileupWindows.find(vol_hash); - if (it != l.fVolumePileupWindows.end()) { + auto it = windows.find(vol_hash); + if (it != windows.end()) { // Return a reference to the existing PileupWindow object for the volume. return it->second; } else { @@ -120,7 +121,8 @@ GateDigitizerPileupActor::GetPileupWindowForVolume( // Create an iterator to be used when digis will be combined into one digi, // due to pile-up. window.digiIter = window.digis->NewIterator(); - window.digiIter.TrackAttribute("TotalEnergyDeposit", &l.edep); + window.digiIter.TrackAttribute("TotalEnergyDeposit", + &fThreadLocalData.Get().edep); // Create a filler to copy all digi attributes from the input collection // into the collection of the window. window.fillerIn = std::make_unique( @@ -135,8 +137,8 @@ GateDigitizerPileupActor::GetPileupWindowForVolume( window.digis, fOutputDigiCollection, filler_out_attributes); // Store the PileupWindow in the map and return a reference. - l.fVolumePileupWindows[vol_hash] = std::move(window); - return l.fVolumePileupWindows[vol_hash]; + windows[vol_hash] = std::move(window); + return windows[vol_hash]; } } diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index 608fb4da8..7ce78e7bb 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -66,7 +66,10 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { std::unique_ptr fillerOut; }; - PileupWindow &GetPileupWindowForVolume(GateUniqueVolumeID::Pointer *volume); + PileupWindow & + GetPileupWindowForCurrentVolume(GateUniqueVolumeID::Pointer *volume, + std::map &windows); + void ProcessPileupWindow(PileupWindow &window); struct threadLocalT { From cbe871b5f3d8337e6883e69c7fdb932e3d75da40 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 21 Dec 2025 18:54:30 +0000 Subject: [PATCH 11/26] improve use of optional --- .../opengate_lib/digitizer/GateDigitizerPileupActor.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 0e39506bc..c06dd126b 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -104,7 +104,7 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( const auto vol_hash = volume->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth); - // Look up the window based + // Look up the window based on volume hash auto it = windows.find(vol_hash); if (it != windows.end()) { // Return a reference to the existing PileupWindow object for the volume. @@ -156,7 +156,7 @@ void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { while (!window.digiIter.IsAtEnd()) { const auto current_edep = *l.edep; // Remember the value and index of the highest deposited energy so far. - if (!highest_edep || current_edep > highest_edep) { + if (!highest_edep.has_value() || current_edep > highest_edep.value()) { highest_edep = current_edep; highest_edep_index = window.digiIter.fIndex; } From c7bcf44401c0b8ee2779bda0a316a180095ace6f Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Tue, 23 Dec 2025 12:39:03 +0000 Subject: [PATCH 12/26] modify pile-up test to contain non-monotonously increasing GlobalTime --- opengate/tests/src/test097_pileup.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py index b325e6c92..c25032fa0 100755 --- a/opengate/tests/src/test097_pileup.py +++ b/opengate/tests/src/test097_pileup.py @@ -70,13 +70,21 @@ crystal.material = "LYSO" crystal.color = green - source = sim.add_source("GenericSource", "b2b") - source.particle = "back_to_back" - source.activity = 10 * 1e6 * Bq - source.position.type = "sphere" - source.position.radius = 0.5 * 1e-15 * mm - source.direction.theta = [90 * deg, 90 * deg] - source.direction.phi = [0, 360 * deg] + source1 = sim.add_source("GenericSource", "b2b_1") + source1.particle = "back_to_back" + source1.activity = 5 * 1e6 * Bq + source1.position.type = "point" + source1.position.translation = [100 * mm, 0, 0] + source1.direction.theta = [90 * deg, 90 * deg] + source1.direction.phi = [0, 360 * deg] + + source2 = sim.add_source("GenericSource", "b2b_2") + source2.particle = "back_to_back" + source2.activity = 5 * 1e6 * Bq + source1.position.translation = [-100 * mm, 0, 0] + source2.position.type = "point" + source2.direction.theta = [90 * deg, 90 * deg] + source2.direction.phi = [0, 360 * deg] # Physics sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" From 663f896903bc9fcc0c17fc4a3e140efdda28797c Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Tue, 23 Dec 2025 12:39:48 +0000 Subject: [PATCH 13/26] add time sorting in pile-up actor --- .../digitizer/GateDigitizerPileupActor.cpp | 78 ++++++----- .../digitizer/GateDigitizerPileupActor.h | 5 + .../opengate_lib/digitizer/GateTimeSorter.cpp | 131 ++++++++++++++++++ .../opengate_lib/digitizer/GateTimeSorter.h | 60 ++++++++ 4 files changed, 242 insertions(+), 32 deletions(-) create mode 100644 core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp create mode 100644 core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index c06dd126b..3289e4296 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -50,44 +50,26 @@ void GateDigitizerPileupActor::DigitInitialize( auto &lr = fThreadLocalVDigitizerData.Get(); auto &l = fThreadLocalData.Get(); - lr.fInputIter.TrackAttribute("GlobalTime", &l.time); - lr.fInputIter.TrackAttribute("PreStepUniqueVolumeID", &l.volID); + l.fTimeSortedDigis = + GateDigiCollectionManager::GetInstance()->NewDigiCollection(GetName() + + "_sorted"); + l.fTimeSortedDigis->InitDigiAttributesFromCopy(fInputDigiCollection); + l.fTimeSortedDigiIterator = l.fTimeSortedDigis->NewIterator(); + l.fTimeSortedDigiIterator.TrackAttribute("GlobalTime", &l.time); + l.fTimeSortedDigiIterator.TrackAttribute("PreStepUniqueVolumeID", &l.volID); + l.fTimeSorter.Init(GetName(), fInputDigiCollection, l.fTimeSortedDigis); } void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { - auto &lr = fThreadLocalVDigitizerData.Get(); auto &l = fThreadLocalData.Get(); - auto &inputIter = lr.fInputIter; - - inputIter.GoToBegin(); - while (!inputIter.IsAtEnd()) { - // Look up or create the pile-up window object for the volume to which the - // current digi belongs. - auto &window = - GetPileupWindowForCurrentVolume(l.volID, l.fVolumePileupWindows); - - const auto current_time = *l.time; - if (window.digis->GetSize() == 0) { - // The window has no digis yet: make the window start at the time of the - // current digi. - window.startTime = current_time; - } else if (current_time - window.startTime > fTimeWindow) { - // The current digi is beyond the time window: process the digis that are - // currently in the window, then make the window start at the time of the - // current digi. - ProcessPileupWindow(window); - window.startTime = current_time; - } - - // Add the current digi to the window. - window.fillerIn->Fill(inputIter.fIndex); - - inputIter++; - } + l.fTimeSorter.Process(); + ProcessTimeSortedDigis(); } void GateDigitizerPileupActor::EndOfRunAction(const G4Run *) { auto &l = fThreadLocalData.Get(); + l.fTimeSorter.Flush(); + ProcessTimeSortedDigis(); for (auto &[_vol_hash, window] : l.fVolumePileupWindows) { ProcessPileupWindow(window); } @@ -123,10 +105,11 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( window.digiIter = window.digis->NewIterator(); window.digiIter.TrackAttribute("TotalEnergyDeposit", &fThreadLocalData.Get().edep); - // Create a filler to copy all digi attributes from the input collection + // Create a filler to copy all digi attributes from the sorted collection // into the collection of the window. + auto &l = fThreadLocalData.Get(); window.fillerIn = std::make_unique( - fInputDigiCollection, window.digis, + l.fTimeSortedDigis, window.digis, window.digis->GetDigiAttributeNames()); // Create a filler to copy digi attributes from the collection of the window // to the output collection (used for the digis that will result from @@ -142,6 +125,37 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( } } +void GateDigitizerPileupActor::ProcessTimeSortedDigis() { + auto &l = fThreadLocalData.Get(); + auto &iter = l.fTimeSortedDigiIterator; + iter.GoToBegin(); + while (!iter.IsAtEnd()) { + // Look up or create the pile-up window object for the volume to which the + // current digi belongs. + auto &window = + GetPileupWindowForCurrentVolume(l.volID, l.fVolumePileupWindows); + + const auto current_time = *l.time; + if (window.digis->GetSize() == 0) { + // The window has no digis yet: make the window start at the time of the + // current digi. + window.startTime = current_time; + } else if (current_time - window.startTime > fTimeWindow) { + // The current digi is beyond the time window: process the digis that are + // currently in the window, then make the window start at the time of the + // current digi. + ProcessPileupWindow(window); + window.startTime = current_time; + } + + // Add the current digi to the window. + window.fillerIn->Fill(iter.fIndex); + + iter++; + } + l.fTimeSortedDigis->SetBeginOfEventIndex(iter.fIndex); +} + void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { // This function simulates pile-up by combining the digis in the given window // into one digi. diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index 7ce78e7bb..d58338938 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -8,6 +8,7 @@ #ifndef GateDigitizerPileupActor_h #define GateDigitizerPileupActor_h +#include "GateTimeSorter.h" #include "GateVDigitizerWithOutputActor.h" #include #include @@ -70,6 +71,7 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { GetPileupWindowForCurrentVolume(GateUniqueVolumeID::Pointer *volume, std::map &windows); + void ProcessTimeSortedDigis(); void ProcessPileupWindow(PileupWindow &window); struct threadLocalT { @@ -77,6 +79,9 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { double *time; double *edep; + GateTimeSorter fTimeSorter; + GateDigiCollection *fTimeSortedDigis; + GateDigiCollectionIterator fTimeSortedDigiIterator; std::map fVolumePileupWindows; }; G4Cache fThreadLocalData; diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp new file mode 100644 index 000000000..53356ab57 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp @@ -0,0 +1,131 @@ +#include "GateTimeSorter.h" +#include "GateDigiCollection.h" +#include "GateDigiCollectionManager.h" +#include "GateHelpersDigitizer.h" +#include + +void GateTimeSorter::Init(const std::string &name, GateDigiCollection *input, + GateDigiCollection *output) { + if (fInitialized) { + Fatal("Init() cannot be called more than once."); + } + + fInputIter = input->NewIterator(); + fInputIter.TrackAttribute("GlobalTime", &fTime); + + auto *manager = GateDigiCollectionManager::GetInstance(); + const auto attribute_names = input->GetDigiAttributeNames(); + + fCurrentStorage = std::make_unique(); + fCurrentStorage->digis = manager->NewDigiCollection(name + "_collectionA"); + fCurrentStorage->digis->InitDigiAttributesFromCopy(input); + fCurrentStorage->fillerIn = std::make_unique( + input, fCurrentStorage->digis, attribute_names); + fCurrentStorage->fillerOut = std::make_unique( + fCurrentStorage->digis, output, attribute_names); + + fFutureStorage = std::make_unique(); + fFutureStorage->digis = manager->NewDigiCollection(name + "_collectionB"); + fFutureStorage->digis->InitDigiAttributesFromCopy(input); + fFutureStorage->fillerIn = std::make_unique( + input, fFutureStorage->digis, attribute_names); + fFutureStorage->fillerOut = std::make_unique( + fFutureStorage->digis, output, attribute_names); + + fCurrentStorage->fillerSwap = std::make_unique( + fCurrentStorage->digis, fFutureStorage->digis, attribute_names); + fFutureStorage->fillerSwap = std::make_unique( + fFutureStorage->digis, fCurrentStorage->digis, attribute_names); + + fInitialized = true; +} + +void GateTimeSorter::SetDelay(double delay) { + if (fProcessingStarted) { + Fatal("SetDelay() cannot be called after Process() has been called."); + } + fDelay = delay; +} + +void GateTimeSorter::SetMaxSize(size_t maxSize) { + if (fProcessingStarted) { + Fatal("SetMaxSize() cannot be called after Process() has been called."); + } + fMaxSize = maxSize; +} + +void GateTimeSorter::Process() { + if (!fInitialized) { + Fatal("Process() called before Init(). " + "Init() must be called first to initialize the time sorter."); + } + if (fFlushed) { + Fatal("Process called after Flush(). The time sorter must not be used " + "after it was flushed."); + } + fProcessingStarted = true; + + auto &iter = fInputIter; + auto &sortedIndices = fCurrentStorage->sortedIndices; + + iter.GoToBegin(); + while (!iter.IsAtEnd()) { + const size_t digiIndex = fCurrentStorage->digis->GetSize(); + const double digiTime = *fTime; + if (fMostRecentTimeDeparted.has_value() && + (digiTime < *fMostRecentTimeDeparted)) { + // TODO warn that digi was dropped to be able to guarantee time + // monotonicity and that delay parameter should be increased. + } else { + fCurrentStorage->fillerIn->Fill(iter.fIndex); + sortedIndices.push({digiIndex, digiTime}); + + if (!fMostRecentTimeArrived || (digiTime > *fMostRecentTimeArrived)) { + fMostRecentTimeArrived = digiTime; + } + } + iter++; + } + + if (!sortedIndices.empty()) { + size_t num_filled = 0; + double timespan = *fMostRecentTimeArrived - sortedIndices.top().time; + while (!sortedIndices.empty() && timespan > fDelay) { + fCurrentStorage->fillerOut->Fill(sortedIndices.top().index); + fMostRecentTimeDeparted = sortedIndices.top().time; + sortedIndices.pop(); + if (!sortedIndices.empty()) { + timespan = *fMostRecentTimeArrived - sortedIndices.top().time; + } + num_filled++; + } + } + + if (fCurrentStorage->digis->GetSize() > fMaxSize) { + Prune(); + } +} + +void GateTimeSorter::Flush() { + auto &sortedIndices = fCurrentStorage->sortedIndices; + while (sortedIndices.size() > 0) { + fCurrentStorage->fillerOut->Fill(sortedIndices.top().index); + sortedIndices.pop(); + } + Prune(); + fFlushed = true; +} + +void GateTimeSorter::Prune() { + auto &sortedIndices = fCurrentStorage->sortedIndices; + while (!sortedIndices.empty()) { + const auto timed_index = sortedIndices.top(); + sortedIndices.pop(); + const size_t digiIndex = fFutureStorage->digis->GetSize(); + const double digiTime = timed_index.time; + fCurrentStorage->fillerSwap->Fill(timed_index.index); + fFutureStorage->sortedIndices.push({digiIndex, digiTime}); + } + fCurrentStorage->digis->Clear(); + std::swap(fCurrentStorage, fFutureStorage); +} diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h new file mode 100644 index 000000000..b03243492 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h @@ -0,0 +1,60 @@ +#ifndef GateTimeSorter_h +#define GateTimeSorter_h + +#include "GateDigiCollectionIterator.h" +#include +#include +#include + +class GateDigiCollection; +class GateDigiAttributesFiller; + +class GateTimeSorter { +public: + GateTimeSorter() = default; + + void Init(const std::string &name, GateDigiCollection *input, + GateDigiCollection *output); + + void SetDelay(double delay); + void SetMaxSize(size_t size); + + void Process(); + void Flush(); + +private: + void Prune(); + + struct TimedDigiIndex { + size_t index; + double time; + + bool operator>(const TimedDigiIndex &other) const { + return time > other.time; + } + }; + + struct TimeSortedStorage { + GateDigiCollection *digis; + std::priority_queue, + std::greater> + sortedIndices; + std::unique_ptr fillerIn; + std::unique_ptr fillerOut; + std::unique_ptr fillerSwap; + }; + + GateDigiCollectionIterator fInputIter; + double *fTime; + double fDelay{1000.0}; // nanoseconds + size_t fMaxSize{100'000}; + bool fInitialized{false}; + bool fProcessingStarted{false}; + bool fFlushed{false}; + std::optional fMostRecentTimeArrived; + std::optional fMostRecentTimeDeparted; + std::unique_ptr fCurrentStorage; + std::unique_ptr fFutureStorage; +}; + +#endif // GateTimeSorter_h From 59d7cc1a22663e6dd0ecd0c88a27d73fb4f45d69 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Tue, 23 Dec 2025 20:09:52 +0000 Subject: [PATCH 14/26] encapsulate Filler and Iterator inside GateTimeSorter --- .../digitizer/GateDigitizerPileupActor.cpp | 20 ++---- .../digitizer/GateDigitizerPileupActor.h | 2 - .../opengate_lib/digitizer/GateTimeSorter.cpp | 71 ++++++++++++------- .../opengate_lib/digitizer/GateTimeSorter.h | 21 ++++-- 4 files changed, 67 insertions(+), 47 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 3289e4296..2ea6fc42f 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -50,14 +50,10 @@ void GateDigitizerPileupActor::DigitInitialize( auto &lr = fThreadLocalVDigitizerData.Get(); auto &l = fThreadLocalData.Get(); - l.fTimeSortedDigis = - GateDigiCollectionManager::GetInstance()->NewDigiCollection(GetName() + - "_sorted"); - l.fTimeSortedDigis->InitDigiAttributesFromCopy(fInputDigiCollection); - l.fTimeSortedDigiIterator = l.fTimeSortedDigis->NewIterator(); - l.fTimeSortedDigiIterator.TrackAttribute("GlobalTime", &l.time); - l.fTimeSortedDigiIterator.TrackAttribute("PreStepUniqueVolumeID", &l.volID); - l.fTimeSorter.Init(GetName(), fInputDigiCollection, l.fTimeSortedDigis); + l.fTimeSorter.Init(fInputDigiCollection); + l.fTimeSorter.OutputIterator().TrackAttribute("GlobalTime", &l.time); + l.fTimeSorter.OutputIterator().TrackAttribute("PreStepUniqueVolumeID", + &l.volID); } void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { @@ -108,9 +104,7 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( // Create a filler to copy all digi attributes from the sorted collection // into the collection of the window. auto &l = fThreadLocalData.Get(); - window.fillerIn = std::make_unique( - l.fTimeSortedDigis, window.digis, - window.digis->GetDigiAttributeNames()); + window.fillerIn = l.fTimeSorter.CreateFiller(window.digis); // Create a filler to copy digi attributes from the collection of the window // to the output collection (used for the digis that will result from // pile-up). @@ -127,7 +121,7 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( void GateDigitizerPileupActor::ProcessTimeSortedDigis() { auto &l = fThreadLocalData.Get(); - auto &iter = l.fTimeSortedDigiIterator; + auto &iter = l.fTimeSorter.OutputIterator(); iter.GoToBegin(); while (!iter.IsAtEnd()) { // Look up or create the pile-up window object for the volume to which the @@ -153,7 +147,7 @@ void GateDigitizerPileupActor::ProcessTimeSortedDigis() { iter++; } - l.fTimeSortedDigis->SetBeginOfEventIndex(iter.fIndex); + l.fTimeSorter.MarkOutputAsProcessed(); } void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index d58338938..fcd7ec17d 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -80,8 +80,6 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { double *edep; GateTimeSorter fTimeSorter; - GateDigiCollection *fTimeSortedDigis; - GateDigiCollectionIterator fTimeSortedDigiIterator; std::map fVolumePileupWindows; }; G4Cache fThreadLocalData; diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp index 53356ab57..185b6bb3d 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp @@ -4,47 +4,45 @@ #include "GateHelpersDigitizer.h" #include -void GateTimeSorter::Init(const std::string &name, GateDigiCollection *input, - GateDigiCollection *output) { - if (fInitialized) { - Fatal("Init() cannot be called more than once."); - } - - fInputIter = input->NewIterator(); +void GateTimeSorter::Init(GateDigiCollection *input) { + fInputCollection = input; + fInputIter = fInputCollection->NewIterator(); fInputIter.TrackAttribute("GlobalTime", &fTime); auto *manager = GateDigiCollectionManager::GetInstance(); - const auto attribute_names = input->GetDigiAttributeNames(); + const auto attribute_names = fInputCollection->GetDigiAttributeNames(); + + fOutputCollection = + manager->NewDigiCollection(fInputCollection->GetName() + "_sorted"); + fOutputCollection->InitDigiAttributesFromCopy(fInputCollection); + fOutputIter = fOutputCollection->NewIterator(); fCurrentStorage = std::make_unique(); - fCurrentStorage->digis = manager->NewDigiCollection(name + "_collectionA"); - fCurrentStorage->digis->InitDigiAttributesFromCopy(input); + fCurrentStorage->digis = + manager->NewDigiCollection(fInputCollection->GetName() + "_temporaryA"); + fCurrentStorage->digis->InitDigiAttributesFromCopy(fInputCollection); fCurrentStorage->fillerIn = std::make_unique( - input, fCurrentStorage->digis, attribute_names); + fInputCollection, fCurrentStorage->digis, attribute_names); fCurrentStorage->fillerOut = std::make_unique( - fCurrentStorage->digis, output, attribute_names); + fCurrentStorage->digis, fOutputCollection, attribute_names); fFutureStorage = std::make_unique(); - fFutureStorage->digis = manager->NewDigiCollection(name + "_collectionB"); - fFutureStorage->digis->InitDigiAttributesFromCopy(input); + fFutureStorage->digis = + manager->NewDigiCollection(fInputCollection->GetName() + "_temporaryB"); + fFutureStorage->digis->InitDigiAttributesFromCopy(fInputCollection); fFutureStorage->fillerIn = std::make_unique( - input, fFutureStorage->digis, attribute_names); + fInputCollection, fFutureStorage->digis, attribute_names); fFutureStorage->fillerOut = std::make_unique( - fFutureStorage->digis, output, attribute_names); - - fCurrentStorage->fillerSwap = std::make_unique( - fCurrentStorage->digis, fFutureStorage->digis, attribute_names); - fFutureStorage->fillerSwap = std::make_unique( - fFutureStorage->digis, fCurrentStorage->digis, attribute_names); + fFutureStorage->digis, fOutputCollection, attribute_names); fInitialized = true; } -void GateTimeSorter::SetDelay(double delay) { +void GateTimeSorter::SetSortingWindow(double duration) { if (fProcessingStarted) { Fatal("SetDelay() cannot be called after Process() has been called."); } - fDelay = delay; + fSortingWindow = duration; } void GateTimeSorter::SetMaxSize(size_t maxSize) { @@ -54,11 +52,21 @@ void GateTimeSorter::SetMaxSize(size_t maxSize) { fMaxSize = maxSize; } -void GateTimeSorter::Process() { +std::unique_ptr +GateTimeSorter::CreateFiller(GateDigiCollection *destination) { if (!fInitialized) { - Fatal("Process() called before Init(). " - "Init() must be called first to initialize the time sorter."); + Fatal("CreateFiller() called before Init()."); } + return std::make_unique( + fOutputCollection, destination, + fOutputCollection->GetDigiAttributeNames()); +} + +GateDigiCollection::Iterator &GateTimeSorter::OutputIterator() { + return fOutputIter; +} + +void GateTimeSorter::Process() { if (fFlushed) { Fatal("Process called after Flush(). The time sorter must not be used " "after it was flushed."); @@ -90,7 +98,7 @@ void GateTimeSorter::Process() { if (!sortedIndices.empty()) { size_t num_filled = 0; double timespan = *fMostRecentTimeArrived - sortedIndices.top().time; - while (!sortedIndices.empty() && timespan > fDelay) { + while (!sortedIndices.empty() && timespan > fSortingWindow) { fCurrentStorage->fillerOut->Fill(sortedIndices.top().index); fMostRecentTimeDeparted = sortedIndices.top().time; sortedIndices.pop(); @@ -106,6 +114,15 @@ void GateTimeSorter::Process() { } } +void GateTimeSorter::MarkOutputAsProcessed() { + if (fOutputCollection->GetSize() <= fMaxSize) { + fOutputCollection->SetBeginOfEventIndex(fOutputIter.fIndex); + } else { + fOutputCollection->Clear(); + fOutputIter.Reset(); + } +} + void GateTimeSorter::Flush() { auto &sortedIndices = fCurrentStorage->sortedIndices; while (sortedIndices.size() > 0) { diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h index b03243492..c19687a24 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h @@ -13,13 +13,16 @@ class GateTimeSorter { public: GateTimeSorter() = default; - void Init(const std::string &name, GateDigiCollection *input, - GateDigiCollection *output); + void Init(GateDigiCollection *input); - void SetDelay(double delay); + void SetSortingWindow(double duration); void SetMaxSize(size_t size); + std::unique_ptr + CreateFiller(GateDigiCollection *destination); + GateDigiCollection::Iterator &OutputIterator(); void Process(); + void MarkOutputAsProcessed(); void Flush(); private: @@ -44,15 +47,23 @@ class GateTimeSorter { std::unique_ptr fillerSwap; }; + double fSortingWindow{1000.0}; // nanoseconds + size_t fMaxSize{100'000}; + + GateDigiCollection *fInputCollection; GateDigiCollectionIterator fInputIter; + + GateDigiCollection *fOutputCollection; + GateDigiCollectionIterator fOutputIter; + double *fTime; - double fDelay{1000.0}; // nanoseconds - size_t fMaxSize{100'000}; + bool fInitialized{false}; bool fProcessingStarted{false}; bool fFlushed{false}; std::optional fMostRecentTimeArrived; std::optional fMostRecentTimeDeparted; + std::unique_ptr fCurrentStorage; std::unique_ptr fFutureStorage; }; From 896bc8cb97944b881535812f69c9fadcab695e1f Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Wed, 24 Dec 2025 07:53:07 +0000 Subject: [PATCH 15/26] add small improvements and comments to GateTimeSorter --- .../opengate_lib/digitizer/GateTimeSorter.cpp | 134 ++++++++++++++---- .../opengate_lib/digitizer/GateTimeSorter.h | 4 + 2 files changed, 107 insertions(+), 31 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp index 185b6bb3d..645393780 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp @@ -4,7 +4,34 @@ #include "GateHelpersDigitizer.h" #include +GateTimeSorter::TimeSortedStorage::TimeSortedStorage( + GateDigiCollection *input, GateDigiCollection *output, + const std::string &name_suffix) { + + auto *manager = GateDigiCollectionManager::GetInstance(); + const auto attribute_names = input->GetDigiAttributeNames(); + + // GateDigiCollection for temporary storage + digis = manager->NewDigiCollection(input->GetName() + "_" + name_suffix); + digis->InitDigiAttributesFromCopy(input); + + // Filler to copy from input collection to temporary collection + fillerIn = + std::make_unique(input, digis, attribute_names); + + // Filler to copy from temporary to output collection + fillerOut = std::make_unique(digis, output, + attribute_names); + + // The GateDigiCollection contains the digis in the order in which they were + // added using fillerIn, but the actual sorting happens in the priority queue + // sortedIndices. +} + void GateTimeSorter::Init(GateDigiCollection *input) { + + // Create an iterator for the input collection, tracking the GlobalTime of + // digis to be able to sort them in time. fInputCollection = input; fInputIter = fInputCollection->NewIterator(); fInputIter.TrackAttribute("GlobalTime", &fTime); @@ -12,33 +39,36 @@ void GateTimeSorter::Init(GateDigiCollection *input) { auto *manager = GateDigiCollectionManager::GetInstance(); const auto attribute_names = fInputCollection->GetDigiAttributeNames(); + // Create an output collection that will receive the time-sorted digis, and an + // iterator for the output collection. fOutputCollection = manager->NewDigiCollection(fInputCollection->GetName() + "_sorted"); fOutputCollection->InitDigiAttributesFromCopy(fInputCollection); fOutputIter = fOutputCollection->NewIterator(); - fCurrentStorage = std::make_unique(); - fCurrentStorage->digis = - manager->NewDigiCollection(fInputCollection->GetName() + "_temporaryA"); - fCurrentStorage->digis->InitDigiAttributesFromCopy(fInputCollection); - fCurrentStorage->fillerIn = std::make_unique( - fInputCollection, fCurrentStorage->digis, attribute_names); - fCurrentStorage->fillerOut = std::make_unique( - fCurrentStorage->digis, fOutputCollection, attribute_names); - - fFutureStorage = std::make_unique(); - fFutureStorage->digis = - manager->NewDigiCollection(fInputCollection->GetName() + "_temporaryB"); - fFutureStorage->digis->InitDigiAttributesFromCopy(fInputCollection); - fFutureStorage->fillerIn = std::make_unique( - fInputCollection, fFutureStorage->digis, attribute_names); - fFutureStorage->fillerOut = std::make_unique( - fFutureStorage->digis, fOutputCollection, attribute_names); + // Create a TimeSortedStorage object for sorting digis. + fCurrentStorage = std::make_unique( + fInputCollection, fOutputCollection, "temporaryA"); + + // Create a second TimeSortedStorage object that will be used later, when the + // first one's memory needs to be freed. + fFutureStorage = std::make_unique( + fInputCollection, fOutputCollection, "temporaryB"); fInitialized = true; } void GateTimeSorter::SetSortingWindow(double duration) { + // The sorting window, specified in nanoseconds, is the difference in + // GlobalTime between the oldest and newest digi being stored in the + // TimeSortedStorage object, before they are transferred to the output + // collection. Since the GlobalTime of incoming digis is not guaranteed to + // increase monotonically, the sorting window must be large enough to ensure + // that all digis can be properly time-sorted. + // Non-monotonicity of GlobalTime can be caused by the time difference between + // a particle's generation and its interaction with a detector. It can also be + // caused by a DigitizerBlurringActor with blur_attribute "GlobalTime". + if (fProcessingStarted) { Fatal("SetDelay() cannot be called after Process() has been called."); } @@ -54,6 +84,9 @@ void GateTimeSorter::SetMaxSize(size_t maxSize) { std::unique_ptr GateTimeSorter::CreateFiller(GateDigiCollection *destination) { + // Creates a GateDigiAttributesFiller from the TimeSorter's output collection + // into the given destination GateDigiCollection. + if (!fInitialized) { Fatal("CreateFiller() called before Init()."); } @@ -63,10 +96,15 @@ GateTimeSorter::CreateFiller(GateDigiCollection *destination) { } GateDigiCollection::Iterator &GateTimeSorter::OutputIterator() { + // Provides access to the GateTimeSorter's output iterator. return fOutputIter; } void GateTimeSorter::Process() { + // Processes all digis from the input collection by copying and sorting them + // according to GlobalTime. Next, copies the oldest sorted digis to the output + // collection. + if (fFlushed) { Fatal("Process called after Flush(). The time sorter must not be used " "after it was flushed."); @@ -76,18 +114,32 @@ void GateTimeSorter::Process() { auto &iter = fInputIter; auto &sortedIndices = fCurrentStorage->sortedIndices; + // Iterate over the input collection and sort the digis. iter.GoToBegin(); while (!iter.IsAtEnd()) { const size_t digiIndex = fCurrentStorage->digis->GetSize(); const double digiTime = *fTime; if (fMostRecentTimeDeparted.has_value() && (digiTime < *fMostRecentTimeDeparted)) { - // TODO warn that digi was dropped to be able to guarantee time - // monotonicity and that delay parameter should be increased. + // The digi is dropped, in order to be able to guarantee monotonous + // GlobalTime in the output collection. + if (!fSortingWindowWarningIssued) { + std::cout << "The digis in " << fInputCollection->GetName() + << " have non-monotonicities in the GlobalTime attribute " + "that exceed " + "the sorting window (" + << fSortingWindow + << " ns). Please increase the sorting window to avoid " + "dropped digis"; + fSortingWindowWarningIssued = true; + } } else { + // Copy the digi into the temporary digi collection. fCurrentStorage->fillerIn->Fill(iter.fIndex); + // Keep a time-sorted list of indices into the temporary digi collection. sortedIndices.push({digiIndex, digiTime}); + // Keep track of the highest GlobalTime observed so far. if (!fMostRecentTimeArrived || (digiTime > *fMostRecentTimeArrived)) { fMostRecentTimeArrived = digiTime; } @@ -95,26 +147,33 @@ void GateTimeSorter::Process() { iter++; } - if (!sortedIndices.empty()) { - size_t num_filled = 0; - double timespan = *fMostRecentTimeArrived - sortedIndices.top().time; - while (!sortedIndices.empty() && timespan > fSortingWindow) { - fCurrentStorage->fillerOut->Fill(sortedIndices.top().index); - fMostRecentTimeDeparted = sortedIndices.top().time; - sortedIndices.pop(); - if (!sortedIndices.empty()) { - timespan = *fMostRecentTimeArrived - sortedIndices.top().time; - } - num_filled++; - } + // Copy the oldest digis from the sorted temporary storage into the output + // collection. Continue as long as the newest digi is at least fSortingWindow + // more recent than the oldest digi. + while ( + !sortedIndices.empty() && + (*fMostRecentTimeArrived - sortedIndices.top().time > fSortingWindow)) { + // Copy oldest digi into the output collection. + fCurrentStorage->fillerOut->Fill(sortedIndices.top().index); + // Keep track of the GlobalTime of the last digi that was copied. + fMostRecentTimeDeparted = sortedIndices.top().time; + // Remove the time-sorted index of the digi. + sortedIndices.pop(); } + // The temporary digi collection keeps growing as more digis are processed. + // The digis that have already been copied to the output must be removed once + // in a while to limit memory usage. if (fCurrentStorage->digis->GetSize() > fMaxSize) { Prune(); } } void GateTimeSorter::MarkOutputAsProcessed() { + // Modifies the start of event index, to ensure that future use of the + // iterator starts with the output digis that have not been processed yet. + // Once in a while, the output collection is cleared to reduce memory + // consumption. if (fOutputCollection->GetSize() <= fMaxSize) { fOutputCollection->SetBeginOfEventIndex(fOutputIter.fIndex); } else { @@ -124,6 +183,11 @@ void GateTimeSorter::MarkOutputAsProcessed() { } void GateTimeSorter::Flush() { + // Copies all remaining sorted digis into the output collection. + // This method is intended to be called once at the end of time sorting, when + // it is known that no more digis will be processed from the input. + // As a consequence, the sorting window does not have to be taken into account + // while flushing. auto &sortedIndices = fCurrentStorage->sortedIndices; while (sortedIndices.size() > 0) { fCurrentStorage->fillerOut->Fill(sortedIndices.top().index); @@ -134,6 +198,14 @@ void GateTimeSorter::Flush() { } void GateTimeSorter::Prune() { + // Frees memory that is currently occupied by digis that have already been + // copied into the output collection: + // 1. The digis that have not been copied to the output yet, are copied into + // the second instance of TimeSortedStorage. + // 2. The digi collection in the first instance of TimeSortedStorage is + // cleared. + // 3. The two instances of TimeSortedStorage are swapped. + auto &sortedIndices = fCurrentStorage->sortedIndices; while (!sortedIndices.empty()) { const auto timed_index = sortedIndices.top(); diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h index c19687a24..8760c7d0d 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h @@ -38,6 +38,9 @@ class GateTimeSorter { }; struct TimeSortedStorage { + TimeSortedStorage(GateDigiCollection *input, GateDigiCollection *output, + const std::string &name_suffix); + GateDigiCollection *digis; std::priority_queue, std::greater> @@ -61,6 +64,7 @@ class GateTimeSorter { bool fInitialized{false}; bool fProcessingStarted{false}; bool fFlushed{false}; + bool fSortingWindowWarningIssued{false}; std::optional fMostRecentTimeArrived; std::optional fMostRecentTimeDeparted; From 5f6d239c59d63b0017320f727b42d2f9d01bfcf0 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 07:32:22 +0000 Subject: [PATCH 16/26] rename parameter to pile-up time --- .../digitizer/GateDigitizerPileupActor.cpp | 15 +++++++-------- .../digitizer/GateDigitizerPileupActor.h | 2 +- opengate/tests/src/test097_pileup.py | 7 +++++-- opengate/tests/src/test097_pileup_helpers.py | 15 +++++++++++---- 4 files changed, 24 insertions(+), 15 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 2ea6fc42f..6753cd03c 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -22,10 +22,10 @@ GateDigitizerPileupActor::~GateDigitizerPileupActor() = default; void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); - // Get time window parameter in ns. - fTimeWindow = 0.0; // default value, no pile-up - if (py::len(user_info) > 0 && user_info.contains("time_window")) { - fTimeWindow = DictGetDouble(user_info, "time_window"); + // Get pile-up time parameter in ns. + fPileupTime = 0.0; // default value, no pile-up + if (py::len(user_info) > 0 && user_info.contains("pileup_time")) { + fPileupTime = DictGetDouble(user_info, "pileup_time"); } fGroupVolumeDepth = -1; fInputDigiCollectionName = DictGetStr(user_info, "input_digi_collection"); @@ -91,6 +91,7 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( // A PileupWindow object does not yet exist for this volume: create one. PileupWindow window; const auto vol_id = volume->get()->GetIdUpToDepth(fGroupVolumeDepth); + auto &l = fThreadLocalData.Get(); // Create a GateDigiCollection for this volume, as a temporary storage for // digis that belong to the same time window (the name must be unique). window.digis = GateDigiCollectionManager::GetInstance()->NewDigiCollection( @@ -99,11 +100,9 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( // Create an iterator to be used when digis will be combined into one digi, // due to pile-up. window.digiIter = window.digis->NewIterator(); - window.digiIter.TrackAttribute("TotalEnergyDeposit", - &fThreadLocalData.Get().edep); + window.digiIter.TrackAttribute("TotalEnergyDeposit", &l.edep); // Create a filler to copy all digi attributes from the sorted collection // into the collection of the window. - auto &l = fThreadLocalData.Get(); window.fillerIn = l.fTimeSorter.CreateFiller(window.digis); // Create a filler to copy digi attributes from the collection of the window // to the output collection (used for the digis that will result from @@ -134,7 +133,7 @@ void GateDigitizerPileupActor::ProcessTimeSortedDigis() { // The window has no digis yet: make the window start at the time of the // current digi. window.startTime = current_time; - } else if (current_time - window.startTime > fTimeWindow) { + } else if (current_time - window.startTime > fPileupTime) { // The current digi is beyond the time window: process the digis that are // currently in the window, then make the window start at the time of the // current digi. diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index fcd7ec17d..afaa4098f 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -44,8 +44,8 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { const std::vector &attributes_not_in_filler) override; // User parameters + double fPileupTime; int fGroupVolumeDepth; - double fTimeWindow; // Output attribute pointer GateVDigiAttribute *fOutputEdepAttribute{}; diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py index c25032fa0..d2b24f95d 100755 --- a/opengate/tests/src/test097_pileup.py +++ b/opengate/tests/src/test097_pileup.py @@ -22,6 +22,7 @@ # Units mm = gate.g4_units.mm sec = gate.g4_units.s + ns = gate.g4_units.ns keV = gate.g4_units.keV Bq = gate.g4_units.Bq gcm3 = gate.g4_units.g_cm3 @@ -118,7 +119,9 @@ pu.authorize_repeated_volumes = True pu.input_digi_collection = sc.name pu.output_filename = sc.output_filename - pu.time_window = 1000.0 + pu.pileup_time = 2000.0 * ns + pu.clear_every = 1e3 + # pu.skip_attributes = ["PreStepUniqueVolumeID"] # Timing sim.run_timing_intervals = [[0, 0.001 * sec]] @@ -130,6 +133,6 @@ sc.output_filename, "Singles_before_pileup", "Singles_after_pileup", - pu.time_window, + pu.pileup_time, ) utility.test_ok(all_match) diff --git a/opengate/tests/src/test097_pileup_helpers.py b/opengate/tests/src/test097_pileup_helpers.py index ab5933cec..1ecb88f77 100644 --- a/opengate/tests/src/test097_pileup_helpers.py +++ b/opengate/tests/src/test097_pileup_helpers.py @@ -1,5 +1,6 @@ import pandas as pd import uproot +import numpy as np def pileup(singles_before_pileup: pd.DataFrame, time_window: float): @@ -112,9 +113,15 @@ def check_gate_pileup( continue expected_values = expected_singles[attr].values actual_values = actual_singles[attr].values - if not all(expected_values == actual_values): - print(f"Volume {volume_id}: Attribute {attr} does not match") - all_match = False - break + if np.issubdtype(expected_values.dtype, np.floating): + if not np.allclose(expected_values, actual_values, rtol=1e-9): + print(f"Volume {volume_id}: Attribute {attr} does not match") + all_match = False + break + else: + if not all(expected_values == actual_values): + print(f"Volume {volume_id}: Attribute {attr} does not match") + all_match = False + break return all_match From 5d0b27fb2ed8cddd25525677cbda7ec354f1d976 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 07:36:29 +0000 Subject: [PATCH 17/26] fix issue in GateTimeSorter --- .../digitizer/GateDigitizerPileupActor.cpp | 1 + .../opengate_lib/digitizer/GateTimeSorter.cpp | 11 ++++++++--- .../opengate_lib/digitizer/GateTimeSorter.h | 1 - opengate/tests/src/test097_pileup.py | 3 +-- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 6753cd03c..0cf6dde43 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -54,6 +54,7 @@ void GateDigitizerPileupActor::DigitInitialize( l.fTimeSorter.OutputIterator().TrackAttribute("GlobalTime", &l.time); l.fTimeSorter.OutputIterator().TrackAttribute("PreStepUniqueVolumeID", &l.volID); + l.fTimeSorter.SetMaxSize(fClearEveryNEvents); } void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp index 645393780..bd921a7e1 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp @@ -126,8 +126,7 @@ void GateTimeSorter::Process() { if (!fSortingWindowWarningIssued) { std::cout << "The digis in " << fInputCollection->GetName() << " have non-monotonicities in the GlobalTime attribute " - "that exceed " - "the sorting window (" + "that exceed the sorting window (" << fSortingWindow << " ns). Please increase the sorting window to avoid " "dropped digis"; @@ -206,15 +205,21 @@ void GateTimeSorter::Prune() { // cleared. // 3. The two instances of TimeSortedStorage are swapped. + // Step 1 + GateDigiAttributesFiller transferFiller( + fCurrentStorage->digis, fFutureStorage->digis, + fCurrentStorage->digis->GetDigiAttributeNames()); auto &sortedIndices = fCurrentStorage->sortedIndices; while (!sortedIndices.empty()) { const auto timed_index = sortedIndices.top(); sortedIndices.pop(); const size_t digiIndex = fFutureStorage->digis->GetSize(); const double digiTime = timed_index.time; - fCurrentStorage->fillerSwap->Fill(timed_index.index); + transferFiller.Fill(timed_index.index); fFutureStorage->sortedIndices.push({digiIndex, digiTime}); } + // Step 2 fCurrentStorage->digis->Clear(); + // Step 3 std::swap(fCurrentStorage, fFutureStorage); } diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h index 8760c7d0d..841eb6693 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h @@ -47,7 +47,6 @@ class GateTimeSorter { sortedIndices; std::unique_ptr fillerIn; std::unique_ptr fillerOut; - std::unique_ptr fillerSwap; }; double fSortingWindow{1000.0}; // nanoseconds diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py index d2b24f95d..14e802bbe 100755 --- a/opengate/tests/src/test097_pileup.py +++ b/opengate/tests/src/test097_pileup.py @@ -120,8 +120,7 @@ pu.input_digi_collection = sc.name pu.output_filename = sc.output_filename pu.pileup_time = 2000.0 * ns - pu.clear_every = 1e3 - # pu.skip_attributes = ["PreStepUniqueVolumeID"] + pu.clear_every = 1e4 # Timing sim.run_timing_intervals = [[0, 0.001 * sec]] From c16d6839f65598872ce535175ac92d2278ed28f5 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 08:55:39 +0000 Subject: [PATCH 18/26] make sorting time configurable and warn when digis are dropped while sorting --- .../digitizer/GateDigitizerPileupActor.cpp | 8 +++--- .../digitizer/GateDigitizerPileupActor.h | 1 + .../opengate_lib/digitizer/GateTimeSorter.cpp | 11 ++++++-- .../opengate_lib/digitizer/GateTimeSorter.h | 1 + opengate/actors/digitizers.py | 26 +++++++++---------- 5 files changed, 29 insertions(+), 18 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 0cf6dde43..c07e23907 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -22,10 +22,11 @@ GateDigitizerPileupActor::~GateDigitizerPileupActor() = default; void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); - // Get pile-up time parameter in ns. - fPileupTime = 0.0; // default value, no pile-up if (py::len(user_info) > 0 && user_info.contains("pileup_time")) { - fPileupTime = DictGetDouble(user_info, "pileup_time"); + fPileupTime = DictGetDouble(user_info, "pileup_time"); // nanoseconds + } + if (py::len(user_info) > 0 && user_info.contains("sorting_time")) { + fSortingTime = DictGetDouble(user_info, "sorting_time"); // nanoseconds } fGroupVolumeDepth = -1; fInputDigiCollectionName = DictGetStr(user_info, "input_digi_collection"); @@ -54,6 +55,7 @@ void GateDigitizerPileupActor::DigitInitialize( l.fTimeSorter.OutputIterator().TrackAttribute("GlobalTime", &l.time); l.fTimeSorter.OutputIterator().TrackAttribute("PreStepUniqueVolumeID", &l.volID); + l.fTimeSorter.SetSortingWindow(fSortingTime); l.fTimeSorter.SetMaxSize(fClearEveryNEvents); } diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index afaa4098f..2fe3a1f05 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -45,6 +45,7 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { // User parameters double fPileupTime; + double fSortingTime; int fGroupVolumeDepth; // Output attribute pointer diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp index bd921a7e1..29b22037f 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp @@ -123,13 +123,14 @@ void GateTimeSorter::Process() { (digiTime < *fMostRecentTimeDeparted)) { // The digi is dropped, in order to be able to guarantee monotonous // GlobalTime in the output collection. + ++fNumDroppedDigi; if (!fSortingWindowWarningIssued) { std::cout << "The digis in " << fInputCollection->GetName() << " have non-monotonicities in the GlobalTime attribute " - "that exceed the sorting window (" + "that exceed the sorting time (" << fSortingWindow << " ns). Please increase the sorting window to avoid " - "dropped digis"; + "dropped digis\n"; fSortingWindowWarningIssued = true; } } else { @@ -194,6 +195,12 @@ void GateTimeSorter::Flush() { } Prune(); fFlushed = true; + if (fNumDroppedDigi > 0) { + std::cout << fNumDroppedDigi + << " digis have been dropped while time-sorting. Please increase " + "the sorting time to a value higher than " + << fSortingWindow << " ns\n"; + } } void GateTimeSorter::Prune() { diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h index 841eb6693..6f274578c 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h @@ -64,6 +64,7 @@ class GateTimeSorter { bool fProcessingStarted{false}; bool fFlushed{false}; bool fSortingWindowWarningIssued{false}; + size_t fNumDroppedDigi{}; std::optional fMostRecentTimeArrived; std::optional fMostRecentTimeDeparted; diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index ce2227990..bb5a098bb 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -540,40 +540,40 @@ class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor) """ user_info_defaults = { - "attributes": ( - [], - { - "doc": "Attributes to be considered. ", - }, - ), "input_digi_collection": ( "Singles", { - "doc": "Digi collection to be used as input. ", + "doc": "Digi collection to be used as input.", }, ), "skip_attributes": ( [], { - "doc": "FIXME", + "doc": "Attributes to be omitted from the output.", }, ), "clear_every": ( 1e5, { - "doc": "FIXME", + "doc": "The memory consumed by the actor is minimized after having processed the specified amount of digis", }, ), - "time_window": ( - None, + "pileup_time": ( + 0, { - "doc": "FIXME", + "doc": "Time interval during which consecutive digis are piled up into one digi.", }, ), "group_volume": ( None, { - "doc": "FIXME", + "doc": "Name of the volume in which digis are piled up.", + }, + ), + "sorting_time": ( + 1e3, + { + "doc": "Time interval during which digis are buffered for time-sorting", }, ), } From fa6b6739d0fcc5151385b32c8815c05dde9aabaf Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 09:47:58 +0000 Subject: [PATCH 19/26] initial pile-up policy: non-paralyzing, resulting digi gets time of first digi, sum of energy of all digis, energy-weighted position of all digis, other attributes from highest-energy digi --- .../digitizer/GateDigitizerPileupActor.cpp | 27 ++++++++++++++++++- .../digitizer/GateDigitizerPileupActor.h | 3 +++ opengate/tests/src/test097_pileup_helpers.py | 17 ++++++++++-- 3 files changed, 44 insertions(+), 3 deletions(-) diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index c07e23907..ab98f87d2 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -40,12 +40,16 @@ void GateDigitizerPileupActor::DigitInitialize( const std::vector &attributes_not_in_filler) { auto a = attributes_not_in_filler; + a.push_back("GlobalTime"); a.push_back("TotalEnergyDeposit"); + a.push_back("PostPosition"); GateVDigitizerWithOutputActor::DigitInitialize(a); // Get output attribute pointer + fOutputTimeAttribute = fOutputDigiCollection->GetDigiAttribute("GlobalTime"); fOutputEdepAttribute = fOutputDigiCollection->GetDigiAttribute("TotalEnergyDeposit"); + fOutputPosAttribute = fOutputDigiCollection->GetDigiAttribute("PostPosition"); // Set up pointers to track specific attributes auto &lr = fThreadLocalVDigitizerData.Get(); @@ -103,7 +107,9 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( // Create an iterator to be used when digis will be combined into one digi, // due to pile-up. window.digiIter = window.digis->NewIterator(); + window.digiIter.TrackAttribute("GlobalTime", &l.time); window.digiIter.TrackAttribute("TotalEnergyDeposit", &l.edep); + window.digiIter.TrackAttribute("PostPosition", &l.pos); // Create a filler to copy all digi attributes from the sorted collection // into the collection of the window. window.fillerIn = l.fTimeSorter.CreateFiller(window.digis); @@ -111,7 +117,9 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( // to the output collection (used for the digis that will result from // pile-up). auto filler_out_attributes = window.digis->GetDigiAttributeNames(); + filler_out_attributes.erase("GlobalTime"); filler_out_attributes.erase("TotalEnergyDeposit"); + filler_out_attributes.erase("PostPosition"); window.fillerOut = std::make_unique( window.digis, fOutputDigiCollection, filler_out_attributes); @@ -157,14 +165,23 @@ void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { // into one digi. auto &l = fThreadLocalData.Get(); + std::optional first_time{}; std::optional highest_edep{}; double total_edep = 0.0; size_t highest_edep_index = 0; + G4ThreeVector weighted_position; // Iterate over all digis in the window from the beginning. + window.digiIter.Reset(); while (!window.digiIter.IsAtEnd()) { const auto current_edep = *l.edep; + const auto current_time = *l.time; + const auto current_pos = *l.pos; + // Remember the time of the first digi. + if (!first_time) { + first_time = current_time; + } // Remember the value and index of the highest deposited energy so far. if (!highest_edep.has_value() || current_edep > highest_edep.value()) { highest_edep = current_edep; @@ -172,11 +189,19 @@ void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { } // Accumulate all deposited energy values. total_edep += current_edep; + // Accumulate the energy-weighted position. + weighted_position += current_pos * current_edep; window.digiIter++; } + weighted_position /= total_edep; - // The resulting pile-up digi gets the total edep value as its edep value. + // The resulting pile-up digi gets: + // - the time of the first contributing digi. + fOutputTimeAttribute->FillDValue(*first_time); + // - the total edep value. fOutputEdepAttribute->FillDValue(total_edep); + // - the energy-weighted position. + fOutputPosAttribute->Fill3Value(weighted_position); // All the other attribute values are taken from the digi with the highest // edep. window.fillerOut->Fill(highest_edep_index); diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h index 2fe3a1f05..04a08d274 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -49,7 +49,9 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { int fGroupVolumeDepth; // Output attribute pointer + GateVDigiAttribute *fOutputTimeAttribute{}; GateVDigiAttribute *fOutputEdepAttribute{}; + GateVDigiAttribute *fOutputPosAttribute{}; // Struct for storing digis in one particular volume which belong to the same // time window. @@ -79,6 +81,7 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { GateUniqueVolumeID::Pointer *volID; double *time; double *edep; + G4ThreeVector *pos; GateTimeSorter fTimeSorter; std::map fVolumePileupWindows; diff --git a/opengate/tests/src/test097_pileup_helpers.py b/opengate/tests/src/test097_pileup_helpers.py index 1ecb88f77..2b0462bbf 100644 --- a/opengate/tests/src/test097_pileup_helpers.py +++ b/opengate/tests/src/test097_pileup_helpers.py @@ -9,8 +9,11 @@ def pileup(singles_before_pileup: pd.DataFrame, time_window: float): The singles_before_pileup are in a pandas DataFrame. It returns a dict in which the keys are volume IDs, and the values are a pandas DataFrame representing the singles for that volume ID after pile-up. - Each single after pile-up has a TotalEnergyDeposit equal to the sum of all singles in the same time window. - The other attributes of a single after pile-up are taken from the single in the time window that has the + Each single after pile-up has: + - GlobalTime taken from the first contributing single. + - TotalEnergyDeposit equal to the sum of all singles in the same time window. + - PostPosition equal to the energy-weighted position of all contributing singles. + The other attributes are taken from the single in the time window that has the highest TotalEnergyDeposit. """ df = singles_before_pileup # df for DataFrame @@ -42,6 +45,16 @@ def pileup(singles_before_pileup: pd.DataFrame, time_window: float): pileup_row["TotalEnergyDeposit"] = group_slice[ "TotalEnergyDeposit" ].sum() + pileup_row["GlobalTime"] = group_slice["GlobalTime"].iloc[0] + pileup_row["PostPosition_X"] = ( + group_slice["PostPosition_X"] * group_slice["TotalEnergyDeposit"] + ).sum() / group_slice["TotalEnergyDeposit"].sum() + pileup_row["PostPosition_Y"] = ( + group_slice["PostPosition_Y"] * group_slice["TotalEnergyDeposit"] + ).sum() / group_slice["TotalEnergyDeposit"].sum() + pileup_row["PostPosition_Z"] = ( + group_slice["PostPosition_Z"] * group_slice["TotalEnergyDeposit"] + ).sum() / group_slice["TotalEnergyDeposit"].sum() # Add the combined single to the output. singles_after_pileup.setdefault(volume_id, []).append(pileup_row) # Update current and next for the next time window. From 5ca9424c554a5742a242eb03f27f79894d55aa44 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 09:51:50 +0000 Subject: [PATCH 20/26] update comments --- opengate/tests/src/test097_pileup_helpers.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/opengate/tests/src/test097_pileup_helpers.py b/opengate/tests/src/test097_pileup_helpers.py index 2b0462bbf..b1ac81820 100644 --- a/opengate/tests/src/test097_pileup_helpers.py +++ b/opengate/tests/src/test097_pileup_helpers.py @@ -40,12 +40,15 @@ def pileup(singles_before_pileup: pd.DataFrame, time_window: float): group_slice = group.iloc[current:next_single] max_energy_idx = group_slice["TotalEnergyDeposit"].idxmax() # Create a single with the attribute values from the max energy single, - # except for the TotalEnergyDeposit, take the sum over all singles in the time window. + # except for the TotalEnergyDeposit, GlobalTime and PostPosition. pileup_row = group.loc[max_energy_idx].copy() + # Energy is the sum of energies of all contributing singles. pileup_row["TotalEnergyDeposit"] = group_slice[ "TotalEnergyDeposit" ].sum() + # Time is taken from the first contributing single. pileup_row["GlobalTime"] = group_slice["GlobalTime"].iloc[0] + # PostPosition is the energy-weighted sum of positions from all contributing singles. pileup_row["PostPosition_X"] = ( group_slice["PostPosition_X"] * group_slice["TotalEnergyDeposit"] ).sum() / group_slice["TotalEnergyDeposit"].sum() From 56be3649e4c060ee27b55cf5b559ca631446b5db Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 14:35:10 +0000 Subject: [PATCH 21/26] ensure that setting group_volume works as intended --- opengate/actors/digitizers.py | 19 +++++++------------ opengate/tests/src/test097_pileup.py | 6 +++--- 2 files changed, 10 insertions(+), 15 deletions(-) diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index bb5a098bb..49b237294 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -546,18 +546,6 @@ class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor) "doc": "Digi collection to be used as input.", }, ), - "skip_attributes": ( - [], - { - "doc": "Attributes to be omitted from the output.", - }, - ), - "clear_every": ( - 1e5, - { - "doc": "The memory consumed by the actor is minimized after having processed the specified amount of digis", - }, - ), "pileup_time": ( 0, { @@ -570,6 +558,12 @@ class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor) "doc": "Name of the volume in which digis are piled up.", }, ), + "clear_every": ( + 1e5, + { + "doc": "The memory consumed by the actor is minimized after having processed the specified amount of digis", + }, + ), "sorting_time": ( 1e3, { @@ -601,6 +595,7 @@ def set_group_by_depth(self): def StartSimulationAction(self): DigitizerBase.StartSimulationAction(self) + self.set_group_by_depth() g4.GateDigitizerPileupActor.StartSimulationAction(self) def EndSimulationAction(self): diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py index 14e802bbe..0d376a7a3 100755 --- a/opengate/tests/src/test097_pileup.py +++ b/opengate/tests/src/test097_pileup.py @@ -115,12 +115,12 @@ # Pile-up pu = sim.add_actor("DigitizerPileupActor", "Singles_after_pileup") - pu.attached_to = hc.attached_to - pu.authorize_repeated_volumes = True pu.input_digi_collection = sc.name - pu.output_filename = sc.output_filename + pu.group_volume = crystal.name + pu.authorize_repeated_volumes = True pu.pileup_time = 2000.0 * ns pu.clear_every = 1e4 + pu.output_filename = sc.output_filename # Timing sim.run_timing_intervals = [[0, 0.001 * sec]] From 4a993ac01158d274930305e3ce5f4cd5e5656666 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 14:35:45 +0000 Subject: [PATCH 22/26] add documentation for GateDigitizerPileupActor --- .../user_guide_reference_actors.rst | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/docs/source/user_guide/user_guide_reference_actors.rst b/docs/source/user_guide/user_guide_reference_actors.rst index 5ae9c49a2..b9cdeb84b 100644 --- a/docs/source/user_guide/user_guide_reference_actors.rst +++ b/docs/source/user_guide/user_guide_reference_actors.rst @@ -671,6 +671,51 @@ Reference .. autoclass:: opengate.actors.digitizers.DigitizerEfficiencyActor +DigitizerPileupActor +-------------------- + +Description +~~~~~~~~~~~ + +The combines singles + +Pile‑up occurs when multiple detector interactions happen within a time interval shorter than the resolving/shaping time, +so their pulses overlap and cannot be distinguished as separate events. +The :class:`~.opengate.actors.digitizers.DigitizerPileupActor` simulates this by combining singles in the same volume +(set by `group_volume`) if they occur in a time interval set by `pileup_time`. + +Currently, the following rules apply: + +* The first single occurring after a time period of at least `pileup_time` without singles opens a time window of duration `pileup_time`. +* Any additional single occurring in that time window is merged with the other singles in that window. + These additional singles to not alter the end of the time window (non-paralyzable behavior). +* The resulting combined single has the following attribute values: + + * `GlobalTime`: same value as the first single in the window + * `TotalEnergyDeposit`: sum of the values of all singles in the window + * `PostPosition`: energy-weighted sum of positions of all singles in the window + * remaining attributes values are taken from the single with the highest deposited energy value + +In the future, policies may be added to modify the way in which the attributes of the resulting single are calculated. + +.. code-block:: python + + pu = sim.add_actor("DigitizerPileupActor", "Singles_with_pileup") + pu.input_digi_collection = "Singles" + pu.group_volume = crystal.name + pu.authorize_repeated_volumes = True + pu.pileup_time = 100.0 * ns + pu.output_filename = "singles.root" + +Refer to test097 for an example. + +Reference +~~~~~~~~~ + +.. autoclass:: opengate.actors.digitizers.DigitizerPileupActor + + + Coincidences Sorter ------------------- From e30da48090295bcf7e4459afd0bb77fe2defc5ad Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 16:24:30 +0000 Subject: [PATCH 23/26] fix missing skip_attributes --- opengate/actors/digitizers.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index 49b237294..5a24e24e6 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -570,6 +570,12 @@ class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor) "doc": "Time interval during which digis are buffered for time-sorting", }, ), + "skip_attributes": ( + [], + { + "doc": "Attributes to be omitted from the output.", + }, + ), } def __init__(self, *args, **kwargs): From e1c67740268ff01bfe4e2ccb58fef74bc3e97ea2 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Sun, 28 Dec 2025 16:37:19 +0000 Subject: [PATCH 24/26] add docstring --- opengate/actors/digitizers.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index 5a24e24e6..89e87f6e8 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -536,7 +536,8 @@ def EndSimulationAction(self): class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor): """ - TODO documentation + Dititizer module for simulating pile-up + (combining singles that occur close to each other in time, in the same volume). """ user_info_defaults = { From 1e2806e864d88522f9ffe7026ca69135729a3a09 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Mon, 29 Dec 2025 12:17:28 +0000 Subject: [PATCH 25/26] update GateDigitizerPileupActor documentation --- docs/source/user_guide/user_guide_reference_actors.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/user_guide/user_guide_reference_actors.rst b/docs/source/user_guide/user_guide_reference_actors.rst index b9cdeb84b..73625c014 100644 --- a/docs/source/user_guide/user_guide_reference_actors.rst +++ b/docs/source/user_guide/user_guide_reference_actors.rst @@ -677,8 +677,6 @@ DigitizerPileupActor Description ~~~~~~~~~~~ -The combines singles - Pile‑up occurs when multiple detector interactions happen within a time interval shorter than the resolving/shaping time, so their pulses overlap and cannot be distinguished as separate events. The :class:`~.opengate.actors.digitizers.DigitizerPileupActor` simulates this by combining singles in the same volume @@ -698,6 +696,8 @@ Currently, the following rules apply: In the future, policies may be added to modify the way in which the attributes of the resulting single are calculated. +The :class:`~.opengate.actors.digitizers.DigitizerPileupActor` can currently only be used in single-threaded simulations. + .. code-block:: python pu = sim.add_actor("DigitizerPileupActor", "Singles_with_pileup") From f055adf6a560fae1135f751cbc9bc60dcae56392 Mon Sep 17 00:00:00 2001 From: Gert Van Hoey Date: Tue, 6 Jan 2026 07:42:16 +0000 Subject: [PATCH 26/26] fix typos in documentation --- docs/source/user_guide/user_guide_reference_actors.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/user_guide/user_guide_reference_actors.rst b/docs/source/user_guide/user_guide_reference_actors.rst index 73625c014..7bba46714 100644 --- a/docs/source/user_guide/user_guide_reference_actors.rst +++ b/docs/source/user_guide/user_guide_reference_actors.rst @@ -677,8 +677,8 @@ DigitizerPileupActor Description ~~~~~~~~~~~ -Pile‑up occurs when multiple detector interactions happen within a time interval shorter than the resolving/shaping time, -so their pulses overlap and cannot be distinguished as separate events. +Pile‑up occurs when multiple detector interactions happen within a time interval shorter than the resolving/shaping time. As a result, +their pulses overlap and cannot be distinguished as separate events. The :class:`~.opengate.actors.digitizers.DigitizerPileupActor` simulates this by combining singles in the same volume (set by `group_volume`) if they occur in a time interval set by `pileup_time`. @@ -686,13 +686,13 @@ Currently, the following rules apply: * The first single occurring after a time period of at least `pileup_time` without singles opens a time window of duration `pileup_time`. * Any additional single occurring in that time window is merged with the other singles in that window. - These additional singles to not alter the end of the time window (non-paralyzable behavior). + These additional singles do not alter the end of the time window (non-paralyzable behavior). * The resulting combined single has the following attribute values: * `GlobalTime`: same value as the first single in the window * `TotalEnergyDeposit`: sum of the values of all singles in the window * `PostPosition`: energy-weighted sum of positions of all singles in the window - * remaining attributes values are taken from the single with the highest deposited energy value + * remaining attribute values are taken from the single with the highest deposited energy value In the future, policies may be added to modify the way in which the attributes of the resulting single are calculated.