eic-opticks/qudarap/qpmt.h

0001 #pragma once
0002 /**
0003 qpmt.h
0004 =======
0005
0006
0007 **/
0008
0009 #if defined(__CUDACC__) || defined(__CUDABE__)
0010    #define QPMT_METHOD __device__
0011 #else
0012    #define QPMT_METHOD
0013 #endif
0014
0015
0016
0017
0018
0019 #if defined(__CUDACC__) || defined(__CUDABE__)
0020 #else
0021 #include "QUDARAP_API_EXPORT.hh"
0022 #endif
0023
0024
0025 template <typename T> struct qprop ;
0026
0027 #include "scuda.h"
0028 #include "squad.h"
0029 #include "qprop.h"
0030
0031 #ifdef WITH_CUSTOM4
0032 #include "C4MultiLayrStack.h"
0033 #endif
0034
0035
0036 #include "s_pmt.h"
0037
0038
0039 template<typename F>
0040 struct qpmt
0041 {
0042     enum { L0, L1, L2, L3 } ;
0043
0044     static constexpr const F hc_eVnm = 1239.84198433200208455673  ;
0045     static constexpr const F zero = 0. ;
0046     static constexpr const F one = 1. ;
0047     // constexpr should mean any double conversions happen at compile time ?
0048
0049     qprop<F>* rindex_prop ;
0050     qprop<F>* qeshape_prop ;
0051     qprop<F>* cetheta_prop ;
0052     qprop<F>* cecosth_prop ;
0053
0054     F*        thickness ;
0055     F*        lcqs ;
0056     int*      i_lcqs ;  // int* "view" of lcqs memory
0057
0058     qprop<F>* s_qeshape_prop ;
0059     F*        s_qescale ;
0060
0061
0062 #if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
0063     // loosely follow SPMT.h
0064     QPMT_METHOD int  get_lpmtcat_from_lpmtid(  int lpmtid  ) const  ;
0065     QPMT_METHOD int  get_lpmtcat_from_lpmtidx( int lpmtidx ) const  ;
0066     QPMT_METHOD F    get_qescale_from_lpmtid(  int lpmtid  ) const  ;
0067     QPMT_METHOD F    get_qescale_from_lpmtidx( int lpmtidx ) const  ;
0068
0069     QPMT_METHOD F    get_s_qescale_from_spmtid(  int spmtid  ) const  ;
0070
0071
0072     QPMT_METHOD F    get_lpmtcat_qe( int pmtcat, F energy_eV ) const ;
0073     QPMT_METHOD F    get_lpmtcat_ce( int pmtcat, F theta ) const ;
0074
0075     QPMT_METHOD F    get_lpmtcat_rindex(    int lpmtcat, int layer, int prop, F energy_eV ) const ;
0076     QPMT_METHOD F    get_lpmtcat_rindex_wl( int lpmtcat, int layer, int prop, F wavelength_nm ) const ;
0077
0078
0079     QPMT_METHOD void get_lpmtcat_stackspec( F* spec16, int pmtcat, F energy_eV ) const ;
0080
0081     QPMT_METHOD void get_lpmtid_stackspec(          F* spec16, int lpmtid, F energy_eV ) const ;
0082
0083
0084 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0085     QPMT_METHOD void get_lpmtid_stackspec_ce_acosf( F* spec16, int lpmtid, F energy_eV, F lposcost, unsigned pidx, bool pidx_debug ) const ;
0086 #else
0087     QPMT_METHOD void get_lpmtid_stackspec_ce_acosf( F* spec16, int lpmtid, F energy_eV, F lposcost ) const ;
0088 #endif
0089
0090     QPMT_METHOD void get_lpmtid_stackspec_ce(       F* spec15, int lpmtid, F energy_eV, F lposcost ) const ;
0091
0092
0093 #ifdef WITH_CUSTOM4
0094     QPMT_METHOD void get_lpmtid_SPEC(   F* spec16 , int lpmtid, F wavelength_nm ) const ;
0095
0096 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0097     QPMT_METHOD void get_lpmtid_SPEC_ce(F* spec16 , int lpmtid, F wavelength_nm, F lposcost, unsigned pidx, bool pidx_debug ) const ;
0098 #else
0099     QPMT_METHOD void get_lpmtid_SPEC_ce(F* spec16 , int lpmtid, F wavelength_nm, F lposcost ) const ;
0100 #endif
0101
0102
0103     QPMT_METHOD void get_lpmtid_LL(  F* ll_128  , int lpmtid, F wavelength_nm, F minus_cos_theta, F dot_pol_cross_mom_nrm ) const ;
0104     QPMT_METHOD void get_lpmtid_COMP(F* comp_32 , int lpmtid, F wavelength_nm, F minus_cos_theta, F dot_pol_cross_mom_nrm ) const ;
0105     QPMT_METHOD void get_lpmtid_ART( F* art_16  , int lpmtid, F wavelength_nm, F minus_cos_theta, F dot_pol_cross_mom_nrm ) const ;
0106     QPMT_METHOD void get_lpmtid_ARTE(  F* ARTE  , int lpmtid, F wavelength_nm, F minus_cos_theta, F dot_pol_cross_mom_nrm ) const ;
0107
0108 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0109     QPMT_METHOD void get_lpmtid_ATQC(  F* ATQC, int lpmtid, F wavelength_nm, F minus_cos_theta, F dot_pol_cross_mom_nrm, F lposcost, unsigned pidx, bool pidx_debug ) const ;
0110 #else
0111     QPMT_METHOD void get_lpmtid_ATQC(  F* ATQC  , int lpmtid, F wavelength_nm, F minus_cos_theta, F dot_pol_cross_mom_nrm, F lposcost ) const ;
0112 #endif
0113
0114 // end WITH_CUSTOM3
0115 #endif
0116
0117 // end __CUDACC__ etc..
0118 #endif
0119 };
0120
0121 #if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
0122
0123 template<typename F>
0124 inline QPMT_METHOD int qpmt<F>::get_lpmtcat_from_lpmtid( int lpmtid ) const
0125 {
0126     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
0127     return lpmtidx < s_pmt::NUM_LPMTIDX && lpmtidx > -1 ? i_lcqs[lpmtidx*2+0] : -2 ;
0128     // extended from former NUM_CD_LPMT_AND_WP after s_pmt.h overhaul for WP_ATM_LPMT WP_WAL_PMT
0129 }
0130
0131 template<typename F>
0132 inline QPMT_METHOD int qpmt<F>::get_lpmtcat_from_lpmtidx( int lpmtidx ) const
0133 {
0134     return lpmtidx < s_pmt::NUM_LPMTIDX && lpmtidx > -1 ? i_lcqs[lpmtidx*2+0] : -2 ;
0135 }
0136
0137 template<typename F>
0138 inline QPMT_METHOD F qpmt<F>::get_qescale_from_lpmtid( int lpmtid ) const
0139 {
0140     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
0141     return lpmtidx < s_pmt::NUM_LPMTIDX && lpmtidx > -1 ? lcqs[lpmtidx*2+1] : -2.f ;
0142 }
0143 template<typename F>
0144 inline QPMT_METHOD F qpmt<F>::get_qescale_from_lpmtidx( int lpmtidx ) const
0145 {
0146     return lpmtidx < s_pmt::NUM_LPMTIDX && lpmtidx > -1 ? lcqs[lpmtidx*2+1] : -2.f ;
0147 }
0148
0149
0150 template<typename F>
0151 inline QPMT_METHOD F qpmt<F>::get_s_qescale_from_spmtid( int spmtid ) const
0152 {
0153     int spmtidx = s_pmt::spmtidx_from_pmtid(spmtid);
0154     return spmtidx < s_pmt::NUM_SPMT && spmtidx > -1 ? s_qescale[spmtidx] : -2.f ;
0155 }
0156
0157
0158
0159
0160
0161
0162
0163 template<typename F>
0164 inline QPMT_METHOD F qpmt<F>::get_lpmtcat_qe( int lpmtcat, F energy_eV ) const
0165 {
0166     return lpmtcat > -1 && lpmtcat < s_pmt::NUM_CAT ? qeshape_prop->interpolate( lpmtcat, energy_eV ) : -1.f ;
0167 }
0168
0169 /**
0170 qpmt::get_lpmtcat_ce
0171 ---------------------
0172
0173 theta_radians range 0->pi/2
0174
0175
0176               cos(th)=1
0177                 th=0
0178
0179                   |
0180                +  |  +
0181            +      |      +
0182                   |
0183          +        |         +
0184         +---------+----------+--   th = pi/2
0185                                 cos(th)=0
0186
0187 **/
0188
0189 template<typename F>
0190 inline QPMT_METHOD F qpmt<F>::get_lpmtcat_ce( int lpmtcat, F theta_radians ) const
0191 {
0192     //return lpmtcat > -1 && lpmtcat < qpmt_NUM_CAT ? cetheta_prop->interpolate( lpmtcat, theta_radians ) : -1.f ;
0193     return lpmtcat > -1 && lpmtcat < s_pmt::NUM_CAT ? cetheta_prop->interpolate( lpmtcat, theta_radians ) : -1.f ;
0194 }
0195 template<typename F>
0196 inline QPMT_METHOD F qpmt<F>::get_lpmtcat_rindex( int lpmtcat, int layer, int prop, F energy_eV ) const
0197 {
0198     //const unsigned idx = lpmtcat*qpmt_NUM_LAYR*qpmt_NUM_PROP + layer*qpmt_NUM_PROP + prop ;
0199     const unsigned idx = lpmtcat*s_pmt::NUM_LAYR*s_pmt::NUM_PROP + layer*s_pmt::NUM_PROP + prop ;
0200     return rindex_prop->interpolate( idx, energy_eV ) ;
0201 }
0202 template<typename F>
0203 inline QPMT_METHOD F qpmt<F>::get_lpmtcat_rindex_wl( int lpmtcat, int layer, int prop, F wavelength_nm ) const
0204 {
0205     const F energy_eV = hc_eVnm/wavelength_nm ;
0206     return get_lpmtcat_rindex(lpmtcat, layer, prop, energy_eV );
0207 }
0208
0209
0210
0211 template<typename F>
0212 inline QPMT_METHOD void qpmt<F>::get_lpmtcat_stackspec( F* spec, int lpmtcat, F energy_eV ) const
0213 {
0214     //const unsigned idx = lpmtcat*qpmt_NUM_LAYR*qpmt_NUM_PROP ;
0215     //const unsigned idx0 = idx + L0*qpmt_NUM_PROP ;
0216     //const unsigned idx1 = idx + L1*qpmt_NUM_PROP ;
0217     //const unsigned idx2 = idx + L2*qpmt_NUM_PROP ;
0218
0219     const unsigned idx = lpmtcat*s_pmt::NUM_LAYR*s_pmt::NUM_PROP ;
0220
0221     const unsigned idx0 = idx + L0*s_pmt::NUM_PROP ;
0222     const unsigned idx1 = idx + L1*s_pmt::NUM_PROP ;
0223     const unsigned idx2 = idx + L2*s_pmt::NUM_PROP ;
0224
0225
0226     spec[0*4+0] = rindex_prop->interpolate( idx0+0u, energy_eV );
0227     spec[0*4+1] = zero ;
0228     spec[0*4+2] = zero ;
0229
0230     spec[1*4+0] = rindex_prop->interpolate( idx1+0u, energy_eV );
0231     spec[1*4+1] = rindex_prop->interpolate( idx1+1u, energy_eV );
0232     spec[1*4+2] = thickness[lpmtcat*s_pmt::NUM_LAYR+L1] ;
0233
0234     spec[2*4+0] = rindex_prop->interpolate( idx2+0u, energy_eV );
0235     spec[2*4+1] = rindex_prop->interpolate( idx2+1u, energy_eV );
0236     spec[2*4+2] = thickness[lpmtcat*s_pmt::NUM_LAYR+L2] ;
0237
0238     spec[3*4+0] = one ;  // Vacuum RINDEX
0239     spec[3*4+1] = zero ;
0240     spec[3*4+2] = zero ;
0241
0242     // "4th" column untouched, as pmtid info goes in there
0243 }
0244
0245
0246
0247 template<typename F>
0248 inline QPMT_METHOD void qpmt<F>::get_lpmtid_stackspec( F* spec, int lpmtid, F energy_eV ) const
0249 {
0250     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
0251     const int& lpmtcat = i_lcqs[lpmtidx*2+0] ;
0252     if(lpmtcat < 0) printf("//qpmt::get_lpmtidx_stackspec lpmtid %d lpmtcat %d \n", lpmtid, lpmtcat );
0253
0254     const F& qe_scale = lcqs[lpmtidx*2+1] ;
0255     const F qe_shape = qeshape_prop->interpolate( lpmtcat, energy_eV ) ;
0256     const F qe = qe_scale*qe_shape ;
0257
0258     spec[0*4+3] = lpmtcat ;  // profligate int to float, s.spec[:,:,0,3].astype(np.int32)
0259     spec[1*4+3] = qe_scale ;
0260     spec[2*4+3] = qe_shape ;
0261     spec[3*4+3] = qe ;
0262
0263     get_lpmtcat_stackspec( spec, lpmtcat, energy_eV );
0264 }
0265
0266
0267 /**
0268 get_lpmtid_stackspec_ce_acosf (see alternative get_lpmtid_stackspec_ce)
0269 -------------------------------------------------------------------------
0270
0271 This uses cetheta_prop interpolation forcing use of acosf to get theta
0272
0273 lposcost
0274     local position cosine theta,
0275     expected range 1->0 (as front of PMT is +Z)
0276     so theta_radians expected 0->pi/2
0277
0278 Currently called by qpmt::get_lpmtid_ATQC
0279
0280 **/
0281
0282
0283 template<typename F>
0284 inline QPMT_METHOD void qpmt<F>::get_lpmtid_stackspec_ce_acosf(
0285     F* spec
0286     , int lpmtid
0287     , F energy_eV
0288     , F lposcost
0289 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0290     , unsigned pidx
0291     , bool pidx_debug
0292 #endif
0293     ) const
0294 {
0295     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);  // lpmtidx:-1 FOR SPMT, BUT THAT SHOULD NEVER HAPPEN
0296     const int& lpmtcat = i_lcqs[lpmtidx*2+0] ;
0297
0298     const F& qe_scale = lcqs[lpmtidx*2+1] ;
0299     const F qe_shape = qeshape_prop->interpolate( lpmtcat, energy_eV ) ;
0300     const F qe = qe_scale*qe_shape ;
0301
0302     const F theta_radians = acosf(lposcost);
0303     const F ce = cetheta_prop->interpolate( lpmtcat, theta_radians );
0304
0305     spec[0*4+3] = lpmtcat ;       //  3
0306     spec[1*4+3] = ce ;            //  7
0307     spec[2*4+3] = qe_shape ;      // 11
0308     spec[3*4+3] = qe ;            // 15
0309
0310     get_lpmtcat_stackspec( spec, lpmtcat, energy_eV );
0311
0312
0313 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0314     if(pidx_debug) printf("//qpmt::get_lpmtid_stackspec_ce_acosf pidx %7d lpmtid %6d lpmtidx %6d lpmtcat %2d qe_scale %8.4f energy_eV %8.4f qe_shape %8.4f qe(qe_scale*qe_shape) %8.4f spec[15] %8.4f lposcost %8.4f theta_radians %8.4f ce %8.4f \n",
0315                                                                       pidx, lpmtid, lpmtidx, lpmtcat,  qe_scale, energy_eV, qe_shape, qe, spec[15], lposcost, theta_radians, ce );
0316 #endif
0317
0318
0319
0320 }
0321
0322
0323
0324 /**
0325 get_lpmtid_stackspec_ce
0326 -------------------------
0327
0328 This uses cecosth_prop interpolation avoiding use of acosf
0329 as directly interpolate in the cosine.
0330
0331 lposcost
0332     local position cosine theta,
0333     expected range 1->0 (as front of PMT is +Z)
0334     so theta_radians expected 0->pi/2
0335
0336 Potentially called by qpmt::get_lpmtid_ATQC
0337
0338 **/
0339
0340
0341 template<typename F>
0342 inline QPMT_METHOD void qpmt<F>::get_lpmtid_stackspec_ce( F* spec, int lpmtid, F energy_eV, F lposcost ) const
0343 {
0344     int lpmtidx = s_pmt::lpmtidx_from_pmtid(lpmtid);
0345     const int& lpmtcat = i_lcqs[lpmtidx*2+0] ;
0346     const F& qe_scale = lcqs[lpmtidx*2+1] ;
0347
0348     const F qe_shape = qeshape_prop->interpolate( lpmtcat, energy_eV ) ;
0349     const F qe = qe_scale*qe_shape ;
0350
0351     const F ce = cecosth_prop->interpolate( lpmtcat, lposcost );
0352
0353     spec[0*4+3] = lpmtcat ;       //  3
0354     spec[1*4+3] = ce ;            //  7
0355     spec[2*4+3] = qe_shape ;      // 11
0356     spec[3*4+3] = qe ;            // 15
0357
0358     get_lpmtcat_stackspec( spec, lpmtcat, energy_eV );
0359
0360     //printf("//qpmt::get_lpmtid_stackspec_ce lpmtid %d lpmtidx %d lpmtcat %d lposcost %7.3f  ce %7.3f spec[7] %7.3f \n", lpmtid, lpmtidx, lpmtcat, lposcost, ce, spec[7] );
0361 }
0362
0363
0364
0365
0366
0367
0368
0369 #ifdef WITH_CUSTOM4
0370
0371 /**
0372 qpmt::get_lpmtid_SPEC
0373 -----------------------
0374
0375 Gets the inputs to the TMM stack calculation
0376 customized for particular lpmtid.
0377
0378 1. HMM: has to convert wavelength_nm to energy_eV
0379    for every photon... could convert to wavelength domain CPU side
0380    to avoid the need to do that
0381
0382 **/
0383
0384 template<typename F>
0385 inline QPMT_METHOD void qpmt<F>::get_lpmtid_SPEC(
0386     F* spec_16,
0387     int lpmtid,
0388     F wavelength_nm ) const
0389 {
0390     const F energy_eV = hc_eVnm/wavelength_nm ;
0391     get_lpmtid_stackspec( spec_16, lpmtid, energy_eV );
0392 }
0393
0394 template<typename F>
0395 inline QPMT_METHOD void qpmt<F>::get_lpmtid_SPEC_ce(
0396      F* spec_16
0397    , int lpmtid
0398    , F wavelength_nm
0399    , F lposcost
0400 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0401    , unsigned pidx
0402    , bool pidx_debug
0403 #endif
0404    ) const
0405 {
0406     const F energy_eV = hc_eVnm/wavelength_nm ;
0407 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0408     get_lpmtid_stackspec_ce_acosf( spec_16, lpmtid, energy_eV, lposcost, pidx, pidx_debug );
0409 #else
0410     get_lpmtid_stackspec_ce_acosf( spec_16, lpmtid, energy_eV, lposcost );
0411 #endif
0412 }
0413
0414
0415
0416 /**
0417 qpmt::get_lpmtid_LL
0418 ----------------------
0419
0420 Full details of all the TMM layers for debugging.
0421
0422 **/
0423
0424 template<typename F>
0425 inline QPMT_METHOD void qpmt<F>::get_lpmtid_LL(
0426     F* ll_128,
0427     int lpmtid,
0428     F wavelength_nm,
0429     F minus_cos_theta,
0430     F dot_pol_cross_mom_nrm ) const
0431 {
0432     const F energy_eV = hc_eVnm/wavelength_nm ;
0433
0434     F spec[16] ;
0435     get_lpmtid_stackspec( spec, lpmtid, energy_eV );
0436
0437     Stack<F,4> stack(wavelength_nm, minus_cos_theta, dot_pol_cross_mom_nrm, spec, 16u );
0438     const F* stack_ll = stack.ll[0].cdata() ;
0439
0440     for(int i=0 ; i < 128 ; i++ ) ll_128[i] = stack_ll[i] ;
0441 }
0442
0443 /**
0444 qpmt::get_lpmtid_COMP
0445 ----------------------
0446
0447 Full details of the composite TMM layer for debugging.
0448
0449 **/
0450
0451 template<typename F>
0452 inline QPMT_METHOD void qpmt<F>::get_lpmtid_COMP(
0453     F* comp_32,
0454     int lpmtid,
0455     F wavelength_nm,
0456     F minus_cos_theta,
0457     F dot_pol_cross_mom_nrm ) const
0458 {
0459     const F energy_eV = hc_eVnm/wavelength_nm ;
0460
0461     F spec[16] ;
0462     get_lpmtid_stackspec( spec, lpmtid, energy_eV );
0463
0464     Stack<F,4> stack(wavelength_nm, minus_cos_theta, dot_pol_cross_mom_nrm, spec, 16u );
0465     const F* stack_comp = stack.comp.cdata() ;
0466     for(int i=0 ; i < 32 ; i++ ) comp_32[i] = stack_comp[i] ;
0467 }
0468
0469 /**
0470 qpmt::get_lpmtid_ART
0471 ----------------------
0472
0473 4x4 ART matrix with A,R,T for S/P/av/actual polarizations
0474
0475 **/
0476
0477
0478 template<typename F>
0479 inline QPMT_METHOD void qpmt<F>::get_lpmtid_ART(
0480     F* art_16,
0481     int lpmtid,
0482     F wavelength_nm,
0483     F minus_cos_theta,
0484     F dot_pol_cross_mom_nrm ) const
0485 {
0486     const F energy_eV = hc_eVnm/wavelength_nm ;
0487
0488     F spec[16] ;
0489     get_lpmtid_stackspec( spec, lpmtid, energy_eV );
0490
0491     Stack<F,4> stack(wavelength_nm, minus_cos_theta, dot_pol_cross_mom_nrm, spec, 16u );
0492     const F* stack_art = stack.art.cdata() ;
0493
0494     for(int i=0 ; i < 16 ; i++ ) art_16[i] = stack_art[i] ;
0495 }
0496
0497 /**
0498 qpmt::get_lpmtid_ARTE
0499 -----------------------
0500
0501 lpmtid and polarization customized TMM calc of::
0502
0503    A : theAbsorption
0504    R : theReflectivity
0505    T : theTransmittance
0506    E : theEfficiency (more specifically QE)
0507
0508 Q: Does theReflectivity+theTransmittace = 1.f ?
0509 A: YES, by construction because A+R+T = one (triplet prob)
0510
0511    so : R+T=one-A
0512
0513    R+T
0514    ----- = one
0515    one-A
0516
0517
0518 **/
0519
0520 template<typename F>
0521 inline QPMT_METHOD void qpmt<F>::get_lpmtid_ARTE(
0522     F* ARTE,
0523     int lpmtid,
0524     F wavelength_nm,
0525     F minus_cos_theta,
0526     F dot_pol_cross_mom_nrm ) const
0527 {
0528     const F energy_eV = hc_eVnm/wavelength_nm ;
0529
0530     F spec[16] ;
0531     get_lpmtid_stackspec( spec, lpmtid, energy_eV );
0532
0533     const F* ss = spec ;
0534     const F& _qe = spec[15] ;
0535
0536 #ifdef MOCK_CURAND_DEBUG
0537     printf("//qpmt.get_lpmtid_ARTE lpmtid %d energy_eV %7.3f _qe %7.3f \n", lpmtid, energy_eV, _qe );
0538 #endif
0539
0540
0541     Stack<F,4> stack ;
0542
0543     if( minus_cos_theta < zero )
0544     {
0545         stack.calc(wavelength_nm, -one, zero, ss, 16u );
0546         ARTE[3] = _qe/stack.art.A ;
0547
0548 #ifdef MOCK_CURAND_DEBUG
0549         printf("//qpmt.get_lpmtid_ARTE stack.art.A %7.3f _qe/stack.art.A %7.3f \n", stack.art.A, ARTE[3] );
0550 #endif
0551     }
0552     else
0553     {
0554         ARTE[3] = zero ;
0555     }
0556
0557     stack.calc(wavelength_nm, minus_cos_theta, dot_pol_cross_mom_nrm, ss, 16u );
0558
0559     const F& A = stack.art.A ;
0560     const F& R = stack.art.R ;
0561     const F& T = stack.art.T ;
0562
0563     ARTE[0] = A ;         // aka theAbsorption
0564     ARTE[1] = R/(one-A) ; // aka theReflectivity
0565     ARTE[2] = T/(one-A) ; // aka theTransmittance
0566
0567 #ifdef MOCK_CURAND_DEBUG
0568     std::cout << "//qpmt.get_lpmtid_ARTE stack.calc.2 " << std::endl ;
0569     std::cout << "//qpmt.get_lpmtid_ARTE wavelength_nm " << wavelength_nm << std::endl ;
0570     std::cout << "//qpmt.get_lpmtid_ARTE minus_cos_theta " << minus_cos_theta << std::endl ;
0571     std::cout << "//qpmt.get_lpmtid_ARTE dot_pol_cross_mom_nrm " << dot_pol_cross_mom_nrm << std::endl ;
0572     std::cout << "//qpmt.get_lpmtid_ARTE " << stack << std::endl ;
0573 #endif
0574
0575 }
0576
0577
0578 /**
0579 qpmt::get_lpmtid_ATQC
0580 ------------------------
0581
0582 Invoked from qsim::propagate_at_surface_CustomART
0583
0584
0585 ATQC[0] A
0586    absorption
0587
0588 ATQC[1] T
0589    transmission, scaled by 1/(1-A) to make R+T = 1
0590
0591 ATQC[2] Q
0592    QE quantum efficiency, depending on photon energy
0593    and PMT TMM parameters
0594
0595 ATQC[3] C
0596    CE collection efficiency, depending on theta of local position on PMT surface
0597    obtained by interpolation over theta domain (OR maybe in future costheta domain)
0598
0599
0600 TODO: compare between the alternates::
0601
0602     get_lpmtid_stackspec_ce_acosf   // interpolates ce in theta of local position in PMT frame
0603     get_lpmtid_stackspec_ce         // interpolates ce in costheta of local position in PMT frame
0604
0605 **/
0606
0607 template<typename F>
0608 inline QPMT_METHOD void qpmt<F>::get_lpmtid_ATQC(
0609       F* ATQC
0610     , int lpmtid
0611     , F wavelength_nm
0612     , F minus_cos_theta
0613     , F dot_pol_cross_mom_nrm
0614     , F lposcost
0615 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0616     , unsigned pidx
0617     , bool pidx_debug
0618 #endif
0619     ) const
0620 {
0621     const F energy_eV = hc_eVnm/wavelength_nm ;
0622
0623     F spec[16] ;
0624
0625
0626 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0627     get_lpmtid_stackspec_ce_acosf( spec, lpmtid, energy_eV, lposcost, pidx, pidx_debug );
0628 #else
0629     get_lpmtid_stackspec_ce_acosf( spec, lpmtid, energy_eV, lposcost );
0630 #endif
0631
0632
0633     const F* ss = spec ;
0634     const F& ce = spec[7] ;
0635     ATQC[3] = ce ;
0636
0637     const F& _qe = spec[15] ;
0638
0639 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0640     //if(pidx_debug) printf("//qpmt.get_lpmtid_ATQC pidx %7d lpmtid %d energy_eV %8.4f _qe(spec[15]) %8.4f lposcost %8.4f ce %8.4f  \n", pidx, lpmtid, energy_eV, _qe, lposcost, ce );
0641 #endif
0642
0643     Stack<F,4> stack ;
0644
0645     if( minus_cos_theta < zero )
0646     {
0647         stack.calc(wavelength_nm, -one, zero, ss, 16u );
0648         ATQC[2] = _qe/stack.art.A ;
0649
0650 #if !defined(PRODUCTION) && defined(DEBUG_PIDX)
0651         if(pidx_debug) printf("//qpmt.get_lpmtid_ATQC pidx %7d lpmtid %d energy_eV %8.4f _qe(spec[15]) %8.4f lposcost %8.4f ce %8.4f stack.art.A %8.4f _qe/stack.art.A=ATQC[2] %8.4f  ce=ATQC[3] %8.4f  \n",
0652                                                       pidx, lpmtid, energy_eV, _qe, lposcost, ce, stack.art.A, ATQC[2], ATQC[3] );
0653 #endif
0654     }
0655     else
0656     {
0657         ATQC[2] = zero ;
0658     }
0659
0660     stack.calc(wavelength_nm, minus_cos_theta, dot_pol_cross_mom_nrm, ss, 16u );
0661
0662     const F& A = stack.art.A ;
0663     const F& T = stack.art.T ;
0664
0665     ATQC[0] = A ;         // aka theAbsorption
0666     ATQC[1] = T/(one-A) ; // aka theTransmittance
0667
0668
0669 #ifdef MOCK_CURAND_DEBUG
0670     std::cout << "//qpmt.get_lpmtid_ATQC stack.calc.2 " << std::endl ;
0671     std::cout << "//qpmt.get_lpmtid_ATQC wavelength_nm " << wavelength_nm << std::endl ;
0672     std::cout << "//qpmt.get_lpmtid_ATQC minus_cos_theta " << minus_cos_theta << std::endl ;
0673     std::cout << "//qpmt.get_lpmtid_ATQC dot_pol_cross_mom_nrm " << dot_pol_cross_mom_nrm << std::endl ;
0674     std::cout << "//qpmt.get_lpmtid_ATQC " << stack << std::endl ;
0675 #endif
0676
0677 }
0678
0679
0680
0681
0682
0683 // end WITH_CUSTOM4
0684 #endif
0685
0686 // end CUDA/MOCK_CUDA
0687 #endif
0688
0689
0690 #if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_CURAND ) || defined(MOCK_CUDA)
0691 #else
0692 template struct QUDARAP_API qpmt<float>;
0693 //template struct QUDARAP_API qpmt<double>;
0694 #endif
0695
0696
0697