qsim.h

/**
qsim.h : GPU side struct prepared CPU side by QSim.hh
========================================================

qsim.h replaces the OptiX 6 context in a CUDA-centric way.
Canonical use is from CSGOptiX/CSGOptiX7.cu:simulate

* qsim.h instance is uploaded once only at CSGOptiX instanciation
  as this encompasses the physics not the event-by-event info.

* qsim encompasses global info relevant to all photons, meaning that any changes
  made to the qsim instance from single photon threads must be into thread-owned "idx"
  slots into arrays to avoid interference

* temporary working state local to each photon is held in *sctx*
  and passed around using reference arguments

TODO:

1. get more of the below to work on CPU with mocked curand (and in future mocked tex2D and cudaTextureObject_t )
   NB must move decl and implementation to do this

**/

#if defined(__CUDACC__) || defined(__CUDABE__)
   #define QSIM_METHOD __device__
#else
   #define QSIM_METHOD
#endif

#include "OpticksGenstep.h"
#include "OpticksPhoton.h"

#include "sflow.h"
#include "sqat4.h"
#include "sc4u.h"
#include "sxyz.h"
#include "sphoton.h"

#include "storch.h"
#include "scarrier.h"
#include "sevent.h"
#include "sstate.h"
#include "smatsur.h"


#ifndef PRODUCTION
#include "srec.h"
#include "sseq.h"
#include "stag.h"
#ifdef DEBUG_LOGF
#define KLUDGE_FASTMATH_LOGF(u) (u < 0.998f ? __logf(u) : __logf(u) - 0.46735790f*1e-7f )
#endif
#endif

#include "sctx.h"

#include "qrng.h"
#include "qbase.h"
#include "qprop.h"
#include "qmultifilm.h"
#include "qbnd.h"
#include "qscint.h"
#include "qcerenkov.h"
#include "qpmt.h"
#include "tcomplex.h"


struct qcerenkov ;

struct qsim
{
    qbase*              base ;
    sevent*             evt ;
    qrng<RNG>*          rng ;
    qbnd*               bnd ;
    qmultifilm*         multifilm;
    qcerenkov*          cerenkov ;
    qscint*             scint ;
    qpmt<float>*        pmt ;

#if defined(__CUDACC__) || defined(__CUDABE__)
#else
    qsim(); // instanciated on CPU (see QSim::init_sim) and copied to device so no ctor in device code
#endif

