Implement G4Trap and G4Trd support in geometries through convexpolyhedron by ggalgoczi · Pull Request #331 · BNLNPPS/simphony

ggalgoczi · 2026-05-15T17:14:14Z

Summary

Adds support for G4Trap and G4Trd solids by converting them to CSG_CONVEXPOLYHEDRON . 8 vertices computed from G4 parameters at conversion time, 6 outward face planes, routed through the existing GPU intersection kernel.

Also adds validation scripts (tests/g4trap_validation.sh) covering both solids with isotropic torch sources and Cherenkov+Scintillation physics from 5 GeV electrons.

Changes

Geometry support (CSG + sysrap + u4)

sysrap/sn.h — ConvexPolyhedron factory, planes field on sn, getPlanes() accessor, CSG_CONVEXPOLYHEDRON case in setAABB_LeafFrame
CSG/CSGImport.cc — plane pass-through to foundry; planeIdx is 0-based (matches csg_intersect_leaf_convexpolyhedron.h kernel); bypasses
ExpectExternalBBox assertion when planes are present
u4/U4Solid.h — _G4Trap / _G4Trd enum, init_Trap() (handles all 11 G4Trap params including alpha shear + symAxis tilt), init_Trd()
(5-param degenerate case), shared _setRoot_FromVertices() helper

G4 hit writer for GPURaytrace (`src/GPURaytrace.h`)

