Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion core/opengate_core/opengate_lib/GateDebugSource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion core/opengate_core/opengate_lib/GateDebugSource.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
5 changes: 3 additions & 2 deletions core/opengate_core/opengate_lib/GateGenericSource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
2 changes: 1 addition & 1 deletion core/opengate_core/opengate_lib/GateGenericSource.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
5 changes: 3 additions & 2 deletions core/opengate_core/opengate_lib/GateLastVertexSource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion core/opengate_core/opengate_lib/GateLastVertexSource.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
101 changes: 86 additions & 15 deletions core/opengate_core/opengate_lib/GateSourceManager.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#if defined(_MSC_VER)
#pragma warning(pop)
#endif
#include <G4Exception.hh>
#include <G4MTRunManager.hh>
#include <G4RunManager.hh>
#include <G4TransportationManager.hh>
Expand All @@ -25,6 +26,7 @@
#include <G4UnitsTable.hh>
#include <cmath>
#include <iostream>
#include <limits>
#include <pybind11/numpy.h>

#ifdef USE_GDML
Expand All @@ -34,6 +36,9 @@
/* There will be one SourceManager per thread */

bool GateSourceManager::fRunTerminationFlag = false;
std::atomic<std::uint64_t> GateSourceManager::fGeneratedPrimariesThisRun{0};
std::atomic<bool> GateSourceManager::fPrimaryLimitWarningIssued{false};
std::uint64_t GateSourceManager::fMaxPrimariesPerRun = INT32_MAX;

GateSourceManager::GateSourceManager() {
fUIEx = nullptr;
Expand Down Expand Up @@ -66,6 +71,50 @@ void GateSourceManager::SetRunTerminationFlag(bool flag) {
fRunTerminationFlag = flag;
}

std::uint64_t GateSourceManager::GetPlatformMaxPrimariesPerRun() {
return std::numeric_limits<G4int>::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;
Expand All @@ -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();

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)) {
Expand Down Expand Up @@ -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)
}
}

Expand Down
10 changes: 10 additions & 0 deletions core/opengate_core/opengate_lib/GateSourceManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@
#include <G4UIsession.hh>
#include <G4VUserPrimaryGeneratorAction.hh>
#include <G4VisExecutive.hh>
#include <atomic>
#include <cstdint>

using namespace indicators;

Expand Down Expand Up @@ -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<std::uint64_t> fGeneratedPrimariesThisRun;
static std::atomic<bool> fPrimaryLimitWarningIssued;
static std::uint64_t fMaxPrimariesPerRun;
bool fVisualizationFlag;
bool fVisualizationVerboseFlag;
std::string fVisualizationType;
Expand Down
5 changes: 2 additions & 3 deletions core/opengate_core/opengate_lib/GateTreatmentPlanPBSource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
2 changes: 1 addition & 1 deletion core/opengate_core/opengate_lib/GateVSource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion core/opengate_core/opengate_lib/GateVSource.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 2 additions & 0 deletions core/opengate_core/opengate_lib/pyGateSourceManager.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
20 changes: 20 additions & 0 deletions docs/source/developer_guide/developer_guide_notes.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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.
7 changes: 7 additions & 0 deletions docs/source/user_guide/user_guide_sources.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
===================================
Expand Down
8 changes: 5 additions & 3 deletions opengate/engines.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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 = (
Expand Down
8 changes: 8 additions & 0 deletions opengate/managers.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
{
Expand Down
Loading
Loading