eic-opticks/qudarap/qcerenkov.h

0001 #pragma once
0002 /**
0003 qcerenkov.h
0004 ==============
0005
0006 **/
0007
0008 #if defined(__CUDACC__) || defined(__CUDABE__)
0009    #define QCERENKOV_METHOD __device__
0010 #else
0011    #define QCERENKOV_METHOD
0012 #endif
0013
0014 #include "scurand.h"
0015 #include "smath.h"
0016 #include "OpticksPhoton.h"
0017
0018 #include "qrng.h"
0019
0020
0021 struct scerenkov ;
0022 struct qbase ;
0023 struct qbnd ;
0024 struct quad6 ;
0025 template <typename T> struct qprop ;
0026
0027 /**
0028
0029 HMM: the qsim member causes chicken-egg setup problem,
0030 as want a cerenkov member of qsim
0031
0032 * using it for sim->boundary_lookup, sim->prop
0033   so need to break things up at finer level for easier reuse
0034
0035 **/
0036
0037 struct qcerenkov
0038 {
0039     qbase* base ;
0040     qbnd*  bnd ;
0041     qprop<float>*  prop ;
0042
0043 #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
0044     QCERENKOV_METHOD void generate( sphoton& p,  RNG& rng, const quad6& gs, unsigned long long idx, int genstep_id ) const ;
0045
0046     template<typename T>
0047     QCERENKOV_METHOD void wavelength_sampled_enprop( float& wavelength, float& cosTheta, float& sin2Theta, RNG& rng, const scerenkov& gs, unsigned long long idx, int genstep_id ) const ;
0048     QCERENKOV_METHOD void wavelength_sampled_bndtex( float& wavelength, float& cosTheta, float& sin2Theta, RNG& rng, const scerenkov& gs, unsigned long long idx, int genstep_id ) const ;
0049
0050     QCERENKOV_METHOD void fraction_sampled(float& fraction, float& delta, RNG& rng, const scerenkov& gs, unsigned long long idx, int gsid ) const ;
0051 #endif
0052
0053 };
0054
0055 #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
0056 inline QCERENKOV_METHOD void qcerenkov::generate( sphoton& p, RNG& rng, const quad6& _gs, unsigned long long idx, int gsid ) const
0057 {
0058     //printf("//qcerenkov::generate idx %d \n", idx );
0059
0060     const scerenkov& gs = (const scerenkov&)_gs ;
0061     const float3 p0 = normalize(gs.DeltaPosition) ;   // TODO: see if can normalize inside the genstep at collection
0062
0063     float wavelength ;
0064     float cosTheta ;
0065     float sin2Theta ;
0066
0067     //wavelength = 500.f ; cosTheta = 0.70710678f ; sin2Theta = 0.5f ;
0068
0069     wavelength_sampled_bndtex(wavelength, cosTheta, sin2Theta, rng, gs, idx, gsid) ;
0070     //wavelength_sampled_enprop<float>(wavelength, cosTheta, sin2Theta, rng, gs, idx, gsid) ;
0071     //wavelength_sampled_enprop<double>(wavelength, cosTheta, sin2Theta, rng, gs, idx, gsid) ;
0072
0073     float sinTheta = sqrtf(sin2Theta);
0074
0075     // Generate random position of photon on cone surface
0076     // defined by Theta
0077
0078     float u0 = curand_uniform(&rng);
0079     float phi = 2.f*M_PIf*u0 ;
0080     float sinPhi = sin(phi);
0081     float cosPhi = cos(phi);
0082
0083     // calculate x,y, and z components of photon energy
0084     // (in coord system with primary particle direction
0085     //  aligned with the z axis)
0086
0087     p.mom.x = sinTheta*cosPhi ;
0088     p.mom.y = sinTheta*sinPhi ;
0089     p.mom.z = cosTheta ;
0090
0091     // Rotate momentum direction back to global reference system
0092     smath::rotateUz(p.mom, p0 );
0093
0094     // Determine polarization of new photon
0095
0096     p.pol.x = cosTheta*cosPhi ;
0097     p.pol.y = cosTheta*sinPhi ;
0098     p.pol.z = -sinTheta ;
0099
0100
0101     // Rotate back to original coord system
0102     smath::rotateUz(p.pol, p0 );
0103
0104     p.wavelength = wavelength ;
0105
0106     float fraction ;
0107     float delta ;
0108
0109     fraction_sampled( fraction, delta, rng, gs, idx, gsid );
0110
0111     float midVelocity = gs.preVelocity + fraction*( gs.postVelocity - gs.preVelocity )*0.5f ;
0112
0113     p.time = gs.time + delta / midVelocity ;
0114     p.pos = gs.pos + fraction * gs.DeltaPosition ;   // NB here gs.DeltaPosition must not be normalized
0115
0116     p.zero_flags();
0117     p.set_flag(CERENKOV) ;
0118
0119 }
0120
0121 /**
0122 qcerenkov::fraction_sampled
0123 ------------------------------
0124
0125 Note that N and NumberOfPhotons are never used below.
0126 The point of the the rejection sampling loop is to come up with a
0127 *fraction* and *delta* that fulfils the theoretical constraint.
0128 This *fraction* controls where the photon gets generated along the
0129 genstep line segment with the rejection sampling serving
0130 to get the appropriate distribution of generation along that line.
0131
0132 Consider the case of DeltaN zero due to gs.MeanNumberOfPhotons1
0133 and gs.MeanNumberOfPhotons2 being equal. The NumberOfPhotons
0134 gets set to MeanNumberOfPhotons1. N will be less that that
0135 so the loop will only turn once with fraction being random.
0136
0137 **/
0138
0139 inline QCERENKOV_METHOD void qcerenkov::fraction_sampled(float& fraction, float& delta, RNG& rng, const scerenkov& gs, unsigned long long idx, int gsid ) const
0140 {
0141     float NumberOfPhotons ;
0142     float N ;
0143     float u ;
0144
0145     float MeanNumberOfPhotonsMax = fmaxf( gs.MeanNumberOfPhotons1, gs.MeanNumberOfPhotons2 );
0146     float DeltaN = (gs.MeanNumberOfPhotons1-gs.MeanNumberOfPhotons2) ;
0147     do
0148     {
0149         fraction = curand_uniform(&rng) ;
0150
0151         delta = fraction * gs.step_length ;
0152
0153         NumberOfPhotons = gs.MeanNumberOfPhotons1 - fraction * DeltaN  ;
0154
0155         u = curand_uniform(&rng) ;
0156
0157         N = u * MeanNumberOfPhotonsMax ;
0158
0159     } while (N > NumberOfPhotons);
0160
0161     //printf("//qcerenkov::fraction_sampled fraction %10.4f delta %10.4f \n", fraction, delta );
0162
0163 }
0164
0165
0166
0167 /**
0168 qcerenkov::wavelength_sampled_bndtex
0169 -----------------------------------------------------
0170
0171 wavelength between Wmin and Wmax is uniform-reciprocal-sampled
0172 to mimic uniform energy range sampling without taking reciprocals
0173 twice
0174
0175 g4-cls G4Cerenkov::
0176
0177     251   G4double Pmin = Rindex->GetMinLowEdgeEnergy();
0178     252   G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
0179     253   G4double dp = Pmax - Pmin;
0180     254
0181     255   G4double nMax = Rindex->GetMaxValue();
0182     256
0183     257   G4double BetaInverse = 1./beta;
0184     258
0185     259   G4double maxCos = BetaInverse / nMax;
0186     260   G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
0187     261
0188     ...
0189     270   for (G4int i = 0; i < fNumPhotons; i++) {
0190     271
0191     272       // Determine photon energy
0192     273
0193     274       G4double rand;
0194     275       G4double sampledEnergy, sampledRI;
0195     276       G4double cosTheta, sin2Theta;
0196     277
0197     278       // sample an energy
0198     279
0199     280       do {
0200     281          rand = G4UniformRand();
0201     282          sampledEnergy = Pmin + rand * dp;                // linear energy sample in Pmin -> Pmax
0202     283          sampledRI = Rindex->Value(sampledEnergy);
0203     284          cosTheta = BetaInverse / sampledRI;              // what cone angle that energy sample corresponds to
0204     285
0205     286          sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
0206     287          rand = G4UniformRand();
0207     288
0208     289         // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
0209     290       } while (rand*maxSin2 > sin2Theta);                // constrain
0210     291
0211
0212 ::
0213
0214
0215
0216                   \    B    /
0217               \    .   |   .    /                            AC     ct / n          1         i       BetaInverse
0218           \    C       |       C    /             cos th =  ---- =  --------   =  ------ =   ---  =  -------------
0219       \    .    \      |      /    .     /                   AB       bct           b n       n        sampledRI
0220        .         \    bct    /          .
0221                   \    |    /                                  BetaInverse
0222                    \   |   /  ct                  maxCos  =  -----------------
0223                     \  |th/  ----                                nMax
0224                      \ | /    n
0225                       \|/
0226                        A
0227
0228     Particle travels AB, light travels AC,  ACB is right angle
0229
0230
0231      Only get Cerenkov radiation when
0232
0233             cos th <= 1 ,
0234
0235             beta >= beta_min = 1/n        BetaInverse <= BetaInverse_max = n
0236
0237
0238      At the small beta threshold AB = AC,   beta = beta_min = 1/n     eg for n = 1.333, beta_min = 0.75
0239
0240             cos th = 1,  th = 0         light is in direction of the particle
0241
0242
0243      For ultra relativistic particle beta = 1, there is a maximum angle
0244
0245             th = arccos( 1/n )
0246
0247     In [5]: np.arccos(0.75)*180./np.pi
0248     Out[5]: 41.40962210927086
0249
0250
0251      So the beta range to have Cerenkov is  :
0252
0253                 1/n       slowest, cos th = 1, th = 0
0254
0255           ->    1         ultra-relativistic, maximum cone angle th  arccos(1/n)
0256
0257
0258      Corresponds to BetaInverse range
0259
0260            BetaInverse =  n            slowest, cos th = 1, th = 0    cone in particle direction
0261
0262            BetaInverse  = 1
0263
0264
0265
0266
0267      The above considers a fixed refractive index.
0268      Actually refractive index varies with wavelength resulting in
0269      a range of cone angles for a fixed particle beta.
0270
0271
0272     * https://www2.physics.ox.ac.uk/sites/default/files/2013-08-20/external_pergamon_jelley_pdf_18410.pdf
0273     * ~/opticks_refs/external_pergamon_jelley_pdf_18410.pdf
0274
0275
0276 HMM: when there is no permissable cerenkov for the material RINDEX this is just going to infinite spin ?
0277 Added the loop count to prevent that.  BUT whats the proper way to avoid that ?
0278
0279 The numphotons should be zero in that situation : so this is an issue of testing with
0280 fabricated gensteps that should not be an issue with real gensteps.
0281
0282
0283 **/
0284
0285 inline QCERENKOV_METHOD void qcerenkov::wavelength_sampled_bndtex(float& wavelength, float& cosTheta, float& sin2Theta, RNG& rng, const scerenkov& gs, unsigned long long idx, int gsid ) const
0286 {
0287     //printf("//qcerenkov::wavelength_sampled_bndtex bnd %p gs.matline %d \n", bnd, gs.matline );
0288     float u0 ;
0289     float u1 ;
0290     float w ;
0291     float sampledRI ;
0292     float u_maxSin2 ;
0293
0294     unsigned count = 0 ;
0295
0296     do {
0297         u0 = curand_uniform(&rng) ;
0298
0299         w = gs.Wmin + u0*(gs.Wmax - gs.Wmin) ;
0300
0301         wavelength = gs.Wmin*gs.Wmax/w ; // reciprocalization : arranges flat energy distribution, expressed in wavelength
0302
0303         float4 props = bnd->boundary_lookup(wavelength, gs.matline, 0u);
0304
0305         sampledRI = props.x ;
0306
0307         //printf("//qcerenkov::wavelength_sampled_bndtex count %d wavelength %10.4f sampledRI %10.4f \n", count, wavelength, sampledRI );
0308
0309         cosTheta = gs.BetaInverse / sampledRI ;
0310
0311         sin2Theta = fmaxf( 0.f, (1.f - cosTheta)*(1.f + cosTheta));
0312
0313         u1 = curand_uniform(&rng) ;
0314
0315         u_maxSin2 = u1*gs.maxSin2 ;
0316
0317         count += 1 ;
0318
0319     } while ( u_maxSin2 > sin2Theta && count < 100 );
0320
0321 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0322     if(count > 50)
0323     printf("//qcerenkov::wavelength_sampled_bndtex idx %6lld sampledRI %7.3f cosTheta %7.3f sin2Theta %7.3f wavelength %7.3f count %d matline %d \n",
0324               idx , sampledRI, cosTheta, sin2Theta, wavelength, count, gs.matline );
0325 #endif
0326
0327 }
0328
0329
0330 /**
0331 qcerenkov::wavelength_sampled_enprop
0332 --------------------------------------
0333
0334 template type controls the type used for the rejection sampling, not the return type
0335
0336
0337 **/
0338 template<typename T>
0339 inline QCERENKOV_METHOD void qcerenkov::wavelength_sampled_enprop(float& f_wavelength, float& f_cosTheta, float& f_sin2Theta, RNG& rng, const scerenkov& gs, unsigned long long idx, int gsid ) const
0340 {
0341     T u0 ;
0342     T u1 ;
0343     T energy ;
0344     T sampledRI ;
0345     T cosTheta ;
0346     T sin2Theta ;
0347     T u_mxs2_s2 ;
0348
0349     T one(1.) ;
0350     T zero(0.) ;
0351
0352     T pmin = gs.Pmin() ;
0353     T pmax = gs.Pmax() ;
0354
0355     unsigned loop = 0u ;
0356
0357     do {
0358
0359         u0 = scurand<T>::uniform(&rng) ;
0360
0361         energy = pmin + u0*(pmax - pmin) ;
0362
0363         sampledRI = prop->interpolate( 0u, energy );
0364
0365         cosTheta = gs.BetaInverse / sampledRI ;
0366
0367         sin2Theta = (one - cosTheta)*(one + cosTheta);
0368
0369         u1 = scurand<T>::uniform(&rng) ;
0370
0371         u_mxs2_s2 = u1*gs.maxSin2 - sin2Theta ;
0372
0373         loop += 1 ;
0374
0375         if( idx == base->pidx )
0376         {
0377             printf("//qcerenkov::cerenkov_generate_enprop idx %d loop %3d u0 %10.5f ri %10.5f ct %10.5f s2 %10.5f u_mxs2_s2 %10.5f \n", idx, loop, u0, sampledRI, cosTheta, sin2Theta, u_mxs2_s2 );
0378         }
0379
0380
0381     } while ( u_mxs2_s2 > zero );
0382
0383     f_wavelength = smath::hc_eVnm/energy ;
0384     f_cosTheta = cosTheta ;
0385     f_sin2Theta = sin2Theta ;
0386 }
0387
0388 #endif
0389
0390
0391