Adds per-event G4 hit writing to g4_hits_output.txt (same format as simphony's writer, for diff-able comparison)
Mutex-guarded with G4AutoLock on a dedicated g4_hits_file_mutex — safe under MT G4 (G4TaskRunManager default)

Test infrastructure (`tests/`)

g4trap_validation.sh + g4trap_validation.README.md — driver script with per-test PASS/FAIL based on count rel-diff and per-axis χ² .
Auto-installs torch configs to share/simphony/config/ and uses macros that already set /process optical/boundary/setInvokeSD true
geom/test_trap.gdml, geom/test_trd.gdml — iso-source geometries (Quartz, flat n=1.46, ABSLENGTH=100m)
geom/test_trap_scint.gdml, geom/test_trd_scint.gdml — full-physics geometries ( dispersive RINDEX + synthetic scintillation, yield=10 γ/MeV, τ=2 ns, ABSLENGTH=100m).
run_validate.mac, run_validate_5evt_1t.mac — G4 macros

Validation results

All 4 tests pass at seed=42, OPTICKS_PROPAGATE_EPSILON=0.001, OPTICKS_MAX_BOUNCE=10000:

Test	Solid	G4 hits	Simphony hits	rel diff	max pos+dir χ²/ndf
`trap` iso	G4Trap	11083	11087	0.04%	< 0.2
`trd` iso	G4Trd	7687	7677	0.13%	< 0.2
`scintillation_trap`	G4Trap	472693	472017	0.14%	< 1.5
`scintillation_trd`	G4Trd	892626	893156	0.06%	< 1.1

Implements G4Trap by computing the 8 vertices and 6 face planes from the 11 G4Trap parameters at conversion time, then routing through the existing CSG_CONVEXPOLYHEDRON GPU intersection (no new device code). Adds plane-array plumbing to the sn -> CSGNode pipeline that was previously only stubbed: - sysrap/sn.h: std::vector<double>* planes field on sn (nullable), sn::ConvexPolyhedron factory, getPlanes() accessor, CSG_CONVEXPOLYHEDRON case in setAABB_LeafFrame - CSG/CSGImport.cc: importNode now passes node planes through to CSGFoundry::addPlan and sets planeIdx/planeNum on the CSGNode; bypasses the ExpectExternalBBox assertion when planes are present - u4/U4Solid.h: _G4Trap enum, Type/Tag/Constituents, init_Trap() (computes 8 verts with theta/phi/alpha shear, then 6 face planes via cross product on vertex triples) - tests/geom/test_trap.gdml: minimal test geometry with a trap solid - tests/visualize_trap.py: round-trip plane->vertex visualization

CSGImport::importNode was setting planeIdx = 1 + fd->plan.size() (1-based), but the GPU intersection at csg_intersect_leaf_convexpolyhedron.h reads plan[planeIdx + i] directly (0-based), as does CSGFoundry::addNode at line 1839 which uses setPlaneIdx(plan.size()) before adding planes.

The dy*alpha contribution to vertex x-coordinates had inverted sign compared to G4's G4Trap::SetAllParameters convention. This mirrored the trapezoid in x, causing an 80mm hit displacement on DIRC prism geometries. Verified against G4 on y_pixel_detector.gdml: GPU hit centroid now matches G4 to 0.01mm.

Extract the 6-face plane + AABB computation from init_Trap into a private _setRoot_FromVertices(v[8][3]) helper, then add init_Trd that fills its 8 vertices and reuses the same path. G4Trd is the {alpha=theta=phi=0, dx1=dx2, dx3=dx4} degenerate case of G4Trap, very common in calorimeter geometries. Trap output is bit-identical to before the refactor (6 planes matching the analytic prism with alpha=-15deg on test_trap.gdml). Trd output on a new test_trd.gdml (x1=200,x2=60,y1=300,y2=80,z=400) gives the expected outward-normal face set tilted toward +z because the +z face is narrower than the -z face.

…rison The default G4 behaviour is to NOT invoke the sensitive detector when a photon hits an EFFICIENCY>0 surface (the boundary process kills the photon at the interface, before the SD volume step would otherwise fire). Setting /process/optical/boundary/setInvokeSD true makes G4 call the SD on Detection, so trap/trd hit counts can be compared directly between G4 and Opticks. Without this flag, G4 records 0 hits while Opticks records ~10000.

Same line format as the existing opticks_hits_output.txt, so the two files can be diffed directly for G4-vs-GPU validation. Previously GPURaytrace only counted G4 hits and wrote nothing about them, so any distributional comparison required adding a writer. Caveat: ofstream open+append is NOT thread-safe. When G4 runs MT (the default), worker threads race on the file and ~12% of lines are lost. The hit count from the master-thread accumulator (fTotalG4Hits) is correct under MT; only the file contents are racy. For clean per-photon output use the new tests/run_validate_10evt_1t.mac which sets /run/numberOfThreads 1. A G4AutoLock-guarded writer would be a proper fix but is out of scope for this branch.

Adds the artefacts needed for a reviewer to reproduce the branch's validation: GDMLs (tests/geom/): - test_trap_side.gdml — detector at +X to probe slanted-wall refraction - test_trap_absorbing.gdml — finite ABSLENGTH eliminates bounce-budget effect - test_trap_dispersive.gdml — Quartz Sellmeier RINDEX, used by Cherenkov test - test_trap_scint.gdml — Quartz with synthetic scintillation properties Macros (tests/): - run_validate_10evt.mac — beamOn 10 with boundary InvokeSD on - run_validate_10evt_1t.mac — same, single-thread to avoid the G4 writer race Script: - g4trap_validation.sh + README — runs all 7 tests, prints PASS/FAIL per case, applies pre-existing torch source configs if missing in share/, uses separate tolerances for torch sources (5%) vs cherenkov (5%/8) vs scint (10%/50) reflecting the natural RNG widths. Reviewer entry point: ./tests/g4trap_validation.sh Per-test: ./tests/g4trap_validation.sh {trap|trd|single_photon|cherenkov|scintillation}

…ent macro 5 GeV electron in 40 cm Quartz at SCINTILLATIONYIELD=1000 produces ~170k detector hits per event from scintillation alone, on top of Cherenkov. G4 single-thread propagation of that many photons is the bottleneck (several minutes per event). 100 photons/MeV cuts to ~17k hits/event, keeping enough statistics for distributional comparison without making the validation script too slow for reviewers to run. run_validate_5evt_1t.mac is also dropped to 5 events for the same reason.

… 5-event macro

The prior writer used 'eventID==0 ? trunc : app' which races under MT — multiple worker threads each saw eventID==0 for their first event and each truncated the file. ~12% of hit lines were lost in 36-thread runs. Switch to a process-wide 'static bool first_event' under a dedicated G4AutoLock, so the file is truncated EXACTLY ONCE across all threads, then appended to (atomically) for every subsequent event. Hit count and hit file now agree under any thread count.

…utex-guarded

…ff, chi2/ndf<5)

Opticks's U4Scint::Create at u4/U4Scint.h:97 requires the OLD G4 10.x property names FASTCOMPONENT, SLOWCOMPONENT, REEMISSIONPROB to recognize a material as a scintillator. With only the new G4 11.x names (SCINTILLATIONCOMPONENT1 etc.) U4Scint returns nullptr, the GPU wavelength ICDF is never built, and scintillation gensteps get sampled with invalid wavelengths -- photons are silently lost during propagation. Symptom before fix: with Cerenkov enabled, all 534664 Opticks hits were tagged CreationProcessID=0 (Cerenkov), zero scintillation hits, while G4 produced ~5k extra hits in the 350-450 nm scintillation peak. With Cerenkov disabled, Opticks generated 1982 photons but recorded 0 hits. Fix: declare both old and new property names pointing to the same matrices. G4 reads the new names, U4Scint reads the old ones.