    QSIM_METHOD void    generate_photon_dummy( sphoton& p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const ;
    QSIM_METHOD static float3 uniform_sphere(const float u0, const float u1);
    QSIM_METHOD static float RandGaussQ_shoot( RNG& rng, float mean, float stdDev );
    QSIM_METHOD static void SmearNormal_SigmaAlpha( RNG& rng, float3* smeared_normal, const float3* direction, const float3* normal, float sigma_alpha, const sctx& ctx );
    QSIM_METHOD static void SmearNormal_Polish(     RNG& rng, float3* smeared_normal, const float3* direction, const float3* normal, float polish     , const sctx& ctx );

#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
    QSIM_METHOD static float3 uniform_sphere(RNG& rng);
#endif

#if defined(__CUDACC__) || defined(__CUDABE__)
    QSIM_METHOD float4  multifilm_lookup(unsigned pmtType, float nm, float aoi);
#endif

#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND )  || defined(MOCK_CUDA)
    QSIM_METHOD static void lambertian_direction(float3* dir, const float3* normal, float orient, RNG& rng, sctx& ctx );
    QSIM_METHOD static void random_direction_marsaglia(float3* dir, RNG& rng, sctx& ctx );
    QSIM_METHOD void rayleigh_scatter(RNG& rng, sctx& ctx );
    QSIM_METHOD int     propagate_to_boundary( unsigned& flag, RNG& rng, sctx& ctx );
#endif

#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
    QSIM_METHOD int     propagate_at_boundary(        unsigned& flag, RNG& rng, sctx& ctx, float theTransmittance=-1.f ) const ;
    QSIM_METHOD int     propagate_at_boundary_with_T( unsigned& flag, RNG& rng, sctx& ctx, float theTransmittance ) const ;
#endif

#if defined(__CUDACC__) || defined(__CUDABE__)
    QSIM_METHOD int     propagate_at_surface_MultiFilm(unsigned& flag, RNG& rng, sctx& ctx );
#endif

#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
    QSIM_METHOD int     propagate_at_surface(           unsigned& flag, RNG& rng, sctx& ctx );
    QSIM_METHOD int     propagate_at_surface_Detect(    unsigned& flag, RNG& rng, sctx& ctx ) const ;
#if defined( WITH_CUSTOM4 )
    QSIM_METHOD int     propagate_at_surface_CustomART( unsigned& flag, RNG& rng, sctx& ctx ) const ;
#endif
#endif

#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
    QSIM_METHOD void    reflect_diffuse(                       RNG& rng, sctx& ctx );
    QSIM_METHOD void    reflect_specular(                      RNG& rng, sctx& ctx );

    QSIM_METHOD void    fake_propagate( sphoton& p, const quad2* mock_prd, RNG& rng, unsigned long long idx );
    QSIM_METHOD int     propagate(const int bounce, RNG& rng, sctx& ctx );

    QSIM_METHOD void    hemisphere_polarized( unsigned polz, bool inwards, RNG& rng, sctx& ctx );
    QSIM_METHOD void    generate_photon_simtrace(         quad4&   p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const ;
    QSIM_METHOD void    generate_photon_simtrace_frame(   quad4&   p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const ;
    QSIM_METHOD void    generate_photon(                  sphoton& p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const ;
#endif
};

// CTOR
#if defined(__CUDACC__) || defined(__CUDABE__)
#else
inline qsim::qsim()    // instanciated on CPU (see QSim::init_sim) and copied to device so no ctor in device code
        :
        base(nullptr),
        evt(nullptr),
        rng(nullptr),
        bnd(nullptr),
        multifilm(nullptr),
        cerenkov(nullptr),
        scint(nullptr),
        pmt(nullptr)
    {
    }
#endif

inline QSIM_METHOD void qsim::generate_photon_dummy(sphoton& p_, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const
{
    quad4& p = (quad4&)p_ ;
#ifndef PRODUCTION
    printf("//qsim::generate_photon_dummy  photon_id %3lld genstep_id %3d  gs.q0.i ( gencode:%3d %3d %3d %3d ) \n",
       photon_id,
       genstep_id,
       gs.q0.i.x,
       gs.q0.i.y,
       gs.q0.i.z,
       gs.q0.i.w
      );
#endif
    p.q0.i.x = 1 ; p.q0.i.y = 2 ; p.q0.i.z = 3 ; p.q0.i.w = 4 ;
    p.q1.i.x = 1 ; p.q1.i.y = 2 ; p.q1.i.z = 3 ; p.q1.i.w = 4 ;
    p.q2.i.x = 1 ; p.q2.i.y = 2 ; p.q2.i.z = 3 ; p.q2.i.w = 4 ;
    p.q3.i.x = 1 ; p.q3.i.y = 2 ; p.q3.i.z = 3 ; p.q3.i.w = 4 ;

    p.set_flag(TORCH);
}

inline QSIM_METHOD float3 qsim::uniform_sphere(const float u0, const float u1)
{
    float phi = u0*2.f*M_PIf;
    float cosTheta = 2.f*u1 - 1.f ; // -1.f -> 1.f
    float sinTheta = sqrtf(1.f-cosTheta*cosTheta);
    return make_float3(cosf(phi)*sinTheta, sinf(phi)*sinTheta, cosTheta);
}


#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
/**
qsim::uniform_sphere
---------------------

**/
inline QSIM_METHOD float3 qsim::uniform_sphere(RNG& rng)
{
    float phi = curand_uniform(&rng)*2.f*M_PIf;
    float cosTheta = 2.f*curand_uniform(&rng) - 1.f ; // -1.f -> 1.f
    float sinTheta = sqrtf(1.f-cosTheta*cosTheta);
    return make_float3(cosf(phi)*sinTheta, sinf(phi)*sinTheta, cosTheta);
}

/**
qsim::RandGaussQ_shoot
------------------------

See::

    sysrap/tests/erfcinvf_Test.sh
    sysrap/tests/S4MTRandGaussQTest.sh

    g4-cls RandGaussQ
    g4-cls G4MTRandGaussQ

**/
inline QSIM_METHOD float qsim::RandGaussQ_shoot( RNG& rng, float mean, float stdDev )
{
    float u2 = 2.f*curand_uniform(&rng) ;
    float v = -M_SQRT2f*erfcinvf(u2)*stdDev + mean ;
    //printf("//qsim.RandGaussQ_shoot mean %10.5f stdDev %10.5f u2 %10.5f v %10.5f \n", mean, stdDev, u2, v  ) ;
    return v ;
}


/**
qsim::SmearNormal_SigmaAlpha
------------------------------

CAUTION : THIS CURRENTLY NOT USED BY ANYTHING OTHER THAN TESTS

The reason is that a Geant4 simulation with some sigma_alpha/polish
surfaces will not actually use G4OpBoundaryProcess::GetFacetNormal unless
a bunch of conditions are satisfied such as having non-zero values of some
of prob_sl/prob_ss/prob_bs that induces G4OpBoundaryProcess::ChooseReflection
to return something other that LambertianReflection.

**THIS MAKES ME SUSPECT ACTUAL USE OF NORMAL SMEARING IS VERY RARE**

+----------------------+-------------------------------+------------------------------+----------------------+
| G4OpBoundaryProcess  | G4MaterialPropertiesIndex.hh  | G4MaterialPropertiesTable.cc | ChooseReflection     |
+======================+===============================+==============================+======================+
| prob_sl              | kSPECULARLOBECONSTANT         | SPECULARLOBECONSTANT         | LobeReflection       |
+----------------------+-------------------------------+------------------------------+----------------------+
| prob_ss              | kSPECULARSPIKECONSTANT        | SPECULARSPIKECONSTANT        | SpikeReflection      |
+----------------------+-------------------------------+------------------------------+----------------------+
| prob_bs              | kBACKSCATTERCONSTANT          | BACKSCATTERCONSTANT          | BackScattering       |
+----------------------+-------------------------------+------------------------------+----------------------+
| all zero =>          |                               |                              | LambertianReflection |
+----------------------+-------------------------------+------------------------------+----------------------+

TODO: full simulation run with breakpoint "BP=C4OpBoundaryProcess::GetFacetNormal"

* C4 (not G4) as Custom4 is in use for the boundary process

**/

inline QSIM_METHOD void qsim::SmearNormal_SigmaAlpha(
    RNG& rng,
    float3* smeared_normal,
    const float3* direction,
    const float3* normal,
    float sigma_alpha,
    const sctx& ctx
   )
{
#if !defined(PRODUCTION) && defined(MOCK_CUDA_DEBUG)
    bool dump = ctx.pidx == -1 ;
#endif

    if(sigma_alpha == 0.f)
    {
        *smeared_normal = *normal ;
        return ;
    }
    float f_max = fminf(1.f,4.f*sigma_alpha);

#if !defined(PRODUCTION) && defined(MOCK_CUDA_DEBUG)
    if(dump) printf("//qsim::SmearNormal_SigmaAlpha.MOCK_CUDA_DEBUG sigma_alpha %10.5f f_max %10.5f  \n", sigma_alpha, f_max );
#endif

    float alpha, sin_alpha, phi, u0, u1, u2 ;
    bool reject_alpha ;
    bool reject_dir ;

    do {
        do {
            //alpha = RandGaussQ_shoot(rng, 0.f, sigma_alpha );  // mean:0.f stdDev:sigma_alpha
            u0 = curand_uniform(&rng) ;
            alpha = -M_SQRT2f*erfcinvf(2.f*u0)*sigma_alpha ;

            sin_alpha = sinf(alpha);
            u1 = curand_uniform(&rng) ;
            reject_alpha = alpha >= M_PIf/2.f || (u1*f_max > sin_alpha) ;

#if !defined(PRODUCTION) && defined(MOCK_CUDA_DEBUG)
            if(dump) printf("//qsim::SmearNormal_SigmaAlpha.MOCK_CUDA_DEBUG u0 %10.5f alpha %10.5f sin_alpha %10.5f u1 %10.5f u1*f_max %10.5f  (u1*f_max > sin_alpha) %d reject_alpha %d  \n",
               u0, alpha, sin_alpha, u1, u1*f_max, (u1*f_max > sin_alpha), reject_alpha );
            // theres lots of alpha rejected : eg all -ve sin_alpha
#endif

        } while( reject_alpha ) ;

        u2 = curand_uniform(&rng) ;
        phi = u2*M_PIf*2.f ;

        smeared_normal->x = sin_alpha * cosf(phi) ;
        smeared_normal->y = sin_alpha * sinf(phi) ;
        smeared_normal->z = cosf(alpha) ;

        smath::rotateUz(*smeared_normal, *normal);
        reject_dir = dot(*smeared_normal, *direction ) >= 0.f ;
        // reject smears that move the normal into same hemi as direction

#if !defined(PRODUCTION) && defined(MOCK_CUDA_DEBUG)
        if(dump) printf("//qsim::SmearNormal_SigmaAlpha.MOCK_CUDA_DEBUG u2 %10.5f phi %10.5f smeared_normal ( %10.5f, %10.5f, %10.5f)  reject_dir %d  \n",
               u2, phi, smeared_normal->x, smeared_normal->y, smeared_normal->z, reject_dir );
#endif


    } while( reject_dir ) ;
}

/**
qsim::SmearNormal_Polish
------------------------------

CAUTION : THIS CURRENTLY NOT USED BY ANYTHING OTHER THAN TESTS : SEE DETAILS ABOVE

**/

inline QSIM_METHOD void qsim::SmearNormal_Polish(
    RNG& rng,
    float3* smeared_normal,
    const float3* direction,
    const float3* normal,
    float polish,
    const sctx& ctx
    )
{
#if !defined(PRODUCTION) && defined(MOCK_CUDA_DEBUG)
    bool dump = ctx.pidx == -1 ;
#endif

    if(polish == 1.f)
    {
        *smeared_normal = *normal ;
        return ;
    }

    float u0, u1, u2 ;
    float3 smear ;
    bool reject_mag ;
    bool reject_dir ;

    do {
        do {
            u0 = curand_uniform(&rng);
            u1 = curand_uniform(&rng) ;
            u2 = curand_uniform(&rng) ;
            smear.x = 2.f*u0 - 1.f ;
            smear.y = 2.f*u1 - 1.f ;
            smear.z = 2.f*u2 - 1.f ;
            reject_mag = length(smear) > 1.f  ;   // HMM: could this use just dot(smear, smear) ?
       }
       while( reject_mag );

       *smeared_normal = *normal + (1.f-polish)*smear;
       reject_dir = dot(*smeared_normal, *direction) >= 0.f ;
    }
    while( reject_dir );
    *smeared_normal = normalize(*smeared_normal);
}





#endif

#if defined(__CUDACC__) || defined(__CUDABE__)
/*
qsim::multifilm_lookup
-----------------------

*/

inline QSIM_METHOD float4 qsim::multifilm_lookup(unsigned pmtType, float nm, float aoi)
{
    float4 value = multifilm->lookup(pmtType, nm, aoi);
    return value;
}

#endif

#if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)

/**
qsim::lambertian_direction following G4LambertianRand
--------------------------------------------------------

g4-cls G4RandomTools::

     59 inline G4ThreeVector G4LambertianRand(const G4ThreeVector& normal)
     60 {
     61   G4ThreeVector vect;
     62   G4double ndotv;
     63   G4int count=0;
     64   const G4int max_trials = 1024;
     65
     66   do
     67   {
     68     ++count;
     69     vect = G4RandomDirection();
     70     ndotv = normal * vect;
     71
     72     if (ndotv < 0.0)
     73     {
     74       vect = -vect;
     75       ndotv = -ndotv;
     76     }
     77
     78   } while (!(G4UniformRand() < ndotv) && (count < max_trials));
     79
     80   return vect;
     81 }


NB: potentially bad for performance for dir pointer to be into global mem
as opposed to local stack float3 : as this keeps changing the dir before
arriving at the final one

**/
inline  QSIM_METHOD void qsim::lambertian_direction(float3* dir, const float3* normal, float orient, RNG& rng, sctx& ctx )
{
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    unsigned long long PIDX = 0xffffffffff ;
    if(ctx.pidx == PIDX )
    {
        printf("//qsim.lambertian_direction.head pidx %7lld : normal = np.array([%10.5f,%10.5f,%10.5f]) ; orient = %10.5f  \n",
            ctx.pidx, normal->x, normal->y, normal->z, orient  );
    }
#endif

    float ndotv ;
    int count = 0 ;
    float u ;
    do
    {
        count++ ;
        random_direction_marsaglia(dir, rng, ctx); // sets dir to random point on unit sphere
        ndotv = dot( *dir, *normal )*orient ;
        if( ndotv < 0.f )
        {
            *dir = -1.f*(*dir) ;
            ndotv = -1.f*ndotv ;
        }
        // when random dir is in opposite hemisphere to oriented normal
        // flip the dir into same hemi and ndotv

        u = curand_uniform(&rng) ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)

        if(ctx.pidx == PIDX)
        {
            printf("//qsim.lambertian_direction.loop pidx %7lld : dir = np.array([%10.5f,%10.5f,%10.5f]) ; count = %d ; ndotv = %10.5f ; u = %10.5f \n",
                ctx.pidx, dir->x, dir->y, dir->z, count, ndotv, u   );

        }
#endif
    }
    while (!(u < ndotv) && (count < 1024)) ;
    // distribution looks pretty similar without the while loop


#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == PIDX)
    {
        printf("//qsim.lambertian_direction.tail pidx %7lld : dir = np.array([%10.5f,%10.5f,%10.5f]) ; count = %d ; ndotv = %10.5f \n",
            ctx.pidx, dir->x, dir->y, dir->z, count, ndotv  );

    }
#endif


}

/**
qsim::random_direction_marsaglia following G4RandomDirection
-------------------------------------------------------------

* https://mathworld.wolfram.com/SpherePointPicking.html

::

                          v
                          |
              +---------.-|- -----------+
              |      .    |     .       |
              |   .       |          .  |
              |           |             |
              | .         |            .|
              |           |             |
          ----+-----------0-------------+---- u
              |.          |             |
              |           |            .|
              | .         |             |
              |           |          .  |
              |    .      |      .      |
              +--------.--|--.----------+
                          |

Marsaglia (1972) derived an elegant method that consists of picking u and v from independent
uniform distributions on (-1,1) and rejecting points for which uu+vv >=1.
So are picking points from within the (u,v) disc.

For those remaining random points on the 2D (u,v) disc the below (u,v) to (x,y,z)
mapping is used to give a 3D position on the unit sphere,

    x=2*u*sqrt(1-(uu+vv))
    y=2*v*sqrt(1-(uu+vv))
    z=1-2(uu+vv)

Checking normalization, it reduces to 1::

   xx + yy + zz =
         4uu (1-(uu+vv))
         4vv (1-(uu+vv)) +
        1 -4(uu+vv) + 4(uu+vv)(uu+vv)
                = 1

So that means the random 3D (x,y,z) points are on the unit sphere.

::

     g4-cls G4RandomDirection

     58 // G.Marsaglia (1972) method
     59 inline G4ThreeVector G4RandomDirection()
     60 {
     61   G4double u, v, b;
     62   do {
     63     u = 2.*G4UniformRand() - 1.;
     64     v = 2.*G4UniformRand() - 1.;
     65     b = u*u + v*v;
     66   } while (b > 1.);
     67   G4double a = 2.*std::sqrt(1. - b);
     68   return G4ThreeVector(a*u, a*v, 2.*b - 1.);
     69 }

**/


inline QSIM_METHOD void qsim::random_direction_marsaglia(float3* dir,  RNG& rng, sctx& ctx  )
{
    // NB: no use of ctx.tagr so this has not been random aligned
    float u0, u1 ;
    float u, v, b, a  ;
    do
    {
        u0 = curand_uniform(&rng);
        u1 = curand_uniform(&rng);
        //if( idx == 0u ) printf("//qsim.random_direction_marsaglia pidx %7lld u0 %10.4f u1 %10.4f \n", ctx.pidx, u0, u1 );
        u = 2.f*u0 - 1.f ;
        v = 2.f*u1 - 1.f ;
        b = u*u + v*v ;
    }
    while( b > 1.f ) ;

    a = 2.f*sqrtf( 1.f - b );

    dir->x = a*u ;
    dir->y = a*v ;
    dir->z = 2.f*b - 1.f ;
}



/**
qsim::rayleigh_scatter
------------------------------

Following G4OpRayleigh::PostStepDoIt

* https://bugzilla-geant4.kek.jp/show_bug.cgi?id=207 Xin Qian patch


Transverse wave nature means::

   dot(p_direction, p_polarization)  = 0
   dot(direction,   polarization)  = 0

*constant* and normalized direction retains transversality thru the candidate scatter::

    pol = p_pol + constant*dir

    dot(pol, dir) = dot(p_pol, dir) + constant* dot(dir, dir)
                  = dot(p_pol, dir) + constant* 1.
                  = dot(p_pol, dir) - dot(p_pol, dir)
                  = 0.

**/

inline QSIM_METHOD void qsim::rayleigh_scatter(RNG& rng, sctx& ctx )
{
    sphoton& p = ctx.p ;
    float3 direction ;
    float3 polarization ;

    bool looping(true) ;
    do
    {
        float u0 = curand_uniform(&rng) ;
        float u1 = curand_uniform(&rng) ;
        float u2 = curand_uniform(&rng) ;
        float u3 = curand_uniform(&rng) ;
        float u4 = curand_uniform(&rng) ;

#if !defined(PRODUCTION) && defined(DEBUG_TAG)
        stagr& tagr = ctx.tagr ;  // UNTESTED
        tagr.add(stag_sc, u0);
        tagr.add(stag_sc, u1);
        tagr.add(stag_sc, u2);
        tagr.add(stag_sc, u3);
        tagr.add(stag_sc, u4);
#endif
        float cosTheta = u0 ;
        float sinTheta = sqrtf(1.0f-u0*u0);
        if(u1 < 0.5f ) cosTheta = -cosTheta ;
        // could use uniform_sphere here : but not doing so to follow G4OpRayleigh more closely

        float sinPhi ;
        float cosPhi ;

#if defined(MOCK_CURAND ) || defined(MOCK_CUDA)
        //__sincosf(2.f*M_PIf*u2,&sinPhi,&cosPhi);   // apple extension
        float phi = 2.f*M_PIf*u2 ;
        sinPhi = sinf(phi);
        cosPhi = cosf(phi);
#else
        sincosf(2.f*M_PIf*u2,&sinPhi,&cosPhi);
#endif

        direction.x = sinTheta * cosPhi;
        direction.y = sinTheta * sinPhi;
        direction.z = cosTheta ;

        smath::rotateUz(direction, p.mom );

        float constant = -dot(direction, p.pol );

        polarization.x = p.pol.x + constant*direction.x ;
        polarization.y = p.pol.y + constant*direction.y ;
        polarization.z = p.pol.z + constant*direction.z ;

        if(dot(polarization, polarization) == 0.f )
        {

#if defined( MOCK_CURAND ) || defined(MOCK_CUDA)
            //__sincosf(2.f*M_PIf*u3,&sinPhi,&cosPhi);
            phi = 2.f*M_PIf*u3 ;
            sinPhi = sinf(phi);
            cosPhi = cosf(phi);
#else
            sincosf(2.f*M_PIf*u3,&sinPhi,&cosPhi);
#endif

            polarization.x = cosPhi ;
            polarization.y = sinPhi ;
            polarization.z = 0.f ;

            smath::rotateUz(polarization, direction);
        }
        else
        {
            // There are two directions which are perpendicular
            // to the new momentum direction
            if(u3 < 0.5f) polarization = -polarization ;
        }
        polarization = normalize(polarization);

        // simulate according to the distribution cos^2(theta)
        // where theta is the angle between old and new polarizations
        float doCosTheta = dot(polarization, p.pol ) ;
        float doCosTheta2 = doCosTheta*doCosTheta ;
        looping = doCosTheta2 < u4 ;

    } while ( looping ) ;

    p.mom = direction ;
    p.pol = polarization ;
}


/**
qsim::propagate_to_boundary
------------------------------

+---------------------+------------------+---------------------------------------------------------+-------------------------------------------------------+
| flag                |   command        |  changed                                                |  note                                                 |
+=====================+==================+=========================================================+=======================================================+
|   BULK_REEMIT       |   CONTINUE       |  time, position, direction, polarization, wavelength    | advance to reemit position with everything changed    |
+---------------------+------------------+---------------------------------------------------------+-------------------------------------------------------+
|   BULK_SCATTER      |   CONTINUE       |  time, position, direction, polarization                | advance to scatter position, new dir+pol              |
+---------------------+------------------+---------------------------------------------------------+-------------------------------------------------------+
|   BULK_ABSORB       |   BREAK          |  time, position                                         | advance to absorption position, dir+pol unchanged     |
+---------------------+------------------+---------------------------------------------------------+-------------------------------------------------------+
|   not set "SAIL"    |   BOUNDARY       |  time, position                                         | advanced to border position, dir+pol unchanged        |
+---------------------+------------------+---------------------------------------------------------+-------------------------------------------------------+


TODO:
   whilst in measurement iteration try changing the four "return"
   in the below into a single return: by setting command variable
   and returning that

**/



inline QSIM_METHOD int qsim::propagate_to_boundary(unsigned& flag, RNG& rng, sctx& ctx)
{
    sphoton& p = ctx.p ;
    const sstate& s = ctx.s ;

    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 = ctx.current_group_velocity;
    const float& distance_to_boundary = ctx.prd->q0.f.w ;


#if !defined(PRODUCTION) && defined(DEBUG_TAG)
    float u_to_sci = curand_uniform(&rng) ;  // purely for alignment with G4
    float u_to_bnd = curand_uniform(&rng) ;  // purely for alignment with G4
#endif
    float u_scattering = curand_uniform(&rng) ;
    float u_absorption = curand_uniform(&rng) ;

#if !defined(PRODUCTION) && defined(DEBUG_TAG)
    stagr& tagr = ctx.tagr ;
    tagr.add( stag_to_sci, u_to_sci);
    tagr.add( stag_to_bnd, u_to_bnd);
    tagr.add( stag_to_sca, u_scattering);
    tagr.add( stag_to_abs, u_absorption);
#endif


#if !defined(PRODUCTION) && defined(DEBUG_LOGF)
    // see notes/issues/U4LogTest_maybe_replacing_G4Log_G4UniformRand_in_Absorption_and_Scattering_with_float_version_will_avoid_deviations.rst
    float scattering_distance = -scattering_length*KLUDGE_FASTMATH_LOGF(u_scattering);
    float absorption_distance = -absorption_length*KLUDGE_FASTMATH_LOGF(u_absorption);
#else
    float scattering_distance = -scattering_length*logf(u_scattering);
    float absorption_distance = -absorption_length*logf(u_absorption);
#endif

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)

    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.propagate_to_boundary.head pidx %7lld : u_absorption %10.8f logf(u_absorption) %10.8f absorption_length %10.4f absorption_distance %10.6f \n",
        ctx.pidx, u_absorption, logf(u_absorption), absorption_length, absorption_distance );

    printf("//qsim.propagate_to_boundary.head pidx %7lld : post = np.array([%10.5f,%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx, p.pos.x, p.pos.y, p.pos.z, p.time );

    printf("//qsim.propagate_to_boundary.head pidx %7lld : distance_to_boundary %10.4f absorption_distance %10.4f scattering_distance %10.4f \n",
             ctx.pidx, distance_to_boundary, absorption_distance, scattering_distance );

    printf("//qsim.propagate_to_boundary.head pidx %7lld : u_scattering %10.4f u_absorption %10.4f \n",
             ctx.pidx, u_scattering, u_absorption  );

    }
#endif





    if (absorption_distance <= scattering_distance)
    {
        if (absorption_distance <= distance_to_boundary)
        {
            p.time += absorption_distance/group_velocity ;
            p.pos  += absorption_distance*(p.mom) ;


#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
            float absorb_time_delta = absorption_distance/group_velocity ;
            if( ctx.pidx == base->pidx )
            {
            printf("//qsim.propagate_to_boundary.body.BULK_ABSORB pidx %7lld : post = np.array([%10.5f,%10.5f,%10.5f,%10.5f]) ; absorb_time_delta = %10.8f   \n",
                    ctx.pidx, p.pos.x, p.pos.y, p.pos.z, p.time, absorb_time_delta  );

            }
#endif

            float u_reemit = reemission_prob == 0.f ? 2.f : curand_uniform(&rng);  // avoid consumption at absorption when not scintillator


#if !defined(PRODUCTION) && defined(DEBUG_TAG)
            if( u_reemit != 2.f ) tagr.add( stag_to_ree, u_reemit) ;
#endif


            if (u_reemit < reemission_prob)
            {
                float u_re_wavelength = curand_uniform(&rng);
                float u_re_mom_ph = curand_uniform(&rng);
                float u_re_mom_ct = curand_uniform(&rng);
                float u_re_pol_ph = curand_uniform(&rng);
                float u_re_pol_ct = curand_uniform(&rng);

                p.wavelength = scint->wavelength(u_re_wavelength);
                p.mom = uniform_sphere(u_re_mom_ph, u_re_mom_ct);
                p.pol = normalize(cross(uniform_sphere(u_re_pol_ph, u_re_pol_ct), p.mom));

#if !defined(PRODUCTION) && defined(DEBUG_TAG)
                tagr.add( stag_re_wl, u_re_wavelength);
                tagr.add( stag_re_mom_ph, u_re_mom_ph);
                tagr.add( stag_re_mom_ct, u_re_mom_ct);
                tagr.add( stag_re_pol_ph, u_re_pol_ph);
                tagr.add( stag_re_pol_ct, u_re_pol_ct);
#endif

                flag = BULK_REEMIT ;
                return CONTINUE;
            }
            else
            {
                flag = BULK_ABSORB ;
                return BREAK;
            }
        }
        //  otherwise sail to boundary
    }
    else
    {
        if (scattering_distance <= distance_to_boundary)
        {
            p.time += scattering_distance/group_velocity ;
            p.pos  += scattering_distance*(p.mom) ;

            rayleigh_scatter(rng, ctx);  // changes dir and pol, consumes 5u at each turn of rejection sampling loop

            flag = BULK_SCATTER;

            return CONTINUE;
        }
          //  otherwise sail to boundary
    }     // if scattering_distance < absorption_distance



    p.pos  += distance_to_boundary*(p.mom) ;
    p.time += distance_to_boundary/group_velocity   ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    float sail_time_delta = distance_to_boundary/group_velocity ;
    if( ctx.pidx == base->pidx ) printf("//qsim.propagate_to_boundary.tail.SAIL pidx %7lld : post = np.array([%10.5f,%10.5f,%10.5f,%10.5f]) ;  sail_time_delta = %10.5f   \n",
          ctx.pidx, p.pos.x, p.pos.y, p.pos.z, p.time, sail_time_delta  );
#endif

    return BOUNDARY ;
}
#endif
#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
/**
qsim::propagate_at_boundary
------------------------------------------

This was brought over from oxrap/cu/propagate.h:propagate_at_boundary_geant4_style
See env-/g4op-/G4OpBoundaryProcess.cc annotations to follow this
and compare the Opticks and Geant4 implementations.

Input:

* p.direction
* p.polarization
* s.material1.x    : refractive index
* s.material2.x    : refractive index
* prd.normal

Consumes one random deciding between BOUNDARY_REFLECT and BOUNDARY_TRANSMIT

+---------------------+------------------+-------------------------------+----------+
| output flag         |   command        |  changed                      |  note    |
+=====================+==================+===============================+==========+
|   BOUNDARY_REFLECT  |    -             |  direction, polarization      |          |
+---------------------+------------------+-------------------------------+----------+
|   BOUNDARY_TRANSMIT |    -             |  direction, polarization      |          |
+---------------------+------------------+-------------------------------+----------+

Notes:

* when geometry and refractive indices dictates TIR there is no dependence on u_reflect and always get reflection


::

                    s1
                  +----+
                   \   .   /      ^
              c1   i\  .  / r    /|\
                     \ . /        |
        material1     \./         | n
        ---------------+----------+----------
        material2      .\
                       . \
                  c2   .  \ t
                       .   \
                       +----+
                         s2

Snells law::

     s1    n2
    --- = ---         s1 n1 = s2 n2         eta = n1/n2     s1 eta = s2
     s2    n1

     s1.s1 = 1 - c1.c1   # trig identity

     s2.s2 = 1 - c2.c2

    s1 eta = s2          # snell

    s1s1 eta eta = s2s2

    ( 1.f - c1c1 ) eta eta = 1.f - c2c2

     c2c2 = 1.f - eta eta ( 1.f - c1c1 )    # snell and trig identity

Because the electromagnetic wave is transverse, the field incident onto the
interface can be decomposed into two polarization components, one P-polarized,
i.e., with the electric field vector inside the plane of incidence, and the
other one S-polarized, i.e., orthogonal to that plane.


inconsistent normal definitions, c1 is expected to be +ve and normal needs to be oriented against initial direction
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This is apparent from reflected direction vector::

      *direction + 2.0f*c1*surface_normal

The standard normal vector at an intersection position on the surface of a shape
is defined to be rigidly oriented outwards away from the shape.
This definition is used by *fill_state* in order to determine properties
of this material m1 and the next material m2 on the other side of the boundary.

The below math assumes that the photon direction is always against the normal
such that the sign of c1 is +ve. Having -ve c1 leads to non-sensical -ve TranCoeff
which results in always relecting.

So what about photons going in the other direction ?
Surface normal is used in several places in the below so presumably must
arrange to have an oriented normal that is flipped appropriately OR perhaps can change the math ?

In summary this is a case of inconsistent definitions of the normal,
that will need to be oriented ~half the time.

TODO: try avoiding "float3 oriented_normal" instead just use "bool orient"
      and multiply prd.normal by 1.f or -1.f depending on orient at every use


random aligned matching with examples/Geant/BoundaryStandalone
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* S/P/"X"-polarized + TIR + normal_incidence all now matching
* noted a 1-in-a-million deviant from TransCoeff cut edge float/double for S and P

* initially had two more deviants at very close to normal incidence that were aligned by changing
  the criteria to match Geant4 "sint1 == 0." better::

    //const bool normal_incidence = fabs(c1) > 0.999999f ;
    const bool normal_incidence = fabs(c1) == 1.f ;

* see notes/issues/QSimTest_propagate_at_boundary_vs_BoundaryStandalone_G4OpBoundaryProcessTest.rst



**Normal Incidence Special Case**

Judging normal_incidence based on absolete dot product being exactly unity "c1 == 1.f" is problematic
as when very near to normal incidence there are vectors for which the absolute dot product
is not quite 1.f but the cross product does give an exactly zero vector which gives
A_trans (nan, nan, nan) from the normalize doing : (zero,zero,zero)/zero.

Solution is to judge normal incidence based on trans_length as that is what the
calculation actually needs to be non-zero in order to be able to normalize trans to give A_trans.

However using "bool normal_incidence = trans_length == 0.f" also problematic
as it means would be using very small trans vectors to define A_trans and this
would cause a difference with double precision Geant4 and float precision Opticks.
So try using a cutoff "trans_length < 1e-6f" below which to special case normal
incidence.

**/

inline QSIM_METHOD int qsim::propagate_at_boundary(unsigned& flag, RNG& rng, sctx& ctx, float theTransmittance ) const
{
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    printf("//propagate_at_boundary.DEBUG_PIDX ctx.pidx %7lld base %p base.pidx %7lld \n", ctx.pidx, base, base->pidx  );
#endif

#if !defined(PRODUCTION) && defined(DEBUG_TAG)
    if(ctx.pidx == base->pidx)
    printf("//propagate_at_boundary.DEBUG_TAG ctx.pidx %7lld base %p base.pidx %7lld \n", ctx.pidx, base, base->pidx  );
#endif
    // stray "return 0;" left here 2024-12-14 caused : ~/j/issues/jok-tds-missing-BR-BT-on-A-side.rst

    sphoton& p = ctx.p ;
    const sstate& s = ctx.s ;

    const float& n1 = s.material1.x ;
    const float& n2 = s.material2.x ;
    const float eta = n1/n2 ;

    const float3* normal = (float3*)&ctx.prd->q0.f.x ;     // geometrical outwards normal

    const float _c1 = -dot(p.mom, *normal );                // _c1 : cos(angle_of_incidence) not yet oriented
    const float3 oriented_normal = _c1 < 0.f ? -(*normal) : (*normal) ; // oriented against incident p.mom
    const float3 trans = cross(p.mom, oriented_normal) ;   // perpendicular to plane of incidence, S-pol direction
    const float trans_length = length(trans) ;             // same as sin(theta), as p.mom and oriented_normal are unit vectors
    const bool normal_incidence = trans_length < 1e-6f  ;  // p.mom parallel/anti-parallel to oriented_normal
    const float3 A_trans = normal_incidence ? p.pol : trans/trans_length ; // normalized unit vector : perpendicular to plane of incidence
    const float E1_perp = dot(p.pol, A_trans);     // amplitude of polarization in direction perpendicular to plane of incidence, ie S polarization

    const float c1 = fabs(_c1) ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.propagate_at_boundary.head pidx %7lld : theTransmittance = %10.8f \n", ctx.pidx, theTransmittance  );
    printf("//qsim.propagate_at_boundary.head pidx %7lld : nrm = np.array([%10.8f,%10.8f,%10.8f]) ; lnrm = %10.8f  \n",
         ctx.pidx, oriented_normal.x, oriented_normal.y, oriented_normal.z, length(oriented_normal) );
    printf("//qsim.propagate_at_boundary.head pidx %7lld : pos = np.array([%10.5f,%10.5f,%10.5f]) ; lpos = %10.8f \n",
         ctx.pidx, p.pos.x, p.pos.y, p.pos.z, length(p.pos) );
    printf("//qsim.propagate_at_boundary.head pidx %7lld : mom0 = np.array([%10.8f,%10.8f,%10.8f]) ; lmom0 = %10.8f \n",
         ctx.pidx, p.mom.x, p.mom.y, p.mom.z, length(p.mom)  );
    printf("//qsim.propagate_at_boundary.head pidx %7lld : pol0 = np.array([%10.8f,%10.8f,%10.8f]) ; lpol0 = %10.8f \n",
          ctx.pidx, p.pol.x, p.pol.y, p.pol.z, length(p.pol)  );
    printf("//qsim.propagate_at_boundary.head pidx %7lld : n1,n2,eta = (%10.8f,%10.8f,%10.8f) \n", ctx.pidx, n1, n2, eta );
    printf("//qsim.propagate_at_boundary.head pidx %7lld : c1 = %10.8f ; normal_incidence = %d \n", ctx.pidx, c1, normal_incidence );
    }
#endif

    const float c2c2 = 1.f - eta*eta*(1.f - c1 * c1 ) ;   // Snells law and trig identity
    bool tir = c2c2 < 0.f ;
    const float EdotN = dot(p.pol, oriented_normal ) ;  // used for TIR polarization
    const float c2 = tir ? 0.f : sqrtf(c2c2) ;   // c2 chosen +ve, set to 0.f for TIR => reflection_coefficient = 1.0f : so will always reflect
    const float n1c1 = n1*c1 ;
    const float n2c2 = n2*c2 ;
    const float n2c1 = n2*c1 ;
    const float n1c2 = n1*c2 ;

    const float2 E1   = normal_incidence ? make_float2( 0.f, 1.f) : make_float2( E1_perp , length( p.pol - (E1_perp*A_trans) ) );
    const float2 E2_t = make_float2(  2.f*n1c1*E1.x/(n1c1+n2c2), 2.f*n1c1*E1.y/(n2c1+n1c2) ) ;  // ( S:perp, P:parl )
    const float2 E2_r = make_float2( E2_t.x - E1.x             , (n2*E2_t.y/n1) - E1.y     ) ;  // ( S:perp, P:parl )
    const float2 RR = normalize(E2_r) ;
    const float2 TT = normalize(E2_t) ;
    const float TransCoeff = theTransmittance >= 0.f ?
                                                           theTransmittance
                                                     :
                                                           ( tir || n1c1 == 0.f ? 0.f : n2c2*dot(E2_t,E2_t)/n1c1 )
                                                     ;

    /*
    E1, E2_t, E2_t: incident, transmitted and reflected amplitudes in S and P directions
    RR, TT: normalized
    */

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.propagate_at_boundary.body pidx %7lld : TransCoeff = %10.8f ; n1c1 = %10.8f ; n2c2 = %10.8f \n",
            ctx.pidx, TransCoeff, n1c1, n2c2 );

