eic-opticks/qudarap/qcerenkov_dev.h

0001 #pragma once
0002
0003
0004 #if defined(__CUDACC__) || defined(__CUDABE__)
0005    #define QCERENKOV_DEV_METHOD __device__
0006 #else
0007    #define QCERENKOV_DEV_METHOD
0008 #endif
0009
0010
0011 #include "qrng.h"
0012 #include "scerenkov.h"
0013
0014
0015 struct qcerenkov_dev
0016 {
0017 #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
0018
0019     QCERENKOV_DEV_METHOD static void generate(              qsim* sim, quad4& q, unsigned id, RNG& rng, const scerenkov& gs );
0020
0021     template<typename T>
0022     QCERENKOV_DEV_METHOD static void generate_enprop(       qsim* sim, quad4& q, unsigned id, RNG& rng, const scerenkov& gs );
0023     QCERENKOV_DEV_METHOD static void generate_expt_double(  qsim* sim, quad4& q, unsigned id, RNG& rng );
0024
0025
0026     QCERENKOV_DEV_METHOD static float wavelength_rejection_sampled(qsim* sim, unsigned id, RNG& rng ) ;
0027     QCERENKOV_DEV_METHOD static void  generate(qsim* sim, quad4& p, unsigned id, RNG& rng );
0028
0029     template<typename T>
0030     QCERENKOV_DEV_METHOD static void generate_enprop(qsim* sim, quad4& p, unsigned id, RNG& rng);
0031
0032
0033 #endif
0034
0035 };
0036
0037
0038
0039
0040 inline QCERENKOV_DEV_METHOD void qcerenkov_dev::generate(qsim* sim, quad4& q, unsigned id, RNG& rng, const scerenkov& gs )
0041 {
0042     float u0 ;
0043     float u1 ;
0044
0045
0046     float w_linear ;
0047     float wavelength ;
0048
0049     float sampledRI ;
0050     float cosTheta ;
0051     float sin2Theta ;
0052     float u_mxs2_s2 ;
0053
0054     unsigned line = gs.matline ; //   line :  4*boundary_idx + OMAT/IMAT (0/3)
0055
0056     unsigned loop = 0u ;
0057
0058     do {
0059
0060 #ifdef FLIP_RANDOM
0061         u0 = 1.f - curand_uniform(&rng) ;
0062 #else
0063         u0 = curand_uniform(&rng) ;
0064 #endif
0065
0066         w_linear = gs.Wmin + u0*(gs.Wmax - gs.Wmin) ;
0067
0068         wavelength = gs.Wmin*gs.Wmax/w_linear ;
0069
0070         float4 props = sim->boundary_lookup( wavelength, line, 0u);
0071
0072         sampledRI = props.x ;
0073
0074         cosTheta = gs.BetaInverse / sampledRI ;
0075
0076         sin2Theta = (1.f - cosTheta)*(1.f + cosTheta);
0077
0078 #ifdef FLIP_RANDOM
0079         u1 = 1.f - curand_uniform(&rng) ;
0080 #else
0081         u1 = curand_uniform(&rng) ;
0082 #endif
0083
0084         u_mxs2_s2 = u1*gs.maxSin2 - sin2Theta ;
0085
0086         loop += 1 ;
0087
0088         if( id == sim->base->pidx )
0089         {
0090             printf("//qcerenkov_dev::cerenkov_generate id %d loop %3d u0 %10.5f ri %10.5f ct %10.5f s2 %10.5f u_mxs2_s2 %10.5f \n", id, loop, u0, sampledRI, cosTheta, sin2Theta, u_mxs2_s2 );
0091         }
0092
0093
0094     } while ( u_mxs2_s2 > 0.f );
0095
0096     float energy = smath::hc_eVnm/wavelength ;
0097
0098     q.q0.f.x = energy ;
0099     q.q0.f.y = wavelength ;
0100     q.q0.f.z = sampledRI ;
0101     q.q0.f.w = cosTheta ;
0102
0103     q.q1.f.x = sin2Theta ;
0104     q.q1.u.y = 0u ;
0105     q.q1.u.z = 0u ;
0106     q.q1.f.w = gs.BetaInverse ;
0107
0108     q.q2.f.x = w_linear ;    // linear sampled wavelenth
0109     q.q2.f.y = wavelength ;  // reciprocalized trick : does it really work
0110     q.q2.f.z = u0 ;
0111     q.q2.f.w = u1 ;
0112
0113     q.q3.u.x = line ;
0114     q.q3.u.y = loop ;
0115     q.q3.f.z = 0.f ;
0116     q.q3.f.w = 0.f ;
0117 }
0118
0119
0120
0121
0122
0123
0124 /**
0125 qcerenkov_dev::generate_enprop
0126 --------------------------------
0127
0128 Variation assuming Wmin, Wmax contain Pmin Pmax and using qprop::interpolate
0129 to sample the RINDEX
0130
0131 **/
0132
0133
0134
0135 template<typename T>
0136 inline QCERENKOV_DEV_METHOD void qcerenkov_dev::generate_enprop(qsim* sim, quad4& q, unsigned id, RNG& rng, const scerenkov& gs )
0137 {
0138     T u0 ;
0139     T u1 ;
0140     T energy ;
0141     T sampledRI ;
0142     T cosTheta ;
0143     T sin2Theta ;
0144     T u_mxs2_s2 ;
0145
0146     T one(1.) ;
0147
0148     unsigned line = gs.matline ; //   line :  4*boundary_idx + OMAT/IMAT (0/3)
0149     unsigned loop = 0u ;
0150
0151     do {
0152
0153         u0 = scurand<T>::uniform(&rng) ;
0154
0155         energy = gs.Wmin + u0*(gs.Wmax - gs.Wmin) ;
0156
0157         sampledRI = sim->prop->interpolate( 0u, energy );
0158
0159         cosTheta = gs.BetaInverse / sampledRI ;
0160
0161         sin2Theta = (one - cosTheta)*(one + cosTheta);
0162
0163         u1 = scurand<T>::uniform(&rng) ;
0164
0165         u_mxs2_s2 = u1*gs.maxSin2 - sin2Theta ;
0166
0167         loop += 1 ;
0168
0169         if( id == sim->base->pidx )
0170         {
0171             printf("//qcerenkov_dev::generate_enprop id %d loop %3d u0 %10.5f ri %10.5f ct %10.5f s2 %10.5f u_mxs2_s2 %10.5f \n", id, loop, u0, sampledRI, cosTheta, sin2Theta, u_mxs2_s2 );
0172         }
0173
0174
0175     } while ( u_mxs2_s2 > 0.f );
0176
0177
0178     float wavelength = smath::hc_eVnm/energy ;
0179
0180
0181
0182     q.q0.f.x = energy ;
0183     q.q0.f.y = wavelength ;
0184     q.q0.f.z = sampledRI ;
0185     q.q0.f.w = cosTheta ;
0186
0187     q.q1.f.x = sin2Theta ;
0188     q.q1.u.y = 0u ;
0189     q.q1.u.z = 0u ;
0190     q.q1.f.w = gs.BetaInverse ;
0191
0192     q.q2.f.x = 0.f ;
0193     q.q2.f.y = 0.f ;
0194     q.q2.f.z = u0 ;
0195     q.q2.f.w = u1 ;
0196
0197     q.q3.u.x = line ;
0198     q.q3.u.y = loop ;
0199     q.q3.f.z = 0.f ;
0200     q.q3.f.w = 0.f ;
0201 }
0202
0203
0204
0205
0206
0207
0208
0209 /**
0210 qcerenkov_dev::generate_expt_double
0211 -------------------------------------
0212
0213 This does the sampling all in double, narrowing to
0214 float just for the photon output.
0215
0216 Note that this is not using a genstep.
0217
0218 Which things have most need to be  double to make any difference ?
0219
0220 **/
0221
0222 inline QCERENKOV_DEV_METHOD void qcerenkov_dev::generate_expt_double(qsim* sim, quad4& q, unsigned id, RNG& rng )
0223 {
0224     double BetaInverse = 1.5 ;
0225     double Pmin = 1.55 ;
0226     double Pmax = 15.5 ;
0227     double nMax = 1.793 ;
0228
0229     //double maxOneMinusCosTheta = (nMax - BetaInverse) / nMax;
0230
0231     double maxCos = BetaInverse / nMax;
0232     double maxSin2 = ( 1. - maxCos )*( 1. + maxCos );
0233     double cosTheta ;
0234     double sin2Theta ;
0235
0236     double reject ;
0237     double u0 ;
0238     double u1 ;
0239     double energy ;
0240     double sampledRI ;
0241
0242     //double oneMinusCosTheta ;
0243
0244     unsigned loop = 0u ;
0245
0246     do {
0247         u0 = curand_uniform_double(&rng) ;
0248         u1 = curand_uniform_double(&rng) ;
0249         energy = Pmin + u0*(Pmax - Pmin) ;
0250         sampledRI = sim->prop->interpolate( 0u, energy );
0251         //oneMinusCosTheta = (sampledRI - BetaInverse) / sampledRI ;
0252         //reject = u1*maxOneMinusCosTheta - oneMinusCosTheta ;
0253         loop += 1 ;
0254
0255         cosTheta = BetaInverse / sampledRI ;
0256         sin2Theta = (1. - cosTheta)*(1. + cosTheta);
0257         reject = u1*maxSin2 - sin2Theta ;
0258
0259     } while ( reject > 0. );
0260
0261
0262     // narrowing for output
0263     q.q0.f.x = energy ;
0264     q.q0.f.y = smath::hc_eVnm/energy ;
0265     q.q0.f.z = sampledRI ;
0266     //p.q0.f.w = 1. - oneMinusCosTheta ;
0267     q.q0.f.w = cosTheta ;
0268
0269     q.q1.f.x = sin2Theta ;
0270     q.q1.u.y = 0u ;
0271     q.q1.u.z = 0u ;
0272     q.q1.f.w = BetaInverse ;
0273
0274     q.q2.f.x = reject ;
0275     q.q2.f.y = 0.f ;
0276     q.q2.f.z = u0 ;
0277     q.q2.f.w = u1 ;
0278
0279     q.q3.f.x = 0.f ;
0280     q.q3.u.y = loop ;
0281     q.q3.f.z = 0.f ;
0282     q.q3.f.w = 0.f ;
0283 }
0284
0285
0286
0287
0288 /**
0289 qcerenkov_dev:cerenkov_wavelength_rejection_sampled
0290 --------------------------------------------
0291
0292 HUH: this is using a GPU fabricated genstep everytime : that is kinda crazy approach.
0293 Makes much more sense to fabricate genstep on CPU and upload it.
0294
0295 **/
0296
0297
0298 inline QCERENKOV_DEV_METHOD float qcerenkov_dev::wavelength_rejection_sampled(qsim* sim, unsigned id, RNG& rng )
0299 {
0300
0301     int matline = 0u ;
0302     int numphoton_per_genstep = 10u ;
0303
0304     // MAKES MORE SENSE TO PREP GS ON CPU ?
0305     scerenkov gs ;
0306     scerenkov::FillGenstep(gs, matline, numphoton_per_genstep, false );
0307
0308     float wavelength = wavelength_rejection_sampled(sim, id, rng );
0309     return wavelength ;
0310 }
0311
0312 inline QCERENKOV_DEV_METHOD void qcerenkov_dev::generate(qsim* sim, quad4& p, unsigned id, RNG& rng )
0313 {
0314     int matline = 0u ;
0315     int numphoton_per_genstep = 10u ;
0316
0317     scerenkov gs ;
0318     scerenkov::FillGenstep(gs, matline, numphoton_per_genstep, false );
0319
0320     generate(sim, p, id, rng, gs );
0321 }
0322
0323 template<typename T>
0324 inline QCERENKOV_DEV_METHOD void qcerenkov_dev::generate_enprop(qsim* sim, quad4& p, unsigned id, RNG& rng)
0325 {
0326     int matline = 0u ;
0327     int numphoton_per_genstep = 10u ;
0328
0329     scerenkov gs ;
0330     scerenkov::FillGenstep(gs, matline, numphoton_per_genstep, false );
0331
0332     generate_enprop<T>(sim, p, id, rng, gs );
0333 }
0334