eic-opticks/sysrap/SCenterExtentGenstep.cc

0001 #include "scuda.h"
0002 #include "squad.h"
0003 #include "stran.h"
0004
0005 #include "SSys.hh"
0006 #include "SPath.hh"
0007 #include "SEvent.hh"
0008
0009 #include "SGenstep.h"
0010 #include "SFrameGenstep.hh"
0011 #include "SCenterExtentGenstep.hh"
0012
0013 #include "NP.hh"
0014
0015 #include "SLOG.hh"
0016
0017 const plog::Severity SCenterExtentGenstep::LEVEL = SLOG::EnvLevel("SCenterExtentGenstep", "DEBUG" );
0018
0019 const char* SCenterExtentGenstep::BASE = "$TMP/SCenterExtentGenstep"  ;
0020
0021 void SCenterExtentGenstep::save() const
0022 {
0023     save(BASE);
0024 }
0025 void SCenterExtentGenstep::save(const char* dir) const
0026 {
0027     const char* fold = SPath::Resolve(dir, DIRPATH);
0028
0029     LOG(LEVEL) << " saving to " << fold ;
0030
0031     gs->save(fold, "gs.npy" );
0032     save_vec(fold, "isect.npy", ii );
0033     save_vec(fold, "photons.npy", pp );
0034     NP::Write(fold, "peta.npy", (float*)(&peta->q0.f.x), 1, 4, 4 );
0035 }
0036
0037 void SCenterExtentGenstep::save_vec(const char* dir, const char* name, const std::vector<quad4>& vv ) const
0038 {
0039     if(vv.size() == 0)
0040     {
0041         LOG(LEVEL) << " skip as no vec entries for " << name  ;
0042         return ;
0043     }
0044     LOG(LEVEL) << "[ " << name << " size " << vv.size() ;
0045     NP* arr = NP::Make<float>(vv.size(), 4, 4);
0046     LOG(LEVEL) << arr->sstr() ;
0047     arr->read<float>((float*)vv.data());
0048     arr->save(dir, name);
0049     LOG(LEVEL) << "]" ;
0050 }
0051
0052 template<typename T> void SCenterExtentGenstep::set_meta(const char* key, T value )
0053 {
0054     assert(gs);
0055     gs->set_meta<T>(key, value) ;
0056 }
0057
0058
0059 SCenterExtentGenstep::SCenterExtentGenstep(const float4* ce_)
0060     :
0061     gs(nullptr),
0062     gridscale(SSys::getenvfloat("GRIDSCALE", 1.0 )),
0063     peta(new quad4),
0064     ce( ce_ ? *ce_ : make_float4(0.f, 0.f, 0.f, 100.f ))
0065 {
0066     init();
0067 }
0068
0069 /**
0070 SCenterExtentGenstep::init
0071 -----------------------------
0072
0073 Interprets the CEGS envvar defining a 2D plane or 3D grid
0074 of integer genstep positions. The integral coordinates get scaled
0075 by GRIDSCALE to give the float positions.
0076
0077 **/
0078
0079 void SCenterExtentGenstep::init()
0080 {
0081     peta->zero();
0082
0083     LOG(LEVEL) << "[ gridscale " << gridscale  ;
0084
0085     SSys::getenvintvec("CEGS", cegs, ':', "16:0:9:10" );
0086     // input CEGS are 4 or 7 ints delimited by colon nx:ny:nz:num_pho OR nx:px:ny:py:nz:py:num_pho
0087
0088
0089     SFrameGenstep::StandardizeCEGS(cegs);
0090     assert( cegs.size() == 8 );
0091
0092     int ix0 = cegs[0] ;
0093     int ix1 = cegs[1] ;
0094     int iy0 = cegs[2] ;
0095     int iy1 = cegs[3] ;
0096     int iz0 = cegs[4] ;
0097     int iz1 = cegs[5] ;
0098     int photons_per_genstep = cegs[6] ;
0099     int high = cegs[7] ;
0100
0101
0102     nx = (ix1 - ix0)/2 ;
0103     ny = (iy1 - iy0)/2 ;
0104     nz = (iz1 - iz0)/2 ;
0105     int gridaxes = SGenstep::GridAxes(nx, ny, nz);
0106
0107     LOG(LEVEL)
0108         << " nx " << nx
0109         << " ny " << ny
0110         << " nz " << nz
0111         << " GridAxes " << gridaxes
0112         << " GridAxesName " << SGenstep::GridAxesName(gridaxes)
0113         << " high " << high
0114         ;
0115
0116     peta->q0.i.x = ix0 ;
0117     peta->q0.i.y = ix1 ;
0118     peta->q0.i.z = iy0 ;
0119     peta->q0.i.w = iy1 ;
0120
0121     peta->q1.i.x = iz0 ;
0122     peta->q1.i.y = iz1 ;
0123     peta->q1.i.z = photons_per_genstep ;
0124     peta->q1.f.w = gridscale ;
0125
0126
0127     //const Tran<double>* geotran = Tran<double>::make_identity();
0128     const Tran<double>* geotran = Tran<double>::make_translate(ce.x, ce.y, ce.z);
0129
0130
0131     LOG(LEVEL) << "ce (" << ce.x << " " << ce.y << " " << ce.z << " " << ce.w << ")" ;
0132
0133     peta->q2.f.x = ce.x ;   // moved from q1
0134     peta->q2.f.y = ce.y ;
0135     peta->q2.f.z = ce.z ;
0136     peta->q2.f.w = ce.w ;
0137
0138     bool ce_scale = true ;
0139
0140     float3 offset = make_float3(0.f, 0.f, 0.f ) ;
0141     std::vector<float3> ce_offset ;
0142     ce_offset.push_back(offset);
0143
0144     gs = SFrameGenstep::MakeCenterExtentGenstep(ce, cegs, gridscale, geotran, ce_offset, ce_scale, nullptr );
0145
0146     const char* topline = SSys::getenvvar("TOPLINE", "SCenterExtentGenstep.topline") ;
0147     const char* botline = SSys::getenvvar("BOTLINE", "SCenterExtentGenstep.botline" ) ;
0148     set_meta<std::string>("TOPLINE", topline );
0149     set_meta<std::string>("BOTLINE", botline );
0150
0151     SFrameGenstep::GenerateCenterExtentGenstepPhotons( pp, gs, gridscale );
0152
0153     LOG(LEVEL) << "]" ;
0154 }
0155
0156 const char* SCenterExtentGenstep::desc() const
0157 {
0158     std::stringstream ss ;
0159     ss << " CEGS (" ;
0160     for(unsigned i=0 ; i < cegs.size() ; i++ ) ss << cegs[i] << " " ;
0161     ss << ")" ;
0162     ss << " nx " << nx ;
0163     ss << " ny " << ny ;
0164     ss << " nz " << nz ;
0165     ss << " GRIDSCALE " << gridscale ;
0166     ss << " CE ("
0167        << ce.x << " "
0168        << ce.y << " "
0169        << ce.z << " "
0170        << ce.w
0171        << ") "
0172        ;
0173
0174     ss << " gs " << gs->sstr() ;
0175     ss << " pp " << pp.size() ;
0176     ss << " ii " << ii.size() ;
0177
0178     std::string s = ss.str();
0179     return strdup(s.c_str());
0180 }
0181
0182
0183 template void     SCenterExtentGenstep::set_meta<int>(const char*, int );
0184 template void     SCenterExtentGenstep::set_meta<unsigned>(const char*, unsigned );
0185 template void     SCenterExtentGenstep::set_meta<float>(const char*, float );
0186 template void     SCenterExtentGenstep::set_meta<double>(const char*, double );
0187 template void     SCenterExtentGenstep::set_meta<std::string>(const char*, std::string );
0188
0189