From 51cb3e796136c88e57f9190c6d96b4231558c5ec Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Mon, 13 Apr 2026 14:47:11 +0000 Subject: [PATCH 01/25] feat(CSG): add G4Trap support via CSG_CONVEXPOLYHEDRON plane storage 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* 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 --- CSG/CSGImport.cc | 24 +++++++++- sysrap/sn.h | 21 +++++++++ tests/geom/test_trap.gdml | 75 ++++++++++++++++++++++++++++++ u4/U4Solid.h | 96 ++++++++++++++++++++++++++++++++++++++- 4 files changed, 213 insertions(+), 3 deletions(-) create mode 100644 tests/geom/test_trap.gdml diff --git a/CSG/CSGImport.cc b/CSG/CSGImport.cc index de7f21e6d0..d68bb1fd87 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 = 1 + fd->plan.size() ; // 1-based + 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/sysrap/sn.h b/sysrap/sn.h index 327400d59b..45f56fd5ca 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" ; @@ -1120,6 +1123,7 @@ inline sn::sn(int typecode_, sn* left_, sn* right_) xform(nullptr), param(nullptr), aabb(nullptr), + planes(nullptr), parent(nullptr), #ifdef WITH_CHILD #else @@ -1173,6 +1177,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) @@ -3259,6 +3264,18 @@ inline sn* sn::Box3(double fx, double fy, double fz ) // static } +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 ; } @@ -5505,6 +5522,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/geom/test_trap.gdml b/tests/geom/test_trap.gdml new file mode 100644 index 0000000000..fe0cb2ef0d --- /dev/null +++ b/tests/geom/test_trap.gdml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/u4/U4Solid.h b/u4/U4Solid.h index 9a47619d58..59c0e7d3ac 100644 --- a/u4/U4Solid.h +++ b/u4/U4Solid.h @@ -48,6 +48,7 @@ npy/NNodeUncoincide npy/NNodeNudger #include "G4Hype.hh" #include "G4MultiUnion.hh" #include "G4Torus.hh" +#include "G4Trap.hh" #include "G4UnionSolid.hh" #include "G4IntersectionSolid.hh" #include "G4SubtractionSolid.hh" @@ -77,7 +78,8 @@ enum { _G4IntersectionSolid, _G4SubtractionSolid, _G4DisplacedSolid, - _G4CutTubs + _G4CutTubs, + _G4Trap }; struct U4Solid @@ -97,6 +99,7 @@ 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* _U4Solid__IsFlaggedLVID = "U4Solid__IsFlaggedLVID" ; static const int IsFlaggedLVID_ ; @@ -144,6 +147,7 @@ struct U4Solid void init_IntersectionSolid(); void init_SubtractionSolid(); void init_CutTubs(); + void init_Trap(); sn* init_Sphere_(char layer); sn* init_Cons_(char layer); @@ -293,6 +297,7 @@ 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 ; return type ; } @@ -316,6 +321,7 @@ 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 ; } return tag ; } @@ -413,6 +419,7 @@ 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 ; } if(!root) std::cerr << "U4Solid::init_Constituents UNHANDLED SOLID TYPE " << type << "\n" << desc() << "\n" ; @@ -790,6 +797,93 @@ inline void U4Solid::init_Box() } +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) + 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 ; + + // Compute 6 face planes from vertex triples (outward normals, CCW winding) + // Same face ordering as legacy Python (analytic/prism.py): + // +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] + + struct { int a, b, c; } faces[6] = { + {3,7,5}, {0,4,6}, {2,6,7}, {1,5,4}, {5,7,6}, {3,1,0} + }; + + double planes[24] ; // 6 planes * 4 doubles (nx,ny,nz,d) + 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++) + { + double* A = v[faces[f].a] ; + double* B = v[faces[f].b] ; + double* C = v[faces[f].c] ; + // BA = B - A, CA = C - A + 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] } ; + // n = BA x CA + 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() { const G4Tubs* tubs = dynamic_cast(solid); From 1eaf66067a1eb5aa388c4d4bff4a1decf8f39321 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Tue, 14 Apr 2026 16:37:45 +0000 Subject: [PATCH 02/25] fix(CSG): 0-based planeIdx for ConvexPolyhedron import 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. --- CSG/CSGImport.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CSG/CSGImport.cc b/CSG/CSGImport.cc index d68bb1fd87..0edc34b29d 100644 --- a/CSG/CSGImport.cc +++ b/CSG/CSGImport.cc @@ -569,7 +569,7 @@ CSGNode* CSGImport::importNode(int nodeOffset, int partIdx, const snode& node, c // ConvexPolyhedron: store planes in foundry, set planeIdx/planeNum on node const std::vector* pl = nd->getPlanes() ; unsigned num_planes = pl->size() / 4 ; - unsigned planeIdx = 1 + fd->plan.size() ; // 1-based + 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] ) ; From 2915c3d43c21359bb4cb246621642f6aa9d276a5 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 16 Apr 2026 15:18:28 +0000 Subject: [PATCH 03/25] fix(u4): correct G4Trap alpha-shear sign in init_Trap vertex computation 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. --- u4/U4Solid.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/u4/U4Solid.h b/u4/U4Solid.h index 59c0e7d3ac..a5a7992c8e 100644 --- a/u4/U4Solid.h +++ b/u4/U4Solid.h @@ -823,16 +823,16 @@ inline void U4Solid::init_Trap() // 8 vertices (G4Trap convention, CCW when viewed from outside) double v[8][3] ; - // Bottom face (z = -dz) - 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 ; + // 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 ; + 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 ; // Compute 6 face planes from vertex triples (outward normals, CCW winding) // Same face ordering as legacy Python (analytic/prism.py): From e5a9bc9551eb433f452361a9c2e47548a476b020 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 13:01:45 +0000 Subject: [PATCH 04/25] feat(u4): add G4Trd via shared ConvexPolyhedron plane helper 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. --- tests/geom/test_trd.gdml | 76 ++++++++++++++++++++++++++++++++++++++++ u4/U4Solid.h | 64 ++++++++++++++++++++++++++------- 2 files changed, 127 insertions(+), 13 deletions(-) create mode 100644 tests/geom/test_trd.gdml diff --git a/tests/geom/test_trd.gdml b/tests/geom/test_trd.gdml new file mode 100644 index 0000000000..1c5e8aa20d --- /dev/null +++ b/tests/geom/test_trd.gdml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/u4/U4Solid.h b/u4/U4Solid.h index a5a7992c8e..c436e0b5ae 100644 --- a/u4/U4Solid.h +++ b/u4/U4Solid.h @@ -49,6 +49,7 @@ npy/NNodeUncoincide npy/NNodeNudger #include "G4MultiUnion.hh" #include "G4Torus.hh" #include "G4Trap.hh" +#include "G4Trd.hh" #include "G4UnionSolid.hh" #include "G4IntersectionSolid.hh" #include "G4SubtractionSolid.hh" @@ -79,7 +80,8 @@ enum { _G4SubtractionSolid, _G4DisplacedSolid, _G4CutTubs, - _G4Trap + _G4Trap, + _G4Trd }; struct U4Solid @@ -100,6 +102,7 @@ struct U4Solid 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_ ; @@ -148,6 +151,8 @@ struct U4Solid 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); @@ -298,6 +303,7 @@ inline int U4Solid::Type(const char* name) // static 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 ; } @@ -322,6 +328,7 @@ inline const char* U4Solid::Tag(int type) // static case _G4DisplacedSolid: tag = G4DisplacedSolid_ ; break ; case _G4CutTubs: tag = G4CutTubs_ ; break ; case _G4Trap: tag = G4Trap_ ; break ; + case _G4Trd: tag = G4Trd_ ; break ; } return tag ; } @@ -420,6 +427,7 @@ inline void U4Solid::init_Constituents() 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" ; @@ -834,17 +842,49 @@ inline void U4Solid::init_Trap() 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 ; - // Compute 6 face planes from vertex triples (outward normals, CCW winding) - // Same face ordering as legacy Python (analytic/prism.py): - // +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] + _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 ; - struct { int a, b, c; } faces[6] = { + // 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] ; // 6 planes * 4 doubles (nx,ny,nz,d) + double planes[24] ; double bbmin[3] = { 1e30, 1e30, 1e30 } ; double bbmax[3] = {-1e30, -1e30, -1e30 } ; @@ -857,13 +897,11 @@ inline void U4Solid::init_Trap() for(int f = 0 ; f < 6 ; f++) { - double* A = v[faces[f].a] ; - double* B = v[faces[f].b] ; - double* C = v[faces[f].c] ; - // BA = B - A, CA = C - A + 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] } ; - // n = BA x CA 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] ; From 8e07f3c1e2b4556e0a72f03b8ee7797a14360fb7 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 13:21:09 +0000 Subject: [PATCH 05/25] test: add run_validate.mac with boundary InvokeSD for G4-vs-GPU comparison 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. --- tests/run_validate.mac | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 tests/run_validate.mac 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 From 5a17104774987853b9b1dd91b8f875b25f11f13f Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:25:19 +0000 Subject: [PATCH 06/25] test(GPURaytrace): write per-event G4 photon hits to g4_hits_output.txt 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. --- src/GPURaytrace.h | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index afd7018538..e335309d8d 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -328,6 +328,11 @@ struct EventAction : G4UserEventAction G4HCofThisEvent *hce = event->GetHCofThisEvent(); if (hce) { + std::ofstream g4OutFile; + int eventID = event->GetEventID(); + g4OutFile.open("g4_hits_output.txt", + eventID == 0 ? std::ios::out : std::ios::app); + for (G4int i = 0; i < hce->GetNumberOfCollections(); i++) { G4VHitsCollection *hc = hce->GetHC(i); @@ -335,7 +340,22 @@ 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(); } } From 74dab0d93809233fe8f0650ce936ae42c94ff973 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:25:36 +0000 Subject: [PATCH 07/25] test: G4Trap / G4Trd validation suite (G4 vs Opticks) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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} --- tests/g4trap_validation.sh | 230 +++++++++++++++++++++++++++ tests/geom/test_trap_absorbing.gdml | 75 +++++++++ tests/geom/test_trap_dispersive.gdml | 75 +++++++++ tests/geom/test_trap_scint.gdml | 91 +++++++++++ tests/geom/test_trap_side.gdml | 75 +++++++++ tests/run_validate_10evt.mac | 5 + tests/run_validate_10evt_1t.mac | 7 + 7 files changed, 558 insertions(+) create mode 100755 tests/g4trap_validation.sh create mode 100644 tests/geom/test_trap_absorbing.gdml create mode 100644 tests/geom/test_trap_dispersive.gdml create mode 100644 tests/geom/test_trap_scint.gdml create mode 100644 tests/geom/test_trap_side.gdml create mode 100644 tests/run_validate_10evt.mac create mode 100644 tests/run_validate_10evt_1t.mac diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh new file mode 100755 index 0000000000..cbc0c95a19 --- /dev/null +++ b/tests/g4trap_validation.sh @@ -0,0 +1,230 @@ +#!/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: build the branch in /opt/eic-opticks/build (or set +# EIC_OPTICKS_BIN/EIC_OPTICKS_CFG). + +set -e + +EIC_OPTICKS_BIN=${EIC_OPTICKS_BIN:-/opt/eic-opticks/bin} +EIC_OPTICKS_CFG=${EIC_OPTICKS_CFG:-/opt/eic-opticks/share/eic-opticks/config} +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +GEOM_DIR="${SCRIPT_DIR}/geom" +OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation} + +export CUDA_VISIBLE_DEVICES=${CUDA_VISIBLE_DEVICES:-1} +export OPTICKS_PROPAGATE_EPSILON=0.0001 +export OPTICKS_PROPAGATE_EPSILON0=0.0001 +export OPTICKS_MAX_BOUNCE=${OPTICKS_MAX_BOUNCE:-100} + +mkdir -p "${OUT_DIR}" + +# ------------------------------------------------------------------ +# torch source configs (auto-installed to share/ on first run) +# ------------------------------------------------------------------ +ensure_torch_config () { + local name=$1 + local content=$2 + if [ ! -f "${EIC_OPTICKS_CFG}/${name}.json" ]; then + echo "${content}" > "${EIC_OPTICKS_CFG}/${name}.json" + fi +} + +ensure_torch_config trap_disc '{ + "torch": { + "gentype": "TORCH", "trackid": 0, "matline": 0, + "numphoton": 10000, + "pos": [0.0, 0.0, -150.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": 30.0, "distance": 0.0, "mode": 255, "type": "disc" + }, + "event": {"mode": "DebugLite", "maxslot": 1000000} +}' + +ensure_torch_config trap_iso '{ + "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} +}' + +ensure_torch_config single_photon_xside '{ + "torch": { + "gentype": "TORCH", "trackid": 0, "matline": 0, + "numphoton": 1, + "pos": [0.0, 0.0, 0.0], "time": 0.0, + "mom": [0.894427, 0.447214, 0.0], "weight": 0.0, + "pol": [0.0, 0.0, 1.0], "wavelength": 420.0, + "zenith": [0.0, 0.0], "azimuth": [0.0, 0.0], + "radius": 0.0, "distance": 0.0, "mode": 255, "type": "point" + }, + "event": {"mode": "DebugLite", "maxslot": 1000000} +}' + +# ------------------------------------------------------------------ +# 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_10EVT="${SCRIPT_DIR}/run_validate_10evt_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 + "${EIC_OPTICKS_BIN}/GPUPhotonSource" -g "${gdml}" -c "${cfg}" -m "${G4_MACRO}" -s ${seed} > run.log 2>&1 +} + +run_cherenkov_or_scint () { + local case=$1; local gdml=$2; local seed=${3:-42} + local outd="${OUT_DIR}/${case}" + mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt + "${EIC_OPTICKS_BIN}/GPURaytrace" -g "${gdml}" -m "${G4_MACRO_10EVT}" -s ${seed} > run.log 2>&1 +} + +# ------------------------------------------------------------------ +# tests +# ------------------------------------------------------------------ +test_trap_disc () { + echo + echo "----- Test: trap, disc source straight +Z -----" + run_torch_test trap_disc "${GEOM_DIR}/test_trap.gdml" trap_disc + python3 "${COMPARE_PY}" "${OUT_DIR}/trap_disc/g4_hits_output.txt" "${OUT_DIR}/trap_disc/opticks_hits_output.txt" "trap disc" +} + +test_trd_disc () { + echo + echo "----- Test: trd, disc source straight +Z -----" + run_torch_test trd_disc "${GEOM_DIR}/test_trd.gdml" trap_disc + python3 "${COMPARE_PY}" "${OUT_DIR}/trd_disc/g4_hits_output.txt" "${OUT_DIR}/trd_disc/opticks_hits_output.txt" "trd disc" +} + +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" +} + +test_single_photon () { + echo + echo "----- Test: single photon normal-incidence at trap +X wall -----" + echo " expected hit: (249.5, 124.634, 0) -- both binaries must agree" + run_torch_test single_photon "${GEOM_DIR}/test_trap_side.gdml" single_photon_xside + head -2 "${OUT_DIR}/single_photon/g4_hits_output.txt" "${OUT_DIR}/single_photon/opticks_hits_output.txt" +} + +test_cherenkov () { + echo + echo "----- Test: Cherenkov from 10 x 5 GeV electrons (single-thread G4) -----" + run_cherenkov_or_scint cherenkov "${GEOM_DIR}/test_trap_dispersive.gdml" + python3 "${COMPARE_PY}" "${OUT_DIR}/cherenkov/g4_hits_output.txt" "${OUT_DIR}/cherenkov/opticks_hits_output.txt" "cherenkov" 5 8 +} + +test_scintillation () { + echo + echo "----- Test: Scintillation from 10 x 5 GeV electrons -----" + run_cherenkov_or_scint scintillation "${GEOM_DIR}/test_trap_scint.gdml" + python3 "${COMPARE_PY}" "${OUT_DIR}/scintillation/g4_hits_output.txt" "${OUT_DIR}/scintillation/opticks_hits_output.txt" "scint" 10 50 +} + +# ------------------------------------------------------------------ +# dispatch +# ------------------------------------------------------------------ +case "${1:-all}" in + trap) test_trap_disc; test_trap_iso ;; + trd) test_trd_disc; test_trd_iso ;; + single_photon|sp) test_single_photon ;; + cherenkov|ck) test_cherenkov ;; + scintillation|sc) test_scintillation ;; + all|*) + test_trap_disc + test_trd_disc + test_trap_iso + test_trd_iso + test_single_photon + test_cherenkov + test_scintillation + ;; +esac diff --git a/tests/geom/test_trap_absorbing.gdml b/tests/geom/test_trap_absorbing.gdml new file mode 100644 index 0000000000..7455cde252 --- /dev/null +++ b/tests/geom/test_trap_absorbing.gdml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/geom/test_trap_dispersive.gdml b/tests/geom/test_trap_dispersive.gdml new file mode 100644 index 0000000000..473592f3d4 --- /dev/null +++ b/tests/geom/test_trap_dispersive.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..d7ac6d3240 --- /dev/null +++ b/tests/geom/test_trap_scint.gdml @@ -0,0 +1,91 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/geom/test_trap_side.gdml b/tests/geom/test_trap_side.gdml new file mode 100644 index 0000000000..49b5e0daa9 --- /dev/null +++ b/tests/geom/test_trap_side.gdml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/run_validate_10evt.mac b/tests/run_validate_10evt.mac new file mode 100644 index 0000000000..80eae40c25 --- /dev/null +++ b/tests/run_validate_10evt.mac @@ -0,0 +1,5 @@ +# Validation macro for multi-event Cherenkov tests. +/run/verbose 1 +/process/optical/boundary/setInvokeSD true +/run/initialize +/run/beamOn 10 diff --git a/tests/run_validate_10evt_1t.mac b/tests/run_validate_10evt_1t.mac new file mode 100644 index 0000000000..f231760e25 --- /dev/null +++ b/tests/run_validate_10evt_1t.mac @@ -0,0 +1,7 @@ +# Single-thread validation macro: avoids race in the GPURaytrace G4 hit writer +# (the writer is not MT-safe and otherwise loses ~12% of hits to the race). +/run/verbose 1 +/run/numberOfThreads 1 +/process/optical/boundary/setInvokeSD true +/run/initialize +/run/beamOn 10 From 2e02cf2acf2944fc1a0ed28ed927c4c30a12c5cf Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:27:24 +0000 Subject: [PATCH 08/25] test(scint): drop synthetic scint yield to 100 photons/MeV + add 5-event 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. --- tests/geom/test_trap_scint.gdml | 2 +- tests/run_validate_5evt_1t.mac | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) create mode 100644 tests/run_validate_5evt_1t.mac diff --git a/tests/geom/test_trap_scint.gdml b/tests/geom/test_trap_scint.gdml index d7ac6d3240..6af7efa587 100644 --- a/tests/geom/test_trap_scint.gdml +++ b/tests/geom/test_trap_scint.gdml @@ -12,7 +12,7 @@ constants accessed via GetConstProperty. --> - + 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 From f31170eef83363a50e6c8e0b50c4152bda8c2e95 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:28:03 +0000 Subject: [PATCH 09/25] test(script): pass macro path to run_cherenkov_or_scint so scint uses 5-event macro --- tests/g4trap_validation.sh | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index cbc0c95a19..416575fd95 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -141,6 +141,7 @@ PYEOF # ------------------------------------------------------------------ G4_MACRO="${SCRIPT_DIR}/run_validate.mac" G4_MACRO_10EVT="${SCRIPT_DIR}/run_validate_10evt_1t.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} @@ -150,10 +151,10 @@ run_torch_test () { } run_cherenkov_or_scint () { - local case=$1; local gdml=$2; local seed=${3:-42} + local case=$1; local gdml=$2; local macro=$3; local seed=${4:-42} local outd="${OUT_DIR}/${case}" mkdir -p "${outd}"; cd "${outd}"; rm -f *_hits_output.txt - "${EIC_OPTICKS_BIN}/GPURaytrace" -g "${gdml}" -m "${G4_MACRO_10EVT}" -s ${seed} > run.log 2>&1 + "${EIC_OPTICKS_BIN}/GPURaytrace" -g "${gdml}" -m "${macro}" -s ${seed} > run.log 2>&1 } # ------------------------------------------------------------------ @@ -198,14 +199,14 @@ test_single_photon () { test_cherenkov () { echo echo "----- Test: Cherenkov from 10 x 5 GeV electrons (single-thread G4) -----" - run_cherenkov_or_scint cherenkov "${GEOM_DIR}/test_trap_dispersive.gdml" + run_cherenkov_or_scint cherenkov "${GEOM_DIR}/test_trap_dispersive.gdml" "${G4_MACRO_10EVT}" python3 "${COMPARE_PY}" "${OUT_DIR}/cherenkov/g4_hits_output.txt" "${OUT_DIR}/cherenkov/opticks_hits_output.txt" "cherenkov" 5 8 } test_scintillation () { echo - echo "----- Test: Scintillation from 10 x 5 GeV electrons -----" - run_cherenkov_or_scint scintillation "${GEOM_DIR}/test_trap_scint.gdml" + echo "----- Test: Scintillation from 5 x 5 GeV electrons (single-thread G4) -----" + run_cherenkov_or_scint scintillation "${GEOM_DIR}/test_trap_scint.gdml" "${G4_MACRO_5EVT}" python3 "${COMPARE_PY}" "${OUT_DIR}/scintillation/g4_hits_output.txt" "${OUT_DIR}/scintillation/opticks_hits_output.txt" "scint" 10 50 } From 73d23fe7a70f839af66a42c30d1cb0ea4473e72e Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:30:03 +0000 Subject: [PATCH 10/25] fix(GPURaytrace): mutex-guard G4 hits file writer for MT runs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- src/GPURaytrace.h | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index e335309d8d..15807484f5 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,10 +329,13 @@ struct EventAction : G4UserEventAction G4HCofThisEvent *hce = event->GetHCofThisEvent(); if (hce) { + G4AutoLock lock(&g4_hits_file_mutex); + + static bool first_event = true; std::ofstream g4OutFile; - int eventID = event->GetEventID(); g4OutFile.open("g4_hits_output.txt", - eventID == 0 ? std::ios::out : std::ios::app); + first_event ? std::ios::out : std::ios::app); + first_event = false; for (G4int i = 0; i < hce->GetNumberOfCollections(); i++) { From f1f446b98bb7e79d2fc1c3e23567b340d0e4d83c Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:34:36 +0000 Subject: [PATCH 11/25] test(script): cherenkov uses MT macro now that the G4 hit writer is mutex-guarded --- tests/g4trap_validation.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index 416575fd95..2018dd9416 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -198,7 +198,7 @@ test_single_photon () { test_cherenkov () { echo - echo "----- Test: Cherenkov from 10 x 5 GeV electrons (single-thread G4) -----" + echo "----- Test: Cherenkov from 10 x 5 GeV electrons -----" run_cherenkov_or_scint cherenkov "${GEOM_DIR}/test_trap_dispersive.gdml" "${G4_MACRO_10EVT}" python3 "${COMPARE_PY}" "${OUT_DIR}/cherenkov/g4_hits_output.txt" "${OUT_DIR}/cherenkov/opticks_hits_output.txt" "cherenkov" 5 8 } From b2e0c443223026d3fcceb17955d2eb84bea2690b Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:45:40 +0000 Subject: [PATCH 12/25] test(scint): final yield=10 photons/MeV (validated G4-vs-GPU 0.91% diff, chi2/ndf<5) --- tests/geom/test_trap_scint.gdml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/geom/test_trap_scint.gdml b/tests/geom/test_trap_scint.gdml index 6af7efa587..629dc579ad 100644 --- a/tests/geom/test_trap_scint.gdml +++ b/tests/geom/test_trap_scint.gdml @@ -12,7 +12,7 @@ constants accessed via GetConstProperty. --> - + From 7aa279440dbf69273ba58b2fc70b705a5dbf74f1 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:59:31 +0000 Subject: [PATCH 13/25] fix(scint gdml): add Geant4 10.x property aliases for Opticks U4Scint 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. --- tests/geom/test_trap_scint.gdml | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/tests/geom/test_trap_scint.gdml b/tests/geom/test_trap_scint.gdml index 629dc579ad..4388d6aa8e 100644 --- a/tests/geom/test_trap_scint.gdml +++ b/tests/geom/test_trap_scint.gdml @@ -16,6 +16,7 @@ + @@ -37,12 +38,25 @@ - + + + + + + + + + From f7389989771d1e5fdc3559abb55b9dbaee9ba026 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 15:43:50 +0000 Subject: [PATCH 14/25] test(eps): change defaults to eps=0.001 + MAX_BOUNCE=10000, ABSLENGTH=100m MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- tests/g4trap_validation.sh | 14 +++++++++++--- tests/geom/test_trap_scint.gdml | 5 ++++- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index 2018dd9416..94c43053d6 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -27,9 +27,17 @@ GEOM_DIR="${SCRIPT_DIR}/geom" OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation} export CUDA_VISIBLE_DEVICES=${CUDA_VISIBLE_DEVICES:-1} -export OPTICKS_PROPAGATE_EPSILON=0.0001 -export OPTICKS_PROPAGATE_EPSILON0=0.0001 -export OPTICKS_MAX_BOUNCE=${OPTICKS_MAX_BOUNCE:-100} +# 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}" diff --git a/tests/geom/test_trap_scint.gdml b/tests/geom/test_trap_scint.gdml index 4388d6aa8e..2e31983f00 100644 --- a/tests/geom/test_trap_scint.gdml +++ b/tests/geom/test_trap_scint.gdml @@ -5,7 +5,10 @@ - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/geom/test_trap_dispersive.gdml b/tests/geom/test_trap_dispersive.gdml deleted file mode 100644 index 473592f3d4..0000000000 --- a/tests/geom/test_trap_dispersive.gdml +++ /dev/null @@ -1,75 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/geom/test_trap_side.gdml b/tests/geom/test_trap_side.gdml deleted file mode 100644 index 49b5e0daa9..0000000000 --- a/tests/geom/test_trap_side.gdml +++ /dev/null @@ -1,75 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/run_validate_10evt.mac b/tests/run_validate_10evt.mac deleted file mode 100644 index 80eae40c25..0000000000 --- a/tests/run_validate_10evt.mac +++ /dev/null @@ -1,5 +0,0 @@ -# Validation macro for multi-event Cherenkov tests. -/run/verbose 1 -/process/optical/boundary/setInvokeSD true -/run/initialize -/run/beamOn 10 diff --git a/tests/run_validate_10evt_1t.mac b/tests/run_validate_10evt_1t.mac deleted file mode 100644 index f231760e25..0000000000 --- a/tests/run_validate_10evt_1t.mac +++ /dev/null @@ -1,7 +0,0 @@ -# Single-thread validation macro: avoids race in the GPURaytrace G4 hit writer -# (the writer is not MT-safe and otherwise loses ~12% of hits to the race). -/run/verbose 1 -/run/numberOfThreads 1 -/process/optical/boundary/setInvokeSD true -/run/initialize -/run/beamOn 10 From 18ed85f203a9ff88bf950b55687eae9df57d36dc Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 16:19:48 +0000 Subject: [PATCH 16/25] test: add full-physics scintillation test for trd 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. --- tests/g4trap_validation.sh | 36 ++++++++----- tests/geom/test_trd_scint.gdml | 97 ++++++++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 12 deletions(-) create mode 100644 tests/geom/test_trd_scint.gdml diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index 252bacc8e5..115a225fb9 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -153,29 +153,41 @@ test_trd_iso () { } # (2) Full-physics test: 5 GeV electrons in synthetic-scintillator Quartz -# trap with dispersive Sellmeier index + finite ABSLENGTH=100m. Folds -# Cherenkov + Scintillation + dispersion + absorption + slanted walls -# + multi-bounce together. Single-thread G4 for deterministic file -# output; ~3 min wall time. -test_scintillation () { +# 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 from 5 x 5 GeV electrons -----" - local outd="${OUT_DIR}/scintillation" + 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 "${EIC_OPTICKS_BIN}/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" 10 50 + 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 + "${EIC_OPTICKS_BIN}/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|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 + test_scintillation_trap + test_scintillation_trd ;; esac 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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 437212eede134c2b29d8fd3ac4ae99fcc0e66267 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:00:29 +0000 Subject: [PATCH 17/25] test(script): update config path to share/simphony/config 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". --- tests/g4trap_validation.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index 115a225fb9..5d86891972 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -21,7 +21,7 @@ set -e EIC_OPTICKS_BIN=${EIC_OPTICKS_BIN:-/opt/eic-opticks/bin} -EIC_OPTICKS_CFG=${EIC_OPTICKS_CFG:-/opt/eic-opticks/share/eic-opticks/config} +EIC_OPTICKS_CFG=${EIC_OPTICKS_CFG:-/opt/eic-opticks/share/simphony/config} SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" GEOM_DIR="${SCRIPT_DIR}/geom" OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation} From 0fcbe7f66472394ad8f150e80c7892cacb46f599 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:05:03 +0000 Subject: [PATCH 18/25] test: ABSLENGTH=100m in iso GDMLs too (was 1km; matches scint GDMLs) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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). --- tests/geom/test_trap.gdml | 2 +- tests/geom/test_trd.gdml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/geom/test_trap.gdml b/tests/geom/test_trap.gdml index fe0cb2ef0d..ba7502b505 100644 --- a/tests/geom/test_trap.gdml +++ b/tests/geom/test_trap.gdml @@ -5,7 +5,7 @@ - + diff --git a/tests/geom/test_trd.gdml b/tests/geom/test_trd.gdml index 1c5e8aa20d..21230d85aa 100644 --- a/tests/geom/test_trd.gdml +++ b/tests/geom/test_trd.gdml @@ -5,7 +5,7 @@ - + From df9af9d272771f178dc6c24f4a7260a54f044bca Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:17:19 +0000 Subject: [PATCH 19/25] style: clang-format touched lines (Microsoft style, repo .clang-format) --- CSG/CSGImport.cc | 20 ++-- src/GPURaytrace.h | 9 +- sysrap/sn.h | 33 ++++--- u4/U4Solid.h | 227 +++++++++++++++++++++++++++------------------- 4 files changed, 165 insertions(+), 124 deletions(-) diff --git a/CSG/CSGImport.cc b/CSG/CSGImport.cc index 0edc34b29d..a608707640 100644 --- a/CSG/CSGImport.cc +++ b/CSG/CSGImport.cc @@ -534,13 +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 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 // Allow CSG_CONVEXPOLYHEDRON when sn node carries plane data - bool expect = external_bbox_is_expected == false || has_planes ; + 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) @@ -564,23 +564,23 @@ CSGNode* CSGImport::importNode(int nodeOffset, int partIdx, const snode& node, c n->setComplement( nd ? nd->complement : false ); n->setTransform(tranIdx); - if( has_planes ) + 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++) + 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) ; + 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->setParam_Narrow(nd ? nd->getPA_data() : nullptr); } n->setAABB_Narrow(aabb ? aabb : nullptr ); diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index 15807484f5..09c2caaebd 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -331,7 +331,7 @@ struct EventAction : G4UserEventAction { G4AutoLock lock(&g4_hits_file_mutex); - static bool first_event = true; + static bool first_event = true; std::ofstream g4OutFile; g4OutFile.open("g4_hits_output.txt", first_event ? std::ios::out : std::ios::app); @@ -344,12 +344,12 @@ struct EventAction : G4UserEventAction { fTotalG4Hits += hc->GetSize(); } - PhotonHitsCollection *phc = dynamic_cast(hc); + PhotonHitsCollection* phc = dynamic_cast(hc); if (phc && g4OutFile.is_open()) { for (size_t j = 0; j < phc->entries(); j++) { - const PhotonHit *p = (*phc)[j]; + const PhotonHit* p = (*phc)[j]; g4OutFile << p->ftime << " " << 1239.84 / p->fenergy << " " << "(" << p->fposition.x() << ", " << p->fposition.y() << ", " << p->fposition.z() << ") " @@ -359,7 +359,8 @@ struct EventAction : G4UserEventAction } } } - if (g4OutFile.is_open()) g4OutFile.close(); + if (g4OutFile.is_open()) + g4OutFile.close(); } } diff --git a/sysrap/sn.h b/sysrap/sn.h index 45f56fd5ca..395475fcca 100644 --- a/sysrap/sn.h +++ b/sysrap/sn.h @@ -156,7 +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) + 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 @@ -472,7 +472,7 @@ struct SYSRAP_API sn 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 ; + 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" ; @@ -1115,8 +1115,7 @@ 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), @@ -1177,7 +1176,7 @@ inline sn::~sn() delete xform ; delete param ; delete aabb ; - delete planes ; + delete planes; // parent is not deleted : as it is regarded as weakly linked (ie not owned by this node) @@ -3263,21 +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 + 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 ; + 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 ; } - - +inline const std::vector* sn::getPlanes() const +{ + return planes; +} /** sn::Torus @@ -5522,7 +5521,7 @@ inline void sn::setAABB_LeafFrame() setBB( pmin.x, pmin.y, -rmax, pmax.x, pmax.y, +rmax ); } - else if( typecode == CSG_CONVEXPOLYHEDRON ) + else if (typecode == CSG_CONVEXPOLYHEDRON) { // AABB already set by ConvexPolyhedron factory; keep it } diff --git a/u4/U4Solid.h b/u4/U4Solid.h index c436e0b5ae..58b68590eb 100644 --- a/u4/U4Solid.h +++ b/u4/U4Solid.h @@ -37,22 +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" @@ -101,8 +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* G4Trap_ = "Tra"; + static constexpr const char* G4Trd_ = "Trd"; static constexpr const char* _U4Solid__IsFlaggedLVID = "U4Solid__IsFlaggedLVID" ; static const int IsFlaggedLVID_ ; @@ -302,8 +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 ; + if (strcmp(name, "G4Trap") == 0) + type = _G4Trap; + if (strcmp(name, "G4Trd") == 0) + type = _G4Trd; return type ; } @@ -327,9 +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 ; - } + case _G4Trap: + tag = G4Trap_; + break; + case _G4Trd: + tag = G4Trd_; + break; + } return tag ; } @@ -426,9 +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 ; - } + 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); @@ -804,74 +814,103 @@ 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) + 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 + 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] ; + 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 ; + 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) ; + 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 ; + 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) ; + 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. @@ -881,47 +920,49 @@ inline void U4Solid::init_Trd() 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} - }; + {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 } ; + 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++) + 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] ; + 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++) + 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 ; + 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]) ; + bbmin[0], bbmin[1], bbmin[2], + bbmax[0], bbmax[1], bbmax[2]); } - inline void U4Solid::init_Tubs() { const G4Tubs* tubs = dynamic_cast(solid); From 38e97b079948ea4bd51aa71f94559ecafabe9d9f Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:49:54 +0000 Subject: [PATCH 20/25] test(scint trap): tilt symmetry axis (theta=10deg, phi=30deg) 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). --- tests/geom/test_trap_scint.gdml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/geom/test_trap_scint.gdml b/tests/geom/test_trap_scint.gdml index 2e31983f00..aa7727aef4 100644 --- a/tests/geom/test_trap_scint.gdml +++ b/tests/geom/test_trap_scint.gdml @@ -65,7 +65,7 @@ - From 552d9e71ed2a78e82deb6a276703f95fd7bd4e2b Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 20:42:33 +0000 Subject: [PATCH 21/25] ci: run g4trap_validation.sh in PR test matrix --- .github/workflows/build-pull-request.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/build-pull-request.yaml b/.github/workflows/build-pull-request.yaml index 3cb278d5e6..ed9baa08dc 100644 --- a/.github/workflows/build-pull-request.yaml +++ b/.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() }} From 0ae5148a5da661e1f5c7d2a07731e457e2dc85d4 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Wed, 20 May 2026 14:30:27 +0000 Subject: [PATCH 22/25] chore(config): ship trap_iso torch config from source, drop runtime install The g4trap validation script previously wrote /opt/eic-opticks/share/simphony/config/trap_iso.json on the fly via an ensure_torch_config helper. That assumes the script can write into the install prefix, which fails in CI when the container user can't write share/simphony/config and aborts the script with set -e. Move the JSON to config/trap_iso.json so the existing install(DIRECTORY ./config ...) rule in the top-level CMakeLists.txt installs it alongside dev.json / wls_test.json / sphere_leak.json / 8x8SiPM_crystal.json. Drop the ensure_torch_config helper and the EIC_OPTICKS_CFG variable from tests/g4trap_validation.sh. --- config/trap_iso.json | 30 ++++++++++++++++++++++++++++++ tests/g4trap_validation.sh | 28 +++------------------------- 2 files changed, 33 insertions(+), 25 deletions(-) create mode 100644 config/trap_iso.json 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/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index 5d86891972..96757e492a 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -16,12 +16,11 @@ # ./tests/g4trap_validation.sh scintillation # Scint from e- # # Pre-requisites: build the branch in /opt/eic-opticks/build (or set -# EIC_OPTICKS_BIN/EIC_OPTICKS_CFG). +# EIC_OPTICKS_BIN). set -e EIC_OPTICKS_BIN=${EIC_OPTICKS_BIN:-/opt/eic-opticks/bin} -EIC_OPTICKS_CFG=${EIC_OPTICKS_CFG:-/opt/eic-opticks/share/simphony/config} SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" GEOM_DIR="${SCRIPT_DIR}/geom" OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation} @@ -41,29 +40,8 @@ export OPTICKS_MAX_BOUNCE=${OPTICKS_MAX_BOUNCE:-10000} mkdir -p "${OUT_DIR}" -# ------------------------------------------------------------------ -# torch source configs (auto-installed to share/ on first run) -# ------------------------------------------------------------------ -ensure_torch_config () { - local name=$1 - local content=$2 - if [ ! -f "${EIC_OPTICKS_CFG}/${name}.json" ]; then - echo "${content}" > "${EIC_OPTICKS_CFG}/${name}.json" - fi -} - -ensure_torch_config trap_iso '{ - "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} -}' +# 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 From 5055e30a8aab6daa8e458ace8fe201ecad1f1af5 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Wed, 20 May 2026 14:59:20 +0000 Subject: [PATCH 23/25] fix(tests): use unqualified binary names so g4trap_validation works in CI image The script hardcoded EIC_OPTICKS_BIN=/opt/eic-opticks/bin and invoked ${EIC_OPTICKS_BIN}/GPUPhotonSource and ${EIC_OPTICKS_BIN}/GPURaytrace via absolute path. The develop docker image sets OPTICKS_PREFIX=/opt/simphony (see Dockerfile line 67) and adds ${OPTICKS_PREFIX}/bin to PATH (line 72), so the absolute path resolves to a non-existent directory and bash exits with 127 ("command not found"). The other test scripts (test_GPURaytrace.sh, test_simg4ox.sh, ...) already rely on PATH, so do the same here: drop the EIC_OPTICKS_BIN variable and invoke GPUPhotonSource / GPURaytrace unqualified. --- tests/g4trap_validation.sh | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index 96757e492a..b81954cae5 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -15,12 +15,12 @@ # ./tests/g4trap_validation.sh cherenkov # Cherenkov from e- # ./tests/g4trap_validation.sh scintillation # Scint from e- # -# Pre-requisites: build the branch in /opt/eic-opticks/build (or set -# EIC_OPTICKS_BIN). +# 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 -EIC_OPTICKS_BIN=${EIC_OPTICKS_BIN:-/opt/eic-opticks/bin} SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" GEOM_DIR="${SCRIPT_DIR}/geom" OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation} @@ -106,7 +106,7 @@ 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 - "${EIC_OPTICKS_BIN}/GPUPhotonSource" -g "${gdml}" -c "${cfg}" -m "${G4_MACRO}" -s ${seed} > run.log 2>&1 + GPUPhotonSource -g "${gdml}" -c "${cfg}" -m "${G4_MACRO}" -s ${seed} > run.log 2>&1 } # ------------------------------------------------------------------ @@ -140,7 +140,7 @@ test_scintillation_trap () { 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 - "${EIC_OPTICKS_BIN}/GPURaytrace" -g "${GEOM_DIR}/test_trap_scint.gdml" -m "${G4_MACRO_5EVT}" -s 42 > run.log 2>&1 + 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 } @@ -149,7 +149,7 @@ test_scintillation_trd () { 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 - "${EIC_OPTICKS_BIN}/GPURaytrace" -g "${GEOM_DIR}/test_trd_scint.gdml" -m "${G4_MACRO_5EVT}" -s 42 > run.log 2>&1 + 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 } From 3f5799232e87ff1877f301ce35027f2a4a9821b9 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Wed, 20 May 2026 15:44:10 +0000 Subject: [PATCH 24/25] fix(tests): don't force CUDA_VISIBLE_DEVICES=1 in g4trap_validation.sh CI launches the test container with `docker run --gpus 'device=1'`, which exposes a single GPU as device 0 inside the container. The script's default `CUDA_VISIBLE_DEVICES=${CUDA_VISIBLE_DEVICES:-1}` then selects a non-existent device, causing CUDA init failure and SIGABRT (exit 134). Drop the default; let the environment provide CUDA_VISIBLE_DEVICES when needed. Matches what every other test script in tests/ does. --- tests/g4trap_validation.sh | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index b81954cae5..2846260184 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -25,7 +25,6 @@ SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" GEOM_DIR="${SCRIPT_DIR}/geom" OUT_DIR=${OUT_DIR:-/tmp/g4trap_validation} -export CUDA_VISIBLE_DEVICES=${CUDA_VISIBLE_DEVICES:-1} # 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 From 6f6585f3c3d6a2e334aa61d9073c71756e0c1546 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Sat, 6 Jun 2026 14:41:42 +0000 Subject: [PATCH 25/25] refactor(tests): address g4trap review feedback - Use named constant kHC_eVnm for eV-to-nm conversion in GPURaytrace hit writer - Fix validation script usage header to match actual test modes - Move propagation parameters (epsilon, max_bounce) into trap_iso.json config - Extract inline comparison Python into tests/g4trap_compare.py --- config/trap_iso.json | 5 ++- src/GPURaytrace.h | 4 +- tests/g4trap_compare.py | 59 ++++++++++++++++++++++++++++ tests/g4trap_validation.sh | 78 +++----------------------------------- 4 files changed, 72 insertions(+), 74 deletions(-) create mode 100644 tests/g4trap_compare.py diff --git a/config/trap_iso.json b/config/trap_iso.json index 9aba18443b..96eef08962 100644 --- a/config/trap_iso.json +++ b/config/trap_iso.json @@ -25,6 +25,9 @@ "event": { "mode": "DebugLite", - "maxslot": 1000000 + "maxslot": 1000000, + "max_bounce": 10000, + "propagate_epsilon": 0.001, + "propagate_epsilon0": 0.001 } } diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index 09c2caaebd..14606bfe34 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -311,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; @@ -351,7 +353,7 @@ struct EventAction : G4UserEventAction { const PhotonHit* p = (*phc)[j]; g4OutFile << p->ftime << " " - << 1239.84 / p->fenergy << " " + << 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() << ")" diff --git a/tests/g4trap_compare.py b/tests/g4trap_compare.py new file mode 100644 index 0000000000..fa1cb92b9a --- /dev/null +++ b/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") diff --git a/tests/g4trap_validation.sh b/tests/g4trap_validation.sh index 2846260184..043c897830 100755 --- a/tests/g4trap_validation.sh +++ b/tests/g4trap_validation.sh @@ -9,11 +9,11 @@ # # 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- +# ./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 @@ -25,75 +25,9 @@ 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 +COMPARE_PY="${SCRIPT_DIR}/g4trap_compare.py" # ------------------------------------------------------------------ # run helpers