From e4f225239294aa9d37b1f092916cbf5c65f6b653 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Mon, 13 Apr 2026 14:47:11 +0000 Subject: [PATCH 01/24] 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 f9a7d3ec5795bd08210558ce1042a8fecbf81792 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Tue, 14 Apr 2026 16:37:45 +0000 Subject: [PATCH 02/24] 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 1be643f74054a48d08cd91a1b8ff8bf79932d339 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 16 Apr 2026 15:18:28 +0000 Subject: [PATCH 03/24] 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 0e549a6baf4eb2e522cc77716a6b5a2c3bd2a4ae Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 13:01:45 +0000 Subject: [PATCH 04/24] 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 83b48bc32a28ccbe12ffa8f179528d5e4df718a2 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 13:21:09 +0000 Subject: [PATCH 05/24] 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 3e21db8e38a463153513fcd58a543506b82d3c7d Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:25:19 +0000 Subject: [PATCH 06/24] 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 96a4f62f41bc5ad6a76629657f166f68501eb068 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:25:36 +0000 Subject: [PATCH 07/24] 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 479d143ccc979cac88b2ab652ac2ca00a0a2883b Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:27:24 +0000 Subject: [PATCH 08/24] 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 62a9f7da4ba7ecfd7de8821a2dc3db43e30ef044 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:28:03 +0000 Subject: [PATCH 09/24] 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 8eed994cb26b07c1c7a17ecb3cc88ad789b24177 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:30:03 +0000 Subject: [PATCH 10/24] 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 80fbde87060e7a856f414a7d4890fe3924b746da Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:34:36 +0000 Subject: [PATCH 11/24] 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 c4bce927a821227b9eb2b622eadc0df88a5750ab Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:45:40 +0000 Subject: [PATCH 12/24] 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 8a3829a9babc5b841272987fbf890e4f32736458 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 14:59:31 +0000 Subject: [PATCH 13/24] 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 91efc62cc11229356257c6d5f21b00a6e4ae55d0 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 15:43:50 +0000 Subject: [PATCH 14/24] 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 1e624cb3178915cb263e9b75bbdbb4d96a67ea60 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 16:19:48 +0000 Subject: [PATCH 16/24] 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 dbead88f79857fb191fd15e3a27dfe54ebb5ccd5 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:00:29 +0000 Subject: [PATCH 17/24] 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 1d5b2391917c49d4eea97c1b40698fef0e35f301 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:05:03 +0000 Subject: [PATCH 18/24] 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 133ab63a6b117eb5d1d3f8c0e51a516757eb148e Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:17:19 +0000 Subject: [PATCH 19/24] 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 d79ed3cc43300fdd398ce7c6f2132c8a9b717e6e Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 17:49:54 +0000 Subject: [PATCH 20/24] 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 8dde97e9a9dc24e8fcc5ff196dc5687f2a72b283 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 15 May 2026 20:42:33 +0000 Subject: [PATCH 21/24] 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 22ac71cd2f..990998cdc9 100644 --- a/.github/workflows/build-pull-request.yaml +++ b/.github/workflows/build-pull-request.yaml @@ -187,6 +187,7 @@ jobs: docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonFileSource.sh docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonSource_8x8SiPM.sh docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_wavelength_shifting.sh + docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/g4trap_validation.sh - name: Cleanup local test image if: ${{ success() }} From 2fb2d446c97e9d3ea1652426c7ed6e7f9946b5c1 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Wed, 20 May 2026 14:30:27 +0000 Subject: [PATCH 22/24] 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 f0800021550c010d56d0b6cc2f0bdcddc9bab128 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Wed, 20 May 2026 14:59:20 +0000 Subject: [PATCH 23/24] 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 8db7324f385098a2078862e2b964f89026d95770 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Wed, 20 May 2026 15:44:10 +0000 Subject: [PATCH 24/24] 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