eic-opticks/sysrap/SEvent.cc

0001 #include <csignal>
0002 #include "NP.hh"
0003
0004 #include "scuda.h"
0005 #include "squad.h"
0006 #include "sqat4.h"
0007 #include "stran.h"
0008 #include "sframe.h"
0009
0010 #include "srngcpu.h"
0011 using RNG = srngcpu ;
0012
0013 #include "storch.h"
0014 #include "scerenkov.h"
0015 #include "sscint.h"
0016 #include "scarrier.h"
0017
0018 #include "sxyz.h"
0019 #include "ssys.h"
0020 #include "sstr.h"
0021
0022 #include "NPX.h"
0023 #include "SLOG.hh"
0024
0025 #include "OpticksGenstep.h"
0026
0027 #include "SEvt.hh"
0028 #include "SGenstep.h"
0029 #include "SEvent.hh"
0030
0031 const plog::Severity SEvent::LEVEL = SLOG::EnvLevel("SEvent", "DEBUG") ;
0032
0033
0034
0035
0036 NP* SEvent::GENSTEP = nullptr ;
0037 NP* SEvent::GetGENSTEP()
0038 {
0039     LOG_IF(info, SEvt::LIFECYCLE ) <<  " GENSTEP " << ( GENSTEP ? GENSTEP->sstr() : "-" );
0040     return GENSTEP ;
0041 }
0042
0043 void SEvent::SetGENSTEP(NP* gs)
0044 {
0045     GENSTEP = gs ;
0046     LOG_IF(info, SEvt::LIFECYCLE ) <<  " GENSTEP " << ( GENSTEP ? GENSTEP->sstr() : "-" );
0047 }
0048
0049 bool SEvent::HasGENSTEP()
0050 {
0051     return GENSTEP != nullptr ;
0052 }
0053
0054 NP* SEvent::HIT = nullptr ;
0055 NP* SEvent::GetHIT()
0056 {
0057     return HIT ;
0058 }
0059 void SEvent::SetHIT(NP* gs)
0060 {
0061     HIT = gs ;
0062 }
0063 bool SEvent::HasHIT()
0064 {
0065     return HIT != nullptr ;
0066 }
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077 NP* SEvent::MakeDemoGenstep(const char* config, int idx)
0078 {
0079     NP* gs = nullptr ;
0080     if(     sstr::StartsWith(config, "count")) gs = MakeCountGenstep(config) ;
0081     else if(sstr::StartsWith(config, "torch")) gs = MakeTorchGenstep(idx) ;
0082     assert(gs);
0083     LOG(LEVEL)
0084        << " config " << ( config ? config :  "-" )
0085        << " gs " << ( gs ? gs->sstr() : "-" )
0086        ;
0087
0088     return gs ;
0089 }
0090
0091
0092
0093
0094 #ifdef WITH_OLD_FRAME
0095 NP* SEvent::MakeInputPhotonGenstep(const NP* input_photon, int gentype, const sframe* fr )
0096 {
0097     assert( gentype == OpticksGenstep_INPUT_PHOTON || gentype == OpticksGenstep_INPUT_PHOTON_SIMTRACE );
0098     std::vector<quad6> qgs(1) ;
0099     qgs[0].zero() ;
0100     qgs[0] = MakeInputPhotonGenstep_(input_photon, gentype, fr );
0101     NP* ipgs = NPX::ArrayFromVec<float,quad6>( qgs, 6, 4) ;
0102     return ipgs ;
0103 }
0104 #else
0105 NP* SEvent::MakeInputPhotonGenstep(const NP* input_photon, int gentype, const sfr* fr )
0106 {
0107     assert( gentype == OpticksGenstep_INPUT_PHOTON || gentype == OpticksGenstep_INPUT_PHOTON_SIMTRACE );
0108     std::vector<quad6> qgs(1) ;
0109     qgs[0].zero() ;
0110     qgs[0] = MakeInputPhotonGenstep_(input_photon, gentype, fr );
0111     NP* ipgs = NPX::ArrayFromVec<float,quad6>( qgs, 6, 4) ;
0112     return ipgs ;
0113 }
0114 #endif
0115
0116
0117
0118 /**
0119 SEvent::MakeInputPhotonGenstep_
0120 ---------------------------------
0121
0122 Now called from SEvt::addInputGenstep (formerly from SEvt::addFrameGenstep and before that SEvt::setFrame)
0123 Note that the only thing taken from the *input_photon* is the
0124 number of photons so this can work with either local or
0125 transformed *input_photon*.
0126
0127 The m2w transform from the frame is copied into the genstep.
0128 HMM: is that actually used ? Because the frame is also persisted.
0129
0130 **/
0131
0132 #ifdef WITH_OLD_FRAME
0133 quad6 SEvent::MakeInputPhotonGenstep_(const NP* input_photon, int gentype, const sframe* fr )
0134 {
0135     LOG(LEVEL) << " input_photon " << NP::Brief(input_photon) ;
0136
0137     quad6 ipgs ;
0138     ipgs.zero();
0139     ipgs.set_gentype( gentype );
0140     ipgs.set_numphoton(  input_photon->shape[0]  );
0141
0142     if(fr) fr->m2w.write(ipgs); // copy fr->m2w into ipgs.q2,q3,q4,q5
0143
0144     return ipgs ;
0145 }
0146 #else
0147 quad6 SEvent::MakeInputPhotonGenstep_(const NP* input_photon, int gentype, const sfr* fr )
0148 {
0149     LOG(LEVEL) << " input_photon " << NP::Brief(input_photon) ;
0150
0151     quad6 ipgs ;
0152     ipgs.zero();
0153     ipgs.set_gentype( gentype );
0154     ipgs.set_numphoton(  input_photon->shape[0]  );
0155
0156     assert(fr);
0157     bool skip_col3 = true ;
0158     ipgs.read_transform( glm::value_ptr(fr->m2w), skip_col3 );
0159
0160     return ipgs ;
0161 }
0162 #endif
0163
0164
0165
0166
0167
0168 /**
0169 SEvent::MakeTorchGenstep
0170 --------------------------
0171
0172 Canonically invoked from SEvt::AddTorchGenstep
0173 which seems to need user code invokation.
0174 HMM: perhaps SEventConfig to do this in standardized place ?
0175
0176 Default idx_arg for these is -1 indicating with_index:false
0177
0178 **/
0179
0180 NP* SEvent::MakeTorchGenstep(   int idx_arg){ return MakeGenstep( OpticksGenstep_TORCH, idx_arg ) ; }
0181 NP* SEvent::MakeCerenkovGenstep(int idx_arg){ return MakeGenstep( OpticksGenstep_CERENKOV, idx_arg ) ; }
0182 NP* SEvent::MakeScintGenstep(   int idx_arg){ return MakeGenstep( OpticksGenstep_SCINTILLATION, idx_arg ) ; }
0183 NP* SEvent::MakeCarrierGenstep( int idx_arg){ return MakeGenstep( OpticksGenstep_CARRIER, idx_arg ) ; }
0184
0185
0186 /**
0187 SEvent::MakeGenstep
0188 ---------------------
0189
0190 NB index_arg is userspace 0-based index, that is not the same as the internal SEvt::index
0191 which may be offset by OPTICKS_START_INDEX
0192
0193 **/
0194
0195
0196 NP* SEvent::MakeGenstep( int gentype, int index_arg )
0197 {
0198     bool with_index = index_arg != -1 ;
0199     if(with_index) assert( index_arg >= 0 );  // index_arg is 0-based
0200
0201
0202     int64_t num_ph = with_index ? SEventConfig::NumPhoton(index_arg)  : ssys::getenvint64spec("SEvent__MakeGenstep_num_ph", "100" ) ;
0203     int64_t num_gs = with_index ? SEventConfig::NumGenstep(index_arg) : ssys::getenvint64spec("SEvent__MakeGenstep_num_gs", "1"   ) ;
0204
0205     bool dump = ssys::getenvbool("SEvent_MakeGenstep_dump");
0206     const int M = 1000000 ;
0207
0208     LOG(LEVEL)
0209         << " gentype " << gentype
0210         << " index_arg " << index_arg
0211         << " with_index " << ( with_index ? "YES" : "NO " )
0212         << " num_ph " << num_ph
0213         << " num_ph/M " << num_ph/M
0214         << " num_gs " << num_gs
0215         << " dump " << dump
0216         ;
0217
0218     assert( num_gs > 0 );
0219
0220     NP* gs = NP::Make<float>(num_gs, 6, 4 );
0221     gs->set_meta<std::string>("creator", "SEvent::MakeGenstep" );
0222     gs->set_meta<int64_t>("num_ph", num_ph );
0223     gs->set_meta<int64_t>("num_gs", num_gs );
0224     gs->set_meta<int>("index_arg",  index_arg );
0225
0226
0227     int gs_start = 0 ;
0228     int gs_stop = num_gs ;
0229     size_t _gs_ph = num_ph/num_gs ; // divide the num_ph equally between the num_gs
0230     assert( _gs_ph <= SGenstep::MAX_SLOT_PER_SLICE );   // fits in int32_t : so can narrow
0231     int gs_ph(_gs_ph) ;
0232
0233     switch(gentype)
0234     {
0235         case  OpticksGenstep_TORCH:         FillGenstep<storch>(   gs, gs_start, gs_stop, gs_ph, dump) ; break ;
0236         case  OpticksGenstep_CERENKOV:      FillGenstep<scerenkov>(gs, gs_start, gs_stop, gs_ph, dump) ; break ;
0237         case  OpticksGenstep_SCINTILLATION: FillGenstep<sscint>(   gs, gs_start, gs_stop, gs_ph, dump) ; break ;
0238         case  OpticksGenstep_CARRIER:       FillGenstep<scarrier>( gs, gs_start, gs_stop, gs_ph, dump) ; break ;
0239     }
0240     return gs ;
0241 }
0242
0243 template<typename T>
0244 void SEvent::FillGenstep( NP* gs, int gs_start, int gs_stop, int numphoton_per_genstep, bool dump )
0245 {
0246     int num_gs = gs ? gs->shape[0] : 0 ;
0247     assert( gs_start <  num_gs );
0248     assert( gs_stop  <= num_gs );
0249
0250     T* tt = (T*)gs->bytes() ;
0251     for(int i=gs_start ; i < gs_stop ; i++ )
0252     {
0253         unsigned genstep_id = i ;
0254         T::FillGenstep( tt[i], genstep_id, numphoton_per_genstep, dump ) ;
0255     }
0256 }
0257
0258 template void SEvent::FillGenstep<storch>(    NP* gs, int gs_start, int gs_stop, int numphoton_per_genstep, bool dump );
0259 template void SEvent::FillGenstep<scerenkov>( NP* gs, int gs_start, int gs_stop, int numphoton_per_genstep, bool dump );
0260 template void SEvent::FillGenstep<sscint>(    NP* gs, int gs_start, int gs_stop, int numphoton_per_genstep, bool dump );
0261 template void SEvent::FillGenstep<scarrier>(  NP* gs, int gs_start, int gs_stop, int numphoton_per_genstep, bool dump );
0262
0263
0264 /**
0265 SEvent::MakeSeed
0266 ------------------
0267
0268 Creates seed array which provides genstep reference indices
0269 for each photon slot.
0270
0271 Normally this is done on device using involved thrust,
0272 here is a simple CPU implementation of that.
0273
0274 **/
0275
0276 NP* SEvent::MakeSeed( const NP* gs )
0277 {
0278     assert( gs->has_shape(-1,6,4) );
0279     int num_gs = gs->shape[0] ;
0280     const storch* tt = (storch*)gs->bytes() ;
0281     // storch type is used only to access numphoton in the quad6 (0,3) position
0282     // which all genstep types place at the same offset
0283     // HMM: using quad6 would be clearer
0284
0285     // vector of numphoton in each genstep
0286     std::vector<int> gsp(num_gs) ;
0287     for(int i=0 ; i < num_gs ; i++ ) gsp[i] = tt[i].numphoton ;
0288
0289     int tot_photons = 0 ;
0290     for(int i=0 ; i < num_gs ; i++ ) tot_photons += gsp[i] ;
0291
0292     NP* se = NP::Make<int>( tot_photons );
0293     int* sev = se->values<int>();
0294
0295     // duplicate genstep reference index into the seed array
0296     // which is the same length as total photons
0297     int offset = 0 ;
0298     for(int i=0 ; i < num_gs ; i++)
0299     {
0300         int np = gsp[i] ;
0301         for(int p=0 ; p < np ; p++) sev[offset+p] = i ;
0302         offset += np ;
0303     }
0304     return se ;
0305 }
0306
0307
0308
0309
0310
0311 NP* SEvent::MakeCountGenstep(const char* config, int* total ) // static
0312 {
0313     std::vector<int>* photon_counts_per_genstep = nullptr ;
0314     if( config == nullptr )
0315     {
0316         (*photon_counts_per_genstep) = { 3, 5, 2, 0, 1, 3, 4, 2, 4 };
0317     }
0318     return MakeCountGenstep(*photon_counts_per_genstep, total);
0319 }
0320
0321 /**
0322 SEvent::MakeCountGenstep
0323 ---------------------------
0324
0325 Used by qudarap/tests/QEvtTest.cc
0326
0327 HMM: MAYBE MOVE THIS TO SGenstep ?
0328
0329
0330 **/
0331
0332
0333 NP* SEvent::MakeCountGenstep(const std::vector<int>& counts, int* total ) // static
0334 {
0335     int gencode = OpticksGenstep_TORCH ;
0336     std::vector<quad6> gensteps ;
0337
0338     if(total) *total = 0 ;
0339     for(unsigned i=0 ; i < counts.size() ; i++)
0340     {
0341         quad6 gs ; gs.zero();
0342
0343         int gridaxes = XYZ ;
0344         int gsid = 0 ;
0345         int photons_per_genstep = counts[i];
0346
0347         if(total) *total += photons_per_genstep ;
0348
0349         SGenstep::ConfigureGenstep(gs, gencode, gridaxes, gsid, photons_per_genstep );
0350
0351         gs.q1.f.x = 0.f ;  gs.q1.f.y = 0.f ;  gs.q1.f.z = 0.f ;  gs.q1.f.w = 0.f ;
0352
0353         // identity transform to avoid nan
0354         gs.q2.f.x = 1.f ;  gs.q2.f.y = 0.f ;  gs.q2.f.z = 0.f ;  gs.q2.f.w = 0.f ;
0355         gs.q3.f.x = 0.f ;  gs.q3.f.y = 1.f ;  gs.q3.f.z = 0.f ;  gs.q3.f.w = 0.f ;
0356         gs.q4.f.x = 0.f ;  gs.q4.f.y = 0.f ;  gs.q4.f.z = 1.f ;  gs.q4.f.w = 0.f ;
0357         gs.q5.f.x = 0.f ;  gs.q5.f.y = 0.f ;  gs.q5.f.z = 0.f ;  gs.q5.f.w = 1.f ;
0358
0359         gensteps.push_back(gs);
0360     }
0361     return SGenstep::MakeArray(gensteps);
0362 }
0363
0364 unsigned SEvent::SumCounts(const std::vector<int>& counts) // static
0365 {
0366     unsigned total = 0 ;
0367     for(unsigned i=0 ; i < counts.size() ; i++) total += counts[i] ;
0368     return total ;
0369 }
0370
0371
0372 /**
0373 SEvent::ExpectedSeeds
0374 ----------------------
0375
0376 From a vector of counts populate the vector of seeds by simple CPU side duplication.
0377
0378 **/
0379
0380 void SEvent::ExpectedSeeds(std::vector<int>& seeds, const std::vector<int>& counts ) // static
0381 {
0382     unsigned total = SumCounts(counts);
0383     unsigned ni = counts.size();
0384     for(unsigned i=0 ; i < ni ; i++)
0385     {
0386         int np = counts[i] ;
0387         for(int p=0 ; p < np ; p++) seeds.push_back(i) ;
0388     }
0389
0390     bool seeds_expect = seeds.size() == total ;
0391     if(!seeds_expect) std::raise(SIGINT);
0392     assert( seeds_expect );
0393 }
0394
0395 int SEvent::CompareSeeds( const int* seeds, const int* xseeds, int num_seed ) // static
0396 {
0397     int mismatch = 0 ;
0398     for(int i=0 ; i < num_seed ; i++) if( seeds[i] != xseeds[i] ) mismatch += 1 ;
0399     return mismatch ;
0400 }
0401
0402
0403 std::string SEvent::DescSeed( const int* seed, int num_seed, int edgeitems )  // static
0404 {
0405     std::stringstream ss ;
0406     ss << "SEvent::DescSeed num_seed " << num_seed << " (" ;
0407
0408     for(int i=0 ; i < num_seed ; i++)
0409     {
0410         if( i < edgeitems || i > num_seed - edgeitems ) ss << seed[i] << " " ;
0411         else if( i == edgeitems )  ss << "... " ;
0412     }
0413     ss << ")"  ;
0414     std::string s = ss.str();
0415     return s ;
0416 }
0417
0418
0419
0420