Implement G4Trap and G4Trd support in geometries through convexpolyhedron

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

@@ -189,6 +189,7 @@ jobs:
           docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonFileSource.sh
           docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonSource_8x8SiPM.sh
           docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_wavelength_shifting.sh
+          docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/g4trap_validation.sh
 
       - name: Cleanup local test image
         if: ${{ success() }}

CSG/CSGImport.cc

@@ -534,11 +534,13 @@ CSGNode* CSGImport::importNode(int nodeOffset, int partIdx, const snode& node, c
 
     int  typecode = nd ? nd->typecode : CSG_ZERO ;
     bool leaf = CSG::IsLeaf(typecode) ;
+    bool has_planes = nd && nd->getPlanes() && nd->getPlanes()->size() > 0;
 
     bool external_bbox_is_expected = CSG::ExpectExternalBBox(typecode);
     // CSG_CONVEXPOLYHEDRON, CSG_CONTIGUOUS, CSG_DISCONTIGUOUS, CSG_OVERLAP
 
-    bool expect = external_bbox_is_expected == false ;
+    // Allow CSG_CONVEXPOLYHEDRON when sn node carries plane data
+    bool expect = external_bbox_is_expected == false || has_planes;
     LOG_IF(fatal, !expect)
         << " NOT EXPECTING LEAF WITH EXTERNAL BBOX EXPECTED "
         << " for node of type " << CSG::Name(typecode)
@@ -561,7 +563,25 @@ CSGNode* CSGImport::importNode(int nodeOffset, int partIdx, const snode& node, c
     n->setBoundary(node.boundary);
     n->setComplement( nd ? nd->complement : false );
     n->setTransform(tranIdx);
-    n->setParam_Narrow( nd ? nd->getPA_data() : nullptr );
+
+    if (has_planes)
+    {
+        // 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
+        for (unsigned i = 0; i < num_planes; i++)
+        {
+            float4 plane = make_float4((*pl)[i * 4 + 0], (*pl)[i * 4 + 1], (*pl)[i * 4 + 2], (*pl)[i * 4 + 3]);
+            fd->addPlan(plane);
+        }
+        n->setPlaneIdx(planeIdx);
+        n->setPlaneNum(num_planes);
+    }
+    else
+    {
+        n->setParam_Narrow(nd ? nd->getPA_data() : nullptr);
+    }
     n->setAABB_Narrow(aabb ? aabb : nullptr  );
 
     return n ;

config/trap_iso.json

@@ -0,0 +1,33 @@
+{
+  "torch": {
+    "gentype": "TORCH",
+    "trackid": 0,
+    "matline": 0,
+    "numphoton": 50000,
+
+    "pos": [20.0, 20.0, -50.0],
+    "time": 0.0,
+
+    "mom": [0.0, 0.0, 1.0],
+    "weight": 0.0,
+
+    "pol": [1.0, 0.0, 0.0],
+    "wavelength": 420.0,
+
+    "zenith":  [0.0, 1.0],
+    "azimuth": [0.0, 1.0],
+
+    "radius": 0.1,
+    "distance": 0.0,
+    "mode": 255,
+    "type": "sphere_marsaglia"
+  },
+
+  "event": {
+    "mode": "DebugLite",
+    "maxslot": 1000000,
+    "max_bounce": 10000,
+    "propagate_epsilon": 0.001,
+    "propagate_epsilon0": 0.001
+  }
+}

src/GPURaytrace.h

@@ -47,6 +47,7 @@
 namespace
 {
 G4Mutex genstep_mutex = G4MUTEX_INITIALIZER;
+G4Mutex g4_hits_file_mutex = G4MUTEX_INITIALIZER;
 }
 
 bool IsSubtractionSolid(G4VSolid *solid)
@@ -310,6 +311,8 @@ struct PrimaryGenerator : G4VUserPrimaryGeneratorAction
     }
 };
 