…=100m eps-scan on the scintillation test geometry found eps=0.001 (1 µm boundary stepping) is the sweet spot vs G4: eps=1e-5 -> 467k hits (count diff +1.17%, all chi2/ndf 50-200) -- photons stick at boundaries from float32 ulp issues at sub-µm steps eps=1e-4 -> 474k hits (count diff -0.40%, dirY chi2/ndf 6.58, t 80.7) eps=1e-3 -> 472k hits (count diff +0.14%, all spatial chi2/ndf < 1.5, t 5.22) eps=1e-2 -> 472k hits (count diff +0.08%, but t chi2/ndf 14.4) eps=5e-2 -> 474k hits (Opticks default; count diff -0.22%, t 51.0) -- photons leak through boundaries before refraction is computed accurately The Quartz trap with ABSLENGTH=100m (down from the previous 1km test value) makes absorption the natural cut-off, so OPTICKS_MAX_BOUNCE=10000 doesn't inflate the late tail any more (the G4 looper-kill artefact from earlier runs is now insignificant). All chi^2/ndf well within tolerance.

Drop: - test_trap_disc, test_trd_disc, single_photon, cherenkov sub-tests - test_trap_side.gdml (single-photon, redundant with iso tests) - test_trap_absorbing.gdml (intermediate, was for ABSLENGTH study) - test_trap_dispersive.gdml (cherenkov; subsumed by test_trap_scint) - run_validate_10evt[_1t].mac (cherenkov macros) - trap_disc / single_photon_xside torch configs Keep: - iso_trap (50k γ isotropic in trap) -- simple, validates G4Trap - iso_trd (50k γ isotropic in trd) -- simple, validates G4Trd - scintillation (5x 5 GeV e- in synthetic scintillator) -- folds Cerenkov + Scintillation + dispersive RINDEX + ABSLENGTH=100m + slanted walls + multi-bounce in one realistic test Two iso runs cover both new solids under heavy TIR; the scintillation run covers every optical physics process together. Total wall time ~5 min. Reviewer entry point: ./tests/g4trap_validation.sh Sub-tests: ./tests/g4trap_validation.sh {trap|trd|scintillation}

Adds test_trd_scint.gdml (same scintillating-Quartz setup as test_trap_scint.gdml but with a G4Trd solid instead of G4Trap) and splits the scintillation sub-test into: scintillation_trap -- 5 GeV e- in trap_scint gdml (472693 G4 hits) scintillation_trd -- 5 GeV e- in trd_scint gdml (892626 G4 hits) Both folded under the dispatcher case "scintillation" so existing invocations keep working. Empirical result on trd full-physics: count: 892626 vs 893156 = 0.06% diff x,y: chi2/ndf 0.96, 0.58 dx,dy,dz: chi2/ndf 0.89, 0.98, 0.68 All pos/dir chi2/ndf well under 1.1.

