Skip to content

Geant4 features as uLib objects

Every Geant4 feature below is a uLib::Object: its parameters are discoverable properties, so it appears in the Properties panel and the node graph, it saves into a .gcx like anything else, and it connects to other objects through ports.

This page walks through the five examples in src/HEP/Geant/examples/, explaining the code that matters in each. The examples are built and run by ctest, so an example that stopped being true fails the suite rather than misleading the next reader. Each also writes the scene it built; the copies in assets/examples/geant/ are those same files.

Example Shows Scene
01 recording every step in a volume 01_sensitive_detector.gcx
02 accumulating totals, and narrowing what counts 02_scoring_with_filters.gcx
03 bending tracks, capping steps 03_fields_and_limits.gcx
04 placing many copies cheaply 04_replicas_and_assemblies.gcx
05 ghost geometry, mirrored geometry 05_parallel_and_reflection.gcx
06 picking a material, and building one 06_materials_and_compounds.gcx
07 what happens when the numbers describe nothing 07_invalid_geometry.gcx
08 naming the physics a run was done with 08_choosing_physics.gcx
09 scoring on a grid rather than per volume 09_dose_map.gcx
10 choosing which tracks are worth following 10_track_filtering.gcx
11 optical properties, and Birks quenching 11_optical_scintillator.gcx
12 truth hits into what an instrument would report 12_digitizing_a_readout.gcx
13 isotopes that are not the natural mix 13_enriched_materials.gcx
14 slicing a mother, and how that differs from a replica 14_dividing_a_volume.gcx
15 run statistics, and one line per track 15_watching_a_run.gcx
16 an image as geometry, and measuring it back out 16_a_voxel_phantom.gcx
17 a plane that watches muons, a recipe, a spectrum 17_a_plane_probe.gcx
19 muons streamed into an image as they are measured 19_streaming_poca.gcx
20 the reconstruction waits, and the waiting starts the run 20_driven_stream.gcx
21 a run's muons into a TFile, with the scene beside them 21_muons_to_root.gcx

Run one directly to see its narration, and point it at a path to write the scene somewhere else:

./build/clang-ninja/src/HEP/Geant/examples/01_sensitive_detector /tmp/mine.gcx

The shape they all share

Three conventions recur, and knowing them makes every example shorter to read.

Ports that grow. Anything applying to a set of volumes — a region, a field, a detector — has one port per attached object plus one spare on the end. The spare is how the next object is attached; without it something holding nothing would have no port to wire to at all.

Attach(detector, "Volume", 0, planes[0]);   // fills Volume0, creates Volume1
Attach(detector, "Volume", 1, planes[1]);   // fills Volume1, creates Volume2

Attach is example plumbing (ExampleSupport.h); underneath it is an ordinary reference port:

auto *port = dynamic_cast<ReferencePropertyBase *>(
    owner->GetProperty("Volume0"));
port->SetReferencedObject(target);

Off is negative, not zero. Every optional cap — a production cut, a step limit — is off when negative. Zero is a real instruction that would stop a track dead, so it cannot double as "unset", and making it do so would turn the value a user is most likely to type by accident into one that silently does nothing.

A carrier is a volume or a region. Geant4 lets both G4LogicalVolume and G4Region hold a field manager and a user limit set, so fields and limits accept either. Sensitive detectors and scoring detectors take volumes only — Geant4 has no notion of a sensitive region.


01 — Sensitive detectors

A sensitive detector is told about every step taken inside the volumes attached to it. Use it when the individual interactions matter: where a muon crossed a plane, how much a strip saw, in what order.

It sits beside the existing DetectorSteppingAction, which reconstructs muon tracks across the chamber planes. That is a different job from recording per-volume hits, and both can run in one simulation.

Attaching is the whole of the wiring — the volume becomes sensitive as a consequence:

auto *detector = new Geant::SensitiveDetector("tracker");
scene.AddObject(detector);
for (int i = 0; i < 3; ++i) Attach(detector, "Volume", i, planes[i]);

Why the rules are public

The recording policy lives in RecordHit(), deliberately apart from the G4Step that normally supplies the numbers:

Geant::Hit hit;
hit.energyDeposit = 1.8;      // MeV
hit.pdgCode = 13;             // a muon
detector->RecordHit(hit);

ProcessHits() is only the part that pulls those numbers off a step. Splitting them means the thresholds, the particle selection and the hit cap can be exercised — and shown — without standing up a beam. The example does exactly that, and so does SensitiveTest.cpp.

What the thresholds do

detector->MinEnergyDeposit.Set(0.5);    // MeV
detector->RecordHit(crossing(0.05, 13, 1));   // dropped
detector->RecordHit(crossing(1.80, 13, 1));   // kept

Rejected hits are counted, not forgotten: GetRejectedCount() is what tells "my detector saw nothing" apart from "my detector was never reached".

Nothing is filtered by default, including a step that deposits zero — a particle passing through is still a crossing, and whether that matters belongs to the consumer rather than the detector.

MaxHits caps memory on a long run; past the cap, arrivals are counted as rejected rather than silently lost.

Lifetime

A G4LogicalVolume keeps a bare pointer to its detector and does not own it, so the detector has to come off before it dies. It detaches through the G4 pointers it recorded, not by re-resolving the uLib objects:

m_Applied.ApplySensitiveDetector(nullptr);

A G4LogicalVolume belongs to G4LogicalVolumeStore and outlives the LogicalVolume that made it, so this still works during scene teardown — when the objects beside it may already be gone. An earlier version re-resolved the carriers with dynamic_cast and segfaulted for exactly that reason.


02 — Scoring with filters

Where a sensitive detector keeps a record per step, a scorer keeps a running total per volume and throws the steps away. For a dose map the total is the answer, at a double per volume rather than a record per step.

Three objects make a setup: a ScoringDetector holding any number of Scorers, each with an optional filter.

One enum instead of nine classes

auto *dose = new Geant::Scorer("dose");
dose->Quantity.Set(static_cast<int>(Geant::ScoredQuantity::DoseDeposit));

Geant4 offers a class per scorer. Here the quantity is a property, so the catalogue is one drop-down to read rather than nine classes to find. ScoringTest.cpp walks every label and builds each one, so no entry in that drop-down can become a dead option.

Changing the quantity discards the old primitive rather than reconfiguring it: a total that switched quantity halfway through a run would mean nothing.

One volume, one detector — many scorers

Attach(detector, "Scorer", 0, dose);
Attach(detector, "Scorer", 1, path);
Attach(detector, "Volume", 0, block);

A volume can carry only one sensitive detector, so several numbers on one volume means several scorers on one G4MultiFunctionalDetector. That is what the detector is for.

Filters are objects

auto *muons = new Geant::ParticleFilter("muons_only");
muons->Particles.Set("mu-,mu+");
dose->SetFilter(...);      // or wire the "Filter" port

Making the filter its own object means one can serve several scorers — "muons only" applied to the dose and the path length at once, edited in one place.

A name no physics list ever built is a fatal G4Exception inside G4SDParticleFilter::add(), which by default calls abort(). Since this is a string a user types, the library catches it (Guard.h), names the offending entry and skips it. Example 02 deliberately includes a typo so you can see the refused line and observe that the program continues.

Reading the numbers

std::map<std::int32_t, double> scores = dose->ReadScores(event);

Scores live in the event, not in the scorer: a primitive scorer fills a hits collection the event owns and keeps nothing of its own between events. Any other way of reading either reaches into G4VPrimitiveScorer's protected state or returns whichever event happened to be last. Before a run, ReadScores(nullptr) is empty — a meaningful answer, not an error.

Scores are keyed by copy number, which is what separates one replica slice from another; they all share a name, so the name cannot.


03 — Fields and limits

A field is more than a field

A field on its own moves nothing. Geant4 needs an equation of motion for it, a stepper to integrate that, a driver to choose step sizes, and a chord finder to decide when a curve is straight enough to treat as a segment — all held by a G4FieldManager. Field owns that chain; subclasses supply the parts that depend on the kind of field.

