diff --git a/core/opengate_core/opengate_lib/digitizer/GateCoincidenceSorterActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateCoincidenceSorterActor.cpp index 299b625e04..4f7e1d60dc 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateCoincidenceSorterActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateCoincidenceSorterActor.cpp @@ -176,7 +176,7 @@ void GateCoincidenceSorterActor::BeginOfRunActionMasterThread(int run_id) { // digis from all threads, so that the coincidence actor can identify // coincidences between singles, irrespective of the thread in which they were // simulated (prompt as well as random coincidences). - fTimeSorter = std::make_unique(); + fTimeSorter = std::make_unique(fOutputDigiCollectionName); fTimeSorter->Init(fInputDigiCollection); fTimeSorter->SetSortingWindow(fSortingTime); fTimeSorter->SetMaxSize(fClearEveryNEvents); diff --git a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp index 426e81f51d..ca11e69302 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.cpp @@ -1,7 +1,7 @@ /* -------------------------------------------------- Copyright (C): OpenGATE Collaboration This software is distributed under the terms - of the GNU Lesser General Public Licence (LGPL) + of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details -------------------------------------------------- */ @@ -12,7 +12,7 @@ #include GateDigitizerPileupActor::GateDigitizerPileupActor(py::dict &user_info) - : GateVDigitizerWithOutputActor(user_info, false) { + : GateVDigitizerWithOutputActor(user_info, true) { fActions.insert("EndOfEventAction"); fActions.insert("EndOfRunAction"); } @@ -69,6 +69,22 @@ void GateDigitizerPileupActor::InitializeUserInfo(py::dict &user_info) { fInputDigiCollectionName = DictGetStr(user_info, "input_digi_collection"); } +void GateDigitizerPileupActor::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("TotalEnergyDeposit", &fTimeSorterOutputEdep); + outputIter.TrackAttribute("PreStepUniqueVolumeID", &fTimeSorterOutputVolID); + + fVolumePileupWindows.clear(); + fWindowExpiry = std::queue(); +} + void GateDigitizerPileupActor::SetGroupVolumeDepth(const int depth) { fGroupVolumeDepth = depth; } @@ -82,41 +98,28 @@ void GateDigitizerPileupActor::DigitInitialize( 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("TotalEnergyDeposit", &l.edep); - l.fTimeSorter.OutputIterator().TrackAttribute("PreStepUniqueVolumeID", - &l.volID); - l.fTimeSorter.SetSortingWindow(fSortingTime); - l.fTimeSorter.SetMaxSize(fClearEveryNEvents); + fOutputDigiCollection->RootInitializeTupleForWorker(); } void GateDigitizerPileupActor::EndOfEventAction(const G4Event *) { - auto &l = fThreadLocalData.Get(); - l.fTimeSorter.Ingest(); - l.fTimeSorter.Process(); - ProcessTimeSortedDigis(); + + fTimeSorter->OnEndOfEventAction([this]() { 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); + + fTimeSorter->OnEndOfRunAction( + [this]() { fOutputDigiCollection->FillToRootIfNeeded(true); }, + [this]() { + ProcessTimeSortedDigis(); + // Process all pile-up windows which still have an expiry item. + while (fWindowExpiry.size() > 0) { + auto &window = + fVolumePileupWindows.at(fWindowExpiry.front().volumeHash); + ProcessPileupWindow(window); + fWindowExpiry.pop(); + } + }); } GateDigitizerPileupActor::PileupWindow & @@ -128,7 +131,7 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( const auto vol_hash = volume->get()->GetIdUpToDepthAsHash(fGroupVolumeDepth); - // Look up the window based on volume hash + // 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. @@ -136,23 +139,25 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( } else { // A PileupWindow object does not yet exist for this volume: create one. PileupWindow window; + window.hash = vol_hash; 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); + window.digis->SetSharedStorage( + true); // Same storage is accessible by all threads. // 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("PostPosition", &l.pos); + window.digiIter.TrackAttribute("TotalEnergyDeposit", &fPileupWindowEdep); + window.digiIter.TrackAttribute("PostPosition", &fPileupWindowPos); // Create a filler to copy all digi attributes from the sorted collection // into the collection of the window. window.fillerIn = std::make_unique( - l.fTimeSorter.OutputCollection(), window.digis, - l.fTimeSorter.OutputCollection()->GetDigiAttributeNames()); + fTimeSorter->OutputCollection(), window.digis, + fTimeSorter->OutputCollection()->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). @@ -169,42 +174,51 @@ GateDigitizerPileupActor::GetPileupWindowForCurrentVolume( } void GateDigitizerPileupActor::ProcessTimeSortedDigis() { - auto &l = fThreadLocalData.Get(); - auto &iter = l.fTimeSorter.OutputIterator(); + + auto &iter = 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); + auto &window = GetPileupWindowForCurrentVolume(fTimeSorterOutputVolID, + fVolumePileupWindows); + + const auto current_time = *fTimeSorterOutputTime; + const auto current_edep = *fTimeSorterOutputEdep; + + // The GlobalTime of digis provided by the time sorter are guaranteed to be + // monotonically increasing. As a consequence, all pile-up windows that have + // expired by the time of the most recently arrived digi can be handled. + ProcessPileupWindows(current_time); - const auto current_time = *l.time; - const auto current_edep = *l.edep; if (window.digis->GetSize() == 0) { - // The window has no digis yet: make the window start at the time of the - // current digi and remember the current edep as highest. + // The window was empty: the newly arrived digi will open it. window.startTime = current_time; + fWindowExpiry.push({window.hash, window.startTime + fTimeWindow}); window.highestEdep = current_edep; - } 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 and remember the current edep as highest. - ProcessPileupWindow(window); - window.startTime = current_time; - window.highestEdep = current_edep; - } - - if (fTimeWindowPolicy == TimeWindowPolicy::Paralyzable) { - // In Paralyzable mode, extend the time window to start at the time of - // the current digi. - window.startTime = current_time; - } else if (fTimeWindowPolicy == TimeWindowPolicy::EnergyWinnerParalyzable) { - // In EnergyWinnerParalyzable mode, extend the time window only if - // the current digi has higher energy than the highest-energy - // digi so far. - if (current_edep > window.highestEdep) { + } else { + // The window was already opened: update the window depending on the + // policy. + switch (fTimeWindowPolicy) { + case TimeWindowPolicy::NonParalyzable: + // window start time remains the same. + break; + case TimeWindowPolicy::Paralyzable: + // The current digi moves the start time forward. window.startTime = current_time; - window.highestEdep = current_edep; + fWindowExpiry.push({window.hash, window.startTime + fTimeWindow}); + break; + case TimeWindowPolicy::EnergyWinnerParalyzable: + // The current digi moves the start time forward if its energy is higher + // than previous energies. + if (current_edep > window.highestEdep) { + window.startTime = current_time; + fWindowExpiry.push({window.hash, window.startTime + fTimeWindow}); + window.highestEdep = current_edep; + } + break; + default: + Fatal("Unknown time window policy"); } } @@ -213,13 +227,11 @@ void GateDigitizerPileupActor::ProcessTimeSortedDigis() { iter++; } - l.fTimeSorter.MarkOutputAsProcessed(); + 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 highest_edep{}; double total_edep = 0.0; @@ -229,22 +241,26 @@ void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { G4ThreeVector weighted_position; G4ThreeVector highest_edep_position; - // Iterate over all digis in the window from the beginning. + if (window.digis->GetSize() == 0) { + return; + } - window.digiIter.Reset(); - while (!window.digiIter.IsAtEnd()) { - const auto current_edep = *l.edep; - const auto current_pos = *l.pos; + // Iterate over all digis in the window from the beginning. + auto &iter = window.digiIter; + iter.GoToBegin(); + while (!iter.IsAtEnd()) { + const auto current_edep = *fPileupWindowEdep; + const auto current_pos = *fPileupWindowPos; // Remember the index of the first digi. if (!first_index) { - first_index = window.digiIter.fIndex; + first_index = iter.fIndex; } // Remember the index of the last digi. - last_index = window.digiIter.fIndex; + last_index = iter.fIndex; // 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; + highest_edep_index = iter.fIndex; highest_edep_position = current_pos; } // Accumulate all deposited energy values. @@ -255,22 +271,28 @@ void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { // Accumulate the energy-weighted position. weighted_position += current_pos * current_edep; } - window.digiIter++; + iter++; } if (fPositionAttributePolicy == PositionAttributePolicy::EnergyWeightedCentroid) { weighted_position /= total_edep; } + // Get output attribute pointer. + auto outputEdepAttribute = + fOutputDigiCollection->GetDigiAttribute("TotalEnergyDeposit"); + auto outputPosAttribute = + fOutputDigiCollection->GetDigiAttribute("PostPosition"); + // The resulting pile-up digi gets: // - the total edep value. - fOutputEdepAttribute->FillDValue(total_edep); + outputEdepAttribute->FillDValue(total_edep); // - the position according to the position attribute policy. if (fPositionAttributePolicy == PositionAttributePolicy::EnergyWinner) { - fOutputPosAttribute->Fill3Value(highest_edep_position); + outputPosAttribute->Fill3Value(highest_edep_position); } else if (fPositionAttributePolicy == PositionAttributePolicy::EnergyWeightedCentroid) { - fOutputPosAttribute->Fill3Value(weighted_position); + outputPosAttribute->Fill3Value(weighted_position); } // All the other attribute values are according to the attribute policy. if (fAttributePolicy == AttributePolicy::First) { @@ -283,3 +305,18 @@ void GateDigitizerPileupActor::ProcessPileupWindow(PileupWindow &window) { // Remove all processed digis from the window. window.digis->Clear(); } + +void GateDigitizerPileupActor::ProcessPileupWindows(double currentTime) { + + // Process the expiry items for which the expiry time is before currentTime. + while (fWindowExpiry.size() > 0 && + currentTime > fWindowExpiry.front().expiryTime) { + auto &window = fVolumePileupWindows.at(fWindowExpiry.front().volumeHash); + // Check again whether the window is actually expired, because its expiry + // time may have been updated in a later expiry item. + if (currentTime > window.startTime + fTimeWindow) { + ProcessPileupWindow(window); + } + fWindowExpiry.pop(); + } +} \ 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 c0fc798fed..e63089d881 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h +++ b/core/opengate_core/opengate_lib/digitizer/GateDigitizerPileupActor.h @@ -1,13 +1,14 @@ /* -------------------------------------------------- Copyright (C): OpenGATE Collaboration This software is distributed under the terms - of the GNU Lesser General Public Licence (LGPL) + of the GNU Lesser General Public Licence (LGPL) See LICENSE.md for further details -------------------------------------------------- */ #ifndef GateDigitizerPileupActor_h #define GateDigitizerPileupActor_h +#include "GateDigiCollection.h" #include "GateTimeSorter.h" #include "GateVDigitizerWithOutputActor.h" #include @@ -15,6 +16,7 @@ #include #include #include +#include namespace py = pybind11; @@ -31,10 +33,10 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { void InitializeUserInfo(py::dict &user_info) override; - // Called every time an Event ends (all threads) + void BeginOfRunActionMasterThread(int run_id) override; + void EndOfEventAction(const G4Event *event) override; - // Called every time a Run ends (all threads) void EndOfRunAction(const G4Run *run) override; void SetGroupVolumeDepth(int depth); @@ -59,14 +61,11 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { 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 { + // Hash of the corresponding volume. + uint64_t hash{}; // Time at which the time window opens. double startTime{}; // Higehst energy deposit in the window. @@ -83,23 +82,32 @@ class GateDigitizerPileupActor : public GateVDigitizerWithOutputActor { std::unique_ptr fillerOut; }; + // Struct that represents when a pile-up window expires. + struct volumeWindowExpiry { + uint64_t volumeHash; + double expiryTime; + }; + + std::unique_ptr fTimeSorter; + std::map fVolumePileupWindows; + std::queue fWindowExpiry; + + // Tracking pointers used by GateTimeSorter output iterator. + GateUniqueVolumeID::Pointer *fTimeSorterOutputVolID{}; + double *fTimeSorterOutputTime{}; + double *fTimeSorterOutputEdep{}; + + // Tracking pointers used by PileupWindow digi iterator. + double *fPileupWindowEdep{}; + G4ThreeVector *fPileupWindowPos{}; + 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; + void ProcessPileupWindows(double currentTime); }; #endif // GateDigitizerPileupActor_h diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp index 7dcfe1559c..6edf83f7b1 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.cpp @@ -6,7 +6,7 @@ #include #include -GateTimeSorter::GateTimeSorter() { +GateTimeSorter::GateTimeSorter(const std::string &name) : fName(name) { fNumWorkingThreads = std::max(1, G4Threading::GetNumberOfRunningWorkerThreads()); fNumActiveWorkingThreads.store(fNumWorkingThreads); @@ -326,6 +326,7 @@ void GateTimeSorter::Process() { while (!iter.IsAtEnd()) { const size_t digiIndex = fSortedCollectionA->GetSize(); const double digiTime = *t; + fNumDigi++; // If a digi is older (lower GlobalTime value) than the newest digi that has // already been transferred to the output collection, then this digi must be // dropped to be able to guarantee time-monotonicity in the output @@ -337,13 +338,14 @@ void GateTimeSorter::Process() { if (fMostRecentTimeDeparted.has_value() && (digiTime < *fMostRecentTimeDeparted)) { ++fNumDroppedDigi; + fMaxDropDelta = + std::max(fMaxDropDelta, *fMostRecentTimeDeparted - digiTime); if (!fSortingWindowWarningIssued) { - std::cout << "The digis in " << fInputCollection->GetName() - << " have non-monotonicities in the GlobalTime attribute " - "that exceed the sorting time (" - << fMinimumSortingWindow - << " ns). Please increase the sorting window to avoid " - "dropped digis\n"; + std::cout << "The digis output by actor '" + << fInputCollection->GetName() + << "' have non-monotonicities in the GlobalTime attribute " + "that exceed the sorting time in actor '" + << fName << "' (" << fMinimumSortingWindow << " ns).\n"; fSortingWindowWarningIssued = true; } } else { @@ -414,10 +416,14 @@ 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 " - << fMinimumSortingWindow << " ns\n"; + const auto percentage = + static_cast(fNumDroppedDigi) / fNumDigi * 100; + std::cout << fNumDroppedDigi << " digis (" << percentage + << " %) have been dropped while time-sorting in actor '" << fName + << "'. Please increase the sorting time to a value higher than " + "the current " + << fMinimumSortingWindow << " ns (suggestion: > " << fMaxDropDelta + << "ns)\n"; } } diff --git a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h index 340ac08f0a..7ffb6c3214 100644 --- a/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h +++ b/core/opengate_core/opengate_lib/digitizer/GateTimeSorter.h @@ -15,7 +15,7 @@ class GateDigiAttributesFiller; class GateTimeSorter { public: - GateTimeSorter(); + GateTimeSorter(const std::string &name); void Init(GateDigiCollection *input); @@ -30,11 +30,11 @@ class GateTimeSorter { GateDigiCollection::Iterator &OutputIterator(); void MarkOutputAsProcessed(); +private: bool Ingest(); void Process(); void Flush(); -private: void IdentifyFastestThread(); void MarkThreadAsFinished(int threadId); void Prune(); @@ -93,11 +93,14 @@ class GateTimeSorter { // GateTimeSorter internal state + std::string fName{}; bool fInitialized{false}; bool fProcessingStarted{false}; bool fFlushed{false}; bool fSortingWindowWarningIssued{false}; size_t fNumDroppedDigi{}; + size_t fNumDigi{}; + double fMaxDropDelta{}; std::optional fMostRecentTimeArrived; std::optional fMostRecentTimeDeparted; }; diff --git a/docs/source/user_guide/user_guide_reference_actors.rst b/docs/source/user_guide/user_guide_reference_actors.rst index cf50175179..b72aad8c18 100644 --- a/docs/source/user_guide/user_guide_reference_actors.rst +++ b/docs/source/user_guide/user_guide_reference_actors.rst @@ -915,8 +915,6 @@ To obtain the same pile-up behavior as in GATE 9, set the following options: position_attribute_policy = "EnergyWinner" attribute_policy = "EnergyWinner" -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") diff --git a/opengate/tests/src/actors/test097_pileup.py b/opengate/tests/src/actors/test097_pileup.py index ac6992971f..3107365a5e 100755 --- a/opengate/tests/src/actors/test097_pileup.py +++ b/opengate/tests/src/actors/test097_pileup.py @@ -9,126 +9,12 @@ PositionAttributePolicy, AttributePolicy, ) - -green = [0, 1, 0, 1] -gray = [0.5, 0.5, 0.5, 1] -white = [1, 1, 1, 0.8] +from test097_pileup_simulation import create_simulation 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.time_window = 2000.0 * ns - pu.clear_every = 1e4 - pu.output_filename = sc.output_filename - - # Timing - sim.run_timing_intervals = [[0, 0.001 * sec]] + sim, pu, root_filename = create_simulation(paths, num_threads=1) test_all_parameter_combinations = False @@ -177,7 +63,7 @@ sim.run(start_new_process=True) all_match = check_gate_pileup( - sc.output_filename, + root_filename, "Singles_before_pileup", "Singles_after_pileup", pu.time_window, diff --git a/opengate/tests/src/actors/test097_pileup_helpers.py b/opengate/tests/src/actors/test097_pileup_helpers.py index 137223ba3f..c39e2e50a1 100644 --- a/opengate/tests/src/actors/test097_pileup_helpers.py +++ b/opengate/tests/src/actors/test097_pileup_helpers.py @@ -183,9 +183,13 @@ def check_gate_pileup( # 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) + actual_singles = ( + actual_singles_after_pileup[ + actual_singles_after_pileup["PreStepUniqueVolumeID"] == volume_id + ] + .sort_values("GlobalTime") + .reset_index(drop=True) + ) # Compare the number of singles if len(expected_singles) != len(actual_singles): @@ -204,11 +208,23 @@ def check_gate_pileup( 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 diff --git a/opengate/tests/src/actors/test097_pileup_mt.py b/opengate/tests/src/actors/test097_pileup_mt.py new file mode 100755 index 0000000000..489e302a90 --- /dev/null +++ b/opengate/tests/src/actors/test097_pileup_mt.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from opengate.tests import utility +from test097_pileup_helpers import ( + check_gate_pileup, + TimeWindowPolicy, + PositionAttributePolicy, + AttributePolicy, +) +from test097_pileup_simulation import create_simulation + +if __name__ == "__main__": + paths = utility.get_default_test_paths(__file__, "gate_test097", "test097") + + sim, pu, root_filename = create_simulation(paths, num_threads=2) + + test_all_parameter_combinations = False + + if test_all_parameter_combinations: + tested_parameter_combinations = [] + for twp in [ + TimeWindowPolicy.NonParalyzable, + TimeWindowPolicy.Paralyzable, + TimeWindowPolicy.EnergyWinnerParalyzable, + ]: + for pap in [ + PositionAttributePolicy.EnergyWeightedCentroid, + PositionAttributePolicy.EnergyWinner, + ]: + for ap in [ + AttributePolicy.First, + AttributePolicy.EnergyWinner, + AttributePolicy.Last, + ]: + tested_parameter_combinations.append((twp, pap, ap)) + else: + tested_parameter_combinations = [ + # Default + ( + TimeWindowPolicy.NonParalyzable, + PositionAttributePolicy.EnergyWeightedCentroid, + AttributePolicy.First, + ), + # GATE 9 behavior + ( + TimeWindowPolicy.EnergyWinnerParalyzable, + PositionAttributePolicy.EnergyWinner, + AttributePolicy.EnergyWinner, + ), + ] + + all_tests_ok = True + + for c in tested_parameter_combinations: + twp, pap, ap = c + + pu.time_window_policy = twp.name + pu.position_attribute_policy = pap.name + pu.attribute_policy = ap.name + + sim.run(start_new_process=True) + + all_match = check_gate_pileup( + root_filename, + "Singles_before_pileup", + "Singles_after_pileup", + pu.time_window, + twp, + pap, + ap, + ) + + if not all_match: + print(f"Pileup test failed for {twp}, {pap}, {ap}") + all_tests_ok = False + + utility.test_ok(all_tests_ok) diff --git a/opengate/tests/src/actors/test097_pileup_simulation.py b/opengate/tests/src/actors/test097_pileup_simulation.py new file mode 100644 index 0000000000..8a332c81ee --- /dev/null +++ b/opengate/tests/src/actors/test097_pileup_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_pileup") + 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 + + # 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.time_window = 2000.0 * ns + pu.clear_every = 1e4 + pu.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, pu, root_filename)