From 73df732c28c37c4b9d7b06f8dcc35c1556adcfac Mon Sep 17 00:00:00 2001 From: KAVYA WADHWA Date: Sat, 17 Oct 2026 16:47:39 +0530 Subject: [PATCH] Restore atomic de-excitation; make Generator work multithreaded The custom physics list registered G4RadioactiveDecayPhysics before the EM constructor. It then creates the atomic de-excitation module and initializes it while the physics list is built, before any material-cuts couple exists, so de-excitation stayed inactive in every volume: no fluorescence x-rays or Auger electrons, while printParameters reported them enabled. Register it after G4EmStandardPhysics_option4, as Geant4 requires. In multithreaded mode the master has no primary generator or stacking action, so app.generator raised "Primary generator action is not available". Their state is static, so return the bound classes instead. /gun/ and /gps/ commands issued before their messengers exist (before initialize() in MT) are now kept and applied right after initialization. --- .../include/geant4_application/Application.h | 5 +- src/geant4_application/src/Application.cpp | 50 +++++++++------- src/geant4_application/src/PhysicsList.cpp | 15 +++-- .../src/PrimaryGeneratorAction.cpp | 16 ++--- src/python/module.cpp | 7 ++- tests/test_physics.py | 60 +++++++++++++++++++ tests/test_run.py | 40 +++++++++++++ 7 files changed, 151 insertions(+), 42 deletions(-) create mode 100644 tests/test_physics.py diff --git a/src/geant4_application/include/geant4_application/Application.h b/src/geant4_application/include/geant4_application/Application.h index 371cdf2..fd98adf 100644 --- a/src/geant4_application/include/geant4_application/Application.h +++ b/src/geant4_application/include/geant4_application/Application.h @@ -35,6 +35,9 @@ class Application { static std::unordered_set eventFieldsComplete; static std::unordered_set eventFields; + // /gun/ and /gps/ commands issued before their messengers exist, applied by Initialize + static std::vector pendingGeneratorCommands; + public: Application(); ~Application(); @@ -67,8 +70,6 @@ class Application { static void ListCommands(const std::string& directory); static G4UImanager* GetUIManager(); - const PrimaryGeneratorAction& GetPrimaryGeneratorAction() const; - const StackingAction& GetStackingAction() const; const DetectorConstruction& GetDetectorConstruction() const; static std::filesystem::path GetTemporaryApplicationDirectory(); diff --git a/src/geant4_application/src/Application.cpp b/src/geant4_application/src/Application.cpp index 38d70ce..c9f021b 100644 --- a/src/geant4_application/src/Application.cpp +++ b/src/geant4_application/src/Application.cpp @@ -232,6 +232,12 @@ void Application::Initialize() { runManager->Initialize(); isInitialized = true; + + const auto pending = std::move(pendingGeneratorCommands); + pendingGeneratorCommands.clear(); + for (const auto& command: pending) { + Command(command); + } } py::list Application::Run(const py::object& primaries) { @@ -333,8 +339,29 @@ G4UImanager* Application::GetUIManager() { return ui; } +// The particle gun and the GPS, and so their /gun/ and /gps/ commands, are created +// with the user actions. In multithreaded mode those exist only on the worker threads, +// which Geant4 starts during initialization (the master thread has no primary +// generator). Commands for them issued earlier are kept and applied, in order, right +// after initialization; the master then forwards them to every worker. +static bool IsPendingGeneratorCommand(G4UImanager* ui, const string& command) { + const auto start = command.find_first_not_of(' '); + if (start == string::npos) { + return false; + } + const auto path = command.substr(start, command.find(' ', start) - start); + if (path.rfind("/gun/", 0) != 0 && path.rfind("/gps/", 0) != 0) { + return false; + } + return ui->FindCommand(path.c_str()) == nullptr; +} + void Application::Command(const string& command) { auto ui = GetUIManager(); + if ((pInstance == nullptr || !pInstance->IsInitialized()) && IsPendingGeneratorCommand(ui, command)) { + pendingGeneratorCommands.push_back(command); + return; + } int code = ui->ApplyCommand(command); if (code != 0) { throw runtime_error("Command '" + command + "' failed with code " + to_string(code)); @@ -346,28 +373,6 @@ void Application::ListCommands(const string& directory) { ui->ListCommands(directory.c_str()); } -const PrimaryGeneratorAction& Application::GetPrimaryGeneratorAction() const { - if (runManager == nullptr) { - throw runtime_error("RunManager needs to be set up first"); - } - const auto primaryGeneratorAction = dynamic_cast(runManager->GetUserPrimaryGeneratorAction()); - if (primaryGeneratorAction == nullptr) { - throw runtime_error("Primary generator action is not available"); - } - return *primaryGeneratorAction; -} - -const StackingAction& Application::GetStackingAction() const { - if (runManager == nullptr) { - throw runtime_error("RunManager needs to be set up first"); - } - const auto stackingAction = dynamic_cast(runManager->GetUserStackingAction()); - if (stackingAction == nullptr) { - throw runtime_error("Stacking action is not available"); - } - return *stackingAction; -} - const DetectorConstruction& Application::GetDetectorConstruction() const { if (runManager == nullptr) { throw runtime_error("RunManager needs to be set up first"); @@ -397,3 +402,4 @@ void Application::SetEventFields(const unordered_set& fields) { } unordered_set Application::eventFields = {}; +vector Application::pendingGeneratorCommands = {}; diff --git a/src/geant4_application/src/PhysicsList.cpp b/src/geant4_application/src/PhysicsList.cpp index d94b2cd..57ed723 100644 --- a/src/geant4_application/src/PhysicsList.cpp +++ b/src/geant4_application/src/PhysicsList.cpp @@ -19,17 +19,20 @@ using namespace geant4_app; PhysicsList::PhysicsList(bool optical) : G4VModularPhysicsList() { SetVerboseLevel(1); - // Decay physics - RegisterPhysics(new G4DecayPhysics()); - - // Radioactive decay - RegisterPhysics(new G4RadioactiveDecayPhysics()); - // Electromagnetic physics // RegisterPhysics(new G4EmStandardPhysics()); // RegisterPhysics(new G4EmExtraPhysics()); RegisterPhysics(new G4EmStandardPhysics_option4()); + // Decay physics + RegisterPhysics(new G4DecayPhysics()); + + // Radioactive decay. Must come after the EM constructor: registered first, it + // creates the atomic de-excitation module and initializes it before the + // material-cuts couples exist, which leaves fluorescence and Auger emission + // off in every volume (while /process/em/printParameters reports them on). + RegisterPhysics(new G4RadioactiveDecayPhysics()); + // Hadronic physics RegisterPhysics(new G4HadronElasticPhysics()); RegisterPhysics(new G4HadronPhysicsFTFP_BERT()); diff --git a/src/geant4_application/src/PrimaryGeneratorAction.cpp b/src/geant4_application/src/PrimaryGeneratorAction.cpp index a8d5890..02db982 100644 --- a/src/geant4_application/src/PrimaryGeneratorAction.cpp +++ b/src/geant4_application/src/PrimaryGeneratorAction.cpp @@ -125,25 +125,21 @@ void PrimaryGeneratorAction::SetGeneratorType(const string& type) { } void PrimaryGeneratorAction::SetEnergy(double energy) { - auto ui = Application::GetUIManager(); - ui->ApplyCommand("/gun/energy " + to_string(energy) + " MeV"); + Application::Command("/gun/energy " + to_string(energy) + " MeV"); } void PrimaryGeneratorAction::SetPosition(const array& position) { - auto ui = Application::GetUIManager(); - ui->ApplyCommand("/gun/position " + to_string(position[0]) + " " + to_string(position[1]) + " " + - to_string(position[2]) + " cm"); + Application::Command("/gun/position " + to_string(position[0]) + " " + to_string(position[1]) + " " + + to_string(position[2]) + " cm"); } void PrimaryGeneratorAction::SetDirection(const array& direction) { - auto ui = Application::GetUIManager(); - ui->ApplyCommand("/gun/direction " + to_string(direction[0]) + " " + to_string(direction[1]) + " " + - to_string(direction[2])); + Application::Command("/gun/direction " + to_string(direction[0]) + " " + to_string(direction[1]) + " " + + to_string(direction[2])); } void PrimaryGeneratorAction::SetParticle(const string& particle) { - auto ui = Application::GetUIManager(); - ui->ApplyCommand("/gun/particle " + particle); + Application::Command("/gun/particle " + particle); } string PrimaryGeneratorAction::generatorType = "gun"; diff --git a/src/python/module.cpp b/src/python/module.cpp index aa6441b..596ce53 100644 --- a/src/python/module.cpp +++ b/src/python/module.cpp @@ -38,9 +38,12 @@ PYBIND11_MODULE(_geant4_application, m) { .def_static("get_event_fields_complete", &Application::GetEventFieldsComplete) .def_static("command", &Application::Command, py::arg("command")) .def_static("list_commands", &Application::ListCommands, py::arg("directory") = "/") - .def_property_readonly("generator", &Application::GetPrimaryGeneratorAction, py::return_value_policy::reference_internal) + // The generator and stacking actions belong to the worker threads in + // multithreaded mode (the master has none), and all their state is static, + // so these return the classes rather than an instance. + .def_property_readonly("generator", [](const Application&) { return py::type::of(); }) .def_property_readonly("detector", &Application::GetDetectorConstruction, py::return_value_policy::reference_internal) - .def_property_readonly("stacking", &Application::GetStackingAction, py::return_value_policy::reference_internal); + .def_property_readonly("stacking", [](const Application&) { return py::type::of(); }); py::class_(m, "PrimaryGeneratorAction") .def_property_static( diff --git a/tests/test_physics.py b/tests/test_physics.py new file mode 100644 index 0000000..6d2a6d0 --- /dev/null +++ b/tests/test_physics.py @@ -0,0 +1,60 @@ +from __future__ import annotations + +import awkward as ak +import numpy as np +import pytest + +from geant4_python_application import Application + +gold_slab_gdml = """ + + + + + + + + + + + + + + + + + + + + + + +""" + + +@pytest.mark.parametrize("n_threads", [0, 2]) +def test_gold_fluorescence(n_threads): + # The default physics list enables fluorescence; it was silently inactive while + # G4RadioactiveDecayPhysics was registered before the EM constructor. + if n_threads > 0 and not Application.multithreading_available(): + pytest.skip("linked Geant4 runtime was built without multithreading") + with Application(gdml=gold_slab_gdml, seed=1234, n_threads=n_threads) as app: + app.generator.particle("gamma").energy(100, "keV") + app.generator.position(0, 0, -10, "cm").direction(0, 0, 1) + events = app.run(500) + + particle = ak.flatten(events.track.particle) + creator = ak.flatten(events.track.creator_process) + energy = np.asarray(ak.flatten(events.track.initial_energy)) # keV + + photoelectrons = ak.sum((particle == "e-") & (creator == "phot")) + fluorescence = energy[np.asarray((particle == "gamma") & (creator == "phot"))] + # Gold K-alpha lines are at 66.99 and 68.80 keV (K-shell binding 80.7 keV). + k_alpha = fluorescence[(fluorescence > 66) & (fluorescence < 70)] + + # About 60% of the photons interact in 0.1 mm of gold, almost all by + # photoabsorption, 80% of it in the K shell, with a 96% fluorescence yield. + assert photoelectrons > 250 + assert len(k_alpha) > 100 + assert np.all(fluorescence < 100) diff --git a/tests/test_run.py b/tests/test_run.py index eac8267..8de16f0 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -47,6 +47,46 @@ def test_general_particle_source(): assert np.allclose(ak.flatten(events.primaries.energy), 2000) +@pytest.mark.parametrize("initialize_first", [False, True]) +def test_generator_multithreaded(initialize_first): + # In multithreaded mode only the workers have a primary generator; commands + # issued before they exist are applied once the application is initialized. + if not Application.multithreading_available(): + pytest.skip("linked Geant4 runtime was built without multithreading") + with Application(gdml=basic_gdml, n_threads=2, seed=1) as app: + if initialize_first: + app.initialize() + assert app.generator.type == "gun" + app.generator.particle("proton").energy(50, "MeV") + app.generator.position(0, 0, -10, "cm").direction(0, 0, 1) + events = app.run(10) + assert len(events) == 10 + assert set(ak.flatten(events.primaries.particle)) == {"proton"} + assert np.allclose(ak.flatten(events.primaries.energy), 50000) + + +def test_general_particle_source_multithreaded(): + if not Application.multithreading_available(): + pytest.skip("linked Geant4 runtime was built without multithreading") + with Application(gdml=basic_gdml, n_threads=2, seed=1) as app: + app.generator.use_gps().particle("geantino").energy(2, "MeV") + app.generator.position(0, 0, -10, "cm").direction(0, 0, 1) + assert app.generator.type == "gps" + events = app.run(4) + assert len(events) == 4 + assert set(ak.flatten(events.primaries.particle)) == {"geantino"} + assert np.allclose(ak.flatten(events.primaries.energy), 2000) + + +def test_generator_command_typo_multithreaded(): + if not Application.multithreading_available(): + pytest.skip("linked Geant4 runtime was built without multithreading") + with Application(gdml=basic_gdml, n_threads=2) as app: + app.command("/gun/energyy 1 MeV") # kept until the generator exists + with pytest.raises(RuntimeError, match="/gun/energyy"): + app.initialize() + + def test_multiple_primary_vertices_per_event(): particle = lambda name, x: { "particle": name,