diff --git a/core/opengate_core/opengate_core.cpp b/core/opengate_core/opengate_core.cpp index f11b57132..22017c541 100644 --- a/core/opengate_core/opengate_core.cpp +++ b/core/opengate_core/opengate_core.cpp @@ -451,6 +451,8 @@ void init_GateHitsCollectionActor(py::module &); void init_GateHitsAdderActor(py::module &); +void init_GateDigitizerDeadTimeActor(py::module &); + void init_GateDigitizerPileupActor(py::module &); void init_GateDigitizerReadoutActor(py::module &m); @@ -765,6 +767,7 @@ PYBIND11_MODULE(opengate_core, m) { init_GateDigiAttributeManager(m); init_GateVDigiAttribute(m); init_GateHitsAdderActor(m); + init_GateDigitizerDeadTimeActor(m); init_GateDigitizerPileupActor(m); init_GateDigitizerReadoutActor(m); init_GateDigitizerBlurringActor(m); diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerDeadTimeActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerDeadTimeActor.cpp new file mode 100644 index 000000000..41530e208 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerDeadTimeActor.cpp @@ -0,0 +1,123 @@ +/* -------------------------------------------------- + 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 "GateDigitizerDeadTimeActor.h" +#include "../GateHelpersDict.h" +#include "GateHelpersDigitizer.h" +#include + +GateDigitizerDeadTimeActor::GateDigitizerDeadTimeActor(py::dict &user_info) + : GateVDigitizerWithOutputActor(user_info, true) { + fActions.insert("EndOfEventAction"); + fActions.insert("EndOfRunAction"); +} + +GateDigitizerDeadTimeActor::~GateDigitizerDeadTimeActor() = default; + +void GateDigitizerDeadTimeActor::InitializeUserInfo(py::dict &user_info) { + + GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); + if (py::len(user_info) > 0 && user_info.contains("dead_time")) { + fDeadTime = DictGetDouble(user_info, "dead_time"); // nanoseconds + } + if (py::len(user_info) > 0 && user_info.contains("policy")) { + const auto policy_str = DictGetStr(user_info, "policy"); + if (policy_str == "NonParalyzable") { + fPolicy = DeadTimePolicy::NonParalyzable; + } else if (policy_str == "Paralyzable") { + fPolicy = DeadTimePolicy::Paralyzable; + } else { + Fatal("Unknown dead time policy '" + policy_str + "'"); + } + } + + 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"); +} + +void GateDigitizerDeadTimeActor::BeginOfRunActionMasterThread(int run_id) { + + fTimeSorter = std::make_unique(fOutputDigiCollectionName); + fTimeSorter->Init(fInputDigiCollection); + fTimeSorter->SetSortingWindow(fSortingTime); + fTimeSorter->SetMaxSize(fClearEveryNEvents); + + auto &outputIter = fTimeSorter->OutputIterator(); + outputIter.TrackAttribute("GlobalTime", &fTimeSorterOutputTime); + outputIter.TrackAttribute("PreStepUniqueVolumeID", &fTimeSorterOutputVolID); + + fillerOut = std::make_unique( + fTimeSorter->OutputCollection(), fOutputDigiCollection, + fTimeSorter->OutputCollection()->GetDigiAttributeNames()); + + fVolumeEndOfDeadTimeInterval.clear(); +} + +void GateDigitizerDeadTimeActor::SetGroupVolumeDepth(const int depth) { + fGroupVolumeDepth = depth; +} + +void GateDigitizerDeadTimeActor::DigitInitialize( + const std::vector &attributes_not_in_filler) { + + auto a = attributes_not_in_filler; + GateVDigitizerWithOutputActor::DigitInitialize(a); + + fOutputDigiCollection->RootInitializeTupleForWorker(); +} + +void GateDigitizerDeadTimeActor::EndOfEventAction(const G4Event *) { + + fTimeSorter->OnEndOfEventAction([this]() { ProcessTimeSortedDigis(); }); +} + +void GateDigitizerDeadTimeActor::EndOfRunAction(const G4Run *) { + + fTimeSorter->OnEndOfRunAction( + [this]() { fOutputDigiCollection->FillToRootIfNeeded(true); }, + [this]() { ProcessTimeSortedDigis(); }); +} + +void GateDigitizerDeadTimeActor::ProcessTimeSortedDigis() { + + auto &iter = fTimeSorter->OutputIterator(); + iter.GoToBegin(); + while (!iter.IsAtEnd()) { + + const auto volHash = + fTimeSorterOutputVolID->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth); + std::optional endTime{}; + auto it = fVolumeEndOfDeadTimeInterval.find(volHash); + if (it != fVolumeEndOfDeadTimeInterval.end()) { + endTime = it->second; + } + + const auto currentTime = *fTimeSorterOutputTime; + if (!endTime.has_value() || currentTime > *endTime) { + // Digi goes to the output. + fillerOut->Fill(iter.fIndex); + // Update end of dead time interval. + fVolumeEndOfDeadTimeInterval[volHash] = currentTime + fDeadTime; + } else { + // Digi is dropped because it arrived during a dead time interval. + if (fPolicy == DeadTimePolicy::NonParalyzable) { + // End of dead time interval does not change. + } else if (fPolicy == DeadTimePolicy::Paralyzable) { + // End of dead time interval is updated. + fVolumeEndOfDeadTimeInterval[volHash] = currentTime + fDeadTime; + } else { + Fatal("Unknown dead time policy"); + } + } + + iter++; + } + fTimeSorter->MarkOutputAsProcessed(); +} diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerDeadTimeActor.h b/core/opengate_core/opengate_lib/digitizer/GateDigitizerDeadTimeActor.h new file mode 100644 index 000000000..b1ce464ec --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerDeadTimeActor.h @@ -0,0 +1,71 @@ +/* -------------------------------------------------- + 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 GateDigitizerDeadTimeActor_h +#define GateDigitizerDeadTimeActor_h + +#include "GateDigiCollection.h" +#include "GateTimeSorter.h" +#include "GateVDigitizerWithOutputActor.h" +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; + +/* + * Digitizer module for dead time. + */ + +class GateDigitizerDeadTimeActor : public GateVDigitizerWithOutputActor { + +public: + explicit GateDigitizerDeadTimeActor(py::dict &user_info); + + ~GateDigitizerDeadTimeActor() override; + + void InitializeUserInfo(py::dict &user_info) override; + + void BeginOfRunActionMasterThread(int run_id) override; + + void EndOfEventAction(const G4Event *event) override; + + void EndOfRunAction(const G4Run *run) override; + + void SetGroupVolumeDepth(int depth); + +protected: + enum class DeadTimePolicy { + NonParalyzable, + Paralyzable, + }; + + void DigitInitialize( + const std::vector &attributes_not_in_filler) override; + + // User parameters + double fDeadTime; + DeadTimePolicy fPolicy; + double fSortingTime; + int fGroupVolumeDepth; + + std::unique_ptr fTimeSorter; + std::map fVolumeEndOfDeadTimeInterval; + + // Tracking pointers used by GateTimeSorter output iterator. + GateUniqueVolumeID::Pointer *fTimeSorterOutputVolID{}; + double *fTimeSorterOutputTime{}; + + std::unique_ptr fillerOut; + + void ProcessTimeSortedDigis(); +}; + +#endif // GateDigitizerDeadTimeActor_h diff --git a/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerDeadTimeActor.cpp b/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerDeadTimeActor.cpp new file mode 100644 index 000000000..be9cca630 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerDeadTimeActor.cpp @@ -0,0 +1,23 @@ +/* -------------------------------------------------- + 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 "GateDigitizerDeadTimeActor.h" + +void init_GateDigitizerDeadTimeActor(py::module &m) { + + py::class_, + GateVDigitizerWithOutputActor>(m, "GateDigitizerDeadTimeActor") + .def(py::init()) + .def("SetGroupVolumeDepth", + &GateDigitizerDeadTimeActor::SetGroupVolumeDepth); +} diff --git a/docs/source/user_guide/user_guide_reference_actors.rst b/docs/source/user_guide/user_guide_reference_actors.rst index b72aad8c1..04062d32b 100644 --- a/docs/source/user_guide/user_guide_reference_actors.rst +++ b/docs/source/user_guide/user_guide_reference_actors.rst @@ -878,6 +878,49 @@ Reference .. autoclass:: opengate.actors.digitizers.DigitizerEfficiencyActor +DigitizerDeadTimeActor +---------------------- + +Description +~~~~~~~~~~~ + +Dead time is the period following the detection of a single during which a detector cannot register a new event. +The :class:`~.opengate.actors.digitizers.DigitizerDeadTimeActor` simulates this effect by dropping singles that arrive within the dead time interval +after a previously accepted single, on a per-volume basis (set by `group_volume`). + +Two policies are available, controlled by the `policy` parameter: + +* **NonParalyzable** (default): the dead time interval is fixed in duration, starting from each accepted single. +Singles arriving during that interval are discarded, but do not extend it. This models, for example, +a detector that keeps accepting events after the dead time has elapsed regardless of what happened during the interval. +* **Paralyzable**: every single that arrives during the current dead time interval extends it by `dead_time` from the time of that single. +If the rate is very high, the detector can become fully paralyzed because each new single keeps resetting the interval. + +The actor internally time-sorts its input digis using a buffer window controlled by `sorting_time` before applying the dead time logic. +This is necessary to correctly handle events from different threads or runs that may arrive out of order. + +.. code-block:: python + + ns = gate.g4_units.ns + + dt = sim.add_actor("DigitizerDeadTimeActor", "Singles_after_deadtime") + dt.attached_to = crystal.name + dt.authorize_repeated_volumes = True + dt.input_digi_collection = "Singles" + dt.group_volume = crystal.name + dt.dead_time = 1.0 * ns + dt.policy = "NonParalyzable" # or "Paralyzable" + dt.output_filename = "singles.root" + +Refer to `test100 `_ for an example. +The multi-threaded variant is in `test100_deadtime_actor_mt `_. + +Reference +~~~~~~~~~ + +.. autoclass:: opengate.actors.digitizers.DigitizerDeadTimeActor + + DigitizerPileupActor -------------------- diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index c551eea9d..f1750075e 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -538,6 +538,93 @@ def EndSimulationAction(self): g4.GateDigitizerBlurringActor.EndSimulationAction(self) +class DigitizerDeadTimeActor(DigitizerWithRootOutput, g4.GateDigitizerDeadTimeActor): + """ + Dititizer module for simulating dead time + (drops singles that occur briefly after a previous one, in the same volume). + """ + + user_info_defaults = { + "input_digi_collection": ( + "Singles", + { + "doc": "Digi collection to be used as input.", + }, + ), + "dead_time": ( + 0, + { + "doc": "Time interval after one digi during which following digis are dropped", + }, + ), + "policy": ( + "NonParalyzable", + { + "doc": "Policy controlling whether the dead time interval is extended when new digis occur. " + + "NonParalyzable: the dead time interval is fixed and does not change when new digis occur during that interval. " + + "Paralyzable: when a new digi arrives during the current dead time interval, the interval is extended to last until dead_time after the new digi. ", + "allowed_values": ( + "NonParalyzable", + "Paralyzable", + ), + }, + ), + "group_volume": ( + None, + { + "doc": "Name of the volume in which the dead time effect takes place", + }, + ), + "clear_every": ( + 1e5, + { + "doc": "The memory consumed by the actor is minimized after having processed the specified amount of digis", + }, + ), + "sorting_time": ( + 1e3, + { + "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): + DigitizerBase.__init__(self, *args, **kwargs) + self.__initcpp__() + + def __initcpp__(self): + g4.GateDigitizerDeadTimeActor.__init__(self, self.user_info) + self.AddActions({"StartSimulationAction", "EndSimulationAction"}) + + def initialize(self): + DigitizerBase.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) + self.set_group_by_depth() + g4.GateDigitizerDeadTimeActor.StartSimulationAction(self) + + def EndSimulationAction(self): + g4.GateDigitizerDeadTimeActor.EndSimulationAction(self) + + class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor): """ Dititizer module for simulating pile-up @@ -1581,6 +1668,7 @@ def volume_name(self, value): process_cls(DigitizerWithRootOutput) process_cls(DigitizerAdderActor) process_cls(DigitizerBlurringActor) +process_cls(DigitizerDeadTimeActor) process_cls(DigitizerPileupActor) process_cls(DigitizerSpatialBlurringActor) process_cls(DigitizerEfficiencyActor) diff --git a/opengate/managers.py b/opengate/managers.py index 0dae8d5c3..d5242eded 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -76,6 +76,7 @@ DigitizerEfficiencyActor, DigitizerEnergyWindowsActor, DigitizerHitsCollectionActor, + DigitizerDeadTimeActor, DigitizerPileupActor, DigitizerProjectionActor, DigitizerReadoutActor, @@ -171,6 +172,7 @@ "DigitizerProjectionActor": DigitizerProjectionActor, "DigitizerEnergyWindowsActor": DigitizerEnergyWindowsActor, "DigitizerHitsCollectionActor": DigitizerHitsCollectionActor, + "DigitizerDeadTimeActor": DigitizerDeadTimeActor, "DigitizerPileupActor": DigitizerPileupActor, "CoincidenceSorterActor": CoincidenceSorterActor, "DigiAttributeProcessDefinedStepInVolumeActor": DigiAttributeProcessDefinedStepInVolumeActor, diff --git a/opengate/tests/src/actors/test100_deadtime_actor.py b/opengate/tests/src/actors/test100_deadtime_actor.py new file mode 100644 index 000000000..83b269f9f --- /dev/null +++ b/opengate/tests/src/actors/test100_deadtime_actor.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from opengate.tests import utility +from test100_deadtime_helpers import ( + check_gate_deadtime, + DeadTimePolicy, +) +from test100_deadtime_simulation import create_simulation + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "gate_test100", "test100") + + sim, dt, root_filename = create_simulation(paths, num_threads=1) + + test_all_parameter_combinations = False + + all_tests_ok = True + + for policy in [ + DeadTimePolicy.NonParalyzable, + DeadTimePolicy.Paralyzable, + ]: + dt.policy = policy.name + + sim.run(start_new_process=True) + + all_match = check_gate_deadtime( + root_filename, + "Singles_before_deadtime", + "Singles_after_deadtime", + dt.dead_time, + policy, + ) + + if not all_match: + print(f"Dead time test failed for {policy}") + all_tests_ok = False + + utility.test_ok(all_tests_ok) diff --git a/opengate/tests/src/actors/test100_deadtime_actor_mt.py b/opengate/tests/src/actors/test100_deadtime_actor_mt.py new file mode 100644 index 000000000..40b5c02f4 --- /dev/null +++ b/opengate/tests/src/actors/test100_deadtime_actor_mt.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from opengate.tests import utility +from test100_deadtime_helpers import ( + check_gate_deadtime, + DeadTimePolicy, +) +from test100_deadtime_simulation import create_simulation + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "gate_test100", "test100") + + sim, dt, root_filename = create_simulation(paths, num_threads=2) + + test_all_parameter_combinations = False + + all_tests_ok = True + + for policy in [ + DeadTimePolicy.NonParalyzable, + DeadTimePolicy.Paralyzable, + ]: + dt.policy = policy.name + + sim.run(start_new_process=True) + + all_match = check_gate_deadtime( + root_filename, + "Singles_before_deadtime", + "Singles_after_deadtime", + dt.dead_time, + policy, + ) + + if not all_match: + print(f"Dead time test failed for {policy}") + all_tests_ok = False + + utility.test_ok(all_tests_ok) diff --git a/opengate/tests/src/actors/test100_deadtime_helpers.py b/opengate/tests/src/actors/test100_deadtime_helpers.py new file mode 100644 index 000000000..cde581eea --- /dev/null +++ b/opengate/tests/src/actors/test100_deadtime_helpers.py @@ -0,0 +1,133 @@ +import pandas as pd +import uproot +import numpy as np + +from enum import Enum, auto + + +class DeadTimePolicy(Enum): + NonParalyzable = auto() + Paralyzable = auto() + + +def deadtime( + singles_before_deadtime: pd.DataFrame, + dead_time: float, + policy: DeadTimePolicy = DeadTimePolicy.NonParalyzable, +) -> dict: + """ + This function simulates dead time (given in ns). + The singles_before_deadtime 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 dead time. + """ + df = singles_before_deadtime # df for DataFrame + grouped = df.groupby("PreStepUniqueVolumeID") + + singles_after_deadtime = {} + 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 + + if len(times) == 0: + continue + + next_single = 0 + while next_single < len(times): + current = next_single # index of a single that opens the current dead time interval + next_single = current + 1 + begin = times[current] + singles_after_deadtime.setdefault(volume_id, []).append(group.iloc[current]) + # Increment next until it points to the single that opens the next dead time interval. + while next_single < len(times) and times[next_single] <= begin + dead_time: + if policy == DeadTimePolicy.Paralyzable: + begin = times[next_single] + next_single += 1 + + return singles_after_deadtime + + +def check_gate_deadtime( + root_file_path: str, + name_before_deadtime: str, + name_after_deadtime: str, + dead_time: float, + policy: DeadTimePolicy, +): + """ + This function checks that the singles generated by GateDigitizerDeadTimeActor are identical + to the ones generated by the Python dead time implementation in the deadtime() function above. + """ + with uproot.open(root_file_path) as root_file: + singles_before_deadtime = root_file[name_before_deadtime].arrays(library="pd") + actual_singles_after_deadtime = root_file[name_after_deadtime].arrays( + library="pd" + ) + + expected_singles_after_deadtime = deadtime( + singles_before_deadtime, + dead_time, + policy, + ) + num_expected_singles = sum( + len(entries) for entries in expected_singles_after_deadtime.values() + ) + + print(f"Singles before dead time: {len(singles_before_deadtime)}") + print( + f"Expected singles after dead time: {num_expected_singles} ({num_expected_singles / len(singles_before_deadtime) * 100:.01f}%)" + ) + print(f"Actual singles after dead time: {len(actual_singles_after_deadtime)}") + + all_match = True + for volume_id in expected_singles_after_deadtime.keys(): + + # Get the expected singles (Python) and actual singles (GateDigitizerDeadTimeActor) for the current volume. + expected_singles = pd.DataFrame(expected_singles_after_deadtime[volume_id]) + actual_singles = ( + actual_singles_after_deadtime[ + actual_singles_after_deadtime["PreStepUniqueVolumeID"] == volume_id + ] + .sort_values("GlobalTime") + .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 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") + mismatch_indices = np.where( + ~np.isclose(expected_values, actual_values, rtol=1e-9) + )[0] + first_mismatch_idx = mismatch_indices[0] + print( + f" Single index {first_mismatch_idx}: expected={expected_values[first_mismatch_idx]}, actual={actual_values[first_mismatch_idx]}" + ) + all_match = False + break + else: + if not all(expected_values == actual_values): + print(f"Volume {volume_id}: Attribute {attr} does not match") + mismatch_indices = np.where(expected_values != actual_values)[0] + first_mismatch_idx = mismatch_indices[0] + print( + f" Single index {first_mismatch_idx}: expected={expected_values[first_mismatch_idx]}, actual={actual_values[first_mismatch_idx]}" + ) + all_match = False + break + + return all_match diff --git a/opengate/tests/src/actors/test100_deadtime_simulation.py b/opengate/tests/src/actors/test100_deadtime_simulation.py new file mode 100644 index 000000000..6b9fa504a --- /dev/null +++ b/opengate/tests/src/actors/test100_deadtime_simulation.py @@ -0,0 +1,116 @@ +import opengate as gate + + +def create_simulation(paths, num_threads): + + root_filename = paths.output / "output_singles.root" + + sim = gate.Simulation() + + sim.visu = False + sim.random_seed = 1234 + sim.number_of_threads = num_threads + + # Units + mm = gate.g4_units.mm + sec = gate.g4_units.s + ns = gate.g4_units.ns + Bq = gate.g4_units.Bq + gcm3 = gate.g4_units.g_cm3 + deg = gate.g4_units.deg + + # 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.material = "G4_AIR" + + # Block + block = sim.add_volume("Box", "block") + block.mother = pet.name + block.size = [60 * mm, 20 * mm, 20 * mm] + translations_ring, rotations_ring = gate.geometry.utility.get_circular_repetition( + 40, [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" + + # Crystal + crystal = sim.add_volume("Box", "crystal") + crystal.mother = block.name + crystal.size = [60 * mm, 20 * mm, 20 * mm] + crystal.material = "LYSO" + + source1 = sim.add_source("GenericSource", "b2b_1") + source1.particle = "back_to_back" + source1.activity = 5 * 1e6 * Bq / num_threads + 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 / num_threads + 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" + + 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_deadtime") + sc.attached_to = hc.attached_to + sc.authorize_repeated_volumes = True + sc.input_digi_collection = hc.name + sc.policy = "EnergyWinnerPosition" + sc.output_filename = root_filename + + # Dead time + dt = sim.add_actor("DigitizerDeadTimeActor", "Singles_after_deadtime") + dt.input_digi_collection = sc.name + dt.group_volume = crystal.name + dt.authorize_repeated_volumes = True + dt.dead_time = 2000.0 * ns + dt.clear_every = 1e4 + dt.output_filename = root_filename + + # Timing + run_duration = 0.0005 * sec + num_runs = 2 + sim.run_timing_intervals = [ + [2 * r * run_duration, (2 * r + 1) * run_duration] for r in range(num_runs) + ] + + return (sim, dt, root_filename)