From 085ffebd10e27beaef84d81c6b629a26c4912121 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 29 May 2026 13:41:06 +0000 Subject: [PATCH 1/9] qbnd/qsim: carry source-medium matline across CSG sibling-pair faces GBndLib does not enumerate boundary entries for CSG sibling-pair faces (two daughters of one parent sharing a face), so OptiX can pick a boundary whose m1/m2 do not reference the photon's current medium -> absorption/scattering sampled in the wrong material1. Carry the photon's matline (seeded from the genstep at birth in CSGOptiX7.cu, updated each propagate_at_boundary) in sctx::current_matline and use it as m1 in qbnd::fill_state when valid. No-op for correctly-enumerated boundaries (carried == the boundary's own m1_line). --- CSGOptiX/CSGOptiX7.cu | 8 ++++++++ qudarap/qbnd.h | 24 ++++++++++++++++++++---- qudarap/qsim.h | 15 ++++++++++++++- sysrap/sctx.h | 1 + 4 files changed, 43 insertions(+), 5 deletions(-) diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index 5821cbaf7..957a18cdc 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -439,6 +439,14 @@ static __forceinline__ __device__ void simulate( const uint3& launch_idx, const #endif sim->generate_photon(ctx.p, rng, gs, photon_idx, genstep_idx ); + // F7: seed carried matline from the genstep's source-medium bnd row + // (gs.q0.u.z == scerenkov/sscint matline). Updated through + // propagate_at_boundary. Guard the -1 sentinel (0xFFFFFFFF) that + // SEvt::setGenstep emits for bad_ck/unmapped mtindex. + { + const unsigned gm = gs.q0.u.z; + ctx.current_matline = (gm == 0xFFFFFFFFu) ? 0u : gm; + } FlowAction command = FlowAction::Start; int bounce = 0 ; diff --git a/qudarap/qbnd.h b/qudarap/qbnd.h index d6c5bdd9d..fe9637a85 100644 --- a/qudarap/qbnd.h +++ b/qudarap/qbnd.h @@ -38,7 +38,7 @@ struct qbnd #if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_TEXTURE) || defined(MOCK_CUDA) QBND_METHOD float4 boundary_lookup( unsigned ix, unsigned iy ); QBND_METHOD float4 boundary_lookup( float nm, unsigned line, unsigned k ); - QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx ); + QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline = 0u); #endif }; @@ -194,14 +194,30 @@ s.optical.x **/ -inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx ) +inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline) { const int line = boundary*_BOUNDARY_NUM_MATSUR ; // now that are not signing boundary use 0-based - const int m1_line = cosTheta > 0.f ? line + IMAT : line + OMAT ; + int m1_line = cosTheta > 0.f ? line + IMAT : line + OMAT; const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT ; const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR ; + // F7 sibling-pair material override. + // GBndLib does not enumerate boundary entries for CSG sibling-pair faces + // (two daughters of one parent sharing a face), so OptiX may pick a boundary + // whose m1/m2 do not reference the photon's actual medium -> wrong material1. + // qsim::propagate passes the photon's carried matline (seeded from the + // genstep, updated each propagate_at_boundary); when it is valid we use it + // directly as m1 so absorption/scattering sample the true current medium. + // (material2 is left as the boundary's other side.) + // Guards: 0 = "no carry" (legacy callers); 0xFFFFFFFF = the lookup_mtline=-1 + // sentinel from SEvt::setGenstep, not a valid table row; and the carried + // slot must be a material row (OMAT/IMAT), not a surface row (OSUR/ISUR). + const unsigned matline_slot = carried_matline & 0x3u; + if (carried_matline != 0u && carried_matline != 0xFFFFFFFFu && (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT))) + { + m1_line = int(carried_matline); + } s.material1 = boundary_lookup( wavelength, m1_line, 0); // refractive_index, absorption_length, scattering_length, reemission_prob s.m1group2 = boundary_lookup( wavelength, m1_line, 1); // x: material1 group_velocity, y: material2 group_velocity, z/w unused @@ -211,7 +227,7 @@ inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wav s.optical = optical[su_line].u ; // 1-based-surface-index-0-meaning-boundary/type/finish/value (type,finish,value not used currently) - s.index.x = optical[m1_line].u.x ; // m1 index (1-based, see sstandard::make_optical) + s.index.x = optical[m1_line].u.x; // m1 index (1-based, see sstandard::make_optical) -- reflects override if any s.index.y = optical[m2_line].u.x ; // m2 index (1-based, see sstandard::make_optical) s.index.z = optical[su_line].u.x ; // su index (1-based, see sstandard::make_optical) s.index.w = 0u ; // avoid undefined memory comparison issues diff --git a/qudarap/qsim.h b/qudarap/qsim.h index 1e110505b..c25e06518 100644 --- a/qudarap/qsim.h +++ b/qudarap/qsim.h @@ -1261,6 +1261,19 @@ inline QSIM_METHOD FlowAction qsim::propagate_at_boundary(unsigned& flag, RNG& r ctx.current_material_index = reflect ? s.index.x : s.index.y; ctx.current_group_velocity = reflect ? s.material1_group_velocity() : s.material2_group_velocity(); + // F7 SIBLING-PAIR carry: update current_matline so propagate() at the NEXT + // boundary knows the photon's actual medium even when that boundary doesn't + // reference it (the sibling-pair case). _c1 = -dot(p.mom,normal) here, so + // _c1 < 0 maps to cosTheta > 0 in qbnd::fill_state (m1=IMAT, m2=OMAT). + { + const unsigned bnd_idx = ctx.prd->boundary(); + const unsigned imat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + IMAT; + const unsigned omat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + OMAT; + const unsigned this_m1_line = (_c1 < 0.f) ? imat_line : omat_line; + const unsigned this_m2_line = (_c1 < 0.f) ? omat_line : imat_line; + ctx.current_matline = reflect ? this_m1_line : this_m2_line; + } + #if !defined(PRODUCTION) && defined(DEBUG_TAG) if( flag == BOUNDARY_REFLECT ) { @@ -2149,7 +2162,7 @@ QSIM_FORCEINLINE_METHOD FlowAction qsim::propagate(const int bounce, RNG& rng, s // copy geometry info into the sphoton struct ctx.p.set_prd(boundary, identity, cosTheta, iindex ); // HMM: lposcost not passed along - bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx ); + bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx, ctx.current_matline); // F7: pass carried matline for sibling-pair override if (ctx.current_material_index == 0u) { diff --git a/sysrap/sctx.h b/sysrap/sctx.h index 431975c54..ce4124357 100644 --- a/sysrap/sctx.h +++ b/sysrap/sctx.h @@ -79,6 +79,7 @@ struct sctx sstate s ; float current_group_velocity = 0.f; unsigned current_material_index = 0u; + unsigned current_matline = 0u; // F7 sibling-pair: photon's authoritative current-medium matline, seeded at birth, updated on transmit/reflect #ifndef PRODUCTION srec rec ; From 756f4cde2ba9c1104d08287af817f4e4eb43290a Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Sat, 30 May 2026 11:19:11 -0400 Subject: [PATCH 2/9] feat(CSGOptiX): index-contrast-gated exit-bias for CSG sibling-coincident faces At a face shared by two sibling primitives (segmented radiator bars + glue, multi-layer CSG lens) OptiX sorts the two coplanar intersections inconsistently; the "exiting" hit (cosI>0) carries the parent medium as m2, wrong for a sibling crossing. Bias the t reported to OptiX's closest-hit sort by +10 um on the exit side so the entering sibling wins. The real unbiased distance is preserved in the PRD, so photon position advancement is unperturbed. Gate the bias per boundary by index contrast |n_m1 - n_m2| > 0.1 (max over wavelength), computed at translation in CSGOptiX::initSimulate from the bnd table and read in __intersection__is. High-contrast faces (e.g. quartz/air ~0.47) need the entering-sibling selection; low-contrast thin gaps (e.g. aerogel/air ~0.025) are harmed by it and gate off. No per-detector hardcode; nullptr => always-on fallback. Stacks on the matline-carry m1 fix (carry source-medium matline across CSG sibling-pair faces): that fixes *which material* m1 is sampled; this fixes *which face* OptiX selects. The two are complementary - on a TIR light-guide either alone collapses the GPU hit count. Verified ePIC hpDIRC mu- 3 GeV: GPU/CPU 0.987, bit-identical to the ungated bias (all glued-sibling faces clear the 0.1 threshold). Low-index-contrast detectors are unaffected (gate stays off at their thin gaps). --- CSGOptiX/CSGOptiX.cc | 46 +++++++++++++++++++++++++++++++++++++++++++ CSGOptiX/CSGOptiX7.cu | 15 ++++++++++++-- CSGOptiX/Params.cc | 1 + CSGOptiX/Params.h | 1 + 4 files changed, 61 insertions(+), 2 deletions(-) diff --git a/CSGOptiX/CSGOptiX.cc b/CSGOptiX/CSGOptiX.cc index 49c942e04..38b47c237 100644 --- a/CSGOptiX/CSGOptiX.cc +++ b/CSGOptiX/CSGOptiX.cc @@ -595,6 +595,52 @@ void CSGOptiX::initSimulate() params->PropagateRefineDistance = SEventConfig::PropagateRefineDistance(); // approx distance beyond which to refine intersect with 2nd trace params->tmin = SEventConfig::PropagateEpsilon() ; // eg 0.1 0.05 to avoid self-intersection off boundaries + + // Per-boundary face-bias gate for the sibling-coincident t-bias (see __intersection__is). + // bit=1 only where the boundary's own index contrast |n_m1 - n_m2| (max over wavelength) + // exceeds a TIR-risk threshold: high-contrast faces (quartz/air ~0.47) NEED the entering- + // sibling bias; low-contrast thin gaps (aerogel/air ~0.025) are harmed and gate OFF. + // Computed from the bnd table; no per-detector hardcode. Allocated once, reused per launch. + { + const SSim* ssim = foundry ? foundry->getSim() : nullptr ; + const NP* bnd = ssim ? ssim->get_bnd() : nullptr ; + if( bnd != nullptr && bnd->shape.size() == 5 ) + { + const unsigned num_bnd = bnd->shape[0] ; + const unsigned num_wl = bnd->shape[3] ; + const double* bv = bnd->cvalues() ; // stree/standard bnd is bits(num_bnd, 0u) ; + for(unsigned b=0 ; b < num_bnd ; b++) + { + double dnmax = 0. ; + for(unsigned w=0 ; w < num_wl ; w++) + { + const double n_om = bv[ (((b*4u + LMAT_OMAT)*2u + 0u)*num_wl + w)*4u + 0u ] ; // RINDEX of omat + const double n_im = bv[ (((b*4u + LMAT_IMAT)*2u + 0u)*num_wl + w)*4u + 0u ] ; // RINDEX of imat + double dn = n_om - n_im ; if(dn < 0.) dn = -dn ; + if(dn > dnmax) dnmax = dn ; + } + bits[b] = (dnmax > DN_THRESHOLD) ? 1u : 0u ; + LOG(LEVEL) << "boundary_face_bias bnd " << b << " dnmax " << dnmax << " gate " << unsigned(bits[b]) ; + } + static unsigned char* s_face_bias = nullptr ; + static unsigned s_face_bias_num = 0u ; + if( s_face_bias == nullptr ){ CUDA_CHECK( cudaMalloc((void**)&s_face_bias, num_bnd) ); s_face_bias_num = num_bnd ; } + if( s_face_bias_num == num_bnd ) + { + CUDA_CHECK( cudaMemcpy(s_face_bias, bits.data(), num_bnd, cudaMemcpyHostToDevice) ); + params->boundary_face_bias = s_face_bias ; + } + LOG(LEVEL) << "boundary_face_bias gate computed, num_bnd " << num_bnd ; + } + else + { + LOG(LEVEL) << "boundary_face_bias: no bnd table (rank != 5) -> gate left null (t-bias always-on)" ; + } + } + params->tmax = 1000000.f ; params->max_time = SEventConfig::MaxTime() ; diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index 957a18cdc..a3bfc0662 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -910,8 +910,19 @@ extern "C" __global__ void __intersection__is() const unsigned boundary = node->boundary() ; // all CSGNode in the tree for one CSGPrim tree have same boundary const unsigned globalPrimIdx_boundary = (( globalPrimIdx & 0xffffu ) << 16 ) | ( boundary & 0xffffu ) ; + // Sibling-touching coincident-face fix: at a face shared by two sibling primitives OptiX + // sorts the two coplanar intersections inconsistently; the exiting hit (cosI>0) carries the + // parent medium as m2, wrong for a sibling crossing. Bias the t REPORTED to OptiX's closest- + // hit sort by +10 um on the exit side so the entering sibling wins; the real unbiased isect.w + // still goes to the PRD/attributes so position advancement is unperturbed. Gated per boundary + // by index contrast (params.boundary_face_bias from CSGOptiX::initSimulate): high-contrast + // faces (quartz/air ~0.47) need it; low-contrast thin gaps gate off; nullptr => always-on. + const float cosI = dot(ray_direction, make_float3(isect.x, isect.y, isect.z)); + const bool face_bias_on = (params.boundary_face_bias == nullptr) || (params.boundary_face_bias[boundary] != 0u) ; + const float t_report = (cosI > 0.f && face_bias_on) ? isect.w + 1.e-2f : isect.w ; + #ifdef WITH_PRD - if(optixReportIntersection( isect.w, hitKind)) + if(optixReportIntersection( t_report, hitKind)) { quad2* prd = SOPTIX_getPRD(); // access prd addr from RG program prd->q0.f = isect ; // .w:distance and .xyz:normal which starts as the local frame one @@ -927,7 +938,7 @@ extern "C" __global__ void __intersection__is() a3 = __float_as_uint( isect.w ) ; a4 = globalPrimIdx_boundary ; a5 = __float_as_uint( lposcost ); - optixReportIntersection( isect.w, hitKind, a0, a1, a2, a3, a4, a5 ); + optixReportIntersection( t_report, hitKind, a0, a1, a2, a3, a4, a5 ); // IS:optixReportIntersection writes the attributes that can be read in CH and AH programs // max 8 attribute registers, see PIP::PIP, communicate to __closesthit__ch diff --git a/CSGOptiX/Params.cc b/CSGOptiX/Params.cc index cc3a8b7c5..df23bb9eb 100644 --- a/CSGOptiX/Params.cc +++ b/CSGOptiX/Params.cc @@ -147,6 +147,7 @@ Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) origin_y(0), tmin(0.f), tmin0(0.f), + boundary_face_bias(nullptr), PropagateEpsilon0Mask(0u), tmax(0.f), vizmask(0xff), diff --git a/CSGOptiX/Params.h b/CSGOptiX/Params.h index 54a93faea..14f313404 100644 --- a/CSGOptiX/Params.h +++ b/CSGOptiX/Params.h @@ -64,6 +64,7 @@ struct Params float tmin ; float tmin0 ; + const unsigned char* boundary_face_bias ; // per-boundary gate (1=apply sibling-coincident t-bias) keyed on |n_m1-n_m2|>0.1; nullptr => always-on unsigned PropagateEpsilon0Mask ; // default from SEventConfig TO,CK,SI,SC,RE float PropagateRefineDistance ; bool PropagateRefine ; From d9502a7b14ca25a6f28fa0241076d8ff8f27733a Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 17 Jul 2026 11:06:51 -0400 Subject: [PATCH 3/9] fix(qudarap): only update carried matline on transmit Reflection never changes the photon's medium, so recomputing the carry from this boundary's own m1 line on the reflect branch could clobber a sibling-pair override from qbnd::fill_state with the parent medium's line (e.g. after a Fresnel/TIR reflection at a sibling-coincident face), poisoning m1 lookups at every subsequent boundary. Leave the carry untouched on reflect; transmit keeps the entered-medium (m2) line as before. --- qudarap/qsim.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/qudarap/qsim.h b/qudarap/qsim.h index c25e06518..53446601e 100644 --- a/qudarap/qsim.h +++ b/qudarap/qsim.h @@ -1265,13 +1265,15 @@ inline QSIM_METHOD FlowAction qsim::propagate_at_boundary(unsigned& flag, RNG& r // boundary knows the photon's actual medium even when that boundary doesn't // reference it (the sibling-pair case). _c1 = -dot(p.mom,normal) here, so // _c1 < 0 maps to cosTheta > 0 in qbnd::fill_state (m1=IMAT, m2=OMAT). + // Only transmit changes medium: on reflect the carry is left untouched, + // preserving any sibling override fill_state applied (recomputing m1 from + // this boundary's own lines would clobber it with the wrong medium). + if (!reflect) { const unsigned bnd_idx = ctx.prd->boundary(); const unsigned imat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + IMAT; const unsigned omat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + OMAT; - const unsigned this_m1_line = (_c1 < 0.f) ? imat_line : omat_line; - const unsigned this_m2_line = (_c1 < 0.f) ? omat_line : imat_line; - ctx.current_matline = reflect ? this_m1_line : this_m2_line; + ctx.current_matline = (_c1 < 0.f) ? omat_line : imat_line; // the m2 (entered-medium) line } #if !defined(PRODUCTION) && defined(DEBUG_TAG) From adad7055466c4589f011a168251734789a89af27 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Sat, 30 May 2026 15:57:56 +0000 Subject: [PATCH 4/9] style: clang-format (v20) PR-changed lines for cpp-linter (lines-changed-only) --- CSGOptiX/CSGOptiX.cc | 79 ++++++++++++++++++++----------------------- CSGOptiX/CSGOptiX7.cu | 8 ++--- CSGOptiX/Params.cc | 6 ++-- CSGOptiX/Params.h | 12 +++---- qudarap/qbnd.h | 35 +++++++++---------- qudarap/qsim.h | 4 +-- sysrap/sctx.h | 6 ++-- 7 files changed, 70 insertions(+), 80 deletions(-) diff --git a/CSGOptiX/CSGOptiX.cc b/CSGOptiX/CSGOptiX.cc index 38b47c237..9a5f5be6d 100644 --- a/CSGOptiX/CSGOptiX.cc +++ b/CSGOptiX/CSGOptiX.cc @@ -592,9 +592,9 @@ void CSGOptiX::initSimulate() params->PropagateEpsilon0Mask = SEventConfig::PropagateEpsilon0Mask(); // eg from CK|SI|TO|SC|RE params->PropagateRefine = SEventConfig::PropagateRefine(); - params->PropagateRefineDistance = SEventConfig::PropagateRefineDistance(); // approx distance beyond which to refine intersect with 2nd trace + params->PropagateRefineDistance = SEventConfig::PropagateRefineDistance(); // approx distance beyond which to refine intersect with 2nd trace - params->tmin = SEventConfig::PropagateEpsilon() ; // eg 0.1 0.05 to avoid self-intersection off boundaries + params->tmin = SEventConfig::PropagateEpsilon(); // eg 0.1 0.05 to avoid self-intersection off boundaries // Per-boundary face-bias gate for the sibling-coincident t-bias (see __intersection__is). // bit=1 only where the boundary's own index contrast |n_m1 - n_m2| (max over wavelength) @@ -602,61 +602,56 @@ void CSGOptiX::initSimulate() // sibling bias; low-contrast thin gaps (aerogel/air ~0.025) are harmed and gate OFF. // Computed from the bnd table; no per-detector hardcode. Allocated once, reused per launch. { - const SSim* ssim = foundry ? foundry->getSim() : nullptr ; - const NP* bnd = ssim ? ssim->get_bnd() : nullptr ; - if( bnd != nullptr && bnd->shape.size() == 5 ) + const SSim* ssim = foundry ? foundry->getSim() : nullptr; + const NP* bnd = ssim ? ssim->get_bnd() : nullptr; + if (bnd != nullptr && bnd->shape.size() == 5) { - const unsigned num_bnd = bnd->shape[0] ; - const unsigned num_wl = bnd->shape[3] ; - const double* bv = bnd->cvalues() ; // stree/standard bnd is bits(num_bnd, 0u) ; - for(unsigned b=0 ; b < num_bnd ; b++) + const unsigned num_bnd = bnd->shape[0]; + const unsigned num_wl = bnd->shape[3]; + const double* bv = bnd->cvalues(); // stree/standard bnd is bits(num_bnd, 0u); + for (unsigned b = 0; b < num_bnd; b++) { - double dnmax = 0. ; - for(unsigned w=0 ; w < num_wl ; w++) + double dnmax = 0.; + for (unsigned w = 0; w < num_wl; w++) { - const double n_om = bv[ (((b*4u + LMAT_OMAT)*2u + 0u)*num_wl + w)*4u + 0u ] ; // RINDEX of omat - const double n_im = bv[ (((b*4u + LMAT_IMAT)*2u + 0u)*num_wl + w)*4u + 0u ] ; // RINDEX of imat - double dn = n_om - n_im ; if(dn < 0.) dn = -dn ; - if(dn > dnmax) dnmax = dn ; + const double n_om = bv[(((b * 4u + LMAT_OMAT) * 2u + 0u) * num_wl + w) * 4u + 0u]; // RINDEX of omat + const double n_im = bv[(((b * 4u + LMAT_IMAT) * 2u + 0u) * num_wl + w) * 4u + 0u]; // RINDEX of imat + double dn = n_om - n_im; + if (dn < 0.) + dn = -dn; + if (dn > dnmax) + dnmax = dn; } - bits[b] = (dnmax > DN_THRESHOLD) ? 1u : 0u ; - LOG(LEVEL) << "boundary_face_bias bnd " << b << " dnmax " << dnmax << " gate " << unsigned(bits[b]) ; + bits[b] = (dnmax > DN_THRESHOLD) ? 1u : 0u; + LOG(LEVEL) << "boundary_face_bias bnd " << b << " dnmax " << dnmax << " gate " << unsigned(bits[b]); } - static unsigned char* s_face_bias = nullptr ; - static unsigned s_face_bias_num = 0u ; - if( s_face_bias == nullptr ){ CUDA_CHECK( cudaMalloc((void**)&s_face_bias, num_bnd) ); s_face_bias_num = num_bnd ; } - if( s_face_bias_num == num_bnd ) + static unsigned char* s_face_bias = nullptr; + static unsigned s_face_bias_num = 0u; + if (s_face_bias == nullptr) { - CUDA_CHECK( cudaMemcpy(s_face_bias, bits.data(), num_bnd, cudaMemcpyHostToDevice) ); - params->boundary_face_bias = s_face_bias ; + CUDA_CHECK(cudaMalloc((void**)&s_face_bias, num_bnd)); + s_face_bias_num = num_bnd; } - LOG(LEVEL) << "boundary_face_bias gate computed, num_bnd " << num_bnd ; + if (s_face_bias_num == num_bnd) + { + CUDA_CHECK(cudaMemcpy(s_face_bias, bits.data(), num_bnd, cudaMemcpyHostToDevice)); + params->boundary_face_bias = s_face_bias; + } + LOG(LEVEL) << "boundary_face_bias gate computed, num_bnd " << num_bnd; } else { - LOG(LEVEL) << "boundary_face_bias: no bnd table (rank != 5) -> gate left null (t-bias always-on)" ; + LOG(LEVEL) << "boundary_face_bias: no bnd table (rank != 5) -> gate left null (t-bias always-on)"; } } - params->tmax = 1000000.f ; - params->max_time = SEventConfig::MaxTime() ; - - + params->tmax = 1000000.f; + params->max_time = SEventConfig::MaxTime(); } - - - - - - - - - - /** CSGOptiX::initRender -------------------------- diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index a3bfc0662..27a3b9294 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -918,11 +918,11 @@ extern "C" __global__ void __intersection__is() // by index contrast (params.boundary_face_bias from CSGOptiX::initSimulate): high-contrast // faces (quartz/air ~0.47) need it; low-contrast thin gaps gate off; nullptr => always-on. const float cosI = dot(ray_direction, make_float3(isect.x, isect.y, isect.z)); - const bool face_bias_on = (params.boundary_face_bias == nullptr) || (params.boundary_face_bias[boundary] != 0u) ; - const float t_report = (cosI > 0.f && face_bias_on) ? isect.w + 1.e-2f : isect.w ; + const bool face_bias_on = (params.boundary_face_bias == nullptr) || (params.boundary_face_bias[boundary] != 0u); + const float t_report = (cosI > 0.f && face_bias_on) ? isect.w + 1.e-2f : isect.w; #ifdef WITH_PRD - if(optixReportIntersection( t_report, hitKind)) + if (optixReportIntersection(t_report, hitKind)) { quad2* prd = SOPTIX_getPRD(); // access prd addr from RG program prd->q0.f = isect ; // .w:distance and .xyz:normal which starts as the local frame one @@ -938,7 +938,7 @@ extern "C" __global__ void __intersection__is() a3 = __float_as_uint( isect.w ) ; a4 = globalPrimIdx_boundary ; a5 = __float_as_uint( lposcost ); - optixReportIntersection( t_report, hitKind, a0, a1, a2, a3, a4, a5 ); + optixReportIntersection(t_report, hitKind, a0, a1, a2, a3, a4, a5); // IS:optixReportIntersection writes the attributes that can be read in CH and AH programs // max 8 attribute registers, see PIP::PIP, communicate to __closesthit__ch diff --git a/CSGOptiX/Params.cc b/CSGOptiX/Params.cc index df23bb9eb..fe026ef4d 100644 --- a/CSGOptiX/Params.cc +++ b/CSGOptiX/Params.cc @@ -123,12 +123,10 @@ std::string Params::detail() const << std::setw(20) << " evt " << std::setw(10) << evt << std::endl ; std::string s = ss.str(); - return s ; + return s; } - -Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) - : +Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) : raygenmode(SRG_RENDER), node(nullptr), plan(nullptr), diff --git a/CSGOptiX/Params.h b/CSGOptiX/Params.h index 14f313404..1bdf84d86 100644 --- a/CSGOptiX/Params.h +++ b/CSGOptiX/Params.h @@ -61,13 +61,13 @@ struct Params float3 V ; float3 W ; float3 WNORM ; - float tmin ; + float tmin; - float tmin0 ; - const unsigned char* boundary_face_bias ; // per-boundary gate (1=apply sibling-coincident t-bias) keyed on |n_m1-n_m2|>0.1; nullptr => always-on - unsigned PropagateEpsilon0Mask ; // default from SEventConfig TO,CK,SI,SC,RE - float PropagateRefineDistance ; - bool PropagateRefine ; + float tmin0; + const unsigned char* boundary_face_bias; // per-boundary gate (1=apply sibling-coincident t-bias) keyed on |n_m1-n_m2|>0.1; nullptr => always-on + unsigned PropagateEpsilon0Mask; // default from SEventConfig TO,CK,SI,SC,RE + float PropagateRefineDistance; + bool PropagateRefine; float tmax ; float4 ZPROJ ; diff --git a/qudarap/qbnd.h b/qudarap/qbnd.h index fe9637a85..3154af4b1 100644 --- a/qudarap/qbnd.h +++ b/qudarap/qbnd.h @@ -35,16 +35,13 @@ struct qbnd unsigned boundary_tex_MaterialLine_LS ; quad* optical ; -#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_TEXTURE) || defined(MOCK_CUDA) - QBND_METHOD float4 boundary_lookup( unsigned ix, unsigned iy ); - QBND_METHOD float4 boundary_lookup( float nm, unsigned line, unsigned k ); +#if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_TEXTURE) || defined(MOCK_CUDA) + QBND_METHOD float4 boundary_lookup(unsigned ix, unsigned iy); + QBND_METHOD float4 boundary_lookup(float nm, unsigned line, unsigned k); QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline = 0u); #endif - }; - - #if defined( MOCK_TEXTURE) || defined(MOCK_CUDA) #include "stexture.h" #endif @@ -196,11 +193,11 @@ s.optical.x inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline) { - const int line = boundary*_BOUNDARY_NUM_MATSUR ; // now that are not signing boundary use 0-based + const int line = boundary * _BOUNDARY_NUM_MATSUR; // now that are not signing boundary use 0-based int m1_line = cosTheta > 0.f ? line + IMAT : line + OMAT; - const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT ; - const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR ; + const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT; + const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR; // F7 sibling-pair material override. // GBndLib does not enumerate boundary entries for CSG sibling-pair faces @@ -219,21 +216,21 @@ inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wav m1_line = int(carried_matline); } - s.material1 = boundary_lookup( wavelength, m1_line, 0); // refractive_index, absorption_length, scattering_length, reemission_prob - s.m1group2 = boundary_lookup( wavelength, m1_line, 1); // x: material1 group_velocity, y: material2 group_velocity, z/w unused - s.material2 = boundary_lookup( wavelength, m2_line, 0); // refractive_index, (absorption_length, scattering_length, reemission_prob) only m2:refractive index actually used - s.surface = boundary_lookup( wavelength, su_line, 0); // detect, , absorb , (reflect_specular), reflect_diffuse [they add to 1. so one not used] + s.material1 = boundary_lookup(wavelength, m1_line, 0); // refractive_index, absorption_length, scattering_length, reemission_prob + s.m1group2 = boundary_lookup(wavelength, m1_line, 1); // x: material1 group_velocity, y: material2 group_velocity, z/w unused + s.material2 = boundary_lookup(wavelength, m2_line, 0); // refractive_index, (absorption_length, scattering_length, reemission_prob) only m2:refractive index actually used + s.surface = boundary_lookup(wavelength, su_line, 0); // detect, , absorb , (reflect_specular), reflect_diffuse [they add to 1. so one not used] s.set_material2_group_velocity(boundary_lookup(wavelength, m2_line, 1).x); - s.optical = optical[su_line].u ; // 1-based-surface-index-0-meaning-boundary/type/finish/value (type,finish,value not used currently) + s.optical = optical[su_line].u; // 1-based-surface-index-0-meaning-boundary/type/finish/value (type,finish,value not used currently) - s.index.x = optical[m1_line].u.x; // m1 index (1-based, see sstandard::make_optical) -- reflects override if any - s.index.y = optical[m2_line].u.x ; // m2 index (1-based, see sstandard::make_optical) - s.index.z = optical[su_line].u.x ; // su index (1-based, see sstandard::make_optical) - s.index.w = 0u ; // avoid undefined memory comparison issues + s.index.x = optical[m1_line].u.x; // m1 index (1-based, see sstandard::make_optical) -- reflects override if any + s.index.y = optical[m2_line].u.x; // m2 index (1-based, see sstandard::make_optical) + s.index.z = optical[su_line].u.x; // su index (1-based, see sstandard::make_optical) + s.index.w = 0u; // avoid undefined memory comparison issues #if !defined(PRODUCTION) && defined(DEBUG_PIDX) - if( idx == base_pidx ) + if (idx == base_pidx) { printf("//qbnd.fill_state idx %7lld boundary %d line %d wavelength %10.4f m1_line %d m2_line %d su_line %d s.optical.x %d \n", idx, boundary, line, wavelength, m1_line, m2_line, su_line, s.optical.x ); printf("//qbnd.fill_state idx %7lld boundary %d [s.index.x-1](m1_index) %d [s.index.y-1](m2_index) %d [s.index.z-1](su_index) %d \n", idx, boundary, s.index.x-1u, s.index.y-1u, s.index.z-1u ); diff --git a/qudarap/qsim.h b/qudarap/qsim.h index 53446601e..e2027307d 100644 --- a/qudarap/qsim.h +++ b/qudarap/qsim.h @@ -1277,7 +1277,7 @@ inline QSIM_METHOD FlowAction qsim::propagate_at_boundary(unsigned& flag, RNG& r } #if !defined(PRODUCTION) && defined(DEBUG_TAG) - if( flag == BOUNDARY_REFLECT ) + if (flag == BOUNDARY_REFLECT) { const float u_br_align_0 = curand_uniform(&rng) ; const float u_br_align_1 = curand_uniform(&rng) ; @@ -2162,7 +2162,7 @@ QSIM_FORCEINLINE_METHOD FlowAction qsim::propagate(const int bounce, RNG& rng, s #endif // copy geometry info into the sphoton struct - ctx.p.set_prd(boundary, identity, cosTheta, iindex ); // HMM: lposcost not passed along + ctx.p.set_prd(boundary, identity, cosTheta, iindex); // HMM: lposcost not passed along bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx, ctx.current_matline); // F7: pass carried matline for sibling-pair override diff --git a/sysrap/sctx.h b/sysrap/sctx.h index ce4124357..2dc78f62a 100644 --- a/sysrap/sctx.h +++ b/sysrap/sctx.h @@ -76,13 +76,13 @@ struct sctx #endif sphoton p ; - sstate s ; - float current_group_velocity = 0.f; + sstate s; + float current_group_velocity = 0.f; unsigned current_material_index = 0u; unsigned current_matline = 0u; // F7 sibling-pair: photon's authoritative current-medium matline, seeded at birth, updated on transmit/reflect #ifndef PRODUCTION - srec rec ; + srec rec; sseq seq ; stagr tagr ; quad4 aux ; From db69d08541a75518f76894f74ef288a7af200020 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 17 Jul 2026 11:36:15 -0400 Subject: [PATCH 5/9] fix(examples): point validate_drich.py at committed geometry The script referenced epic_drich_simple.xml / drich_1sector_simple.xml, which were never committed, and hardcoded /tmp/simphony as the script directory. Use the bundled single-sector compact via DRICH_1SECTOR_XML, matching benchmark_drich.py, and derive the script directory from __file__ so the example runs from any checkout location. Note: constructing the 1-sector geometry against epic >= 26.06 still fails until the bundled radiator XML gains the element that newer epic DRICH plugins require; that applies equally to benchmark_drich.py and is left to a separate change. --- dd4hepplugins/examples/validate_drich.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/dd4hepplugins/examples/validate_drich.py b/dd4hepplugins/examples/validate_drich.py index 687208cc3..677a49c8e 100644 --- a/dd4hepplugins/examples/validate_drich.py +++ b/dd4hepplugins/examples/validate_drich.py @@ -1,8 +1,8 @@ #!/usr/bin/env python3 -"""Run GPU and CPU with simplified dRICH geometry (no non-optical PDU components).""" +"""Run GPU and CPU with the bundled single-sector dRICH geometry.""" import math, os, sys, numpy as np -_SCRIPT_DIR = "/tmp/simphony/dd4hepplugins/examples" +_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) _GEOM_DIR = os.path.join(_SCRIPT_DIR, "geometry") sys.path.insert(0, _SCRIPT_DIR) @@ -17,9 +17,9 @@ def setup_kernel(mode): cppyy.include("G4OpticalParameters.hh") from cppyy.gbl import G4OpticalParameters - compact = os.path.join(_GEOM_DIR, "epic_drich_simple.xml") - if "DRICH_SIMPLE_XML" not in os.environ: - os.environ["DRICH_SIMPLE_XML"] = os.path.join(_GEOM_DIR, "drich_1sector_simple.xml") + compact = os.path.join(_GEOM_DIR, "epic_drich_1sector.xml") + if "DRICH_1SECTOR_XML" not in os.environ: + os.environ["DRICH_1SECTOR_XML"] = os.path.join(_GEOM_DIR, "drich_1sector.xml") if mode == "gpu": os.environ["OPTICKS_EVENT_MODE"] = "Minimal" From 57340fc40211308faadc6ddf5d644362b2500a6b Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Wed, 29 Jul 2026 22:34:38 +0000 Subject: [PATCH 6/9] fix(CSGOptiX): safely resolve coincident sibling faces Seed the carried material line only for Cerenkov and scintillation gensteps, whose q0.z field has that contract. Treat UINT_MAX as the unavailable sentinel so material line zero remains valid, and validate carried rows against the uploaded boundary table before property or optical-buffer lookup. Replace the fixed 10 um, refractive-index-gated exit bias with a simulation-only one-ULP tie-break. This lets an entering sibling win exact coplanar ties without crossing a representable physical gap, while preserving the unbiased intersection for position advancement. Removing the old gate also eliminates its static GPU allocation and per-boundary debug logging. Add regression coverage for genstep material-line ownership; valid line-zero and inner-material carries; sentinel, surface-row, and out-of-range fallbacks; and the updated Params layout. Validated with a Release CUDA/OptiX build of CSGOptiX, GenstepTest, qbnd_boundary_lookup_test, and ParamsTest. All three focused CTest cases pass. Refs: BNLNPPS/simphony#346 Refs: https://github.com/BNLNPPS/simphony/pull/346#discussion_r3389469269 Refs: https://github.com/BNLNPPS/simphony/pull/346#discussion_r3389469313 Refs: https://github.com/BNLNPPS/simphony/pull/346#discussion_r3389469340 Refs: https://github.com/BNLNPPS/simphony/pull/346#discussion_r3389469368 --- CSGOptiX/CSGOptiX.cc | 52 ---------------------- CSGOptiX/CSGOptiX7.cu | 28 +++++------- CSGOptiX/Params.cc | 1 - CSGOptiX/Params.h | 9 ++-- qudarap/qbnd.h | 14 +++--- qudarap/tests/qbnd_boundary_lookup_test.cc | 32 +++++++++++++ sysrap/OpticksGenstep.h | 30 +++++++++++++ sysrap/sctx.h | 2 +- sysrap/tests/CMakeLists.txt | 1 + sysrap/tests/OpticksGenstepTest.cc | 28 ++++++++++++ 10 files changed, 116 insertions(+), 81 deletions(-) create mode 100644 sysrap/tests/OpticksGenstepTest.cc diff --git a/CSGOptiX/CSGOptiX.cc b/CSGOptiX/CSGOptiX.cc index 9a5f5be6d..e95166a34 100644 --- a/CSGOptiX/CSGOptiX.cc +++ b/CSGOptiX/CSGOptiX.cc @@ -596,58 +596,6 @@ void CSGOptiX::initSimulate() params->tmin = SEventConfig::PropagateEpsilon(); // eg 0.1 0.05 to avoid self-intersection off boundaries - // Per-boundary face-bias gate for the sibling-coincident t-bias (see __intersection__is). - // bit=1 only where the boundary's own index contrast |n_m1 - n_m2| (max over wavelength) - // exceeds a TIR-risk threshold: high-contrast faces (quartz/air ~0.47) NEED the entering- - // sibling bias; low-contrast thin gaps (aerogel/air ~0.025) are harmed and gate OFF. - // Computed from the bnd table; no per-detector hardcode. Allocated once, reused per launch. - { - const SSim* ssim = foundry ? foundry->getSim() : nullptr; - const NP* bnd = ssim ? ssim->get_bnd() : nullptr; - if (bnd != nullptr && bnd->shape.size() == 5) - { - const unsigned num_bnd = bnd->shape[0]; - const unsigned num_wl = bnd->shape[3]; - const double* bv = bnd->cvalues(); // stree/standard bnd is bits(num_bnd, 0u); - for (unsigned b = 0; b < num_bnd; b++) - { - double dnmax = 0.; - for (unsigned w = 0; w < num_wl; w++) - { - const double n_om = bv[(((b * 4u + LMAT_OMAT) * 2u + 0u) * num_wl + w) * 4u + 0u]; // RINDEX of omat - const double n_im = bv[(((b * 4u + LMAT_IMAT) * 2u + 0u) * num_wl + w) * 4u + 0u]; // RINDEX of imat - double dn = n_om - n_im; - if (dn < 0.) - dn = -dn; - if (dn > dnmax) - dnmax = dn; - } - bits[b] = (dnmax > DN_THRESHOLD) ? 1u : 0u; - LOG(LEVEL) << "boundary_face_bias bnd " << b << " dnmax " << dnmax << " gate " << unsigned(bits[b]); - } - static unsigned char* s_face_bias = nullptr; - static unsigned s_face_bias_num = 0u; - if (s_face_bias == nullptr) - { - CUDA_CHECK(cudaMalloc((void**)&s_face_bias, num_bnd)); - s_face_bias_num = num_bnd; - } - if (s_face_bias_num == num_bnd) - { - CUDA_CHECK(cudaMemcpy(s_face_bias, bits.data(), num_bnd, cudaMemcpyHostToDevice)); - params->boundary_face_bias = s_face_bias; - } - LOG(LEVEL) << "boundary_face_bias gate computed, num_bnd " << num_bnd; - } - else - { - LOG(LEVEL) << "boundary_face_bias: no bnd table (rank != 5) -> gate left null (t-bias always-on)"; - } - } - params->tmax = 1000000.f; params->max_time = SEventConfig::MaxTime(); } diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index 27a3b9294..c567846d4 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -439,14 +439,10 @@ static __forceinline__ __device__ void simulate( const uint3& launch_idx, const #endif sim->generate_photon(ctx.p, rng, gs, photon_idx, genstep_idx ); - // F7: seed carried matline from the genstep's source-medium bnd row - // (gs.q0.u.z == scerenkov/sscint matline). Updated through - // propagate_at_boundary. Guard the -1 sentinel (0xFFFFFFFF) that - // SEvt::setGenstep emits for bad_ck/unmapped mtindex. - { - const unsigned gm = gs.q0.u.z; - ctx.current_matline = (gm == 0xFFFFFFFFu) ? 0u : gm; - } + // Only Cerenkov and scintillation gensteps store a source-material line in + // q0.u.z. TORCH/FRAME use that slot for a genstep id, so leave the carry + // invalid for all genstep types without the material-line contract. + ctx.current_matline = OpticksGenstep_UsesMaterialLine(gs.q0.u.x) ? gs.q0.u.z : 0xFFFFFFFFu; FlowAction command = FlowAction::Start; int bounce = 0 ; @@ -910,16 +906,14 @@ extern "C" __global__ void __intersection__is() const unsigned boundary = node->boundary() ; // all CSGNode in the tree for one CSGPrim tree have same boundary const unsigned globalPrimIdx_boundary = (( globalPrimIdx & 0xffffu ) << 16 ) | ( boundary & 0xffffu ) ; - // Sibling-touching coincident-face fix: at a face shared by two sibling primitives OptiX - // sorts the two coplanar intersections inconsistently; the exiting hit (cosI>0) carries the - // parent medium as m2, wrong for a sibling crossing. Bias the t REPORTED to OptiX's closest- - // hit sort by +10 um on the exit side so the entering sibling wins; the real unbiased isect.w - // still goes to the PRD/attributes so position advancement is unperturbed. Gated per boundary - // by index contrast (params.boundary_face_bias from CSGOptiX::initSimulate): high-contrast - // faces (quartz/air ~0.47) need it; low-contrast thin gaps gate off; nullptr => always-on. + // Coincident sibling faces can produce equal entering and exiting hits. + // Move only the reported exiting distance to the next representable float + // during photon simulation so the entering sibling wins an exact tie. + // This one-ULP ordering change cannot overstep a representable physical + // gap. Keep the unbiased isect.w in the PRD for position advancement. const float cosI = dot(ray_direction, make_float3(isect.x, isect.y, isect.z)); - const bool face_bias_on = (params.boundary_face_bias == nullptr) || (params.boundary_face_bias[boundary] != 0u); - const float t_report = (cosI > 0.f && face_bias_on) ? isect.w + 1.e-2f : isect.w; + const bool break_exit_tie = params.raygenmode == SRG_SIMULATE && cosI > 0.f; + const float t_report = break_exit_tie ? nextafterf(isect.w, CUDART_INF_F) : isect.w; #ifdef WITH_PRD if (optixReportIntersection(t_report, hitKind)) diff --git a/CSGOptiX/Params.cc b/CSGOptiX/Params.cc index fe026ef4d..3d388a2fe 100644 --- a/CSGOptiX/Params.cc +++ b/CSGOptiX/Params.cc @@ -145,7 +145,6 @@ Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) origin_y(0), tmin(0.f), tmin0(0.f), - boundary_face_bias(nullptr), PropagateEpsilon0Mask(0u), tmax(0.f), vizmask(0xff), diff --git a/CSGOptiX/Params.h b/CSGOptiX/Params.h index 1bdf84d86..7b927e1fa 100644 --- a/CSGOptiX/Params.h +++ b/CSGOptiX/Params.h @@ -63,11 +63,10 @@ struct Params float3 WNORM ; float tmin; - float tmin0; - const unsigned char* boundary_face_bias; // per-boundary gate (1=apply sibling-coincident t-bias) keyed on |n_m1-n_m2|>0.1; nullptr => always-on - unsigned PropagateEpsilon0Mask; // default from SEventConfig TO,CK,SI,SC,RE - float PropagateRefineDistance; - bool PropagateRefine; + float tmin0; + unsigned PropagateEpsilon0Mask; // default from SEventConfig TO,CK,SI,SC,RE + float PropagateRefineDistance; + bool PropagateRefine; float tmax ; float4 ZPROJ ; diff --git a/qudarap/qbnd.h b/qudarap/qbnd.h index 3154af4b1..a554c454b 100644 --- a/qudarap/qbnd.h +++ b/qudarap/qbnd.h @@ -38,7 +38,7 @@ struct qbnd #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_TEXTURE) || defined(MOCK_CUDA) QBND_METHOD float4 boundary_lookup(unsigned ix, unsigned iy); QBND_METHOD float4 boundary_lookup(float nm, unsigned line, unsigned k); - QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline = 0u); + QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline = 0xFFFFFFFFu); #endif }; @@ -207,11 +207,15 @@ inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wav // genstep, updated each propagate_at_boundary); when it is valid we use it // directly as m1 so absorption/scattering sample the true current medium. // (material2 is left as the boundary's other side.) - // Guards: 0 = "no carry" (legacy callers); 0xFFFFFFFF = the lookup_mtline=-1 - // sentinel from SEvt::setGenstep, not a valid table row; and the carried - // slot must be a material row (OMAT/IMAT), not a surface row (OSUR/ISUR). + // UINT_MAX means that no source material is available. Line zero remains a + // valid OMAT row. Also reject surface rows and values beyond the boundary + // texture before using the carry for texture or optical-buffer indexing. + const unsigned num_matline = boundary_meta->q0.u.y / _BOUNDARY_NUM_FLOAT4; const unsigned matline_slot = carried_matline & 0x3u; - if (carried_matline != 0u && carried_matline != 0xFFFFFFFFu && (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT))) + const bool carried_matline_valid = + carried_matline < num_matline && + (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT)); + if (carried_matline_valid) { m1_line = int(carried_matline); } diff --git a/qudarap/tests/qbnd_boundary_lookup_test.cc b/qudarap/tests/qbnd_boundary_lookup_test.cc index 77cb99deb..e0e134854 100644 --- a/qudarap/tests/qbnd_boundary_lookup_test.cc +++ b/qudarap/tests/qbnd_boundary_lookup_test.cc @@ -102,6 +102,10 @@ int main() qbnd qb; qb.boundary_tex = texObj; qb.boundary_meta = &meta; + quad optical[ni * nj] = {}; + for (unsigned line = 0; line < ni * nj; line++) + optical[line].u.x = 1000u + line; + qb.optical = optical; printf("qbnd_boundary_lookup_test\n"); @@ -160,6 +164,34 @@ int main() check("row0 clamp above: .x == 30 (bin 3)", feq(p.x, 30.f)); } + // --- Test 4: Carried material-line validation --- + // Boundary 1 entered from outside normally uses line 4 (OMAT). + // A valid carry overrides only material1; invalid, surface, and out-of-range + // rows must fall back to line 4 without indexing outside the buffers. + { + const unsigned boundary = 1u; + const float cosTheta = -1.f; + sstate s = {}; + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 0u); + check("carried line 0 remains a valid OMAT override", feq(s.material1.x, 0.f)); + check("carried line 0 selects optical material index", s.index.x == 1000u); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 7u); + check("carried line 7 remains a valid IMAT override", feq(s.material1.x, 1400.f)); + check("carried line 7 selects optical material index", s.index.x == 1007u); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 0xFFFFFFFFu); + check("invalid sentinel falls back to boundary OMAT", feq(s.material1.x, 800.f)); + check("invalid sentinel keeps boundary material index", s.index.x == 1004u); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 5u); + check("surface row carry falls back to boundary OMAT", feq(s.material1.x, 800.f)); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, ni * nj); + check("out-of-range carry falls back to boundary OMAT", feq(s.material1.x, 800.f)); + } + printf("qbnd_boundary_lookup_test: %s (%d failure%s)\n", g_fail == 0 ? "PASS" : "FAIL", g_fail, g_fail == 1 ? "" : "s"); return g_fail == 0 ? 0 : 1; diff --git a/sysrap/OpticksGenstep.h b/sysrap/OpticksGenstep.h index 9f61e8b4b..5ba589286 100644 --- a/sysrap/OpticksGenstep.h +++ b/sysrap/OpticksGenstep.h @@ -38,6 +38,36 @@ enum OpticksGenstep_NumType = 21 }; +#if defined(__CUDACC__) || defined(__CUDABE__) +#define OPTICKS_GENSTEP_METHOD __forceinline__ __device__ +#else +#define OPTICKS_GENSTEP_METHOD inline +#endif + +/** + * Only Cerenkov and scintillation gensteps use q0.u.z for the source material + * line. Other gensteps reuse the slot for unrelated data, such as the packed + * genstep id used by TORCH and FRAME. + */ +OPTICKS_GENSTEP_METHOD bool OpticksGenstep_UsesMaterialLine(unsigned gentype) +{ + switch (gentype) + { + case OpticksGenstep_G4Cerenkov_1042: + case OpticksGenstep_G4Scintillation_1042: + case OpticksGenstep_DsG4Cerenkov_r3971: + case OpticksGenstep_DsG4Scintillation_r3971: + case OpticksGenstep_DsG4Scintillation_r4695: + case OpticksGenstep_CERENKOV: + case OpticksGenstep_SCINTILLATION: + case OpticksGenstep_G4Cerenkov_modified: + return true; + default: + return false; + } +} + +#undef OPTICKS_GENSTEP_METHOD #if defined(__CUDACC__) || defined(__CUDABE__) #else diff --git a/sysrap/sctx.h b/sysrap/sctx.h index 2dc78f62a..88ef6650b 100644 --- a/sysrap/sctx.h +++ b/sysrap/sctx.h @@ -79,7 +79,7 @@ struct sctx sstate s; float current_group_velocity = 0.f; unsigned current_material_index = 0u; - unsigned current_matline = 0u; // F7 sibling-pair: photon's authoritative current-medium matline, seeded at birth, updated on transmit/reflect + unsigned current_matline = 0xFFFFFFFFu; // sibling-pair current-medium line; UINT_MAX means unavailable #ifndef PRODUCTION srec rec; diff --git a/sysrap/tests/CMakeLists.txt b/sysrap/tests/CMakeLists.txt index 6c4ae88c5..4ae07743f 100644 --- a/sysrap/tests/CMakeLists.txt +++ b/sysrap/tests/CMakeLists.txt @@ -85,6 +85,7 @@ set(TEST_SOURCES SFrameGenstep_MakeCenterExtentGensteps_Test.cc SEventTest.cc + OpticksGenstepTest.cc SThetaCutTest.cc SPhiCutTest.cc scanvasTest.cc diff --git a/sysrap/tests/OpticksGenstepTest.cc b/sysrap/tests/OpticksGenstepTest.cc new file mode 100644 index 000000000..e5481aa10 --- /dev/null +++ b/sysrap/tests/OpticksGenstepTest.cc @@ -0,0 +1,28 @@ +#include + +#include "OpticksGenstep.h" + +static int fail = 0; + +static void check(const char* label, bool expected, bool actual) +{ + const bool pass = expected == actual; + std::printf(" %-44s %s\n", label, pass ? "PASS" : "*** FAIL ***"); + if (!pass) + fail++; +} + +int main() +{ + check("Cerenkov carries a material line", true, OpticksGenstep_UsesMaterialLine(OpticksGenstep_CERENKOV)); + check("scintillation carries a material line", true, OpticksGenstep_UsesMaterialLine(OpticksGenstep_SCINTILLATION)); + check("modified G4 Cerenkov carries a material line", true, OpticksGenstep_UsesMaterialLine(OpticksGenstep_G4Cerenkov_modified)); + check("TORCH does not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_TORCH)); + check("FRAME does not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_FRAME)); + check("input photons do not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_INPUT_PHOTON)); + check("carrier gensteps do not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_CARRIER)); + + std::printf("OpticksGenstepTest: %s (%d failure%s)\n", + fail == 0 ? "PASS" : "FAIL", fail, fail == 1 ? "" : "s"); + return fail == 0 ? 0 : 1; +} From 7796d9c7334c8415b75eeb3a66a934abe0581518 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 30 Jul 2026 17:09:22 +0000 Subject: [PATCH 7/9] fix(sysrap): remove retired material-line gensteps Align OpticksGenstep_UsesMaterialLine with the current genstep protocol after BNLNPPS/simphony#429 retired the five legacy suffixed identifiers. Only the generic Cerenkov, scintillation, and modified G4 Cerenkov producers retain the q0.z material-line contract. Cover numeric protocol slots 1 through 5 as retired non-material-line gensteps so a future reintroduction cannot silently restore the stale classification. Validated by rebuilding CSGOptiX7_ptx and OpticksGenstepTest and running SysRapTest.OpticksGenstepTest plus SysRapTest.GenstepTest. Refs: BNLNPPS/simphony#346 Refs: BNLNPPS/simphony#429 --- sysrap/OpticksGenstep.h | 5 ----- sysrap/tests/OpticksGenstepTest.cc | 2 ++ 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/sysrap/OpticksGenstep.h b/sysrap/OpticksGenstep.h index 5ba589286..314e37b9a 100644 --- a/sysrap/OpticksGenstep.h +++ b/sysrap/OpticksGenstep.h @@ -53,11 +53,6 @@ OPTICKS_GENSTEP_METHOD bool OpticksGenstep_UsesMaterialLine(unsigned gentype) { switch (gentype) { - case OpticksGenstep_G4Cerenkov_1042: - case OpticksGenstep_G4Scintillation_1042: - case OpticksGenstep_DsG4Cerenkov_r3971: - case OpticksGenstep_DsG4Scintillation_r3971: - case OpticksGenstep_DsG4Scintillation_r4695: case OpticksGenstep_CERENKOV: case OpticksGenstep_SCINTILLATION: case OpticksGenstep_G4Cerenkov_modified: diff --git a/sysrap/tests/OpticksGenstepTest.cc b/sysrap/tests/OpticksGenstepTest.cc index e5481aa10..c84032be4 100644 --- a/sysrap/tests/OpticksGenstepTest.cc +++ b/sysrap/tests/OpticksGenstepTest.cc @@ -22,6 +22,8 @@ int main() check("input photons do not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_INPUT_PHOTON)); check("carrier gensteps do not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_CARRIER)); + for (unsigned type = 1; type <= 5; type++) + check("retired genstep does not carry a material line", false, OpticksGenstep_UsesMaterialLine(type)); std::printf("OpticksGenstepTest: %s (%d failure%s)\n", fail == 0 ? "PASS" : "FAIL", fail, fail == 1 ? "" : "s"); return fail == 0 ? 0 : 1; From 72ddeaacbf9c124190102339e0e9884a70d341b2 Mon Sep 17 00:00:00 2001 From: Dmitri Smirnov Date: Thu, 30 Jul 2026 13:37:14 -0400 Subject: [PATCH 8/9] fix --- CSGOptiX/CSGOptiX7.cu | 11 ++++++--- dd4hepplugins/examples/validate_drich.py | 10 ++++---- qudarap/qbnd.h | 29 ++++++++++++----------- qudarap/qsim.h | 22 ++++++++++------- sysrap/OpticksGenstep.h | 25 -------------------- sysrap/tests/CMakeLists.txt | 1 - sysrap/tests/OpticksGenstepTest.cc | 30 ------------------------ 7 files changed, 41 insertions(+), 87 deletions(-) delete mode 100644 sysrap/tests/OpticksGenstepTest.cc diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index c567846d4..0a9002140 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -440,9 +440,14 @@ static __forceinline__ __device__ void simulate( const uint3& launch_idx, const sim->generate_photon(ctx.p, rng, gs, photon_idx, genstep_idx ); // Only Cerenkov and scintillation gensteps store a source-material line in - // q0.u.z. TORCH/FRAME use that slot for a genstep id, so leave the carry - // invalid for all genstep types without the material-line contract. - ctx.current_matline = OpticksGenstep_UsesMaterialLine(gs.q0.u.x) ? gs.q0.u.z : 0xFFFFFFFFu; + // q0.u.z. TORCH and FRAME use that slot for a genstep id, so leave the + // carried line invalid for all other genstep types. + const unsigned gentype = gs.q0.u.x; + const bool carries_matline = + gentype == OpticksGenstep_CERENKOV || + gentype == OpticksGenstep_SCINTILLATION || + gentype == OpticksGenstep_G4Cerenkov_modified; + ctx.current_matline = carries_matline ? gs.q0.u.z : 0xFFFFFFFFu; FlowAction command = FlowAction::Start; int bounce = 0 ; diff --git a/dd4hepplugins/examples/validate_drich.py b/dd4hepplugins/examples/validate_drich.py index 677a49c8e..687208cc3 100644 --- a/dd4hepplugins/examples/validate_drich.py +++ b/dd4hepplugins/examples/validate_drich.py @@ -1,8 +1,8 @@ #!/usr/bin/env python3 -"""Run GPU and CPU with the bundled single-sector dRICH geometry.""" +"""Run GPU and CPU with simplified dRICH geometry (no non-optical PDU components).""" import math, os, sys, numpy as np -_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) +_SCRIPT_DIR = "/tmp/simphony/dd4hepplugins/examples" _GEOM_DIR = os.path.join(_SCRIPT_DIR, "geometry") sys.path.insert(0, _SCRIPT_DIR) @@ -17,9 +17,9 @@ def setup_kernel(mode): cppyy.include("G4OpticalParameters.hh") from cppyy.gbl import G4OpticalParameters - compact = os.path.join(_GEOM_DIR, "epic_drich_1sector.xml") - if "DRICH_1SECTOR_XML" not in os.environ: - os.environ["DRICH_1SECTOR_XML"] = os.path.join(_GEOM_DIR, "drich_1sector.xml") + compact = os.path.join(_GEOM_DIR, "epic_drich_simple.xml") + if "DRICH_SIMPLE_XML" not in os.environ: + os.environ["DRICH_SIMPLE_XML"] = os.path.join(_GEOM_DIR, "drich_1sector_simple.xml") if mode == "gpu": os.environ["OPTICKS_EVENT_MODE"] = "Minimal" diff --git a/qudarap/qbnd.h b/qudarap/qbnd.h index a554c454b..8dd2597fa 100644 --- a/qudarap/qbnd.h +++ b/qudarap/qbnd.h @@ -199,22 +199,23 @@ inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wav const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT; const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR; - // F7 sibling-pair material override. - // GBndLib does not enumerate boundary entries for CSG sibling-pair faces - // (two daughters of one parent sharing a face), so OptiX may pick a boundary - // whose m1/m2 do not reference the photon's actual medium -> wrong material1. - // qsim::propagate passes the photon's carried matline (seeded from the - // genstep, updated each propagate_at_boundary); when it is valid we use it - // directly as m1 so absorption/scattering sample the true current medium. - // (material2 is left as the boundary's other side.) - // UINT_MAX means that no source material is available. Line zero remains a - // valid OMAT row. Also reject surface rows and values beyond the boundary - // texture before using the carry for texture or optical-buffer indexing. + // A boundary normally describes a daughter volume relative to its parent. + // When two daughter volumes touch, however, a photon can move directly from + // one sibling to the other. The boundary selected by OptiX may then name the + // common parent instead of the material the photon is actually leaving. + // + // qsim carries the material-table row for the photon's current volume across + // interactions. Use that row as material1 when it is valid so absorption and + // scattering are sampled from the correct medium. Material2 and the surface + // still come from the selected boundary. + // + // UINT_MAX means that no current-material row is available. Line zero is a + // valid material row. Surface rows and rows beyond the uploaded boundary + // table must not be used for property or optical-buffer lookups. const unsigned num_matline = boundary_meta->q0.u.y / _BOUNDARY_NUM_FLOAT4; const unsigned matline_slot = carried_matline & 0x3u; - const bool carried_matline_valid = - carried_matline < num_matline && - (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT)); + const bool carried_matline_valid = carried_matline < num_matline && (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT)); + if (carried_matline_valid) { m1_line = int(carried_matline); diff --git a/qudarap/qsim.h b/qudarap/qsim.h index e2027307d..68de1fc34 100644 --- a/qudarap/qsim.h +++ b/qudarap/qsim.h @@ -1261,19 +1261,21 @@ inline QSIM_METHOD FlowAction qsim::propagate_at_boundary(unsigned& flag, RNG& r ctx.current_material_index = reflect ? s.index.x : s.index.y; ctx.current_group_velocity = reflect ? s.material1_group_velocity() : s.material2_group_velocity(); - // F7 SIBLING-PAIR carry: update current_matline so propagate() at the NEXT - // boundary knows the photon's actual medium even when that boundary doesn't - // reference it (the sibling-pair case). _c1 = -dot(p.mom,normal) here, so - // _c1 < 0 maps to cosTheta > 0 in qbnd::fill_state (m1=IMAT, m2=OMAT). - // Only transmit changes medium: on reflect the carry is left untouched, - // preserving any sibling override fill_state applied (recomputing m1 from - // this boundary's own lines would clobber it with the wrong medium). + // Remember the material entered by a transmitted photon so the next + // boundary lookup knows the photon's actual current medium. This matters + // when touching sibling volumes share a face because the selected boundary + // may describe one sibling relative to their common parent and omit the + // material the photon is leaving. + // + // Reflection does not change the current medium, so preserve the existing + // material row in that case. For transmission, _c1 identifies which side of + // the selected boundary is the material the photon enters. if (!reflect) { const unsigned bnd_idx = ctx.prd->boundary(); const unsigned imat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + IMAT; const unsigned omat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + OMAT; - ctx.current_matline = (_c1 < 0.f) ? omat_line : imat_line; // the m2 (entered-medium) line + ctx.current_matline = (_c1 < 0.f) ? omat_line : imat_line; // material entered by this transmission } #if !defined(PRODUCTION) && defined(DEBUG_TAG) @@ -2164,7 +2166,9 @@ QSIM_FORCEINLINE_METHOD FlowAction qsim::propagate(const int bounce, RNG& rng, s // copy geometry info into the sphoton struct ctx.p.set_prd(boundary, identity, cosTheta, iindex); // HMM: lposcost not passed along - bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx, ctx.current_matline); // F7: pass carried matline for sibling-pair override + // Supply the current-material row so fill_state can correct material1 when + // the selected boundary does not describe the photon's actual current medium. + bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx, ctx.current_matline); if (ctx.current_material_index == 0u) { diff --git a/sysrap/OpticksGenstep.h b/sysrap/OpticksGenstep.h index 314e37b9a..9f61e8b4b 100644 --- a/sysrap/OpticksGenstep.h +++ b/sysrap/OpticksGenstep.h @@ -38,31 +38,6 @@ enum OpticksGenstep_NumType = 21 }; -#if defined(__CUDACC__) || defined(__CUDABE__) -#define OPTICKS_GENSTEP_METHOD __forceinline__ __device__ -#else -#define OPTICKS_GENSTEP_METHOD inline -#endif - -/** - * Only Cerenkov and scintillation gensteps use q0.u.z for the source material - * line. Other gensteps reuse the slot for unrelated data, such as the packed - * genstep id used by TORCH and FRAME. - */ -OPTICKS_GENSTEP_METHOD bool OpticksGenstep_UsesMaterialLine(unsigned gentype) -{ - switch (gentype) - { - case OpticksGenstep_CERENKOV: - case OpticksGenstep_SCINTILLATION: - case OpticksGenstep_G4Cerenkov_modified: - return true; - default: - return false; - } -} - -#undef OPTICKS_GENSTEP_METHOD #if defined(__CUDACC__) || defined(__CUDABE__) #else diff --git a/sysrap/tests/CMakeLists.txt b/sysrap/tests/CMakeLists.txt index 4ae07743f..6c4ae88c5 100644 --- a/sysrap/tests/CMakeLists.txt +++ b/sysrap/tests/CMakeLists.txt @@ -85,7 +85,6 @@ set(TEST_SOURCES SFrameGenstep_MakeCenterExtentGensteps_Test.cc SEventTest.cc - OpticksGenstepTest.cc SThetaCutTest.cc SPhiCutTest.cc scanvasTest.cc diff --git a/sysrap/tests/OpticksGenstepTest.cc b/sysrap/tests/OpticksGenstepTest.cc deleted file mode 100644 index c84032be4..000000000 --- a/sysrap/tests/OpticksGenstepTest.cc +++ /dev/null @@ -1,30 +0,0 @@ -#include - -#include "OpticksGenstep.h" - -static int fail = 0; - -static void check(const char* label, bool expected, bool actual) -{ - const bool pass = expected == actual; - std::printf(" %-44s %s\n", label, pass ? "PASS" : "*** FAIL ***"); - if (!pass) - fail++; -} - -int main() -{ - check("Cerenkov carries a material line", true, OpticksGenstep_UsesMaterialLine(OpticksGenstep_CERENKOV)); - check("scintillation carries a material line", true, OpticksGenstep_UsesMaterialLine(OpticksGenstep_SCINTILLATION)); - check("modified G4 Cerenkov carries a material line", true, OpticksGenstep_UsesMaterialLine(OpticksGenstep_G4Cerenkov_modified)); - check("TORCH does not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_TORCH)); - check("FRAME does not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_FRAME)); - check("input photons do not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_INPUT_PHOTON)); - check("carrier gensteps do not carry a material line", false, OpticksGenstep_UsesMaterialLine(OpticksGenstep_CARRIER)); - - for (unsigned type = 1; type <= 5; type++) - check("retired genstep does not carry a material line", false, OpticksGenstep_UsesMaterialLine(type)); - std::printf("OpticksGenstepTest: %s (%d failure%s)\n", - fail == 0 ? "PASS" : "FAIL", fail, fail == 1 ? "" : "s"); - return fail == 0 ? 0 : 1; -} From ab6c8b14a812208f671070ee49fba38855316b34 Mon Sep 17 00:00:00 2001 From: Dmitri Smirnov Date: Thu, 30 Jul 2026 19:24:10 +0000 Subject: [PATCH 9/9] test(CSGOptiX): verify sibling material propagation Exercise the carried material-line lifecycle directly through qsim, checking that transmission records the entered material while reflection preserves the current one. Share the CUDA-local genstep classification and reported-distance rules with an on-device test covering TORCH/Cerenkov handling and exact entering-versus-exiting intersection ties. Add a compact touching-sibling GDML fixture and simg4ox CPU/GPU hit comparison so the production propagation path must cross the shared face and retain the correct medium properties to reach the sensor. Refs: #346 --- CSGOptiX/CMakeLists.txt | 1 + CSGOptiX/CSGOptiX7.cu | 9 +- CSGOptiX/CSGOptiXHelpers.h | 30 ++++++ CSGOptiX/tests/CMakeLists.txt | 3 + CSGOptiX/tests/CSGOptiXHelpersTest.cu | 92 ++++++++++++++++++ config/sibling_pair.json | 34 +++++++ qudarap/tests/CMakeLists.txt | 5 + qudarap/tests/QSim_MaterialCarryTest.cc | 121 ++++++++++++++++++++++++ tests/CMakeLists.txt | 1 + tests/geom/sibling_pair.gdml | 97 +++++++++++++++++++ tests/test_simg4ox.sh | 7 ++ 11 files changed, 394 insertions(+), 6 deletions(-) create mode 100644 CSGOptiX/CSGOptiXHelpers.h create mode 100644 CSGOptiX/tests/CSGOptiXHelpersTest.cu create mode 100644 config/sibling_pair.json create mode 100644 qudarap/tests/QSim_MaterialCarryTest.cc create mode 100644 tests/geom/sibling_pair.gdml diff --git a/CSGOptiX/CMakeLists.txt b/CSGOptiX/CMakeLists.txt index d68ac98f7..4d0ea3540 100644 --- a/CSGOptiX/CMakeLists.txt +++ b/CSGOptiX/CMakeLists.txt @@ -22,6 +22,7 @@ list(APPEND SOURCES list(APPEND HEADERS CSGOptiX.h CSGOptiXService.h + CSGOptiXHelpers.h CSGOPTIX_API_EXPORT.hh CSGOPTIX_LOG.hh ) diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index 0a9002140..5c2d93b9e 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -77,6 +77,7 @@ __intersection__is #include "Binding.h" #include "Params.h" +#include "CSGOptiXHelpers.h" #ifdef WITH_PRD #include "scuda_pointer.h" @@ -443,10 +444,7 @@ static __forceinline__ __device__ void simulate( const uint3& launch_idx, const // q0.u.z. TORCH and FRAME use that slot for a genstep id, so leave the // carried line invalid for all other genstep types. const unsigned gentype = gs.q0.u.x; - const bool carries_matline = - gentype == OpticksGenstep_CERENKOV || - gentype == OpticksGenstep_SCINTILLATION || - gentype == OpticksGenstep_G4Cerenkov_modified; + const bool carries_matline = CSGOptiX7_GenstepCarriesMaterialLine(gentype); ctx.current_matline = carries_matline ? gs.q0.u.z : 0xFFFFFFFFu; FlowAction command = FlowAction::Start; @@ -917,8 +915,7 @@ extern "C" __global__ void __intersection__is() // This one-ULP ordering change cannot overstep a representable physical // gap. Keep the unbiased isect.w in the PRD for position advancement. const float cosI = dot(ray_direction, make_float3(isect.x, isect.y, isect.z)); - const bool break_exit_tie = params.raygenmode == SRG_SIMULATE && cosI > 0.f; - const float t_report = break_exit_tie ? nextafterf(isect.w, CUDART_INF_F) : isect.w; + const float t_report = CSGOptiX7_ReportedIntersectionDistance(isect.w, params.raygenmode, cosI); #ifdef WITH_PRD if (optixReportIntersection(t_report, hitKind)) diff --git a/CSGOptiX/CSGOptiXHelpers.h b/CSGOptiX/CSGOptiXHelpers.h new file mode 100644 index 000000000..cc943db7b --- /dev/null +++ b/CSGOptiX/CSGOptiXHelpers.h @@ -0,0 +1,30 @@ +#pragma once + +#include + +#include "OpticksGenstep.h" +#include "SRG.h" + +#if defined(__CUDACC__) || defined(__CUDABE__) +#define CSGOPTIX7_HELPER_METHOD __device__ __forceinline__ +#else +#define CSGOPTIX7_HELPER_METHOD inline +#endif + +CSGOPTIX7_HELPER_METHOD bool CSGOptiX7_GenstepCarriesMaterialLine(unsigned gentype) +{ + return gentype == OpticksGenstep_CERENKOV || + gentype == OpticksGenstep_SCINTILLATION || + gentype == OpticksGenstep_G4Cerenkov_modified; +} + +CSGOPTIX7_HELPER_METHOD float CSGOptiX7_ReportedIntersectionDistance( + float distance, + unsigned raygenmode, + float direction_dot_normal) +{ + const bool defer_exit = raygenmode == SRG_SIMULATE && direction_dot_normal > 0.f; + return defer_exit ? nextafterf(distance, INFINITY) : distance; +} + +#undef CSGOPTIX7_HELPER_METHOD diff --git a/CSGOptiX/tests/CMakeLists.txt b/CSGOptiX/tests/CMakeLists.txt index d5dbfe09b..3789957c1 100644 --- a/CSGOptiX/tests/CMakeLists.txt +++ b/CSGOptiX/tests/CMakeLists.txt @@ -10,6 +10,7 @@ or doing geometry specific things. set(SINGLE_SOURCES + CSGOptiXHelpersTest.cu CSGOptiXDescTest.cc CSGOptiXSimtraceTest.cc CSGOptiXSimulateTest.cc @@ -49,4 +50,6 @@ foreach(SRC ${TEST_SOURCES}) ) endforeach() +add_test(NAME ${name}.CSGOptiXHelpersTest COMMAND $) + set_tests_properties(${name}.CSGOptiXRenderTest PROPERTIES DISABLED TRUE) diff --git a/CSGOptiX/tests/CSGOptiXHelpersTest.cu b/CSGOptiX/tests/CSGOptiXHelpersTest.cu new file mode 100644 index 000000000..5e28e15ce --- /dev/null +++ b/CSGOptiX/tests/CSGOptiXHelpersTest.cu @@ -0,0 +1,92 @@ +#include +#include + +#include + +#include "CSGOptiXHelpers.h" + +namespace +{ + void require(bool condition, const char* expression, int line) + { + if (condition) return; + std::fprintf(stderr, "%s:%d requirement failed: %s\n", __FILE__, line, expression); + std::exit(EXIT_FAILURE); + } +} + +#define REQUIRE(expression) require((expression), #expression, __LINE__) + +struct CSGOptiXHelpersTestResult +{ + bool torch_carries_matline; + bool frame_carries_matline; + bool cerenkov_carries_matline; + bool scintillation_carries_matline; + bool modified_cerenkov_carries_matline; + float entering_distance; + float exiting_distance; + float non_simulation_exit_distance; +}; + +__global__ void CSGOptiXHelpersTest_kernel(CSGOptiXHelpersTestResult* result) +{ + result->torch_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_TORCH); + result->frame_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_FRAME); + result->cerenkov_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_CERENKOV); + result->scintillation_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_SCINTILLATION); + result->modified_cerenkov_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_G4Cerenkov_modified); + + const float coincident_distance = 10.f; + result->entering_distance = CSGOptiX7_ReportedIntersectionDistance( + coincident_distance, + SRG_SIMULATE, + -1.f); + result->exiting_distance = CSGOptiX7_ReportedIntersectionDistance( + coincident_distance, + SRG_SIMULATE, + 1.f); + result->non_simulation_exit_distance = CSGOptiX7_ReportedIntersectionDistance( + coincident_distance, + SRG_SIMTRACE, + 1.f); +} + +int main() +{ + CSGOptiXHelpersTestResult* device_result = nullptr; + REQUIRE(cudaMalloc(&device_result, sizeof(CSGOptiXHelpersTestResult)) == cudaSuccess); + + CSGOptiXHelpersTest_kernel<<<1, 1>>>(device_result); + REQUIRE(cudaGetLastError() == cudaSuccess); + REQUIRE(cudaDeviceSynchronize() == cudaSuccess); + + CSGOptiXHelpersTestResult result = {}; + REQUIRE(cudaMemcpy( + &result, + device_result, + sizeof(CSGOptiXHelpersTestResult), + cudaMemcpyDeviceToHost) == cudaSuccess); + REQUIRE(cudaFree(device_result) == cudaSuccess); + + REQUIRE(!result.torch_carries_matline); + REQUIRE(!result.frame_carries_matline); + REQUIRE(result.cerenkov_carries_matline); + REQUIRE(result.scintillation_carries_matline); + REQUIRE(result.modified_cerenkov_carries_matline); + + REQUIRE(result.entering_distance == 10.f); + REQUIRE(result.exiting_distance == nextafterf(10.f, INFINITY)); + REQUIRE(result.entering_distance < result.exiting_distance); + REQUIRE(result.non_simulation_exit_distance == 10.f); + + std::puts("CSGOptiXHelpersTest: PASS"); + return 0; +} + +#undef REQUIRE diff --git a/config/sibling_pair.json b/config/sibling_pair.json new file mode 100644 index 000000000..97ca98e8e --- /dev/null +++ b/config/sibling_pair.json @@ -0,0 +1,34 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 2000, + + "pos": [-10.0, 0.0, -50.1], + "time": 0.0, + + "mom": [0.3, 0.0, 1.0], + "weight": 0.0, + + "pol": [0.0, 1.0, 0.0], + "wavelength": 420.0, + + "zenith": [-1.0, 1.0], + "azimuth": [0.25, 0.25], + + "radius": 0.05, + "distance": 0.0, + "mode": 255, + "type": "disc" + }, + + "event": { + "mode": "HitPhoton", + "maxslot": 10000, + "max_bounce": 8, + "propagate_epsilon": 0.05, + "propagate_epsilon0": 0.05, + "propagate_epsilon0_mask": "TO,CK,SI,SC,RE" + } +} diff --git a/qudarap/tests/CMakeLists.txt b/qudarap/tests/CMakeLists.txt index 5deea0e09..2a5e173e2 100644 --- a/qudarap/tests/CMakeLists.txt +++ b/qudarap/tests/CMakeLists.txt @@ -18,6 +18,7 @@ set(TEST_SOURCES QTexMakerTest.cc QTexLookupTest.cc qbnd_boundary_lookup_test.cc + QSim_MaterialCarryTest.cc QMultiFilmTest.cc @@ -91,6 +92,10 @@ foreach(SRC ${TEST_SOURCES}) add_executable(${TGT} ${SRC}) target_link_libraries(${TGT} QUDARap CUDA::cudart) + if(TGT STREQUAL QSim_MaterialCarryTest) + target_compile_definitions(${TGT} PRIVATE MOCK_CUDA MOCK_CURAND) + endif() + if(TGT STREQUAL QEvt_Lifecycle_Test) foreach(INPUT_PHOTON_NAME ${QUDARAP_INPUT_PHOTON_NAMES}) set(TEST_NAME ${name}.${TGT}.${INPUT_PHOTON_NAME}) diff --git a/qudarap/tests/QSim_MaterialCarryTest.cc b/qudarap/tests/QSim_MaterialCarryTest.cc new file mode 100644 index 000000000..34b237b34 --- /dev/null +++ b/qudarap/tests/QSim_MaterialCarryTest.cc @@ -0,0 +1,121 @@ +#include +#include + +#include "scerenkov.h" +#include "scuda.h" +#include "smath.h" +#include "srec.h" +#include "srngcpu.h" +using RNG = srngcpu; + +#include "stexture.h" + +#include "OpticksPhoton.h" +#include "qsim.h" + +namespace +{ +void require(bool condition, const char* expression, int line) +{ + if (condition) + return; + std::fprintf(stderr, "%s:%d requirement failed: %s\n", __FILE__, line, expression); + std::exit(EXIT_FAILURE); +} + +#define REQUIRE(expression) require((expression), #expression, __LINE__) + +constexpr unsigned BOUNDARY_INDEX = 2u; +constexpr unsigned OUTER_MATERIAL_LINE = + BOUNDARY_INDEX * _BOUNDARY_NUM_MATSUR + OMAT; +constexpr unsigned INNER_MATERIAL_LINE = + BOUNDARY_INDEX * _BOUNDARY_NUM_MATSUR + IMAT; + +void initialize(sctx& ctx, quad2& prd) +{ + prd.zero(); + prd.q0.f.z = 1.f; + prd.q1.u.w = BOUNDARY_INDEX; + + ctx = {}; + ctx.prd = &prd; + ctx.pidx = 0u; + ctx.p.zero(); + ctx.p.mom = make_float3(0.f, 0.f, -1.f); + ctx.p.pol = make_float3(1.f, 0.f, 0.f); + ctx.p.wavelength = 420.f; + + ctx.s.material1 = make_float4(1.f, 0.f, 0.f, 0.f); + ctx.s.material2 = make_float4(1.f, 0.f, 0.f, 0.f); + ctx.s.m1group2 = make_float4(100.f, 200.f, 0.f, 0.f); + ctx.s.index = make_uint4(10u, 20u, 0u, 0u); +} + +void test_transmission_updates_carried_material(const qsim& sim, RNG& rng) +{ + sctx ctx = {}; + quad2 prd = {}; + initialize(ctx, prd); + ctx.current_matline = 99u; + + unsigned flag = 0u; + const FlowAction action = sim.propagate_at_boundary(flag, rng, ctx, 1.f); + + REQUIRE(action == FlowAction::Continue); + REQUIRE(flag == BOUNDARY_TRANSMIT); + REQUIRE(ctx.current_matline == INNER_MATERIAL_LINE); + REQUIRE(ctx.current_material_index == 20u); + REQUIRE(ctx.current_group_velocity == 200.f); +} + +void test_reflection_preserves_carried_material(const qsim& sim, RNG& rng) +{ + sctx ctx = {}; + quad2 prd = {}; + initialize(ctx, prd); + ctx.current_matline = 99u; + + unsigned flag = 0u; + const FlowAction action = sim.propagate_at_boundary(flag, rng, ctx, 0.f); + + REQUIRE(action == FlowAction::Continue); + REQUIRE(flag == BOUNDARY_REFLECT); + REQUIRE(ctx.current_matline == 99u); + REQUIRE(ctx.current_material_index == 10u); + REQUIRE(ctx.current_group_velocity == 100.f); +} + +void test_outward_transmission_selects_outer_material(const qsim& sim, RNG& rng) +{ + sctx ctx = {}; + quad2 prd = {}; + initialize(ctx, prd); + ctx.p.mom = make_float3(0.f, 0.f, 1.f); + + unsigned flag = 0u; + const FlowAction action = sim.propagate_at_boundary(flag, rng, ctx, 1.f); + + REQUIRE(action == FlowAction::Continue); + REQUIRE(flag == BOUNDARY_TRANSMIT); + REQUIRE(ctx.current_matline == OUTER_MATERIAL_LINE); +} +} // namespace + +int main() +{ + qbase base = {}; + base.pidx = ~0ull; + + qsim sim; + sim.base = &base; + + RNG rng; + rng.set_fake(0.5); + + test_transmission_updates_carried_material(sim, rng); + test_reflection_preserves_carried_material(sim, rng); + test_outward_transmission_selects_outer_material(sim, rng); + + std::puts("QSim_MaterialCarryTest: PASS"); + return 0; +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6ad1f41ef..725d4167f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -29,6 +29,7 @@ add_simg4ox_integration_test(dune_mock_wls_detector_box) add_simg4ox_integration_test(8x8SiPM_w_CSI_optial_grease) add_simg4ox_integration_test(opticks_two_spheres) add_simg4ox_integration_test(drich) +add_simg4ox_integration_test(sibling_pair) set(SIMPHONY_SIMG4OX_MULTIEVENT_TEST_WORKDIR "${CMAKE_CURRENT_BINARY_DIR}/simg4ox_multievent") file(MAKE_DIRECTORY "${SIMPHONY_SIMG4OX_MULTIEVENT_TEST_WORKDIR}") diff --git a/tests/geom/sibling_pair.gdml b/tests/geom/sibling_pair.gdml new file mode 100644 index 000000000..05c01350d --- /dev/null +++ b/tests/geom/sibling_pair.gdml @@ -0,0 +1,97 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_simg4ox.sh b/tests/test_simg4ox.sh index 648a51fed..cf33613c5 100755 --- a/tests/test_simg4ox.sh +++ b/tests/test_simg4ox.sh @@ -22,6 +22,7 @@ usage() { echo " 8x8SiPM_w_CSI_optial_grease" >&2 echo " opticks_two_spheres" >&2 echo " drich" >&2 + echo " sibling_pair" >&2 } if [[ $# -gt 1 ]]; then @@ -91,6 +92,12 @@ case "${TEST_CASE}" in # Invoke Geant4 SDs for the EFFICIENCY=1 SiPM skin surfaces. run_hit_validation drich drich "${REPO_DIR}/tests/run_validate.mac" ;; + sibling_pair) + # The beam can reach the sensor only by crossing the exact shared face + # from the left sibling into the right sibling while retaining the + # correct current-medium properties for the segment ending at that face. + run_hit_validation sibling_pair sibling_pair "${REPO_DIR}/tests/run_validate.mac" + ;; *) echo "Unknown simg4ox test case: ${TEST_CASE}" >&2 usage