diff --git a/core/opengate_core/opengate_lib/GateDebugSource.cpp b/core/opengate_core/opengate_lib/GateDebugSource.cpp index 668b42c6b..483d1326f 100644 --- a/core/opengate_core/opengate_lib/GateDebugSource.cpp +++ b/core/opengate_core/opengate_lib/GateDebugSource.cpp @@ -36,7 +36,7 @@ void GateDebugSource::UpdateActivity(const double time) { } double GateDebugSource::PrepareNextTime(const double current_simulation_time, - double NumberOfGeneratedEvents) { + unsigned long NumberOfGeneratedEvents) { DDD("GateDebugSource::PrepareNextTime", " ", G4BestUnit(current_simulation_time, "Time"), " ", NumberOfGeneratedEvents); diff --git a/core/opengate_core/opengate_lib/GateDebugSource.h b/core/opengate_core/opengate_lib/GateDebugSource.h index 2f46ffd6d..d02529560 100644 --- a/core/opengate_core/opengate_lib/GateDebugSource.h +++ b/core/opengate_core/opengate_lib/GateDebugSource.h @@ -27,7 +27,7 @@ class GateDebugSource : public GateVSource { void UpdateActivity(const double time) override; double PrepareNextTime(double current_simulation_time, - double NumberOfGeneratedEvents) override; + unsigned long NumberOfGeneratedEvents) override; void PrepareNextRun() override; diff --git a/core/opengate_core/opengate_lib/GateGenericSource.cpp b/core/opengate_core/opengate_lib/GateGenericSource.cpp index b2785b682..1cf076c6e 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.cpp +++ b/core/opengate_core/opengate_lib/GateGenericSource.cpp @@ -139,8 +139,9 @@ void GateGenericSource::UpdateActivityWithTAC(const double time) { fActivity = fTAC_Activities[i] * w1 + fTAC_Activities[i + 1] * w2; } -double GateGenericSource::PrepareNextTime(const double current_simulation_time, - double NumberOfGeneratedEvents) { +double +GateGenericSource::PrepareNextTime(const double current_simulation_time, + unsigned long NumberOfGeneratedEvents) { // initialization of the effective event time (it can be in the // future according to the current_simulation_time) if (fEffectiveEventTime < current_simulation_time) { diff --git a/core/opengate_core/opengate_lib/GateGenericSource.h b/core/opengate_core/opengate_lib/GateGenericSource.h index 2984d1d49..bc506d004 100644 --- a/core/opengate_core/opengate_lib/GateGenericSource.h +++ b/core/opengate_core/opengate_lib/GateGenericSource.h @@ -32,7 +32,7 @@ class GateGenericSource : public GateVSource { void InitializeUserInfo(py::dict &user_info) override; double PrepareNextTime(double current_simulation_time, - double NumberOfGeneratedEvents) override; + unsigned long NumberOfGeneratedEvents) override; void PrepareNextRun() override; diff --git a/core/opengate_core/opengate_lib/GateLastVertexSource.cpp b/core/opengate_core/opengate_lib/GateLastVertexSource.cpp index ab61b05c6..fc81b1ff6 100644 --- a/core/opengate_core/opengate_lib/GateLastVertexSource.cpp +++ b/core/opengate_core/opengate_lib/GateLastVertexSource.cpp @@ -19,8 +19,9 @@ void GateLastVertexSource::InitializeUserInfo(py::dict &user_info) { fN = 0; } -double GateLastVertexSource::PrepareNextTime(double current_simulation_time, - double NumberOfGeneratedEvents) { +double +GateLastVertexSource::PrepareNextTime(double current_simulation_time, + unsigned long NumberOfGeneratedEvents) { /* // If all N events have been generated, we stop (negative time) diff --git a/core/opengate_core/opengate_lib/GateLastVertexSource.h b/core/opengate_core/opengate_lib/GateLastVertexSource.h index f7f5aecff..ffa90d0de 100644 --- a/core/opengate_core/opengate_lib/GateLastVertexSource.h +++ b/core/opengate_core/opengate_lib/GateLastVertexSource.h @@ -32,7 +32,7 @@ class GateLastVertexSource : public GateVSource { void InitializeUserInfo(py::dict &user_info) override; double PrepareNextTime(double current_simulation_time, - double NumberOfGeneratedEvents) override; + unsigned long NumberOfGeneratedEvents) override; void PrepareNextRun() override; diff --git a/core/opengate_core/opengate_lib/GateSourceManager.cpp b/core/opengate_core/opengate_lib/GateSourceManager.cpp index 38d7bd59b..8b407ba22 100644 --- a/core/opengate_core/opengate_lib/GateSourceManager.cpp +++ b/core/opengate_core/opengate_lib/GateSourceManager.cpp @@ -17,6 +17,7 @@ #if defined(_MSC_VER) #pragma warning(pop) #endif +#include #include #include #include @@ -25,6 +26,7 @@ #include #include #include +#include #include #ifdef USE_GDML @@ -34,6 +36,9 @@ /* There will be one SourceManager per thread */ bool GateSourceManager::fRunTerminationFlag = false; +std::atomic GateSourceManager::fGeneratedPrimariesThisRun{0}; +std::atomic GateSourceManager::fPrimaryLimitWarningIssued{false}; +std::uint64_t GateSourceManager::fMaxPrimariesPerRun = INT32_MAX; GateSourceManager::GateSourceManager() { fUIEx = nullptr; @@ -66,6 +71,50 @@ void GateSourceManager::SetRunTerminationFlag(bool flag) { fRunTerminationFlag = flag; } +std::uint64_t GateSourceManager::GetPlatformMaxPrimariesPerRun() { + return std::numeric_limits::max(); +} + +void GateSourceManager::SetMaxPrimariesPerRun(std::uint64_t value) { + if (value == 0) { + fMaxPrimariesPerRun = GetPlatformMaxPrimariesPerRun(); + return; + } + fMaxPrimariesPerRun = value; +} + +void GateSourceManager::ResetPrimaryCounterForRun() { + fGeneratedPrimariesThisRun.store(0, std::memory_order_relaxed); + fPrimaryLimitWarningIssued.store(false, std::memory_order_relaxed); +} + +bool GateSourceManager::TryReservePrimarySlot() { + auto current = fGeneratedPrimariesThisRun.load(std::memory_order_relaxed); + while (current < fMaxPrimariesPerRun) { + if (fGeneratedPrimariesThisRun.compare_exchange_weak( + current, current + 1, std::memory_order_relaxed, + std::memory_order_relaxed)) { + return true; + } + } + return false; +} + +void GateSourceManager::WarnPrimaryLimitReached() { + bool warning_already_issued = + fPrimaryLimitWarningIssued.exchange(true, std::memory_order_relaxed); + if (warning_already_issued) { + return; + } + + G4ExceptionDescription ed; + ed << "The current Geant4 BeamOn() call reached the maximum supported " + << "number of primaries (" << fMaxPrimariesPerRun << "). " + << "The run is being stopped before generating additional primaries."; + G4Exception("GateSourceManager::GeneratePrimaries", "GateSourceManager0001", + JustWarning, ed); +} + void GateSourceManager::Initialize(const TimeIntervals &simulation_times, py::dict &options) { fSimulationTimes = simulation_times; @@ -86,6 +135,11 @@ void GateSourceManager::Initialize(const TimeIntervals &simulation_times, fVisCommands = DictGetVecStr(options, "visu_commands"); fVerboseLevel = DictGetInt(options, "running_verbose_level"); fProgressBarFlag = DictGetBool(options, "progress_bar"); + if (options.contains("max_primaries_per_run")) { + SetMaxPrimariesPerRun(DictGetInt(options, "max_primaries_per_run")); + } else { + SetMaxPrimariesPerRun(GetPlatformMaxPrimariesPerRun()); + } InstallSignalHandler(); InitializeProgressBar(); @@ -158,6 +212,7 @@ void GateSourceManager::StartMasterThread() { auto &l = fThreadLocalData.Get(); l.fStartNewRun = true; for (size_t run_id = 0; run_id < fSimulationTimes.size(); run_id++) { + ResetPrimaryCounterForRun(); // Start Begin Of Run for MasterThread // (both for multi-thread and mono-thread app) // The conventional (threaded) BeginOfRun will be called @@ -264,7 +319,8 @@ void GateSourceManager::PrepareNextSource() const { // Ask all sources their next time, keep the closest one for (auto *source : fSources) { - G4int numberOfSimulatedEvents = source->GetNumberOfSimulatedEvents(); + unsigned long numberOfSimulatedEvents = + source->GetNumberOfSimulatedEvents(); auto t = source->PrepareNextTime(l.fCurrentSimulationTime, numberOfSimulatedEvents); if ((t >= min_time) && (t < max_time)) { @@ -319,20 +375,35 @@ void GateSourceManager::GeneratePrimaries(G4Event *event) { event->AddPrimaryVertex(vertex); // (allocated memory is leaked) } else { - // shoot particle - l.fNextActiveSource->GeneratePrimaries(event, l.fCurrentSimulationTime); - // log (after particle creation) - if (LogLevel_EVENT <= GateSourceManager::fVerboseLevel) { - const auto *prim = event->GetPrimaryVertex(0)->GetPrimary(0); - std::string t = G4BestUnit(l.fCurrentSimulationTime, "Time"); - std::string e = G4BestUnit(prim->GetKineticEnergy(), "Energy"); - std::string s = l.fNextActiveSource->fName; - Log(LogLevel_EVENT, fVerboseLevel, - "Event {} {} {} {} {:.2f} {:.2f} {:.2f} ({})\n", event->GetEventID(), - t, prim->GetParticleDefinition()->GetParticleName(), e, - event->GetPrimaryVertex(0)->GetPosition()[0], - event->GetPrimaryVertex(0)->GetPosition()[1], - event->GetPrimaryVertex(0)->GetPosition()[2], s); + if (TryReservePrimarySlot()) { + // shoot particle + l.fNextActiveSource->GeneratePrimaries(event, l.fCurrentSimulationTime); + // log (after particle creation) + if (LogLevel_EVENT <= GateSourceManager::fVerboseLevel) { + const auto *prim = event->GetPrimaryVertex(0)->GetPrimary(0); + std::string t = G4BestUnit(l.fCurrentSimulationTime, "Time"); + std::string e = G4BestUnit(prim->GetKineticEnergy(), "Energy"); + std::string s = l.fNextActiveSource->fName; + Log(LogLevel_EVENT, fVerboseLevel, + "Event {} {} {} {} {:.2f} {:.2f} {:.2f} ({})\n", + event->GetEventID(), t, + prim->GetParticleDefinition()->GetParticleName(), e, + event->GetPrimaryVertex(0)->GetPosition()[0], + event->GetPrimaryVertex(0)->GetPosition()[1], + event->GetPrimaryVertex(0)->GetPosition()[2], s); + } + } else { + WarnPrimaryLimitReached(); + fRunTerminationFlag = true; + + auto *particle_table = G4ParticleTable::GetParticleTable(); + const auto *particle_def = particle_table->FindParticle("geantino"); + auto *particle = new G4PrimaryParticle(particle_def); + const auto p = G4ThreeVector(); + auto *vertex = new G4PrimaryVertex(p, l.fCurrentSimulationTime); + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); + // (allocated memory is leaked) } } diff --git a/core/opengate_core/opengate_lib/GateSourceManager.h b/core/opengate_core/opengate_lib/GateSourceManager.h index 49fea47f1..6751b7bad 100644 --- a/core/opengate_core/opengate_lib/GateSourceManager.h +++ b/core/opengate_core/opengate_lib/GateSourceManager.h @@ -27,6 +27,8 @@ #include #include #include +#include +#include using namespace indicators; @@ -108,9 +110,17 @@ class GateSourceManager : public G4VUserPrimaryGeneratorAction { void ComputeExpectedNumberOfEvents(); static void SetRunTerminationFlag(bool flag); + static void ResetPrimaryCounterForRun(); + static bool TryReservePrimarySlot(); + static void WarnPrimaryLimitReached(); + static void SetMaxPrimariesPerRun(std::uint64_t value); + static std::uint64_t GetPlatformMaxPrimariesPerRun(); // fRunTerminationFlag should not be thread local static bool fRunTerminationFlag; + static std::atomic fGeneratedPrimariesThisRun; + static std::atomic fPrimaryLimitWarningIssued; + static std::uint64_t fMaxPrimariesPerRun; bool fVisualizationFlag; bool fVisualizationVerboseFlag; std::string fVisualizationType; diff --git a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp index 3a07c8e91..9f1fef3d9 100644 --- a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp +++ b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp @@ -90,9 +90,8 @@ double GateTreatmentPlanPBSource::CalcNextTime(double current_simulation_time) { return next_time; } -double -GateTreatmentPlanPBSource::PrepareNextTime(double current_simulation_time, - double NumberOfGeneratedEvents) { +double GateTreatmentPlanPBSource::PrepareNextTime( + double current_simulation_time, unsigned long NumberOfGeneratedEvents) { if (current_simulation_time < fStartTime) { return fStartTime; diff --git a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h index ceee46dc8..0b11e6d95 100644 --- a/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h +++ b/core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.h @@ -26,7 +26,7 @@ class GateTreatmentPlanPBSource : public GateVSource { void InitializeUserInfo(py::dict &user_info) override; void GeneratePrimaries(G4Event *event, double time) override; double PrepareNextTime(double current_simulation_time, - double NumberOfGeneratedEvents) override; + unsigned long NumberOfGeneratedEvents) override; void PrepareNextRun() override; double CalcNextTime(double current_simulation_time) override; diff --git a/core/opengate_core/opengate_lib/GateVSource.cpp b/core/opengate_core/opengate_lib/GateVSource.cpp index 55d5ee488..5c5b5fec7 100644 --- a/core/opengate_core/opengate_lib/GateVSource.cpp +++ b/core/opengate_core/opengate_lib/GateVSource.cpp @@ -69,7 +69,7 @@ void GateVSource::PrepareNextRun() { } double GateVSource::PrepareNextTime(double current_simulation_time, - double numberOfGeneratedEvents) { + unsigned long numberOfGeneratedEvents) { UpdateActivity(current_simulation_time); if ((fMaxN <= 0) || ((fMaxN > numberOfGeneratedEvents) && (fMaxN > 0))) { if (current_simulation_time < fStartTime) diff --git a/core/opengate_core/opengate_lib/GateVSource.h b/core/opengate_core/opengate_lib/GateVSource.h index 0f6996c69..d121bb92d 100644 --- a/core/opengate_core/opengate_lib/GateVSource.h +++ b/core/opengate_core/opengate_lib/GateVSource.h @@ -39,7 +39,7 @@ class GateVSource { virtual void PrepareNextRun(); virtual double PrepareNextTime(double current_simulation_time, - double NumberOfGeneratedEvents); + unsigned long NumberOfGeneratedEvents); virtual void GeneratePrimaries(G4Event *event, double current_simulation_time); diff --git a/core/opengate_core/opengate_lib/pyGateSourceManager.cpp b/core/opengate_core/opengate_lib/pyGateSourceManager.cpp index 81cd3abbc..978715dd4 100644 --- a/core/opengate_core/opengate_lib/pyGateSourceManager.cpp +++ b/core/opengate_core/opengate_lib/pyGateSourceManager.cpp @@ -20,6 +20,8 @@ void init_GateSourceManager(py::module &m) { .def("SetActors", &GateSourceManager::SetActors) .def("GetExpectedNumberOfEvents", &GateSourceManager::GetExpectedNumberOfEvents) + .def_static("GetPlatformMaxPrimariesPerRun", + &GateSourceManager::GetPlatformMaxPrimariesPerRun) .def_readwrite("fUserEventInformationFlag", &GateSourceManager::fUserEventInformationFlag) .def("StartMasterThread", &GateSourceManager::StartMasterThread, diff --git a/docs/source/developer_guide/developer_guide_notes.rst b/docs/source/developer_guide/developer_guide_notes.rst index a70988fa7..24f30a13d 100644 --- a/docs/source/developer_guide/developer_guide_notes.rst +++ b/docs/source/developer_guide/developer_guide_notes.rst @@ -246,3 +246,23 @@ Geant4’s MT mechanism based on locks needs to be bound to python with a “py::gil_scoped_release release” statement as above. The serial version G4RunManager::Initialize() does not need this statement (and should not have it) because it does not check locks at any point. + +BeamOn() event-count limit in Gate +---------------------------------- + +For the actual simulation run, Gate currently drives Geant4 by calling +``/run/beamOn INT32_MAX`` once per run timing interval from the source manager. +This is intentional: Gate does not usually know the exact number of primaries +in advance, because many sources sample them stochastically over time. +Instead, Gate enters a large Geant4 event loop and terminates the run early +from inside ``GateSourceManager::GeneratePrimaries()`` when the sources are +exhausted. + +This design means the Geant4 ``BeamOn()`` limit matters during the run, not +only during Python-side setup. Since ``G4RunManager::BeamOn()`` takes its event +count as a ``G4int``, one ``BeamOn()`` call cannot safely exceed the maximum +value of that type. Gate therefore keeps a dedicated per-run primary counter in +``GateSourceManager`` and aborts the run with a warning before generating more +than the allowed number of real primaries. The check must live in the runtime +loop, because pre-computing the exact number of primaries is not reliable for +all source types. diff --git a/docs/source/user_guide/user_guide_sources.rst b/docs/source/user_guide/user_guide_sources.rst index 1aebc98fa..6d0a72640 100644 --- a/docs/source/user_guide/user_guide_sources.rst +++ b/docs/source/user_guide/user_guide_sources.rst @@ -53,6 +53,13 @@ Alternatively, you can specify the fixed number of primary particles to be gener Note, however, that certain functionality regarding simulations with multiple runs is currently not supported when sources specify the fixed number of particles as above. +There is an important Geant4 limitation behind the scenes: one ``BeamOn()`` +call can only handle up to the maximum value of Geant4's ``G4int`` type +(typically ``2147483647``). In OpenGATE, this limit applies per run timing +interval, because each interval is executed as a separate Geant4 run. If a run +would generate more primaries than this limit, OpenGATE stops that run +cleanly and emits a warning instead of letting the simulation crash. + How to position the source in space =================================== diff --git a/opengate/engines.py b/opengate/engines.py index 1434cf7ee..b10f219e6 100644 --- a/opengate/engines.py +++ b/opengate/engines.py @@ -68,9 +68,8 @@ class SourceEngine(EngineBase): Source Engine manages the G4 objects of sources at runtime """ - # G4RunManager::BeamOn takes an int as input. The max cpp int value is currently 2147483647 - # Python manages int differently (no limit), so we need to set the max value here. - max_int = 2147483647 + # One Geant4 BeamOn() call is limited by the compiled G4int width. + max_int = g4.GateSourceManager.GetPlatformMaxPrimariesPerRun() def __init__(self, simulation_engine): super().__init__(simulation_engine) @@ -203,6 +202,9 @@ def create_g4_thread_source_manager(self, append=True): self.source_manager_options["progress_bar"] = ( self.simulation_engine.simulation.progress_bar ) + self.source_manager_options["max_primaries_per_run"] = ( + self.simulation_engine.simulation.max_primaries_per_run + ) ms.Initialize(self.run_timing_intervals, self.source_manager_options) self.expected_number_of_events = ( diff --git a/opengate/managers.py b/opengate/managers.py index 48b6128c8..f5099fbfc 100644 --- a/opengate/managers.py +++ b/opengate/managers.py @@ -1841,6 +1841,14 @@ class Simulation(GateObject): "doc": "Display a progress bar during the simulation", }, ), + "max_primaries_per_run": ( + g4.GateSourceManager.GetPlatformMaxPrimariesPerRun(), + { + "doc": "Maximum number of primaries allowed in a single Geant4 BeamOn() call. " + "Defaults to the platform-dependent maximum of Geant4's G4int type. " + "Primarily useful for testing the overflow guard.", + }, + ), "dyn_geom_open_close": ( True, { diff --git a/opengate/tests/src/source/test098_stop_simulation_at_max_primaries.py b/opengate/tests/src/source/test098_stop_simulation_at_max_primaries.py new file mode 100644 index 000000000..21873aa78 --- /dev/null +++ b/opengate/tests/src/source/test098_stop_simulation_at_max_primaries.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +import opengate_core +from opengate.tests import utility + + +def collect_run_counts(simulation_engine): + current_run = simulation_engine.g4_RunManager.GetCurrentRun() + simulation_engine.user_hook_log.append( + { + "number_of_events": current_run.GetNumberOfEvent(), + "number_of_events_to_be_processed": current_run.GetNumberOfEventToBeProcessed(), + } + ) + + +def main(): + paths = utility.get_default_test_paths(__file__, output_folder="test098") + + sim = gate.Simulation() + sim.g4_verbose = False + sim.visu = False + sim.number_of_threads = 1 + sim.random_seed = 123456 + sim.output_dir = paths.output + sim.max_primaries_per_run = 1000 + sim.user_hook_after_run = collect_run_counts + + m = gate.g4_units.m + Bq = gate.g4_units.Bq + MeV = gate.g4_units.MeV + sec = gate.g4_units.s + + sim.world.size = [1 * m, 1 * m, 1 * m] + sim.world.material = "G4_Galactic" + sim.physics_manager.physics_list_name = "QGSP_BERT_EMZ" + + source = sim.add_source("GenericSource", "too_many_primaries") + source.particle = "gamma" + source.activity = 1.0e6 * Bq + source.energy.mono = 0.1 * MeV + source.position.type = "point" + source.direction.type = "iso" + + sim.run_timing_intervals = [(0, 1 * sec)] + sim.run() + + run_info = sim.user_hook_log[0] + produced_events = run_info["number_of_events"] + requested_events = run_info["number_of_events_to_be_processed"] + platform_limit = opengate_core.GateSourceManager.GetPlatformMaxPrimariesPerRun() + is_ok = True + + utility.print_test( + requested_events == platform_limit, + f"Geant4 BeamOn requests the platform G4int maximum: {requested_events}", + ) + is_ok = requested_events == platform_limit and is_ok + + utility.print_test( + produced_events <= sim.max_primaries_per_run + 1, + f"Run stopped at the configured ceiling: processed {produced_events} events for a " + f"primary limit of {sim.max_primaries_per_run} " + f"(the optional +1 is the dummy termination event)", + ) + is_ok = produced_events <= sim.max_primaries_per_run + 1 and is_ok + + utility.print_test( + produced_events < requested_events, + f"Run terminated early before exhausting Geant4 BeamOn: produced {produced_events}, requested {requested_events}", + ) + is_ok = produced_events < requested_events and is_ok + + utility.test_ok(is_ok) + + +if __name__ == "__main__": + main()