    printf("//qsim.propagate_at_boundary.body pidx %7lld : E2_t = np.array([%10.8f,%10.8f]) ; lE2_t = %10.8f \n",
            ctx.pidx,  E2_t.x, E2_t.y, length(E2_t) );

    printf("//qsim.propagate_at_boundary.body pidx %7lld : A_trans = np.array([%10.8f,%10.8f,%10.8f]) ; lA_trans = %10.8f \n",
            ctx.pidx,  A_trans.x, A_trans.y, A_trans.z, length(A_trans) );

    }
#endif


#if !defined(PRODUCTION) && defined(DEBUG_TAG)
    const float u_boundary_burn = curand_uniform(&rng) ;  // needed for random consumption alignment with Geant4 G4OpBoundaryProcess::PostStepDoIt
#endif
    const float u_reflect = curand_uniform(&rng) ;
    bool reflect = u_reflect > TransCoeff  ;

#if !defined(PRODUCTION) && defined(DEBUG_TAG)
    stagr& tagr = ctx.tagr ;
    tagr.add( stag_at_burn_sf_sd, u_boundary_burn);
    tagr.add( stag_at_ref,  u_reflect);
#endif

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.propagate_at_boundary.body pidx %7lld : u_reflect %10.4f TransCoeff %10.4f reflect %d \n",
              ctx.pidx,  u_reflect, TransCoeff, reflect   );

    printf("//qsim.propagate_at_boundary.body pidx %7lld : mom0 = np.array([%10.8f,%10.8f,%10.8f]) ; lmom0 = %10.8f \n",
               ctx.pidx, p.mom.x, p.mom.y, p.mom.z, length(p.mom) ) ;

    printf("//qsim.propagate_at_boundary.body pidx %7lld : pos = np.array([%10.5f,%10.5f,%10.5f]) ; lpos = %10.8f \n",
               ctx.pidx, p.pos.x, p.pos.y, p.pos.z, length(p.pos)  );

    printf("//qsim.propagate_at_boundary.body pidx %7lld : nrm = np.array([%10.8f,%10.8f,%10.8f]) ; lnrm = %10.8f \n",
               ctx.pidx, oriented_normal.x, oriented_normal.y, oriented_normal.z, length(oriented_normal) ) ;

    printf("//qsim.propagate_at_boundary.body pidx %7lld : n1 = %10.8f ; n2 = %10.8f ; eta = %10.8f  \n",
               ctx.pidx, n1, n2, eta );

    printf("//qsim.propagate_at_boundary.body pidx %7lld : c1 = %10.8f ; eta_c1 = %10.8f ; c2 = %10.8f ; eta_c1__c2 = %10.8f \n",
               ctx.pidx, c1, eta*c1, c2, (eta*c1 - c2) );

    }