Upstream renamed the share dir from 'eic-opticks' to 'simphony' (PR #322 or similar). After rebase the validation script's EIC_OPTICKS_CFG default was pointing at the old path; binaries couldn't find trap_iso.json and aborted with std::runtime_error "Could not find config file".

The two iso-test geometries (test_trap.gdml, test_trd.gdml) had ABSLENGTH=1km while the scint geometries had 100m — inconsistent and exposed Opticks to a small bounce-budget asymmetry vs G4. With both at 100m the test suite is uniform and trap iso rel diff tightens from 0.27% to 0.036% (trd: 0.065% -> 0.13%, also fine).

Adds non-zero theta and phi to the alpha-sheared scint-trap geometry so the test exercises init_Trap's cx_bot/cy_bot/cx_top/cy_top tilt offsets, which the previous theta=phi=0 setup left at zero. Same material + scint + ABSLENGTH config otherwise. Validation result tightens from 0.14% to 0.006% count diff (max pos/dir chi2/ndf <1).

ggalgoczi · 2026-05-15T17:50:46Z

Since then pushed a new, more general trapezoid. Khi2 tests still pass but the numbers will be different from mentioned in the initial PR comment.

plexoos

Thank you for pushing this long-awaited support forward. It’s great to see G4Trap/G4Trd moving ahead.

A couple of immediate follow-ups before merge:

Please add the new validation script to CI so we know it passes in the standard environments, not just locally.
Have you tested this on the original pfRICH geometry with trapezoid volumes? The standalone validation cases look promising, but I’d like to confirm behavior in the motivating geometry as well.

ggalgoczi · 2026-05-15T20:05:34Z

Thank you for pushing this long-awaited support forward. It’s great to see G4Trap/G4Trd moving ahead.

A couple of immediate follow-ups before merge:
* Please add the new validation script to CI so we know it passes in the standard environments, not just locally.

* Have you tested this on the original pfRICH geometry with trapezoid volumes? The standalone validation cases look promising, but I’d like to confirm behavior in the motivating geometry as well.

Thank you for checking. I did not test this with pfRICH, I am testing dRICH now. I am not sure whether the best is to test such a feature on a complex geometry. Do you want to have a pfrich validation in this specific PR? Will add CI.

plexoos · 2026-05-15T20:50:23Z

I did not test this with pfRICH, I am testing dRICH now. I am not sure whether the best is to test such a feature on a complex geometry. Do you want to have a pfrich validation in this specific PR?

I think testing this on dRICH is fine, as long as it includes pyramid mirrors like those in pfRICH. After all, these trapezoid volumes were the primary motivation for this change.

I am not sure why you describe the pfRICH geometry as complex. From my recollection, we used it for the performance studies, and those jobs were relatively short. We have capacity to add more (longer) tests, since the current test stage is still quite quick.

plexoos · 2026-05-15T20:53:45Z

@codex review

ggalgoczi · 2026-05-15T21:13:29Z

I did not test this with pfRICH, I am testing dRICH now. I am not sure whether the best is to test such a feature on a complex geometry. Do you want to have a pfrich validation in this specific PR?

I think testing this on dRICH is fine, as long as it includes pyramid mirrors like those in pfRICH. After all, these trapezoid volumes were the primary motivation for this change.

I am not sure why you describe the pfRICH geometry as complex. From my recollection, we used it for the performance studies, and those jobs were relatively short. We have capacity to add more (longer) tests, since the current test stage is still quite quick.

How do you define complex? pfRICH geom had hundreds of volumes before simplifying. Even after simplifying it had more than a 100. I am not going to add pfRICH or dRICH test in this PR. Can follow up with an other PR later perhaps.

Copilot

Pull request overview

This PR extends the geometry conversion pipeline to support G4Trap and G4Trd by converting them into a plane-based CSG_CONVEXPOLYHEDRON, passing those planes through import into the CSGFoundry, and validating behavior by comparing Geant4 vs GPU hit outputs.

Changes:

Add G4Trap/G4Trd conversion in u4/U4Solid.h by computing 8 vertices, deriving 6 outward planes, and creating an sn::ConvexPolyhedron leaf.
Extend sn and CSGImport to carry/import per-leaf plane data into the foundry plane buffer with planeIdx/planeNum.
Add a new validation suite (GDML + macros + driver script) and run it in the PR CI workflow; extend GPURaytrace to write per-event Geant4 hits for diffable comparisons.

Reviewed changes

Copilot reviewed 12 out of 12 changed files in this pull request and generated 4 comments.

Show a summary per file

File	Description
u4/U4Solid.h	Adds `G4Trap`/`G4Trd` handling and convex-polyhedron plane construction from vertices.
sysrap/sn.h	Adds optional plane payload to `sn` plus a `ConvexPolyhedron` factory and keeps AABB for that leaf type.
CSG/CSGImport.cc	Imports plane payload into `CSGFoundry` and sets `planeIdx/planeNum` on nodes.
src/GPURaytrace.h	Adds mutex-guarded Geant4 hit writing to `g4_hits_output.txt`.
tests/g4trap_validation.sh	New end-to-end Geant4 vs GPU validation driver (torch + physics cases) with an embedded comparison script.
tests/run_validate.mac	New Geant4 macro enabling optical boundary SD invocation and running 1 event.
tests/run_validate_5evt_1t.mac	New Geant4 macro for deterministic single-thread 5-event physics runs.
tests/geom/test_trap.gdml	New simple `G4Trap` test geometry for isotropic torch validation.
tests/geom/test_trd.gdml	New simple `G4Trd` test geometry for isotropic torch validation.
tests/geom/test_trap_scint.gdml	New `G4Trap` test geometry exercising Cherenkov+scintillation in synthetic scintillator setup.
tests/geom/test_trd_scint.gdml	New `G4Trd` test geometry exercising Cherenkov+scintillation in synthetic scintillator setup.
.github/workflows/build-pull-request.yaml	Adds the new validation suite to CI test runs.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

+GEOM_DIR="${SCRIPT_DIR}/geom"
+OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation}
+
+export CUDA_VISIBLE_DEVICES=${CUDA_VISIBLE_DEVICES:-1}


