Skip to content

Commit cfdf88d

Browse files
ggalgocziplexoos
andauthored
Implement G4Trap and G4Trd support in geometries through convexpolyhedron (#331)
## 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 | --------- Co-authored-by: Dmitri Smirnov <dmixsmi@gmail.com>
1 parent 2fcfe1f commit cfdf88d

14 files changed

Lines changed: 824 additions & 17 deletions

.github/workflows/build-pull-request.yaml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,7 @@ jobs:
189189
docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonFileSource.sh
190190
docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonSource_8x8SiPM.sh
191191
docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_wavelength_shifting.sh
192+
docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/g4trap_validation.sh
192193
193194
- name: Cleanup local test image
194195
if: ${{ success() }}

CSG/CSGImport.cc

Lines changed: 22 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -534,11 +534,13 @@ CSGNode* CSGImport::importNode(int nodeOffset, int partIdx, const snode& node, c
534534

535535
int typecode = nd ? nd->typecode : CSG_ZERO ;
536536
bool leaf = CSG::IsLeaf(typecode) ;
537+
bool has_planes = nd && nd->getPlanes() && nd->getPlanes()->size() > 0;
537538

538539
bool external_bbox_is_expected = CSG::ExpectExternalBBox(typecode);
539540
// CSG_CONVEXPOLYHEDRON, CSG_CONTIGUOUS, CSG_DISCONTIGUOUS, CSG_OVERLAP
540541

541-
bool expect = external_bbox_is_expected == false ;
542+
// Allow CSG_CONVEXPOLYHEDRON when sn node carries plane data
543+
bool expect = external_bbox_is_expected == false || has_planes;
542544
LOG_IF(fatal, !expect)
543545
<< " NOT EXPECTING LEAF WITH EXTERNAL BBOX EXPECTED "
544546
<< " for node of type " << CSG::Name(typecode)
@@ -561,7 +563,25 @@ CSGNode* CSGImport::importNode(int nodeOffset, int partIdx, const snode& node, c
561563
n->setBoundary(node.boundary);
562564
n->setComplement( nd ? nd->complement : false );
563565
n->setTransform(tranIdx);
564-
n->setParam_Narrow( nd ? nd->getPA_data() : nullptr );
566+
567+
if (has_planes)
568+
{
569+
// ConvexPolyhedron: store planes in foundry, set planeIdx/planeNum on node
570+
const std::vector<double>* pl = nd->getPlanes();
571+
unsigned num_planes = pl->size() / 4;
572+
unsigned planeIdx = fd->plan.size(); // 0-based, matches csg_intersect_leaf_convexpolyhedron.h
573+
for (unsigned i = 0; i < num_planes; i++)
574+
{
575+
float4 plane = make_float4((*pl)[i * 4 + 0], (*pl)[i * 4 + 1], (*pl)[i * 4 + 2], (*pl)[i * 4 + 3]);
576+
fd->addPlan(plane);
577+
}
578+
n->setPlaneIdx(planeIdx);
579+
n->setPlaneNum(num_planes);
580+
}
581+
else
582+
{
583+
n->setParam_Narrow(nd ? nd->getPA_data() : nullptr);
584+
}
565585
n->setAABB_Narrow(aabb ? aabb : nullptr );
566586

567587
return n ;

config/trap_iso.json

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
{
2+
"torch": {
3+
"gentype": "TORCH",
4+
"trackid": 0,
5+
"matline": 0,
6+
"numphoton": 50000,
7+
8+
"pos": [20.0, 20.0, -50.0],
9+
"time": 0.0,
10+
11+
"mom": [0.0, 0.0, 1.0],
12+
"weight": 0.0,
13+
14+
"pol": [1.0, 0.0, 0.0],
15+
"wavelength": 420.0,
16+
17+
"zenith": [0.0, 1.0],
18+
"azimuth": [0.0, 1.0],
19+
20+
"radius": 0.1,
21+
"distance": 0.0,
22+
"mode": 255,
23+
"type": "sphere_marsaglia"
24+
},
25+
26+
"event": {
27+
"mode": "DebugLite",
28+
"maxslot": 1000000,
29+
"max_bounce": 10000,
30+
"propagate_epsilon": 0.001,
31+
"propagate_epsilon0": 0.001
32+
}
33+
}

src/GPURaytrace.h

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@
4747
namespace
4848
{
4949
G4Mutex genstep_mutex = G4MUTEX_INITIALIZER;
50+
G4Mutex g4_hits_file_mutex = G4MUTEX_INITIALIZER;
5051
}
5152

5253
bool IsSubtractionSolid(G4VSolid *solid)
@@ -310,6 +311,8 @@ struct PrimaryGenerator : G4VUserPrimaryGeneratorAction
310311
}
311312
};
312313

314+
constexpr double kHC_eVnm = 1239.8418754200; // h*c in eV*nm (same as smath::hc_eVnm)
315+
313316
struct EventAction : G4UserEventAction
314317
{
315318
SEvt *sev;
@@ -328,14 +331,38 @@ struct EventAction : G4UserEventAction
328331
G4HCofThisEvent *hce = event->GetHCofThisEvent();
329332
if (hce)
330333
{
334+
G4AutoLock lock(&g4_hits_file_mutex);
335+
336+
static bool first_event = true;
337+
std::ofstream g4OutFile;
338+
g4OutFile.open("g4_hits_output.txt",
339+
first_event ? std::ios::out : std::ios::app);
340+
first_event = false;
341+
331342
for (G4int i = 0; i < hce->GetNumberOfCollections(); i++)
332343
{
333344
G4VHitsCollection *hc = hce->GetHC(i);
334345
if (hc)
335346
{
336347
fTotalG4Hits += hc->GetSize();
337348
}
349+
PhotonHitsCollection* phc = dynamic_cast<PhotonHitsCollection*>(hc);
350+
if (phc && g4OutFile.is_open())
351+
{
352+
for (size_t j = 0; j < phc->entries(); j++)
353+
{
354+
const PhotonHit* p = (*phc)[j];
355+
g4OutFile << p->ftime << " "
356+
<< kHC_eVnm / p->fenergy << " "
357+
<< "(" << p->fposition.x() << ", " << p->fposition.y() << ", " << p->fposition.z() << ") "
358+
<< "(" << p->fdirection.x() << ", " << p->fdirection.y() << ", " << p->fdirection.z() << ") "
359+
<< "(" << p->fpolarization.x() << ", " << p->fpolarization.y() << ", " << p->fpolarization.z() << ")"
360+
<< std::endl;
361+
}
362+
}
338363
}
364+
if (g4OutFile.is_open())
365+
g4OutFile.close();
339366
}
340367
}
341368

sysrap/sn.h

Lines changed: 24 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -156,6 +156,7 @@ struct SYSRAP_API sn
156156
s_tv* xform ;
157157
s_pa* param ;
158158
s_bb* aabb ;
159+
std::vector<double>* planes; // auxiliary plane data for CSG_CONVEXPOLYHEDRON (4 doubles per plane: nx,ny,nz,d)
159160
sn* parent ; // NOT owned by this sn
160161

161162
#ifdef WITH_CHILD
@@ -470,6 +471,8 @@ struct SYSRAP_API sn
470471
static sn* ZSphere(double radius, double z1, double z2);
471472
static sn* Box3(double fullside);
472473
static sn* Box3(double fx, double fy, double fz );
474+
static sn* ConvexPolyhedron(const double* pl, unsigned num_planes, double bbmin_x, double bbmin_y, double bbmin_z, double bbmax_x, double bbmax_y, double bbmax_z);
475+
const std::vector<double>* getPlanes() const;
473476
static sn* Torus(double rmin, double rmax, double rtor, double startPhi_deg, double deltaPhi_deg );
474477

475478
static constexpr const char* sn__PhiCut_PACMAN_ALLOWED = "sn__PhiCut_PACMAN_ALLOWED" ;
@@ -1112,14 +1115,14 @@ as other ctor/dtor can change the pool while this holds on to the old stale pid
11121115
11131116
**/
11141117

1115-
inline sn::sn(int typecode_, sn* left_, sn* right_)
1116-
:
1118+
inline sn::sn(int typecode_, sn* left_, sn* right_) :
11171119
typecode(typecode_),
11181120
complement(0),
11191121
lvid(-1),
11201122
xform(nullptr),
11211123
param(nullptr),
11221124
aabb(nullptr),
1125+
planes(nullptr),
11231126
parent(nullptr),
11241127
#ifdef WITH_CHILD
11251128
#else
@@ -1173,6 +1176,7 @@ inline sn::~sn()
11731176
delete xform ;
11741177
delete param ;
11751178
delete aabb ;
1179+
delete planes;
11761180
// parent is not deleted : as it is regarded as weakly linked (ie not owned by this node)
11771181

11781182

@@ -3258,9 +3262,21 @@ inline sn* sn::Box3(double fx, double fy, double fz ) // static
32583262
return nd ;
32593263
}
32603264