#endif

    p.mom = reflect
                    ?
                       p.mom + 2.0f*c1*oriented_normal
                    :
                       eta*(p.mom) + (eta*c1 - c2)*oriented_normal
                    ;


    // Q: Does the new p.mom need to be normalized ?
    // A: NO, it is inherently normalized as derived in the comment below


    const float3 A_paral = normalize(cross(p.mom, A_trans));   // new P-pol direction

    p.pol =  normal_incidence ?
                                         ( reflect ?  p.pol*(n2>n1? -1.f:1.f) : p.pol )
                                      :
                                         ( reflect ?
                                                   ( tir ?  -p.pol + 2.f*EdotN*oriented_normal : RR.x*A_trans + RR.y*A_paral )

                                                   :
                                                       TT.x*A_trans + TT.y*A_paral

                                                   )
                                      ;



     // Q: Above expression kinda implies A_trans and A_paral are same for reflect and transmit ?
     // A: NO IT DOESNT,
     //    A_trans is the same (except for normal incidence) as there is only one perpendicular
     //    to the plane of incidence which is the same for i,r,t.
     //
     //    A_paral depends on the new p.mom (is has to be orthogonal to p.mom and A_trans)
     //    and p.mom of course is different for r and t
     //    (the reflect bool is used in multiple places, not just here)



#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.propagate_at_boundary.tail pidx %7lld : reflect %d tir %d TransCoeff %10.4f u_reflect %10.4f \n", ctx.pidx, reflect, tir, TransCoeff, u_reflect );
    printf("//qsim.propagate_at_boundary.tail pidx %7lld : mom1 = np.array([%10.8f,%10.8f,%10.8f]) ; lmom1 = %10.8f  \n",
        ctx.pidx, p.mom.x, p.mom.y, p.mom.z, length(p.mom) );
    printf("//qsim.propagate_at_boundary.tail pidx %7lld : pol1 = np.array([%10.8f,%10.8f,%10.8f]) ; lpol1 = %10.8f \n",
        ctx.pidx, p.pol.x, p.pol.y, p.pol.z, length(p.pol) );
    }

    /*
    if(ctx.pidx == 251959)
    {
        printf("//qsim.propagate_at_boundary RR.x %10.4f A_trans (%10.4f %10.4f %10.4f )  RR.y %10.4f  A_paral (%10.4f %10.4f %10.4f ) \n",
              RR.x, A_trans.x, A_trans.y, A_trans.z,
              RR.y, A_paral.x, A_paral.y, A_paral.z );

        printf("//qsim.propagate_at_boundary reflect %d  tir %d polarization (%10.4f, %10.4f, %10.4f) \n", reflect, tir, p.pol.x, p.pol.y, p.pol.z );
    }
    */
#endif

    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 )
    {
        const float u_br_align_0 = curand_uniform(&rng) ;
        const float u_br_align_1 = curand_uniform(&rng) ;
        const float u_br_align_2 = curand_uniform(&rng) ;
        const float u_br_align_3 = curand_uniform(&rng) ;

        // switched below standard tags from stag_br_align_0/1/2/3 to simplify A:tag to B:stack mapping
        tagr.add( stag_to_sci, u_br_align_0 );
        tagr.add( stag_to_bnd, u_br_align_1 );
        tagr.add( stag_to_sca, u_br_align_2 );
        tagr.add( stag_to_abs, u_br_align_3 );
    }
#endif

    return CONTINUE ;
}
/**
qsim::propagate_at_boundary_with_T
------------------------------------

Variant of qsim::propagate_at_boundary where a value of theTransmittance
is passed in as an argument (eg from CustomART calculation) rather then
calculated here.

HMM: was hoping that passing in theTransmittance would afford some
simplifications, but it seems almost no simplification is possible
as need to calculate everything anyhow for the polarization.

Given that this is almost exactly the same as qsim::propagate_at_boundary
and no simplification is afforded, might as well just keep using that
and leave this just for testing.

**/

inline QSIM_METHOD int qsim::propagate_at_boundary_with_T(unsigned& flag, RNG& rng, sctx& ctx, float theTransmittance ) const
{
    sphoton& p = ctx.p ;
    const sstate& s = ctx.s ;

    const float& n1 = s.material1.x ;
    const float& n2 = s.material2.x ;
    const float eta = n1/n2 ;

    const float3* normal = (float3*)&ctx.prd->q0.f.x ;     // geometrical outwards normal

    const float _c1 = -dot(p.mom, *normal );                // _c1 : cos(angle_of_incidence) not yet oriented
    const float3 oriented_normal = _c1 < 0.f ? -(*normal) : (*normal) ; // oriented against incident p.mom
    const float3 trans = cross(p.mom, oriented_normal) ;   // perpendicular to plane of incidence, S-pol direction
    const float trans_length = length(trans) ;             // same as sin(theta), as p.mom and oriented_normal are unit vectors
    const bool normal_incidence = trans_length < 1e-6f  ;  // p.mom parallel/anti-parallel to oriented_normal
    const float3 A_trans = normal_incidence ? p.pol : trans/trans_length ; // normalized unit vector : perpendicular to plane of incidence
    const float E1_perp = dot(p.pol, A_trans);     // amplitude of polarization in direction perpendicular to plane of incidence, ie S polarization

    const float c1 = fabs(_c1) ;

    const float c2c2 = 1.f - eta*eta*(1.f - c1 * c1 ) ;   // Snells law and trig identity
    bool tir = c2c2 < 0.f ;
    const float EdotN = dot(p.pol, oriented_normal ) ;  // used for TIR polarization
    const float c2 = tir ? 0.f : sqrtf(c2c2) ;   // c2 chosen +ve, set to 0.f for TIR => reflection_coefficient = 1.0f : so will always reflect
    const float n1c1 = n1*c1 ;
    const float n2c2 = n2*c2 ;
    const float n2c1 = n2*c1 ;
    const float n1c2 = n1*c2 ;

    const float2 E1   = normal_incidence ? make_float2( 0.f, 1.f) : make_float2( E1_perp , length( p.pol - (E1_perp*A_trans) ) );
    const float2 E2_t = make_float2(  2.f*n1c1*E1.x/(n1c1+n2c2), 2.f*n1c1*E1.y/(n2c1+n1c2) ) ;  // ( S:perp, P:parl )
    const float2 E2_r = make_float2( E2_t.x - E1.x             , (n2*E2_t.y/n1) - E1.y     ) ;  // ( S:perp, P:parl )
    const float2 RR = normalize(E2_r) ;
    const float2 TT = normalize(E2_t) ;


/*
    const float TransCoeff = theTransmittance >= 0.f ?
                                                           theTransmittance
                                                     :
                                                           ( tir || n1c1 == 0.f ? 0.f : n2c2*dot(E2_t,E2_t)/n1c1 )
                                                     ;
*/

    const float& TransCoeff = theTransmittance ;

    const float u_reflect = curand_uniform(&rng) ;
    bool reflect = u_reflect > TransCoeff  ;

    p.mom = reflect
                    ?
                       p.mom + 2.0f*c1*oriented_normal
                    :
                       eta*(p.mom) + (eta*c1 - c2)*oriented_normal
                    ;


    const float3 A_paral = normalize(cross(p.mom, A_trans));   // new P-pol direction

    p.pol =  normal_incidence ?
                                         ( reflect ?  p.pol*(n2>n1? -1.f:1.f) : p.pol )
                                      :
                                         ( reflect ?
                                                   ( tir ?  -p.pol + 2.f*EdotN*oriented_normal : RR.x*A_trans + RR.y*A_paral )

                                                   :
                                                       TT.x*A_trans + TT.y*A_paral

                                                   )
                                      ;

    flag = reflect ? BOUNDARY_REFLECT : BOUNDARY_TRANSMIT ;


    return CONTINUE ;
}