auto *field = new Geant::MagneticField("dipole");
field->FieldValue.Set(Vector3d(0.0, 1.5, 0.0));   // tesla
Attach(field, "Carrier", 0, bore);

A magnetic field does no work, so Geant4's default equation and stepper are correct and G4ChordFinder builds both. An electric field changes a particle's energy, so ElectricField builds the chain explicitly — G4EqMagElectricField over eight variables rather than six.

A zero field is a normal state meaning "nothing to apply", not a failure — it is what you get the moment one is created. IsBuilt() reports false and nothing is attached.

The accuracy properties govern integration; the tightest is DeltaIntersection, because a misplaced boundary crossing puts energy in the wrong volume. Changing any of them rebuilds the chain, and the carriers are re-pointed at the new manager rather than left on the freed one.

Regions and limits

auto *region = new Geant::Region("tracker_region");
region->GammaCut.Set(0.1);          // mm
Attach(region, "Volume", 0, tracker);

auto *limits = new Geant::UserLimits("fine_steps");
limits->MaxStep.Set(1.0);           // mm
Attach(limits, "Carrier", 0, region);

MaxStep is the limit usually wanted. Geant4 will happily cross a metre of scintillator in a single step, and a detector recording per step then records one hit for the whole traversal. Capping the step is what makes the response resolvable along the track.

Attaching the limits to the region rather than to each volume is the point of a region: it covers everything in it, including volumes added later.


04 — Replicas, parameterised volumes and assemblies

A PhysicalVolume places one thing once. At a thousand copies that costs a thousand placements, and Geant4 offers three ways out. Which applies depends on how the copies are arranged.

Requirement Cost
ReplicaVolume slices tile the mother exactly one volume, position computed arithmetically
ParameterisedVolume none — placement is a rule one volume, placement computed per copy
AssemblyVolume none — a group placed repeatedly one placement per part per imprint
slices->Axis.Set(static_cast<int>(Geant::DivisionAxis::Z));
slices->Count.Set(20);
slices->Width.Set(100.0);   // 20 x 100 mm fills the 2000 mm mother

The only-daughter rule

A replica or parameterised volume must be its mother's only daughter. Two consequences, both of which the example demonstrates:

Rebuilding must remove the old placement first. A G4 placement adds itself to its mother's daughter list on construction, so a rebuild that merely forgot the previous pointer would leave two — and Geant4 then refuses the new one. The library removes and deletes the previous placement on every rebuild, so changing Count works:

slices->Count.Set(40);
slices->Width.Set(50.0);
slices->Update();           // mother still has exactly one daughter

Nothing else may share that mother. Example 04 puts its assembly in its own ring volume rather than in the hall that already holds the ladder. Getting this wrong produces the most confusing error in this area, because the refusal names the volume being placed rather than the parameterised one that caused it. The first draft of this example made exactly that mistake and failed the suite.

Assemblies

MakeBoxVolume(*sector->GetChildren(), "scintillator", 100.0, ...);
MakeBoxVolume(*sector->GetChildren(), "lightguide", 60.0, ...);
sector->Imprints.Set(6);

An assembly is not a volume — it is a group placed together, repeatedly. Six imprints of two parts is twelve placements, which is the point of an assembly over placing each by hand. GetParts() picks the LogicalVolumes out of the children and ignores the solids and materials beside them.


05 — Parallel worlds and reflections

Both are narrow tools that are easy to mistake for something broader.

A parallel world asks questions, it does not add geometry

The volumes in it are ghosts: a particle passes through without interacting, but the navigator still reports the crossing, so a scorer or sensitive detector on a ghost measures what went through it.

auto *readout = new Geant::ParallelWorld("readout");
Attach(readout, "Ghost", 0, strip);

The use is a readout finer than the physical structure — a single slab read out in strips. Cutting the slab into strips would change the physics; laying a ghost grid over it does not.

Ghosts are placed by Construct(), which the run manager calls during initialisation and this library does not. So GetG4ParallelWorld() hands the object out to be registered rather than installing it, and GetPlacedCount() answers zero until a run has initialised. That zero is the honest answer, not a failure, and the test asserts it.

A reflection is not a negative scale

Geant4 will not accept a negative scale on a placement. A reflected solid is genuinely a different solid with its own surface normals, and a navigator handed a negative determinant gets inside and outside the wrong way round. G4ReflectionFactory decomposes the transform, builds the mirrored solid, and places what is left.

mirrored->MirrorX.Set(true);
mirrored->Update();
G4LogicalVolume *reflectedLV = mirrored->GetReflectedLogical();

Two traps, both named in the API so they cannot be hit silently:

GetReflectedLogical() is the mirror image — a new G4LogicalVolume over a new solid, asked of the factory's map rather than remembered, so a rebuild cannot leave it stale.

GetPlacementInReflectedMother() is not the mirror image, despite being the second half of what G4ReflectionFactory::Place() returns. It is the twin placement needed when the mother is itself reflected, and is null otherwise. The first draft of ParallelTest.cpp asserted it was the mirror and failed, which is how the distinction earned its name.

Parity, not presence:

mirrored->MirrorX.Set(true);
mirrored->MirrorY.Set(true);
mirrored->IsMirrored();     // false -- this is a rotation about Z

Two mirrors compose to a rotation about the third axis: positive determinant, no mirrored solid, an ordinary placement. IsMirrored() counts parity rather than asking whether any axis was ticked. Ticking two boxes is legal and almost never what was meant.


06 — Materials and compounds

Two ways to say what something is made of.

Picking from the catalogue

A material's Name carries the whole NIST catalogue as its choices:

auto *air = new Geant::Material("G4_AIR");
air->SetName("G4_WATER");

In gcompose that renders as a filterable picker rather than a text field — 309 materials on a typical build, so the field is editable with a contains-match completer: type poly and the list narrows. The catalogue is read from Geant4, not written out in uLib, so it is right for whatever build you are running.

It is a suggestion, not a constraint. A name the catalogue does not have is kept, with a note:

a mistyped name stays put: 'G4_WATTER', with the note:
    'G4_WATTER' is not in the NIST catalogue

Silently substituting G4_AIR — which is what this used to do — is a substitution you do not see, and it surfaces much later as the wrong physics rather than as a bad name.

Building one that is not in the catalogue

A scintillator recipe, an alloy, a borated concrete. Set Source to Compound, give it a density, and attach components:

scintillator->Source.Set(static_cast<int>(Geant::MaterialSource::Compound));
scintillator->SetDensity(1.032);

auto *carbonAtoms = new Geant::MaterialComponent("carbon_atoms");
carbonAtoms->SetIngredient(...carbon...);
carbonAtoms->Measure.Set(static_cast<int>(Geant::ComponentMeasure::AtomCount));
carbonAtoms->AtomCount.Set(9);
Attach(scintillator, "Component", 0, carbonAtoms);

Each ingredient is its own object, because "how much" is a value someone edits and a port cannot carry one. It also means an ingredient that is itself a compound appears in the node graph with its own inputs.

An ingredient may be an Element or another Material. Atom counts only mean something for an element — a count of molecules of a mixture is not a quantity Geant4 has — and MaterialComponent::IsUsable() says so rather than letting G4 refuse it later.

Elements

auto *iron = new Geant::Element("Fe");        // natural isotopic mix

Symbol carries the element catalogue, the same way Name does. Turning FromNist off switches to explicit AtomicNumber and AtomicMass, which is for the case the table cannot express: an enriched or depleted sample, whose whole point is that the mix is not natural.

What a half-finished compound does

Nothing. It is the normal state while someone is typing, so the material keeps working and the note says what is missing:

with the fractions short, the note reads: mass fractions add to 0.8, not 1

GetFractionSum() is there so an editor can say this before the build, while the number can still be fixed. The two measures cannot be mixed in one material — Geant4 needs one or the other, and its own error for trying is about internal bookkeeping rather than about what you did, so that case is caught and named here too.

