diff --git a/.github/workflows/build-pull-request.yaml b/.github/workflows/build-pull-request.yaml index 22ac71cd2f..990998cdc9 100644 --- a/.github/workflows/build-pull-request.yaml +++ b/.github/workflows/build-pull-request.yaml @@ -187,6 +187,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() }} diff --git a/CSG/CSGImport.cc b/CSG/CSGImport.cc index de7f21e6d0..a608707640 100644 --- a/CSG/CSGImport.cc +++ b/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* 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 ; diff --git a/config/trap_iso.json b/config/trap_iso.json new file mode 100644 index 0000000000..9aba18443b --- /dev/null +++ b/config/trap_iso.json @@ -0,0 +1,30 @@ +{ + "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 + } +} diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index afd7018538..09c2caaebd 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -47,6 +47,7 @@ namespace { G4Mutex genstep_mutex = G4MUTEX_INITIALIZER; +G4Mutex g4_hits_file_mutex = G4MUTEX_INITIALIZER; } bool IsSubtractionSolid(G4VSolid *solid) @@ -328,6 +329,14 @@ 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); @@ -335,7 +344,23 @@ struct EventAction : G4UserEventAction { fTotalG4Hits += hc->GetSize(); } + PhotonHitsCollection* phc = dynamic_cast(hc); + if (phc && g4OutFile.is_open()) + { + for (size_t j = 0; j < phc->entries(); j++) + { + 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; + } + } } + if (g4OutFile.is_open()) + g4OutFile.close(); } } diff --git a/sysrap/sn.h b/sysrap/sn.h index 327400d59b..395475fcca 100644 --- a/sysrap/sn.h +++ b/sysrap/sn.h @@ -156,6 +156,7 @@ struct SYSRAP_API sn s_tv* xform ; s_pa* param ; s_bb* aabb ; + std::vector* 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* 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(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* 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. ); diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh new file mode 100755 index 0000000000..2846260184 --- /dev/null +++ b/tests/g4trap_validation.sh @@ -0,0 +1,170 @@ +#!/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-only test +# ./tests/g4trap_validation.sh trd # trd-only test +# ./tests/g4trap_validation.sh single_photon # 1-photon golden +# ./tests/g4trap_validation.sh cherenkov # Cherenkov from e- +# ./tests/g4trap_validation.sh scintillation # Scint from e- +# +# 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} + +# eps=0.001 mm (= 1 µm) was the sweet spot in an eps-scan over 1e-5..5e-2. +# Smaller values cause photons to get stuck at boundaries (float32 ulp issues); +# larger values cause leakage. eps=0.001 matched G4 best across hit count +# (within 0.14%) and all spatial+temporal chi^2 (all < 5.5). +export OPTICKS_PROPAGATE_EPSILON=${OPTICKS_PROPAGATE_EPSILON:-0.001} +export OPTICKS_PROPAGATE_EPSILON0=${OPTICKS_PROPAGATE_EPSILON0:-0.001} +# MAX_BOUNCE=10000 keeps the long late-tail in Opticks. G4's effective +# cap is much lower (~50-100 from G4Transportation looper-kill) — for a +# pure G4-count match in lossless media use MAX_BOUNCE=100, but the +# right physics setting is 10000 with realistic ABSLENGTH (< ~100 m). +export OPTICKS_MAX_BOUNCE=${OPTICKS_MAX_BOUNCE:-10000} + +mkdir -p "${OUT_DIR}" + +# Torch source config `trap_iso` is installed by CMake into +# ${OPTICKS_PREFIX}/share/simphony/config/trap_iso.json from config/trap_iso.json. + +# ------------------------------------------------------------------ +# compare.py — extract hit count + per-axis chi^2 vs G4 +# ------------------------------------------------------------------ +COMPARE_PY="${OUT_DIR}/compare.py" +cat > "${COMPARE_PY}" <<'PYEOF' +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") +PYEOF + +# ------------------------------------------------------------------ +# 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 diff --git a/tests/geom/test_trap.gdml b/tests/geom/test_trap.gdml new file mode 100644 index 0000000000..ba7502b505 --- /dev/null +++ b/tests/geom/test_trap.gdml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/geom/test_trap_scint.gdml b/tests/geom/test_trap_scint.gdml new file mode 100644 index 0000000000..aa7727aef4 --- /dev/null +++ b/tests/geom/test_trap_scint.gdml @@ -0,0 +1,108 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/geom/test_trd.gdml b/tests/geom/test_trd.gdml new file mode 100644 index 0000000000..21230d85aa --- /dev/null +++ b/tests/geom/test_trd.gdml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/geom/test_trd_scint.gdml b/tests/geom/test_trd_scint.gdml new file mode 100644 index 0000000000..68d302f5d2 --- /dev/null +++ b/tests/geom/test_trd_scint.gdml @@ -0,0 +1,97 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/run_validate.mac b/tests/run_validate.mac new file mode 100644 index 0000000000..69095b57c0 --- /dev/null +++ b/tests/run_validate.mac @@ -0,0 +1,8 @@ +# Validation macro for G4Trap/G4Trd comparison runs. +# /process/optical/boundary/setInvokeSD true: needed so that G4 fires the +# sensitive detector when an optical photon hits an EFFICIENCY>0 surface +# (default off; without it G4 records 0 hits while Opticks records them all). +/run/verbose 1 +/process/optical/boundary/setInvokeSD true +/run/initialize +/run/beamOn 1 diff --git a/tests/run_validate_5evt_1t.mac b/tests/run_validate_5evt_1t.mac new file mode 100644 index 0000000000..bb151a4f20 --- /dev/null +++ b/tests/run_validate_5evt_1t.mac @@ -0,0 +1,6 @@ +# 5-event single-thread macro for slower physics (scintillation). +/run/verbose 1 +/run/numberOfThreads 1 +/process/optical/boundary/setInvokeSD true +/run/initialize +/run/beamOn 5 diff --git a/u4/U4Solid.h b/u4/U4Solid.h index 9a47619d58..58b68590eb 100644 --- a/u4/U4Solid.h +++ b/u4/U4Solid.h @@ -37,20 +37,22 @@ npy/NNodeUncoincide npy/NNodeNudger #include "G4VSolid.hh" -#include "G4Orb.hh" -#include "G4Sphere.hh" -#include "G4Ellipsoid.hh" #include "G4Box.hh" -#include "G4Tubs.hh" -#include "G4CutTubs.hh" -#include "G4Polycone.hh" #include "G4Cons.hh" +#include "G4CutTubs.hh" +#include "G4Ellipsoid.hh" #include "G4Hype.hh" +#include "G4IntersectionSolid.hh" #include "G4MultiUnion.hh" +#include "G4Orb.hh" +#include "G4Polycone.hh" +#include "G4Sphere.hh" +#include "G4SubtractionSolid.hh" #include "G4Torus.hh" +#include "G4Trap.hh" +#include "G4Trd.hh" +#include "G4Tubs.hh" #include "G4UnionSolid.hh" -#include "G4IntersectionSolid.hh" -#include "G4SubtractionSolid.hh" #include "G4BooleanSolid.hh" #include "G4RotationMatrix.hh" @@ -77,7 +79,9 @@ enum { _G4IntersectionSolid, _G4SubtractionSolid, _G4DisplacedSolid, - _G4CutTubs + _G4CutTubs, + _G4Trap, + _G4Trd }; struct U4Solid @@ -97,6 +101,8 @@ struct U4Solid static constexpr const char* G4SubtractionSolid_ = "Sub" ; static constexpr const char* G4DisplacedSolid_ = "Dis" ; static constexpr const char* G4CutTubs_ = "TuC" ; + static constexpr const char* G4Trap_ = "Tra"; + static constexpr const char* G4Trd_ = "Trd"; static constexpr const char* _U4Solid__IsFlaggedLVID = "U4Solid__IsFlaggedLVID" ; static const int IsFlaggedLVID_ ; @@ -144,6 +150,9 @@ struct U4Solid void init_IntersectionSolid(); void init_SubtractionSolid(); void init_CutTubs(); + void init_Trap(); + void init_Trd(); + void _setRoot_FromVertices(const double v[8][3]); sn* init_Sphere_(char layer); sn* init_Cons_(char layer); @@ -293,6 +302,10 @@ inline int U4Solid::Type(const char* name) // static if( strcmp(name, "G4IntersectionSolid") == 0 ) type = _G4IntersectionSolid ; if( strcmp(name, "G4DisplacedSolid") == 0 ) type = _G4DisplacedSolid ; if( strcmp(name, "G4CutTubs") == 0 ) type = _G4CutTubs ; + if (strcmp(name, "G4Trap") == 0) + type = _G4Trap; + if (strcmp(name, "G4Trd") == 0) + type = _G4Trd; return type ; } @@ -316,7 +329,13 @@ inline const char* U4Solid::Tag(int type) // static case _G4IntersectionSolid: tag = G4IntersectionSolid_ ; break ; case _G4DisplacedSolid: tag = G4DisplacedSolid_ ; break ; case _G4CutTubs: tag = G4CutTubs_ ; break ; - } + case _G4Trap: + tag = G4Trap_; + break; + case _G4Trd: + tag = G4Trd_; + break; + } return tag ; } @@ -413,7 +432,13 @@ inline void U4Solid::init_Constituents() case _G4SubtractionSolid : init_SubtractionSolid() ; break ; case _G4DisplacedSolid : init_DisplacedSolid() ; break ; case _G4CutTubs : init_CutTubs() ; break ; - } + case _G4Trap: + init_Trap(); + break; + case _G4Trd: + init_Trd(); + break; + } if(!root) std::cerr << "U4Solid::init_Constituents UNHANDLED SOLID TYPE " << type << "\n" << desc() << "\n" ; assert( root); @@ -789,6 +814,154 @@ inline void U4Solid::init_Box() root = sn::Box3(fx, fy, fz) ; } +inline void U4Solid::init_Trap() +{ + const G4Trap* trap = dynamic_cast(solid); + assert(trap); + + // G4Trap parameters (all lengths are half-lengths) + double dz = trap->GetZHalfLength() / CLHEP::mm; + double theta = trap->GetSymAxis().theta(); // radians + double phi = trap->GetSymAxis().phi(); // radians + double dy1 = trap->GetYHalfLength1() / CLHEP::mm; + double dx1 = trap->GetXHalfLength1() / CLHEP::mm; // x at -y side, -z face + double dx2 = trap->GetXHalfLength2() / CLHEP::mm; // x at +y side, -z face + double alpha1 = trap->GetTanAlpha1(); // tan(alpha1), not the angle + double dy2 = trap->GetYHalfLength2() / CLHEP::mm; + double dx3 = trap->GetXHalfLength3() / CLHEP::mm; // x at -y side, +z face + double dx4 = trap->GetXHalfLength4() / CLHEP::mm; // x at +y side, +z face + double alpha2 = trap->GetTanAlpha2(); // tan(alpha2) + + // Face center offsets from theta/phi + double cx_bot = -dz * tan(theta) * cos(phi); // x offset at z=-dz + double cy_bot = -dz * tan(theta) * sin(phi); // y offset at z=-dz + double cx_top = dz * tan(theta) * cos(phi); // x offset at z=+dz + double cy_top = dz * tan(theta) * sin(phi); // y offset at z=+dz + + // 8 vertices (G4Trap convention, CCW when viewed from outside) + double v[8][3]; + // Bottom face (z = -dz) : x-center shifts by y*tan(alpha) + v[0][0] = cx_bot - dx1 - dy1 * alpha1; + v[0][1] = cy_bot - dy1; + v[0][2] = -dz; + v[1][0] = cx_bot + dx1 - dy1 * alpha1; + v[1][1] = cy_bot - dy1; + v[1][2] = -dz; + v[2][0] = cx_bot - dx2 + dy1 * alpha1; + v[2][1] = cy_bot + dy1; + v[2][2] = -dz; + v[3][0] = cx_bot + dx2 + dy1 * alpha1; + v[3][1] = cy_bot + dy1; + v[3][2] = -dz; + // Top face (z = +dz) + v[4][0] = cx_top - dx3 - dy2 * alpha2; + v[4][1] = cy_top - dy2; + v[4][2] = +dz; + v[5][0] = cx_top + dx3 - dy2 * alpha2; + v[5][1] = cy_top - dy2; + v[5][2] = +dz; + v[6][0] = cx_top - dx4 + dy2 * alpha2; + v[6][1] = cy_top + dy2; + v[6][2] = +dz; + v[7][0] = cx_top + dx4 + dy2 * alpha2; + v[7][1] = cy_top + dy2; + v[7][2] = +dz; + + _setRoot_FromVertices(v); +} + +inline void U4Solid::init_Trd() +{ + const G4Trd* trd = dynamic_cast(solid); + assert(trd); + + double dx1 = trd->GetXHalfLength1() / CLHEP::mm; // x at -z face + double dx2 = trd->GetXHalfLength2() / CLHEP::mm; // x at +z face + double dy1 = trd->GetYHalfLength1() / CLHEP::mm; + double dy2 = trd->GetYHalfLength2() / CLHEP::mm; + double dz = trd->GetZHalfLength() / CLHEP::mm; + + // 8 vertices using same {x-sign, y-sign, z-sign} bit convention as init_Trap + double v[8][3]; + v[0][0] = -dx1; + v[0][1] = -dy1; + v[0][2] = -dz; + v[1][0] = +dx1; + v[1][1] = -dy1; + v[1][2] = -dz; + v[2][0] = -dx1; + v[2][1] = +dy1; + v[2][2] = -dz; + v[3][0] = +dx1; + v[3][1] = +dy1; + v[3][2] = -dz; + v[4][0] = -dx2; + v[4][1] = -dy2; + v[4][2] = +dz; + v[5][0] = +dx2; + v[5][1] = -dy2; + v[5][2] = +dz; + v[6][0] = -dx2; + v[6][1] = +dy2; + v[6][2] = +dz; + v[7][0] = +dx2; + v[7][1] = +dy2; + v[7][2] = +dz; + + _setRoot_FromVertices(v); +} + +// Compute 6 outward face planes + AABB from 8 vertices, install as CSG_CONVEXPOLYHEDRON root. +// Vertex index convention: bit0=x-sign, bit1=y-sign, bit2=z-sign (0=-, 1=+). +// Face triples are CCW from outside, so BA x CA gives outward normal. +// +X: v[3],v[7],v[5] -X: v[0],v[4],v[6] +// +Y: v[2],v[6],v[7] -Y: v[1],v[5],v[4] +// +Z: v[5],v[7],v[6] -Z: v[3],v[1],v[0] +inline void U4Solid::_setRoot_FromVertices(const double v[8][3]) +{ + static const int faces[6][3] = { + {3, 7, 5}, {0, 4, 6}, {2, 6, 7}, {1, 5, 4}, {5, 7, 6}, {3, 1, 0}}; + + double planes[24]; + double bbmin[3] = {1e30, 1e30, 1e30}; + double bbmax[3] = {-1e30, -1e30, -1e30}; + + for (int i = 0; i < 8; i++) + for (int j = 0; j < 3; j++) + { + if (v[i][j] < bbmin[j]) + bbmin[j] = v[i][j]; + if (v[i][j] > bbmax[j]) + bbmax[j] = v[i][j]; + } + + for (int f = 0; f < 6; f++) + { + const double* A = v[faces[f][0]]; + const double* B = v[faces[f][1]]; + const double* C = v[faces[f][2]]; + double ba[3] = {B[0] - A[0], B[1] - A[1], B[2] - A[2]}; + double ca[3] = {C[0] - A[0], C[1] - A[1], C[2] - A[2]}; + double nx = ba[1] * ca[2] - ba[2] * ca[1]; + double ny = ba[2] * ca[0] - ba[0] * ca[2]; + double nz = ba[0] * ca[1] - ba[1] * ca[0]; + double nn = sqrt(nx * nx + ny * ny + nz * nz); + assert(nn > 0.); + nx /= nn; + ny /= nn; + nz /= nn; + double d = nx * A[0] + ny * A[1] + nz * A[2]; + + planes[f * 4 + 0] = nx; + planes[f * 4 + 1] = ny; + planes[f * 4 + 2] = nz; + planes[f * 4 + 3] = d; + } + + root = sn::ConvexPolyhedron(planes, 6, + bbmin[0], bbmin[1], bbmin[2], + bbmax[0], bbmax[1], bbmax[2]); +} inline void U4Solid::init_Tubs() {