aegir
Phlex-based simulation framework for the SHiP experiment.
Loading...
Searching...
No Matches
geant4_integration

Geant4 integration strategy

aegir uses a direct worker integration pattern where Phlex framework threads act as Geant4 worker threads directly, with no separate G4 thread pool or event queue. This is inspired by the CMSSW OscarMTProducer pattern used in CMS.

Why direct workers?

We evaluated five integration strategies in the ship-phlex-sim benchmarking repository:

  1. Synchronous single-threaded – G4RunManager::BeamOn(1) on a dedicated thread.
  2. Asynchronous single-threaded – dedicated G4 thread with an event queue and std::future for synchronisation.
  3. Synchronous multi-threaded (slots) – G4TaskRunManager with per-worker rendezvous slots.
  4. Asynchronous multi-threaded (queue) – G4MTRunManager on a dedicated thread; G4 workers pop from a concurrent queue.
  5. Direct workers – framework threads are G4 workers; no queue, no BeamOn.

Strategies 1–4 all involve cross-thread handoff overhead (queues, promises, condition variables). The direct worker approach eliminates this entirely: the thread that receives the event from the framework is the same thread that processes it through Geant4.

In benchmarks with particle gun events on a 4-core machine, the direct pattern matched the best queue-based strategy (sync MT slots) while being simpler to configure and reason about. It also removes the num_events configuration footgun present in the queue-based MT module, where the configured event count had to exactly match the driver's total.

Architecture

Master thread (std::thread)
|
+-- Owns G4MTRunManager
+-- Runs Initialize() + RunInitialization()
+-- Stores world volume, physics list, detector construction
+-- Blocks until shutdown
Phlex TBB worker threads (N = concurrency)
|
+-- Lazy init: G4WorkerRunManagerKernel per thread
+-- simulate() builds G4Event from MCParticles
+-- Calls G4EventManager::ProcessOneEvent()
+-- Collects results from thread-local storage

Master initialisation

A dedicated std::thread creates a G4MTRunManager and runs the full Geant4 initialisation sequence (geometry, physics, run initialisation). It then stores the world physical volume, physics list, and detector construction pointers for workers to use. The master thread blocks on a shutdown future until the module is destroyed.

Worker initialisation

Each Phlex thread lazily initialises on its first call to simulate():

  1. Assigns a unique G4 thread ID via G4Threading::G4SetThreadId()
  2. Calls G4WorkerThread::BuildGeometryAndPhysicsVector()
  3. Creates a G4WorkerRunManagerKernel and defines the world volume
  4. Initialises physics and sensitive detectors for this thread
  5. Sets user actions (tracking, energy cut) on the event manager

Event processing

Instead of G4RunManager::BeamOn(), the module:

  1. Builds a G4Event directly from the input MCParticle vector — keeping only final-state entries, see Primary selection — creating G4PrimaryVertex and G4PrimaryParticle objects (no G4VUserPrimaryGeneratorAction involved)
  2. Sets the G4 state to G4State_GeomClosed
  3. Calls G4EventManager::ProcessOneEvent(event)
  4. Collects hits and particles from thread-local storage

This bypasses the G4 run loop entirely, giving the framework full control over event scheduling.

Primary selection

SHiP::MCParticle::status carries the HepMC status code — the same convention Pythia8's Particle::statusHepMC() produces and EDM4hep uses for its generatorStatus:

Code Meaning Tracked?
0 Empty entry, no meaningful information no
1 Final state — not decayed by the generator; may still be unstable yes
2 Decayed Standard Model hadron, tau or muon no
3 Documentation entry no
4 Incoming beam particle no
11-200 Intermediate entry, generator-dependent classification no

Only status-1 entries become primaries. Note that "final state" does not mean stable: a final-state K_S is status 1 and is expected to decay in the detector — that is exactly the work being handed over. Everything else is either the generator's own bookkeeping or a particle the generator has already decayed, and tracking it alongside the daughters it already produced would double-count the event.

The predicate lives in src/hepmc_status.hpp, next to the code set. Geant4's own HepMC example uses the same status == 1 test.