A failed build keeps the previous material. Dropping the geometry's material on every keystroke would be worse than holding a stale one until the new one is ready.

One thing to know about rebuilding

A G4Material cannot be edited once it exists — every change means a new one, and G4MaterialTable keeps them all for the life of the process. Since Rebuild() is reached from every property on the object, a signature of everything the material depends on is compared first, and an unchanged rebuild is skipped. Genuine changes still mint a new material (you will see names like my_scintillator_2), which is inherent to Geant4's design rather than something this layer can avoid.

Saved scenes

Material gained five fields, so it carries BOOST_CLASS_VERSION(..., 1) and reads them only when the payload has them:

if (version >= 1) { ar & HRPE("Source", ...); /* ... */ }

Note this is boost's per-class version, not kSceneDocumentFormatVersion. A class that defines serialize() is reached through ar & *self, which hands it the class version; only classes without one are given the document version through SerializeValues(). Reading the new fields unconditionally would take five values too many out of an older payload, and boost's binary primitives decode that misalignment as garbage rather than failing cleanly.


07 — Invalid geometry

Geant4 validates a solid's parameters in its constructor and reports a bad combination through G4Exception at fatal severity — and the default handler prints and calls abort(). Example 07 makes those mistakes deliberately.

ellipsoid->BottomCutZ.Set(100.0);
ellipsoid->TopCutZ.Set(-100.0);    // bottom above top: not a solid
ellipsoid->Update();

Three things happen, and the second matters most:

  • The build is refused, with a line naming the solid and G4's own message.
  • The previous shape is kept. An editor that dropped the geometry on every bad keystroke would blank the viewport while someone was still typing.
  • It recovers: the next sane value builds normally.

BuildGuarded is usable directly for anything that calls into Geant4 with numbers a user chose:

G4VSolid *impossible = Geant::BuildGuarded("hand_rolled", []() -> G4VSolid * {
  return new G4Tubs("bad_tube", 500.0, 100.0, ...);   // inner > outer
});
// null, and the program is still running

08 — Choosing the physics

"What physics list did you use?" is the first question asked of a Geant4 result. Physics answers it:

auto *physics = new Geant::Physics("physics");
physics->Reference.Set("FTFP_BERT");     // 36 available on this build
physics->EmOption.Set("_EMZ");           // 12 EM variants, chosen separately
world->SetPhysics(physics);

Both catalogues are read from G4PhysListFactory, so they are right for whatever Geant4 you are running. The EM half is a separate suffix because 36 × 12 combined names would help nobody.

A name this build does not have is refused, not substituted — a run reporting physics it did not use is a worse failure than one that will not start.

Two things worth knowing:

  • StepLimiter is on by default, and it matters: a UserLimits MaxStep is a request, and this is the process that honours it. Without it, capping the step silently does nothing.
  • The list is fixed once the run manager initialises. Choose before the first run; a later change is recorded but does not take effect.

A scene without a Physics keeps the built-in default, so existing scenes are unaffected.


09 — A dose map

A Scorer on a ScoringDetector answers "how much in that block". A ScoringMesh answers "how much here, and here, and here" — a grid of bins over space, independent of what volumes are there:

auto *mesh = new Geant::ScoringMesh("dose_map");
mesh->Size.Set(Vector3d(1000.0, 1000.0, 1000.0));   // half-lengths
mesh->Segments.Set(Vector3i(4, 4, 8));              // 128 bins
Attach(mesh, "Scorer", 0, edep);                    // an ordinary Scorer

A mesh takes no volumes — it covers whatever is inside it. That's the difference from a ScoringDetector, and it's why a dose map doesn't require cutting the geometry into as many pieces as you want resolution.

The scorer is an ordinary Scorer: the same nine quantities and the same filters. Only where it accumulates differs.

Reading back after a run:

std::map<std::int32_t, double> scores = mesh->ReadScores("edep");
double total = mesh->ReadTotal("edep");

Firing a muon pencil along z through a 4 × 4 × 8 mesh scores 8 bins of 128 — one column, which is the visible sign the mesh really is binning by position.

Two details worth knowing: a bin holds a G4StatDouble, so Geant4 tracks the weighted sum, its square and the count; ReadScores reports sum_wx(), the accumulated quantity, where mean() would divide by the event count and mean something else. And an empty map before a run is an answer, not a failure.

The bins are a picture

VoxImage<Voxel> image;
mesh->ToVoxImage("edep", image);      // dims, spacing and corner from the mesh
4 x 4 x 8 voxels of 500 x 500 x 250 mm
columns with anything in them: 4 of 16
on the beam axis: 9.10 MeV, off it: 0.00 MeV

Geant4 flattens a bin as i*nj*nk + j*nk + k, and BinOf(key) takes it apart again — which is what makes the map an image rather than a column of numbers. A VoxImage is what the rest of this library reconstructs into, so a simulated dose map and a measured one are now the same kind of object.

This did not work until recently, and the reason is worth knowing. The mesh registered plain primitive scorers, and a plain scorer indexes a hit by the copy number of the innermost replica — the z axis alone. Every bin sharing a z merged: a 4 × 4 × 8 mesh reported eight numbers instead of 128, and there was no honest way to map them onto anything. Scorer::MakeG4MeshScorer builds the 3D variants, which read all three replica depths. If you write a scorer of your own for a mesh, it has to be one of those.


10 — Track filtering

A 4 GeV muon crossing a metre of iron makes thousands of secondaries, nearly all of them low-energy electrons that stop within a millimetre. Following them costs most of the run. A TrackFilter is a say over every track before Geant4 follows it:

auto *filter = new Geant::TrackFilter("filter");
filter->PrimariesOnly.Set(true);      // or one of the milder settings below
world->SetTrackFilter(filter);

In the example, ten muons through iron:

hits recorded species seen
following everything ~15 800 5
primaries only ~1 200 1

The muons still cross, still deposit and are still recorded — what is gone is the shower around them. Whether that is acceptable is the physics question; the filter only makes it askable.

It is not a production cut

A cut (Region, Physics::DefaultCut) decides what is created. A filter decides what is followed once created. "Follow muons, ignore what they make" is a filter, and no cut can say it. They compose: a cut keeps the soft tail from existing at all, a filter drops what a cut cannot express.

The settings, in order of how much they change

Property Effect
MinSecondaryEnergy a floor in MeV; drops the soft tail only
DropNeutralSecondaries keeps ionising tracks, drops the gammas
ParticleFilter one PDG code and nothing else — primaries too
PrimariesOnly the bluntest; nothing the beam makes is followed

Two rules hold whatever is set. A primary is never dropped by an energy or charge rule — a beam somebody asked for is not noise, however feeble. And Enabled off beats everything, so a filter that is attached but switched off leaves a simulation exactly as it was, which is also the default state of a freshly created one.

ParticleFilter is the one rule that reaches a primary: a scene set to follow muons follows no electron, whoever launched it. Antiparticles are their own PDG code and are not implied — 13 does not admit -13.

The policy is a pure function

Classify(TrackFacts) takes four numbers — kinetic energy, PDG code, parent id, charge — and returns Track, Postpone or Kill. That is the whole policy, and it is public on purpose: every rule can be exercised without standing up a beam, which is what TrackFilterTest does before it fires anything. The Geant4 side only reads those numbers off a G4Track.

Attaching, changing and detaching between runs

Geant4 builds a worker's user actions once, when the thread starts, not per run. An action bound to whichever filter existed at that moment would freeze the choice for the life of the process. So the stacking action reads the scene's filter per track: a filter can be attached after a run has already happened, exchanged, or removed, and the next run honours it. An empty slot follows everything, which is exactly the behaviour of no stacking action at all.

GetKilledCount() and GetTrackedCount() are atomic and accumulate across workers and across runs; ResetCounts() clears them. They move nowhere but inside ClassifyNewTrack, so a non-zero killed count is proof Geant4 actually consulted the filter.