#endif








/**
Reflected momentum vector
---------------------------

::

     p.mom(new) = p.mom(old) + 2.f*c1*oriented_normal


     PdotN = dot(p.mom(old), oriented_normal) = -c1

     "c1" is +ve, above dot product -ve as incident vector is against the normal

     p.mom(new) = p.mom(old) -2.f dot(p.mom(old), oriented_normal) * oriented_normal

     NewMomentum = OldMomentum - 2*PdotN*oriented_normal



     OA = oriented_normal

     IO = OC = p.mom (old)
          OR = p.mom (new)

          OR = OC + CD + DR  = OC + 2.f*DR
          CD = DR = c1*oriented_normal


                  s1   s1
                I----A----R
                .\   .   /.
            c1  . \  .  / .  c1
                .  \ . /  .
                .   \./   .
       ---------+----O----D------
                :    :.   :
                :    : .  :  c1
           c1   :    :  . :
                :    :   .:
                +....B....C
                 s1    s1


* NB that is inherently normalized, as demonstrated in the below reference


Transmitted momentum vector
-------------------------------

* https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
* ~/opticks_refs/bram_de_greve_reflection_refraction_in_ray_tracing.pdf


::

     eta = n1/n2

     p.mom(new) =  eta*p.mom(old) + (eta*c1 - c2)*oriented_normal


Consider the perp and para components of the incident and transmitted mom vectors::

    i = i_perp + i_para
    t = t_perp + t_para

    i_perp = ( i.n ) n      (component of i in normal direction)
    t_perp = ( t.n ) n      (component of t in normal direction )

    i_para = i - (i.n) n
    t_para = t - (t.n) n

    # HMM caution with backwards oriented normal direction convention
    #     and dot product signs


Snell::

    n1 s1 = n2 s2

    eta s1 = s2      (eta = n1/n2)

    s2 = eta s1

    |i_para| = s1
    |t_para| = s2

    t_para  = eta i_para     (as they point in same direction)

    t_para  = eta ( i - (i.n) n )
            = eta ( i + c1 n )

    t_perp  = - sqrt( 1-|t_para|^2 ) n     (-ve as opposite direction to oriented normal)
            = - c2 n

    t   = eta i + ( eta c1 - c2 ) n


Now is that inherently normalized ? YES

   t.t =  eta*eta*i.i  + (eta c1 - c2 )*(eta c1 - c2 )*n.n + 2 eta (eta c1 - c2) i.n

         (i.i = 1, n.n = 1, i.n = -c1 )

       = eta*eta + (eta c1 - c2)*(eta c1 - c2) + 2 eta(eta c1 - c2) (-c1 )

       = eta eta  +  eta eta c1 c1 + c2 c2 - 2 eta c1 c2  - 2 eta eta c1 c1 + 2 eta c2 c1
                                           --------------                    --------------
       =   eta eta  ( 1 - c1 c1 ) + c2 c2

       = 1.f      ( form above   c2c2 = 1.f - eta eta ( 1.f - c1c1 ) )    # snell and trig identity



Compare transmitted vector with G4OpBoundaryProcess::DielectricDielectric
----------------------------------------------------------------------------

::

    1226                 theStatus = FresnelRefraction;
    1227
    1228                 if (sint1 > 0.0) {      // incident ray oblique
    1229
    1230                    alpha = cost1 - cost2*(Rindex2/Rindex1);
    1231                    NewMomentum = OldMomentum + alpha*theFacetNormal;
    1232                    NewMomentum = NewMomentum.unit();


    t = eta i  + (eta c1 - c2 ) n      eta = n1/n2
    t/eta = i + (c1 - c2/eta ) n

    Because Geant4 normalizes NewMomentum it gets away with
    playing fast and loose with factors of 1/eta


**/







/*
G4OpBoundaryProcess::DielectricDielectric


1152               if (sint1 > 0.0) {
1153                  A_trans = OldMomentum.cross(theFacetNormal);
1154                  A_trans = A_trans.unit();
1155                  E1_perp = OldPolarization * A_trans;
1156                  E1pp    = E1_perp * A_trans;
1157                  E1pl    = OldPolarization - E1pp;
1158                  E1_parl = E1pl.mag();
1159               }
1160               else {
1161                  A_trans  = OldPolarization;
1162                  // Here we Follow Jackson's conventions and we set the
1163                  // parallel component = 1 in case of a ray perpendicular
1164                  // to the surface
1165                  E1_perp  = 0.0;
1166                  E1_parl  = 1.0;
1167               }
1168
1169               s1 = Rindex1*cost1;
1170               E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
1171               E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
1172               E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1173               s2 = Rindex2*cost2*E2_total;
1174
1175               if (theTransmittance > 0) TransCoeff = theTransmittance;
1176               else if (cost1 != 0.0) TransCoeff = s2/s1;
1177               else TransCoeff = 0.0;


reflect

1217
1218                        E2_parl   = Rindex2*E2_parl/Rindex1 - E1_parl;
1219                        E2_perp   = E2_perp - E1_perp;
1220                        E2_total  = E2_perp*E2_perp + E2_parl*E2_parl;
1221                        A_paral   = NewMomentum.cross(A_trans);
1222                        A_paral   = A_paral.unit();
1223                        E2_abs    = std::sqrt(E2_total);
1224                        C_parl    = E2_parl/E2_abs;
1225                        C_perp    = E2_perp/E2_abs;
1226
1227                        NewPolarization = C_parl*A_paral + C_perp*A_trans;
1228

transmit

1253                    alpha = cost1 - cost2*(Rindex2/Rindex1);
1254                    NewMomentum = OldMomentum + alpha*theFacetNormal;
1255                    NewMomentum = NewMomentum.unit();
1256 //                   PdotN = -cost2;
1257                    A_paral = NewMomentum.cross(A_trans);
1258                    A_paral = A_paral.unit();

1259                    E2_abs     = std::sqrt(E2_total);
1260                    C_parl     = E2_parl/E2_abs;
1261                    C_perp     = E2_perp/E2_abs;
1262
1263                    NewPolarization = C_parl*A_paral + C_perp*A_trans;

*/


/*
qsim::propagate_at_surface_MultiFilm
-------------------------------

based on https://juno.ihep.ac.cn/trac/browser/offline/trunk/Simulation/DetSimV2/PMTSim/src/junoPMTOpticalModel.cc
which follows a flawed approach to polarization

Rs: s-component reflect probability
Ts: s-component transmit probability
Rp: p-component reflect probability
Tp: p-component reflect probability

TODO:
   access the qe
   access the multifilm catagory
   optical photon polarization

*/

#if defined(__CUDACC__) || defined(__CUDABE__)
inline QSIM_METHOD int qsim::propagate_at_surface_MultiFilm(unsigned& flag, RNG& rng, sctx& ctx )
{

	const float one = 1.0f;
    const sphoton& p = ctx.p ;
    const float3* normal = (float3*)&ctx.prd->q0.f.x ;
    int   lpmtid = ctx.prd->identity() - 1 ;  // identity comes from optixInstance.instanceId where 0 means not-a-sensor

    float minus_cos_theta = dot(p.mom, *normal);
    int pmtcat = pmt->get_lpmtcat_from_lpmtid(lpmtid);
    float wv_nm = p.wavelength;

    float4 RsTsRpTp = multifilm->lookup(pmtcat, wv_nm, minus_cos_theta);


    const float c1 = fabs(minus_cos_theta) ;
    const float s1 = sqrtf(one -c1*c1);

    float EsEs = s1 > 0.f ? dot(p.pol, cross( p.mom, *normal))/s1 : 0.f ;
    EsEs *= EsEs;   //   orienting normal doesnt matter as squared : this is S_vs_P power fraction


    float3 ART ;
    ART.z = RsTsRpTp.y*EsEs + RsTsRpTp.w*(one - EsEs);
    ART.y = RsTsRpTp.x*EsEs + RsTsRpTp.z*(one - EsEs);
    ART.x = one - (ART.y+ART.z);

    const float& A = ART.x ;
    const float& T = ART.z ;


    float4 RsTsRpTpNormal = multifilm->lookup(pmtcat, wv_nm, -one );
    // Normal means the photon incident from glass to vacuum, AOI = 0 deg  cos_theta = -1.f

    float3 ART_normal;
    ART_normal.z = 0.5f*(RsTsRpTpNormal.y + RsTsRpTpNormal.w); // T:0.5f*(Ts+Tp)
    ART_normal.y = 0.5f*(RsTsRpTpNormal.x + RsTsRpTpNormal.z); // R:0.5f*(Rs+Rp)
    ART_normal.x = one -(ART_normal.y + ART_normal.z) ;        // 1.f - (R+T)

    const float& An = ART_normal.x ;
    const float energy_eV = qpmt<float>::hc_eVnm/wv_nm ;
    const float qe_scale = pmt->get_qescale_from_lpmtid(lpmtid);
    const float qe_shape = pmt->get_lpmtcat_qe(pmtcat, energy_eV);

    const float _qe = minus_cos_theta > 0.f ? 0.f : qe_scale * qe_shape;

    const float& theAbsorption = A;
    const float& theTransmittance = T/(one-A);
    const float& theEfficiency = _qe/An;

    float u_theAbsorption = curand_uniform(&rng);
    int action = u_theAbsorption < theAbsorption  ? BREAK : CONTINUE ;

    if( action == BREAK )
    {
        float u_theEfficiency = curand_uniform(&rng) ;
        flag = u_theEfficiency < theEfficiency ? SURFACE_DETECT : SURFACE_ABSORB ;
    }
    else
    {
        propagate_at_boundary( flag, rng, ctx, theTransmittance  );
    }

    //printf("//qsim.propagate_at_surface_MultiFilm pidx %7lld lpmtid %d ART ( %7.3f %7.3f %7.3f ) u_theAbsorption  %7.3f action %d \n",
    //ctx.pidx, lpmtid, ART.x, ART.y, ART.z, u_theAbsorption, action);

    return action ;

}

#endif


#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)

