fix: align post-boundary optical timing with Geant4 11.4+ (#260)

Geant4 11.4 substantially refactored G4OpBoundaryProcess, including the
dielectric-dielectric boundary handling path. In our simg4ox and
GPURaytrace workflows, Geant4 11.4+ is the first release family where
timing-sensitive optical comparisons diverge from the Geant4 11.3
baselines. That makes post-boundary timing correctness explicit in a way
our previous implementation did not handle.

This branch fixes two gaps in our handling of post-boundary optical
timing:

1. The Opticks GPU path advanced the next segment using the pre-boundary
   medium group velocity even after transmission into material2.

2. The lightweight Geant4-side tracking/recording path used by simg4ox
   did not explicitly enable UseGivenVelocity(true) for optical photons,
   so the track was not forced to follow the boundary-selected velocity
   during post-boundary recording.

Fix the GPU side by carrying the active-medium timing state through
qsim:
- store both material1 and material2 group velocities in sstate
- track current_material_index/current_group_velocity in sctx
- use current_group_velocity in propagate_to_boundary
- switch the carried material/velocity in propagate_at_boundary after
  reflection/transmission
- refresh the carried velocity after each boundary lookup when the
  boundary orientation still matches the active medium
- initialize the new fields in the QSim test/debug kernels

Fix the Geant4-side CPU recording path by enabling
UseGivenVelocity(true) for optical photons in g4app tracking action, so
the recorded post-boundary timing follows the velocity selected by
Geant4.

Because hit merging in this codebase is time-bucket based, correcting
post-boundary timing changes timing-sensitive observables and can change
merged hit counts. The follow-up test updates on this branch therefore
make expected residual diffs and hit-count baselines depend on the
Geant4 release family instead of treating the new values as regressions.

---------

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

Author: plexoos

qudarap/QSim.cu

@@ -324,6 +324,8 @@ __global__ void _QSim_propagate_to_boundary( qsim* sim, sphoton* photon, unsigne
     ctx.prd = &dbg->prd ;  // no need for local copy when readonly
     ctx.s = dbg->s ;
     ctx.p = dbg->p ;      // need local copy of photon otherwise will have write interference between threads
+    ctx.current_group_velocity = ctx.s.material1_group_velocity();
+    ctx.current_material_index = ctx.s.index.x;
 
     unsigned flag = 0u ;
     //sim->propagate_to_boundary( flag, p, prd, s, rng, idx, tagr );
@@ -357,6 +359,8 @@ __global__ void _QSim_propagate_at_boundary_generate( qsim* sim, sphoton* photon
     ctx.prd = &dbg->prd ;  // no need for local copy when readonly
     ctx.s = dbg->s ;
     ctx.p = dbg->p ;    // need local copy of photon otherwise will have write interference between threads
+    ctx.current_group_velocity = ctx.s.material1_group_velocity();
+    ctx.current_material_index = ctx.s.index.x;
 
     quad4& q = (quad4&)ctx.p ;
     q.q0.f = q.q1.f ;   // non-standard record initial mom and pol into q0, q3
@@ -402,6 +406,8 @@ __global__ void _QSim_propagate_at_boundary_mutate( qsim* sim, sphoton* photon,
     ctx.p = photon[idx] ;
     ctx.s = dbg->s ;
     ctx.prd = &dbg->prd ;
+    ctx.current_group_velocity = ctx.s.material1_group_velocity();
+    ctx.current_material_index = ctx.s.index.x;
 
     quad4&  q  = (quad4&)ctx.p ;
     q.q0.f = q.q1.f ;   // non-standard record initial mom and pol into q0, q3
@@ -441,6 +447,8 @@ __global__ void _QSim_propagate_at_multifilm_mutate( qsim* sim, sphoton* photon,
     ctx.p = photon[idx] ;
     ctx.s = dbg->s ;
     ctx.prd = &dbg->prd ;
+    ctx.current_group_velocity = ctx.s.material1_group_velocity();
+    ctx.current_material_index = ctx.s.index.x;
 
     quad4&  q  = (quad4&)ctx.p ;
 
@@ -820,9 +828,3 @@ extern void QSim_multifilm_lookup_all(dim3 numBlocks, dim3 threadsPerBlock, qsim
     printf("//QSim_multifilm_lookup width %d  height %d \n", width, height );
     _QSim_multifilm_lookup_all<<<numBlocks,threadsPerBlock>>>( sim, sample,result , width, height );
 }
-
-
-
-
-
-

qudarap/QState.hh

@@ -51,6 +51,7 @@ inline sstate QState::Make()
     float m1_scattering_length = 1000.f ; 
     float m1_reemission_prob = 0.f ; 
     float m1_group_velocity = 300.f ; 
+    float m2_group_velocity = 300.f;
 
     float m2_refractive_index = qenvfloat("M2_REFRACTIVE_INDEX", "1.5" ) ; 
     float m2_absorption_length = 1000.f ; 
@@ -65,7 +66,7 @@ inline sstate QState::Make()
     sstate s ; 
     s.material1 = make_float4( m1_refractive_index, m1_absorption_length, m1_scattering_length, m1_reemission_prob ); 
     s.material2 = make_float4( m2_refractive_index, m2_absorption_length, m2_scattering_length, m2_reemission_prob );  
-    s.m1group2  = make_float4( m1_group_velocity, 0.f, 0.f, 0.f ); 
+    s.m1group2  = make_float4(m1_group_velocity, m2_group_velocity, 0.f, 0.f);
     s.surface   = make_float4( su_detect, su_absorb, su_reflect_specular, su_reflect_diffuse ); 
     s.optical   = make_uint4( 0u, 0u, 0u, 0u );  // x/y/z/w index/type/finish/value  
     s.index     = make_uint4( 0u, 0u, 0u, 0u );  // indices of m1/m2/surf/sensor
@@ -144,5 +145,3 @@ inline std::string QState::Desc( const sstate& s )
     std::string repr = ss.str(); 
     return repr ; 
 }
-
-

qudarap/qbnd.h

@@ -166,7 +166,9 @@ Notice that s.optical.x and s.index.z are the same thing.
 So half of s.index is extraneous and the m1 index and m2 index
 is not much used.
 
-Also only one elemnt of m1group2 is actually used
+Also only one element of m1group2 was formerly used. The y slot now carries
+material2 group velocity so transmitted segments can carry the correct active
+medium timing across boundary crossings.
 
 
 s.optical.x
@@ -191,9 +193,10 @@ inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wav
 
 
     s.material1 = boundary_lookup( wavelength, m1_line, 0);   // refractive_index, absorption_length, scattering_length, reemission_prob
-    s.m1group2  = boundary_lookup( wavelength, m1_line, 1);   // group_velocity ,  (unused          , unused           , unused)
+    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)
 

qudarap/qsim.h

@@ -723,7 +723,7 @@ inline QSIM_METHOD int qsim::propagate_to_boundary(unsigned& flag, RNG& rng, sct
     const float& absorption_length = s.material1.y ;
     const float& scattering_length = s.material1.z ;
     const float& reemission_prob = s.material1.w ;
-    const float& group_velocity = s.m1group2.x ;
+    const float& group_velocity = ctx.current_group_velocity;
     const float& distance_to_boundary = ctx.prd->q0.f.w ;
 
 
@@ -1184,6 +1184,8 @@ inline QSIM_METHOD int qsim::propagate_at_boundary(unsigned& flag, RNG& rng, sct
 
     flag = reflect ? BOUNDARY_REFLECT : BOUNDARY_TRANSMIT ;
 
+    ctx.current_material_index = reflect ? s.index.x : s.index.y;
+    ctx.current_group_velocity = reflect ? s.material1_group_velocity() : s.material2_group_velocity();
 
 #if !defined(PRODUCTION) && defined(DEBUG_TAG)
     if( flag ==  BOUNDARY_REFLECT )
@@ -2249,6 +2251,18 @@ inline QSIM_METHOD int qsim::propagate(const int bounce, RNG& rng, sctx& ctx )
 
     bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx );
 
+    if (ctx.current_material_index == 0u)
+    {
+        ctx.current_material_index = ctx.s.index.x;
+        ctx.current_group_velocity = ctx.s.material1_group_velocity();
+    }
+    else if (ctx.s.index.x == ctx.current_material_index)
+    {
+        // When the boundary orientation agrees with the carried medium, refresh
+        // the active velocity from the wavelength-dependent material lookup.
+        ctx.current_group_velocity = ctx.s.material1_group_velocity();
+    }
+
     unsigned flag = 0 ;
 
     int command = propagate_to_boundary( flag, rng, ctx );
@@ -2540,4 +2554,3 @@ inline QSIM_METHOD void qsim::generate_photon(sphoton& p, RNG& rng, const quad6&
     p.set_index(photon_id) ;
 }
 #endif
-

src/g4app.h

@@ -376,6 +376,14 @@ struct TrackingAction : G4UserTrackingAction
 
     void PreUserTrackingAction(const G4Track *track) override
     {
+        if (track->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
+        {
+            // Geant4 boundary updates optical velocity via ProposeVelocity, but the
+            // track must honor the given velocity for post-boundary timing to match.
+            G4Track *mutable_track = const_cast<G4Track *>(track);
+            mutable_track->UseGivenVelocity(true);
+        }
+
         spho *label = STrackInfo::GetRef(track);
 
         if (label == nullptr)

sysrap/sctx.h

@@ -77,6 +77,8 @@ struct sctx
 
     sphoton p ;
     sstate  s ;
+    float current_group_velocity = 0.f;
+    unsigned current_material_index = 0u;
 
 #ifndef PRODUCTION
     srec rec ;
@@ -184,4 +186,3 @@ SCTX_METHOD void sctx::end()
 }
 
 #endif
-

sysrap/sstate.h

@@ -25,7 +25,7 @@ BUT seems no point doing that, can just directly use them from PRD.
 struct sstate
 {
     float4 material1 ;    // refractive_index/absorption_length/scattering_length/reemission_prob
-    float4 m1group2 ;     // group_velocity/spare1/spare2/spare3
+    float4 m1group2 ;     // material1_group_velocity/material2_group_velocity/spare2/spare3
     float4 material2 ;   
     float4 surface ;      // detect/absorb/reflect_specular/reflect_diffuse
 
@@ -42,6 +42,9 @@ struct sstate
     void save(const char* dir) const ; 
 #endif
 
+    SSTATE_METHOD float material1_group_velocity() const { return m1group2.x; }
+    SSTATE_METHOD float material2_group_velocity() const { return m1group2.y; }
+    SSTATE_METHOD void  set_material2_group_velocity(float gv) { m1group2.y = gv; }
 
 };
 
@@ -71,7 +74,7 @@ inline std::ostream& operator<<(std::ostream& os, const sstate& s )
        << " (refractive_index/absorption_length/scattering_length/reemission_prob) " 
        << std::endl 
        << " m1group2 " << s.m1group2
-       << " (group_velocity/spare1/spare2/spare3) "
+       << " (material1_group_velocity/material2_group_velocity/spare2/spare3) "
        << std::endl 
        << " material2 " << s.material2 
        << " (refractive_index/absorption_length/scattering_length/reemission_prob) " 
@@ -134,6 +137,3 @@ But if decide to pursure streamlined qstate with mixed up fields
 the named field access would be helpful to isolate user code from changes to the struct. 
 
 **/
-
-
-

tests/compare_ab.py

@@ -1,14 +1,37 @@
 import numpy as np
 
+from optiphy.geant4_version import detect_geant4_version, geant4_series
+
+
+EXPECTED_DIFF = {
+    "11.3": [14, 22, 32, 34, 40, 81, 85],
+    "11.4+": [0, 30, 32, 34, 42, 69, 78, 85, 86],
+}
+
+
+def expected_diff_for_version(version):
+    return EXPECTED_DIFF[geant4_series(version)]
+
+
+geant4_version = detect_geant4_version()
+expected_diff = expected_diff_for_version(geant4_version)
+
 a = np.load("/tmp/fakeuser/opticks/GEOM/fakegeom/simg4ox/ALL0_no_opticks_event_name/A000/record.npy")
 b = np.load("/tmp/fakeuser/opticks/GEOM/fakegeom/simg4ox/ALL0_no_opticks_event_name/B000/f000/record.npy")
 
 print(a.shape)
 print(b.shape)
+print(f"GEANT4_VERSION={geant4_version}")
+print(f"EXPECTED_DIFF={expected_diff}")
 
 assert a.shape == b.shape
 
-diff = [i for i, (a, b) in enumerate(zip(a[:, 1:], b[:, 0:-1])) if not np.allclose(a, b, rtol=0, atol=1e-5)]
+# Geant4 and Opticks record one-step-shifted sequences for this test geometry,
+# so compare the aligned slices directly, including time.
+a_cmp = a[:, 1:]
+b_cmp = b[:, 0:-1]
+
+diff = [i for i, (ac, bc) in enumerate(zip(a_cmp, b_cmp)) if not np.allclose(ac, bc, rtol=0, atol=1e-5)]
 print(diff)
 
-assert diff == [14, 22, 32, 34, 40, 81, 85]
+assert diff == expected_diff

tests/test_GPURaytrace.sh

@@ -2,11 +2,50 @@
 
 set -e
 
-
 SEED=42
 TOLERANCE=113
 PASS=true
 
+GEANT4_VERSION=$(geant4-version version)
+GEANT4_SERIES=$(geant4-version series)
+
+declare -Ar SUPPORTED_GEANT4_SERIES=(
+    ["11.3"]=1
+    ["11.4+"]=1
+)
+
+declare -Ar EXPECTED_G4_HITS=(
+    ["11.3:cerenkov"]=12672
+    ["11.3:cerenkov_scintillation"]=17473
+    ["11.4+:cerenkov"]=11842
+    ["11.4+:cerenkov_scintillation"]=21411
+)
+
+declare -Ar EXPECTED_OPTICKS_HITS=(
+    ["11.3:cerenkov"]=12664
+    ["11.3:cerenkov_scintillation"]=17443
+    ["11.4+:cerenkov"]=11827
+    ["11.4+:cerenkov_scintillation"]=21390
+)
+
+set_expected_hits() {
+    local test_name=$1
+    local key="${GEANT4_SERIES}:${test_name}"
+
+    if [[ -z ${SUPPORTED_GEANT4_SERIES[$GEANT4_SERIES]+x} ]]; then
+        echo "FAILED: unsupported Geant4 version $GEANT4_VERSION. Add hit-count references for this release."
+        exit 1
+    fi
+
+    if [[ -z ${EXPECTED_G4_HITS[$key]+x} || -z ${EXPECTED_OPTICKS_HITS[$key]+x} ]]; then
+        echo "FAILED: unknown test name: $test_name"
+        exit 1
+    fi
+
+    EXPECTED_G4=${EXPECTED_G4_HITS[$key]}
+    EXPECTED_OPTICKS=${EXPECTED_OPTICKS_HITS[$key]}
+}
+
 check_hits() {
     local label=$1
     local actual=$2
@@ -21,8 +60,12 @@ check_hits() {
     fi
 }
 
+echo "Using Geant4 version $GEANT4_VERSION"
+
 # ---- Test 1: Cerenkov only (opticks_raindrop.gdml) ----
 echo "=== Test 1: Cerenkov only (opticks_raindrop.gdml) ==="
+set_expected_hits cerenkov
+
 echo "Running GPURaytrace with seed $SEED ..."
 OUTPUT=$(USER=fakeuser GEOM=fakegeom GPURaytrace \
     -g "$OPTICKS_HOME/tests/geom/opticks_raindrop.gdml" \
@@ -32,15 +75,17 @@ OUTPUT=$(USER=fakeuser GEOM=fakegeom GPURaytrace \
 G4_HITS=$(echo "$OUTPUT" | grep "Geant4: NumHits:" | awk '{print $NF}')
 OPTICKS_HITS=$(echo "$OUTPUT" | grep "Opticks: NumHits:" | awk '{print $NF}')
 
-echo "Geant4:  NumHits: $G4_HITS  (expected 12672 +/- $TOLERANCE)"
-echo "Opticks: NumHits: $OPTICKS_HITS  (expected 12664 +/- $TOLERANCE)"
+echo "Geant4:  NumHits: $G4_HITS  (expected $EXPECTED_G4 +/- $TOLERANCE)"
+echo "Opticks: NumHits: $OPTICKS_HITS  (expected $EXPECTED_OPTICKS +/- $TOLERANCE)"
 
-check_hits "Geant4 NumHits"  "$G4_HITS"      12672
-check_hits "Opticks NumHits" "$OPTICKS_HITS"  12664
+check_hits "Geant4 NumHits"  "$G4_HITS"      "$EXPECTED_G4"
+check_hits "Opticks NumHits" "$OPTICKS_HITS" "$EXPECTED_OPTICKS"
 
 # ---- Test 2: Cerenkov + Scintillation (opticks_raindrop_with_scintillation.gdml) ----
 echo ""
 echo "=== Test 2: Cerenkov + Scintillation (opticks_raindrop_with_scintillation.gdml) ==="
+set_expected_hits cerenkov_scintillation
+
 echo "Running GPURaytrace with seed $SEED ..."
 OUTPUT=$(USER=fakeuser GEOM=fakegeom GPURaytrace \
     -g "$OPTICKS_HOME/tests/geom/opticks_raindrop_with_scintillation.gdml" \
@@ -50,11 +95,11 @@ OUTPUT=$(USER=fakeuser GEOM=fakegeom GPURaytrace \
 G4_HITS=$(echo "$OUTPUT" | grep "Geant4: NumHits:" | awk '{print $NF}')
 OPTICKS_HITS=$(echo "$OUTPUT" | grep "Opticks: NumHits:" | awk '{print $NF}')
 
-echo "Geant4:  NumHits: $G4_HITS  (expected 17473 +/- $TOLERANCE)"
-echo "Opticks: NumHits: $OPTICKS_HITS  (expected 17443 +/- $TOLERANCE)"
+echo "Geant4:  NumHits: $G4_HITS  (expected $EXPECTED_G4 +/- $TOLERANCE)"
+echo "Opticks: NumHits: $OPTICKS_HITS  (expected $EXPECTED_OPTICKS +/- $TOLERANCE)"
 
-check_hits "Geant4 NumHits"  "$G4_HITS"      17473
-check_hits "Opticks NumHits" "$OPTICKS_HITS"  17443
+check_hits "Geant4 NumHits"  "$G4_HITS"      "$EXPECTED_G4"
+check_hits "Opticks NumHits" "$OPTICKS_HITS" "$EXPECTED_OPTICKS"
 
 # ---- Summary ----
 echo ""