11 — A scintillator that emits light

Physics::Optical registers G4OpticalPhysics, and until materials could carry optical properties that bought nothing: optical physics acts on a material's G4MaterialPropertiesTable, and no material had one. A scintillator could be described as a lump of the right stuff, at the right density, that emitted no light.

An OpticalProperty is one entry of that table, attached to a Material through a port:

auto *rindex = new Geant::OpticalProperty("rindex");
rindex->Quantity.Set((int)Geant::OpticalQuantity::RefractiveIndex);
rindex->PhotonEnergies.Set({2.0, 2.8, 3.2});     // eV
rindex->Values.Set({1.57, 1.58, 1.60});
Attach(plastic, "Optical", 0, rindex);

Curves and constants

Geant4's table has two shapes and so does this. A curve is a value per photon energy — a refractive index, an absorption length, an emission spectrum. A constant is one number for the material — a light yield, a decay time. Which one a quantity is is a fact about the quantity, not a choice: IsCurve() answers it.

Quantity Key Shape Unit
RefractiveIndex RINDEX curve —
AbsorptionLength ABSLENGTH curve mm
RayleighLength RAYLEIGH curve mm
ScintillationSpectrum SCINTILLATIONCOMPONENT1 curve —
Reflectivity REFLECTIVITY curve —
Efficiency EFFICIENCY curve —
ScintillationYield SCINTILLATIONYIELD constant photons/MeV
ResolutionScale RESOLUTIONSCALE constant —
FastTimeConstant SCINTILLATIONTIMECONSTANT1 constant ns
SlowTimeConstant SCINTILLATIONTIMECONSTANT2 constant ns
Custom whatever CustomKey says either —

Energies are in eV throughout, because that is how optical photons are quoted — visible light is 1.7 to 3.3 eV, and a table written in MeV would be all zeroes to five decimal places.

Custom reaches anything the list does not name, including keys this Geant4 build has never heard of, which are created rather than refused.

What optical physics costs

Physics::Optical on is not free: a relativistic muon in a plastic radiator makes hundreds of Cerenkov photons per millimetre, and every one is a track. Five settings on Physics control that, and they reach G4OpticalParameters, which Geant4 reads at the start of each run — so unlike the physics list itself they can be changed between runs.

Property Effect
CerenkovMaxPhotons cap per step; zero leaves Geant4's own. The one that matters
Cerenkov generate Cerenkov photons at all — off times a run without them
Scintillation generate scintillation photons at all
TrackPhotonsFirst follow a step's photons before continuing its own track
ScintillationByParticle yield depends on what deposited the energy

What scintillation actually needs

A spectrum, a yield, a resolution scale and a time constant. All four, or Geant4's scintillation process has nothing to sample. A refractive index alone gives Cerenkov light and no scintillation at all.

Five 500 MeV muons through a 20 mm tile at 40 photons/MeV:

hits recorded: 3580
  pdg -22: 3509 (optical photons)
  pdg 11: 39 (delta electrons)
  pdg 13: 32 (the beam)

The optical photon's PDG code is -22 — not 22, which is a real photon, and not 0. It is a modelling construct with no code in the particle data tables, parked on the unused negative of a real photon.

Wrapping it

The bulk says how light is made and how far it gets. It says nothing about the edge of the volume, and the edge is most of the answer: a bare tile loses light through every face, the same tile in a reflective wrapper keeps it bouncing until it is absorbed or reaches a readout.

auto *wrapping = new Geant::OpticalSurface("wrapper");
wrapping->Kind.Set((int)Geant::SurfaceKind::DielectricMetal);
wrapping->Finish.Set((int)Geant::SurfaceFinish::Polished);
Attach(wrapping, "Property", 0, reflectivity);   // an OpticalProperty
Attach(wrapping, "Volume", 0, tile);             // a LogicalVolume
optical photon steps, bare:    3509
optical photon steps, wrapped: 33498

A surface carries its own table, built from ordinary OpticalProperty entries — Reflectivity is the usual one, Efficiency for a photocathode. Model chooses Unified (micro-facets, SigmaAlpha) or Glisur (one Polish number); Kind says what is on the far side; Finish the state of the surface and any paint.

Attaching a volume gives it a skin surface: the same treatment on every face. Geant4 keys skins by volume in a map, so a second one for the same volume would silently displace the first — editing the surface re-points the existing skin rather than making another.

Between two placements

A skin wraps a volume: every face, whatever is on the other side. A BorderSurface is the interface between two particular placements — glue between a crystal and a light guide, an air gap between two tiles:

auto *joint = new Geant::BorderSurface("tile_to_tile");
joint->SetSurface(glue);      // an OpticalSurface, referenced not owned
joint->SetFrom(first);        // PhysicalVolume
joint->SetTo(second);

It is directional. A photon crossing from From into To sees the surface; one going the other way sees nothing. That is Geant4's rule and not an oversight — a back-painted surface, a dichroic filter and a photocathode all behave differently depending on which side light arrives from. Symmetric on builds the reverse border as well, for the cases where the interface really is the same both ways.

The surface itself is referenced rather than owned, so one definition — "optical glue" — can be used at every border that has it. Geant4 keys borders by the ordered pair, so a rebuild re-points the existing one rather than making a second, the same rule the skins follow.

From and To being the same placement is refused: a volume's interface with itself is a skin surface, and saying so is more use than building something no photon can ever cross.

Birks quenching

plastic->BirksConstant.Set(0.126);    // mm/MeV, polystyrene

A scintillator does not emit light in proportion to energy deposited: a densely ionising track saturates the medium along its path. Without Birks' law a slow proton and a minimum-ionising muon depositing the same energy give the same light, and they do not. It sits on the material rather than on an optical entry, because it is a property of the medium's ionisation and applies whether or not anything optical is attached.

What is refused, and why

Two ways to get a curve wrong, and Geant4 refuses neither: energies out of order, which it interpolates wrongly rather than rejecting, and a value list that does not match the energy list. Both are refused here, with GetBuildError() saying which. A wrong refractive index is the kind of mistake that survives all the way into a result.

A material whose every entry is unusable gets no table at all rather than an empty one — optical physics reads an empty table as "transparent and non-scintillating" rather than as "unspecified", and a half-written entry is the normal state of something being edited.

One thing to know about ownership

G4Material owns the properties table it is given: handing it another deletes the first. A NIST material is also shared — G4NistManager hands the same G4_WATER to everybody — so an optical table set on one is visible from every Material object naming it, and the last writer wins. That is Geant4's model, not this library's.


12 — Digitizing a readout

A Hit is truth: this much energy, here, at this time, by this particle. No instrument reports that. A real channel integrates over a gate, sees a smeared total because its photostatistics are finite, reports nothing below its discriminator, and quantises what is left into counts. Comparing a simulation with data means comparing digis with digis.

auto *digitizer = new Geant::Digitizer("channels");
Attach(digitizer, "Detector", 0, readout);       // a SensitiveDetector
digitizer->Grouping.Set((int)Geant::DigiGrouping::PerVolume);
digitizer->EnergyResolution.Set(0.12);           // 12% at 1 MeV
digitizer->Threshold.Set(0.2);                   // MeV
digitizer->CountsPerMeV.Set(100.0);
digitizer->CountRange.Set(4095);                 // twelve bits
digitizer->Run();

Eight muons through three 20 mm tiles:

hits recorded: 62          (every step in the tiles)
digis: 24                  (three channels, eight events)
one channel: 3.593 MeV over 1 step -> 3.838 MeV = 384 counts at 5.967 ns

What a channel is

Grouping One channel per
PerCopy volume and copy number — a replica slice, a strip
PerVolume volume name, all copies together
Merged everything the detector saw

A hit records the physical volume's name, so three separate placements are three channels under PerVolume even when they share one LogicalVolume. Copy numbers tell replica slices apart, where the name is shared.

Events never merge. Whatever the grouping, hits from different events are different readings: hits accumulate across events, and summing two muons into one reading is energy no channel ever saw at once. This is what Hit::eventId is for.

