eic-opticks/sysrap/storch.h

0001 #pragma once
0002 /**
0003 storch.h : replace (but stay similar to) : npy/NStep.hpp optixrap/cu/torchstep.h
0004 ===================================================================================
0005
0006 NB sizeof storch struct is **CONSTRAINED TO MATCH quad6** like all gensteps
0007
0008 Bringing over some of the old torch genstep generation into the modern workflow
0009 with mocking on CPU and pure-CUDA test cababilities.
0010
0011 Notes
0012 --------
0013
0014 Techniques to implement the spirit of the old torch gensteps in much less code
0015
0016 * sharing types and code between GPU and CPU
0017 * quad6 and NP and casting between them
0018 * union between quad6 and simple torch struct eliminates tedious get/set of NStep.hpp
0019 * macros to use same headers on CPU and GPU, eg enum strings live with enum values in same header
0020   but are hidden from nvcc
0021
0022 **/
0023
0024 #if defined(__CUDACC__) || defined(__CUDABE__)
0025    #define STORCH_METHOD __device__
0026 #else
0027    #define STORCH_METHOD inline
0028 #endif
0029
0030
0031
0032 #include "OpticksGenstep.h"
0033 #include "OpticksPhoton.h"
0034
0035 //#include "scurand.h"
0036 #include "smath.h"
0037 #include "sphoton.h"
0038 #include "storchtype.h"
0039
0040 /**
0041 **/
0042
0043 struct storch
0044 {
0045     // ctrl
0046     unsigned gentype ;  // eg OpticksGenstep_TORCH
0047     unsigned trackid ;
0048     unsigned matline ;
0049     unsigned numphoton ;
0050
0051     float3   pos ;
0052     float    time ;
0053
0054     float3   mom ;
0055     float    weight ;
0056
0057     float3   pol ;
0058     float    wavelength ;
0059
0060     float2  zenith ;  // for T_RECTANGLE : repurposed for the Z values of rect sides
0061     float2  azimuth ; // for T_RECTANGLE : repurposed for the X values of rect sides
0062
0063     // beam
0064     float    radius ;
0065     float    distance ;
0066     unsigned mode ;     // basemode
0067     unsigned type ;     // basetype
0068
0069     // NB : organized into 6 quads : are constained not to change that
0070
0071     STORCH_METHOD static void generate( sphoton& p, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id );
0072
0073
0074 #if defined(__CUDACC__) || defined(__CUDABE__)
0075 #else
0076    float* cdata() const {  return (float*)&gentype ; }
0077    static constexpr const char* storch_FillGenstep_pos        = "storch_FillGenstep_pos" ;
0078    static constexpr const char* storch_FillGenstep_time       = "storch_FillGenstep_time" ;
0079    static constexpr const char* storch_FillGenstep_mom        = "storch_FillGenstep_mom" ;
0080    static constexpr const char* storch_FillGenstep_wavelength = "storch_FillGenstep_wavelength" ;
0081    static constexpr const char* storch_FillGenstep_distance   = "storch_FillGenstep_distance" ;
0082    static constexpr const char* storch_FillGenstep_weight     = "storch_FillGenstep_weight" ;
0083    static constexpr const char* storch_FillGenstep_radius     = "storch_FillGenstep_radius" ;
0084    static constexpr const char* storch_FillGenstep_zenith     = "storch_FillGenstep_zenith" ;
0085    static constexpr const char* storch_FillGenstep_azimuth    = "storch_FillGenstep_azimuth" ;
0086    static constexpr const char* storch_FillGenstep_type       = "storch_FillGenstep_type" ;
0087    static void FillGenstep( storch& gs, int genstep_id, int numphoton_per_genstep, bool dump=false ) ;
0088    std::string desc() const ;
0089 #endif
0090
0091 };
0092
0093
0094 #if defined(__CUDACC__) || defined(__CUDABE__)
0095 #else
0096
0097 /**
0098 storch::FillGenstep
0099 ----------------------
0100
0101 Canonically invoked from SEvent::MakeGensteps
0102
0103 **/
0104 inline void storch::FillGenstep( storch& gs, int genstep_id, int numphoton_per_genstep, bool dump )
0105 {
0106     gs.gentype = OpticksGenstep_TORCH ;
0107     gs.numphoton = numphoton_per_genstep  ;
0108
0109     qvals( gs.pos , storch_FillGenstep_pos , "0,0,-90" );
0110     if(dump) printf("//storch::FillGenstep storch_FillGenstep_pos gs.pos (%10.4f %10.4f %10.4f) \n", gs.pos.x, gs.pos.y, gs.pos.z );
0111
0112     qvals( gs.time, storch_FillGenstep_time, "0.0" );
0113     if(dump) printf("//storch::FillGenstep storch_FillGenstep_time gs.time (%10.4f) \n", gs.time );
0114
0115     qvals( gs.mom , storch_FillGenstep_mom , "0,0,1" );
0116     gs.mom = normalize(gs.mom);  // maybe should skip this float normalized, relying instead on U4VPrimaryGenerator::GetPhotonParam to do the normalize ?
0117     if(dump) printf("//storch::FillGenstep storch_FillGenstep_mom gs.mom (%10.4f %10.4f %10.4f) \n", gs.mom.x, gs.mom.y, gs.mom.z );
0118
0119     qvals( gs.wavelength, storch_FillGenstep_wavelength, "420" );
0120     if(dump) printf("//storch::FillGenstep storch_FillGenstep_wavelength gs.wavelength (%10.4f) \n", gs.wavelength  );
0121
0122     qvals( gs.distance, storch_FillGenstep_distance, "0" );
0123     if(dump) printf("//storch::FillGenstep storch_FillGenstep_distance gs.distance (%10.4f) \n", gs.distance  );
0124
0125     qvals( gs.weight, storch_FillGenstep_weight, "0" );
0126     if(dump) printf("//storch::FillGenstep storch_FillGenstep_weight gs.weight (%10.4f) \n", gs.weight  );
0127
0128     qvals( gs.zenith,  storch_FillGenstep_zenith,  "0,1" );
0129     if(dump) printf("//storch::FillGenstep storch_FillGenstep_zenith gs.zenith (%10.4f,%10.4f) \n", gs.zenith.x, gs.zenith.y  );
0130
0131     qvals( gs.azimuth,  storch_FillGenstep_azimuth,  "0,1" );
0132     if(dump) printf("//storch::FillGenstep storch_FillGenstep_azimuth gs.azimuth (%10.4f,%10.4f) \n", gs.azimuth.x, gs.azimuth.y  );
0133
0134     qvals( gs.radius, storch_FillGenstep_radius, "50" );
0135     if(dump) printf("//storch::FillGenstep storch_FillGenstep_radius gs.radius (%10.4f) \n", gs.radius );
0136
0137     const char* type = qenv(storch_FillGenstep_type, "disc" );
0138     unsigned ttype = storchtype::Type(type) ;
0139     bool ttype_valid = storchtype::IsValid(ttype) ;
0140     if(!ttype_valid) printf("//storch::FillGenstep FATAL TTYPE INVALID %s:[%s][%d] \n", storch_FillGenstep_type, type, ttype ) ;
0141     assert(ttype_valid);
0142
0143     gs.type = ttype ;
0144     if(dump) printf("//storch::FillGenstep storch_FillGenstep_type %s  gs.type %d \n", type, gs.type );
0145     gs.mode = 255 ;    //torchmode::Type("...");
0146 }
0147
0148
0149 inline std::string storch::desc() const
0150 {
0151     std::stringstream ss ;
0152     ss << "storch::desc"
0153        << " gentype " << gentype
0154        << " mode " << mode
0155        << " type " << type
0156        ;
0157     std::string s = ss.str();
0158     return s ;
0159 }
0160 #endif
0161
0162
0163
0164 /**
0165 storch::generate
0166 -----------------
0167
0168 On GPU this is invoked by::
0169
0170    CSGOptiX7.cu:simulate
0171    qsim::generate_photon
0172
0173 On CPU this is invoked using MOCK_CURAND with for example::
0174
0175    G4CXApp::GeneratePrimaries
0176    U4VPrimaryGenerator::GeneratePrimaries
0177    SGenerate::GeneratePhotons
0178    storch::generate
0179
0180
0181 Populate "sphoton& p" as parameterized by "const quad6& gs_" which casts to "const storch& gs",
0182 the photon_id and genstep_id inputs are informational only.
0183
0184 Old workflow equivalent ~/opticks/optixrap/cu/torchstep.h
0185
0186
0187 **/
0188
0189 STORCH_METHOD void storch::generate( sphoton& p, RNG& rng, const quad6& gs_, unsigned long long photon_id, unsigned genstep_id )  // static
0190 {
0191     const storch& gs = (const storch&)gs_ ;   // casting between union-ed types : quad6 and storch
0192
0193 #ifdef STORCH_DEBUG
0194     printf("//storch::generate photon_id %3d genstep_id %3d  gs gentype/trackid/matline/numphoton(%3d %3d %3d %3d) type %d \n",
0195        photon_id,
0196        genstep_id,
0197        gs.gentype,
0198        gs.trackid,
0199        gs.matline,
0200        gs.numphoton,
0201        gs.type
0202       );
0203 #endif
0204     if( gs.type == T_DISC )
0205     {
0206
0207         /**
0208         disc/T_DISC
0209             zenith.x->zenith.y radial range, eg [0. 100.] filled disc, [90., 100.] annulus
0210             azimuth.x->azimuth.y phi segment in fraction of twopi [0,1] for complete segment
0211
0212         **/
0213         //printf("//storch::generate T_DISC gs.type %d gs.mode %d  \n", gs.type, gs.mode );
0214
0215         p.wavelength = gs.wavelength ;
0216         p.time = gs.time ;
0217         p.mom = gs.mom ;
0218
0219         float u_zenith  = gs.zenith.x  + curand_uniform(&rng)*(gs.zenith.y-gs.zenith.x)   ;
0220         float u_azimuth = gs.azimuth.x + curand_uniform(&rng)*(gs.azimuth.y-gs.azimuth.x) ;
0221
0222
0223
0224         float r = gs.radius*u_zenith ;
0225
0226         float phi = 2.f*M_PIf*u_azimuth ;
0227         float sinPhi = sinf(phi);
0228         float cosPhi = cosf(phi);
0229         // __sincosf(phi,&sinPhi,&cosPhi);   // HMM: think thats an apple extension
0230
0231         p.pos.x = r*cosPhi ;
0232         p.pos.y = r*sinPhi ;
0233         p.pos.z = 0.f ;
0234         // 3D rotate the positions to make their disc perpendicular to p.mom for a nice beam
0235         smath::rotateUz(p.pos, p.mom) ;
0236         p.pos = p.pos + gs.pos ; // translate position after orienting the disc
0237
0238         p.pol.x = sinPhi ;
0239         p.pol.y = -cosPhi ;
0240         p.pol.z = 0.f ;
0241         // p.pol.z zero in initial frame, so rotating the frame to arrange
0242         // z to be in p.mom direction makes pol transverse to mom
0243         // NOTICE : ARBITRARY PHASE OF POLARIZATION CHOICE HERE
0244         smath::rotateUz(p.pol, p.mom) ;
0245     }
0246     else if( gs.type == T_SPHERE )
0247     {
0248         /**
0249         T_SPHERE
0250              generates positions on a sphere of gs.radius and radial momentum direction
0251              outwards(inwards) for gs.radius +ve(-ve)
0252         **/
0253
0254
0255         p.wavelength = gs.wavelength ;
0256         p.time = gs.time ;
0257
0258         float u_zenith  = gs.zenith.x  + curand_uniform(&rng)*(gs.zenith.y-gs.zenith.x)   ;
0259         float u_azimuth = gs.azimuth.x + curand_uniform(&rng)*(gs.azimuth.y-gs.azimuth.x) ;
0260
0261         float phi = 2.f*M_PIf*u_azimuth ;  // azimuth range 0->2pi
0262         float sinPhi = sinf(phi);
0263         float cosPhi = cosf(phi);
0264
0265         float cosTheta = 1.f - 2.0f*u_zenith  ;   // polar range 0->pi
0266         float sinTheta = sqrtf( 1.0f - cosTheta*cosTheta );
0267
0268         float flip = copysignf( 1.f, gs.radius );
0269
0270         // gs.radius -ve(+ve) => flip -1.f(+1.f)
0271         // below flips direction and not position for outwards/inwards control
0272
0273         p.mom.x = flip*sinTheta*cosPhi ;
0274         p.mom.y = flip*sinTheta*sinPhi ;
0275         p.mom.z = flip*cosTheta ;
0276
0277         float radius = fabs(gs.radius);
0278         p.pos.x = sinTheta*cosPhi*radius ;
0279         p.pos.y = sinTheta*sinPhi*radius ;
0280         p.pos.z = cosTheta*radius ;
0281
0282         // float frac_twopi = 0.0f ;    // tangent vectors up the sphere (towards +Z pole, increasing theta)
0283         // float frac_twopi = 0.5f ;    // tangent vectors down the sphere (towards -Z pole, decreasing theta)
0284         // float frac_twopi = 0.25f ;   // tangents around the sphere in direction of increasing phi
0285         // float frac_twopi = 0.75f ;      // tangents around the sphere in direction of decreasing phi
0286         // NOTICE : ARBITRARY PHASE OF POLARIZATION CHOICE HERE : VARIOUS TANGENTS TO THE SPHERE
0287         float frac_twopi = gs.distance ;  // repurpose the distance, as frac_twopi
0288
0289         float phase = 2.f*M_PIf*frac_twopi ;
0290         p.pol.x = cosf(phase) ;
0291         p.pol.y = sinf(phase) ;
0292         p.pol.z = 0.f ;
0293         // p.pol.z zero in initial frame, so rotating the frame to arrange
0294         // z to be in p.mom direction makes pol transverse to mom
0295         smath::rotateUz(p.pol, p.mom);
0296     }
0297     else if( gs.type == T_SPHERE_MARSAGLIA )
0298     {
0299         /**
0300         T_SPHERE_MARSAGLIA
0301              generates positions on a sphere of gs.radius and radial momentum direction
0302              outwards(inwards) for gs.radius +ve(-ve)
0303
0304              uses Marsaglia rejection sampling to get points on unit sphere
0305
0306              using zenith/azimuth does restrict the range, but in a funny tent shape
0307
0308         **/
0309         p.wavelength = gs.wavelength ;
0310         p.time = gs.time ;
0311
0312         float u0_zenith, u1_azimuth  ;
0313         float u, v, b, a  ;
0314
0315         do
0316         {
0317             u0_zenith  = gs.zenith.x  + curand_uniform(&rng)*(gs.zenith.y-gs.zenith.x)   ;   // aka polar theta
0318             u1_azimuth = gs.azimuth.x + curand_uniform(&rng)*(gs.azimuth.y-gs.azimuth.x) ;   // aka azimuth phi
0319
0320             u = 2.f*u0_zenith - 1.f ;
0321             v = 2.f*u1_azimuth - 1.f ;
0322             b = u*u + v*v ;
0323         }
0324         while( b > 1.f ) ;
0325         a = 2.f*sqrtf( 1.f - b );
0326
0327 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0328         //printf("//storch::generate T_SPHERE gs.radius %10.4f gs.distance %10.4f \n", gs.radius, gs.distance );
0329 #endif
0330         float radius = fabs(gs.radius) ;
0331         float flip = copysignf( 1.f, gs.radius );
0332
0333         // gs.radius -ve(+ve) => flip -1.f(+1.f)
0334         // want to flip direction but not position and avoid extra storage
0335         p.mom.x = flip*a*u ;
0336         p.mom.y = flip*a*v ;
0337         p.mom.z = flip*(2.f*b - 1.f) ;
0338
0339         p.pos.x = a*u*radius ;
0340         p.pos.y = a*v*radius ;
0341         p.pos.z = (2.f*b-1.f)*radius ;
0342
0343         // float frac_twopi = 0.0f ;    // tangent vectors up the sphere (towards +Z pole, increasing theta)
0344         // float frac_twopi = 0.5f ;    // tangent vectors down the sphere (towards -Z pole, decreasing theta)
0345         // float frac_twopi = 0.25f ;   // tangents around the sphere in direction of increasing phi
0346         // float frac_twopi = 0.75f ;      // tangents around the sphere in direction of decreasing phi
0347         // NOTICE : ARBITRARY PHASE OF POLARIZATION CHOICE HERE : VARIOUS TANGENTS TO THE SPHERE
0348         float frac_twopi = gs.distance ;  // repurpose the distance, as frac_twopi
0349
0350         float phase = 2.f*M_PIf*frac_twopi ;
0351         p.pol.x = cosf(phase) ;
0352         p.pol.y = sinf(phase) ;
0353         p.pol.z = 0.f ;
0354         // p.pol.z zero in initial frame, so rotating the frame to arrange
0355         // z to be in p.mom direction makes pol transverse to mom
0356         smath::rotateUz(p.pol, p.mom);
0357     }
0358     else if( gs.type == T_LINE )
0359     {
0360         /**
0361         T_LINE
0362             photons start at positions (varying by photon_id)
0363             along a line from -gs.radius to +gs.radius
0364             wavelength and time are fixed.
0365             The position and polarization are oriented
0366             such that the local frame z is in the gs.mom direction.
0367
0368         **/
0369         p.wavelength = gs.wavelength ;
0370         p.time = gs.time ;
0371         p.mom = gs.mom ;
0372
0373         float frac = float(photon_id)/float(gs.numphoton) ;  // 0->~1
0374         float sfrac = 2.f*(frac-0.5f) ;     // -1 -> ~1
0375         float r = gs.radius*sfrac ;        // -gs.radius -> gs.radius  (NB gets offset by gs.pos too)
0376
0377         p.pos.x = r ;
0378         p.pos.y = 0.f ;
0379         p.pos.z = 0.f ;
0380
0381         smath::rotateUz(p.pos, p.mom) ;
0382         p.pos = p.pos + gs.pos ; // translate position after orienting the line
0383
0384         p.pol.x = 0.f ;
0385         p.pol.y = -1.f ;
0386         p.pol.z = 0.f ;
0387         smath::rotateUz(p.pol, p.mom) ;
0388     }
0389     else if( gs.type == T_POINT )
0390     {
0391         /**
0392         T_POINT
0393              all photons start at gs.pos, local frame z is rotate into
0394              gs.mom direction as is local -Y polarization direction
0395         **/
0396         p.wavelength = gs.wavelength ;
0397         p.time = gs.time ;
0398         p.mom = gs.mom ;
0399
0400         p.pos.x = 0.f ;
0401         p.pos.y = 0.f ;
0402         p.pos.z = 0.f ;
0403
0404         smath::rotateUz(p.pos, p.mom) ;
0405         p.pos = p.pos + gs.pos ; // translate position after orienting the line
0406
0407         p.pol.x = 0.f ;
0408         p.pol.y = -1.f ;
0409         p.pol.z = 0.f ;
0410         smath::rotateUz(p.pol, p.mom) ;
0411     }
0412     else if( gs.type == T_CIRCLE )
0413     {
0414         /**
0415         T_CIRCLE
0416              phi position around circle of radius |gs.radius| based on photon_id
0417              at position gs.pos
0418
0419              gs.radius>0(<0)
0420                  local mom is radially outwards(inwards) in XZ plane
0421
0422              local -Y pol direction is oriented according to the
0423              local radial mom direction
0424
0425         **/
0426         p.wavelength = gs.wavelength ;
0427         p.time = gs.time ;
0428         float f = float(photon_id)/float(gs.numphoton) ;      // 0->~1
0429         float frac = gs.azimuth.x*(1.f-f) + gs.azimuth.y*(f) ; // gs.azimuth.x -> gs.azimuth.y which defaults to 0->1
0430
0431         float phi = 2.f*M_PIf*frac ;
0432         float sinPhi = sinf(phi);
0433         float cosPhi = cosf(phi);
0434
0435         float r = gs.radius < 0.f ? -gs.radius : gs.radius ;
0436
0437         // -ve radius for inwards rays
0438         // +ve radius or zero for outwards rays
0439
0440         p.mom.x = gs.radius < 0.f ? -cosPhi : cosPhi ;
0441         p.mom.y = 0.f ;
0442         p.mom.z = gs.radius < 0.f ? -sinPhi : sinPhi ;
0443
0444         p.pos.x = r*cosPhi ;
0445         p.pos.y = 0.f ;
0446         p.pos.z = r*sinPhi ;
0447         // smath::rotateUz(p.pos, p.mom) ;  // dont do that that
0448         p.pos = p.pos + gs.pos ; // translate position after orienting the line
0449
0450         /*
0451         printf("// T_CIRCLE frac %10.4f gs.radius %10.4f r %10.4f  p.mom (%10.4f %10.4f %10.4f) p.pos (%10.4f %10.4f %10.4f) \n",
0452            frac, gs.radius, r, p.mom.x, p.mom.y, p.mom.z, p.pos.x, p.pos.y, p.pos.z );
0453         */
0454
0455         p.pol.x = 0.f ;
0456         p.pol.y = -1.f ;
0457         p.pol.z = 0.f ;
0458         smath::rotateUz(p.pol, p.mom) ;
0459     }
0460     else if( gs.type == T_RECTANGLE )
0461     {
0462         /**
0463         DIVIDE TOTAL PHOTON SLOTS INTO FOUR SIDES
0464         ie side is 0,0,0...,1,1,1...,2,2,2,..,3,3,3...
0465
0466               +------3:top-------+- gs.zenith.y
0467               |                  |
0468               |                  |
0469               0:left             1:right
0470               |                  |
0471               |                  |
0472               +---2:bottom-------+- gs.zenith.x
0473               |                  |
0474               gs.azimuth.x       gs.azimuth.y
0475
0476         **/
0477
0478         p.wavelength = gs.wavelength ;
0479         p.time = gs.time ;
0480
0481         int side_size = gs.numphoton/4 ;
0482         int side = photon_id/side_size ;
0483         int side_offset = side*side_size ;
0484         int side_index = photon_id - side_offset ; // index within the side
0485         float frac = float(side_index)/float(side_size) ;  // 0->~1  within the side
0486
0487         if( side == 0 || side == 1 ) // left or right
0488         {
0489             p.pos.x = side == 0 ? gs.azimuth.x : gs.azimuth.y ;
0490             p.pos.y = 0.f ;
0491             p.pos.z = (1.f-frac)*gs.zenith.x + frac*gs.zenith.y ;
0492
0493             p.mom.x = side == 0 ? 1.f : -1.f ;   // inwards
0494             p.mom.y = 0.f ;
0495             p.mom.z = 0.f ;
0496         }
0497         else if( side == 2 || side == 3)  // bottom or top
0498         {
0499             p.pos.x = (1.f-frac)*gs.azimuth.x + frac*gs.azimuth.y ;
0500             p.pos.y = 0.f ;
0501             p.pos.z = side == 2 ? gs.zenith.x : gs.zenith.y ;
0502
0503             p.mom.x = 0.f ;
0504             p.mom.y = 0.f ;
0505             p.mom.z = side == 2 ? 1.f : -1.f ;
0506         }
0507         p.pos = p.pos + gs.pos ; // translate position, often gs.pos is origin anyhow
0508
0509         p.pol.x = 0.f ;
0510         p.pol.y = -1.f ;    // point out the XZ plane, so its transverse
0511         p.pol.z = 0.f ;
0512         smath::rotateUz(p.pol, p.mom) ;
0513     }
0514     p.zero_flags();
0515     p.set_flag(TORCH);
0516 }
0517
0518
0519
0520
0521 /**
0522 * qtorch : union between quad6 and specific genstep types for easy usage and yet no serialize/deserialize needed
0523 **/
0524
0525 union qtorch
0526 {
0527    quad6  q ;
0528    storch t ;
0529 };
0530
0531
0532