Skipped primaries are counted and reported once per event in an aggregated warning; they remain in the output mc_particles product, so the record is not lost.

This filter is a no-op for a generator running in the default ‘record: 'final_state’mode, where every emitted entry is status 1. It matters underrecord: 'full'` — see generator_record.md.

Without it the failure mode is silent and partial rather than loud: Geant4 declines to track short-lived definitions (quarks, gluons, diquarks, strings) on its own, but happily tracks decayed hadrons and beam particles. Setting track_all_primaries on a full record demonstrates this — the primary count rises to exactly the number of status 1, 2 and 4 entries.

Shutdown

The G4MTRunManager is intentionally leaked at shutdown. Its destructor accesses global singletons that may already be torn down during plugin unloading, causing crashes. This is a known Geant4 limitation.

Configuration

Parameter Type Default Description
physics_list string FTFP_BERT Geant4 reference physics list name
concurrency int 1 Number of concurrent worker threads
verbosity int 0 Geant4 run verbosity level
sd_mode string scoring Sensitive detector mode: scoring or crossing
ke_threshold double 0.0 Kinetic energy threshold for crossing SD (GeV)
energy_cut bool false Enable stepping energy cut
energy_cut_threshold double ke_threshold KE below which tracks are killed (GeV)
particle_ke_cut double 0.0 KE below which secondary particles are not recorded (GeV)
track_all_primaries bool false Hand every input MCParticle to Geant4 regardless of its generator status. By default only HepMC status-1 (final-state) entries are tracked; set this only for inputs whose status does not follow the HepMC convention
regions map {} Volume name pattern to production cut (mm) mapping
scoring_planes list [] Scoring planes to add to the geometry, see Scoring planes
export_gdml string *(unset)* Write the constructed geometry to this GDML file after initialisation. Errors if the file exists. Lets external tools (e.g. the GENIE event generator) use exactly the geometry Geant4 tracks in
progress_interval int 100 Log a progress line (event count and average rate) every this many simulated events; 0 disables

For consumers that read the exported file with ROOT's TGeo importer, run scripts/gdml_fix_element_names.py on it first: ROOT confuses materials and elements that share a name after pointer-suffix stripping (Geant4's NIST materials routinely do — material "Iron" made of element "Iron") and silently imports them as empty mixtures. The script renames the colliding elements; the physics is unchanged. scripts/gdml_target_nuclei.py then lists the nuclide PDG codes of the geometry — e.g. the target list neutrino cross-section splines must cover.

Example workflow configuration:

{
modules: {
geant4: {
cpp: 'geant4_module',
physics_list: 'FTFP_BERT',
concurrency: 4,
sd_mode: 'crossing',
ke_threshold: 0.5,
energy_cut: true,
particle_ke_cut: 1.0,
regions: { '/SHiP/target': 50, '/SHiP/muon_shield/magn_absorb': 50 },
},
},
}

Scoring planes

A scoring plane is a thin box or cylindrical shell that records the particles crossing it. It is aegir's counterpart of FairShip's exitHadronAbsorber and is meant for flux studies: the muons leaving the hadron absorber, the neutrinos coming out of the target, and so on.

The planes live in a Geant4 parallel world (aegir_scoring_planes), overlaid on whatever geometry the provider built — builtin, GDML or GeoModel alike. Geant4 tracks every particle through both worlds at once (G4ParallelWorldPhysics, added to the physics list whenever planes are configured). A plane can therefore be any size, anywhere: it can cut through or extend beyond the detector's volumes — a muon-shield container, a magnet, the target — and the parallel world carries no material of its own, so the detector's materials apply everywhere. Planes must not overlap one another (touching is fine) or extend beyond the world volume.