The gate

IntegrationTime measures from the first hit of a reading, not from the previous hit — otherwise a slow drizzle would hold a channel open for ever. Four hits 8 ns apart under a 10 ns gate are two readings of two, not one of four. Zero puts a whole event in one reading.

A muon crossing a tile is over in 70 ps, so a gate changes nothing for a prompt run; it earns its keep when something arrives late, such as a delayed capture.

Resolution, threshold, counts

EnergyResolution is the fraction at 1 MeV and the spread scales as √E, which is photostatistics: 12% at 1 MeV is 3.8% at 10. Threshold is applied after smearing, because that is what a discriminator sees, and what falls below is counted in GetDroppedCount() rather than forgotten. CountsPerMeV quantises; CountRange saturates, counted in GetSaturatedCount().

Every digi carries both energy and trueEnergy (and time / trueTime): the difference between them is the instrument, and a study of resolution or of dynamic range needs both. A saturated reading stops at the top of the range while the truth beside it does not.

counts is -1 when no quantisation was asked for — which is not the same as a reading of zero counts.

Per event, as a real readout works

Run() after a run turns every hit held into digis, which is fine for a short study. PerEvent instead digitises at the end of each event, taking that event's hits from the detectors as it goes:

digitizer->PerEvent.Set(true);
world->SetDigitizer(digitizer);      // a Scene port, read per event

The hits are consumed — the digis become the record. That is the trade, and it is what keeps a long run from holding every step it ever took in memory. A study that wants both keeps PerEvent off.

It is safe with several workers, and the reason is worth stating: a worker asks only for the event it has just finished, so those hits are complete and nobody else is adding to them; TakeHitsForEvent runs under the detector's lock, and the digi list and counters are guarded too. DigitizerThreadedTest runs 40 events across 4 workers and requires exactly 40 digis with 40 distinct event ids and no hits left behind.

The Scene reads its digitizer per event rather than binding once, so one can be attached, swapped or detached between runs — the same arrangement TrackFilter uses, and for the same reason.

Joining digis to events

std::vector<Geant::Digi> readings = digitizer->GetDigisForEvent(id);

A digi carries the event that produced it, so the join is a lookup. There is deliberately no digi collection hanging off GeantEvent: that record does not persist its own muon path, so nothing would travel with it that does not travel already, and every consumer would pay for digis whether it wanted them or not.

Off is a pass-through

Enabled off gives one digi per hit, ungrouped, unsmeared, nothing dropped — the truth in the readout's own shape, so it can be compared with the instrument's version of the same run. It is not a silence.

Why it is not a G4VDigitizerModule

Geant4's digitizer machinery reads G4VHitsCollections out of G4HCofThisEvent. This library's SensitiveDetector keeps flat Hit records instead — which is what lets them serialise, travel over a socket and be reached from the object tree — so the digitisation runs on those. The gain is that Digitize(hits) is a pure function: every rule above is checked by handing it three hits, with no beam at all.


13 — Enriched materials

The NIST catalogue gives every element the isotopic mix a sample dug out of the ground has, and for almost everything that is right. Enriched uranium is the counter-example a muon tomography framework exists to look at, and inventing an element with the right average mass does not describe it: that gets the stopping power roughly right and the nuclear cross-sections wrong, because those belong to the nucleus rather than to an average of nuclei.

auto *u235 = new Geant::Isotope("u235");
u235->AtomicNumber.Set(92); u235->NucleonCount.Set(235);
u235->Abundance.Set(0.20);
Attach(element, "Isotope", 0, u235);      // and 0.80 of U-238 beside it
natural uranium: 238.029 g/mole over 3 isotopes
  U-234 at 0.00005, U-235 at 0.00720, U-238 at 0.99275
20% enriched:    237.450 g/mole over 2 isotopes

AtomicMass at zero leaves the mass to Geant4's own table, which is right for anything that is not exotic.

Isotopes beat FromNist, and are meant to: a mix somebody spelled out is not the natural one, by definition. Nothing has to be switched off first.

Two refusals, both because Geant4's own answer is worse than a refusal. A nucleus with fewer nucleons than protons is a fatal exception there, which is a poor reply to a number somebody is still typing. Abundances that do not add up to 1 are a warning during a build the caller never sees, after which Geant4 carries on with a normalisation nobody asked for — turning a typo into a quietly different material. GetAbundanceSum() and GetBuildError() say so while the number can still be fixed.

Editing an abundance after attaching reaches the element: the element listens to its isotopes, because a port reports attach and detach only.


14 — Dividing a volume

Example 04 places many copies with a ReplicaVolume: you size the daughter and Geant4 repeats it. Getting the arithmetic wrong — ten slices of 100 in a mother of 900 — is fatal, because Geant4 has no licence to change what you built.

A DivisionVolume is told the mother and how to cut it, and sizes the slices itself. The daughter's own dimensions along the divided axis are overwritten, which is the point: "cut this into ten" needs no arithmetic.

Mode Given Geant4 works out
ByCount how many how thick
ByWidth how thick how many fit
ByCountAndWidth both whether they fit

Offset leaves the first stretch undivided — the other thing a replica cannot say at all.

mother 1000 mm, slice as built 7 mm (and it does not matter)
ByCount, 10:                            10 slices
ByWidth, 125 mm:                         8 slices
ByCountAndWidth, 4 x 200 after 200 mm:   4 slices
5 x 300 in a block of 1000:              not built

Geant4 restricts what is divisible: the daughter's solid must be the same kind as the mother's, and not every kind can be cut along every axis. A refusal leaves IsBuilt() false rather than aborting, like every other bad geometry here — see 07.

One thing done differently from Geant4. Its width-only constructor cannot be used in 11.2: it reaches G4PVReplica::CheckAndSetParameters with zero replicas and aborts whatever the numbers. ByWidth therefore works the count out from the mother's bounding extent along the axis and uses the count-and-width form. ByWidth along Rho or Phi is refused, those not being lengths.


15 — Watching a run

Geant4 offers five places to hang behaviour on a run. The stepping action drives the existing muon reconstruction, 10 took the stacking one, and these are the other two.

RunSummary

auto *summary = new Geant::RunSummary("summary");
summary->ClearReadouts.Set(true);
Attach(summary, "Readout", 0, detector);    // a SensitiveDetector or Digitizer
world->SetRunSummary(summary);
runs finished: 8
events in the last one: 6
it took 0.002 s (3135 events/s)

The numbers are the last run's, not a running total: how big was that run is the question a summary answers. They are read on the master thread, where Geant4 has merged the workers' events — a worker's own G4Run knows only what that worker did, and four workers each reporting their share would be four wrong answers rather than one right one.

ClearReadouts empties the attached detectors and digitizers when a run begins. Off by default, because it changes what a second run means and something may be relying on hits accumulating; on is what an interactive session wants, where pressing run twice should give the second run's answer.

Geant4 allows one run action per thread, so this shares the one that re-clones the emitter between runs.

TrackRecorder

One record per track rather than one per step:

auto *recorder = new Geant::TrackRecorder("tracks");
world->SetTrackRecorder(recorder);          // PrimariesOnly is on by default

A TrackRecord holds where the particle started and stopped, the direction at each end, the energy at each end, the path length, the time and the step count — and ScatteringAngle(), the angle between the two directions.

That last one is close to the whole measurement a scattering tomograph makes. Firing 3 GeV muons half at a 300 mm lead block and half at air beside it:

through the lead: 24 muons, mean scattering 49.1 mrad
through the air:  24 muons, mean scattering  0.7 mrad

Two directions per track and a place they crossed is the measurement; the step-by-step record is thousands of numbers saying the same thing.

PrimariesOnly is on by default here, unlike everywhere else, and for a reason worth knowing: a muon through iron makes thousands of secondaries, and recording a line for each costs more than the tracking does. ParticleFilter, MinInitialEnergy and MaxTracks narrow it further; what the filters drop is counted in GetRejectedCount() rather than forgotten.