/**
qsim::propagate_at_surface   (HMM: perhaps propagate_at_simplified_surface )
-------------------------------------------------------------------------------

+---------------------+------------------+---------------------------------------------------------+--------------+
| output flag         |   command        |  changed                                                |  note        |
+=====================+==================+=========================================================+==============+
|   SURFACE_ABSORB    |    BREAK         |                                                         |              |
+---------------------+------------------+---------------------------------------------------------+--------------+
|   SURFACE_DETECT    |    BREAK         |                                                         |              |
+---------------------+------------------+---------------------------------------------------------+--------------+
|   SURFACE_DREFLECT  |    CONTINUE      |   direction, polarization                               |              |
+---------------------+------------------+---------------------------------------------------------+--------------+
|   SURFACE_SREFLECT  |    CONTINUE      |   direction, polarization                               |              |
+---------------------+------------------+---------------------------------------------------------+--------------+

Action depens on s.surface float4 params::

+---------------+---------------------+----------------------------------+
| s.surface.x   | detect fraction     |                                  |
+---------------+---------------------+----------------------------------+
| s.surface.y   | absorb fraction     |                                  |
+---------------+---------------------+----------------------------------+
| s.surface.z   | reflect specular    | NOT USED AS ALL FOUR SUM TO 1.   |
+---------------+---------------------+----------------------------------+
| s.surface.w   | reflect diffuse     |                                  |
+---------------+---------------------+----------------------------------+

Excluding technical alignment burns a single u_surface random is throw
to decide on what happens at the surface according to the float4 probabilities
that are assumed to sum to unity::

   +-------------------+-------------------+---------------------+------------------------+
   |                   |                   |                     |                        |
   |  absorb           |  detect           |  reflect_diffuse_   |  "reflect_specular_"   |
   0--SURFACE_ABSORB---|--SURFACE_DETECT---|--SURFACE_DREFLECT---|--SURFACE_SREFLECT------1
   |  s.surface.y      |  s.surface.x      |  s.surface.w        |  s.surface.z           |
   |                   |                   |                     |                        |
   +-------------------+-------------------+---------------------+------------------------+

Refs::

    GSurfaceLib::propertyName
    GSurfaceLib::createStandardSurface

Different types of surface are modelled by changing the four values,
in such a way that they always sum to unity::

    In [20]: t.oldsur.shape
    Out[20]: (46, 2, 761, 4)

    In [21]: t.oldsur[:,0,:].shape  ## pick payload group 0
    Out[21]: (46, 761, 4)

    assert np.all( np.sum(t.oldsur[:,0,:], axis=2) == 1. )

The s.surface float4 is filled by qbnd::fill_state via::

    const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR ;
    s.surface = boundary_lookup(wavelength,su_line,0) ;
    s.optical = optical[su_line].u ;

**/

inline QSIM_METHOD int qsim::propagate_at_surface(unsigned& flag, RNG& rng, sctx& ctx)
{
    const sstate& s = ctx.s ;
    const float& detect = s.surface.x ;
    const float& absorb = s.surface.y ;
    //const float& reflect_specular_ = s.surface.z ;
    const float& reflect_diffuse_  = s.surface.w ;

    float u_surface = curand_uniform(&rng);

#if !defined(PRODUCTION) && defined(DEBUG_TAG)
    stagr& tagr = ctx.tagr ;
    float u_surface_burn = curand_uniform(&rng);
    tagr.add( stag_at_burn_sf_sd, u_surface);
    tagr.add( stag_sf_burn,       u_surface_burn);
#endif


    int action = u_surface < absorb + detect ? BREAK : CONTINUE  ;

    if( action == BREAK )
    {
#if defined(WITH_CUSTOM4)
        int pmtid = ctx.prd->identity() - 1 ;     // identity comes from optixInstance.instanceId where 0 means not-a-sensor
        float qe = 1.f ;
        float u_qe = curand_uniform(&rng);
        if( s_pmt::is_spmtid(pmtid) )
        {
            const float energy_eV = qpmt<float>::hc_eVnm/ctx.p.wavelength ;
            float qe_shape = pmt->s_qeshape_prop->interpolate( 0, energy_eV );
            float qe_scale = pmt->get_s_qescale_from_spmtid(pmtid);
            qe = qe_shape*qe_scale ;
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
            if(ctx.pidx == base->pidx)
            printf("//qsim.propagate_at_surface.BREAK.is_spmtid pidx %7lld : pmtid %d energy_eV %7.3f qe_shape %7.3f qe_scale %7.3f qe %7.3f detect %7.3f absorb %7.3f reflect_specular %7.3f reflect_diffuse %7.3f \n" ,
                ctx.pidx, pmtid, energy_eV, qe_shape, qe_scale, qe, detect, absorb, s.surface.z, reflect_diffuse_ );
#endif
        }
        flag = u_surface < absorb ?
                                      SURFACE_ABSORB
                                  :
                                      ( u_qe < qe  ? EFFICIENCY_COLLECT : EFFICIENCY_CULL  )
                                  ;
#else
        flag = u_surface < absorb ?
                                      SURFACE_ABSORB
                                  :
                                      SURFACE_DETECT
                                  ;
#endif

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
        if(ctx.pidx == base->pidx)
        printf("//qsim.propagate_at_surface.SA/SD.BREAK pidx %7lld : flag %d \n" , ctx.pidx, flag );
#endif
    }
    else
    {
        flag = u_surface < absorb + detect + reflect_diffuse_ ?  SURFACE_DREFLECT : SURFACE_SREFLECT ;
        switch(flag)
        {
            case SURFACE_DREFLECT: reflect_diffuse( rng, ctx)  ; break ;
            case SURFACE_SREFLECT: reflect_specular(rng, ctx)  ; break ;
        }
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
        if(ctx.pidx == base->pidx)
        printf("//qsim.propagate_at_surface.DR/SR.CONTINUE pidx %7lld : flag %d \n" , ctx.pidx, flag );
#endif
    }
    return action ;
}


inline QSIM_METHOD int qsim::propagate_at_surface_Detect(unsigned& flag, RNG& rng, sctx& ctx) const
{
    float u_surface_burn = curand_uniform(&rng);  // for random alignment
    flag = SURFACE_DETECT ;
    return BREAK ;
}


#if defined(WITH_CUSTOM4)

/**
qsim::propagate_at_surface_CustomART
-------------------------------------

lpmtid:-1
   indicates "not-a-sensor", that occurring at this juncture would
   indicate a sensor_identity issue as CustomART special surfaces
   are expected to always have sensor identities


Where ctx.prd->identity() comes from ? Where is the "+ 1" done ?
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

1. SBT::collectInstances sets OptixInstance::instanceId from sqat4::get_IAS_OptixInstance_instanceId
   aka "sensor_identifier"

2. CSGFoundry::addInstance for firstcall:true does the "+1" as done by CSGFoundry::addInstanceVector

3. original access to the copyno from Geant4 in U4SensorIdentifierDefault::getInstanceIdentity
   which is used from U4Tree::identifySensitiveInstances to populate the stree.h snode::sensor_id

**/

inline QSIM_METHOD int qsim::propagate_at_surface_CustomART(unsigned& flag, RNG& rng, sctx& ctx) const
{

    sphoton& p = ctx.p ;
    const float3* normal = (float3*)&ctx.prd->q0.f.x ;  // geometrical outwards normal
    int lpmtid = ctx.prd->identity() - 1 ;  // identity comes from optixInstance.instanceId where 0 means not-a-sensor
    const float lposcost = ctx.prd->lposcost() ;  // local frame intersect position cosine theta


    float minus_cos_theta = dot(p.mom, *normal);
    float dot_pol_cross_mom_nrm = dot(p.pol,cross(p.mom,*normal)) ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if( ctx.pidx_debug )
    {
    float3 cross_mom_nrm = cross(p.mom, *normal) ;
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld : mom = np.array([%10.8f,%10.8f,%10.8f]) ; lmom = %10.8f \n",
       ctx.pidx, p.mom.x, p.mom.y, p.mom.z, length(p.mom) );
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld : pol = np.array([%10.8f,%10.8f,%10.8f]) ; lpol = %10.8f \n",
       ctx.pidx, p.pol.x, p.pol.y, p.pol.z, length(p.pol) );
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld : nrm = np.array([%10.8f,%10.8f,%10.8f]) ; lnrm = %10.8f \n",
       ctx.pidx, normal->x, normal->y, normal->z, length(*normal) );
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld : cross_mom_nrm = np.array([%10.8f,%10.8f,%10.8f]) ; lcross_mom_nrm = %10.8f  \n",
           ctx.pidx, cross_mom_nrm.x, cross_mom_nrm.y, cross_mom_nrm.z, length(cross_mom_nrm)  );
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld : dot_pol_cross_mom_nrm = %10.8f \n", ctx.pidx, dot_pol_cross_mom_nrm );
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld : minus_cos_theta = %10.8f \n", ctx.pidx, minus_cos_theta );
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld : lposcost = %10.8f (expect 0->1)\n", ctx.pidx, lposcost );
    }
#endif

    // formerly excluded Custom4 hits onto WP PMTs see ~/j/issues/jok-tds-mu-running-NOT-A-SENSOR-warnings.rst
    //if(lpmtid < s_pmt::OFFSET_CD_LPMT || lpmtid >= s_pmt::OFFSET_WP_PMT_END )
    //if(lpmtid < s_pmt::OFFSET_CD_LPMT || lpmtid >= s_pmt::OFFSET_WP_ATM_LPMT_END )
    if(lpmtid < s_pmt::OFFSET_CD_LPMT || lpmtid >=   s_pmt::OFFSET_WP_WAL_PMT_END )
    {
        flag = NAN_ABORT ;
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
        printf("//qsim::propagate_at_surface_CustomART pidx %7lld lpmtid %d : ERROR UNEXPECTED LPMTID : NAN_ABORT \n", ctx.pidx, lpmtid );
#endif
        return BREAK ;
    }

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if( ctx.pidx_debug )
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld lpmtid %d wl %8.4f mct %8.4f dpcmn %8.4f pmt %p pre-ATQC \n",
           ctx.pidx, lpmtid, p.wavelength, minus_cos_theta, dot_pol_cross_mom_nrm, pmt );
#endif

    float ATQC[4] = {} ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(lpmtid > -1 && pmt != nullptr) pmt->get_lpmtid_ATQC(ATQC, lpmtid, p.wavelength, minus_cos_theta, dot_pol_cross_mom_nrm, lposcost, ctx.pidx, ctx.pidx_debug );
#else
    if(lpmtid > -1 && pmt != nullptr) pmt->get_lpmtid_ATQC(ATQC, lpmtid, p.wavelength, minus_cos_theta, dot_pol_cross_mom_nrm, lposcost );
#endif



#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if( ctx.pidx_debug )
    printf("//qsim::propagate_at_surface_CustomART pidx %7lld lpmtid %d wl %8.4f mct %8.4f dpcmn %8.4f lpc %8.4f ATQC ( %8.4f %8.4f %8.4f %8.4f  ) \n",
           ctx.pidx, lpmtid, p.wavelength, minus_cos_theta, dot_pol_cross_mom_nrm, lposcost, ATQC[0], ATQC[1], ATQC[2], ATQC[3] );
#endif


    const float& theAbsorption        = ATQC[0];
    const float& theTransmittance     = ATQC[1];
    const float& theEfficiency        = ATQC[2];
    const float& collectionEfficiency = ATQC[3];

    float u_theAbsorption = curand_uniform(&rng);
    int action = u_theAbsorption < theAbsorption  ? BREAK : CONTINUE ;


#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if( ctx.pidx_debug )
        printf("//qsim.propagate_at_surface_CustomART pidx %7lld lpmtid %d ATQC ( %8.4f %8.4f %8.4f %8.4f ) u_theAbsorption  %7.3f action %d \n",
        ctx.pidx, lpmtid, ATQC[0], ATQC[1], ATQC[2], ATQC[3], u_theAbsorption, action  );
#endif

    if( action == BREAK )
    {
        float u_theEfficiency = curand_uniform(&rng) ;
        float u_collectionEfficiency = curand_uniform(&rng);

        flag = u_theEfficiency < theEfficiency ?
                                                    (  u_collectionEfficiency < collectionEfficiency ? EFFICIENCY_COLLECT : EFFICIENCY_CULL )
                                               :
                                                    SURFACE_ABSORB
                                               ;

        // former SD:SURFACE_DETECT, now becomes EC:EFFICIENCY_COLLECT or EX:EFFICIENCY_CULL depending on collectionEfficiency and random throw

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
        if( ctx.pidx_debug )
            printf("//qsim.propagate_at_surface_CustomART.BREAK.SD/SA EC/EX pidx %7lld lpmtid %d ATQC ( %8.4f %8.4f %8.4f %8.4f ) u_theEfficiency  %8.4f theEfficiency %8.4f flag %d \n",
                                                                    ctx.pidx, lpmtid, ATQC[0],ATQC[1], ATQC[2],ATQC[3],  u_theEfficiency,  theEfficiency, flag );
#endif

    }
    else
    {

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
        if( ctx.pidx_debug )
            printf("//qsim.propagate_at_surface_CustomART.CONTINUE pidx %7lld lpmtid %d ATQC ( %7.3f %7.3f %7.3f %7.3f ) theTransmittance %7.3f  \n",
            ctx.pidx, lpmtid, ATQC[0], ATQC[1], ATQC[2], ATQC[3], theTransmittance  );
#endif

        propagate_at_boundary( flag, rng, ctx, theTransmittance  );
    }
    return action ;
}
#endif