+constexpr double kHC_eVnm = 1239.8418754200; // h*c in eV*nm (same as smath::hc_eVnm)
+
 struct EventAction : G4UserEventAction
 {
     SEvt *sev;
@@ -328,14 +331,38 @@ struct EventAction : G4UserEventAction
         G4HCofThisEvent *hce = event->GetHCofThisEvent();
         if (hce)
         {
+            G4AutoLock lock(&g4_hits_file_mutex);
+
+            static bool   first_event = true;
+            std::ofstream g4OutFile;
+            g4OutFile.open("g4_hits_output.txt",
+                           first_event ? std::ios::out : std::ios::app);
+            first_event = false;
+
             for (G4int i = 0; i < hce->GetNumberOfCollections(); i++)
             {
                 G4VHitsCollection *hc = hce->GetHC(i);
                 if (hc)
                 {
                     fTotalG4Hits += hc->GetSize();
                 }
+                PhotonHitsCollection* phc = dynamic_cast<PhotonHitsCollection*>(hc);
+                if (phc && g4OutFile.is_open())
+                {
+                    for (size_t j = 0; j < phc->entries(); j++)
+                    {
+                        const PhotonHit* p = (*phc)[j];
+                        g4OutFile << p->ftime << " "
+                                  << kHC_eVnm / 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;
+                    }
+                }
             }
+            if (g4OutFile.is_open())
+                g4OutFile.close();
         }
     }
 

sysrap/sn.h

@@ -156,6 +156,7 @@ struct SYSRAP_API sn
     s_tv* xform ;
     s_pa* param ;
     s_bb* aabb  ;
+    std::vector<double>* planes;    // auxiliary plane data for CSG_CONVEXPOLYHEDRON (4 doubles per plane: nx,ny,nz,d)
     sn*   parent ;   // NOT owned by this sn
 
 #ifdef WITH_CHILD
@@ -470,6 +471,8 @@ struct SYSRAP_API sn
     static sn* ZSphere(double radius, double z1, double z2);
     static sn* Box3(double fullside);
     static sn* Box3(double fx, double fy, double fz );
+    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);
+    const std::vector<double>* getPlanes() const;
     static sn* Torus(double rmin, double rmax, double rtor, double startPhi_deg, double deltaPhi_deg );
 
     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
 
 **/
 
-inline sn::sn(int typecode_, sn* left_, sn* right_)
-    :
+inline sn::sn(int typecode_, sn* left_, sn* right_) :
     typecode(typecode_),
     complement(0),
     lvid(-1),
     xform(nullptr),
     param(nullptr),
     aabb(nullptr),
+    planes(nullptr),
     parent(nullptr),
 #ifdef WITH_CHILD
 #else
@@ -1173,6 +1176,7 @@ inline sn::~sn()
     delete xform ;
     delete param ;
     delete aabb  ;
+    delete planes;
     // parent is not deleted : as it is regarded as weakly linked (ie not owned by this node)
 
 
@@ -3258,9 +3262,21 @@ inline sn* sn::Box3(double fx, double fy, double fz )  // static
     return nd ;
 }
 
+inline sn* 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) // static
+{
+    sn* nd = Create(CSG_CONVEXPOLYHEDRON);
+    nd->planes = new std::vector<double>(pl, pl + num_planes * 4);
+    nd->setPA(0., 0., 0., 0., 0., 0.);
+    nd->setBB(bbmin_x, bbmin_y, bbmin_z, bbmax_x, bbmax_y, bbmax_z);
+    return nd;
+}
 