Both objects hang off Scene ports read per run, so either can be attached, swapped or detached between runs — with one caveat: the tracking action is bound to its recorder when Geant4 builds a worker's actions, so detaching a recorder mid-process does not stop an already-built worker recording. The run summary and the track filter read their slot per call and have no such limit.


16 — A voxel phantom

Everywhere else here a VoxImage is a result: what a reconstruction produces, and what 09 converts a scoring mesh into. This is the other direction — an image becomes geometry Geant4 tracks through, voxel by voxel, with each voxel's material decided by its value.

auto *phantom = new Geant::VoxelPhantom("scan");
phantom->SetImage(image);
Attach(phantom, "Material", 0, air);      // band 0: the bottom of the range
Attach(phantom, "Material", 1, lead);     // band 1: the top
container->GetChildren()->AddObject(phantom);   // nested, like a replica

Three levels, not a million placements

A million placed boxes would leave Geant4 navigating them one by one. A nested parameterisation is three levels of repetition — slices in x, columns in y, voxels in z — where only the innermost is parameterised and a voxel's material is computed from the copy numbers of all three. Navigation becomes arithmetic, so an 8 × 8 × 8 phantom costs three divisions rather than 512 placements, and a 256³ one costs the same three.

What a value means

A voxel's value becomes a material by falling into one of the bands between ValueMin and ValueMax — one band per attached Material, in order. Values outside the range clamp to the end bands rather than wrapping: an image is not obliged to stay inside what somebody typed.

That is the honest general answer. What a number in an image means is a fact about where the image came from, not about this class.

The container has to fit

GetRequiredSize() is the image's dimensions times its spacing, and the mother box has to be exactly that. A mother of the wrong size is refused, not stretched: a phantom that quietly does not match its image is a wrong answer nobody would think to look for. MaxVoxels guards against turning an image nobody meant into geometry.

Closing the loop

Eight muon columns across a phantom with a 2 × 2 lead block drawn through its middle:

x column 0: 0.7 mrad
x column 1: 0.5 mrad
x column 2: 0.6 mrad
x column 3: 81.9 mrad ########################################
x column 4: 84.8 mrad ##########################################
x column 5: 0.4 mrad
x column 6: 0.5 mrad
x column 7: 0.5 mrad

The image became geometry, the muons crossed it, and the measurement put the block back where it was drawn.


17 — A plane probe

A sensitive detector needs a volume, which means building geometry for something that is not there — a slab of air the muons have to be told not to scatter off. A PlaneProbe is the other thing you often want: a surface with no substance that says what went through it.

PlaneProbe --(EventTable)--> Recipe --(Histogram)--> a plot

It is a TRS, so it is placed and turned exactly like geometry — drag it in the viewport, type coordinates in the panel, or nest it under something that moves. Its plane is its own local xy at z = 0, bounded by Width and Height, so rotate it 90° about x means what it looks like it means.

It is found by being in the tree. Nest one under a Scene and every step of every track is offered to it; there is no port to wire and no volume to build, the same as for an emitter or a placement.

What it reports

Crossings go into an EventTable on the probe's Events output port: rows of doubles with named columns, so a recipe asks for "momentum" rather than counting along a row.

Report Columns
CrossingPoint event, track, pdg, x, y, z, u, v, dx, dy, dz, momentum, energy, time
WholeTrack those, plus startx/y/z, endx/y/z, initialEnergy, finalEnergy, length

u and v are the crossing in the probe's own frame, which is usually what you want for a hit map; x, y, z is the same point in the world.

To plot two columns against each other — u against v for where they landed, u against momentum for whether the fast ones landed anywhere different — the recipe fills a PointSet instead of a Histogram, and a RootScatter draws it. One line: Crossings.SetPoints(Events.Column("u"), Events.Column("v")). See Histograms and plots.

Everything is interpolated to the crossing itself rather than taken from the start of the step: for a metre-long step in air the pre-step values describe a muon a metre away, which is not the same muon.

ParticleFilter defaults to 13 and PrimariesOnly to on, because a probe reporting every delta electron as well would bury the muons a hundred to one.

What the recipe decides

The probe does no arithmetic beyond finding the crossing. What the numbers mean is the recipe's business:

Events = InputPort("Events")
Spectrum = OutputPort("Spectrum")

def run(ctx):
    Spectrum.Reset()      # optional: ResetOnRefill already does this
    for i in range(Events.GetRowCount()):
        Spectrum.Fill(Events.Value(i, "momentum"))

Wire probe.Events → recipe.Events and recipe.Spectrum → a Histogram, and that is the chain. Histograms and plots is the other end of it: what a Histogram holds, and how to get it onto a canvas. Wanting the scattering angle instead, or the momentum only for tracks arriving off-axis, is an edit to the script — no object changes. EventTable, Histogram and PointSet are all bound for Python, which is what makes those three lines work — and what lets the same recipe fill a PointSet instead, for a scatter rather than a spectrum.

Sixty 3 GeV muons through 600 mm of iron, as the example runs it:

crossings recorded: 60 of 43313 steps offered
momentum spectrum at the plane:
  2150 MeV/c: ### 6
  2250 MeV/c: ############## 22
  2350 MeV/c: #################### 31

Moving it measures again, live

The rows in the table were measured with the plane where it was. Changing anything a crossing depends on makes them out of date:

  • its placement (position, rotation, scale), from the panel or from a drag in the viewport;
  • Width or Height;
  • Enabled, ParticleFilter, PrimariesOnly or MinMomentum;
  • Report.

A rename, a display setting, or a property re-announced with the same value changes nothing. Only the probe's own placement counts, so moving a parent it is nested under does not.

A moved probe replays the last run instead of running again. A probe has no substance, so where it stands changes nothing about what the particles did. During a run, the Scene keeps the steps its probes could ask about. A probe that has changed offers those same steps to its new plane and replaces its table in one go. The result is exactly the table a probe that had been there all along would have recorded: same events, same rows. It takes milliseconds, where a run takes seconds.

  • What is kept: only what the probes present at run time could want. That means primaries only if they all watch primaries only, and one particle if they all watch the same one. For default probes (primary muons) that is a few dozen steps per event.
  • Limit: at most Scene::SetProbeStepLimit() steps (2 million by default, about 200 bytes each). A run over the limit keeps nothing.
  • When a replay is not possible: the probe is changed to want steps that were not kept (every particle, say), the limit was exceeded, or there has been no run yet. The table is then emptied and the probe emits UpdateRequested(). The Scene it is nested in answers by running again, if AutoRun is set on the scene. That run keeps enough for the new settings.

A replay emits Updated() on the table once. A recipe reading that table through an input port, with AutoRun set on the recipe, runs on it and refills its histogram, and a ROOT canvas showing the histogram redraws. In gcompose, dragging the probe therefore gives a live pipeline:

drag ─▶ probe replays ─▶ table ─▶ recipe (AutoRun) ─▶ histogram ─▶ canvas

The replay and the recipe are each throttled to the AutoRun interval (see AutoRun and live updates). The plane itself redraws on every mouse move, and the rest follows a few times a second.

To set it up, starting from example 17:

  1. Run the scene once, so there is something to replay.
  2. Set AutoRun on the recipe.
  3. Switch the histogram's eye on in a ROOT canvas pane.
  4. Drag the probe.

Each run of the recipe replaces the histogram rather than adding to it. The histogram's ResetOnRefill (on by default) empties it when the recipe starts, whether or not the script calls Spectrum.Reset(). Turn it off only to accumulate over runs on purpose.

A Scene run empties all its probes together with its event collection, so a table always holds exactly one run.

One thing that had to be got right

A step that runs along a rotated probe's plane comes out of the single-precision world matrix with its ends at z = +1e-7 and −1e-7. Compared against zero those are opposite sides, and the probe would report a crossing of a plane the particle never went through. The sign test uses a tolerance that scales with the step, because the noise scales with the coordinates.


