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..ab98f87d2 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -0,0 +1,210 @@ +/* -------------------------------------------------- + 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" +#include "GateDigiCollectionManager.h" +#include "GateHelpersDigitizer.h" +#include + +GateDigitizerPileupActor::GateDigitizerPileupActor(py::dict &user_info) + : GateVDigitizerWithOutputActor(user_info, false) { + fActions.insert("EndOfEventAction"); + fActions.insert("EndOfRunAction"); +} + +GateDigitizerPileupActor::~GateDigitizerPileupActor() = default; + +void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { + + GateVDigitizerWithOutputActor::InitializeUserInfo(user_info); + if (py::len(user_info) > 0 && user_info.contains("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"); +} + +void GateDigitizerPileupActor::SetGroupVolumeDepth(const int depth) { + fGroupVolumeDepth = depth; +} + +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(); + auto &l = fThreadLocalData.Get(); + + l.fTimeSorter.Init(fInputDigiCollection); + l.fTimeSorter.OutputIterator().TrackAttribute("GlobalTime", &l.time); + l.fTimeSorter.OutputIterator().TrackAttribute("PreStepUniqueVolumeID", + &l.volID); + l.fTimeSorter.SetSortingWindow(fSortingTime); + l.fTimeSorter.SetMaxSize(fClearEveryNEvents); +} + +void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { + auto &l = fThreadLocalData.Get(); + 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); + } + // Make sure everything is output into the root file. + fOutputDigiCollection->FillToRootIfNeeded(true); +} + +GateDigitizerPileupActor::PileupWindow & +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); + + // 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. + 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); + 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( + 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("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); + // 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("GlobalTime"); + filler_out_attributes.erase("TotalEnergyDeposit"); + filler_out_attributes.erase("PostPosition"); + window.fillerOut = std::make_unique( + window.digis, fOutputDigiCollection, filler_out_attributes); + + // Store the PileupWindow in the map and return a reference. + windows[vol_hash] = std::move(window); + return windows[vol_hash]; + } +} + +void GateDigitizerPileupActor::ProcessTimeSortedDigis() { + auto &l = fThreadLocalData.Get(); + 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 + // 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 > 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. + ProcessPileupWindow(window); + window.startTime = current_time; + } + + // Add the current digi to the window. + window.fillerIn->Fill(iter.fIndex); + + iter++; + } + l.fTimeSorter.MarkOutputAsProcessed(); +} + +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(); + + 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; + highest_edep_index = window.digiIter.fIndex; + } + // 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 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); + // Remove all processed digis from the window. + window.digis->Clear(); +} 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..04a08d274 --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -0,0 +1,92 @@ +/* -------------------------------------------------- + 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 "GateTimeSorter.h" +#include "GateVDigitizerWithOutputActor.h" +#include +#include +#include +#include +#include + +namespace py = pybind11; + +/* + * Digitizer module for pile-up. + */ + +class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { + +public: + explicit GateDigitizerPileupActor(py::dict &user_info); + + ~GateDigitizerPileupActor() override; + + void InitializeUserInfo(py::dict &user_info) override; + + // 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 + double fPileupTime; + double fSortingTime; + 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. + 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 & + GetPileupWindowForCurrentVolume(GateUniqueVolumeID::Pointer *volume, + std::map &windows); + + void ProcessTimeSortedDigis(); + void ProcessPileupWindow(PileupWindow &window); + + struct threadLocalT { + GateUniqueVolumeID::Pointer *volID; + double *time; + double *edep; + G4ThreeVector *pos; + + GateTimeSorter fTimeSorter; + std::map fVolumePileupWindows; + }; + G4Cache fThreadLocalData; +}; + +#endif // GateDigitizerPileupActor_h 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..29b22037f --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp @@ -0,0 +1,232 @@ +#include "GateTimeSorter.h" +#include "GateDigiCollection.h" +#include "GateDigiCollectionManager.h" +#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); + + 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(); + + // 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."); + } + fSortingWindow = duration; +} + +void GateTimeSorter::SetMaxSize(size_t maxSize) { + if (fProcessingStarted) { + Fatal("SetMaxSize() cannot be called after Process() has been called."); + } + fMaxSize = 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()."); + } + return std::make_unique( + fOutputCollection, destination, + fOutputCollection->GetDigiAttributeNames()); +} + +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."); + } + fProcessingStarted = true; + + 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)) { + // 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 time (" + << fSortingWindow + << " ns). Please increase the sorting window to avoid " + "dropped digis\n"; + 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; + } + } + iter++; + } + + // 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 { + fOutputCollection->Clear(); + fOutputIter.Reset(); + } +} + +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); + sortedIndices.pop(); + } + 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() { + // 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. + + // 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; + 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 new file mode 100644 index 000000000..6f274578c --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h @@ -0,0 +1,75 @@ +#ifndef GateTimeSorter_h +#define GateTimeSorter_h + +#include "GateDigiCollectionIterator.h" +#include +#include +#include + +class GateDigiCollection; +class GateDigiAttributesFiller; + +class GateTimeSorter { +public: + GateTimeSorter() = default; + + void Init(GateDigiCollection *input); + + 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: + void Prune(); + + struct TimedDigiIndex { + size_t index; + double time; + + bool operator>(const TimedDigiIndex &other) const { + return time > other.time; + } + }; + + struct TimeSortedStorage { + TimeSortedStorage(GateDigiCollection *input, GateDigiCollection *output, + const std::string &name_suffix); + + GateDigiCollection *digis; + std::priority_queue, + std::greater> + sortedIndices; + std::unique_ptr fillerIn; + std::unique_ptr fillerOut; + }; + + double fSortingWindow{1000.0}; // nanoseconds + size_t fMaxSize{100'000}; + + GateDigiCollection *fInputCollection; + GateDigiCollectionIterator fInputIter; + + GateDigiCollection *fOutputCollection; + GateDigiCollectionIterator fOutputIter; + + double *fTime; + + bool fInitialized{false}; + bool fProcessingStarted{false}; + bool fFlushed{false}; + bool fSortingWindowWarningIssued{false}; + size_t fNumDroppedDigi{}; + std::optional fMostRecentTimeArrived; + std::optional fMostRecentTimeDeparted; + + std::unique_ptr fCurrentStorage; + std::unique_ptr fFutureStorage; +}; + +#endif // GateTimeSorter_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..6a8b7d3db --- /dev/null +++ b/core/opengate_core/opengate_lib/digitizer/pyGateDigitizerPileupActor.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 "GateDigitizerPileupActor.h" + +void init_GateDigitizerPileupActor(py::module &m) { + + py::class_, + GateVDigitizerWithOutputActor>(m, "GateDigitizerPileupActor") + .def(py::init()) + .def("SetGroupVolumeDepth", + &GateDigitizerPileupActor::SetGroupVolumeDepth); +} diff --git a/docs/source/user_guide/user_guide_reference_actors.rst b/docs/source/user_guide/user_guide_reference_actors.rst index 5ae9c49a2..7bba46714 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 +~~~~~~~~~~~ + +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`. + +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 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 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. + +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") + 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 ------------------- diff --git a/opengate/actors/digitizers.py b/opengate/actors/digitizers.py index 3882f7191..89e87f6e8 100644 --- a/opengate/actors/digitizers.py +++ b/opengate/actors/digitizers.py @@ -534,6 +534,81 @@ def EndSimulationAction(self): g4.GateDigitizerBlurringActor.EndSimulationAction(self) +class DigitizerPileupActor(DigitizerWithRootOutput, g4.GateDigitizerPileupActor): + """ + Dititizer module for simulating pile-up + (combining singles that occur close to each other in time, in the same volume). + """ + + user_info_defaults = { + "input_digi_collection": ( + "Singles", + { + "doc": "Digi collection to be used as input.", + }, + ), + "pileup_time": ( + 0, + { + "doc": "Time interval during which consecutive digis are piled up into one digi.", + }, + ), + "group_volume": ( + None, + { + "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, + { + "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.GateDigitizerPileupActor.__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.GateDigitizerPileupActor.StartSimulationAction(self) + + def EndSimulationAction(self): + g4.GateDigitizerPileupActor.EndSimulationAction(self) + + class DigitizerSpatialBlurringActor( DigitizerWithRootOutput, g4.GateDigitizerSpatialBlurringActor ): @@ -1276,6 +1351,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, diff --git a/opengate/tests/src/test097_pileup.py b/opengate/tests/src/test097_pileup.py new file mode 100755 index 000000000..0d376a7a3 --- /dev/null +++ b/opengate/tests/src/test097_pileup.py @@ -0,0 +1,137 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +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] + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "gate_test097", "test097") + + 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 + ns = gate.g4_units.ns + 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, 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" + block.color = white + + # Crystal + crystal = sim.add_volume("Box", "crystal") + crystal.mother = block.name + crystal.size = [60 * mm, 20 * mm, 20 * mm] + crystal.material = "LYSO" + crystal.color = green + + 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" + + 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.input_digi_collection = sc.name + 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]] + + sim.run() + print(stats) + + all_match = check_gate_pileup( + sc.output_filename, + "Singles_before_pileup", + "Singles_after_pileup", + 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 new file mode 100644 index 000000000..b1ac81820 --- /dev/null +++ b/opengate/tests/src/test097_pileup_helpers.py @@ -0,0 +1,143 @@ +import pandas as pd +import uproot +import numpy as np + + +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: + - 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 + 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_single = 1 + while next_single < len(times): + # Increment next until it points to the single that opens the next time window. + 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_single] + max_energy_idx = group_slice["TotalEnergyDeposit"].idxmax() + # Create a single with the attribute values from the max energy single, + # 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() + 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. + 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. + singles_after_pileup.setdefault(volume_id, []).append( + group.iloc[current] + ) + # Update current and next for the next time window. + current += 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( + 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 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