-
-
+inline const std::vector<double>* sn::getPlanes() const
+{
+    return planes;
+}
 
 /**
 sn::Torus
@@ -5505,6 +5521,10 @@ inline void sn::setAABB_LeafFrame()
         setBB( pmin.x, pmin.y, -rmax, pmax.x, pmax.y, +rmax );
 
     }
+    else if (typecode == CSG_CONVEXPOLYHEDRON)
+    {
+        // AABB already set by ConvexPolyhedron factory; keep it
+    }
     else if( typecode == CSG_UNION || typecode == CSG_INTERSECTION || typecode == CSG_DIFFERENCE )
     {
         setBB( 0. );

tests/g4trap_compare.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+"""Compare G4 vs Opticks hit files: hit count + per-axis chi2."""
+import re, sys, numpy as np
+
+PAT = re.compile(r'^\s*([\-\d.eE+]+)\s+([\-\d.eE+]+)\s+\(([^)]+)\)\s+\(([^)]+)\)')
+
+def load(p):
+    rows = []
+    for line in open(p):
+        m = PAT.match(line)
+        if not m: continue
+        rows.append((float(m.group(1)), float(m.group(2)),
+                     *[float(x) for x in m.group(3).split(',')],
+                     *[float(x) for x in m.group(4).split(',')]))
+    return np.array(rows, float) if rows else np.empty((0,8))
+
+def chi2_1d(a, b, bins):
+    ha, _ = np.histogram(a, bins=bins)
+    hb, _ = np.histogram(b, bins=bins)
+    m = (ha + hb) > 0
+    return float(np.sum((ha[m] - hb[m])**2 / (ha[m] + hb[m]))), int(m.sum())
+
+g_path, o_path, label = sys.argv[1], sys.argv[2], sys.argv[3]
+tolerance_count = float(sys.argv[4]) if len(sys.argv) > 4 else 5.0
+tolerance_chi2  = float(sys.argv[5]) if len(sys.argv) > 5 else 5.0
+
+g = load(g_path); o = load(o_path)
+n_g, n_o = len(g), len(o)
+rel = abs(n_o - n_g) / ((n_o + n_g) / 2) * 100 if n_o + n_g > 0 else 0
+
+print(f"=== {label} ===")
+print(f"  G4 hits:      {n_g}")
+print(f"  Opticks hits: {n_o}")
+print(f"  rel diff:     {rel:.3f}%   (tol={tolerance_count}%)")
+
+fail = []
+if rel > tolerance_count:
+    fail.append(f"count rel-diff {rel:.2f}% > {tolerance_count}%")
+if n_g == 0 or n_o == 0:
+    print("  no hits, skip distributions")
+else:
+    for col, name, bins in [
+        (2, 'x',  np.linspace(min(g[:,2].min(), o[:,2].min()), max(g[:,2].max(), o[:,2].max()), 31)),
+        (3, 'y',  np.linspace(min(g[:,3].min(), o[:,3].min()), max(g[:,3].max(), o[:,3].max()), 31)),
+        (5, 'dx', np.linspace(-1, 1, 21)),
+        (6, 'dy', np.linspace(-1, 1, 21)),
+        (7, 'dz', np.linspace(-1, 1, 21)),
+    ]:
+        chi2, ndf = chi2_1d(g[:, col], o[:, col], bins)
+        ratio = chi2 / max(ndf, 1)
+        marker = "FAIL" if ratio > tolerance_chi2 else "ok"
+        print(f"  {name:>2}: chi2/ndf = {chi2:7.2f}/{ndf} = {ratio:5.2f}  {marker}")
+        if ratio > tolerance_chi2:
+            fail.append(f"{name} chi2/ndf {ratio:.2f}")
+
+if fail:
+    print(f"  FAILED: {', '.join(fail)}")
+    sys.exit(1)
+print("  PASS")

tests/g4trap_validation.sh

@@ -0,0 +1,104 @@
+#!/usr/bin/env bash
+#
+# G4Trap / G4Trd geometry validation test suite.
+#
+# Compares Geant4 (CPU) and Opticks (GPU) photon hits on identical inputs,
+# for each of the new convex-polyhedron-based solids. The branch under test
+# is g4trap-to-convexpolyhedron; this script runs every check that was used
+# to validate it.
+#
+# Usage:
+#   ./tests/g4trap_validation.sh                          # all tests
+#   ./tests/g4trap_validation.sh trap                     # trap iso-source only
+#   ./tests/g4trap_validation.sh trd                      # trd iso-source only
+#   ./tests/g4trap_validation.sh scintillation            # scint+Cherenkov (trap & trd)
+#   ./tests/g4trap_validation.sh scintillation_trap       # scint+Cherenkov trap only
+#   ./tests/g4trap_validation.sh scintillation_trd        # scint+Cherenkov trd only
+#
+# Pre-requisites: GPUPhotonSource and GPURaytrace on PATH (any standard
+# install of simphony puts them in OPTICKS_PREFIX/bin which is added to
+# PATH in the Dockerfile and devcontainer).
+
+set -e
+
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+GEOM_DIR="${SCRIPT_DIR}/geom"
+OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation}
+
+mkdir -p "${OUT_DIR}"
+
+COMPARE_PY="${SCRIPT_DIR}/g4trap_compare.py"
+
+# ------------------------------------------------------------------
+# run helpers
+# ------------------------------------------------------------------
+G4_MACRO="${SCRIPT_DIR}/run_validate.mac"
+G4_MACRO_5EVT="${SCRIPT_DIR}/run_validate_5evt_1t.mac"
+
+run_torch_test () {
+    local case=$1; local gdml=$2; local cfg=$3; local seed=${4:-42}
+    local outd="${OUT_DIR}/${case}"
+    mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt
+    GPUPhotonSource -g "${gdml}" -c "${cfg}" -m "${G4_MACRO}" -s ${seed} > run.log 2>&1
+}
+
+# ------------------------------------------------------------------
+# tests
+# ------------------------------------------------------------------
+# (1) Simple test class: isotropic torch source in the new solid, no
+#     scintillation / Cherenkov physics. Validates the convex-polyhedron
+#     conversion under heavy TIR / multi-bounce. Runs on BOTH trap and trd
+#     so each new solid is exercised.
+test_trap_iso () {
+    echo
+    echo "----- Test: trap, isotropic source (probes slanted walls) -----"
+    run_torch_test trap_iso "${GEOM_DIR}/test_trap.gdml" trap_iso
+    python3 "${COMPARE_PY}" "${OUT_DIR}/trap_iso/g4_hits_output.txt" "${OUT_DIR}/trap_iso/opticks_hits_output.txt" "trap iso"
+}
+
+test_trd_iso () {
+    echo
+    echo "----- Test: trd, isotropic source -----"
+    run_torch_test trd_iso "${GEOM_DIR}/test_trd.gdml" trap_iso
+    python3 "${COMPARE_PY}" "${OUT_DIR}/trd_iso/g4_hits_output.txt" "${OUT_DIR}/trd_iso/opticks_hits_output.txt" "trd iso"
+}
+
+# (2) Full-physics test: 5 GeV electrons in synthetic-scintillator Quartz
+#     solid with finite ABSLENGTH=100m. Folds Cerenkov + Scintillation +
+#     absorption + slanted walls + multi-bounce. Run on BOTH trap and trd
+#     so each solid sees the full physics chain. Single-thread G4 for
+#     deterministic file output; ~3 min wall time per solid.
+test_scintillation_trap () {
+    echo
+    echo "----- Test: Scintillation+Cherenkov on trap, 5 x 5 GeV electrons -----"
+    local outd="${OUT_DIR}/scintillation_trap"
+    mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt
+    GPURaytrace -g "${GEOM_DIR}/test_trap_scint.gdml" -m "${G4_MACRO_5EVT}" -s 42 > run.log 2>&1
+    python3 "${COMPARE_PY}" "${outd}/g4_hits_output.txt" "${outd}/opticks_hits_output.txt" "scintillation trap" 10 50
+}
+
+test_scintillation_trd () {
+    echo
+    echo "----- Test: Scintillation+Cherenkov on trd, 5 x 5 GeV electrons -----"
+    local outd="${OUT_DIR}/scintillation_trd"
+    mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt
+    GPURaytrace -g "${GEOM_DIR}/test_trd_scint.gdml" -m "${G4_MACRO_5EVT}" -s 42 > run.log 2>&1
+    python3 "${COMPARE_PY}" "${outd}/g4_hits_output.txt" "${outd}/opticks_hits_output.txt" "scintillation trd" 10 50
+}
+
+# ------------------------------------------------------------------
+# dispatch
+# ------------------------------------------------------------------
+case "${1:-all}" in
+    trap|iso_trap)              test_trap_iso ;;
+    trd|iso_trd)                test_trd_iso ;;
+    scintillation|sc)           test_scintillation_trap; test_scintillation_trd ;;
+    scintillation_trap|sc_trap) test_scintillation_trap ;;
+    scintillation_trd|sc_trd)   test_scintillation_trd ;;
+    all|*)
+        test_trap_iso
+        test_trd_iso
+        test_scintillation_trap
+        test_scintillation_trd
+        ;;
+esac