19 — Tomography from a run that keeps nothing

The usual muon-tomography chain remembers every muon it measures:

Scene --(ScatterCollection)--> MuonScatterCollection --(Events)--> PocaImageAlgorithm

Ten million events is ten million MuonScatter values — around 1.4 GB — held so that an algorithm can read them back once and drop them. Moving one wire skips the remembering entirely:

Scene --(ScatterSink)--> PocaImageAlgorithm --(Image)--> MultiVoxImage

ScatterSink takes anything implementing IMuonScatterSink, which every MuonTomographyAlgorithm does. The scene opens a session when the run starts, hands over each muon at the end of the Geant4 event that measured it, and closes the session when the run ends. The muon is folded into its voxel and forgotten: the run holds one at a time, whatever its length, and the image fills while the beam is still going.

Why a push and not a stream port. A DataBuffer and the algorithm's Stream port already move items one at a time, but they are pull: the consumer asks, the producer waits. A Geant4 run cannot wait in the middle of BeamOn — it is the one driving — so a pull seam needs a thread and a queue between them, and a queue of muons is the collection you were avoiding. Here the producer calls and the consumer folds.

What it costs. No second pass. A collection can be re-read with different cuts, saved, or handed to another algorithm later; a streamed muon cannot, because after it lands in the image it is nowhere else. Both ports can be wired at once, which collects and streams — for a short run worth keeping as well as reconstructing.

LiveUpdateEvery on the algorithm decides how often the image redraws while being fed — every 200 muons by default, zero for once at the end. It changes nothing about the result, only whether a long run shows itself filling.

Muons measured at one end only are streamed too, exactly as they are collected, and a consumer counts them as skipped. That number is worth watching: a run where most muons are one-ended is a lower chamber in the wrong place.


20 — A stream that is waited on

Example 19 pushes: press Run on the scene, and the muons land in the image as they are measured. This is the same three objects with control the other way round, and a stream between them:

Scene --(ScatterSource)--> MuonScatterSource --(Stream)--> PocaImageAlgorithm

Press Run on the algorithm. It reads its stream port, finds nothing, and that ask is what starts the run. Muons arrive while it waits, a batch at a time; each is folded into the image and dropped; when the run reaches NumberOfEvents the stream reports itself finished and the algorithm stops.

Two kinds of producer behind one port

pull push
what it is a collection, a decoded file, anything iterable a detector reading out, a simulation while it runs
asking asking is answering nothing is there to be asked for
an empty answer the end of the stream a moment before the next item
an idle producer there is no such thing the ask starts it

A consumer does not choose between them, and does not need configuring for either: it asks the producer which it is — IStreamWaitProducer::CanPullStream() — and behaves accordingly. So the same PocaImageAlgorithm reads a MuonScatterCollection, a DataBuffer or a live run without knowing the difference, and a stream port means the same thing in every graph it appears in.

MuonScatterSource is the push side. It is a sink to the run and a stream to the algorithm, so neither half knows the other exists, and it holds Capacity muons at most — the whole memory cost of the pipeline, whatever the length of the run. BatchSize is how many one ask is worth: asked for one muon at a time a simulation spends its time starting runs; asked for too many it works far ahead of what is being read. WaitTimeout bounds a single ask, so a consumer waiting on a producer that has gone quiet still notices a Pause or a Stop.

Why the run happens in chunks and not on a thread

Geant4's geometry and run manager belong to the thread that built them, which is why Scene::Start() runs on the caller's thread rather than a background one. So an ask for muons runs a chunk of the simulation — BatchSize events — right there on the asking thread, and returns. The consumer sees muons long before the run is over, and Geant4 stays where it needs to be. A producer that can run asynchronously is free to do so instead: the source waits on it either way, and the consumer cannot tell which it got.

The run ends, and stays ended

A stream that refilled itself whenever it was emptied would have no end, and a consumer reading to the end would never return — which is exactly the loop the first version of this deadlocked in. So a driven run that reaches its event count reports itself finished and goes on reporting it. Scene::ResetStreamProduction() is what asks for another pass; wiring a source to the port calls it for you, so a freshly connected pair always has a run to do.

Note: NotifyEvery vs LiveUpdateEvery on Streamed Events

When streaming events into a reconstruction algorithm (such as PocaImageAlgorithm), two settings manage intermediate cadence and updates, but at fundamentally different architectural levels:

[!NOTE] Key Difference: - NotifyEvery (on MuonScatterSource / StreamOutput): An execution trigger for push/chained pipelines. When the buffer queue reaches \(N\) items, the stream announces itself (Updated()), causing a downstream algorithm with AutoRun enabled to launch a new run to drain and fold the batch into the image. - LiveUpdateEvery (on PocaImageAlgorithm / MuonTomographyAlgorithm): A display redraw throttle within a single ongoing run. Every \(N\) events accumulated, the algorithm emits image->Updated() to refresh 3D viewports and slice views so the UI doesn't look frozen during a long run.

Setting MuonScatterSource.NotifyEvery PocaImageAlgorithm.LiveUpdateEvery
Defined On Stream producer/buffer (MuonScatterSource) Reconstruction algorithm (MuonTomographyAlgorithm)
Pipeline Role Execution trigger (Push / Chained) Display refresh cadence (Internal Redraw)
Trigger Condition Queue length reaches \(N\) waiting items \(N\) events processed within the current algorithm run
Action Taken Calls MuonScatterSource::Updated(), notifying downstream graph Emits image->Updated() and this->Updated()
Effect on Run Triggers a new algorithm execution (if AutoRun is on) Refreshes the viewport mid-run without stopping execution
Active Mode Run simulation (Scene); algorithm has AutoRun=true, DriveStream=false, ResetImage=false Run algorithm directly (or feed it as a direct sink via ScatterSink)

The two ways to observe progress

  1. Chained / Push-driven (NotifyEvery): You start the run on the simulation (Scene). The simulation pushes muons into MuonScatterSource. Every NotifyEvery events (e.g., 20), the source announces a batch. The reconstruction algorithm (AutoRun=true, DriveStream=false, ResetImage=false) wakes up, drains the available queue, folds those events into the existing image without wiping it, and returns. The image grows in discrete steps as the beam continues.
  2. Interactive Redraw (LiveUpdateEvery): Used during a single long-running accumulation (such as direct sink streaming via world->SetScatterSink(poca) in Example 19). The algorithm checks LiveUpdateEvery during AcceptMuonScatter() and emits Updated() so VTK viewports redraw intermediate progress without waiting for all events to finish. (Note: when driving a stream via poca->Run() as in 20_driven_stream.gcx, the consumer pulls from StreamInput::Read(), which by default finalizes the image at the end of the full stream).

21 — Muons and the scene in one TFile

Examples 19 and 20 keep nothing: the muons become an image and are gone. That is right until somebody asks to re-run the reconstruction with a different cut. This is the third destination:

Scene --(ScatterSink)--> RootScatterWriter --> muons.root
                                                    |
                         RootScatterReader <--------+
                                 |
                           (Stream) --> PocaImageAlgorithm --> image

RootScatterWriter is a sink like any other, so the run cannot tell it is writing a file rather than filling an image, and holds one muon either way. RootScatterReader is a pull stream — the muons exist already — so the reconstruction reads the file through the port it reads a live run through, unchanged.

What is in the file

root [0] TFile f("21_muons.root"); f.ls()
  KEY: TTree                  muons              uLib muon scatter events
  KEY: ROOT::Mutom::SceneBlob uLibScene          the scene, header plus bytes
  KEY: TTree                  uLibSceneContents  what that scene holds
root [1] muons->Draw("scatterAngle")

The muon tree holds ROOT::Mutom::MuonScatter (src/Root/RootMuonScatter.h): plain floats, split into leaves, with the PoCA and the scattering angle worked out once and stored beside the tracks. ROOT reads it from the StreamerInfo in the file, so an analysis needs neither uLib nor our dictionary — it tries to load libmutomRootDict first and carries on happily when it cannot.

