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,