3265+
inline sn* sn::ConvexPolyhedron(const double* pl, unsigned num_planes,
3266+
double bbmin_x, double bbmin_y, double bbmin_z,
3267+
double bbmax_x, double bbmax_y, double bbmax_z) // static
3268+
{
3269+
sn* nd = Create(CSG_CONVEXPOLYHEDRON);
3270+
nd->planes = new std::vector<double>(pl, pl + num_planes * 4);
3271+
nd->setPA(0., 0., 0., 0., 0., 0.);
3272+
nd->setBB(bbmin_x, bbmin_y, bbmin_z, bbmax_x, bbmax_y, bbmax_z);
3273+
return nd;
3274+
}
32613275

3262-
3263-
3276+
inline const std::vector<double>* sn::getPlanes() const
3277+
{
3278+
return planes;
3279+
}
32643280

32653281
/**
32663282
sn::Torus
@@ -5505,6 +5521,10 @@ inline void sn::setAABB_LeafFrame()
55055521
setBB( pmin.x, pmin.y, -rmax, pmax.x, pmax.y, +rmax );
55065522

55075523
}
5524+
else if (typecode == CSG_CONVEXPOLYHEDRON)
5525+
{
5526+
// AABB already set by ConvexPolyhedron factory; keep it
5527+
}
55085528
else if( typecode == CSG_UNION || typecode == CSG_INTERSECTION || typecode == CSG_DIFFERENCE )
55095529
{
55105530
setBB( 0. );

tests/g4trap_compare.py

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
#!/usr/bin/env python3
2+
"""Compare G4 vs Opticks hit files: hit count + per-axis chi2."""
3+
import re, sys, numpy as np
4+
5+
PAT = re.compile(r'^\s*([\-\d.eE+]+)\s+([\-\d.eE+]+)\s+\(([^)]+)\)\s+\(([^)]+)\)')
6+
7+
def load(p):
8+
rows = []
9+
for line in open(p):
10+
m = PAT.match(line)
11+
if not m: continue
12+
rows.append((float(m.group(1)), float(m.group(2)),
13+
*[float(x) for x in m.group(3).split(',')],
14+
*[float(x) for x in m.group(4).split(',')]))
15+
return np.array(rows, float) if rows else np.empty((0,8))
16+
17+
def chi2_1d(a, b, bins):
18+
ha, _ = np.histogram(a, bins=bins)
19+
hb, _ = np.histogram(b, bins=bins)
20+
m = (ha + hb) > 0
21+
return float(np.sum((ha[m] - hb[m])**2 / (ha[m] + hb[m]))), int(m.sum())
22+
23+
g_path, o_path, label = sys.argv[1], sys.argv[2], sys.argv[3]
24+
tolerance_count = float(sys.argv[4]) if len(sys.argv) > 4 else 5.0
25+
tolerance_chi2 = float(sys.argv[5]) if len(sys.argv) > 5 else 5.0
26+
27+
g = load(g_path); o = load(o_path)
28+
n_g, n_o = len(g), len(o)
29+
rel = abs(n_o - n_g) / ((n_o + n_g) / 2) * 100 if n_o + n_g > 0 else 0
30+
31+
print(f"=== {label} ===")
32+
print(f" G4 hits: {n_g}")
33+
print(f" Opticks hits: {n_o}")
34+
print(f" rel diff: {rel:.3f}% (tol={tolerance_count}%)")
35+
36+
fail = []
37+
if rel > tolerance_count:
38+
fail.append(f"count rel-diff {rel:.2f}% > {tolerance_count}%")
39+
if n_g == 0 or n_o == 0:
40+
print(" no hits, skip distributions")
41+
else:
42+
for col, name, bins in [
43+
(2, 'x', np.linspace(min(g[:,2].min(), o[:,2].min()), max(g[:,2].max(), o[:,2].max()), 31)),
44+
(3, 'y', np.linspace(min(g[:,3].min(), o[:,3].min()), max(g[:,3].max(), o[:,3].max()), 31)),
45+
(5, 'dx', np.linspace(-1, 1, 21)),
46+
(6, 'dy', np.linspace(-1, 1, 21)),
47+
(7, 'dz', np.linspace(-1, 1, 21)),
48+
]:
49+
chi2, ndf = chi2_1d(g[:, col], o[:, col], bins)
50+
ratio = chi2 / max(ndf, 1)
51+
marker = "FAIL" if ratio > tolerance_chi2 else "ok"
52+
print(f" {name:>2}: chi2/ndf = {chi2:7.2f}/{ndf} = {ratio:5.2f} {marker}")
53+
if ratio > tolerance_chi2:
54+
fail.append(f"{name} chi2/ndf {ratio:.2f}")
55+
56+
if fail:
57+
print(f" FAILED: {', '.join(fail)}")
58+
sys.exit(1)
59+
print(" PASS")

tests/g4trap_validation.sh

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
#!/usr/bin/env bash
2+
#
3+
# G4Trap / G4Trd geometry validation test suite.
4+
#
5+
# Compares Geant4 (CPU) and Opticks (GPU) photon hits on identical inputs,
6+
# for each of the new convex-polyhedron-based solids. The branch under test
7+
# is g4trap-to-convexpolyhedron; this script runs every check that was used
8+
# to validate it.
9+
#
10+
# Usage:
11+
# ./tests/g4trap_validation.sh # all tests
12+
# ./tests/g4trap_validation.sh trap # trap iso-source only
13+
# ./tests/g4trap_validation.sh trd # trd iso-source only
14+
# ./tests/g4trap_validation.sh scintillation # scint+Cherenkov (trap & trd)
15+
# ./tests/g4trap_validation.sh scintillation_trap # scint+Cherenkov trap only
16+
# ./tests/g4trap_validation.sh scintillation_trd # scint+Cherenkov trd only
17+
#
18+
# Pre-requisites: GPUPhotonSource and GPURaytrace on PATH (any standard
19+
# install of simphony puts them in OPTICKS_PREFIX/bin which is added to
20+
# PATH in the Dockerfile and devcontainer).
21+
22+
set -e
23+
24+
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
25+
GEOM_DIR="${SCRIPT_DIR}/geom"
26+
OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation}
27+
28+
mkdir -p "${OUT_DIR}"
29+
30+
COMPARE_PY="${SCRIPT_DIR}/g4trap_compare.py"
31+
32+
# ------------------------------------------------------------------
33+
# run helpers
34+
# ------------------------------------------------------------------
35+
G4_MACRO="${SCRIPT_DIR}/run_validate.mac"
36+
G4_MACRO_5EVT="${SCRIPT_DIR}/run_validate_5evt_1t.mac"
37+
38+
run_torch_test () {
39+
local case=$1; local gdml=$2; local cfg=$3; local seed=${4:-42}
40+
local outd="${OUT_DIR}/${case}"
41+
mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt
42+
GPUPhotonSource -g "${gdml}" -c "${cfg}" -m "${G4_MACRO}" -s ${seed} > run.log 2>&1
43+
}
44+
45+
# ------------------------------------------------------------------
46+
# tests
47+
# ------------------------------------------------------------------
48+
# (1) Simple test class: isotropic torch source in the new solid, no
49+
# scintillation / Cherenkov physics. Validates the convex-polyhedron
50+
# conversion under heavy TIR / multi-bounce. Runs on BOTH trap and trd
51+
# so each new solid is exercised.
52+
test_trap_iso () {
53+
echo
54+
echo "----- Test: trap, isotropic source (probes slanted walls) -----"
55+
run_torch_test trap_iso "${GEOM_DIR}/test_trap.gdml" trap_iso
56+
python3 "${COMPARE_PY}" "${OUT_DIR}/trap_iso/g4_hits_output.txt" "${OUT_DIR}/trap_iso/opticks_hits_output.txt" "trap iso"
57+
}
58+
59+
test_trd_iso () {
60+
echo
61+
echo "----- Test: trd, isotropic source -----"
62+
run_torch_test trd_iso "${GEOM_DIR}/test_trd.gdml" trap_iso
63+
python3 "${COMPARE_PY}" "${OUT_DIR}/trd_iso/g4_hits_output.txt" "${OUT_DIR}/trd_iso/opticks_hits_output.txt" "trd iso"
64+
}
65+
66+
# (2) Full-physics test: 5 GeV electrons in synthetic-scintillator Quartz
67+
# solid with finite ABSLENGTH=100m. Folds Cerenkov + Scintillation +
68+
# absorption + slanted walls + multi-bounce. Run on BOTH trap and trd
69+
# so each solid sees the full physics chain. Single-thread G4 for
70+
# deterministic file output; ~3 min wall time per solid.
71+
test_scintillation_trap () {
72+
echo
73+
echo "----- Test: Scintillation+Cherenkov on trap, 5 x 5 GeV electrons -----"
74+
local outd="${OUT_DIR}/scintillation_trap"
75+
mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt
76+
GPURaytrace -g "${GEOM_DIR}/test_trap_scint.gdml" -m "${G4_MACRO_5EVT}" -s 42 > run.log 2>&1
77+
python3 "${COMPARE_PY}" "${outd}/g4_hits_output.txt" "${outd}/opticks_hits_output.txt" "scintillation trap" 10 50
78+
}
79+
80+
test_scintillation_trd () {
81+
echo
82+
echo "----- Test: Scintillation+Cherenkov on trd, 5 x 5 GeV electrons -----"
83+
local outd="${OUT_DIR}/scintillation_trd"
84+
mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt
85+
GPURaytrace -g "${GEOM_DIR}/test_trd_scint.gdml" -m "${G4_MACRO_5EVT}" -s 42 > run.log 2>&1
86+
python3 "${COMPARE_PY}" "${outd}/g4_hits_output.txt" "${outd}/opticks_hits_output.txt" "scintillation trd" 10 50
87+
}
88+
89+
# ------------------------------------------------------------------
90+
# dispatch
91+
# ------------------------------------------------------------------
92+
case "${1:-all}" in
93+
trap|iso_trap) test_trap_iso ;;
94+
trd|iso_trd) test_trd_iso ;;
95+
scintillation|sc) test_scintillation_trap; test_scintillation_trd ;;
96+
scintillation_trap|sc_trap) test_scintillation_trap ;;
97+
scintillation_trd|sc_trd) test_scintillation_trd ;;
98+
all|*)
99+
test_trap_iso
100+
test_trd_iso
101+
test_scintillation_trap
102+
test_scintillation_trd
103+
;;
104+
esac

0 commit comments

Comments
 (0)