A plane records one SimHit for every track that enters it through its surface, whatever the particle — neutrinos and other neutrals included — and whether or not it deposits energy. It does this through its own sensitive detector, so it works the same with ‘sd_mode: 'scoring’and sd_mode: 'crossing', alongside the geometry's own sensitive volumes, and thesensitive_volumesand fieldvolume_patternpatterns never match the planes. A track entering a plane twice (a cylinder crossed in and out, a particle scattered back) is recorded each time; setkill` to stop a track at the first plane it reaches.

Planes do change how the detector geometry is stepped through. With a parallel world registered, Geant4 uses coupled transportation, and steps in the detector geometry also end at every plane face. Hence:

  • ‘sd_mode: 'scoring’writes one hit per step, so a plane cutting through a sensitive volume splits some of its hits in two; the total energy deposit is unchanged. -sd_mode: 'crossing'` is unaffected, since its first-step-in-volume flag only follows the detector's own boundaries.
  • Runs with and without planes are not identical event by event, even with the same seed.
  • Each step costs some extra CPU time, so only configure the planes a run needs.

The hits are told apart from detector hits by their ids:

SimHit field Scoring-plane value
detectorId 100 (SHiP::g4::kScoringPlaneDetectorId), clear of the SHiP subsystem ids and of the sensitive_volumes indices (which limits a geometry to 100 sensitive_volumes patterns)
geometryNodeId the plane's id
position, momentum, time at the entry point
energyDeposit 0
pathLength track length from the track's origin to the entry point. A ‘sd_mode: 'crossing’` hit gives it at the end of its first step in the volume instead, so the two don't mean quite the same thing

In the output module's validation histograms, plane hits fill their own h_plane_hit_* histograms (hits per event, plane id, momentum, PDG code) and are left out of the h_hit_* ones, which keep describing the detector's hits: with detector id 100 and no energy deposit, plane hits would otherwise land in the overflow of h_hit_detector and pile up at zero in h_hit_edep. h_plane_hit_id covers plane ids 0–99; other ids land in its underflow or overflow bin.

Each entry of scoring_planes takes:

Key Type Default Description
name string — Volume name; must not already exist in the geometry
id int — Written to SimHit::geometryNodeId, telling the planes' hits apart; required and unique
shape string box box or tube (cylindrical shell, axis along z)
half_x, half_y double — Box half-widths (mm); required for a box, rejected for a tube
r_min, r_max double — Tube inner and outer radius (mm); required for a tube, rejected for a box
half_z double — Half-length along z (mm); required for both shapes
position [x, y, z] — Centre in the global frame (mm)
after_volume string — Place the plane downstream of this detector volume, centred transversely on its global bounding box
gap double 0 With after_volume: distance from the volume's downstream end to the plane's upstream face (mm)
ke_threshold double 0 Record only tracks entering with at least this kinetic energy (GeV). With kill, tracks below it are stopped too, unrecorded
kill bool false Stop and kill every track that steps into the plane, as FairShip's exitHadronAbsorber does: recorded or not, including secondaries born inside the plane

Exactly one of position and after_volume must be given. after_volume is an exact logical-volume name, and the volume must be placed exactly once in the geometry — otherwise the plane's position would be ambiguous. Planes are always aligned with the global axes. Unknown keys and inconsistent settings fail the job at configuration time. A plane that overlaps another plane or extends beyond the world fails it at initialisation, before any plane is built. That test is exact for these axis-aligned shapes, so even a small corner overlap between thin planes is caught.

geant4: lib.geant4 {
scoring_planes: [
{
// FairShip's plane 1 cm behind the hadron absorber
name: 'exit_hadron_absorber',
id: 1,
half_x: 3559.0, half_y: 1699.0, half_z: 1.0,
after_volume: '/SHiP/muon_shield/magn_absorb',
gap: 10.0,
kill: true,
},
{
// cylindrical shell around a target centred at the origin
name: 'around_target',
id: 2,
shape: 'tube',
r_min: 80.0, r_max: 81.0, half_z: 500.0,
position: [0.0, 0.0, 0.0],
},
],
}

workflows/gun_st_scoring_planes.jsonnet places planes of each kind in the builtin geometry, one of them cutting through the tungsten target, and scripts/scoring_planes.sh (ctest scoring_planes) checks the recorded hits against the geometry for neutrinos and muons.