#endif


#if defined(__CUDACC__) || defined(__CUDABE__)  || defined(MOCK_CURAND) || defined(MOCK_CUDA)

/**
qsim::reflect_diffuse cf G4OpBoundaryProcess::DoReflection
-----------------------------------------------------------

* LobeReflection is not yet implemnented in qsim.h

::

    355 inline
    356 void G4OpBoundaryProcess_MOCK::DoReflection()
    357 {
    358         if ( theStatus == LambertianReflection ) {
    359
    360           NewMomentum = G4LambertianRand(theGlobalNormal);
    361           theFacetNormal = (NewMomentum - OldMomentum).unit();
    362
    363         }
    364         else if ( theFinish == ground ) {
    365
    366           theStatus = LobeReflection;
    367           if ( PropertyPointer1 && PropertyPointer2 ){
    368           } else {
    369              theFacetNormal =
    370                  GetFacetNormal(OldMomentum,theGlobalNormal);
    371           }
    372           G4double PdotN = OldMomentum * theFacetNormal;
    373           NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
    374
    375         }
    376         else {
    377
    378           theStatus = SpikeReflection;
    379           theFacetNormal = theGlobalNormal;
    380           G4double PdotN = OldMomentum * theFacetNormal;
    381           NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
    382
    383         }
    384         G4double EdotN = OldPolarization * theFacetNormal;
    385         NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
    386 }




Orient the normal depending on incident momentum to handle inside/outside
--------------------------------------------------------------------------

HMM: Geant4 is very flippy with normals whereas Opticks uses fixed
geometrical normals that may then need to be oriented for a specfic task,
such as diffuse reflection from the shape.


::

        geometrical normal
         :
       \ : /  diffuse reflection from outside in same hemi as geometrical normal : orient 1.f
        \:/   incident momentum is against the geometrical normal   dot( old_mom, *normal ) < 0.f
       --*------+
      /         |
     /          |
     \          /
      \ \   /  /   diffuse reflection from inside needs to use a flipped geometrical normal : orient -1.f
       \ \ /  /    incident moment is with the geometrical normal   dot( old_mom, *normal ) > 0.f
        \_*__/
          :
          :
          geometrical normal


orient is used to flip the reflection normal back against the incident direction

**/

inline QSIM_METHOD void qsim::reflect_diffuse( RNG& rng, sctx& ctx )
{
    sphoton& p = ctx.p ;

    float3 old_mom = p.mom ;

    const float3* normal = ctx.prd->normal() ;    // geometrical outwards normal

    //const float orient = -1.f ; // BUG : FIXED ORIENT FLIP CANNOT BE CORRECT
    const float orient = dot( old_mom, *normal ) > 0.f ? -1.f : 1.f ;

    lambertian_direction( &p.mom, normal, orient, rng, ctx );


    float3 facet_normal = normalize( p.mom - old_mom );
    const float EdotN = dot( p.pol, facet_normal );
    p.pol = -1.f*(p.pol) + 2.f*EdotN*facet_normal ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.reflect_diffuse pidx %7lld : old_mom = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  old_mom.x, old_mom.y, old_mom.z ) ;

    printf("//qsim.reflect_diffuse pidx %7lld : normal0 = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  normal->x, normal->y, normal->z ) ;

    printf("//qsim.reflect_diffuse pidx %7lld : p.mom = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  p.mom.x, p.mom.y, p.mom.z ) ;

    printf("//qsim.reflect_diffuse pidx %7lld : facet_normal = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  facet_normal.x, facet_normal.y, facet_normal.z ) ;
    }
#endif

}

/**
qsim::reflect_specular
-----------------------

TODO: check the fixed orient flip with specular reflections
from inside and outside

Its not so obvious here because the oriented normal is part of the
mom calc so it would be easy to double flip.

G4OpBoundaryProcess::PostStepDoIt does an early flip of theGlobalNormal
but the G4Navigator does flips before that too... so Geant4 is too flippy
to be helpful.

**/

inline QSIM_METHOD void qsim::reflect_specular( RNG& rng, sctx& ctx )
{
    sphoton& p = ctx.p ;
    const float3* normal = ctx.prd->normal() ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.reflect_specular.head pidx %7lld : normal0 = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  normal->x, normal->y, normal->z ) ;

    printf("//qsim.reflect_specular.head pidx %7lld : mom0 = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  p.mom.x, p.mom.y, p.mom.z ) ;

    printf("//qsim.reflect_specular.head pidx %7lld : pol0 = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  p.pol.x, p.pol.y, p.pol.z ) ;
    }
#endif

#ifdef WITH_ORIENT
    //const float orient = -1.f ;
    const float orient = 1.f ;
    // because orient appears twice in the below p.mom p.pol calcs
    // it being +1.f or -1.f makes no difference

    const float PdotN = dot( p.mom, *normal )*orient ;
    p.mom = p.mom - 2.f*PdotN*(*normal)*orient ;

    const float EdotN = dot( p.pol, *normal )*orient ;
    p.pol = -1.f*(p.pol) + 2.f*EdotN*(*normal)*orient  ;
#else
    // removed orient as does not effect calc, hence confusing and pointless
    const float PdotN = dot( p.mom, *normal ) ;
    p.mom = p.mom - 2.f*PdotN*(*normal) ;

    const float EdotN = dot( p.pol, *normal ) ;
    p.pol = -1.f*(p.pol) + 2.f*EdotN*(*normal)  ;
#endif

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if(ctx.pidx == base->pidx)
    {
    printf("//qsim.reflect_specular.tail pidx %7lld : mom1 = np.array([%10.5f,%10.5f,%10.5f]) ; PdotN = %10.5f ; EdotN = %10.5f \n",
        ctx.pidx,  p.mom.x, p.mom.y, p.mom.z, PdotN, EdotN  ) ;

    printf("//qsim.reflect_specular.tail pidx %7lld : pol1 = np.array([%10.5f,%10.5f,%10.5f]) \n",
        ctx.pidx,  p.pol.x, p.pol.y, p.pol.z ) ;

    }
#endif


}

/**
qsim::fake_propagate
-----------------------

This uses mock input prd (quad2) to provide a CUDA only (no OptiX, no geometry)
test of qsim propagation machinery.

* NB: qsim::mock_propagate is intended to be very similar to CSGOptiX/CSGOptiX7.cu::simulate

TODO
~~~~~

* compare with cx/CSGOptiX7.cu::simulate and find users of this to see if it could be made more similar to cx::simulate
* can sstate be slimmed : seems not very easily
* simplify sstate persisting (quad6?quad5?)

Stages within bounce loop
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

1. mock call to OptiX trace : doing "geometry lookup"
   photon position + direction -> surface normal + distance and identity, boundary

   * cosTheta sign gives boundary orientation

2. lookup material/surface properties using boundary and orientation (cosTheta) from geometry lookup

3. mutate photon and set flag using material properties

  * note that photons that SAIL to boundary are mutated twice within the while loop (by propagate_to_boundary and propagate_at_boundary/surface)

Thoughts
~~~~~~~~~~~

NB: the different collection streams (photon, record, rec, seq) must stay independent
so can switch them off easily in production running


**/

inline QSIM_METHOD void qsim::fake_propagate( sphoton& p, const quad2* mock_prd, RNG& rng, unsigned long long idx )
{
    p.set_flag(TORCH);  // setting initial flag : in reality this should be done by generation

    qsim* sim = this ;

    sctx ctx = {} ;
    ctx.p = p ;     // Q: Why is this different from CSGOptiX7.cu:simulate ? A: Presumably due to input photon.
    ctx.evt = evt ;
    ctx.pidx = idx ;

    int command = START ;
    int bounce = 0 ;
#ifndef PRODUCTION
    ctx.point(bounce);
#endif

    while( bounce < evt->max_bounce )
    {
        ctx.prd = mock_prd + (evt->max_bounce*idx+bounce) ;
        if( ctx.prd->boundary() == 0xffffu ) break ;   // SHOULD NEVER HAPPEN : propagate can do nothing meaningful without a boundary
#ifndef PRODUCTION
        ctx.trace(bounce);
#endif

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
        if(idx == base->pidx)
        printf("//qsim.fake_propagate pidx %7lld bounce %d evt.max_bounce %d prd.q0.f.xyzw (%10.4f %10.4f %10.4f %10.4f) \n",
             idx, bounce, evt->max_bounce, ctx.prd->q0.f.x, ctx.prd->q0.f.y, ctx.prd->q0.f.z, ctx.prd->q0.f.w );
#endif
        command = sim->propagate(bounce, rng, ctx );
        bounce++;
#ifndef PRODUCTION
        ctx.point(bounce);
#endif
        if(command == BREAK) break ;
    }
#ifndef PRODUCTION
    ctx.end();
#endif
    evt->photon[idx] = ctx.p ;
}

/**
qsim::propagate : one "bounce" propagate_to_boundary/propagate_at_boundary
-----------------------------------------------------------------------------

This is canonically invoked from within the bounce loop of CSGOptiX/OptiX7Test.cu:simulate
after tracing (or mock tracing) has populated ctx.prd

MISSING intersects
~~~~~~~~~~~~~~~~~~~

Formerly thought "MISSING needs to return BREAK" but MISSING
intersects only happen with "broken" test geometry so no need to be
concerned with that.

z-plus special sensor surfaces
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The two kinds of z-plus sensor surfaces need to fallback to ordinary
surface handling for z-minus, when lposcost < 0.f at lower z-hemi.
Initially thought would need to use optical.z to refer to
the fallback and call qbnd::fill_state again.
But of course that is wrong, the osur/isur surface references
already provide the z-minus fallback.
Hence can implement the fallback using standard surface branch
with no change to qbnd::fill_state.
For z-plus with lposcost > 0.f upper z-hemi simply ignore
the s.surface

traditional POM special surface : smatsur_Surface_zplus_sensor_A
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

HMM: could use standard surface handling for this but with perfect
detector surface swapped in.


TMM multilayer POM special surface : smatsur_Surface_zplus_sensor_CustomART
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Enabling special surface handling is controlled by ctx.optical.y "ems" enum
which wheels in the qpmt.h stack calc following C4CustomART::doIt
for CustomART special surfaces.

Prior to supporting special surfaces, within the command == BOUNDARY used::

    command = ctx.s.optical.x == 0 ?
                                      propagate_at_boundary( flag, rng, ctx )
                                  :
                                      propagate_at_surface( flag, rng, ctx )
                                  ;

**/

inline QSIM_METHOD int qsim::propagate(const int bounce, RNG& rng, sctx& ctx )  // ::simulate
{
    const unsigned boundary = ctx.prd->boundary() ;
    const unsigned identity = ctx.prd->identity() ; // sensor_identifier+1, 0:not-a-sensor
    const unsigned iindex = ctx.prd->iindex() ;
    const float lposcost = ctx.prd->lposcost() ;  // local frame intersect position cosine theta

    const float3* normal = ctx.prd->normal();
    float cosTheta = dot(ctx.p.mom, *normal ) ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if( ctx.pidx == base->pidx )
    {
    printf("\n//qsim.propagate.head pidx %7lld : ctx.evt.index %d evt.index %d \n", ctx.pidx, ctx.evt->index, evt->index );

    printf("\n//qsim.propagate.head pidx %7lld : bnc %d boundary %d cosTheta %10.8f \n", ctx.pidx, bounce, boundary, cosTheta );

    printf("//qsim.propagate.head pidx %7lld : mom = np.array([%10.8f,%10.8f,%10.8f]) ; lmom = %10.8f  \n",
                 ctx.pidx, ctx.p.mom.x, ctx.p.mom.y, ctx.p.mom.z, length(ctx.p.mom) ) ;

    printf("//qsim.propagate.head pidx %7lld : pos = np.array([%10.5f,%10.5f,%10.5f]) ; lpos = %10.8f \n",
                 ctx.pidx, ctx.p.pos.x, ctx.p.pos.y, ctx.p.pos.z, length(ctx.p.pos) ) ;

    printf("//qsim.propagate.head pidx %7lld : nrm = np.array([(%10.8f,%10.8f,%10.8f]) ; lnrm = %10.8f  \n",
                 ctx.pidx, normal->x, normal->y, normal->z, length(*normal) );

    }
#endif

    // 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 );

    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 );
    /**
    command possibilities:

    1. CONTINUE (eg scatter)
    2. BREAK (eg absorb)
    3. BOUNDARY

    **/

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if( ctx.pidx == base->pidx )
    printf("//qsim.propagate.body pidx %7lld bounce %d command %d flag %d s.optical.x %d s.optical.y %d \n",
          ctx.pidx, bounce, command, flag, ctx.s.optical.x, ctx.s.optical.y );
#endif

    if( command == BOUNDARY )
    {
        const int& ems = ctx.s.optical.y ;

#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
        if( ctx.pidx == base->pidx )
        {
#if defined(WITH_CUSTOM4)
            printf("//qsim.propagate.body.WITH_CUSTOM4 pidx %7lld  BOUNDARY ems %d lposcost %7.3f \n", ctx.pidx, ems, lposcost );
#else
            printf("//qsim.propagate.body.NOT:WITH_CUSTOM4 pidx %7lld BOUNDARY ems %d lposcost %7.3f \n", ctx.pidx, ems, lposcost);
#endif
        }
#endif

        if( ems == smatsur_NoSurface )
        {
            command = propagate_at_boundary( flag, rng, ctx ) ;
        }
        else if( ems == smatsur_Surface )
        {
            command = propagate_at_surface( flag, rng, ctx ) ;
        }
        else if( lposcost < 0.f )  // could combine with prior, but handy for debug to keep separate
        {
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
            if( ctx.pidx == base->pidx )
                printf("//qsim.propagate.body (lposcost < 0.f) pidx %7lld bounce %d command %d flag %d ems %d \n",
                ctx.pidx, bounce, command, flag, ems  );
#endif
            command = propagate_at_surface( flag, rng, ctx ) ;
        }
        else if( ems == smatsur_Surface_zplus_sensor_A )
        {
            command = propagate_at_surface_Detect( flag, rng, ctx ) ;
        }
        else if( ems == smatsur_Surface_zplus_sensor_CustomART )
        {
#if defined(WITH_CUSTOM4)
            command = propagate_at_surface_CustomART( flag, rng, ctx ) ;
            //command = base->custom_lut == 0u ? propagate_at_surface_CustomART( flag, rng, ctx ) : propagate_at_surface_MultiFilm(flag, rng, ctx );

#endif
        }
    }
    ctx.p.set_flag(flag);
    // Q: Does flag need to be single bit at this point OR can multiple "flags" be OR-ed together here ?
    // A: Decided to keep the flag as single bitted, and directly set EFFICENCY_COLLECT/CULL into ctx.p.flagmask


#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
    if( ctx.pidx == base->pidx )
    printf("//qsim.propagate.tail pidx %7lld bounce %d command %d flag %d ctx.s.optical.y(ems) %d \n",
             ctx.pidx, bounce, command, flag, ctx.s.optical.y  );
#endif

    return command ;
}
/**
Q: Where does ctx.s.optical come from ?
A: Populated in qbnd::fill_state based on boundary and cosTheta sign to get the
   su_line index of the optical buffer which distinguishes surface from boundary,
   as acted upon above. See also GBndLib::createOpticalBuffer

Q: Can the big "if" above be changed into a switch ?
A: YES, but non-trivially and probably confusingly. This is because
   the special surface handling condition needs to be after the possibiliy
   of fallback to ordinary surface handling for "lposcost < 0.f" and also
   ordinary NoSurface must be possible no matter the lposcost value.

**/