+#   ./tests/g4trap_validation.sh single_photon            # 1-photon golden
+#   ./tests/g4trap_validation.sh cherenkov                # Cherenkov from e-


+                    {
+                        const PhotonHit* p = (*phc)[j];
+                        g4OutFile << p->ftime << " "
+                                  << 1239.84 / p->fenergy << "  "
+                                  << "(" << p->fposition.x() << ", " << p->fposition.y() << ", " << p->fposition.z() << ")  "
+                                  << "(" << p->fdirection.x() << ", " << p->fdirection.y() << ", " << p->fdirection.z() << ")  "
+                                  << "(" << p->fpolarization.x() << ", " << p->fpolarization.y() << ", " << p->fpolarization.z() << ")"
+                                  << std::endl;


+        // ConvexPolyhedron: store planes in foundry, set planeIdx/planeNum on node
+        const std::vector<double>* pl = nd->getPlanes();
+        unsigned                   num_planes = pl->size() / 4;
+        unsigned                   planeIdx = fd->plan.size(); // 0-based, matches csg_intersect_leaf_convexpolyhedron.h


plexoos · 2026-05-15T22:36:58Z

How do you define complex?

In this case, processing time is a useful proxy for complexity.

I am not going to add pfRICH or dRICH tests in this PR.

NO problem. I'll add them once I figure out how to run those tests.

ggalgoczi added 18 commits May 15, 2026 16:47

test(script): pass macro path to run_cherenkov_or_scint so scint uses…

62a9f7d

… 5-event macro

test(script): cherenkov uses MT macro now that the G4 hit writer is m…

80fbde8

…utex-guarded

test(scint): final yield=10 photons/MeV (validated G4-vs-GPU 0.91% di…

c4bce92

…ff, chi2/ndf<5)

ggalgoczi self-assigned this May 15, 2026

ggalgoczi added the enhancement New feature or request label May 15, 2026

ggalgoczi added 2 commits May 15, 2026 17:17

style: clang-format touched lines (Microsoft style, repo .clang-format)

133ab63

ggalgoczi force-pushed the g4trap-to-convexpolyhedron branch 2 times, most recently from 1053771 to d79ed3c Compare May 15, 2026 17:58

ggalgoczi requested a review from plexoos May 15, 2026 18:02

ggalgoczi changed the title ~~G4trap and G4trd to convexpolyhedron~~ Implement G4Trap and G4Trd support in geometries through convexpolyhedron May 15, 2026

plexoos requested changes May 15, 2026

View reviewed changes

ci: run g4trap_validation.sh in PR test matrix

8dde97e

plexoos requested a review from Copilot May 15, 2026 22:28

Copilot started reviewing on behalf of plexoos May 15, 2026 22:29 View session

Copilot AI reviewed May 15, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement G4Trap and G4Trd support in geometries through convexpolyhedron#331

Implement G4Trap and G4Trd support in geometries through convexpolyhedron#331
ggalgoczi wants to merge 21 commits into
mainfrom
g4trap-to-convexpolyhedron

ggalgoczi commented May 15, 2026 •

edited

Loading

Uh oh!

ggalgoczi commented May 15, 2026

Uh oh!

plexoos left a comment

Uh oh!

ggalgoczi commented May 15, 2026

Uh oh!

plexoos commented May 15, 2026

Uh oh!

plexoos commented May 15, 2026

Uh oh!

ggalgoczi commented May 15, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

plexoos commented May 15, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

		# ./tests/g4trap_validation.sh single_photon # 1-photon golden
		# ./tests/g4trap_validation.sh cherenkov # Cherenkov from e-

Conversation

ggalgoczi commented May 15, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Changes

Geometry support (CSG + sysrap + u4)

G4 hit writer for GPURaytrace (src/GPURaytrace.h)

Test infrastructure (tests/)

Validation results

Uh oh!

ggalgoczi commented May 15, 2026

Uh oh!

plexoos left a comment

Choose a reason for hiding this comment

Uh oh!

ggalgoczi commented May 15, 2026

Uh oh!

plexoos commented May 15, 2026

Uh oh!

plexoos commented May 15, 2026

Uh oh!

ggalgoczi commented May 15, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

plexoos commented May 15, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

ggalgoczi commented May 15, 2026 •

edited

Loading

G4 hit writer for GPURaytrace (`src/GPURaytrace.h`)

Test infrastructure (`tests/`)