eic-opticks/qudarap/qscint.h

0001 #pragma once
0002 /**
0003 qscint.h
0004 ==================
0005
0006 **/
0007
0008 #if defined(__CUDACC__) || defined(__CUDABE__)
0009    #define QSCINT_METHOD __device__
0010 #else
0011    #define QSCINT_METHOD
0012 #endif
0013
0014
0015 #include "qrng.h"
0016
0017 struct quad4 ;
0018 struct quad6 ;
0019 struct sphoton ;
0020
0021 #include "OpticksPhoton.h"
0022
0023 struct qscint
0024 {
0025     cudaTextureObject_t scint_tex ;
0026     quad4*              scint_meta ; // HUH: not used ?
0027     unsigned            hd_factor ;
0028
0029 #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
0030     QSCINT_METHOD void    generate( sphoton& p, RNG& rng, const quad6& gs, unsigned long long photon_id, int genstep_id ) const ;
0031     QSCINT_METHOD void    reemit(   sphoton& p, RNG& rng, float scintillationTime) const ;
0032     QSCINT_METHOD void    mom_pol_wavelength(sphoton& p, RNG& rng) const ;
0033     // sets direction, polarization and wavelength as needed by both generate and reemit
0034
0035     QSCINT_METHOD float   wavelength(     const float& u0) const ;
0036     QSCINT_METHOD float   wavelength_hd0( const float& u0) const ;
0037     QSCINT_METHOD float   wavelength_hd10(const float& u0) const ;
0038     QSCINT_METHOD float   wavelength_hd20(const float& u0) const ;
0039
0040 #endif
0041
0042 };
0043
0044
0045 #if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
0046
0047 #include "sscint.h"
0048
0049 /**
0050 qscint::generate_photon
0051 ------------------------
0052
0053 **/
0054
0055 inline QSCINT_METHOD void qscint::generate(
0056     sphoton& p,
0057     RNG& rng,
0058     const quad6& _gs,
0059     unsigned long long photon_id,
0060     int genstep_id ) const
0061 {
0062     mom_pol_wavelength(p, rng );
0063
0064     const sscint& gs = (const sscint&)_gs ;
0065
0066     float fraction = gs.charge == 0.f  ? 1.f : curand_uniform(&rng) ;
0067     p.pos = gs.pos + fraction*gs.DeltaPosition ;
0068
0069     float u4 = curand_uniform(&rng) ;
0070     float deltaTime = fraction*gs.step_length/gs.meanVelocity - gs.ScintillationTime*logf(u4) ;
0071
0072     p.time = gs.time + deltaTime ;
0073     p.zero_flags();
0074     p.set_flag(SCINTILLATION) ;
0075 }
0076
0077
0078 /**
0079 qscint::reemit_photon
0080 ------------------------
0081
0082 OLD NOTES IN NEED OF REVIST : HOW TO HANDLE REEMISSION scintillationTime ?
0083
0084 As reemission happens inside scintillators for photons arising from Cerenkov (or Torch)
0085 gensteps need to special case the handing of the reemission scintillationTime somehow
0086 because do not have access to scintillation gensteps when handling cerenkov or torch photons.
0087
0088 Could carry the single float (could be domain compressed, it is eg 1.5 ns) in other gensteps ?
0089 But it is material specific (if you had more than one scintillator)
0090 just like REEMISSIONPROB so its more appropriate
0091 to live in the boundary_tex alongside the REEMISSIONPROB ?
0092
0093 But it could be carried in the genstep(or anywhere) as its use is "gated" by a non-zero REEMISSIONPROB.
0094
0095 Prefer to just hold it in the context, and provide G4Opticks::setReemissionScintillationTime API
0096 for setting it (default 0.) that is used from detector specific code which can read from
0097 the Geant4 properties directly.  What about geocache ? Can hold/persist with GScintillatorLib metadata.
0098
0099 **/
0100
0101 inline QSCINT_METHOD void qscint::reemit(sphoton& p, RNG& rng, float scintillationTime) const
0102 {
0103     mom_pol_wavelength(p, rng);
0104     float u3 = curand_uniform(&rng) ;
0105     p.time += -scintillationTime*logf(u3) ;
0106 }
0107
0108
0109 /**
0110 qscint::mom_pol_wavelength : dir,pol and wavelength do not depend on genstep param
0111 --------------------------------------------------------------------------------
0112
0113 Called from qscint::generate, a translation of "jcv DsG4Scintillation"
0114
0115 **/
0116
0117 inline QSCINT_METHOD void qscint::mom_pol_wavelength(sphoton& p, RNG& rng) const
0118 {
0119     float u0 = curand_uniform(&rng);
0120     float u1 = curand_uniform(&rng);
0121     float u2 = curand_uniform(&rng);
0122     float u3 = curand_uniform(&rng);
0123
0124     float cost = 1.f - 2.f*u0;
0125     float sint = sqrt((1.f-cost)*(1.f+cost));
0126     float phi = 2.f*M_PIf*u1;
0127     float sinp = sin(phi);
0128     float cosp = cos(phi);
0129
0130     p.mom.x = sint*cosp;
0131     p.mom.y = sint*sinp;
0132     p.mom.z = cost ;
0133
0134     // Determine polarization of new photon
0135     p.pol.x = cost*cosp ;
0136     p.pol.y = cost*sinp ;
0137     p.pol.z = -sint ;
0138
0139     phi = 2.f*M_PIf*u2 ;
0140     sinp = sin(phi);
0141     cosp = cos(phi);
0142
0143     p.pol = normalize( cosp*p.pol + sinp*cross(p.mom, p.pol) ) ;
0144     p.wavelength = wavelength(u3);
0145 }
0146
0147
0148
0149 inline QSCINT_METHOD float qscint::wavelength(const float& u0) const
0150 {
0151     float wl ;
0152     switch(hd_factor)
0153     {
0154         case 0:  wl = wavelength_hd0(u0)  ; break ;
0155         case 10: wl = wavelength_hd10(u0) ; break ;
0156         case 20: wl = wavelength_hd20(u0) ; break ;
0157         default: wl = 0.f ;
0158     }
0159     //printf("//qscint::wavelength wl %10.4f hd %d \n", wl, hd_factor );
0160     return wl ;
0161 }
0162
0163
0164 inline QSCINT_METHOD float qscint::wavelength_hd0(const float& u0) const
0165 {
0166     constexpr float y0 = 0.5f/3.f ;
0167     return tex2D<float>(scint_tex, u0, y0 );
0168 }
0169
0170 /**
0171 qscint::wavelength_hd10
0172 --------------------------------------------------
0173
0174 Idea is to improve handling of extremes by throwing ten times the bins
0175 at those regions, using simple and cheap linear mappings.
0176
0177 TODO: move hd "layers" into float4 payload so the 2d cerenkov and 1d scint
0178 icdf texture can share some of teh implementation
0179
0180 **/
0181
0182 inline QSCINT_METHOD float qscint::wavelength_hd10(const float& u0) const
0183 {
0184     float wl ;
0185
0186     constexpr float y0 = 0.5f/3.f ;
0187     constexpr float y1 = 1.5f/3.f ;
0188     constexpr float y2 = 2.5f/3.f ;
0189
0190     if( u0 < 0.1f )
0191     {
0192         wl = tex2D<float>(scint_tex, u0*10.f , y1 );
0193     }
0194     else if ( u0 > 0.9f )
0195     {
0196         wl = tex2D<float>(scint_tex, (u0 - 0.9f)*10.f , y2 );
0197     }
0198     else
0199     {
0200         wl = tex2D<float>(scint_tex, u0,  y0 );
0201     }
0202     return wl ;
0203 }
0204
0205
0206
0207 inline QSCINT_METHOD float qscint::wavelength_hd20(const float& u0) const
0208 {
0209     float wl ;
0210
0211     constexpr float y0 = 0.5f/3.f ;
0212     constexpr float y1 = 1.5f/3.f ;
0213     constexpr float y2 = 2.5f/3.f ;
0214
0215     if( u0 < 0.05f )
0216     {
0217         wl = tex2D<float>(scint_tex, u0*20.f , y1 );
0218     }
0219     else if ( u0 > 0.95f )
0220     {
0221         wl = tex2D<float>(scint_tex, (u0 - 0.95f)*20.f , y2 );
0222     }
0223     else
0224     {
0225         wl = tex2D<float>(scint_tex, u0,  y0 );
0226     }
0227     return wl ;
0228 }
0229
0230
0231 #endif
0232
0233