eic-opticks/sysrap/scerenkov.h

0001 #pragma once
0002 /**
0003 scerenkov.h : replace (but stay similar to) : npy/NStep.hpp optixrap/cu/cerenkovstep.h
0004 ========================================================================================
0005
0006 * FOLLOWING PATTERN OF storch.h
0007
0008 **/
0009
0010
0011 #if defined(__CUDACC__) || defined(__CUDABE__)
0012    #define SCERENKOV_METHOD __device__
0013 #else
0014    #define SCERENKOV_METHOD
0015 #endif
0016
0017 #include "OpticksGenstep.h"
0018 #include "OpticksPhoton.h"
0019
0020 #include "smath.h"
0021 #include "scuda.h"
0022 #include "squad.h"
0023
0024 #if defined(__CUDACC__) || defined(__CUDABE__)
0025 #else
0026 #include "NP.hh"
0027 #endif
0028
0029 struct scerenkov
0030 {
0031     // ctrl
0032     unsigned gentype ;   // formerly Id
0033     unsigned trackid ;   // formerly ParentId
0034     unsigned matline ;   // formerly MaterialIndex, used by qbnd::boundary_lookup
0035     unsigned numphoton ; // formerly NumPhotons
0036
0037     float3   pos ;  // formerly x0
0038     float    time ; // formerly t0
0039
0040     float3 DeltaPosition ;  // aka p0 in G4Cerenkov,  G4Step::GetDeltaPosition is not normalized
0041     float  step_length ;
0042
0043     int code;
0044     float charge ;
0045     float weight ;
0046     float preVelocity ;
0047
0048     /// the above first 4 quads are common to both CerenkovStep and ScintillationStep
0049
0050     float BetaInverse ;
0051     float Wmin ;
0052     float Wmax ;
0053     float maxCos ;
0054
0055     float maxSin2 ;
0056     float MeanNumberOfPhotons1 ;
0057     float MeanNumberOfPhotons2 ;
0058     float postVelocity ;
0059
0060     SCERENKOV_METHOD float Pmin() const { return smath::hc_eVnm/Wmax ; }
0061     SCERENKOV_METHOD float Pmax() const { return smath::hc_eVnm/Wmin ; }
0062
0063    static void FillGenstep( scerenkov& gs, int matline, int numphoton_per_genstep, bool dump ) ;
0064
0065 #if defined(__CUDACC__) || defined(__CUDABE__)
0066 #else
0067    float* cdata() const {  return (float*)&gentype ; }
0068    std::string desc() const ;
0069    static void MinMaxPost( float* mn, float* mx, const NP* genstep);
0070
0071    template<typename T>
0072    static bool IsGenstepArray( const NP* a );
0073
0074 #endif
0075
0076 };
0077
0078
0079
0080 /**
0081 scerenkov::FillGenstep : fabricate some values for a demo genstep
0082 ---------------------------------------------------------------------
0083
0084 Uses hard coded values depending on RINDEX of material (LS) that will be used : ie fixing the cone angle.
0085 A better way of doing this would use the MaterialLine as input and obtain the values from the RINDEX,
0086
0087 * (as this code is only needed on CPU only it would be perfectly feasible to read RINDEX arrays)
0088 * however this is debugging code so non-general techniques are acceptable
0089 * NB matline is crucial as that determines which materials RINDEX is used
0090
0091 **/
0092
0093 inline void scerenkov::FillGenstep( scerenkov& gs, int matline, int numphoton_per_genstep, bool dump )
0094 {
0095     float nMax = 1.793f ;
0096     float BetaInverse = 1.500f ;
0097     float maxCos = BetaInverse / nMax;
0098     float maxSin2 = (1.f - maxCos) * (1.f + maxCos) ;
0099
0100     float Pmin = 1.55f ;    // eV
0101     float Pmax = 15.5f ;
0102     float Wmin = smath::hc_eVnm/Pmax ; // close to: 1240./15.5 = 80.
0103     float Wmax = smath::hc_eVnm/Pmin ;  // close to: 1240./1.55 = 800.
0104     // NB in reality this should not use standard domain, but rather the domain of the RINDEX property
0105
0106     gs.gentype = OpticksGenstep_CERENKOV ;
0107     gs.trackid = 0u ;
0108     gs.matline = matline ;
0109     gs.numphoton = numphoton_per_genstep  ;
0110
0111     gs.pos.x = 100.f ;
0112     gs.pos.y = 100.f ;
0113     gs.pos.z = 100.f ;
0114     gs.time = 20.f ;
0115
0116     gs.DeltaPosition.x = 1000.f ;    // aka p0
0117     gs.DeltaPosition.y = 1000.f ;
0118     gs.DeltaPosition.z = 1000.f ;
0119     gs.step_length = 1000.f ;
0120
0121     gs.code = 1 ;
0122     gs.charge = 1.f ;
0123     gs.weight = 1.f ;
0124     gs.preVelocity = 10.f ;
0125
0126     gs.BetaInverse = BetaInverse ;
0127     gs.Wmin = Wmin ;
0128     gs.Wmax = Wmax ;
0129     gs.maxCos = maxCos ;  // NOT USED ?
0130
0131     gs.maxSin2 = maxSin2 ;              // constrains cone angle rejection sampling
0132     gs.MeanNumberOfPhotons1 = 100.f ;   // used for profile sampling to decide where along the step
0133     gs.MeanNumberOfPhotons2 = 200.f ;
0134     gs.postVelocity = 20.f ;
0135
0136 }
0137
0138
0139
0140 #if defined(__CUDACC__) || defined(__CUDABE__)
0141 #else
0142
0143 inline std::string scerenkov::desc() const
0144 {
0145     std::stringstream ss ;
0146     ss << "scerenkov::desc"
0147        << " gentype " << gentype
0148        ;
0149     std::string s = ss.str();
0150     return s ;
0151 }
0152
0153
0154 inline void scerenkov::MinMaxPost( float* mn, float* mx, const NP* _a )
0155 {
0156     NP* a = const_cast<NP*>(_a);
0157
0158     bool is_f = IsGenstepArray<float>(a);
0159
0160     if(!is_f) std::cerr
0161         << "scerenkov::MinMaxPost FATAL UNEXPECTED ARRAY "
0162         << " is_f " << ( is_f ? "YES" : "NO " )
0163         << "\n"
0164         ;
0165
0166     bool expect = is_f ;
0167     assert(expect);
0168     if(!expect) std::raise(SIGINT);
0169
0170     int ni = a->num_items() ;
0171     int nj = 4 ;
0172
0173     float MAX = std::numeric_limits<float>::max() ;
0174     for(int j=0 ; j < nj ; j++)
0175     {
0176         mn[j] = MAX ;
0177         mx[j] = -MAX ;
0178     }
0179
0180     for(int i=0 ; i < ni ; i++)
0181     {
0182         scerenkov* gs = reinterpret_cast<scerenkov*>(a->bytes() + i*a->item_bytes());
0183         assert( gs );
0184
0185         //const float3& dpos = gs->DeltaPosition ;
0186         float4 _p0 = make_float4( gs->pos.x, gs->pos.y, gs->pos.z, gs->time );
0187
0188         float* p0 = &_p0.x ;
0189
0190         for(int j=0 ; j < 4 ; j++)
0191         {
0192             float vj = p0[j];
0193             if( vj < mn[j] ) mn[j] = vj ;
0194             if( vj > mx[j] ) mx[j] = vj ;
0195         }
0196     }
0197
0198
0199
0200 }
0201
0202 template<typename T>
0203 inline bool scerenkov::IsGenstepArray( const NP* a )
0204 {
0205     return a && a->has_shape(-1,6,4) && a->ebyte == sizeof(T) ;
0206 }
0207
0208
0209 #endif
0210