qbnd/qsim: carry source-medium matline across CSG sibling-pair faces

CSGOptiX/CSGOptiX.cc

@@ -609,6 +609,59 @@ 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<double>();    // stree/standard bnd is <f8 (double)
+            const unsigned             LMAT_OMAT = 0u, LMAT_IMAT = 3u; // sstandard mat/sur line order {OMAT,OSUR,ISUR,IMAT}
+            const float                DN_THRESHOLD = 0.1f;
+            std::vector<unsigned char> 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() ;
 

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;
+    }
 
     int command = START ;
     int bounce = 0 ;
@@ -906,8 +914,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<quad2>(); // access prd addr from RG program
             prd->q0.f = isect ;  // .w:distance and .xyz:normal which starts as the local frame one
@@ -923,7 +942,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

CSGOptiX/Params.cc

@@ -126,9 +126,7 @@ std::string Params::detail() const
     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),
@@ -147,6 +145,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),

CSGOptiX/Params.h

@@ -1,3 +1,3 @@
 #pragma once
 
 /**
@@ -64,6 +64,7 @@
     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 ;

qudarap/qbnd.h

@@ -1,3 +1,3 @@
 #pragma once
 /**
 qbnd.h
@@ -38,7 +38,7 @@
 #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
 
 };
@@ -183,14 +183,30 @@
 
 **/
 
-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
@@ -200,7 +216,7 @@
 
     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

qudarap/qsim.h

@@ -1263,6 +1263,19 @@ inline QSIM_METHOD int qsim::propagate_at_boundary(unsigned& flag, RNG& rng, sct
     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 )
     {
@@ -2153,7 +2166,7 @@ inline QSIM_METHOD int qsim::propagate(const int bounce, RNG& rng, sctx& ctx )
     // 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)
     {

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 ;