/**
qsim::hemisphere_polarized
------------------------------


          direction      | surface_normal
               \         |
                \        |
              .  \       |
          .       \      |
      .            \     |
     within         \    |
                     \   |
                      \  |
                       \ |
                        \|
           --------------+------------


*plane of incidence*
    plane containing *surface_normal* *direction* and *within* vectors

*transverse*
    vector transverse to the plane of incidence (S polarized)

*within*
    vector within the plane of incidence and perpendicular to *direction* (P polarized)


A +ve Z upper hemisphere of *direction* is generated and then rotateUz oriented
to adopt the *surface_normal* vector as its Z direction.

For inwards=true the normal direction is flipped to orient all the directions
inwards.

**/

inline QSIM_METHOD void qsim::hemisphere_polarized( unsigned polz, bool inwards, RNG& rng, sctx& ctx )
{
    sphoton& p = ctx.p ;
    const float3* normal = ctx.prd->normal() ;

    //printf("//qsim.hemisphere_polarized polz %d normal (%10.4f, %10.4f, %10.4f) \n", polz, normal->x, normal->y, normal->z );

    float u_hemipol_phi = curand_uniform(&rng) ;
    float phi = u_hemipol_phi*2.f*M_PIf;  // 0->2pi
    float cosTheta = curand_uniform(&rng) ;      // 0->1

#if !defined(PRODUCTION) && defined(DEBUG_TAG)
    stagr& tagr = ctx.tagr ;
    tagr.add( stag_hp_ph, u_hemipol_phi );
    tagr.add( stag_hp_ph, cosTheta );    // trying to reduce stag::BITS from 5 to 4, so change stag_hp_ct to stag_hp_ph
#endif

    float sinTheta = sqrtf(1.f-cosTheta*cosTheta);

    p.mom.x = cosf(phi)*sinTheta ;
    p.mom.y = sinf(phi)*sinTheta ;
    p.mom.z = cosTheta ;

    smath::rotateUz( p.mom, (*normal) * ( inwards ? -1.f : 1.f ));

    //printf("//qsim.hemisphere_polarized polz %d p.mom (%10.4f, %10.4f, %10.4f) \n", polz, p.mom.x, p.mom.y, p.mom.z );

    // what about normal incidence ?
    const float3 transverse = normalize(cross(p.mom, (*normal) * ( inwards ? -1.f : 1.f )  )) ; // perpendicular to plane of incidence
    const float3 within = normalize( cross(p.mom, transverse) );  //   within plane of incidence and perpendicular to direction

    switch(polz)
    {
        case 0: p.pol = transverse ; break ;   // S-polarizatiom
        case 1: p.pol = within     ; break ;   // P-polarization
        case 2: p.pol = normalize( 0.5f*transverse + (1.f-0.5f)*within )  ; break ;  // equal admixture
    }
}




/**
qsim::generate_photon_simtrace
--------------------------------

Canonical invokation from CSGOptiX7.cu:simtrace

* NB simtrace cxs center-extent-genstep are very different to standard Cerenkov/Scintillation gensteps

These gensteps are for example created in SEvent::MakeCenterExtentGensteps the below generation
should be comparable to the CPU implementation::

    SFrameGenstep::GenerateCenterExtentGenstepPhotons
    SFrameGenstep::GenerateSimtracePhotons

HMM: note that the photon_id is global to the launch making it potentially a very large number
but for identication purposes having a local to the photons of the genstep index would
be more useful and would allow storage within much less bits.

TODO: implement local index by including photon_id offset with the gensteps

* NB the sxyz.h enum order is different to the python one  eg XYZ=0

SEE : sevent::add_simtrace

**/

inline QSIM_METHOD void qsim::generate_photon_simtrace(quad4& p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const
{
    const int& gencode = gs.q0.i.x ;
    switch(gencode)
    {
        case OpticksGenstep_FRAME:                   generate_photon_simtrace_frame(p, rng, gs, photon_id, genstep_id ); break ;
        case OpticksGenstep_INPUT_PHOTON_SIMTRACE:   { p = (quad4&)evt->simtrace[photon_id] ; }                          ; break ;
    }
}

inline QSIM_METHOD void qsim::generate_photon_simtrace_frame(quad4& p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const
{
    C4U gsid ;

    //int gencode          = gs.q0.i.x ;
    int gridaxes           = gs.q0.i.y ;  // { XYZ, YZ, XZ, XY }
    gsid.u                 = gs.q0.i.z ;
    //unsigned num_photons = gs.q0.u.w ;

    p.q0.f.x = gs.q1.f.x ;   // start with genstep local frame position, typically origin  (0,0,0)
    p.q0.f.y = gs.q1.f.y ;
    p.q0.f.z = gs.q1.f.z ;
    p.q0.f.w = 1.f ;

    //printf("//qsim.generate_photon_simtrace_frame gridaxes %d gs.q1 (%10.4f %10.4f %10.4f %10.4f) \n", gridaxes, gs.q1.f.x, gs.q1.f.y, gs.q1.f.z, gs.q1.f.w );

    float u0 = curand_uniform(&rng);
    float sinPhi, cosPhi;
#if defined(MOCK_CURAND) || defined(MOCK_CUDA)
    __sincosf(2.f*M_PIf*u0,&sinPhi,&cosPhi);
#else
    sincosf(2.f*M_PIf*u0,&sinPhi,&cosPhi);
#endif

    float u1 = curand_uniform(&rng);
    float cosTheta = 2.f*u1 - 1.f ;
    float sinTheta = sqrtf(1.f-cosTheta*cosTheta) ;

    //printf("//qsim.generate_photon_simtrace_frame u0 %10.4f sinPhi   %10.4f cosPhi   %10.4f \n", u0, sinPhi, cosPhi );
    //printf("//qsim.generate_photon_simtrace_frame u1 %10.4f sinTheta %10.4f cosTheta %10.4f \n", u1, sinTheta, cosTheta );
    //printf("//qsim.generate_photon_simtrace_frame  u0 %10.4f sinPhi   %10.4f cosPhi   %10.4f u1 %10.4f sinTheta %10.4f cosTheta %10.4f \n",  u0, sinPhi, cosPhi, u1, sinTheta, cosTheta );

    switch( gridaxes )
    {
        case YZ:  { p.q1.f.x = 0.f    ;  p.q1.f.y = cosPhi ;  p.q1.f.z = sinPhi ;  p.q1.f.w = 0.f ; } ; break ;
        case XZ:  { p.q1.f.x = cosPhi ;  p.q1.f.y = 0.f    ;  p.q1.f.z = sinPhi ;  p.q1.f.w = 0.f ; } ; break ;
        case XY:  { p.q1.f.x = cosPhi ;  p.q1.f.y = sinPhi ;  p.q1.f.z = 0.f    ;  p.q1.f.w = 0.f ; } ; break ;
        case XYZ: { p.q1.f.x = sinTheta*cosPhi ;
                    p.q1.f.y = sinTheta*sinPhi ;
                    p.q1.f.z = cosTheta        ;
                    p.q1.f.w = 0.f ; } ; break ;   // previously used XZ
    }


    qat4 qt(gs) ; // copy 4x4 transform from last 4 quads of genstep
    qt.right_multiply_inplace( p.q0.f, 1.f );   // position
    qt.right_multiply_inplace( p.q1.f, 0.f );   // direction


    unsigned char ucj = (photon_id < 255 ? photon_id : 255 ) ;
    gsid.c4.w = ucj ;
    p.q3.u.w = gsid.u ;   // WARNING : THIS GSID LOOKS TO BE STOMPED ON BY sevent::add_simtrace
}

/**
qsim::generate_photon
----------------------

Moved non-standard center-extent (aka frame) gensteps to use qsim::generate_photon_simtrace not this

**/

inline QSIM_METHOD void qsim::generate_photon(sphoton& p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const
{
    const int& gencode = gs.q0.i.x ;
    switch(gencode)
    {
        case OpticksGenstep_CARRIER:         scarrier::generate(     p, rng, gs, photon_id, genstep_id)  ; break ;
        case OpticksGenstep_TORCH:           storch::generate(       p, rng, gs, photon_id, genstep_id ) ; break ;

        case OpticksGenstep_G4Cerenkov_modified:
        case OpticksGenstep_CERENKOV:
                                              cerenkov->generate(    p, rng, gs, photon_id, genstep_id ) ; break ;

        case OpticksGenstep_DsG4Scintillation_r4695:
        case OpticksGenstep_SCINTILLATION:
                                              scint->generate(        p, rng, gs, photon_id, genstep_id ) ; break ;

        case OpticksGenstep_INPUT_PHOTON:    { p = evt->photon[photon_id] ; p.set_flag(TORCH) ; }        ; break ;
        default:                             generate_photon_dummy(  p, rng, gs, photon_id, genstep_id)  ; break ;
    }
    p.set_index(photon_id) ;
}
#endif