eic-opticks/sysrap/SPMT.h

0001 #pragma once
0002 /**
0003 SPMT.h : summarize PMTSimParamData NPFold into the few arrays needed on GPU
0004 ============================================================================
0005
0006 Usage sketch
0007 -------------
0008
0009 ::
0010
0011     SSim::get_spmt
0012     QSim::UploadComponents
0013     QPMT
0014     qpmt
0015
0016
0017 WITH_CUSTOM4
0018 --------------
0019
0020 Custom4 is used only for access to the header-only
0021 TMM complex calculation.
0022 Lack of the WITH_CUSTOM4 flag in sysrap IS A problem as
0023 SPMT.h despite being header only is not directly included
0024 into QPMT it gets compiled into SSim
0025
0026 Aims
0027 ----
0028
0029 1. replace and extend from PMTSim/JPMT.h based off the complete
0030    serialized PMT info rather than the partial NP_PROP_BASE rindex,
0031    thickness info loaded by JPMT.h.
0032
0033    * DONE : (rindex,thickness) matches JPMT after ordering and scale fixes
0034    * TODO : include some pmtcat names in source metadata (in _PMTSimParamData?)
0035             to makes the order more explicit
0036    * DONE : compared get_stackspec scan from JPMT and SPMT
0037
0038 2. DONE :  (17612,2) PMT info arrays with [pmtcat 0/1/2, qescale]
0039 3. DONE : update QPMT.hh/qpmt.h to upload the SPMT.h arrays and test them on device
0040
0041
0042
0043 Related developments
0044 ---------------------
0045
0046 * jcv PMTSimParamSvc     # junosw code that populates the below data core
0047 * jcv PMTSimParamData    # PMT data core knocked out from PMTSimParamSvc for low dependency PMT data access
0048 * jcv _PMTSimParamData   # code that serializes the PMT data core into NPFold
0049 * SPMT.h                 # summarize PMT data into arrays for GPU upload
0050 * QPMT.hh                # CPU QProp setup, uploads
0051 * qpmt.h                 # on GPU use of PMT data
0052 * j/Layr/JPMT.h          # earlier incarnation using "dirty" NP_PROP_BASE approach
0053 * j/Layr/JPMTTest.sh
0054
0055 * Simulation/SimSvc/PMTSimParamSvc/PMTSimParamSvc/tests/PMTSimParamData.sh
0056
0057   * python load the persisted PMTSimParamData
0058
0059 * Simulation/SimSvc/PMTSimParamSvc/PMTSimParamSvc/tests/PMTSimParamData_test.sh
0060
0061   * _PMTSimParamData::Load from "$HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/extra/jpmt/PMTSimParamData"
0062   * test a few simple queries against the loaded PMTSimParamData
0063   * does _PMTSimParamData::Scan_pmtid_qe
0064
0065 * Simulation/SimSvc/PMTSimParamSvc/PMTSimParamSvc/tests/PMTAccessor_test.sh
0066
0067   * PMTAccessor::Load from "$HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/extra/jpmt"
0068   * standalone CPU use of PMTAccessor to do the stack calc
0069
0070 * qudarap/tests/QPMTTest.sh
0071
0072   * formerly used JPMT NP_PROP_BASE loading rindex and thickness, not qeshape and lcqs
0073   * DONE: add in SPMT.h rather than JPMT and extend to include qeshape and lcqs
0074   * DONE: on GPU interpolation check using QPMT
0075
0076
0077 **/
0078
0079 #include <cstdio>
0080 #include <csignal>
0081
0082 #include "NPFold.h"
0083 #include "NPX.h"
0084 #include "scuda.h"
0085 #include "squad.h"
0086 #include "ssys.h"
0087 #include "sproc.h"
0088 #include "spath.h"
0089 #include "s_pmt.h"
0090
0091
0092 #ifdef WITH_CUSTOM4
0093 #include "C4MultiLayrStack.h"
0094 #endif
0095
0096
0097 struct SPMT_Total
0098 {
0099     static constexpr const int FIELDS = 7 ;
0100
0101     int CD_LPMT ;
0102     int SPMT ;
0103     int WP ;
0104     int WP_ATM_LPMT ;
0105     int WP_ATM_MPMT ;
0106     int WP_WAL_PMT ;
0107     int ALL ;
0108
0109     int SUM() const ;
0110     std::string desc() const ;
0111 };
0112
0113
0114 inline int SPMT_Total::SUM() const
0115 {
0116     return CD_LPMT + SPMT + WP + WP_ATM_LPMT + WP_ATM_MPMT + WP_WAL_PMT ;
0117 }
0118 inline std::string SPMT_Total::desc() const
0119 {
0120     std::stringstream ss ;
0121     ss
0122        << "[SPMT_Total.desc\n"
0123 #ifdef WITH_MPMT
0124        << " WITH_MPMT : YES\n"
0125 #else
0126        << " WITH_MPMT : NO\n"
0127 #endif
0128        << std::setw(16) << "CD_LPMT"     << std::setw(7) << CD_LPMT << "\n"
0129        << std::setw(16) << "SPMT"        << std::setw(7) << SPMT << "\n"
0130        << std::setw(16) << "WP"          << std::setw(7) << WP << "\n"
0131        << std::setw(16) << "WP_ATM_LPMT" << std::setw(7) << WP_ATM_LPMT << "\n"
0132        << std::setw(16) << "WP_ATM_MPMT" << std::setw(7) << WP_ATM_MPMT << "\n"
0133        << std::setw(16) << "WP_WAL_PMT"  << std::setw(7) << WP_WAL_PMT << "\n"
0134        << std::setw(16) << "ALL"         << std::setw(7) << ALL << "\n"
0135        << std::setw(16) << "SUM:"        << std::setw(7) << SUM() << "\n"
0136        << std::setw(16) << "SUM()==ALL"  << std::setw(7) << ( SUM() == ALL ? "YES" : "NO " ) << "\n"
0137        << std::setw(16) << "SUM()-ALL"   << std::setw(7) << ( SUM() - ALL ) << "\n"
0138        << "]SPMT_Total.desc\n"
0139        ;
0140     std::string str = ss.str() ;
0141     return str ;
0142 }
0143
0144
0145 struct SPMT_Num
0146 {
0147     static constexpr const int FIELDS = 7 ;
0148
0149     int m_nums_cd_lpmt ;
0150     int m_nums_cd_spmt ;
0151     int m_nums_wp_lpmt ;
0152     int m_nums_wp_atm_lpmt ;
0153     int m_nums_wp_wal_pmt ;
0154     int m_nums_wp_atm_mpmt ;
0155     int ALL ;
0156
0157     int SUM() const ;
0158     std::string desc() const ;
0159 };
0160
0161 inline int SPMT_Num::SUM() const
0162 {
0163     return m_nums_cd_lpmt + m_nums_cd_spmt + m_nums_wp_lpmt + m_nums_wp_atm_lpmt + m_nums_wp_wal_pmt + m_nums_wp_atm_mpmt ;
0164 }
0165
0166 inline std::string SPMT_Num::desc() const
0167 {
0168     std::stringstream ss ;
0169     ss
0170        << "[SPMT_Num.desc\n"
0171 #ifdef WITH_MPMT
0172        << " WITH_MPMT : YES\n"
0173 #else
0174        << " WITH_MPMT : NO\n"
0175 #endif
0176        << std::setw(20) << "m_nums_cd_lpmt     :" << std::setw(7) << m_nums_cd_lpmt << "\n"
0177        << std::setw(20) << "m_nums_cd_spmt     :" << std::setw(7) << m_nums_cd_spmt << "\n"
0178        << std::setw(20) << "m_nums_wp_lpmt     :" << std::setw(7) << m_nums_wp_lpmt << "\n"
0179        << std::setw(20) << "m_nums_wp_atm_lpmt :" << std::setw(7) << m_nums_wp_atm_lpmt << "\n"
0180        << std::setw(20) << "m_nums_wp_wal_pmt  :" << std::setw(7) << m_nums_wp_wal_pmt << "\n"
0181        << std::setw(20) << "m_nums_wp_atm_mpmt :" << std::setw(7) << m_nums_wp_atm_mpmt << "\n"
0182        << std::setw(20) << "ALL                :" << std::setw(7) << ALL << "\n"
0183        << std::setw(20) << "SUM                :" << std::setw(7) << SUM() << "\n"
0184        << std::setw(20) << "SUM==ALL           :" << std::setw(7) << ( SUM() == ALL ? "YES" : "NO " ) << "\n"
0185        << "]SPMT_Num.desc\n"
0186        ;
0187     std::string str = ss.str() ;
0188     return str ;
0189 }
0190
0191
0192
0193
0194
0195 struct SPMT
0196 {
0197     static constexpr const char* _level = "SPMT__level" ;
0198     static const int level ;
0199     static constexpr const float hc_eVnm = 1239.84198433200208455673  ;
0200
0201     enum { L0, L1, L2, L3 } ;
0202     enum { RINDEX, KINDEX } ;
0203
0204     struct LCQS { int lc ; float qs ; } ;
0205
0206
0207     static constexpr const float EN0 = 1.55f ;
0208     static constexpr const float EN1 = 4.20f ;  // 15.5
0209     static constexpr const int   N_EN = 420 - 155 + 1 ;
0210
0211     static constexpr const float WL0 = 440.f ;
0212     static constexpr const float WL1 = 440.f ;
0213     static constexpr const int   N_WL = 1 ;
0214
0215
0216     static constexpr const char* PATH = "$CFBaseFromGEOM/CSGFoundry/SSim/extra/jpmt" ;
0217
0218     // TODO: get these from s_pmt.h also
0219     static constexpr int NUM_PMTCAT = 3 ; // (NNVT, HAMA, NNVT_HiQE)
0220     static constexpr int NUM_LAYER = 4 ;  // (Pyrex, ARC, PHC, Vacuum)
0221     static constexpr int NUM_PROP = 2 ;   // (RINDEX,KINDEX) real and imaginary parts of the index
0222
0223     // below three can be changed via envvars
0224     static const int N_LPMT ;   // N_LPMT must be less than or equal to NUM_CD_LPMT
0225     static const int N_MCT ;
0226     static const int N_SPOL ;
0227
0228     static constexpr const char* QE_shape_PMTCAT_NAMES = "QEshape_NNVT.npy,QEshape_R12860.npy,QEshape_NNVT_HiQE.npy" ;
0229     // follows PMTCategory kPMT enum order but using QEshape array naming convention
0230     // because there is no consistently used naming convention have to do these dirty things
0231
0232     static constexpr const char* QE_shape_S_PMTCAT_NAMES = "QEshape_HZC.npy" ;
0233
0234
0235     static constexpr const char* CE_theta_PMTCAT_NAMES = "CE_NNVTMCP.npy,CE_R12860.npy,CE_NNVTMCP_HiQE.npy" ;
0236     static constexpr const char* CE_costh_PMTCAT_NAMES = "CECOS_NNVTMCP.npy,CECOS_R12860.npy,CECOS_NNVTMCP_HiQE.npy" ;
0237
0238
0239     static const NPFold* CreateFromJPMTAndSerialize(const char* path=nullptr);
0240     static SPMT* CreateFromJPMT(const char* path=nullptr);
0241     SPMT(const NPFold* jpmt);
0242     double get_pmtid_qescale( int pmtid ) const ;
0243
0244     void init();
0245
0246     void init_total();
0247
0248     static NP* FAKE_WHEN_MISSING_PMTParamData_serialize_pmtNum();
0249     static constexpr const char* SPMT__init_pmtNum_FAKE_WHEN_MISSING = "SPMT__init_pmtNum_FAKE_WHEN_MISSING" ;
0250     static int init_pmtNum_FAKE_WHEN_MISSING ;
0251     void init_pmtNum();
0252
0253     void init_pmtCat();
0254
0255     void init_lpmtCat();
0256     void init_qeScale();
0257
0258     void init_rindex_thickness();
0259     static NP* MakeCatPropArrayFromFold( const NPFold* fold, const char* _names, std::vector<const NP*>& v_prop, double domain_scale );
0260     void init_qeshape();
0261     void init_s_qeshape();
0262     void init_cetheta();
0263     void init_cecosth();
0264
0265
0266
0267     void init_lcqs();
0268     void init_s_qescale();
0269
0270     static int TranslateCat(int lpmtcat);
0271
0272     std::string descDetail() const ;
0273     std::string desc() const ;
0274
0275     bool is_complete() const ;
0276     NPFold* serialize() const ;
0277     NPFold* serialize_() const ;
0278
0279     template<typename T>
0280     static T GetFrac(int i, int ni );
0281
0282     template<typename T>
0283     static T GetValueInRange(int j, int nj, T X0, T X1 );
0284
0285
0286     float get_energy(int j, int nj) const ;
0287     float get_wavelength(int j, int nj) const ;
0288
0289     float get_minus_cos_theta_linear_angle(int k, int nk ) const ;
0290     float get_minus_cos_theta_linear_cosine(int k, int nk ) const ;
0291
0292     float get_rindex(int cat, int layr, int prop, float energy_eV) const ;
0293     NP* get_rindex() const ;
0294
0295     float get_qeshape(int cat, float energy_eV) const ;
0296     float get_s_qeshape(int cat, float energy_eV) const ;
0297     NP* get_qeshape() const ;
0298     NP* get_s_qeshape() const ;
0299
0300     float get_thickness_nm(int cat, int layr) const ;
0301     NP* get_thickness_nm() const ;
0302
0303
0304     void get_lpmtid_stackspec( quad4& spec, int lpmtid, float energy_eV) const ;  // EXPT
0305
0306     static constexpr const char* LPMTID_LIST = "0,10,55,98,100,137,1000,10000,17611,50000,51000,52000,52399" ;
0307     static const NP* Make_LPMTID_LIST();
0308
0309
0310 #ifdef WITH_CUSTOM4
0311
0312     /**
0313     SPMTData : Full calculation details for debugging
0314     ----------------------------------------------------
0315     **/
0316     struct SPMTData
0317     {
0318         float4             args ;         //  (1,4)
0319         float4             ARTE ;         //  (1,4)
0320         float4             extra ;        //  (1,4)
0321         quad4              spec ;         //  (4,4)
0322         Stack<float,4>     stack ;        // (44,4)
0323                                           // ------
0324                                           // (51,4)
0325
0326         const float* cdata() const { return &args.x ; }
0327     };
0328
0329     void annotate( NP* art ) const ;
0330
0331     void get_ARTE(SPMTData& pd, int lpmtid, float wavelength_nm, float minus_cos_theta, float dot_pol_cross_mom_nrm ) const ;
0332
0333
0334     NPFold* make_c4scan() const ;
0335 #endif
0336
0337     void get_stackspec( quad4& spec, int cat, float energy_eV) const ;
0338     NP*  get_stackspec() const ;
0339
0340     // lpmtidx : contiguous for CD_LPMT + WP
0341     // lpmtid : non-contiguous CD_LPMT, SPMT, WP
0342
0343     int  get_lpmtcat_from_lpmtidx( int lpmtidx ) const ;
0344     NP*  get_lpmtcat_from_lpmtidx() const ;
0345
0346     int  get_lpmtcat_from_lpmtid(  int lpmtid  ) const ;
0347     int  get_lpmtcat_from_lpmtid(  int* lpmtcat, const int* lpmtid , int num ) const ;
0348     NP*  get_lpmtcat_from_lpmtid(const NP* lpmtid) const ;
0349
0350
0351
0352
0353     float get_qescale_from_lpmtid(int lpmtid) const ;
0354     int   get_copyno_from_lpmtid(int lpmtid) const  ;
0355
0356     float get_qescale_from_lpmtidx(int lpmtidx) const ;
0357     int   get_copyno_from_lpmtidx(int lpmtidx) const ;
0358
0359     NP*   get_qescale_from_lpmtidx() const ;
0360     NP*   get_copyno_from_lpmtidx() const ;
0361
0362     void  get_lcqs_from_lpmtid(  int& lc, float& qs, int lpmtid ) const ;
0363     void  get_lcqs_from_lpmtidx( int& lc, float& qs, int lpmtidx ) const ;
0364     NP*   get_lcqs_from_lpmtidx() const ;
0365
0366     float get_s_qescale_from_spmtid( int pmtid ) const ;
0367     NP*   get_s_qescale_from_spmtid() const ;
0368
0369
0370     float get_pmtcat_qe(int cat, float energy_eV) const ;
0371     NP*   get_pmtcat_qe() const ;
0372
0373     float get_lpmtidx_qe(int lpmtidx, float energy_eV) const ;
0374     NP*   get_lpmtidx_qe() const ;
0375
0376     NPFold* make_testfold() const ;
0377
0378
0379
0380     const char* ExecutableName ;
0381
0382     const NPFold* jpmt ;
0383     const NPFold* PMT_RINDEX ;       // PyrexRINDEX and VacuumRINDEX
0384     const NPFold* PMTSimParamData ;  // lpmtCat for 17612 and many other arrays
0385     const NPFold* PMTParamData ;     // only pmtCat for 45612
0386
0387     SPMT_Num num = {} ;
0388     SPMT_Total total = {} ;
0389
0390     const NPFold* MPT ;              // ARC_RINDEX, ARC_KINDEX, PHC_RINDEX, PHC_KINDEX
0391     const NPFold* CONST ;            // ARC_THICKNESS PHC_THICKNESS
0392     const NPFold* QE_shape ;
0393     const NPFold* CE_theta ;
0394     const NPFold* CE_costh ;
0395
0396     std::vector<const NP*> v_rindex ;
0397     std::vector<const NP*> v_qeshape ;
0398     std::vector<const NP*> v_cetheta ;
0399     std::vector<const NP*> v_cecosth ;
0400
0401     std::vector<const NP*> v_s_qeshape ;
0402
0403
0404     std::vector<LCQS>      v_lcqs ;    // NUM_CD_LPMT + NUM_WP
0405
0406     NP* rindex ;    // (NUM_PMTCAT, NUM_LAYER, NUM_PROP, N_EN, 2:[energy,value] )
0407     NP* qeshape ;   // (NUM_PMTCAT, NUM_SAMPLES~44, 2:[energy,value] )
0408     NP* cetheta ;   // (NUM_PMTCAT, NUM_SAMPLES~9, 2:[angle_radians,value] )
0409     NP* cecosth ;   // (NUM_PMTCAT, NUM_SAMPLES~9, 2:[cosine_angle,value] )
0410     NP* lcqs ;      // (NUM_CD_LPMT + NUM_WP, 2)
0411
0412     NP* thickness ; // (NUM_PMTCAT, NUM_LAYER, 1:value )
0413     float* tt ;
0414
0415     const NP* pmtCat ;
0416     int pmtCat_ni ;
0417     const int* pmtCat_v ;
0418     const NP* pmtNum ;
0419
0420     const NP* lpmtCat ;
0421     const int* lpmtCat_v ;
0422
0423     const NP* qeScale ;
0424     const double* qeScale_v ;
0425
0426     NP* s_qeshape ;  // (NUM_PMTCAT:1, NUM_SAMPLES~60, 2:[energy,value] )
0427     NP* s_qescale ;  // (NUM_SPMT, 1)
0428     float* s_qescale_v ;
0429
0430 };
0431
0432
0433 const int SPMT::level  = ssys::getenvint(_level, 0);
0434 const int SPMT::N_LPMT = ssys::getenvint("N_LPMT", 1 ); // 10 LPMT default for fast scanning
0435 const int SPMT::N_MCT  = ssys::getenvint("N_MCT",  180 );  // "AOI" (actually mct) scan points from -1. to 1.
0436 const int SPMT::N_SPOL = ssys::getenvint("N_SPOL", 1 ); // polarization scan points from S-pol to P-pol
0437
0438 inline const NPFold* SPMT::CreateFromJPMTAndSerialize(const char* path) // static
0439 {
0440     SPMT* spmt = SPMT::CreateFromJPMT(path);
0441     const NPFold* fold = spmt ? spmt->serialize() : nullptr ;
0442     return fold ;
0443 }
0444
0445 /**
0446 SPMT::CreateFromJPMT  (formerly SPMT::Load)
0447 ---------------------------------------------
0448
0449 Default path is::
0450
0451     $CFBaseFromGEOM/CSGFoundry/SSim/extra/jpmt
0452
0453 Which is a directory expected to contain sub-dirs::
0454
0455      PMT_RINDEX
0456      PMTSimParamData
0457
0458      PMTParamData   # this also there, but not used ?
0459
0460
0461 **/
0462
0463
0464 inline SPMT* SPMT::CreateFromJPMT(const char* path_)
0465 {
0466     if(level > 0) printf("[SPMT::CreateFromJPMT [%s]\n", ( path_ == nullptr ? "path_-null" : path_ ));
0467
0468     const char* path = spath::Resolve( path_ != nullptr ? path_ : PATH ) ;
0469     bool unresolved = sstr::StartsWith(path,"CFBaseFromGEOM");
0470     if(unresolved) printf("-SPMT::CreateFromJPMT unresolved path[%s]\n", path) ;
0471     NPFold* fold = NPFold::LoadIfExists(path) ;
0472     if(level > 0) printf("-SPMT::CreateFromJPMT path %s \n", ( path == nullptr ? "path-null" : path ) );
0473     if(level > 0) printf("-SPMT::CreateFromJPMT fold %s \n", ( fold == nullptr ? "fold-null" : "fold-ok" ) );
0474     SPMT* sp = fold ? new SPMT(fold) : nullptr ;
0475
0476     if(level > 0) printf("]SPMT::CreateFromJPMT sp[%s]\n", ( sp ? "YES" : "NO " ) );
0477     return sp ;
0478 }
0479
0480 inline SPMT::SPMT(const NPFold* jpmt_)
0481     :
0482     ExecutableName(sproc::ExecutableName()),
0483     jpmt(jpmt_),
0484     PMT_RINDEX(     jpmt ? jpmt->get_subfold("PMT_RINDEX")      : nullptr ),
0485     PMTSimParamData(jpmt ? jpmt->get_subfold("PMTSimParamData") : nullptr ),
0486     PMTParamData(   jpmt ? jpmt->get_subfold("PMTParamData")    : nullptr ),
0487     MPT(            PMTSimParamData ? PMTSimParamData->get_subfold("MPT")   : nullptr ),
0488     CONST(          PMTSimParamData ? PMTSimParamData->get_subfold("CONST") : nullptr ),
0489     QE_shape(       PMTSimParamData ? PMTSimParamData->get_subfold("QEshape") : nullptr ),
0490     CE_theta(       PMTSimParamData ? PMTSimParamData->get_subfold("CEtheta") : nullptr ),
0491     CE_costh(       PMTSimParamData ? PMTSimParamData->get_subfold("CECOStheta") : nullptr ),
0492     v_lcqs(0),
0493     rindex(nullptr),
0494     qeshape(nullptr),
0495     cetheta(nullptr),
0496     cecosth(nullptr),
0497     lcqs(nullptr),
0498     thickness(NP::Make<float>(NUM_PMTCAT, NUM_LAYER, 1)),
0499     tt(thickness->values<float>()),
0500     pmtCat( PMTParamData ? PMTParamData->get("pmtCat") : nullptr ),
0501     pmtCat_ni( pmtCat ? pmtCat->shape[0] : 0 ),
0502     pmtCat_v( pmtCat ? pmtCat->cvalues<int>() : nullptr ),
0503     pmtNum( PMTParamData ? PMTParamData->get("pmtNum") : nullptr ),
0504     lpmtCat( PMTSimParamData ? PMTSimParamData->get("lpmtCat")  : nullptr ),
0505     lpmtCat_v( lpmtCat ? lpmtCat->cvalues<int>() : nullptr ),
0506     qeScale( PMTSimParamData ? PMTSimParamData->get("qeScale")  : nullptr ),
0507     qeScale_v( qeScale ? qeScale->cvalues<double>() : nullptr ),
0508     s_qeshape(nullptr),
0509     s_qescale(NP::Make<float>(s_pmt::NUM_SPMT,1)),
0510     s_qescale_v( s_qescale ? s_qescale->values<float>() : nullptr )
0511 {
0512     init();
0513 }
0514
0515
0516 /**
0517 SPMT::get_pmtid_qescale
0518 ------------------------
0519
0520 qeScale from  _PMTSimParamData::serialize uses s_pmt.h contiguousidx order : CD_LPMT, WP, SPMT
0521
0522 **/
0523
0524 inline double SPMT::get_pmtid_qescale( int pmtid ) const
0525 {
0526     int contiguousidx = s_pmt::contiguousidx_from_pmtid( pmtid );
0527     float qesc = qeScale_v[contiguousidx] ;
0528     return qesc ;
0529 }
0530
0531
0532
0533
0534 /**
0535 SPMT::init
0536 -----------
0537
0538 Converts PMTSimParamData NPFold into arrays ready
0539 for upload to GPU (by QPMT) with various changes:
0540
0541 1. only 3 LPMT categories selected
0542 2. energy domain scaled to use eV
0543 3. layer thickness changed to use nm
0544 4. arrays are narrowed from double to float
0545
0546 **/
0547
0548 inline void SPMT::init()
0549 {
0550     init_total();
0551
0552     init_pmtNum();
0553     init_pmtCat();
0554
0555     init_lpmtCat();
0556     init_qeScale();
0557
0558     init_rindex_thickness();
0559     init_qeshape();
0560     init_s_qeshape();
0561     init_cetheta();
0562     init_cecosth();
0563
0564     init_lcqs();
0565     init_s_qescale();
0566 }
0567
0568
0569
0570 /**
0571 SPMT::init_total
0572 -----------------
0573
0574 OLD::
0575
0576     SPMT::init_total SPMT_Total CD_LPMT: 17612 SPMT: 25600 WP:  2400 ALL: 45612
0577
0578 SPMT_test.sh::
0579
0580     In [21]: s.jpmt.PMTSimParamData.pmtTotal
0581     Out[21]: array([17612, 25600,  2400,   348,     0,     5, 45965], dtype=int32)
0582
0583     In [22]: s.jpmt.PMTSimParamData.pmtTotal_names
0584     Out[22]: array(['PmtTotal', 'PmtTotal_SPMT', 'PmtTotal_WP', 'PmtTotal_WP_ATM_LPMT', 'PmtTotal_WP_ATM_MPMT', 'PmtTotal_WP_WAL_PMT', 'PmtTotal_ALL'], dtype='<U20')
0585
0586
0587     In [8]: 17612 + 25600 + 2400 + 348 + 5
0588     Out[8]: 45965
0589
0590
0591 The pmtTotal do not yet have the 600 MPMT
0592
0593 Provenance of the pmtTotal::
0594
0595     jcv PMTSimParamData
0596     jcv _PMTSimParamData
0597     jcv PMTSimParamSvc
0598
0599
0600 **/
0601
0602 inline void SPMT::init_total()
0603 {
0604     if( PMTSimParamData == nullptr )
0605     {
0606         std::cout << "SPMT::init_total exit as PMTSimParamData is NULL\n" ;
0607         return ;
0608     }
0609     assert( PMTSimParamData );
0610     const NP* pmtTotal = PMTSimParamData->get("pmtTotal") ; // (7,)   formerly (4,)
0611
0612     if(level > 0) std::cout
0613        << "SPMT::init_total"
0614        << " pmtTotal.sstr " << ( pmtTotal ? pmtTotal->sstr() : "-" )
0615        << "\n"
0616        ;
0617
0618     assert( pmtTotal && pmtTotal->uifc == 'i' && pmtTotal->ebyte == 4 );
0619     assert( pmtTotal->shape.size() == 1 && pmtTotal->shape[0] == SPMT_Total::FIELDS );
0620     assert( pmtTotal->names.size() == SPMT_Total::FIELDS );
0621     [[maybe_unused]] const int* pmtTotal_v = pmtTotal->cvalues<int>();
0622
0623     std::vector<std::string> xnames = {
0624          "PmtTotal",
0625          "PmtTotal_SPMT",
0626          "PmtTotal_WP",
0627          "PmtTotal_WP_ATM_LPMT",
0628          "PmtTotal_WP_ATM_MPMT",
0629          "PmtTotal_WP_WAL_PMT",
0630          "PmtTotal_ALL"
0631        };
0632     assert( pmtTotal->names == xnames );
0633
0634
0635     total.CD_LPMT     = pmtTotal->get_named_value<int>("PmtTotal",      -1) ;
0636     total.SPMT        = pmtTotal->get_named_value<int>("PmtTotal_SPMT", -1) ;
0637     total.WP          = pmtTotal->get_named_value<int>("PmtTotal_WP",   -1) ;
0638     total.WP_ATM_LPMT = pmtTotal->get_named_value<int>("PmtTotal_WP_ATM_LPMT",   -1) ;
0639     total.WP_ATM_MPMT = pmtTotal->get_named_value<int>("PmtTotal_WP_ATM_MPMT",   -1) ;
0640     total.WP_WAL_PMT  = pmtTotal->get_named_value<int>("PmtTotal_WP_WAL_PMT",   -1) ;
0641     total.ALL         = pmtTotal->get_named_value<int>("PmtTotal_ALL",  -1) ;
0642
0643     if(level > 0) std::cout << "SPMT::init_total " << total.desc() << "\n" ;
0644
0645     bool x_total_ALL = total.ALL == total.SUM() ;
0646     if(!x_total_ALL) std::cerr
0647         << "SPMT::init_total"
0648         << " x_total_ALL " << ( x_total_ALL ? "YES" : "NO " )
0649         << total.desc()
0650         << "\n"
0651         ;
0652
0653     assert( x_total_ALL );
0654
0655     assert( pmtTotal_v[0]     == total.CD_LPMT );
0656     assert( total.CD_LPMT     == s_pmt::NUM_CD_LPMT );
0657     assert( total.SPMT        == s_pmt::NUM_SPMT ) ;
0658     assert( total.WP          == s_pmt::NUM_WP  ) ;
0659     assert( total.WP_ATM_LPMT == s_pmt::NUM_WP_ATM_LPMT ) ;
0660     assert( total.WP_WAL_PMT  == s_pmt::NUM_WP_WAL_PMT ) ;
0661
0662 #ifdef WITH_MPMT
0663     bool x_MPMT = total.WP_ATM_MPMT == s_pmt::NUM_WP_ATM_MPMT ;
0664 #else
0665     bool x_MPMT = total.WP_ATM_MPMT == 0 ;
0666 #endif
0667     if(!x_MPMT) std::cerr
0668         << "SPMT::init_total"
0669 #ifdef WITH_MPMT
0670         << " WITH_MPMT "
0671 #else
0672         << " NOT:WITH_MPMT "
0673 #endif
0674         << " x_MPMT " << ( x_MPMT ? "YES" : "NO " )
0675         << " total.WP_ATM_MPMT " << total.WP_ATM_MPMT
0676         << " s_pmt::NUM_WP_ATM_MPMT " << s_pmt::NUM_WP_ATM_MPMT
0677         << total.desc()
0678         ;
0679     assert( x_MPMT ) ;
0680
0681 }
0682
0683
0684
0685 /**
0686 SPMT::FAKE_WHEN_MISSING_PMTParamData_serialize_pmtNum
0687 -------------------------------------------------------
0688
0689 Attempt to give backward compatibility for older junosw branch
0690 that lacks the pmtNum by faking it.
0691
0692 **/
0693
0694 inline NP* SPMT::FAKE_WHEN_MISSING_PMTParamData_serialize_pmtNum() // static
0695 {
0696     NPX::KV<int> kv ;
0697     kv.add("m_nums_cd_lpmt", s_pmt::NUM_CD_LPMT );
0698     kv.add("m_nums_cd_spmt", s_pmt::NUM_SPMT );
0699     kv.add("m_nums_wp_lpmt", s_pmt::NUM_WP );
0700     kv.add("m_nums_wp_atm_lpmt", s_pmt::NUM_WP_ATM_LPMT );
0701     kv.add("m_nums_wp_wal_pmt", s_pmt::NUM_WP_WAL_PMT );
0702     kv.add("m_nums_wp_atm_mpmt",  s_pmt::NUM_WP_ATM_MPMT ); // HMM _ALREADY ?
0703     return kv.values() ;
0704 }
0705
0706
0707 /**
0708 SPMT::init_pmtNum
0709 ------------------
0710
0711 pmtNum and pmtCat info both come from the serialization of PMTParamData
0712 so they are required to correspond to each other
0713
0714 **/
0715
0716 inline int SPMT::init_pmtNum_FAKE_WHEN_MISSING = ssys::getenvint(SPMT__init_pmtNum_FAKE_WHEN_MISSING, 1 );  // default to doing this
0717
0718 inline void SPMT::init_pmtNum()
0719 {
0720     std::cout
0721        << "SPMT::init_pmtNum"
0722        << " pmtNum.sstr " <<  (pmtNum ? pmtNum->sstr() : "-"  )
0723        << " init_pmtNum_FAKE_WHEN_MISSING " << init_pmtNum_FAKE_WHEN_MISSING
0724        << "\n"
0725        ;
0726     if(!pmtNum)
0727     {
0728         if( init_pmtNum_FAKE_WHEN_MISSING == 0 )
0729         {
0730             return ;
0731         }
0732         else
0733         {
0734             pmtNum = FAKE_WHEN_MISSING_PMTParamData_serialize_pmtNum();
0735         }
0736     }
0737
0738     std::vector<std::string> xnames = {
0739          "m_nums_cd_lpmt",
0740          "m_nums_cd_spmt",
0741          "m_nums_wp_lpmt",
0742          "m_nums_wp_atm_lpmt",
0743          "m_nums_wp_wal_pmt",
0744          "m_nums_wp_atm_mpmt"
0745        };
0746     assert( pmtNum->names == xnames );
0747
0748     num.m_nums_cd_lpmt      = pmtNum->get_named_value<int>("m_nums_cd_lpmt", -1) ;
0749     num.m_nums_cd_spmt      = pmtNum->get_named_value<int>("m_nums_cd_spmt", -1) ;
0750     num.m_nums_wp_lpmt      = pmtNum->get_named_value<int>("m_nums_wp_lpmt", -1) ;
0751     num.m_nums_wp_atm_lpmt  = pmtNum->get_named_value<int>("m_nums_wp_atm_lpmt", -1) ;
0752     num.m_nums_wp_wal_pmt   = pmtNum->get_named_value<int>("m_nums_wp_wal_pmt", -1) ;
0753     num.m_nums_wp_atm_mpmt  = pmtNum->get_named_value<int>("m_nums_wp_atm_mpmt", -1) ;
0754     num.ALL = num.SUM();
0755
0756     std::cout
0757         << "[SPMT::init_pmtNum\n"
0758         << num.desc()
0759         << "]SPMT::init_pmtNum\n"
0760         ;
0761
0762     assert( num.m_nums_cd_lpmt     == s_pmt::NUM_CD_LPMT );
0763     assert( num.m_nums_cd_spmt     == s_pmt::NUM_SPMT ) ;
0764     assert( num.m_nums_wp_lpmt     == s_pmt::NUM_WP ) ;
0765     assert( num.m_nums_wp_atm_lpmt == s_pmt::NUM_WP_ATM_LPMT ) ;
0766     assert( num.m_nums_wp_wal_pmt  == s_pmt::NUM_WP_WAL_PMT ) ;
0767     //assert( num.m_nums_wp_atm_mpmt == s_pmt::NUM_WP_ATM_MPMT_ALREADY ) ;
0768
0769 }
0770
0771
0772 /**
0773 SPMT::init_pmtCat
0774 --------------------
0775
0776 Comes from _PMTParamData::serialize::
0777
0778     f->add("pmtCat", NPX::ArrayFromDiscoMapUnordered<int>(data.m_pmt_categories));
0779
0780 ::
0781
0782     In [2]: f.pmtCat[:17612]    ## CD_LPMT
0783     Out[2]:
0784     array([[    0,     1],
0785            [    1,     3],
0786            [    2,     1],
0787            [    3,     3],
0788            [    4,     1],
0789            ...,
0790            [17607,     1],
0791            [17608,     3],
0792            [17609,     1],
0793            [17610,     1],
0794            [17611,     1]], shape=(17612, 2), dtype=int32)
0795
0796     In [7]: f.pmtCat[17612:17612+25600]    ## S_PMT
0797     Out[7]:
0798     array([[20000,     2],
0799            [20001,     2],
0800            [20002,     2],
0801            [20003,     2],
0802            [20004,     2],
0803            ...,
0804            [45595,     2],
0805            [45596,     2],
0806            [45597,     2],
0807            [45598,     2],
0808            [45599,     2]], shape=(25600, 2), dtype=int32)
0809
0810     In [8]: f.pmtCat[17612+25600:17612+25600+2400]   ## WP_PMT
0811     Out[8]:
0812     array([[50000,     3],
0813            [50001,     0],
0814            [50002,     0],
0815            [50003,     3],
0816            [50004,     0],
0817            ...,
0818            [52395,     0],
0819            [52396,     3],
0820            [52397,     3],
0821            [52398,     3],
0822            [52399,     0]], shape=(2400, 2), dtype=int32)
0823
0824
0825
0826 2025/06 WP_ATM and WP_WAL have been added to pmtCat
0827 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0828
0829 PMTParamSvc::init_IDService_WP_ATM_LPMT
0830 PMTParamSvc::init_IDService_WP_WAL_PMT
0831
0832
0833 From SPMT_test.sh::
0834
0835     In [16]: s.jpmt.PMTParamData.pmtCat[17612+25600+2400:17612+25600+2400+348]   ## WP_ATM_LPMT
0836     Out[16]:
0837     array([[52400,     3],
0838            [52401,     0],
0839            [52402,     0],
0840            [52403,     3],
0841            [52404,     3],
0842            ...
0843            [52744,     0],
0844            [52745,     3],
0845            [52746,     0],
0846            [52747,     0]], dtype=int32)
0847
0848
0849     In [18]: s.jpmt.PMTParamData.pmtCat[17612+25600+2400+348:17612+25600+2400+348+5]   ## WP_WAL_PMT
0850     Out[18]:
0851     array([[54000,     3],
0852            [54001,     0],
0853            [54002,     0],
0854            [54003,     0],
0855            [54004,     0]], dtype=int32)
0856
0857
0858
0859 2026/01 yupd_bottompipe_adjust
0860 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0861
0862 In older branch without MPMT fully impl they are already in pmtCat::
0863
0864     In [27]: np.all( f.pmtCat[17612+25600+2400+348:17612+25600+2400+348+600][:,1] == 4 )
0865     Out[27]: np.True_
0866
0867     In [28]: f.pmtCat[17612+25600+2400+348+600:17612+25600+2400+348+600+5]
0868     Out[28]:
0869     array([[54000,     3],
0870            [54001,     0],
0871            [54002,     0],
0872            [54003,     0],
0873            [54004,     0]], dtype=int32)
0874
0875 So the implicit contiguous index of pmtCat follows the order, and included MPMT::
0876
0877     CD-LPMT:17612
0878     CD-SPMT:25600
0879     WP-LPMT:2400
0880     WP-Atmosphere-LPMT:348
0881     WP-Atmosphere-MPMT:600
0882     WP-Water-attenuation-length:5
0883
0884
0885 **/
0886
0887
0888 void SPMT::init_pmtCat()
0889 {
0890     bool expected_type  = pmtCat && pmtCat->uifc == 'i' && pmtCat->ebyte == 4 ;
0891     bool expected_shape = pmtCat && pmtCat->shape.size() == 2 && pmtCat->shape[0] == num.ALL && pmtCat->shape[1] == 2 ;
0892
0893     if(!expected_shape || !expected_type) std::cerr
0894        << "SPMT::init_pmtCat"
0895        << " expected_type " << ( expected_type ? "YES" : "NO " )
0896        << " expected_shape[first-dim-is-num.ALL] " << ( expected_shape ? "YES" : "NO " )
0897        << " pmtCat.shape[0] " << pmtCat->shape[0]
0898        << " num.ALL " << num.ALL
0899        << " pmtCat " << ( pmtCat ? pmtCat->sstr() : "-" )
0900        << " num " << num.desc()
0901        << "\n"
0902        ;
0903
0904     if(!pmtCat) return ;
0905     assert( expected_type );
0906     assert( expected_shape );
0907     assert( pmtCat_v );
0908 }
0909
0910
0911
0912
0913
0914
0915 void SPMT::init_lpmtCat()
0916 {
0917     if(!lpmtCat) return ;
0918     assert( lpmtCat && lpmtCat->uifc == 'i' && lpmtCat->ebyte == 4 );
0919     assert( lpmtCat->shape[0] == s_pmt::NUM_CD_LPMT );
0920     assert( lpmtCat_v );
0921 }
0922
0923 void SPMT::init_qeScale()
0924 {
0925     if(!qeScale) return ;
0926     assert( qeScale && qeScale->uifc == 'f' && qeScale->ebyte == 8 );
0927     assert( qeScale_v );
0928
0929 #ifdef WITH_MPMT
0930     int expect_qeScale_items = s_pmt::NUM_ALL  ;
0931 #else
0932     assert( s_pmt::NUM_CD_LPMT + s_pmt::NUM_SPMT + s_pmt::NUM_WP + s_pmt::NUM_WP_ATM_LPMT + s_pmt::NUM_WP_WAL_PMT  == s_pmt::NUM_ALL_EXCEPT_MPMT );
0933     int expect_qeScale_items = s_pmt::NUM_ALL_EXCEPT_MPMT ;
0934 #endif
0935     bool qeScale_shape_expect = qeScale->shape[0] == expect_qeScale_items ;
0936
0937     if(!qeScale_shape_expect)
0938     std::cerr
0939         << "SPMT::init_qeScale"
0940         << " qeScale.sstr " << ( qeScale ? qeScale->sstr() : "-" )
0941         << " qeScale_shape_expect " << ( qeScale_shape_expect ? "YES" : "NO " )
0942         << " expect_qeScale_items " << expect_qeScale_items
0943         << "\n"
0944         << s_pmt::desc()
0945         ;
0946
0947    assert( qeScale_shape_expect  );
0948 }
0949
0950
0951
0952 /**
0953 SPMT::init_rindex_thickness
0954 ------------------------------
0955
0956 1. this is similar to JPMT::init_rindex_thickness
0957 2. energy domain is scaled to use eV units prior to narrowing from double to float
0958 3. thickness values are scaled to use nm, not m
0959 4. originally there was dumb vacuum rindex with (4,2) for a constant.
0960    Fixed this to (2,2) with NP::MakePCopyNotDumb. Before the fix::
0961
0962     In [19]: t.rindex[0,3]
0963     Out[19]:
0964     array([[[ 1.55,  1.  ],
0965             [ 6.2 ,  1.  ],
0966             [10.33,  1.  ],
0967             [15.5 ,  1.  ],
0968             [ 0.  ,  0.  ],
0969             [ 0.  ,  0.  ],
0970
0971 **/
0972
0973 inline void SPMT::init_rindex_thickness()
0974 {
0975     if(MPT == nullptr || CONST == nullptr) return ;
0976     int MPT_sub = MPT->get_num_subfold() ;
0977     int CONST_items = CONST->num_items() ;
0978
0979     bool MPT_sub_expect = MPT_sub == NUM_PMTCAT  ;
0980     if(!MPT_sub_expect) std::raise(SIGINT);
0981     assert( MPT_sub_expect );    // NUM_PMTCAT:3
0982
0983     bool CONST_items_expect = CONST_items == NUM_PMTCAT  ;
0984
0985     if(!CONST_items_expect) std::cerr
0986         << "SPMT::init_rindex_thickness"
0987         << " CONST_items_expect " << ( CONST_items_expect ? "YES" : "NO " )
0988         << " CONST_items " << CONST_items
0989         << " NUM_PMTCAT " << NUM_PMTCAT
0990         << " MPT_sub " << MPT_sub
0991         << "\n"
0992         ;
0993
0994     if(!CONST_items_expect) std::raise(SIGINT);
0995     assert( CONST_items_expect );
0996
0997     double dscale = 1e-6 ;  // make energy domain scale consistent
0998
0999     for(int i=0 ; i < NUM_PMTCAT ; i++)
1000     {
1001         //const char* name = MPT->get_subfold_key(i) ;
1002         NPFold* pmtcat = MPT->get_subfold(i);
1003         const NP* pmtconst = CONST->get_array(i);
1004
1005         for(int j=0 ; j < NUM_LAYER ; j++)  // NUM_LAYER:4
1006         {
1007             for(int k=0 ; k < NUM_PROP ; k++)   // NUM_PROP:2 (RINDEX,KINDEX)
1008             {
1009                 const NP* v = nullptr ;
1010                 switch(j)
1011                 {
1012                    case 0: v = (k == 0 ? PMT_RINDEX->get("PyrexRINDEX")  : NP::ZEROProp<double>(dscale) ) ; break ;
1013                    case 1: v = pmtcat->get( k == 0 ? "ARC_RINDEX" : "ARC_KINDEX"   )              ; break ;
1014                    case 2: v = pmtcat->get( k == 0 ? "PHC_RINDEX" : "PHC_KINDEX"   )              ; break ;
1015                    case 3: v = (k == 0 ? PMT_RINDEX->get("VacuumRINDEX") : NP::ZEROProp<double>(dscale) ) ; break ;
1016                 }
1017
1018                 NP* vc = NP::MakePCopyNotDumb<double>(v);  // avoid dumb constants with > 2 domain values
1019                 vc->pscale(1e6,0);
1020
1021                 NP* vn = NP::MakeWithType<float>(vc);   // narrow
1022
1023                 v_rindex.push_back(vn) ;
1024             }
1025
1026             double d = 0. ;
1027             //double scale = 1e9 ; // express thickness in nm (not meters)
1028             double scale = 1e6 ;  // HUH: why source scale changed ?
1029             switch(j)
1030             {
1031                case 0: d = 0. ; break ;
1032                case 1: d = scale*pmtconst->get_named_value<double>("ARC_THICKNESS", -1) ; break ;
1033                case 2: d = scale*pmtconst->get_named_value<double>("PHC_THICKNESS", -1) ; break ;
1034                case 3: d = 0. ; break ;
1035             }
1036             tt[i*NUM_LAYER + j] = float(d) ;
1037         }
1038     }
1039     rindex = NP::Combine(v_rindex);
1040     const std::vector<NP::INT>& shape = rindex->shape ;
1041     assert( shape.size() == 3 );
1042     rindex->change_shape( NUM_PMTCAT, NUM_LAYER, NUM_PROP, shape[shape.size()-2], shape[shape.size()-1] );
1043 }
1044
1045
1046
1047 /**
1048 SPMT::MakeCatPropArrayFromFold
1049 --------------------------------
1050
1051 1. selects just the LPMT relevant 3 PMTCAT
1052 2. domain scaling is done whilst still in double precision
1053    and before combination to avoid stomping on last column int
1054    (HMM NP::Combine may have special handling to preserve that now?)
1055
1056 **/
1057
1058 inline NP* SPMT::MakeCatPropArrayFromFold( const NPFold* fold, const char* _names, std::vector<const NP*>& v_prop, double domain_scale ) // static
1059 {
1060     if(fold == nullptr) return nullptr ;
1061     int num_items = fold->num_items();
1062     if(level > 0) std::cout << "SPMT::MakeCatPropArrayFromFold num_items : " << num_items << "\n" ;
1063
1064     std::vector<std::string> names ;
1065     U::Split(_names, ',', names );
1066
1067     for(unsigned i=0 ; i < names.size() ; i++)
1068     {
1069         const char* k = names[i].c_str();
1070         const NP* v = fold->get(k);
1071         if(level > 0) std::cout << std::setw(20) << k << " : " << ( v ? v->sstr() : "-" ) << std::endl ;
1072
1073         NP* vc = v->copy();
1074         vc->pscale(domain_scale, 0);
1075         NP* vn = NP::MakeWithType<float>(vc);   // narrow
1076
1077         v_prop.push_back(vn) ;
1078     }
1079     NP* catprop = NP::Combine(v_prop);
1080     catprop->set_names(names);
1081     return catprop ;
1082 }
1083
1084
1085 inline void SPMT::init_qeshape()
1086 {
1087     double domain_scale = 1e6 ; // MeV to eV
1088     qeshape = MakeCatPropArrayFromFold( QE_shape, QE_shape_PMTCAT_NAMES, v_qeshape, domain_scale );
1089 }
1090 inline void SPMT::init_s_qeshape()
1091 {
1092     double domain_scale = 1e6 ; // MeV to eV
1093     s_qeshape = MakeCatPropArrayFromFold( QE_shape, QE_shape_S_PMTCAT_NAMES, v_s_qeshape, domain_scale );
1094 }
1095
1096
1097
1098
1099 inline void SPMT::init_cetheta()
1100 {
1101     double domain_scale = 1. ;
1102     cetheta = MakeCatPropArrayFromFold( CE_theta, CE_theta_PMTCAT_NAMES, v_cetheta, domain_scale );
1103 }
1104 inline void SPMT::init_cecosth()
1105 {
1106     double domain_scale = 1. ;
1107     cecosth = MakeCatPropArrayFromFold( CE_costh, CE_costh_PMTCAT_NAMES, v_cecosth, domain_scale );
1108 }
1109
1110
1111
1112
1113
1114
1115
1116
1117 /**
1118 SPMT::init_lcqs
1119 -----------------
1120
1121 1. get lpmtCat, qeScale arrays from PMTSimParamData NPFold
1122 2. check appropriate sizes with info for all NUM_CD_LPMT 17612
1123 3. populate v_lcqs vector of LCQS struct holding int:lc
1124    "local 0/1/2 pmtcat" and float:qeScale
1125 4. convert the vector of LCQS struct into lcqs array
1126
1127
1128 s.lcqs
1129 ~~~~~~~
1130
1131 ::
1132
1133     In [3]: s.lcqs[:,1].view(np.float32)
1134     Out[3]:
1135     array([1.025, 1.201, 1.238, 1.172, 1.137, 1.177, 1.098, 1.03 , 1.245, 1.026, 1.162, 1.105, 1.046, 1.042, 1.186, 1.059, ..., 0.919, 0.972, 1.005, 0.932, 0.993, 1.011, 1.019, 1.023, 1.026, 1.104,
1136            0.991, 1.006, 0.98 , 1.053, 1.003, 1.009], shape=(20965,), dtype=float32)
1137
1138     In [4]: s.lcqs[:,1].view(np.float32).min()
1139     Out[4]: np.float32(0.38050073)
1140
1141     In [5]: s.lcqs[:,1].view(np.float32).max()
1142     Out[5]: np.float32(1.5862682)
1143
1144     In [6]: s.lcqs[:,0]
1145     Out[6]:
1146     array([  1,   2,   1,   2,   1,   2,   1,   1,   2,   1,   2,   1,   1,   1,   2,   1, ..., -99, -99, -99, -99, -99, -99, -99, -99, -99, -99, -99,   2,   0,   0,   0,   0],
1147           shape=(20965,), dtype=int32)
1148
1149     In [7]: 17612 + 2400 + 348 + 600 + 5
1150     Out[7]: 20965
1151
1152
1153     In [9]: np.c_[np.unique(s.lcqs[0:17612,0],return_counts=True)]
1154     Out[9]:
1155     array([[   0, 2738],
1156            [   1, 4955],
1157            [   2, 9919]])
1158
1159     In [10]: np.c_[np.unique(s.lcqs[17612:17612+2400,0],return_counts=True)]
1160     Out[10]:
1161     array([[   0, 1836],
1162            [   2,  564]])
1163
1164     In [11]: np.c_[np.unique(s.lcqs[17612+2400:17612+2400+348,0],return_counts=True)]
1165     Out[11]:
1166     array([[  0, 132],
1167            [  2, 216]])
1168
1169     In [14]: np.c_[np.unique(s.lcqs[17612+2400+348:17612+2400+348+600,0],return_counts=True)]
1170     Out[14]: array([[-99, 600]])   ## ALL 600 MPMT with cat -99
1171
1172     In [15]: np.c_[np.unique(s.lcqs[17612+2400+348+600:17612+2400+348+600+5,0],return_counts=True)]
1173     Out[15]:
1174     array([[0, 4],
1175            [2, 1]])
1176
1177
1178 Q: Where is lcqs used ? A: qpmt.h
1179 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1180
1181 ::
1182
1183     qpmt<F>::get_lpmtcat_from_lpmtid
1184     qpmt<F>::get_lpmtcat_from_lpmtidx
1185     qpmt<F>::get_qescale_from_lpmtidx
1186
1187
1188 pmtCat
1189 ~~~~~~~~
1190
1191 pmtCat_v comes from pmtCat array serialized from PMTParamData::m_pmt_categories
1192 whereas most everything else comes from PMTSimParamData
1193
1194 pmtCat is a serialization of m_pmt_categories unordered_map done by NPX::ArrayFromDiscoMapUnordered
1195 which orders the (key,val) array in numerically ascending absolute pmt identifier order.
1196
1197 Ordering is::
1198
1199     CD_LPMT        17612          OFFSET_CD_LPMT:OFFSET_CD_LPMT_END               0:17612
1200     CD_SPMT        25600          OFFSET_CD_SPMT:OFFSET_CD_SPMT_END           20000:45600
1201     WP_LPMT         2400          OFFSET_WP_PMT:OFFSET_WP_PMT_END             50000:52400
1202     WP_ATM_LPMT      348          OFFSET_WP_ATM_LPMT:OFFSET_WP_ATM_LPMT_END   52400:52748
1203     WP_ATM_MPMT      600          OFFSET_WP_ATM_MPMT:OFFSET_WP_ATM_MPMT_END   53000:53600
1204     WP_ATM_WAL         5          OFFSET_WP_WAL_PMT:OFFSET_WP_WAL_PMT_END     54000:54005
1205
1206 Note that MPMT is included even prior to MPMT being fully implemented.
1207
1208 Meanings of quantities used below
1209
1210 lpmtidx
1211    index in range 0:s_pmt::NUM_LPMTIDX
1212
1213 pmtid
1214    id obtained from s_pmt::pmtid_from_lpmtidx
1215
1216
1217 **/
1218
1219 inline void SPMT::init_lcqs()
1220 {
1221
1222     // 0. assert on inputs
1223
1224     assert( PMTSimParamData );
1225     assert( PMTParamData );
1226     assert( s_pmt::NUM_CD_LPMT == total.CD_LPMT );
1227     assert( s_pmt::NUM_WP   == total.WP );
1228
1229     // 1. check consistency between start of pmtCat_v and lpmtCat_v
1230     for(int i=0 ; i < s_pmt::NUM_CD_LPMT ; i++)  // only to 17612
1231     {
1232         assert( pmtCat_v[2*i+1] == lpmtCat_v[i] ) ;
1233     }
1234
1235     // 2. populate v_lcqs array excluding SPMT
1236
1237     v_lcqs.resize(s_pmt::NUM_LPMTIDX);
1238     assert( pmtCat_ni == s_pmt::NUM_OLDCONTIGUOUSIDX ); // MPMT:600 included in pmtCat even prior to full impl
1239
1240     for(int i=0 ; i < s_pmt::NUM_LPMTIDX ; i++ )
1241     {
1242         // outcome depends crucially on the implicit order of lpmtidx from s_pmt.h
1243         int lpmtidx = i ;
1244         int pmtid = s_pmt::pmtid_from_lpmtidx( lpmtidx );
1245         int lpmtidx2 = s_pmt::lpmtidx_from_pmtid( pmtid );
1246         assert( lpmtidx == lpmtidx2 );
1247
1248         // *ix:oldcontiguousidx* follows absolute pmtid numerical order but removes the gaps
1249         // [making it a well founded ordering that applies to all PMTs]
1250         int oldcontiguousidx = s_pmt::oldcontiguousidx_from_pmtid( pmtid ); // offsets and num
1251         assert( oldcontiguousidx < s_pmt::NUM_OLDCONTIGUOUSIDX );
1252         assert( oldcontiguousidx < pmtCat_ni );
1253
1254
1255         int copyno = pmtCat_v[2*oldcontiguousidx+0] ;
1256         int cat    = pmtCat_v[2*oldcontiguousidx+1] ;
1257         bool copyno_pmtid_consistent = copyno == pmtid ;
1258
1259         //const char* delim = "\n" ;
1260         const char* delim = "" ;
1261         if(!copyno_pmtid_consistent)
1262         std::cerr
1263            << "SPMT::init_lcqs" << delim
1264            << " i/lpmtidx "                   << std::setw(6) << lpmtidx << delim
1265            << " lpmtidx2 "                    << std::setw(6) << lpmtidx2 << delim
1266            << " oldcontiguousidx "            << std::setw(6) << oldcontiguousidx << delim
1267            << " s_pmt::NUM_OLDCONTIGUOUSIDX " << std::setw(6) << s_pmt::NUM_OLDCONTIGUOUSIDX << delim
1268            << " s_pmt::NUM_LPMTIDX "          << std::setw(6) << s_pmt::NUM_LPMTIDX << delim
1269            << " pmtid "                       << std::setw(6) << pmtid << delim
1270            << " copyno[pmtCat_v] "            << std::setw(6) << copyno << delim
1271            << " cat[pmtCat_v] "               << std::setw(6) << cat << delim
1272            << " copyno_pmtid_consistent "     << ( copyno_pmtid_consistent ? "YES" : "NO " ) << delim
1273            << "\n"
1274            ;
1275
1276         assert( copyno_pmtid_consistent );
1277
1278         float qesc = get_pmtid_qescale( pmtid );
1279         v_lcqs[lpmtidx] = { TranslateCat(cat), qesc } ;
1280     }
1281     lcqs = NPX::ArrayFromVec<int,LCQS>( v_lcqs ) ;
1282
1283     if(level > 0) std::cout
1284        << "SPMT::init_lcqs" << std::endl
1285        << " NUM_CD_LPMT " << s_pmt::NUM_CD_LPMT << std::endl
1286        << " pmtCat " << ( pmtCat ? pmtCat->sstr() : "-" ) << std::endl
1287        << " lpmtCat " << ( lpmtCat ? lpmtCat->sstr() : "-" ) << std::endl
1288        << " qeScale " << ( qeScale ? qeScale->sstr() : "-" ) << std::endl
1289        << " lcqs " << ( lcqs ? lcqs->sstr() : "-" ) << std::endl
1290        ;
1291
1292     assert( lcqs->shape[0] == s_pmt::NUM_LPMTIDX );
1293
1294     assert( s_pmt::NUM_CD_LPMT == 17612 );
1295     assert( s_pmt::NUM_WP == 2400 );
1296     assert( s_pmt::NUM_WP_ATM_LPMT == 348 );
1297     assert( s_pmt::NUM_WP_WAL_PMT == 5 );
1298
1299 }
1300
1301
1302
1303
1304
1305
1306 /**
1307 SPMT::init_s_qescale
1308 ----------------------
1309
1310 pmtCat
1311     created by _PMTParamData uses "oldcontiguous" order : CD_LPMT, SPMT, WP
1312
1313 qeScale
1314     created by _PMTSimParamData::serialize::
1315
1316         NP* qeScale = NPX::ArrayFromVec<double,double>(data.m_all_pmtID_qe_scale) ;
1317
1318     ordering used by data.m_all_pmtID_qe_scale is SPMT at end (aka contiguousidx)
1319
1320 **/
1321
1322
1323 inline void SPMT::init_s_qescale()
1324 {
1325     int ni = s_qescale ? s_qescale->shape[0] : 0 ;
1326     int nj = s_qescale ? s_qescale->shape[1] : 0 ;
1327     assert( ni == s_pmt::NUM_SPMT );
1328     assert( nj == 1 );
1329     assert( s_qescale_v );
1330
1331     for(int i=0 ; i < ni ; i++ )
1332     {
1333         int spmtidx = i ;
1334         int pmtid = s_pmt::pmtid_from_spmtidx( spmtidx );
1335
1336         // pmtCat from _PMTParamData uses s_pmt.h "oldcontiguous" order : CD_LPMT, SPMT, WP
1337         int oldcontiguousidx = s_pmt::oldcontiguousidx_from_pmtid( pmtid );
1338         int copyno = pmtCat_v[2*oldcontiguousidx+0] ;
1339         int cat    = pmtCat_v[2*oldcontiguousidx+1] ;
1340         assert( copyno == pmtid );
1341         assert( cat == 2 );
1342
1343         double qesc = get_pmtid_qescale( pmtid );
1344
1345         s_qescale_v[spmtidx*nj + 0] = qesc  ;
1346     }
1347 }
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360 /**
1361 SPMT::TranslateCat : 0,1,3 => 0,1,2
1362 --------------------------------------
1363
1364 ::
1365
1366     In [10]: np.c_[np.unique( t.lpmtCat[:,0], return_counts=True )]
1367     Out[10]:
1368     array([[   0, 2720],
1369            [   1, 4997],
1370            [   3, 9895]])
1371
1372
1373     In [4]: t.lcqs[:,0]
1374     Out[4]: array([1, 1, 2, 1, 2, ..., 1, 1, 1, 1, 1], dtype=int32)
1375
1376     In [5]: np.unique(t.lcqs[:,0], return_counts=True)
1377     Out[5]: (array([0, 1, 2], dtype=int32), array([2720, 4997, 9895]))
1378
1379
1380 Including WP should make no difference::
1381
1382     In [19]: np.unique(f.PMTParamData.pmtCat[:17612,1], return_counts=True)
1383     Out[19]: (array([0, 1, 3], dtype=int32), array([2738, 4955, 9919]))
1384
1385     In [20]: np.unique(f.PMTParamData.pmtCat[17612+25600:17612+25600+2400,1],return_counts=True)
1386     Out[20]: (array([0, 3], dtype=int32), array([1836,  564]))
1387
1388
1389 **/
1390
1391 inline int SPMT::TranslateCat(int lpmtcat) // static
1392 {
1393     int lcat = -1 ;
1394     switch( lpmtcat )
1395     {
1396         case -1: lcat = -99  ; break ;   // kPMT_Unknown     see "jcv PMTCategory"
1397         case  0: lcat =  0   ; break ;   // kPMT_NNVT
1398         case  1: lcat =  1   ; break ;   // kPMT_Hamamatsu
1399         case  2: lcat =  -99 ; break ;   // kPMT_HZC
1400         case  3: lcat =  2   ; break ;   // kPMT_NNVT_HighQE
1401         case  4: lcat = -99  ; break ;   // ?MPMT:3?
1402         default: lcat = -99  ; break ;
1403     }
1404     //assert( lcat >= 0 );
1405     return lcat ;
1406 }
1407
1408
1409
1410 inline std::string SPMT::descDetail() const
1411 {
1412     std::stringstream ss ;
1413     ss << "SPMT::descDetail" << std::endl ;
1414     ss << ( jpmt ? jpmt->desc() : "NO_JPMT" ) << std::endl ;
1415     ss << "PMT_RINDEX.desc " << std::endl ;
1416     ss << ( PMT_RINDEX ? PMT_RINDEX->desc() : "NO_PMT_RINDEX" ) << std::endl ;
1417     ss << "PMTSimParamData.desc " << std::endl ;
1418     ss << ( PMTSimParamData ? PMTSimParamData->desc() : "NO_PMTSimParamData" ) << std::endl ;
1419     ss << "jpmt/PMTSimParamData/MPT " << std::endl << std::endl ;
1420     ss << ( MPT ? MPT->desc() : "NO_MPT" ) << std::endl ;
1421     std::string str = ss.str();
1422     return str ;
1423 }
1424
1425 inline std::string SPMT::desc() const
1426 {
1427     std::stringstream ss ;
1428     ss << "SPMT::desc" << std::endl ;
1429     ss << "jpmt.loaddir " << ( jpmt && jpmt->loaddir ? jpmt->loaddir : "NO_LOAD" ) << std::endl ;
1430     ss << "rindex " << ( rindex ? rindex->sstr() : "-" ) << std::endl ;
1431     ss << "thickness " << ( thickness ? thickness->sstr() : "-" ) << std::endl ;
1432     ss << "qeshape " << ( qeshape ? qeshape->sstr() : "-" ) << std::endl ;
1433     ss << "s_qeshape " << ( s_qeshape ? s_qeshape->sstr() : "-" ) << std::endl ;
1434     ss << "s_qescale " << ( s_qescale ? s_qescale->sstr() : "-" ) << std::endl ;
1435     ss << "cetheta " << ( cetheta ? cetheta->sstr() : "-" ) << std::endl ;
1436     ss << "cecosth " << ( cecosth ? cecosth->sstr() : "-" ) << std::endl ;
1437     ss << "lcqs " << ( lcqs ? lcqs->sstr() : "-" ) << std::endl ;
1438
1439     std::string str = ss.str();
1440     return str ;
1441 }
1442
1443 inline bool SPMT::is_complete() const
1444 {
1445     return
1446        rindex != nullptr &&
1447        thickness != nullptr &&
1448        qeshape != nullptr &&
1449        s_qeshape != nullptr &&
1450        s_qescale != nullptr &&
1451        cetheta != nullptr &&
1452        cecosth != nullptr &&
1453        lcqs != nullptr
1454        ;
1455 }
1456
1457 inline NPFold* SPMT::serialize() const   // formerly get_fold
1458 {
1459     return is_complete() ? serialize_() : nullptr ;
1460 }
1461
1462 inline NPFold* SPMT::serialize_() const   // formerly get_fold
1463 {
1464     NPFold* fold = new NPFold ;
1465
1466     if(jpmt) fold->add_subfold("jpmt", const_cast<NPFold*>(jpmt) ) ;
1467
1468     if(rindex) fold->add("rindex", rindex) ;
1469     if(thickness) fold->add("thickness", thickness) ;
1470     if(qeshape) fold->add("qeshape", qeshape) ;
1471     if(s_qeshape) fold->add("s_qeshape", s_qeshape) ;
1472     if(s_qescale) fold->add("s_qescale", s_qescale) ;
1473     if(cetheta) fold->add("cetheta", cetheta) ;
1474     if(cecosth) fold->add("cecosth", cecosth) ;
1475     if(lcqs) fold->add("lcqs", lcqs) ;
1476     return fold ;
1477 }
1478
1479
1480 template<typename T>
1481 inline T SPMT::GetFrac(int i, int ni )  // static
1482 {
1483    return ni == 1 ? 0.f : T(i)/T(ni-1) ;
1484 }
1485
1486 template<typename T>
1487 inline T SPMT::GetValueInRange(int j, int nj, T X0, T X1 )  // static
1488 {
1489     assert( j < nj );
1490     T one(1.);
1491     T fr = GetFrac<T>(j, nj);
1492     T x = X0*(one-fr) + X1*fr ;
1493     return x ;
1494 }
1495
1496 inline float SPMT::get_energy(int j, int nj) const
1497 {
1498     return GetValueInRange<float>(j, nj, EN0, EN1) ;
1499 }
1500 inline float SPMT::get_wavelength(int j, int nj) const
1501 {
1502     return GetValueInRange<float>(j, nj, WL0, WL1) ;
1503 }
1504
1505 inline float SPMT::get_minus_cos_theta_linear_angle(int k, int nk) const
1506 {
1507     assert( k < nk );
1508     float theta_frac = GetFrac<float>(k, nk);
1509     float max_theta_pi = 1.f ;
1510     float theta = theta_frac*max_theta_pi*M_PI ;
1511     float mct = -cos(theta) ;
1512     return mct ;
1513 }
1514
1515 inline float SPMT::get_minus_cos_theta_linear_cosine(int k, int nk) const
1516 {
1517     assert( k < nk );
1518     float fr = GetFrac<float>(k, nk);
1519     float MCT0 = -1.f ;
1520     float MCT1 =  1.f ;
1521     float mct = MCT0*(1.f-fr) + MCT1*fr ;
1522     return mct ;
1523 }
1524
1525 inline float SPMT::get_rindex(int cat, int layr, int prop, float energy_eV) const
1526 {
1527     assert( cat == 0 || cat == 1 || cat == 2 );
1528     assert( layr == 0 || layr == 1 || layr == 2 || layr == 3 );
1529     assert( prop == 0 || prop == 1 );
1530     return rindex->combined_interp_5( cat, layr, prop,  energy_eV ) ;
1531 }
1532
1533 inline NP* SPMT::get_rindex() const
1534 {
1535     std::cout << "SPMT::get_rindex " << std::endl ;
1536
1537     int ni = 3 ;  // pmtcat [0,1,2]
1538     int nj = 4 ;  // layers [0,1,2,3]
1539     int nk = 2 ;  // props [0,1] (RINDEX,KINDEX)
1540     int nl = N_EN ; // energies [0..N_EN-1]
1541     int nn = 2 ;   // payload [energy_eV,rindex_value]
1542
1543     NP* a = NP::Make<float>(ni,nj,nk,nl,nn) ;
1544     float* aa = a->values<float>();
1545
1546     for(int i=0 ; i < ni ; i++)
1547     for(int j=0 ; j < nj ; j++)
1548     for(int k=0 ; k < nk ; k++)
1549     for(int l=0 ; l < nl ; l++)
1550     {
1551         float en = get_energy(l, nl );
1552         float ri = get_rindex(i, j, k, en) ;
1553         int idx = i*nj*nk*nl*nn+j*nk*nl*nn+k*nl*nn+l*nn ;
1554         aa[idx+0] = en ;
1555         aa[idx+1] = ri ;
1556     }
1557     return a ;
1558 }
1559
1560
1561 inline float SPMT::get_qeshape(int cat, float energy_eV) const
1562 {
1563     assert( cat == 0 || cat == 1 || cat == 2 );
1564     return qeshape->combined_interp_3( cat, energy_eV ) ;
1565 }
1566
1567 inline float SPMT::get_s_qeshape(int cat, float energy_eV) const
1568 {
1569     assert( cat == 0  );
1570     return s_qeshape->combined_interp_3( cat, energy_eV ) ;
1571 }
1572
1573
1574
1575
1576
1577 inline NP* SPMT::get_qeshape() const
1578 {
1579     std::cout << "SPMT::get_qeshape " << std::endl ;
1580
1581     int ni = 3 ;   // pmtcat [0,1,2]
1582     int nj = N_EN ; // energies [0..N_EN-1]
1583     int nk = 2 ;   // payload [energy_eV,qeshape_value]
1584
1585     NP* a = NP::Make<float>(ni,nj,nk) ;
1586     float* aa = a->values<float>();
1587
1588     for(int i=0 ; i < ni ; i++)
1589     for(int j=0 ; j < nj ; j++)
1590     {
1591         float en = get_energy(j, nj );
1592         float qe = get_qeshape(i, en) ;
1593         int idx = i*nj*nk+j*nk ;
1594         aa[idx+0] = en ;
1595         aa[idx+1] = qe ;
1596     }
1597     return a ;
1598 }
1599
1600 inline NP* SPMT::get_s_qeshape() const
1601 {
1602     std::cout << "SPMT::get_s_qeshape " << std::endl ;
1603
1604     int ni = 1 ;   // only cat 0 for small PMT
1605     int nj = N_EN ;
1606     int nk = 2 ;   // payload [energy_eV,qeshape_value]
1607
1608     NP* a = NP::Make<float>(ni,nj,nk) ;
1609     float* aa = a->values<float>();
1610
1611     for(int i=0 ; i < ni ; i++)
1612     for(int j=0 ; j < nj ; j++)
1613     {
1614         float en = get_energy(j, nj );
1615         float qe = get_s_qeshape(i, en) ;
1616         int idx = i*nj*nk+j*nk ;
1617         aa[idx+0] = en ;
1618         aa[idx+1] = qe ;
1619     }
1620     return a ;
1621 }
1622
1623
1624
1625
1626
1627 float SPMT::get_thickness_nm(int cat, int lay) const
1628 {
1629     assert( cat == 0 || cat == 1 || cat == 2 );
1630     assert( lay == 0 || lay == 1 || lay == 2 || lay == 3 );
1631     return tt[cat*NUM_LAYER+lay] ;
1632 }
1633
1634 NP* SPMT::get_thickness_nm() const
1635 {
1636     std::cout << "SPMT::get_thickness_nm " << std::endl ;
1637     int ni = NUM_PMTCAT ;
1638     int nj = NUM_LAYER ;
1639     int nk = 1 ;
1640     NP* a = NP::Make<float>(ni, nj, nk );
1641     float* aa = a->values<float>();
1642     for(int i=0 ; i < ni ; i++)
1643     for(int j=0 ; j < nj ; j++)
1644     {
1645         int idx = i*NUM_LAYER + j ;
1646         aa[idx] = get_thickness_nm(i, j);
1647     }
1648     return a ;
1649 }
1650
1651 /**
1652 SPMT::get_pmtid_stackspec
1653 ---------------------------
1654
1655 Expt with using "spare" fourth column .w spec slots to hold onto calculation
1656 intermediates to provide a working context, avoiding API
1657 contortions or repeated lookups.
1658
1659    +----+-------------+----------+------------+-------------+-------------+
1660    |    |    x        |   y      |  z         |  w          |  Notes      |
1661    +====+=============+==========+============+=============+=============+
1662    | q0 | rindex      |  0.f     | 0.f        | i:pmtcat    | Pyrex       |
1663    +----+-------------+----------+------------+-------------+-------------+
1664    | q1 | rindex      | kindex   | thickness  | f:qe_scale  | ARC         |
1665    +----+-------------+----------+------------+-------------+-------------+
1666    | q2 | rindex      | kindex   | thickness  | f:qe        | PHC         |
1667    +----+-------------+----------+------------+-------------+-------------+
1668    | q3 | rindex 1.f  |  0.f     | 0.f        | f:_qe       | Vacuum      |
1669    +----+-------------+----------+------------+-------------+-------------+
1670
1671
1672 1. lookup pmtcat, qe_scale for the pmtid
1673 2. using pmtcat and energy_eV do rindex,kindex interpolation and thickness lookups
1674
1675 Note that the rindex, kindex and thickness only depend on the pmtcat.
1676 However the fourth column qe related values do depend on pmtid with
1677 differences coming in via the qe_scale.
1678
1679 **/
1680 void SPMT::get_lpmtid_stackspec( quad4& spec, int lpmtid, float energy_eV) const
1681 {
1682     spec.zero();
1683
1684     int& cat = spec.q0.i.w ;
1685     float& qe_scale = spec.q1.f.w ;
1686     float& qe = spec.q2.f.w ;
1687     float& _qe = spec.q3.f.w ;
1688     // above are refs to locations currently all holding zero
1689
1690     int lpmtidx = s_pmt::lpmtidx_from_pmtid( lpmtid );
1691     get_lcqs_from_lpmtidx(cat, qe_scale, lpmtidx);
1692
1693     assert( cat > -1 && cat < NUM_PMTCAT );
1694     assert( qe_scale > 0.f );
1695
1696     qe = get_pmtcat_qe(cat, energy_eV ) ; // qeshape interpolation
1697     _qe = qe*qe_scale ;
1698
1699     bool expected_range = _qe > 0.f && _qe < 1.f ;
1700
1701     if(!expected_range) std::cout
1702         << "SPMT::get_pmtid_stackspec"
1703         << " expected_range " << ( expected_range ? "YES" : "NO " )
1704         << " lpmtid " << lpmtid
1705         << " energy_eV " << energy_eV
1706         << " qe " << qe
1707         << " qe_scale " << qe_scale
1708         << " _qe " << _qe
1709         << std::endl
1710         ;
1711
1712     assert( expected_range );
1713
1714
1715     spec.q0.f.x = get_rindex( cat, L0, RINDEX, energy_eV );
1716
1717     spec.q1.f.x = get_rindex(       cat, L1, RINDEX, energy_eV );
1718     spec.q1.f.y = get_rindex(       cat, L1, KINDEX, energy_eV );
1719     spec.q1.f.z = get_thickness_nm( cat, L1 );
1720
1721     spec.q2.f.x = get_rindex(       cat, L2, RINDEX, energy_eV );
1722     spec.q2.f.y = get_rindex(       cat, L2, KINDEX, energy_eV );
1723     spec.q2.f.z = get_thickness_nm( cat, L2 );
1724
1725     spec.q3.f.x = 1.f ; // Vacuum
1726 }
1727
1728
1729 #ifdef WITH_CUSTOM4
1730 /**
1731 SPMT::get_ARTE : TMM MultiLayerStack calculation using complex rindex, thickness, ...
1732 ---------------------------------------------------------------------------------------
1733
1734 Output:ARTE float4
1735     theAbsorption,theReflectivity,theTransmittance,theEfficiency
1736
1737 pmtid
1738     integer in range 0->N_LPMT-1 eg N_LPMT 17612
1739
1740     * used to lookup pmtcat 0,1,2 and qe_scale
1741
1742 energy_eV
1743     float in range 1.55 to 15.5
1744
1745     * used for rindex, kindex, qeshape interpolated property lookups
1746
1747 minus_cos_theta
1748     obtain from dot(mom,nrm) OR OldPolarization*theRecoveredNormal
1749
1750     * expresses the angle of incidence of the photon onto the surface
1751     * -1.f:forward normal incidence
1752     *  0.f:skimming incidence
1753     * +1.f:backward normal incidence
1754
1755 dot_pol_cross_mom_nrm
1756
1757     * dot(pol,cross(mom,nrm)) where pol, mom and nrm are all normalized float3 vectors
1758     * OldPolarization*OldMomentum.cross(theRecoveredNormal) in G4ThreeVector form
1759     * expresses degree of S vs P polarization of the incident photon
1760     * pol,mom,nrm all expected to be normalized vectors
1761     * cross(mom,nrm) vector is transverse to plane of incidence which contains mom and nrm by definition
1762     * cross(mom,nrm) has magnitude sine(angle_between_mom_and_nrm)
1763     * cross(mom,nrm) becomes zero at normal incidence, but S/P have no meaning in that case anyhow
1764     * dot(pol,cross(mom,nrm)) value is cosine(angle_between_pol_and_transverse_vector)*sine(angle_between_mom_and_nrm)
1765     * NB when testing the value of dot_pol_cross_mom_nrm needs to be consistent with minus_cos_theta
1766
1767 **NB : minus_cos_theta AND dot_pol_cross_mom_nrm ARGS ARE RELATED**
1768
1769 dot_pol_cross_mom_nrm includes cross(mom,nrm) and minus_cos_theta is dot(mom,nrm)
1770 so one incorporates the cross product and the other is the dot product
1771 of the same two vectors. Hence care is needed to prepare correctly
1772 related arguments when test scanning the API.
1773
1774 S_pol incorporation
1775 ~~~~~~~~~~~~~~~~~~~~~~
1776
1777 * stack.art.A,R,T now incorporatws S_pol frac obtained from (minus_cos_theta, dot_pol_cross_mom_nrm )
1778
1779
1780 **/
1781
1782 inline void SPMT::get_ARTE(
1783     SPMTData& pd,
1784     int   lpmtid,
1785     float wavelength_nm,
1786     float minus_cos_theta,
1787     float dot_pol_cross_mom_nrm ) const
1788 {
1789     const float energy_eV = hc_eVnm/wavelength_nm ;
1790     get_lpmtid_stackspec(pd.spec, lpmtid, energy_eV);
1791
1792     const float* ss = pd.spec.cdata() ;
1793     const float& _qe = ss[15] ;
1794
1795     pd.args.x = lpmtid ;
1796     pd.args.y = energy_eV ;
1797     pd.args.z = minus_cos_theta ;
1798     pd.args.w = dot_pol_cross_mom_nrm ;
1799
1800     if( minus_cos_theta < 0.f ) // only ingoing photons
1801     {
1802         pd.stack.calc(wavelength_nm, -1.f, 0.f, ss, 16u );
1803         pd.ARTE.w = _qe/pd.stack.art.A ;  // aka theEfficiency and escape_fac, no mct dep
1804
1805         pd.extra.x = 1.f - (pd.stack.art.T_av + pd.stack.art.R_av ) ;  // old An
1806         pd.extra.y = pd.stack.art.A_av ;
1807         pd.extra.z = pd.stack.art.A   ;
1808         pd.extra.w = pd.stack.art.A_s ;
1809     }
1810     else
1811     {
1812         pd.ARTE.w = 0.f ;
1813     }
1814
1815     pd.stack.calc(wavelength_nm, minus_cos_theta, dot_pol_cross_mom_nrm, ss, 16u );
1816
1817     const float& A = pd.stack.art.A ;
1818     const float& R = pd.stack.art.R ;
1819     const float& T = pd.stack.art.T ;
1820
1821     pd.ARTE.x = A ;         // aka theAbsorption
1822     pd.ARTE.y = R/(1.f-A) ; // aka theReflectivity
1823     pd.ARTE.z = T/(1.f-A) ; // aka theTransmittance
1824 }
1825
1826
1827 /**
1828 stack::calc notes
1829 -------------------
1830
1831 0. backwards _qe is set to zero (+ve minus_cos_theta, ie mom with nrm)
1832 1. DONE : compare An with stackNormal.art.A_av
1833 2. DONE : remove _av as at normal incidence S/P are meaningless and give same values anyhow,
1834 3. DONE : removed stackNormal.calc for minus_cos_theta > 0.f
1835 4. DONE : even better get rid of stackNormal by reusing the one stack instance
1836
1837 **/
1838
1839
1840 /**
1841 SPMT::annotate
1842 ----------------
1843
1844 **/
1845
1846 inline void SPMT::annotate( NP* art ) const
1847 {
1848     std::vector<std::pair<std::string, std::string>> kvs =
1849     {
1850         { "title", "SPMT.title" },
1851         { "brief", "SPMT.brief" },
1852         { "name",  "SPMT.name"  },
1853         { "label", "SPMT.label" },
1854         { "ExecutableName", ExecutableName }
1855     };
1856
1857     art->set_meta_kv<std::string>(kvs);
1858 }
1859
1860
1861
1862 inline const NP* SPMT::Make_LPMTID_LIST()  // static
1863 {
1864     const char* lpmtid_list = ssys::getenvvar("LPMTID_LIST", LPMTID_LIST) ;
1865     const NP* lpmtid = NPX::FromString<int>(lpmtid_list,',');
1866     return lpmtid ;
1867 }
1868
1869
1870
1871
1872
1873 /**
1874 SPMT::make_c4scan
1875 -----------------
1876
1877 Scan over (lpmtid, wl, mct, spol )
1878
1879 Potentially could scan over any of (ni,nj,nk,nl)
1880 so should add arrays of the ranges which will have
1881 only one value when not scanned over.
1882
1883 **/
1884
1885
1886 inline NPFold* SPMT::make_c4scan() const
1887 {
1888     if(level > 0) std::cout << "[SPMT::make_c4scan level " << level << std::endl;
1889     NPFold* fold = new NPFold ;
1890
1891     const NP* lpmtid_domain = Make_LPMTID_LIST() ;
1892     const NP* lpmtcat_domain = get_lpmtcat_from_lpmtid( lpmtid_domain );
1893     const NP* mct_domain = NP::MinusCosThetaLinearAngle<float>( N_MCT );
1894     const NP* st_domain = NP::SqrtOneMinusSquare(mct_domain) ;
1895
1896     if(level > 1) std::cout << " N_MCT " << N_MCT << std::endl ;
1897     if(level > 1) std::cout << " mct_domain.desc " << mct_domain->desc() << std::endl ;
1898
1899     fold->add("lpmtid_domain", lpmtid_domain);
1900     fold->add("lpmtcat_domain", lpmtcat_domain);
1901     fold->add("mct_domain", mct_domain );
1902     fold->add("st_domain", st_domain );
1903
1904     SPMTData pd ;
1905
1906     int ni = lpmtid_domain->shape[0] ;
1907     int nj = N_WL ;
1908     int nk = N_MCT ;
1909     int nl = N_SPOL ;
1910
1911
1912
1913     NP* args  = NP::Make<float>(ni, nj, nk, nl, 4 );
1914     NP* ARTE  = NP::Make<float>(ni, nj, nk, nl, 4 );
1915     NP* extra = NP::Make<float>(ni, nj, nk, nl, 4 );
1916     NP* spec  = NP::Make<float>(ni, nj, nk, nl, 4, 4 );
1917
1918     // Make_ allows arbitrary dimensions
1919     NP* stack = NP::Make_<float>(ni, nj, nk, nl, 44, 4 );
1920     NP* ll    = NP::Make_<float>(ni, nj, nk, nl, 4, 4, 4, 2) ;  // (32,4)   4*4*2 = 32
1921     NP* comp  = NP::Make_<float>(ni, nj, nk, nl, 1, 4, 4, 2) ;  // ( 8,4)   4*2 =  8
1922     NP* art   = NP::Make_<float>(ni, nj, nk, nl, 4, 4) ;        // ( 4,4)
1923                                                                 // --------
1924                                                                 //  (44,4)
1925
1926     // stack is composed of (ll,comp,art) but its not easy to use because
1927     // those have different shapes, hence also split them into separate
1928     // (ll,comp,art) arrays for easier querying
1929
1930     annotate(art);
1931
1932     assert( sizeof(pd.stack)/sizeof(float) == 44*4 );
1933
1934     fold->add("args", args);
1935     fold->add("ARTE", ARTE);
1936     fold->add("extra",extra );
1937     fold->add("spec", spec);
1938
1939     fold->add("stack", stack);
1940     fold->add("ll", ll);
1941     fold->add("comp", comp);
1942     fold->add("art", art);
1943
1944     if(level > 1) std::cout << "SPMT::make_c4scan args.sstr " << args->sstr() << std::endl ;
1945
1946     const int* lpmtid_domain_v = lpmtid_domain->cvalues<int>();
1947     const float* mct_domain_v = mct_domain->cvalues<float>();
1948     const float* st_domain_v = st_domain->cvalues<float>();
1949
1950     float* args_v = args->values<float>();
1951     float* ARTE_v = ARTE->values<float>();
1952     float* extra_v = extra->values<float>() ;
1953     float* spec_v  = spec->values<float>() ;
1954
1955     float* stack_v = stack->values<float>() ;
1956     float* ll_v    = ll->values<float>() ;
1957     float* comp_v  = comp->values<float>() ;
1958     float* art_v   = art->values<float>() ;
1959
1960     for(int i=0 ; i < ni ; i++)
1961     {
1962         int lpmtid = lpmtid_domain_v[i] ;
1963         if( i % 100 == 0 && level > 1) std::cout << "SPMT::make_c4scan lpmtid " << lpmtid << std::endl ;
1964         for(int j=0 ; j < nj ; j++)
1965         {
1966            float wavelength_nm = get_wavelength(j, nj );
1967            for(int k=0 ; k < nk ; k++)
1968            {
1969               float minus_cos_theta = mct_domain_v[k] ;
1970               float sin_theta = st_domain_v[k] ;
1971               {
1972                   float minus_cos_theta_0 = get_minus_cos_theta_linear_angle(k, nk );
1973                   float minus_cos_theta_diff = std::abs( minus_cos_theta - minus_cos_theta_0 );
1974                   bool minus_cos_theta_diff_expect = minus_cos_theta_diff < 1e-6 ;
1975                   assert( minus_cos_theta_diff_expect );
1976                   if(!minus_cos_theta_diff_expect) std::cerr << "SPMT::make_c4scan minus_cos_theta_diff_expect " << std::endl ;
1977                   float sin_theta_0 = sqrt( 1.f - minus_cos_theta*minus_cos_theta );
1978                   float sin_theta_diff = std::abs(sin_theta - sin_theta_0)  ;
1979                   bool sin_theta_expect = sin_theta_diff < 1e-6 ;
1980                   if(!sin_theta_expect) std::cout
1981                       << "SPMT::make_c4scan"
1982                       << " k " << k
1983                       << " minus_cos_theta " << std::setw(10) << std::fixed << std::setprecision(5) << minus_cos_theta
1984                       << " sin_theta " << std::setw(10) << std::fixed << std::setprecision(5) << sin_theta
1985                       << " sin_theta_0 " << std::setw(10) << std::fixed << std::setprecision(5) << sin_theta_0
1986                       << " sin_theta_diff " << std::setw(10) << std::fixed << std::setprecision(5) << sin_theta_diff
1987                       << std::endl
1988                       ;
1989                   assert( sin_theta_expect );
1990               }
1991
1992               for(int l=0 ; l < nl ; l++)
1993               {
1994                   float s_pol_frac = GetFrac<float>(l, nl) ;
1995                   float dot_pol_cross_mom_nrm = sin_theta*s_pol_frac ;
1996
1997                   get_ARTE(pd, lpmtid, wavelength_nm, minus_cos_theta, dot_pol_cross_mom_nrm );
1998
1999                   int idx = i*nj*nk*nl*4 + j*nk*nl*4 + k*nl*4 + l*4 ;
2000
2001                   int args_idx = idx ;
2002                   int ARTE_idx = idx ;
2003                   int extra_idx = idx ;
2004                   int spec_idx = idx*4 ;
2005
2006                   int stack_idx = idx*44 ;
2007                   int ll_idx    = idx*4*4*2 ; // 32
2008                   int comp_idx  = idx*1*4*2 ; //  8
2009                   int art_idx   = idx*4 ;     //  4
2010
2011                   memcpy( args_v + args_idx,   &pd.args.x, sizeof(float)*4 );
2012                   memcpy( ARTE_v + ARTE_idx,   &pd.ARTE.x, sizeof(float)*4 );
2013                   memcpy( extra_v + extra_idx, &pd.extra.x,   sizeof(float)*4 );
2014                   memcpy( spec_v + spec_idx,   pd.spec.cdata(), sizeof(float)*4*4 );
2015
2016                   memcpy( stack_v + stack_idx,  pd.stack.cdata(),        sizeof(float)*44*4 );
2017                   memcpy( ll_v    + ll_idx,     pd.stack.ll[0].cdata(),  sizeof(float)*32*4 );
2018                   memcpy( comp_v  + comp_idx,   pd.stack.comp.cdata(),   sizeof(float)*8*4 );
2019                   memcpy( art_v   + art_idx,    pd.stack.art.cdata(),    sizeof(float)*4*4 );
2020
2021                   if(level > 0) std::cout
2022                       << "SPMT::make_c4scan"
2023                       << " i " << std::setw(5) << i
2024                       << " j " << std::setw(5) << j
2025                       << " k " << std::setw(5) << k
2026                       << " l " << std::setw(5) << l
2027                       << " args_idx " << std::setw(10) << args_idx
2028                       << std::endl
2029                       ;
2030
2031               }
2032            }
2033         }
2034     }
2035     if(level > 0) std::cout << "]SPMT::make_c4scan " << std::endl;
2036     return fold ;
2037 }
2038 #endif
2039
2040
2041
2042
2043
2044 void SPMT::get_stackspec( quad4& spec, int cat, float energy_eV) const
2045 {
2046     spec.zero();
2047     spec.q0.f.x = get_rindex( cat, L0, RINDEX, energy_eV );
2048
2049     spec.q1.f.x = get_rindex(       cat, L1, RINDEX, energy_eV );
2050     spec.q1.f.y = get_rindex(       cat, L1, KINDEX, energy_eV );
2051     spec.q1.f.z = get_thickness_nm( cat, L1 );
2052
2053     spec.q2.f.x = get_rindex(       cat, L2, RINDEX, energy_eV );
2054     spec.q2.f.y = get_rindex(       cat, L2, KINDEX, energy_eV );
2055     spec.q2.f.z = get_thickness_nm( cat, L2 );
2056
2057     spec.q3.f.x = get_rindex( cat, L3, RINDEX, energy_eV );
2058     //spec.q3.f.x = 1.f ; // Vacuum, so could just set to 1.f
2059 }
2060
2061 NP* SPMT::get_stackspec() const
2062 {
2063     int ni = NUM_PMTCAT ;
2064     int nj = N_EN ;
2065     int nk = 4 ;
2066     int nl = 4 ;
2067
2068     NP* a = NP::Make<float>(ni, nj, nk, nl );
2069     std::cout << "[ SPMT::get_stackspec " << a->sstr() << std::endl ;
2070
2071     float* aa = a->values<float>();
2072
2073     quad4 spec ;
2074     for(int i=0 ; i < ni ; i++)
2075     for(int j=0 ; j < nj ; j++)
2076     {
2077        float en = get_energy(j, nj );
2078        get_stackspec(spec, i, en );
2079        int idx = i*nj*nk*nl + j*nk*nl ;
2080        memcpy( aa+idx, spec.cdata(), nk*nl*sizeof(float) );
2081     }
2082     std::cout << "] SPMT::get_stackspec " << a->sstr() << std::endl ;
2083     return a ;
2084 }
2085
2086 /**
2087 SPMT::get_lpmtcat_from_lpmtidx
2088 --------------------------------
2089
2090 For lpmtidx(0->17612+2400-1 ) returns 0, 1 or 2 corresponding to NNVT, HAMA, NNVT_HiQE
2091
2092 **/
2093
2094 inline int SPMT::get_lpmtcat_from_lpmtidx(int lpmtidx) const
2095 {
2096     assert( lpmtidx >= 0 && lpmtidx < s_pmt::NUM_LPMTIDX );
2097     const int* lcqs_i = lcqs->cvalues<int>() ;
2098     return lcqs_i[lpmtidx*2+0] ;
2099 }
2100
2101 inline NP* SPMT::get_lpmtcat_from_lpmtidx() const
2102 {
2103     std::cout << "SPMT::get_lpmtcat_from_lpmtidx " << std::endl ;
2104     NP* a = NP::Make<int>( s_pmt::NUM_LPMTIDX ) ;
2105     int* aa = a->values<int>();
2106     for(int i=0 ; i < a->shape[0] ; i++) aa[i] = get_lpmtcat_from_lpmtidx(i) ;
2107     return a ;
2108 }
2109
2110
2111
2112
2113
2114
2115 inline int SPMT::get_lpmtcat_from_lpmtid(int lpmtid) const
2116 {
2117     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
2118     return get_lpmtcat_from_lpmtidx(lpmtidx);
2119 }
2120
2121 inline int SPMT::get_lpmtcat_from_lpmtid( int* lpmtcat_ , const int* lpmtid_ , int num ) const
2122 {
2123     for(int i=0 ; i < num ; i++)
2124     {
2125         int lpmtid = lpmtid_[i] ;
2126         int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
2127         int lpmtcat = get_lpmtcat_from_lpmtidx(lpmtidx) ;
2128         lpmtcat_[i] = lpmtcat ;
2129     }
2130     return num ;
2131 }
2132
2133
2134 inline NP* SPMT::get_lpmtcat_from_lpmtid(const NP* lpmtid ) const
2135 {
2136     assert( lpmtid->shape.size() == 1 );
2137     int num = lpmtid->shape[0] ;
2138     NP* lpmtcat = NP::Make<int>(num);
2139     int num2 = get_lpmtcat_from_lpmtid( lpmtcat->values<int>(), lpmtid->cvalues<int>(), num ) ;
2140
2141     bool num_expect = num2 == num ;
2142     assert( num_expect );
2143     if(!num_expect) std::raise(SIGINT);
2144
2145     return lpmtcat ;
2146 }
2147
2148
2149 inline float SPMT::get_qescale_from_lpmtid(int lpmtid) const
2150 {
2151     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
2152     return get_qescale_from_lpmtidx(lpmtidx);
2153 }
2154 inline int SPMT::get_copyno_from_lpmtid(int lpmtid) const
2155 {
2156     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
2157     int copyno = get_copyno_from_lpmtidx(lpmtidx);
2158     assert( copyno == lpmtid );
2159     return copyno ;
2160 }
2161
2162
2163
2164
2165 inline float SPMT::get_qescale_from_lpmtidx(int lpmtidx) const
2166 {
2167     assert( lpmtidx >= 0 && lpmtidx < s_pmt::NUM_LPMTIDX );
2168     const float* lcqs_f = lcqs->cvalues<float>() ;
2169     return lcqs_f[lpmtidx*2+1] ;
2170 }
2171 inline int SPMT::get_copyno_from_lpmtidx(int lpmtidx) const
2172 {
2173     assert( lpmtidx >= 0 && lpmtidx < s_pmt::NUM_LPMTIDX );
2174     const int* lcqs_i = lcqs->cvalues<int>() ;
2175     return lcqs_i[lpmtidx*2+0] ;
2176 }
2177
2178
2179
2180
2181 inline NP* SPMT::get_qescale_from_lpmtidx() const
2182 {
2183     std::cout << "SPMT::get_qescale_from_lpmtidx " << std::endl ;
2184     NP* a = NP::Make<float>( s_pmt::NUM_LPMTIDX ) ;
2185     float* aa = a->values<float>();
2186     for(int i=0 ; i < a->shape[0] ; i++) aa[i] = get_qescale_from_lpmtidx(i) ;
2187     return a ;
2188 }
2189 inline NP* SPMT::get_copyno_from_lpmtidx() const
2190 {
2191     std::cout << "SPMT::get_copyno_from_lpmtidx " << std::endl ;
2192     NP* a = NP::Make<int>( s_pmt::NUM_LPMTIDX ) ;
2193     int* aa = a->values<int>();
2194     for(int i=0 ; i < a->shape[0] ; i++) aa[i] = get_copyno_from_lpmtidx(i) ;
2195     return a ;
2196 }
2197
2198
2199
2200
2201
2202 /**
2203 SPMT::get_lcqs
2204 ----------------
2205
2206 Accesses the lcqs array returning:
2207
2208 * lc(int): local (0,1,2) lpmt category
2209 * qs(float):qescale
2210
2211 ::
2212
2213     In [5]: np.c_[np.unique( s.lcqs[:,0], return_counts=True )]
2214     Out[5]:
2215     array([[   0, 2720],
2216            [   1, 4997],
2217            [   2, 9895]])
2218
2219
2220
2221
2222 **/
2223 inline void SPMT::get_lcqs_from_lpmtid( int& lc, float& qs, int lpmtid) const
2224 {
2225     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
2226     return get_lcqs_from_lpmtidx(lc, qs, lpmtidx);
2227 }
2228 inline void SPMT::get_lcqs_from_lpmtidx(int& lc, float& qs, int lpmtidx) const
2229 {
2230     assert( lpmtidx >= 0 && lpmtidx < s_pmt::NUM_LPMTIDX );
2231     const int*   lcqs_i = lcqs->cvalues<int>() ;
2232     const float* lcqs_f = lcqs->cvalues<float>() ;
2233     lc = lcqs_i[lpmtidx*2+0] ;
2234     qs = lcqs_f[lpmtidx*2+1] ;
2235 }
2236 /**
2237 Q: Does not this just duplicate lcqs ?
2238 A: YES, that is the point : to check that it does.
2239 **/
2240 inline NP* SPMT::get_lcqs_from_lpmtidx() const
2241 {
2242     std::cout << "SPMT::get_lcqs_from_lpmtidx " << std::endl ;
2243     int ni = s_pmt::NUM_LPMTIDX;
2244     int nj = 2 ;
2245     NP* a = NP::Make<int>(ni, nj) ;
2246     int* ii   = a->values<int>() ;
2247     float* ff = a->values<float>() ;
2248     for(int i=0 ; i < ni ; i++) get_lcqs_from_lpmtidx( ii[i*nj+0], ff[i*nj+1], i );
2249     return a ;
2250 }
2251
2252
2253 inline float SPMT::get_s_qescale_from_spmtid( int pmtid ) const
2254 {
2255     assert( s_pmt::is_spmtid( pmtid ));
2256     int spmtidx = s_pmt::spmtidx_from_pmtid(pmtid);
2257     assert( s_pmt::is_spmtidx( spmtidx ));
2258     float s_qs = s_qescale_v[spmtidx] ;
2259     return s_qs ;
2260 }
2261
2262 inline NP* SPMT::get_s_qescale_from_spmtid() const
2263 {
2264     int ni = s_pmt::NUM_SPMT;
2265     int nj = 1 ;
2266     NP* a = NP::Make<float>(ni, nj) ;
2267     float* aa = a->values<float>() ;
2268     for(int i=0 ; i < ni ; i++)
2269     {
2270         int spmtidx = i ;
2271         int spmtid = s_pmt::pmtid_from_spmtidx(spmtidx);
2272         aa[nj*i+0] = get_s_qescale_from_spmtid( spmtid );
2273     }
2274     return a ;
2275 }
2276
2277
2278
2279
2280 inline float SPMT::get_pmtcat_qe(int cat, float energy_eV) const
2281 {
2282     assert( cat == 0 || cat == 1 || cat == 2 );
2283     return qeshape->combined_interp_3( cat, energy_eV ) ;
2284 }
2285 inline NP* SPMT::get_pmtcat_qe() const
2286 {
2287     std::cout << "SPMT::get_pmtcat_qe " << std::endl ;
2288     int ni = 3 ;
2289     int nj = N_EN ;
2290     int nk = 2 ;
2291
2292     NP* a = NP::Make<float>(ni, nj, nk) ;
2293     float* aa = a->values<float>();
2294
2295     for(int i=0 ; i < ni ; i++)
2296     {
2297         for(int j=0 ; j < nj ; j++)
2298         {
2299             float en = get_energy(j, nj );
2300             aa[i*nj*nk+j*nk+0] = en ;
2301             aa[i*nj*nk+j*nk+1] = get_pmtcat_qe(i, en) ;
2302         }
2303     }
2304     return a ;
2305 }
2306
2307
2308
2309 /**
2310 SPMT::get_lpmtidx_qe
2311 ---------------------
2312
2313 ::
2314
2315     q = t.test.get_pmtid_qe
2316
2317     In [21]: q.shape
2318     Out[21]: (17612, 266, 2)
2319
2320     In [22]: np.argmax(q[:,:,0], axis=1)
2321     Out[22]: array([265, 265, 265, 265, 265, ..., 265, 265, 265, 265, 265])
2322
2323     In [23]: np.argmax(q[:,:,1], axis=1)
2324     Out[23]: array([163, 163, 163, 163, 163, ..., 163, 163, 163, 163, 163])
2325
2326 * surprised that max qe at same energy for all 17612 pmtid ?
2327
2328 **/
2329
2330 inline float SPMT::get_lpmtidx_qe(int lpmtidx, float energy_eV) const
2331 {
2332     int cat(-1) ;
2333     float qe_scale(-1.f) ;
2334
2335     get_lcqs_from_lpmtidx(cat, qe_scale, lpmtidx);
2336
2337     assert( cat == 0 || cat == 1 || cat == 2 );
2338     assert( qe_scale > 0.f );
2339
2340     float qe = get_pmtcat_qe(cat, energy_eV);
2341     qe *= qe_scale ;
2342
2343     return qe ;
2344 }
2345
2346
2347
2348
2349 inline NP* SPMT::get_lpmtidx_qe() const
2350 {
2351     std::cout << "SPMT::get_lpmtidx_qe " << std::endl ;
2352     int ni = N_LPMT  ;
2353     int nj = N_EN ;
2354     int nk = 2 ;
2355
2356     NP* a = NP::Make<float>(ni, nj, nk) ;
2357     float* aa = a->values<float>();
2358
2359     for(int i=0 ; i < ni ; i++)
2360     for(int j=0 ; j < nj ; j++)
2361     {
2362         float en = get_energy(j, nj );
2363         int idx = i*nj*nk+j*nk ;
2364         aa[idx+0] = en ;
2365         aa[idx+1] = get_lpmtidx_qe(i, en) ;
2366     }
2367     return a ;
2368 }
2369
2370
2371 /**
2372 SPMT::make_testfold
2373 ---------------------
2374
2375 Invoked from main of sysrap/tests/SPMT_test.cc
2376
2377 **/
2378
2379
2380 inline NPFold* SPMT::make_testfold() const
2381 {
2382     std::cout << "[ SPMT::make_testfold " << std::endl ;
2383
2384     NPFold* f = new NPFold ;
2385
2386     f->add("get_lpmtcat_from_lpmtidx", get_lpmtcat_from_lpmtidx() );
2387     f->add("get_qescale_from_lpmtidx", get_qescale_from_lpmtidx() );
2388     f->add("get_lcqs_from_lpmtidx", get_lcqs_from_lpmtidx() );
2389     f->add("get_pmtcat_qe", get_pmtcat_qe() );
2390     f->add("get_lpmtidx_qe", get_lpmtidx_qe() );
2391
2392     f->add("get_rindex", get_rindex() );
2393     f->add("get_qeshape", get_qeshape() );
2394     f->add("get_s_qeshape", get_s_qeshape() );
2395     f->add("get_s_qescale_from_spmtid", get_s_qescale_from_spmtid() );
2396     f->add("get_thickness_nm", get_thickness_nm() );
2397     f->add("get_stackspec", get_stackspec() );
2398
2399 #ifdef WITH_CUSTOM4
2400     f->add_subfold("c4scan", make_c4scan() );
2401 #endif
2402
2403     std::cout << "] SPMT::make_testfold " << std::endl ;
2404     return f ;
2405 }