The result beside the measurement

A muon needed a class of its own to live in a TFile, because ROOT has no idea what a muon is. A histogram does not: ROOT has had TH1 since before uLib existed. So the data classes are converted rather than mirrored (src/Root/RootMathIO.h):

uLib ROOT
Histogram (1, 2 or 3 axes) TH1D / TH2D / TH3D
PointSet TGraph, or TGraphErrors when it has errors
VoxImage / MultiVoxImage TH3D, one per channel, in world coordinates

WriteContextDataToRootFile() walks a context and writes every one it finds, reporting both the keys it wrote and anything it had to skip. So the example's file ends up holding the measurement, the result, a plot of it, and the scene that produced all three:

--- everything in the file ---
  muons                  TTree
  poca_image             TH3D
  scatter_angles         TH1D
  uLibScene              ROOT::Mutom::SceneBlob
  uLibSceneContents      TTree
--- image: 20 x 20 x 20 voxels, 149 total, x in [-500, 500] mm ---
--- angles: 150 entries, mean 0.02247 rad ---

An image as a TH3D is the one worth arguing about: a TH3 is a histogram and an image is not. It wins anyway, because it carries what an image needs — three axes with real extents and a value per cell — and poca_image->Draw("iso") works in a bare ROOT session. The axes are the grid's world coordinates, so the image lands where the scene says it is; the conversion back restores dims, spacing and position, and the test checks the position rather than only the values, because an image half a box away from where it belongs is the failure that looks like success.

A scene in a ROOT file, like a scene in a .gcx

WriteContextToRootFile() puts a whole context into the same file (src/Root/RootSceneFile.h). A scene holds arbitrary uLib classes — whatever was dropped into it, including classes a plugin only adds at run time — so there can be no per-class ROOT dictionary written in advance. What goes in is the same self-describing document a .gcx carries, as bytes, in a dictionary-backed SceneBlob whose header fields (application, date, format version, object count) are readable without decoding anything.

Beside it goes uLibSceneContents: one entry per object, with its class, name, parent and payload size. Nothing restores from that tree — the bytes are the truth — but it is what lets TBrowser or tree->Scan() answer "what is in this scene" to someone who has ROOT and nothing else.

ReadContextFromRootFile() is the other direction, and it goes through the same RestoreContext a .gcx does, so a scene saved either way comes back the same.

One caution

Two TFiles on one path do not coexist: a reader holding the file open READ while the scene is written UPDATE ends up pointing at a tree that no longer exists, and the symptom is an entry count that quietly becomes zero. The example closes the reader before writing the scene, and says so where it does it.


Working in gcompose

Nothing gcompose shows about an object is declared in gcompose. The Add menu comes from the factory category, the panel's sections from the property groups, the tooltips from the descriptions, the dropdowns from the enum labels, and the dots in the node graph from the port directions. All of it lives on the object, and GuiMetadataTest checks that none of it is missing.

Finding a class

The Geant4 menu has sections, one per thing you might be trying to do:

Section What is in it
(top level) Scene, Physics — what a simulation needs before the rest makes sense
Solids thirty shapes, from BoxSolid to ExtrudedSolid
Boolean unions, intersections, subtractions, scaled solids
Volumes logical and physical volumes, replicas, divisions, assemblies, regions, parallel worlds, the voxel phantom
Materials materials, elements, isotopes, compound components
Optical optical properties, skin and border surfaces
Readout sensitive detectors, scorers, the scoring mesh, the digitizer, the three filters
Fields magnetic and electric fields, user limits
Actions the track filter, the run summary, the track recorder
Sources the particle gun, the isotropic source, the three cosmic emitters
Results events and event collections

Reading an object

Two property groups say what happened rather than what to do, and both are read-only and never saved — a run's outcome is not part of the scene:

  • Status — did it build, and why not. Material.Problem says the fractions add to 0.8; VoxelPhantom.Container Size Needed says what the box has to be; BorderSurface.Problem says which end is still unwired. A half-finished object is the normal state of one being edited, and this is where it says so.
  • Results — what a run produced. TrackFilter.Killed, Digitizer.Below Threshold, RunSummary.Rate, SensitiveDetector.Hits. A zero here after a run usually means the object was never reached, which is a different problem from one that reported nothing.

Every property has a tooltip. Some of them are worth reading twice: a solid's Height is the full height where Geant4 takes a half-length, and every angle is in radians.

Wiring one up

Reference ports carry a dot in the node graph, and the ones on Scene are where a simulation is assembled: Physics, Emitter, EventCollection, TrackFilter, Digitizer, RunSummary, TrackRecorder. A LogicalVolume's Solid and Material are ports too — a logical volume is a shape plus what it is made of, and joining those two is the gesture the graph exists for.

A ScoringMesh has an output: its Image port hands out the dose map as a VoxImage, refreshed by Update Image, so a mesh can be wired straight into an image filter or a viewer. ImageScorer picks which scorer feeds it; empty takes the only one there is.

Two references deliberately have no dot. Container and Representation hold a component the object builds for itself, and pointing them elsewhere would destroy it.

Seeing what came out

A Histogram is drawn by a ROOT-canvas pane, and nothing is drawn until the eye next to it in the context tree is switched on — a context full of histograms would otherwise fill every canvas in the window. Canvas Rows and Canvas Columns in the pane's settings say how many plots fit side by side. The axis labels and log scales are display properties of the plot, not of the histogram.

A PointSet is drawn the same way, as a RootScatter. Two spectra that answer the same question belong on the same axes rather than in two pads: nest both plots under a RootPlotGroup and switch the group's eye on, and each member keeps an eye of its own for leaving one out. A scatter and a spectrum cannot share axes and the group says so rather than drawing them over each other. See Histograms and plots.


Surviving a fatal G4Exception

Geant4 reports a bad argument by calling G4Exception, and the default handler for a fatal severity prints and calls abort(). That is defensible for a batch job whose input was built once and correctly. It is not survivable in an editor, where every intermediate state of a drag is a combination the user has not finished specifying.

Geant/Guard.h is the one place that turns any of them into something a caller can catch:

G4VSolid *built = BuildGuarded("my_trap", [&] { return new G4Trap(...); });
if (!built) { /* the previous shape is still valid and on screen */ }

RunGuarded(name.c_str(), [&] { filter->add(name); });

The guard is scoped, not installed globally, and restores the previous handler on the way out including when the throw happens. Outside construction — during a run, in the physics — a fatal G4Exception means something this framework has no recovery for, and converting it to a throw would only move the crash somewhere less informative.

Known offenders: G4Trap with non-planar side faces, G4Ellipsoid cut below its own bottom, G4SDParticleFilter given a particle the physics list never built, and any placement violating the only-daughter rule.


Reference

Header Objects
Geant/Matter.h Element, MaterialComponent, Material, the NIST catalogues
Geant/Solid.h 30+ solids, LogicalVolume, PhysicalVolume
Geant/Volumes.h ReplicaVolume, ParameterisedVolume, AssemblyVolume, Region
Geant/AttachmentPorts.h AttachmentPorts -- the growing port group every object above uses
Geant/Field.h UserLimits, MagneticField, ElectricField, CarrierRecord
Geant/Sensitive.h Hit, SensitiveDetector
Geant/Scoring.h ScoringDetector, Scorer, ParticleFilter, ChargeFilter, EnergyFilter
Geant/ScoringMesh.h ScoringMesh — scoring on a grid over space
Geant/Physics.h Physics — reference lists and EM variants
Geant/Sources.h ParticleGun, IsotropicSource
Geant/Parallel.h ParallelWorld, ReflectedVolume
Geant/Guard.h BuildGuarded, RunGuarded, ExceptionScope

Tests live in src/HEP/Geant/testing/: NewSolidsTest, VolumesTest, FieldTest, SensitiveTest, ScoringTest, ParallelTest, CompoundMaterialTest, RegionSerializeTest.