Restore atomic de-excitation; make Generator work multithreaded

src/geant4_application/include/geant4_application/Application.h

@@ -35,6 +35,9 @@ class Application {
     static std::unordered_set<std::string> eventFieldsComplete;
     static std::unordered_set<std::string> eventFields;
 
+    // /gun/ and /gps/ commands issued before their messengers exist, applied by Initialize
+    static std::vector<std::string> 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();

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<const PrimaryGeneratorAction*>(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<const StackingAction*>(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<string>& fields) {
 }
 
 unordered_set<string> Application::eventFields = {};
+vector<string> Application::pendingGeneratorCommands = {};

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());

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<double, 3>& 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<double, 3>& 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";

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<PrimaryGeneratorAction>(); })
             .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<StackingAction>(); });
 
     py::class_<PrimaryGeneratorAction>(m, "PrimaryGeneratorAction")
             .def_property_static(

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 = """<?xml version="1.0" encoding="UTF-8"?>
+<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+      xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">
+    <solids>
+        <box name="WorldSolid" x="100" y="100" z="100" lunit="cm"/>
+        <box name="SlabSolid" x="10" y="10" z="0.1" lunit="mm"/>
+    </solids>
+    <structure>
+        <volume name="Slab">
+            <materialref ref="G4_Au"/>
+            <solidref ref="SlabSolid"/>
+        </volume>
+        <volume name="World">
+            <materialref ref="G4_Galactic"/>
+            <solidref ref="WorldSolid"/>
+            <physvol name="slab">
+                <volumeref ref="Slab"/>
+            </physvol>
+        </volume>
+    </structure>
+    <setup name="Default" version="1.0">
+        <world ref="World"/>
+    </setup>
+</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)

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,