eic-opticks/sysrap/sphoton.h

0001 #pragma once
0002 /**
0003 sphoton.h
0004 ============
0005
0006 +----+-------------------------+------------------+-----------------+---------------------+------------------------------+
0007 | q  |               x         |        y         |      z          |           w         |  notes                       |
0008 +====+=========================+==================+=================+=====================+==============================+
0009 |    |           pos.x         |    pos.y         |   pos.z         |       time          |                              |
0010 | q0 |                         |                  |                 |                     |                              |
0011 |    |                         |                  |                 |                     |                              |
0012 +----+-------------------------+------------------+-----------------+---------------------+------------------------------+
0013 |    |           mom.x         |    mom.y         |  mom.z          |   u:hitcount_iindex |                              |
0014 | q1 |                         |                  |                 |      16b       16b  |                              |
0015 |    |                         |                  |                 | (1,3)               |                              |
0016 +----+-------------------------+------------------+-----------------+---------------------+------------------------------+
0017 |    |           pol.x         |    pol.y         |   pol.z         |       wavelength    |                              |
0018 | q2 |                         |                  |                 |                     |                              |
0019 |    |                         |                  |                 |                     |                              |
0020 +----+-------------------------+------------------+-----------------+---------------------+------------------------------+
0021 |    | u:orient_boundary_flag  |  u:identity      | u:index         |     u:flagmask      |                              |
0022 | q3 |    1b     15b      16b  |    8b 24b        |                 |                     |                              |
0023 |    |                         |  8b:extend index |                 |                     |                              |
0024 |    | (3,0)                   | (3,1)            | (3,2)           | (3,3)               |                              |
0025 +----+-------------------------+------------------+-----------------+---------------------+------------------------------+
0026
0027
0028
0029 iindex
0030     instance index of intersected geometry
0031     (see squad.h:quad2 for details, and notes/issues/iindex_making_sense_of_it.rst)
0032
0033     In [11]: ii = f.hit[:,1,3].view(np.uint32)
0034
0035     In [12]: ii.min(),ii.max()
0036     Out[12]: (np.uint32(1), np.uint32(43196))
0037
0038     In [2]: ii = f.photon[:,1,3].view(np.uint32)
0039     In [3]: ii.min(),ii.max()
0040     Out[3]: (np.uint32(0), np.uint32(48593))
0041
0042
0043
0044 hitcount (16 bit)
0045     WIP:use for thrust based hit merging, unclear where best to set it
0046
0047 boundary (16 bit)
0048     boundary index of intersected geometry
0049     index corresponds to unique combinations of four material and surface indices
0050     (omat,osur,isur,imat)
0051
0052     [HMM, does not need 16bit]
0053
0054 flag (16 bit)
0055     OpticksPhoton.h history flag enum eg: TO CK SI BT BR SR AB RE ...
0056
0057 identity (32 bit)
0058     currently sensor_identifier+1 (see squad.h:quad2 for details)
0059     Formerly was combination of primIdx and instanceId of intersected geometry.
0060
0061     In [9]: id = f.hit[:,3,1].view(np.uint32)
0062
0063     In [10]: id.min(), id.max()
0064     Out[10]: (np.uint32(1), np.uint32(45600))
0065
0066
0067     In [3]: id = f.photon.view(np.uint32)[:,3,1]
0068
0069     In [4]: id.min(), id.max()
0070     Out[4]: (np.uint32(0), np.uint32(45600))
0071
0072     In [5]: np.c_[np.unique(id, return_counts=True)]
0073     Out[5]:
0074     array([[     0, 418619],    ## zero means did not end on a sensor
0075            [     1,     22],
0076            [     2,     22],
0077            [     3,     20],
0078            [     4,     18],
0079            ...,
0080            [ 45596,      4],
0081            [ 45597,      1],
0082            [ 45598,      1],
0083            [ 45599,      1],
0084            [ 45600,      2]], shape=(34580, 2))
0085
0086
0087     In [11]: 0xffff,0xffffff   ## 16/24 bits would be enough for identity, sparing 16/8 bits
0088     Out[11]: (65535, 16777215)
0089
0090
0091 orient (1 bit)
0092     set according to the sign of cosTheta
0093     (dot product of surface normal and momentum at last intersect)
0094
0095 index (32 bit)
0096     photon index, always exists even before any intersect
0097     (NB this is limited to 32 bits, maximum (unsigned)0xffffffff is 4.29 billion,
0098     so it will be clocked when simulating events with more photons than that.
0099
0100     As such large simulations are rare and the index is informational and can easily
0101     be post-corrected using uint64 arrays this clocking of the index will not be fixed.
0102
0103     HMM sparing 8 bits from another element would get to more than a trillion::
0104
0105             In [3]: (0x1 << (32+8))/1e9
0106             Out[3]: 1099.511627776
0107
0108
0109
0110 flagmask (32 bit)
0111     bitwise-OR of step point flag enumeration
0112     (see sysrap/OpticksPhoton.h sysrap/OpticksPhoton.hh for details)
0113     Always exists even before any intersect with the generation flag CK/SI/TO.
0114
0115
0116
0117
0118
0119 NB locations and packing here need to match ana/p.py
0120 ------------------------------------------------------
0121
0122 See ana/p.py for python accessors such as::
0123
0124     boundary_  = lambda p:p.view(np.uint32)[3,0] >> 16
0125     flag_      = lambda p:p.view(np.uint32)[3,0] & 0xffff
0126
0127     identity_ = lambda p:p.view(np.uint32)[3,1]
0128     primIdx_   = lambda p:identity_(p) >> 16
0129     instanceId_  = lambda p:identity_(p) & 0xffff
0130
0131     idx_      = lambda p:p.view(np.uint32)[3,2] & 0x7fffffff
0132     orient_   = lambda p:p.view(np.uint32)[3,2] >> 31
0133
0134     flagmask_ = lambda p:p.view(np.uint32)[3,3]
0135
0136
0137
0138 **/
0139
0140
0141
0142 #if defined(__CUDACC__) || defined(__CUDABE__)
0143 #    define SPHOTON_METHOD __host__ __device__ __forceinline__
0144 #else
0145 #    define SPHOTON_METHOD inline
0146 #endif
0147 #include <cstdint>
0148
0149 #if defined(__CUDACC__) || defined(__CUDABE__)
0150 #else
0151    #include <iostream>
0152    #include <iomanip>
0153    #include <sstream>
0154    #include <vector>
0155    #include <cstring>
0156    #include <csignal>
0157    #include <cassert>
0158    #include <glm/glm.hpp>
0159
0160    #include "scuda.h"
0161    #include "stran.h"
0162    #include "smath.h"
0163    #include "OpticksPhoton.h"
0164
0165    struct NP ;
0166 #endif
0167
0168
0169 // set_orient_iindex removed
0170
0171 struct sphoton
0172 {
0173 #if defined(__CUDACC__) || defined(__CUDABE__)
0174 #else
0175     static constexpr const char* NAME = "sphoton" ;
0176 #endif
0177
0178
0179     float3 pos ;        // 0th quad
0180     float  time ;
0181
0182     float3 mom ;        // 1st quad
0183     unsigned hitcount_iindex ;   // hi16:hitcount, lo16:iindex
0184
0185     float3 pol ;        // 2nd quad
0186     float  wavelength ;
0187
0188                         // 3rd quad
0189     unsigned orient_boundary_flag ;  // hi(1,15):(orient,boundary), lo16:flag
0190     unsigned identity ;              // hi8  extend range of index, lo24:identity
0191     unsigned index ;                 // lower 32 bits of the full index (extended into identity)
0192     unsigned flagmask ;              // full32 for flagmask
0193
0194     SPHOTON_METHOD void set_prd( unsigned  boundary, unsigned  identity, float  orient, unsigned iindex );
0195
0196     // hitcount starts as 1, hmm when to do that ?
0197     SPHOTON_METHOD void set_hitcount_one_iindex( unsigned ii ){ hitcount_iindex = (                   0x00010000u ) | (( 0x0000ffffu & ii ) << 0 ); }
0198     SPHOTON_METHOD void set_iindex__( unsigned ii ){            hitcount_iindex = ( hitcount_iindex & 0xffff0000u ) | (( 0x0000ffffu & ii ) << 0 ); }
0199     SPHOTON_METHOD void set_hitcount( unsigned hc ){            hitcount_iindex = ( hitcount_iindex & 0x0000ffffu ) | (( 0x0000ffffu & hc ) << 16); }
0200
0201     SPHOTON_METHOD unsigned iindex()   const {                 return ( hitcount_iindex & 0x0000ffffu ) >> 0  ; }
0202     SPHOTON_METHOD unsigned hitcount() const {                 return ( hitcount_iindex & 0xffff0000u ) >> 16 ; }
0203
0204     SPHOTON_METHOD unsigned flag() const {     return (orient_boundary_flag & 0x0000ffffu) >>  0 ; } // flag___     = lambda p:(p.view(np.uint32)[...,3,0] & 0xffff)
0205     SPHOTON_METHOD unsigned boundary() const { return (orient_boundary_flag & 0x7fff0000u) >> 16 ; } // boundary___ = lambda p:(p.view(np.uint32)[...,3,0] & 0x7fff0000) >> 16
0206     SPHOTON_METHOD float    orient() const {   return (orient_boundary_flag & 0x80000000u) ? -1.f : 1.f ; }
0207
0208     // bit manipulation approach of setters: clear bit(s) for flag/boundary/orient and then set them
0209     SPHOTON_METHOD void     set_flag(unsigned flag) {         orient_boundary_flag = ( orient_boundary_flag & 0xffff0000u ) | (( flag     & 0xffffu         ) <<  0 ) ; flagmask |= flag ; }
0210     SPHOTON_METHOD void     set_boundary(unsigned boundary) { orient_boundary_flag = ( orient_boundary_flag & 0x8000ffffu ) | (( boundary & 0x7fffu         ) << 16 ) ; }
0211     SPHOTON_METHOD void     set_orient(float orient){         orient_boundary_flag = ( orient_boundary_flag & 0x7fffffffu ) | (( orient < 0.f ? 0x1u : 0x0u ) << 31 ) ; }
0212
0213
0214     SPHOTON_METHOD void     set_flagmask(unsigned _flagmask) { flagmask = _flagmask ; }
0215     SPHOTON_METHOD void     addto_flagmask(unsigned flag) { flagmask |=  flag ; }
0216     SPHOTON_METHOD void     scrub_flagmask(unsigned flag) { flagmask &= ~flag ; }
0217
0218     SPHOTON_METHOD void zero_flags() { orient_boundary_flag = 0u ; identity = 0u ; index = 0u ; flagmask = 0u ; hitcount_iindex = 0u ; }
0219
0220     SPHOTON_METHOD float* data() {               return &pos.x ; }
0221     SPHOTON_METHOD const float* cdata() const {  return &pos.x ; }
0222
0223
0224     SPHOTON_METHOD void set_index(uint64_t full_index)
0225     {
0226         // lower 32 bits of full_index into index
0227         index    = ( static_cast<unsigned>((full_index >> 0 ) & 0xffffffffu) <<  0 ) ;
0228
0229         // upper 8 bits of full_index into upper 8 bits of identity with lower 24 bits of identity preserved
0230         identity = ( static_cast<unsigned>((full_index >> 32) & 0xffu      ) << 24 ) | ( identity & 0xffffffu ) ;
0231     }
0232     SPHOTON_METHOD uint64_t get_index() const {  return (static_cast<uint64_t>(identity >> 24) << 32) | static_cast<uint64_t>(index); }
0233
0234
0235     // upper 8 bits of identity field used to extend range of index beyond 4.29 billion
0236     SPHOTON_METHOD void     set_identity(unsigned id) { identity = ( identity & 0xff000000u ) | ( id & 0x00ffffffu ); }
0237     SPHOTON_METHOD unsigned get_identity() const {            return identity & 0x00ffffffu ; }
0238     SPHOTON_METHOD unsigned pmtid() const
0239     {
0240         unsigned id = get_identity() ;
0241         return id == 0 ? 0xffffffffu : id - 1 ;    // id:0 means not-a-sensor
0242     }
0243
0244
0245     SPHOTON_METHOD void zero()
0246     {
0247        pos.x = 0.f ; pos.y = 0.f ; pos.z = 0.f ; time = 0.f ;
0248        mom.x = 0.f ; mom.y = 0.f ; mom.z = 0.f ; hitcount_iindex = 0u ;
0249        pol.x = 0.f ; pol.y = 0.f ; pol.z = 0.f ; wavelength = 0.f ;
0250        orient_boundary_flag = 0u ; identity = 0u ; index = 0u ; flagmask = 0u ;
0251     }
0252
0253 #ifdef WITH_LOCALIZE
0254     template<typename T>
0255     SPHOTON_METHOD void localize(sphoton& l, const glm::tmat4x4<T>& tr, bool normalize) const
0256     {
0257         // promote float->T for transform then T->float for storage
0258         glm::tvec4<T> p_pos(pos.x, pos.y, pos.z, T(1.f));
0259         glm::tvec4<T> l_pos = tr * p_pos;
0260         l.pos = make_float3(l_pos.x, l_pos.y, l_pos.z);
0261
0262         glm::tvec4<T> p_mom(mom.x, mom.y, mom.z, T(0.f));
0263         glm::tvec4<T> l_mom_ = tr * p_mom;
0264         glm::tvec3<T> l_mom(l_mom_);
0265         if(normalize) l_mom = glm::normalize(l_mom);
0266         l.mom = make_float3(l_mom.x, l_mom.y, l_mom.z);
0267
0268         glm::tvec4<T> p_pol(pol.x, pol.y, pol.z, T(0.f));
0269         glm::tvec4<T> l_pol_ = tr * p_pol;
0270         glm::tvec3<T> l_pol(l_pol_);
0271         if(normalize) l_pol = glm::normalize(l_pol);
0272         l.pol = make_float3(l_pol.x, l_pol.y, l_pol.z);
0273     }
0274 #endif
0275
0276
0277     struct key_functor
0278     {
0279         float    tw;
0280         SPHOTON_METHOD uint64_t operator()(const sphoton& p) const
0281         {
0282             unsigned id = p.get_identity() ;
0283             unsigned bucket = static_cast<unsigned>(p.time / tw);
0284             return (uint64_t(id) << 48) | uint64_t(bucket);
0285         }
0286     };
0287
0288     struct reduce_op
0289     {
0290         SPHOTON_METHOD sphoton operator()(const sphoton& a, const sphoton& b) const
0291         {
0292             bool a_first = fminf(a.time, b.time) == a.time ;
0293             sphoton r   = a_first ? a : b ;
0294             r.flagmask  = a.flagmask | b.flagmask ;     // combined flagmask OR first flagmask ?
0295             r.set_hitcount( a.hitcount() + b.hitcount() );
0296             return r;
0297         }
0298     };
0299
0300     struct select_pred
0301     {
0302         unsigned mask;
0303         SPHOTON_METHOD bool operator()(const sphoton& p) const { return (p.flagmask & mask) != 0; }
0304     };
0305
0306
0307
0308
0309 #if defined(__CUDACC__) || defined(__CUDABE__)
0310 #else
0311
0312     SPHOTON_METHOD bool is_cerenkov()  const { return (flagmask & CERENKOV) != 0 ; }
0313     SPHOTON_METHOD bool is_reemit() const {    return (flagmask & BULK_REEMIT) != 0 ; }
0314
0315     SPHOTON_METHOD unsigned flagmask_count() const ;
0316     SPHOTON_METHOD std::string desc() const ;
0317     SPHOTON_METHOD std::string descBase() const ;
0318     SPHOTON_METHOD std::string descPos() const ;
0319     SPHOTON_METHOD std::string descDir() const ;
0320     SPHOTON_METHOD std::string descDetail() const ;
0321     SPHOTON_METHOD std::string descFlag() const ;
0322     SPHOTON_METHOD std::string descDigest() const ;
0323
0324     SPHOTON_METHOD const char* flag_() const ;
0325     SPHOTON_METHOD const char* abbrev_() const ;
0326     SPHOTON_METHOD std::string flagmask_() const ;
0327
0328
0329     SPHOTON_METHOD void ephoton() ;
0330     SPHOTON_METHOD void normalize_mom_pol();
0331     SPHOTON_METHOD void transverse_mom_pol();
0332     SPHOTON_METHOD void set_polarization(float frac_twopi);
0333     SPHOTON_METHOD static sphoton make_ephoton();
0334     SPHOTON_METHOD static NP* make_ephoton_array(size_t num_photon);
0335     SPHOTON_METHOD static NP* zeros(size_t num_photon) ;
0336     SPHOTON_METHOD static NP* demoarray(size_t num_photon);
0337     SPHOTON_METHOD std::string digest(unsigned numval=16) const  ;
0338     SPHOTON_METHOD static bool digest_match( const sphoton& a, const sphoton& b, unsigned numval=16 ) ;
0339
0340     SPHOTON_METHOD static float  DeltaMax( const float* aa, const float* bb, int num );
0341     SPHOTON_METHOD static float4 DeltaMax( const sphoton& a, const sphoton& b ) ;
0342     SPHOTON_METHOD static bool EqualFlags( const sphoton& a, const sphoton& b) ;
0343
0344     SPHOTON_METHOD static void Get( sphoton& p, const NP* a, unsigned idx );
0345     SPHOTON_METHOD static void Get( std::vector<sphoton>& pp, const NP* a );
0346
0347     SPHOTON_METHOD static void MinMaxPost( float* mn, float* mx, const NP* a, bool skip_flagmask_zero );
0348     SPHOTON_METHOD static std::string DescMinMaxPost( const NP* _a, bool skip_flagmask_zero );
0349     SPHOTON_METHOD static std::string Desc(const NP* a, size_t edge=20);
0350     SPHOTON_METHOD static NP* MockupForMergeTest(size_t ni);
0351
0352     template<typename T>
0353     SPHOTON_METHOD static bool IsPhotonArray( const NP* a );
0354
0355     SPHOTON_METHOD static void ChangeTimeInsitu( NP* a, float t0 );
0356
0357     SPHOTON_METHOD void transform_float( const glm::tmat4x4<float>&  tr, bool normalize=true );  // widens transform and uses below
0358     SPHOTON_METHOD void transform(       const glm::tmat4x4<double>& tr, bool normalize=true );
0359
0360     SPHOTON_METHOD void populate_simtrace_input( quad4& q ) const ;
0361
0362     SPHOTON_METHOD float get_cost() const ;
0363     SPHOTON_METHOD float get_fphi() const ;
0364     SPHOTON_METHOD void set_lpos(float lposcost, float lposfphi); // FOR TEMPLATE COMPATIBILITY WITH sphotonlite used only in tests
0365
0366     SPHOTON_METHOD static size_t GetIndexBeyondCursorContiguousIIndex( const sphoton* pp, size_t num, size_t cursor=0 );
0367     SPHOTON_METHOD static void   GetContiguousIIndexStartIndices( std::vector<size_t>& starts, const sphoton* pp, size_t num );
0368 #endif
0369
0370 };
0371
0372
0373
0374
0375
0376 #if defined(__CUDACC__) || defined(__CUDABE__)
0377 #else
0378
0379 struct sphotond
0380 {
0381     double3 pos ;
0382     double  time ;
0383
0384     double3 mom ;
0385     unsigned long long hitcount_iindex ;
0386
0387     double3 pol ;
0388     double  wavelength ;
0389
0390     unsigned long long orient_boundary_flag ;
0391     unsigned long long identity ;
0392     unsigned long long index ;
0393     unsigned long long flagmask ;
0394
0395
0396    SPHOTON_METHOD static void FromFloat( sphotond& d, const sphoton& s );
0397    SPHOTON_METHOD static void Get( sphotond& p, const NP* a, unsigned idx );
0398    SPHOTON_METHOD void transform_float( const glm::tmat4x4<float>& tr,  bool normalize=true );
0399    SPHOTON_METHOD void transform(       const glm::tmat4x4<double>& tr, bool normalize=true );
0400    SPHOTON_METHOD const double* cdata() const {  return &pos.x ; }
0401
0402
0403 };
0404
0405
0406 SPHOTON_METHOD void sphotond::FromFloat( sphotond& d, const sphoton& s )
0407 {
0408     typedef unsigned long long ull ;
0409
0410     d.pos.x = double(s.pos.x) ;
0411     d.pos.y = double(s.pos.y) ;
0412     d.pos.z = double(s.pos.z) ;
0413     d.time  = double(s.time) ;
0414
0415     d.mom.x = double(s.mom.x) ;
0416     d.mom.y = double(s.mom.y) ;
0417     d.mom.z = double(s.mom.z) ;
0418     d.hitcount_iindex = ull(s.hitcount_iindex) ;
0419
0420     d.pol.x = double(s.pol.x) ;
0421     d.pol.y = double(s.pol.y) ;
0422     d.pol.z = double(s.pol.z) ;
0423     d.wavelength  = double(s.wavelength) ;
0424
0425     d.orient_boundary_flag = ull(s.orient_boundary_flag) ;
0426     d.identity      = ull(s.identity) ;
0427     d.index         = ull(s.index) ;
0428     d.flagmask      = ull(s.flagmask) ;
0429 }
0430
0431 SPHOTON_METHOD void sphotond::Get( sphotond& p, const NP* a, unsigned idx )
0432 {
0433     assert(a && a->has_shape(-1,4,4) && a->ebyte == sizeof(double) && idx < unsigned(a->shape[0]) );
0434     assert( sizeof(sphotond) == sizeof(double)*16 );
0435     memcpy( &p, a->cvalues<double>() + idx*16, sizeof(sphotond) );
0436 }
0437
0438 SPHOTON_METHOD void sphotond::transform_float( const glm::tmat4x4<float>& tr, bool normalize )
0439 {
0440     glm::tmat4x4<double> trd ;
0441     TranConvert(trd, tr);
0442     transform(trd, normalize );
0443 }
0444
0445 SPHOTON_METHOD void sphotond::transform( const glm::tmat4x4<double>& tr, bool normalize )
0446 {
0447     double one(1.);
0448     double zero(0.);
0449
0450     unsigned count = 1 ;
0451     unsigned stride = 4*4 ; // effectively not used as count is 1
0452
0453     assert( sizeof(*this) == sizeof(double)*16 );
0454     double* p0 = (double*)this ;
0455
0456     Tran<double>::Apply( tr, p0, one,  count, stride, 0, false );      // transform pos as position
0457     Tran<double>::Apply( tr, p0, zero, count, stride, 4, normalize );  // transform mom as direction
0458     Tran<double>::Apply( tr, p0, zero, count, stride, 8, normalize );  // transform pol as direction
0459 }
0460
0461 #endif
0462
0463
0464
0465
0466 /**
0467 sphoton::set_prd
0468 -----------------
0469
0470 This is canonically invoked GPU side by qsim::propagate
0471 copying geometry info from the quad2 PRD struct into the sphoton.
0472
0473 TODO: relocate identity - 1 offsetting into here as this
0474 marks the transition from geometry to event information
0475 and would allow the offsetting to be better hidden.
0476
0477 See ~/opticks/notes/issues/sensor_identifier_offset_by_one_wrinkle.rst
0478
0479 **/
0480
0481
0482 SPHOTON_METHOD void sphoton::set_prd( unsigned  boundary_, unsigned  identity_, float  orient_, unsigned iindex_ )
0483 {
0484     set_boundary(boundary_);
0485     set_identity(identity_);  // formerly did identity = identity_ which scrubs hi index bits for huge simulations > 4.29 billion
0486     set_orient( orient_ );
0487     set_hitcount_one_iindex( iindex_ );
0488 }
0489
0490
0491
0492 #if defined(__CUDACC__) || defined(__CUDABE__)
0493 #else
0494
0495 #include <cassert>
0496 #include <bitset>
0497 #include "sdigest.h"
0498 #include "OpticksPhoton.hh"
0499 #include "NP.hh"
0500
0501
0502
0503 SPHOTON_METHOD unsigned sphoton::flagmask_count() const
0504 {
0505     return std::bitset<32>(flagmask).count() ;   // NB counting bits, not nibbles with bits
0506 }
0507
0508 SPHOTON_METHOD std::string sphoton::desc() const
0509 {
0510     std::stringstream ss ;
0511     ss << descBase() << " "<< descDetail() ;
0512     std::string s = ss.str();
0513     return s ;
0514 }
0515
0516 SPHOTON_METHOD std::string sphoton::descPos() const
0517 {
0518     std::stringstream ss ;
0519     ss
0520         << " pos " << pos
0521         << " t  " << std::setw(8) << time
0522         << " "
0523         ;
0524     std::string s = ss.str();
0525     return s ;
0526 }
0527
0528 SPHOTON_METHOD std::string sphoton::descDir() const
0529 {
0530     std::stringstream ss ;
0531     ss
0532         << " mom " << mom
0533         << " iindex " << std::setw(4) << iindex()
0534         << " "
0535         << " pol " << pol
0536         << " wl " << std::setw(8) << wavelength
0537         << " "
0538         ;
0539     std::string s = ss.str();
0540     return s ;
0541 }
0542
0543 SPHOTON_METHOD std::string sphoton::descBase() const
0544 {
0545     std::stringstream ss ;
0546     ss  << descPos()
0547         << descDir()
0548         ;
0549     std::string s = ss.str();
0550     return s ;
0551 }
0552 SPHOTON_METHOD std::string sphoton::descDetail() const
0553 {
0554     std::stringstream ss ;
0555     ss
0556         << " bn " << boundary()
0557         << " fl " << std::hex << flag() << std::dec
0558         << " id " << identity
0559         << " or " << orient()
0560         << " ix " << index
0561         << " fm " << std::hex << flagmask  << std::dec
0562         << " ab " << OpticksPhoton::Abbrev( flag() )
0563         << " ii " << iindex()
0564         ;
0565     std::string s = ss.str();
0566     return s ;
0567 }
0568
0569 SPHOTON_METHOD std::string sphoton::descDigest() const
0570 {
0571     std::stringstream ss ;
0572     ss
0573         << " "
0574         << " digest(16) " << digest(16)
0575         << " "
0576         << " digest(12) " << digest(12)
0577         ;
0578     std::string s = ss.str();
0579     return s ;
0580 }
0581
0582
0583 SPHOTON_METHOD const char* sphoton::flag_() const
0584 {
0585     return OpticksPhoton::Flag(flag()) ;
0586 }
0587
0588 SPHOTON_METHOD const char* sphoton::abbrev_() const
0589 {
0590     return OpticksPhoton::Abbrev(flag()) ;
0591 }
0592 SPHOTON_METHOD std::string sphoton::flagmask_() const
0593 {
0594     return OpticksPhoton::FlagMaskLabel(flagmask) ;
0595 }
0596
0597
0598 SPHOTON_METHOD std::string sphoton::descFlag() const
0599 {
0600     std::stringstream ss ;
0601     ss
0602        << " sphoton idx " << index
0603        << " flag " << flag_()
0604        << " flagmask " << flagmask_()
0605        ;
0606     std::string s = ss.str();
0607     return s ;
0608 }
0609
0610
0611 /**
0612 sphoton::ephoton
0613 ------------------
0614
0615 *ephoton* is used from qudarap/tests/QSimTest generate tests as the initial photon,
0616 which gets persisted to p0.npy
0617 The script qudarap/tests/ephoton.sh sets the envvars defining the photon
0618 depending on the TEST envvar.
0619
0620 **/
0621
0622 SPHOTON_METHOD void sphoton::ephoton()
0623 {
0624     quad4& q = (quad4&)(*this) ;
0625
0626     qvals( q.q0.f ,        "EPHOTON_POST" , "0,0,0,0" );                      // position, time
0627     qvals( q.q1.f, q.q2.f, "EPHOTON_MOMW_POLW", "1,0,0,1,0,1,0,500" );  // direction, weight,  polarization, wavelength
0628     qvals( q.q3.i ,        "EPHOTON_FLAG", "0,0,0,0" );
0629     normalize_mom_pol();
0630     transverse_mom_pol();
0631 }
0632
0633 SPHOTON_METHOD void sphoton::normalize_mom_pol()
0634 {
0635     mom = normalize(mom);
0636     pol = normalize(pol);
0637 }
0638
0639 SPHOTON_METHOD void sphoton::transverse_mom_pol()
0640 {
0641     float mom_pol = fabsf( dot(mom, pol)) ;
0642     float eps = 1e-5 ;
0643     bool is_transverse = mom_pol < eps ;
0644
0645     if(!is_transverse )
0646     {
0647         std::cout
0648              << " sphoton::transverse_mom_pol "
0649              << " FAIL "
0650              << " mom " << mom
0651              << " pol " << pol
0652              << " mom_pol " << mom_pol
0653              << " eps " << eps
0654              << " is_transverse " << is_transverse
0655              << std::endl
0656              ;
0657     }
0658     assert(is_transverse);
0659 }
0660
0661
0662 /**
0663 sphoton::set_polarization
0664 ---------------------------
0665
0666 Start with a vector in XY plane and rotate that using smath::rotateUz
0667 to be transverse to mom.
0668
0669 **/
0670
0671 SPHOTON_METHOD void sphoton::set_polarization(float frac_twopi)
0672 {
0673     float phase = 2.f*M_PIf*frac_twopi ;
0674     pol.x = cosf(phase) ;
0675     pol.y = sinf(phase) ;
0676     pol.z = 0.f ;
0677     smath::rotateUz(pol, mom); // rotate pol to be transverse to mom
0678 }
0679
0680
0681 SPHOTON_METHOD sphoton sphoton::make_ephoton()  // static
0682 {
0683     sphoton p ;
0684     p.ephoton();
0685     return p ;
0686 }
0687
0688 /**
0689 sphoton::make_ephoton_array
0690 -----------------------------
0691
0692 Formerly QSim::duplicate_dbg_ephoton
0693
0694 This is called from QSimTest::mock_propagate
0695
0696 **/
0697
0698 SPHOTON_METHOD NP* sphoton::make_ephoton_array(size_t num_photon) // static
0699 {
0700     sphoton ep ;
0701     ep.ephoton() ;
0702
0703     NP* ph = zeros(num_photon);
0704     sphoton* pp  = (sphoton*)ph->bytes();
0705     for(size_t i=0 ; i < num_photon ; i++)
0706     {
0707         sphoton p = ep  ;  // start from ephoton
0708         p.pos.y = float(i)*100.f ;
0709         pp[i] = p ;
0710     }
0711     return ph ;
0712 }
0713 SPHOTON_METHOD NP* sphoton::zeros(size_t num_photon) // static
0714 {
0715     NP* ph = NP::Make<float>(num_photon, 4, 4 );
0716     return ph ;
0717 }
0718
0719 SPHOTON_METHOD NP* sphoton::demoarray(size_t num_photon) // static
0720 {
0721     NP* ph = zeros(num_photon);
0722     sphoton* pp  = (sphoton*)ph->bytes();
0723     for(size_t i=0 ; i < num_photon ; i++)
0724     {
0725         sphoton& p = pp[i] ;
0726
0727         p.pos = make_float3( 0.1f + i*0.1f, 0.2f + i*0.1f , 0.3f + i*0.1f );
0728         p.time = i*0.1f ;
0729
0730         p.mom = make_float3( 0.1f + i*0.1f, 0.2f + i*0.1f , 0.3f + i*0.1f );
0731         p.set_orient( i % 2 == 0 ? 1.f : -1.f );
0732         p.set_hitcount_one_iindex( i*1000 );
0733
0734         p.pol = make_float3( 0.1f + i*0.1f, 0.2f + i*0.1f , 0.3f + i*0.1f );
0735         p.wavelength = 400.f + 10.f*i ;
0736
0737         p.set_boundary( i*10 );
0738         p.set_flag( i*100 );
0739
0740         uint64_t full_index = 0xffffffff + i ;
0741         p.set_index( full_index );
0742
0743         p.set_identity(  i*200 );
0744         p.set_flagmask( SURFACE_DETECT | EFFICIENCY_COLLECT );
0745
0746     }
0747     return ph ;
0748 }
0749
0750
0751
0752
0753
0754
0755
0756 SPHOTON_METHOD std::string sphoton::digest(unsigned numval) const
0757 {
0758     assert( numval <= 16 );
0759     return sdigest::Buf( (const char*)cdata() , numval*sizeof(float) );
0760 }
0761
0762 SPHOTON_METHOD bool sphoton::digest_match( const sphoton& a, const sphoton& b, unsigned numval )  // static
0763 {
0764     std::string adig = a.digest(numval);
0765     std::string bdig = b.digest(numval);
0766     return strcmp( adig.c_str(), bdig.c_str() ) == 0 ;
0767 }
0768
0769 SPHOTON_METHOD float sphoton::DeltaMax( const float* aa, const float* bb, int num )  // static
0770 {
0771     float dmax = 0.f ;
0772     for(int i=0 ; i < num ; i++)
0773     {
0774         float a = *(aa+i) ;
0775         float b = *(bb+i) ;
0776         float ab = std::abs(a - b);
0777         if(ab > dmax) dmax = ab ;
0778     }
0779     return dmax ;
0780 }
0781
0782 SPHOTON_METHOD float4 sphoton::DeltaMax( const sphoton& a, const sphoton& b )  // static
0783 {
0784     float4 dmax = make_float4( 0.f, 0.f, 0.f, 0.f );
0785     dmax.x = DeltaMax( &a.pos.x, &b.pos.x, 3 );
0786     dmax.y = DeltaMax( &a.mom.x, &b.mom.x, 3 );
0787     dmax.z = DeltaMax( &a.pol.x, &b.pol.x, 3 );
0788     dmax.w = std::max( DeltaMax( &a.time,  &b.time, 1 ), DeltaMax( &a.wavelength,  &b.wavelength, 1 ) ) ;
0789     return dmax ;
0790 }
0791
0792
0793 SPHOTON_METHOD bool sphoton::EqualFlags( const sphoton& a, const sphoton& b) // static
0794 {
0795     return a.hitcount_iindex == b.hitcount_iindex &&
0796            a.orient_boundary_flag == b.orient_boundary_flag &&
0797            a.identity == b.identity &&
0798            a.index == b.index &&
0799            a.flagmask == b.flagmask
0800            ;
0801
0802 }
0803
0804
0805
0806 SPHOTON_METHOD void sphoton::Get( sphoton& p, const NP* a, unsigned idx )
0807 {
0808     bool expected = a && a->has_shape(-1,4,4) && a->ebyte == sizeof(float) && idx < unsigned(a->shape[0]) ;
0809     if(!expected) std::cerr
0810         << "sphoton::Get not expected error "
0811         << " a " << ( a ? "Y" : "N" )
0812         << " a.shape " << ( a ? a->sstr() : "-" )
0813         << " a.ebyte " << ( a ? a->ebyte : -1 )
0814         << " a.shape[0] " << ( a ? a->shape[0] : -1 )
0815         << " idx " << idx
0816         << std::endl
0817         ;
0818
0819     assert( expected  );
0820     assert( sizeof(sphoton) == sizeof(float)*16 );
0821     memcpy( &p, a->cvalues<float>() + idx*16, sizeof(sphoton) );
0822 }
0823
0824 SPHOTON_METHOD void sphoton::Get( std::vector<sphoton>& pp, const NP* a )
0825 {
0826     bool is_f = IsPhotonArray<float>(a);
0827     assert( is_f );
0828     unsigned num = a->shape[0] ;
0829     pp.resize(num);
0830     memcpy( pp.data(), a->cvalues<float>(), sizeof(sphoton)*num );
0831 }
0832
0833 /**
0834 sphoton::MinMaxPost
0835 --------------------
0836
0837 Find minimum and maximum position and time over all
0838 filled sphoton entries in the array.  This works with
0839 record arrays by temporarily reshaping to be photon array
0840 shaped and then reshaping back again.
0841
0842 **/
0843
0844 SPHOTON_METHOD void sphoton::MinMaxPost( float* mn, float* mx, const NP* _a, bool skip_flagmask_zero  )
0845 {
0846     NP* a = const_cast<NP*>(_a);
0847
0848     std::vector<NP::INT> sh = a->shape ;
0849     a->change_shape(-1,4,4);
0850
0851     bool is_f = IsPhotonArray<float>(a);
0852     bool is_d = IsPhotonArray<double>(a);
0853     bool is_f_xor_d = is_f ^ is_d ;
0854
0855     if(!is_f_xor_d) std::cerr
0856         << "sphoton::MinMaxPost FATAL UNEXPECTED ARRAY "
0857         << " is_f " << ( is_f ? "YES" : "NO " )
0858         << " is_d " << ( is_d ? "YES" : "NO " )
0859         << "\n"
0860         ;
0861
0862     bool expect = is_f_xor_d ;
0863     assert(expect);
0864     if(!expect) std::raise(SIGINT);
0865
0866     int ni = a->num_items() ;
0867     int nj = 4 ;
0868     float MAX = std::numeric_limits<float>::max() ;
0869
0870     for(int j=0 ; j < nj ; j++)
0871     {
0872         mn[j] = MAX ;
0873         mx[j] = -MAX ;
0874     }
0875     // float limits are big enough as output is in float anyhow
0876
0877
0878     for(int i=0 ; i < ni ; i++)
0879     {
0880         if( is_f )
0881         {
0882             sphoton* p = reinterpret_cast<sphoton*>(a->bytes() + i*a->item_bytes());
0883             if(skip_flagmask_zero && p->flagmask == 0u) continue ;
0884             const float* xyzt = p->cdata();
0885             for(int j=0 ; j < nj ; j++)
0886             {
0887                 float vj = xyzt[j] ;
0888                 if( vj < mn[j] ) mn[j] = vj ;
0889                 if( vj > mx[j] ) mx[j] = vj ;
0890             }
0891         }
0892         else if( is_d )
0893         {
0894             sphotond* p = reinterpret_cast<sphotond*>(a->bytes() + i*a->item_bytes());
0895             if(skip_flagmask_zero && p->flagmask == 0u) continue ;
0896             const double* xyzt = p->cdata();
0897
0898             for(int j=0 ; j < nj ; j++)
0899             {
0900                 double vj = xyzt[j] ;
0901                 if( vj < double(mn[j]) ) mn[j] = vj ;
0902                 if( vj > double(mx[j]) ) mx[j] = vj ;
0903             }
0904         }
0905     }
0906     a->reshape(sh);
0907 }
0908
0909 SPHOTON_METHOD std::string sphoton::DescMinMaxPost( const NP* _a, bool skip_flagmask_zero )
0910 {
0911     float4 mn = {} ;
0912     float4 mx = {} ;
0913     MinMaxPost( &mn.x , &mx.x, _a, skip_flagmask_zero );
0914
0915     std::stringstream ss ;
0916     ss
0917         << "[sphoton::DescMinMaxPost\n"
0918         << " a " << ( _a ? _a->sstr() : "-" ) << "\n"
0919         << " mn " << mn << "\n"
0920         << " mx " << mx << "\n"
0921         << "]sphoton::DescMinMaxPost\n"
0922         ;
0923     std::string str = ss.str();
0924     return str ;
0925 }
0926
0927 SPHOTON_METHOD std::string sphoton::Desc(const NP* a, size_t edge) // static
0928 {
0929     std::stringstream ss ;
0930     ss << "[sphoton::Desc a.sstr " << ( a ? a->sstr() : "-" ) << "\n" ;
0931     size_t ni = a->num_items();
0932     sphoton* pp = (sphoton*)a->bytes();
0933     for(size_t i=0 ; i < ni ; i++)
0934     {
0935         if(i < edge || i > (ni - edge)) ss << std::setw(6) << i << " : " << pp[i].desc() << "\n" ;
0936         else if( i == edge )  ss << std::setw(6) << "" << " : " << "..." << "\n" ;
0937     }
0938
0939     std::string str = ss.str() ;
0940     return str ;
0941 }
0942
0943 /**
0944 sphoton::MockupForMergeTest
0945 ----------------------------
0946
0947 See sutil::unmerge and SPM_test::merge_with_expectation
0948 for a better way to test merging with known expected result.
0949
0950 **/
0951
0952
0953 SPHOTON_METHOD NP* sphoton::MockupForMergeTest(size_t ni) // static
0954 {
0955     NP* a = NP::Make<float>(ni, 4, 4);
0956     sphoton* aa = (sphoton*)a->bytes();
0957     for(size_t i=0 ; i < ni ; i++)
0958     {
0959         sphoton& p = aa[i];
0960         p.time = float( i % 100 )*0.1f ;
0961         p.set_index(i);
0962         p.set_identity( i % 100 );
0963         p.set_hitcount( 1 );
0964         p.flagmask = i % 5 == 0 ? EFFICIENCY_COLLECT : EFFICIENCY_CULL ;
0965     }
0966     return a ;
0967 }
0968
0969
0970
0971
0972
0973
0974
0975
0976
0977
0978
0979
0980
0981 template<typename T>
0982 SPHOTON_METHOD bool sphoton::IsPhotonArray( const NP* a )
0983 {
0984     assert( sizeof(sphoton) == sizeof(float)*16 );
0985     assert( sizeof(sphotond) == sizeof(double)*16 );
0986     return a && a->has_shape(-1,4,4) && a->ebyte == sizeof(T) ;
0987 }
0988
0989
0990 /**
0991 sphoton::ChangeTimeInsitu
0992 --------------------------
0993
0994 Casts the bytes of each item of the array to an sphoton pointer
0995 and uses that to inplace change all the times.  This lowlevel
0996 kinda thing only works because the sphoton struct is and
0997 will remain very simple.
0998
0999 **/
1000
1001 SPHOTON_METHOD void sphoton::ChangeTimeInsitu( NP* a, float t0 )
1002 {
1003     bool is_f = IsPhotonArray<float>(a);
1004     bool is_d = IsPhotonArray<double>(a);
1005     bool is_f_xor_d = is_f ^ is_d ;
1006
1007     if(!is_f_xor_d) std::cerr
1008         << "sphoton::ChangeTimeInsitu FATAL "
1009         << " is_f " << ( is_f ? "YES" : "NO " )
1010         << " is_d " << ( is_d ? "YES" : "NO " )
1011         << "\n"
1012         ;
1013
1014     bool expect = is_f_xor_d ;
1015     assert(expect);
1016     if(!expect) std::raise(SIGINT);
1017
1018     for(NP::INT i=0 ; i < a->num_items() ; i++)
1019     {
1020         if( is_f )
1021         {
1022             sphoton* p = reinterpret_cast<sphoton*>(a->bytes() + i*a->item_bytes());
1023             p->time = t0 ;
1024         }
1025         else if( is_d )
1026         {
1027             sphotond* p = reinterpret_cast<sphotond*>(a->bytes() + i*a->item_bytes());
1028             p->time = t0 ;
1029         }
1030     }
1031 }
1032
1033
1034
1035
1036 /**
1037 sphoton::transform_float
1038 --------------------------------
1039
1040 Its better to keep transforms in double and use sphoton::transform but if
1041 double precision transforms are not yet available this will widen
1042 the transform and use that, which is better than using float.
1043
1044 **/
1045
1046 SPHOTON_METHOD void sphoton::transform_float( const glm::tmat4x4<float>& tr, bool normalize )
1047 {
1048     glm::tmat4x4<double> trd ;
1049     TranConvert(trd, tr);
1050     transform(trd, normalize );
1051 }
1052
1053 /**
1054 sphoton::transform
1055 --------------------
1056
1057 Applies a double precision transform to float precision sphoton
1058 with Tran<double>::ApplyToFloat by temporarily widening the sphoton
1059 pos/mom/pol in order to do the matrix multiplication in double precision
1060 and then save back into the sphoton in float.
1061
1062 And alternative to this would be to widen the sphoton to sphotond
1063 and then use that.
1064
1065 **/
1066
1067 SPHOTON_METHOD void sphoton::transform( const glm::tmat4x4<double>& tr, bool normalize )
1068 {
1069     float one(1.);
1070     float zero(0.);
1071
1072     unsigned count = 1 ;
1073     unsigned stride = 4*4 ; // effectively not used as count is 1
1074
1075     assert( sizeof(*this) == sizeof(float)*16 );
1076     float* p0 = (float*)this ;
1077
1078     Tran<double>::ApplyToFloat( tr, p0, one,  count, stride, 0, false );      // transform pos as position
1079     Tran<double>::ApplyToFloat( tr, p0, zero, count, stride, 4, normalize );  // transform mom as direction
1080     Tran<double>::ApplyToFloat( tr, p0, zero, count, stride, 8, normalize );  // transform pol as direction
1081 }
1082
1083
1084
1085 /**
1086 sphoton::populate_simtrace_input
1087 ---------------------------------
1088
1089 **/
1090
1091
1092 SPHOTON_METHOD void sphoton::populate_simtrace_input( quad4& q ) const
1093 {
1094     q.q0.f.x = pos.x ;
1095     q.q0.f.y = pos.y ;
1096     q.q0.f.z = pos.z ;
1097
1098     q.q1.f.x = mom.x ;
1099     q.q1.f.y = mom.y ;
1100     q.q1.f.z = mom.z ;
1101 }
1102
1103
1104
1105
1106 SPHOTON_METHOD float sphoton::get_cost() const
1107 {
1108     return normalize_cost(pos);  // scuda.h
1109 }
1110 SPHOTON_METHOD float sphoton::get_fphi() const
1111 {
1112     return normalize_fphi(pos);  // scuda.h
1113 }
1114 /**
1115 sphoton::set_lpos USED ONLY IN TESTS [ASSUMES UNIT RADIUS]
1116 -----------------------------------------------------------
1117
1118 ONLY FOR TEMPLATE COMPATIBILITY WITH sphotonlite::set_lpos
1119
1120 pos [-1,0,0] is problematic
1121
1122 **/
1123
1124 SPHOTON_METHOD void sphoton::set_lpos(float lposcost, float lposfphi)
1125 {
1126     float cost = lposcost < -1.f ? -1.f : (lposcost > 1.f ? 1.f : lposcost);
1127     float fphi = lposfphi < 0.f ? 0.f : (lposfphi >= 1.f ? 0.999999f : lposfphi);  // avoid exactly 1.0
1128
1129     float theta = acosf(cost);
1130
1131     float phi = (fphi * 2.0f - 1.0f) * M_PIf; // Recover original signed phi in [-π, π)
1132
1133     float sin_theta = sinf(theta);
1134
1135     pos.x = sin_theta * cosf(phi);
1136     pos.y = sin_theta * sinf(phi);
1137     pos.z = cost;
1138 }
1139
1140 /**
1141 sphoton::GetIndexBeyondCursorContiguousIIndex
1142 --------------------------------------------------
1143
1144 The "IIndex" is the instance-index that can be used to
1145 access a transform for the geometry.
1146
1147 1. get *ii0* iindex at starting *cursor* index
1148 2. iterate ahead thru photon array until a different iindex is found
1149    at which point the index of the new start is returned
1150 3. if no different iindex is found return num, which is one beyond
1151    the highest valid index
1152
1153 **/
1154
1155
1156 SPHOTON_METHOD size_t sphoton::GetIndexBeyondCursorContiguousIIndex( const sphoton* pp, size_t num, size_t cursor ) // static
1157 {
1158     assert( cursor < num );
1159     unsigned ii0 = pp[cursor].iindex();
1160     for(size_t i=cursor+1 ; i < num ; i++)
1161     {
1162         unsigned ii = pp[i].iindex();
1163         if( ii != ii0 ) return i ;
1164     }
1165     return num ;
1166 }
1167
1168 /**
1169 sphoton::GetContiguousIIndexStartIndices
1170 -------------------------------------------
1171
1172 Note that there is no need an equivalent sphotonlite implementation
1173 because the localization for sphotonlite is done directly.
1174
1175 **/
1176
1177 SPHOTON_METHOD void sphoton::GetContiguousIIndexStartIndices( std::vector<size_t>& starts, const sphoton* pp, size_t num )
1178 {
1179     size_t cursor = 0 ;
1180     starts.push_back(cursor);
1181     do
1182     {
1183         cursor = sphoton::GetIndexBeyondCursorContiguousIIndex(pp, num, cursor );
1184         starts.push_back(cursor);
1185     }
1186     while( cursor < num );
1187 }
1188
1189
1190
1191
1192
1193 #endif
1194
1195
1196 #if defined(__CUDACC__) || defined(__CUDABE__)
1197 #else
1198 inline std::ostream& operator<<(std::ostream& os, const sphoton& p )
1199 {
1200     os << "sphoton"
1201        << std::endl
1202        << p.desc()
1203        << std::endl
1204        ;
1205     return os;
1206 }
1207 #endif
1208
1209
1210 /*
1211 union qphoton
1212 {
1213     quad4   q ;
1214     sphoton p ;
1215 };
1216 */
1217
1218
1219 struct sphoton_selector
1220 {
1221     uint32_t hitmask ;
1222     sphoton_selector(uint32_t hitmask_) : hitmask(hitmask_) {};
1223     SPHOTON_METHOD bool operator() (const sphoton& p) const { return ( p.flagmask  & hitmask ) == hitmask  ; }   // require all bits of the mask to be set
1224     SPHOTON_METHOD bool operator() (const sphoton* p) const { return ( p->flagmask & hitmask ) == hitmask  ; }   // require all bits of the mask to be set
1225 };
1226
1227