Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-09 07:49:08 

0001 #include <sstream>
0002 #include <algorithm>
0003 #include <csignal>
0004 
0005 #include "SLOG.hh"
0006 #include "spath.h"
0007 #include "scuda.h"
0008 
0009 #include "NP.hh"
0010 
0011 #include "QUDA_CHECK.h"
0012 #include "QRng.hh"
0013 #include "QCerenkovIntegral.hh"
0014 
0015 #include "QCK.hh"
0016 
0017 
0018 const plog::Severity QCerenkovIntegral::LEVEL = SLOG::EnvLevel("QCerenkovIntegral", "DEBUG"); 
0019 
0020 
0021 const unsigned QCerenkovIntegral::SPLITBIN_PAYLOAD_SIZE = 8 ; 
0022 const unsigned QCerenkovIntegral::UPPERCUT_PAYLOAD_SIZE = 3 ; 
0023 
0024 
0025 
0026 
0027 
0028 /**
0029 QCerenkovIntegral::QCerenkovIntegral
0030 ------------------------------------
0031 
0032 Requires path to RINDEX.npy 
0033 Note this is just handling rindex for a single material. 
0034 
0035 
0036 **/
0037 
0038 QCerenkovIntegral::QCerenkovIntegral(const char* path_ )
0039     :
0040     path( spath::Resolve(path_) ),
0041     dsrc( NP::LoadIfExists(path)),
0042     src( nullptr ),
0043     emn(0.),
0044     emx(0.),
0045     rmn(0.),
0046     rmx(0.)
0047 {
0048     init(); 
0049 }
0050 
0051 
0052 void QCerenkovIntegral::init()
0053 {
0054     LOG_IF(error, dsrc == nullptr ) 
0055         << " FAILED to load " 
0056         << " path " << path 
0057         ; 
0058     if(dsrc == nullptr) return ; 
0059 
0060     if(dsrc) src = dsrc->ebyte == 4 ? dsrc : NP::MakeNarrow(dsrc) ; 
0061 
0062     dsrc->pscale<double>(1e6, 0u) ; //  change energy scale from MeV to eV,   1.55 to 15.5 eV
0063     dsrc->minmax<double>(emn, emx, 0u ); 
0064     dsrc->minmax<double>(rmn, rmx, 1u ); 
0065 
0066     LOG(info) 
0067         << " emn " << std::setw(10) << std::fixed << std::setprecision(4) << emn 
0068         << " emx " << std::setw(10) << std::fixed << std::setprecision(4) << emx 
0069         << " rmn " << std::setw(10) << std::fixed << std::setprecision(4) << rmn 
0070         << " rmx " << std::setw(10) << std::fixed << std::setprecision(4) << rmx 
0071         ;
0072 
0073     // cannot makeTex here as the icdf is not yet formed 
0074     // it takes a lot of integrals and CDF inversion to form it 
0075 }
0076 
0077 /**
0078 QCerenkovIntegral::GetAverageNumberOfPhotons_s2
0079 ----------------------------------------
0080 
0081 This was prototyped and tested in::
0082 
0083     ~/opticks/ana/ckn.py 
0084     ~/opticks/examples/Geant4/CerenkovStandalone/G4Cerenkov_modifiedTest.cc
0085 
0086 
0087 "charge" argument is typically 1. or -1, units of Rfact are (eV cm)^-1 so energies must be in eV 
0088 to give Cerenkov photon yield per mm. Essentially Rfact is : "2*pi*fine_structure_constant/hc", approx::
0089 
0090     In [1]: 2.*np.pi*(1./137.)*(1./1240.)*1e6  # hc ~ 1240 eV nm, nm = 1e-6 mm, fine_structure_const ~ 1/137.   
0091     Out[1]: 36.9860213514221
0092 
0093 To understand the more precise Rfact value see ~/opticks/examples/UseGeant4/UseGeant4.cc and google "Frank Tamm Formula" 
0094 
0095 **/
0096 
0097 template <typename T>
0098 NP* QCerenkovIntegral::GetAverageNumberOfPhotons_s2_(T& emin,  T& emax, const T BetaInverse, const T  charge ) const 
0099 {
0100     emin = emx ; // start with inverted range
0101     emax = emn ; 
0102  
0103     const T* vv = dsrc->cvalues<T>(); 
0104     unsigned ni = dsrc->shape[0] ; 
0105     unsigned nj = dsrc->shape[1] ; 
0106 
0107     bool dsrc_expect = nj == 2 && ni > 1 ;
0108     assert(dsrc_expect); 
0109     if(!dsrc_expect) std::raise(SIGINT); 
0110 
0111     T s2integral(0.) ;  
0112 
0113     T en_a = -1. ; 
0114     T ri_a = -1. ; 
0115     T en_b = -1. ; 
0116     T ri_b = -1. ; 
0117 
0118     NP* s2i = NP::Make<T>(ni) ;   // hmm: top value will always be zero 
0119     T* s2i_vv = s2i->values<T>(); 
0120     bool dump = false ; 
0121 
0122     for(unsigned i=0 ; i < ni - 1 ; i++)
0123     {
0124         T en_0 = vv[2*(i+0)+0] ; 
0125         T en_1 = vv[2*(i+1)+0] ; 
0126 
0127         T ri_0 = vv[2*(i+0)+1] ; 
0128         T ri_1 = vv[2*(i+1)+1] ; 
0129 
0130         T bin_integral = charge*charge*GetS2Integral_WithCut<T>( emin, emax, BetaInverse, en_0, en_1, ri_0, ri_1, en_a, ri_a, en_b, ri_b, dump );
0131 
0132         s2i_vv[i] = bin_integral ; 
0133         s2integral += bin_integral ; 
0134     }
0135     return s2i ; 
0136 }
0137 template NP* QCerenkovIntegral::GetAverageNumberOfPhotons_s2_<double>(double& emin,  double& emax, const double BetaInverse, const double  charge ) const ;
0138 template NP* QCerenkovIntegral::GetAverageNumberOfPhotons_s2_<float>( float& emin,   float&  emax, const float  BetaInverse, const float   charge ) const ;
0139 
0140 
0141 
0142 template <typename T>
0143 T QCerenkovIntegral::GetAverageNumberOfPhotons_s2(T& emin,  T& emax, const T BetaInverse, const T  charge ) const 
0144 {
0145     NP* s2i = GetAverageNumberOfPhotons_s2_<T>(emin, emax, BetaInverse, charge ); 
0146     T s2integral = s2i->psum<T>(0u); 
0147     return s2integral ; 
0148 }
0149 
0150 template double QCerenkovIntegral::GetAverageNumberOfPhotons_s2<double>(double& emin,  double& emax, const double BetaInverse, const double  charge ) const ;
0151 template float  QCerenkovIntegral::GetAverageNumberOfPhotons_s2<float>( float& emin,   float&  emax, const float  BetaInverse, const float   charge ) const ;
0152 
0153 
0154 
0155 template <typename T>
0156 NP* QCerenkovIntegral::getAverageNumberOfPhotons_s2( const NP* bis ) const   
0157 {
0158     T emin ; 
0159     T emax ; 
0160     T charge = T(1.) ;  
0161 
0162     unsigned ni = bis->shape[0] ; 
0163     NP* avph = NP::Make<T>(ni, 4 ); 
0164     T* avph_v = avph->values<T>(); 
0165 
0166     for(unsigned i=0 ; i < ni ; i++)
0167     {
0168         const T BetaInverse = bis->get<T>(i) ; 
0169         T avgNumPhotons = GetAverageNumberOfPhotons_s2<T>(emin, emax, BetaInverse, charge ); 
0170 
0171         avph_v[4*i+0] = BetaInverse ; 
0172         avph_v[4*i+1] = emin  ; 
0173         avph_v[4*i+2] = emax ; 
0174         avph_v[4*i+3] = avgNumPhotons ; 
0175     }
0176     return avph ; 
0177 }
0178 
0179 template NP* QCerenkovIntegral::getAverageNumberOfPhotons_s2<float>( const NP* ) const ; 
0180 template NP* QCerenkovIntegral::getAverageNumberOfPhotons_s2<double>( const NP* ) const ; 
0181 
0182 
0183 
0184 
0185 
0186 
0187 
0188 
0189 
0190 
0191 const double QCerenkovIntegral::FINE_STRUCTURE_OVER_HBARC_EVMM = 36.981 ; 
0192 
0193 const char* QCerenkovIntegral::NONE_ = "NONE" ; 
0194 const char* QCerenkovIntegral::UNCUT_ = "UNCUT" ; 
0195 const char* QCerenkovIntegral::SUB_ = "SUB" ; 
0196 const char* QCerenkovIntegral::CUT_ = "CUT" ; 
0197 const char* QCerenkovIntegral::FULL_ = "FULL" ; 
0198 const char* QCerenkovIntegral::PART_ = "PART" ; 
0199 const char* QCerenkovIntegral::ERR_ = "ERR" ; 
0200 
0201 const char* QCerenkovIntegral::State(int state)
0202 {
0203     const char* s = nullptr ; 
0204     switch(state)
0205     {
0206        case NONE:  s = NONE_  ; break ; 
0207        case SUB:   s = SUB_   ; break ; 
0208        case UNCUT: s = UNCUT_ ; break ; 
0209        case CUT:   s = CUT_   ; break ; 
0210        case FULL:  s = FULL_  ; break ; 
0211        case PART:  s = PART_  ; break ; 
0212        case ERR:   s = ERR_   ; break ; 
0213     }
0214     return s ; 
0215 }
0216 
0217 
0218 
0219 /**
0220 QCerenkovIntegral::GetS2Integral_WithCut
0221 ------------------------------------------
0222 
0223 Whether the en_cut and ri_cut are used depends on the en_cut value compared to en_0 and en_1 
0224 
0225 When integrating trapezoids (like B,C,D) cut bin partial integrals
0226 are simply done by replacing en_1 with the somewhat smaller en_cut.
0227 Integrating edge triangle integrals need more care as the result depends 
0228 on how en_cut compares with en_cross.::
0229 
0230                                  *
0231                                 /|\
0232                                / | \
0233                               /  |  \
0234                              /   |   \
0235                             /    |    \
0236              . .           /     |     \   . .
0237              . .          /      |      \  . .
0238              . . *-------*       |       * . .
0239              . ./|       |       |       |\. .
0240              . / |       |       |       | \ .
0241              ./.A|  B    |   C   |   D   |E.\.
0242      -----0--x-.-1-------2-------3-------4-.-y----5----------
0243              . .                           . .
0244              . en_b                        . en_cross
0245              .                             .          
0246              en_cross                      en_b   
0247 
0248 
0249     full A:  (en_1 - en_cross)*s2_1*half
0250 
0251     part A:  0.                                  en_cut <= en_cross 
0252     part A:  (en_cut - en_cross)*s2_cut*half 
0253 
0254     full C:  (en_1 - en_0)*(s2_0 + s2_1)*half 
0255     part C:  (en_cut - en_0)*(s2_0 + s2_cut)*half
0256 
0257     full E:  (en_cross - en_0)*s2_0*half
0258     part E:  (en_cross - en_0)*s2_0*half           en_cut >= en_cross
0259     part E:  (en_cut - en_0)*(s2_0+s2_cut)*half    en_cut < en_cross       *rhs triangle becomes trapezoid* 
0260 
0261 * en_b (formerly known as en_cut) is used for multi-bin integrals up to en_b(en_cut)
0262 * en_b turns triangle E into a trapezoid for en_b < en_cross 
0263 * en_b can only turn triangle A into a smaller triangle  
0264 
0265 * subsequently the en_a argument was added to provide SUB-bin integrals between en_a and en_b 
0266   which are both required to be within en_0 en_1 of a single bin 
0267 
0268 
0269 
0270               en_0                 en_1 
0271                +--------------------+
0272                |                    |
0273          bbbbbb|                    |                     ( en_0 > en_b  && en_1 > en_b )     -> CUT     bin above the cut 
0274                |                    |
0275                |                    |bbbbb                ( en_0 < en_b  && en_1 <= en_b )    -> FULL    bin below the cut                      
0276                |                    |
0277                bbbbbbbbbbbbbbbbbbbbb|                     ( en_0 <= en_b && en_b < en_1 )     -> PART           
0278                |                    |
0279                aaaaaaaabbbbbbbbbbbbbb                     ( en_a >= en_0 && en_b <= en_1 )     -> SUB 
0280                |                    |
0281                |                    |
0282                |                    |
0283                +--------------------+
0284 
0285 **/
0286 
0287 template<typename T>
0288 T QCerenkovIntegral::GetS2Integral_WithCut( T& emin, T& emax, T BetaInverse, T en_0, T en_1 , T ri_0, T ri_1, T en_a, T ri_a, T en_b, T ri_b, bool dump ) // static 
0289 {
0290 
0291     int state = NONE ; 
0292     if(      en_b <= 0.  )                   state = UNCUT ;   
0293     else if( en_a >= en_0 && en_b <= en_1 )  state = SUB ;   // SUB-bin integral with sub-range en_a:en_b entirely within(including edge) range of big bin en_0:en_1
0294     else if( en_0 >  en_b && en_1 >  en_b )  state = CUT ;   
0295     else if( en_0 <  en_b && en_1 <= en_b )  state = FULL ;  // edges of this bin both below en_b (or en_1 == en_b), so full bin integral 
0296     else if( en_0 <= en_b && en_b <  en_1 )  state = PART ;  // edges of this bin straddle ecut, so partial bin integral  
0297     else                                     state = ERR ; 
0298    
0299     bool missed_condition = state == ERR || state == NONE ; 
0300  
0301     LOG_IF(error, missed_condition) << " missed condition ?" << " en_0 " << en_0 << " en_1 " << en_1 << " en_a " << en_a << " en_b " << en_b ;  
0302  
0303     const T zero(0.) ; 
0304     if( state == CUT ) return zero ; 
0305     assert( state == FULL || state == PART || state == UNCUT || state == SUB ); 
0306     if(state == SUB) assert( en_a >= en_0 && en_b <= en_1 ); 
0307 
0308     const T one(1.) ; 
0309     const T half(0.5) ; 
0310     T s2integral(0.) ;  
0311 
0312     T ct_0 = BetaInverse/ri_0 ; 
0313     T ct_1 = BetaInverse/ri_1 ; 
0314     T ct_a = en_a > 0 ? BetaInverse/ri_a : -1 ; 
0315     T ct_b = en_b > 0 ? BetaInverse/ri_b : -1 ; 
0316     T s2_0 = ( one - ct_0 )*( one + ct_0 ); 
0317     T s2_1 = ( one - ct_1 )*( one + ct_1 );
0318     T s2_a = en_a > 0 ? ( one - ct_a )*( one + ct_a ) : -1 ; 
0319     T s2_b = en_b > 0 ? ( one - ct_b )*( one + ct_b ) : -1 ; 
0320 
0321     if(dump) std::cout 
0322         << "QCerenkovIntegral::GetS2Integral_WithCut"
0323         << " en_0 " << en_0 
0324         << " en_1 " << en_1 
0325         << " en_a " << en_a 
0326         << " en_b " << en_b 
0327         << " s2_0 " << s2_0 
0328         << " s2_1 " << s2_1 
0329         << " s2_a " << s2_a
0330         << " s2_b " << s2_b
0331         << " state " << State(state)
0332         << " BetaInverse " << BetaInverse
0333         << std::endl
0334         ;
0335 
0336 
0337     bool cross = s2_0*s2_1 < zero ; 
0338     T en_cross = cross ? en_0 + (BetaInverse - ri_0)*(en_1 - en_0)/(ri_1 - ri_0) : -one  ;
0339 
0340 
0341     if( s2_0 < zero && s2_1 > zero )  // s2 becomes +ve within the bin, eg bin A : lhs triangle 
0342     {    
0343         assert( en_cross > zero ); 
0344 
0345         if( state == FULL || state == UNCUT )
0346         {  
0347             s2integral =  (en_1 - en_cross)*s2_1*half ;     // left triangle 
0348 
0349             emin = std::min<T>( emin, en_cross );  
0350             emax = std::max<T>( emax, en_1 );
0351         }
0352         else if (state == PART )
0353         {
0354             s2integral = en_b <= en_cross ? zero : (en_b - en_cross)*s2_b*half ; 
0355 
0356             emin = std::min<T>( emin, en_b <= en_cross ? emin : en_cross );  
0357             emax = std::max<T>( emax, emax );
0358         }
0359         else if (state == SUB )
0360         {
0361             /**
0362              Three possibilites for en_a and en_b wrt en_cross : both below, straddling, both above 
0363                                   /
0364                  en_a     en_b   / 
0365                +---+-------+----*---------------+         zero        
0366               en_0            en_cross         en_1
0367             
0368                                   /
0369                           en_a   /  en_b                  triangle from en_cross to en_b 
0370                +-----------+----*----+----------+                
0371               en_0             en_cross        en_1
0372             
0373                                   /
0374                                  /  en_a  en_b            trapezoid independent of en_cross
0375                +----------------*----+-----+----+                
0376               en_0           en_cross          en_1
0377             
0378             **/
0379 
0380             if( en_a <= en_cross && en_b <= en_cross )
0381             {
0382                  s2integral = zero ; 
0383             }
0384             else
0385             {
0386                  if( en_a <= en_cross && en_b >= en_cross )
0387                  {
0388                      s2integral = (en_b - en_cross)*s2_b*half ; 
0389                      emin = std::min<T>( emin, en_cross ) ; 
0390                      emax = std::min<T>( emax, en_b ) ; 
0391                  }
0392                  else
0393                  {
0394                      s2integral = (en_b - en_a)*(s2_a + s2_b)*half ; 
0395                      emin = std::min<T>( emin, en_a ) ; 
0396                      emax = std::min<T>( emax, en_b ) ; 
0397                  }
0398             }
0399         }
0400     }    
0401     else if( s2_0 >= zero && s2_1 >= zero )   // s2 +ve across full bin 
0402     {    
0403         if( state == FULL || state == UNCUT )
0404         {  
0405             s2integral = (en_1 - en_0)*(s2_0 + s2_1)*half ;  
0406 
0407             emin = std::min<T>( emin, en_0 );  
0408             emax = std::max<T>( emax, en_1 );
0409         }
0410         else if( state == PART )
0411         {
0412             s2integral = (en_b - en_0)*(s2_0 + s2_b)*half ;  
0413             if(dump) std::cout << " s2 +ve across full bin " << " s2integral " << s2integral << std::endl ; 
0414             emin = std::min<T>( emin, en_0 );  
0415             emax = std::max<T>( emax, en_b );
0416         }
0417         else if (state == SUB )
0418         {
0419             s2integral = (en_b - en_a)*(s2_a + s2_b)*half ;  
0420             emin = std::min<T>( emin, en_a );  
0421             emax = std::max<T>( emax, en_b );
0422         }
0423     }     
0424     else if( s2_0 > zero && s2_1 < zero )  // s2 becomes -ve within the bin, eg bin E :  rhs triangle 
0425     {    
0426         assert( en_cross > zero ); 
0427 
0428         if( state == FULL || state == UNCUT )
0429         {  
0430             s2integral =  (en_cross - en_0)*s2_0*half ;  
0431 
0432             emin = std::min<T>( emin, en_0 );  
0433             emax = std::max<T>( emax, en_cross );
0434 
0435         }
0436         else if( state == PART )
0437         {
0438             s2integral = en_b >= en_cross ? (en_cross - en_0)*s2_0*half : (en_b - en_0)*(s2_0+s2_b)*half ;   
0439 
0440             emin = std::min<T>( emin, en_0 );  
0441             emax = std::max<T>( emax, en_b >= en_cross ? en_cross : en_b );
0442         }
0443         else if (state == SUB )
0444         {
0445             /**
0446              Three possibilites for en_a and en_b wrt en_cross : both below, straddling, both above 
0447 
0448                              \
0449                               \
0450                  en_a     en_b \
0451                +---+-------+----*---------------+         trapezoid independent of en_cross        
0452               en_0            en_cross         en_1
0453             
0454                              \
0455                               \
0456                           en_a \    en_b                  triangle from en_a to en_cross 
0457                +-----------+----*----+----------+                
0458               en_0             en_cross        en_1
0459             
0460                              \
0461                               \
0462                                \    en_a  en_b            zero
0463                +----------------*----+-----+----+                
0464               en_0           en_cross          en_1
0465             
0466             **/
0467             if( en_a <= en_cross && en_b <= en_cross )
0468             {
0469                  s2integral = (en_b - en_a)*(s2_a + s2_b)*half ; 
0470                  emin = std::min<T>( emin, en_a ) ; 
0471                  emax = std::min<T>( emax, en_b ) ; 
0472             }
0473             else
0474             {
0475                  if( en_a <= en_cross && en_b >= en_cross )
0476                  {
0477                      s2integral = (en_cross - en_a)*s2_a*half ; 
0478                      emin = std::min<T>( emin, en_a ) ; 
0479                      emax = std::min<T>( emax, en_cross ) ; 
0480                  }
0481                  else
0482                  {
0483                      s2integral = zero ; 
0484                  }
0485             }
0486         }
0487     }    
0488     const T Rfact = FINE_STRUCTURE_OVER_HBARC_EVMM ;  
0489     return Rfact*s2integral ; 
0490 }
0491 
0492 template float QCerenkovIntegral::GetS2Integral_WithCut( float& , float&,    float,  float,  float,  float,  float,  float,  float, float, float, bool  ); 
0493 template double QCerenkovIntegral::GetS2Integral_WithCut( double&, double&,  double,  double,  double,  double,  double,  double,  double, double, double, bool  ); 
0494 
0495 
0496 /**
0497 QCerenkovIntegral::getS2Integral_WithCut_
0498 ------------------------------------
0499 
0500 For debugging purposes the contributions to the integral from each rindex bin are returned in an array.
0501 
0502 **/
0503 
0504 template <typename T> NP* QCerenkovIntegral::getS2Integral_WithCut_( T& emin, T& emax, T BetaInverse, T en_a, T en_b, bool dump ) const 
0505 {
0506     const T* vv = dsrc->cvalues<T>(); 
0507     unsigned ni = dsrc->shape[0] ; 
0508     unsigned nj = dsrc->shape[1] ; 
0509     bool dsrc_expect = nj == 2 && ni > 1  ;
0510     assert(dsrc_expect ); 
0511     if(!dsrc_expect) std::raise(SIGINT); 
0512 
0513     NP* s2i = NP::Make<T>(ni);
0514     T* s2i_vv = s2i->values<T>();  
0515 
0516     T ri_a = en_a > 0 ? dsrc->interp<T>(en_a) : -1 ; 
0517     T ri_b = en_b > 0 ? dsrc->interp<T>(en_b) : -1 ; 
0518 
0519     for(unsigned i=0 ; i < ni - 1 ; i++)
0520     {
0521         T en_0 = vv[2*(i+0)+0] ; 
0522         T en_1 = vv[2*(i+1)+0] ; 
0523 
0524         T ri_0 = vv[2*(i+0)+1] ; 
0525         T ri_1 = vv[2*(i+1)+1] ; 
0526 
0527         s2i_vv[i] = GetS2Integral_WithCut<T>(emin, emax, BetaInverse, en_0, en_1, ri_0, ri_1, en_a, ri_a, en_b, ri_b, dump );
0528     }
0529     return s2i ; 
0530 }
0531  
0532 template NP* QCerenkovIntegral::getS2Integral_WithCut_( double&, double&, double, double, double, bool  ) const ; 
0533 template NP* QCerenkovIntegral::getS2Integral_WithCut_( float&,  float&,  float,  float,  float,  bool  ) const ; 
0534 
0535 
0536 
0537 template <typename T>
0538 T QCerenkovIntegral::getS2Integral_WithCut( T& emin, T& emax, T BetaInverse, T en_a, T en_b, bool dump ) const 
0539 {
0540     NP* s2i = getS2Integral_WithCut_<T>(emin, emax, BetaInverse, en_a, en_b, dump ); 
0541     T s2integral = s2i->psum<T>(0u); 
0542     return s2integral ; 
0543 }  
0544 template double QCerenkovIntegral::getS2Integral_WithCut( double&, double&, double, double, double, bool  ) const ; 
0545 template float  QCerenkovIntegral::getS2Integral_WithCut( float&,  float&,  float, float, float, bool   ) const ; 
0546 
0547 
0548 /**
0549 QCerenkovIntegral::getS2Integral_SplitBin
0550 -------------------------------------
0551 
0552 Returns *s2c* array of shape (s2_edges, PAYLOAD_SIZE(=8)) 
0553 with the last payload entry being s2integral.
0554 
0555 The s2 value of evaluated at en_cut and s2integral 
0556 is the cumulative integral up to en_cut.  The first 
0557 value of s2integral is zero. 
0558 
0559 * this approach leads to lots of zero bins because of the fixed full energy range 
0560 * would be simple to restrict to big bins with contributions allowing mul to be increased at the same cost 
0561   by greatly reducing zeros
0562 
0563 
0564 Why smaller sub-bins are needed ?
0565 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0566 
0567 Objective is to effectively provide an inverse-CDF such that
0568 a random number lookup yields an energy. 
0569 As s2 is piecewise linear with energy(from the piecewise linear RINDEX) 
0570 the cumulative integral will be piecewise parabolic. 
0571 Hence linear interpolation on that parabolic is not going to be a very 
0572 good approximation, so want to make the binning smaller than the
0573 rindex bins. 
0574 
0575 Simplification : by bin-splitting
0576 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0577 
0578 Avoid complexity and imprecision arising from differences between 
0579 rindex energy edges and the smaller sub bins by adopting approach of ana/edges.py:divide_bins
0580 and more directly NP::MakeDiv 
0581 
0582 * there is no need for the sub-bins to all be equally spaced, the motivation for 
0583   the sub-bins is to profit from the piecewise linear nature of s2 which means that 
0584   s2 can be linearly interpolated with very small errors unlike the parabolic
0585   cumulative integral for which linear interpolation is a poor approximation 
0586 
0587 * splitting the edges into sub-bins allows to simply loop within each big rindex 
0588   bin saving the partials off into the cumulative integral array   
0589 
0590 
0591 not easy to avoid lots of s2c zero bins in this _SplitBin approach
0592 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0593 
0594 * with _UpperCut the number of energy bins is imposed to 
0595   be equal for all BetaInverse and different Cerenkov permissable 
0596   energy ranges are used for each BetaInverse
0597 
0598 * to do something similar here in _SplitBin would need to have the 
0599   number of energy bins different for each BetaInverse, 
0600   for example including only sub-bins within big-bins with contributions
0601   
0602 * would need to NP::Combine ragged individual BetaInverse s2c 
0603 * complexity seems not worth it, especially as s2c/s2cn are just 
0604   intermediaries on way to creating fixed 0:1 domain icdf array 
0605   that is the input to creating the GPU texture 
0606 
0607 * that icdf is inherently regular fixed domain as created by NP::pdomain 
0608   inversion of s2cn across 0->1 range  
0609 
0610 
0611 prev/next cumulative approach of Geant4
0612 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0613 
0614 Prev/next approach commonly used in Geant4 BuildThePhysicsTable 
0615 for cumulative integration could slightly reduce computation as en_1 ri_1 s2_1 
0616 becomes en_0 ri_0 s2_0 for next bin. 
0617 
0618 Maybe that helps avoid non-monotonic ?
0619 
0620 
0621 See Also
0622 ~~~~~~~~~~~~
0623 
0624 * ana/piecewise.py for some whacky parabolic Ideas
0625 
0626 TODO:
0627 
0628 * create the icdf float4 payload array, by baking NP::pdomain results from 0->1 for each BetaInverse 
0629 * check lookups using it still OK chi2 with sampling 
0630 * turn into GPU texture and check again, will probably need the hd_factor trick from 
0631   X4Scintillation::CreateGeant4InterpolatedInverseCDF in the energy extremes where the CDF 
0632   tends to get very flat (and hence ICDF is very steep)
0633   
0634 
0635 **/
0636 
0637 
0638 template <typename T>
0639 NP* QCerenkovIntegral::getS2Integral_SplitBin(const T BetaInverse, unsigned mul, bool dump ) const 
0640 {
0641     T emin, emax ; 
0642     T charge = T(1.) ;  
0643     T avgNumPhotons = GetAverageNumberOfPhotons_s2<T>(emin, emax, BetaInverse, charge ); 
0644     if(avgNumPhotons <= 0. ) return nullptr ; 
0645 
0646     LOG(debug)
0647        << " emin " << std::setw(10) << std::fixed << std::setprecision(4) << emin
0648        << " emax " << std::setw(10) << std::fixed << std::setprecision(4) << emax
0649        ;
0650 
0651     const T* ri_v = dsrc->cvalues<T>(); 
0652 
0653     unsigned ri_ni = dsrc->shape[0] ; 
0654     unsigned ri_nj = dsrc->shape[1] ; 
0655 
0656     bool dsrc_expect = ri_nj == 2 && ri_ni > 1 ; 
0657     assert(dsrc_expect );
0658     if(!dsrc_expect) std::raise(SIGINT); 
0659 
0660     unsigned ns = 1+mul ;                // sub divisions of one bin 
0661     unsigned s2_edges = getNumEdges_SplitBin<T>(mul); 
0662 
0663     NP* s2c = NP::Make<T>(s2_edges, SPLITBIN_PAYLOAD_SIZE) ; // number of values/edges is one more than bins
0664     T* s2c_v = s2c->values<T>(); 
0665 
0666     std::vector<unsigned> idxs ;  
0667     unsigned idx = 0 ; 
0668     for(unsigned p=0 ; p < SPLITBIN_PAYLOAD_SIZE ; p++) s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + p] = ( p == 0 || p == 1) ? emn : 0. ; 
0669     // artifical cumulative integral zero for 1st entry 
0670     idxs.push_back(idx); 
0671 
0672     T s2integral = 0. ; 
0673 
0674     for(unsigned i=0 ; i < ri_ni - 1 ; i++)
0675     {
0676         T en_0 = ri_v[2*(i+0)+0] ; 
0677         T ri_0 = ri_v[2*(i+0)+1] ; 
0678 
0679         T en_1 = ri_v[2*(i+1)+0] ; 
0680         T ri_1 = ri_v[2*(i+1)+1] ; 
0681 
0682         T en_01 = en_1 - en_0 ; 
0683 
0684         if(dump) std::cout 
0685              << "QCerenkovIntegral::getS2Integral_splitbin"
0686              << " i " << std::setw(3) << i 
0687              << " ns " << std::setw(3) << ns 
0688              << " en_0 " << std::setw(10) << std::fixed << std::setprecision(4) << en_0
0689              << " en_1 " << std::setw(10) << std::fixed << std::setprecision(4) << en_1
0690              << std::endl 
0691              ;
0692 
0693         // starting from s=1 skips first sub-edge as that is the same as the last edge from prior bin
0694         for(unsigned s=1 ; s < 1+mul  ; s++)   
0695         { 
0696             T en_a = en_0 + en_01*T(s-1)/T(mul) ;  // mul=1 ->s=1 only : en_a=en_0      
0697             T en_b = en_0 + en_01*T(s+0)/T(mul) ;  // mul=1 ->s=1 only : en_b=en_1     
0698 
0699             T ri_a = en_a > 0. ? dsrc->interp<T>(en_a) : -1 ; 
0700             T ri_b = en_b > 0. ? dsrc->interp<T>(en_b) : -1 ; 
0701 
0702             T s2_a = getS2<T>( BetaInverse, en_a ); 
0703             T s2_b = getS2<T>( BetaInverse, en_b ); 
0704 
0705             T sub = GetS2Integral_WithCut<T>(emin, emax, BetaInverse, en_0, en_1, ri_0, ri_1, en_a, ri_a, en_b, ri_b, dump );
0706 
0707             s2integral += sub  ; 
0708 
0709             unsigned idx = i*mul + s ; 
0710             idxs.push_back(idx); 
0711             assert( idx < s2_edges ); 
0712 
0713             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 0 ] = en_b ;        // 1st payload slot must be "domain" for NP::pdomain lookup ... 
0714             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 1 ] = en_a ;        // en_b en_a inversion needed .. as 1st entry is auto set to emn 
0715             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 2 ] = ri_a ;  
0716             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 3 ] = ri_b ;  
0717             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 4 ] = s2_a ;  
0718             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 5 ] = s2_b ;  
0719             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 6 ] = sub ;  
0720             s2c_v[SPLITBIN_PAYLOAD_SIZE*idx + 7 ] = s2integral ;   // last payload slot must be "value" for NP::pdomain lookup
0721         }
0722     }
0723 
0724     for(unsigned i=0 ; i < idxs.size() ; i++) assert( idxs[i] == i );  // check idx has it covered 
0725     return s2c ; 
0726 } 
0727 
0728 
0729 template NP* QCerenkovIntegral::getS2Integral_SplitBin( const double, unsigned, bool ) const ; 
0730 template NP* QCerenkovIntegral::getS2Integral_SplitBin( const float,  unsigned, bool ) const ; 
0731 
0732 
0733 template<typename T>
0734 unsigned QCerenkovIntegral::getNumEdges_SplitBin(unsigned mul ) const 
0735 {
0736     unsigned ri_ni = dsrc->shape[0] ; 
0737     unsigned ri_bins = ri_ni - 1 ;       // number of bins is one less than number of values 
0738     unsigned s2_bins = ri_bins*mul ; 
0739     return s2_bins + 1 ; 
0740 }
0741 
0742 
0743 
0744 
0745 
0746 template unsigned QCerenkovIntegral::getNumEdges_SplitBin<float>(unsigned ) const ; 
0747 template unsigned QCerenkovIntegral::getNumEdges_SplitBin<double>(unsigned ) const ; 
0748 
0749 /**
0750 QCerenkovIntegral::getS2Integral_SplitBin
0751 -----------------------------------
0752 
0753 Note that this is requiring the same number of s2c s2_edges for all BetaInverse.
0754 Directly using the s2c to yield the icdf would remove that requirement.
0755 
0756 **/
0757 
0758 template <typename T>
0759 NP* QCerenkovIntegral::getS2Integral_SplitBin( const NP* bis, unsigned mul, bool dump) const 
0760 {
0761     unsigned ni = bis->shape[0] ; 
0762     unsigned s2_edges = getNumEdges_SplitBin<T>(mul); 
0763     unsigned nj = s2_edges ; 
0764 
0765     NP* s2c = NP::Make<T>(ni, nj, SPLITBIN_PAYLOAD_SIZE) ; 
0766 
0767     LOG(LEVEL) << "[ creating s2c " << s2c->sstr() ; 
0768 
0769     for(unsigned i=0 ; i < ni ; i++)
0770     {
0771         const T BetaInverse = bis->get<T>(i) ; 
0772         NP* s2c_one = getS2Integral_SplitBin<T>(BetaInverse, mul, dump ); 
0773  
0774         if(s2c_one == nullptr)
0775         {
0776             LOG(info) 
0777                 << " s2c_one NULL " 
0778                 << " i " << i 
0779                 << " ni " << ni 
0780                 << " BetaInverse " << std::setw(10) << std::fixed << std::setprecision(4) << BetaInverse 
0781                 ; 
0782             continue; 
0783         }
0784         unsigned s2c_one_bytes = s2c_one->arr_bytes() ;  
0785         memcpy( s2c->bytes() + i*s2c_one_bytes, s2c_one->bytes(), s2c_one_bytes ); 
0786         s2c_one->clear(); 
0787     }
0788     LOG(LEVEL) << "] creating s2c " << s2c->sstr() ; 
0789     return s2c ; 
0790 }
0791 
0792 template NP* QCerenkovIntegral::getS2Integral_SplitBin<float>(  const NP*, unsigned, bool) const ; 
0793 template NP* QCerenkovIntegral::getS2Integral_SplitBin<double>( const NP*, unsigned, bool) const ; 
0794 
0795 
0796 
0797 
0798 template <typename T>
0799 T QCerenkovIntegral::getS2( const T BetaInverse, const T en ) const 
0800 {
0801     const T one(1.) ; 
0802     T ri = dsrc->interp<T>(en) ;  // linear interpolation of the rindex values at edges 
0803     T ct = BetaInverse/ri ; 
0804     T s2 = ( one - ct )*( one + ct ); 
0805     return s2 ; 
0806 }
0807 
0808 template double QCerenkovIntegral::getS2( const double, const double ) const ; 
0809 template float  QCerenkovIntegral::getS2( const float, const float ) const ; 
0810 
0811 
0812 
0813 /**
0814 QCerenkovIntegral::getS2Integral_UpperCut
0815 -------------------------------------
0816 
0817 Returns array of shape (nx, 2) with energies and 
0818 cululative integrals up to those energies. 
0819 Note that the energy range is adapted to correspond to the
0820 permissable range of Cerenkov for the BetaInverse to 
0821 make best use of the available bins. 
0822 
0823 The cumulative integral values are not normalized. 
0824 
0825 
0826 HMM: getting slightly non-monotonic again in some regions, 
0827 idea to avoid that is to avoid repeating the calc just save the cumulative 
0828 result in the one pass
0829 
0830 **/
0831 
0832 template <typename T>
0833 NP* QCerenkovIntegral::getS2Integral_UpperCut( const T BetaInverse, unsigned nx ) const   
0834 {
0835     T emin ; 
0836     T emax ; 
0837     T charge = T(1.) ;  
0838 
0839     T avgNumPhotons = GetAverageNumberOfPhotons_s2<T>(emin, emax, BetaInverse, charge ); 
0840     if(avgNumPhotons <= 0. ) return nullptr ; 
0841 
0842 
0843     NP* s2c = NP::Make<T>(nx, UPPERCUT_PAYLOAD_SIZE) ; 
0844     T* cc = s2c->values<T>();
0845 
0846     NP* full_s2i = GetAverageNumberOfPhotons_s2_<T>(emin, emax, BetaInverse, charge );
0847 
0848     LOG(debug) 
0849         << " BetaInverse " << std::setw(10) << std::fixed << std::setprecision(4) << BetaInverse
0850         << " emin " << std::setw(10) << std::fixed << std::setprecision(4) << emin
0851         << " emax " << std::setw(10) << std::fixed << std::setprecision(4) << emax
0852         << " avgNumPhotons " << std::setw(10) << std::fixed << std::setprecision(4) << avgNumPhotons
0853         << " avgNumPhotons*1e6 " << std::setw(10) << std::fixed << std::setprecision(4) << avgNumPhotons*1e6
0854         ;
0855 
0856     const NP* edom = NP::Linspace<T>( emin, emax, nx ); 
0857     const T* ee = edom->cvalues<T>(); 
0858 
0859     T last_ecut = 0. ; 
0860     T last_s2integral = 0. ; 
0861     NP* last_s2i = nullptr ; 
0862     bool dump = false ; 
0863 
0864     // This is doing the full integral with increasing ecut=en_b at each turn 
0865     // so it slightly grows, problem with this is that its very repetitive and hence slow
0866     // and also it risks going slightly non-monotonic from numerical imprecision.
0867     //
0868  
0869     for(unsigned i=0 ; i < nx ; i++)
0870     {
0871         T en_a = -1. ; 
0872         T en_b = ee[i] ; 
0873         
0874         NP* s2i = getS2Integral_WithCut_<T>(emin, emax, BetaInverse, en_a, en_b, dump ); 
0875 
0876         T s2integral = s2i->psum<T>(0u); 
0877         T s2_b = getS2<T>( BetaInverse, en_b ); 
0878 
0879         cc[UPPERCUT_PAYLOAD_SIZE*i+0] = en_b ;         // domain qty must be 1st payload entry for NP::pdomain
0880         cc[UPPERCUT_PAYLOAD_SIZE*i+1] = s2_b ; 
0881         cc[UPPERCUT_PAYLOAD_SIZE*i+2] = s2integral ;   // qty to be normalized needs to be in the last payload slot for NP::divide_by_last  
0882 
0883         if( i == nx - 1 )
0884         {
0885             last_ecut = en_b ; 
0886             last_s2integral = s2integral ; 
0887             last_s2i = s2i ; 
0888         }
0889         else
0890         {
0891             s2i->clear(); 
0892         }
0893     }
0894 
0895     T diff = std::abs(avgNumPhotons - last_s2integral); 
0896     bool close = diff < 1e-6 ; 
0897     if(!close)
0898     {
0899         std::cout << "NP::DumpCompare a:full_s2i b:last_s2i " << std::endl;  
0900         NP::DumpCompare<double>( full_s2i, last_s2i, 0u, 0u, 1e-6 ); 
0901 
0902         std::cout 
0903             << " QCerenkovIntegral::getS2Integral_UpperCut "
0904             << " BetaInverse " << std::setw(10) << std::fixed << std::setprecision(4) << BetaInverse
0905             << " avgNumPhotons " << std::setw(10) << std::fixed << std::setprecision(4) << avgNumPhotons   
0906             << " last_ecut " << std::setw(10) << std::fixed << std::setprecision(4) << last_ecut 
0907             << " last_s2integral " << std::setw(10) << std::fixed << std::setprecision(4) << last_s2integral
0908             << " diff " << std::setw(10) << std::fixed << std::setprecision(4) << diff
0909             << " close " << close 
0910             << std::endl 
0911             ;
0912     }
0913     assert(close);   // last cut integral should be close to the avgNumPhotons
0914 
0915     return s2c ; 
0916 }
0917 
0918 template NP* QCerenkovIntegral::getS2Integral_UpperCut( const double BetaInverse, unsigned nx ) const ; 
0919 template NP* QCerenkovIntegral::getS2Integral_UpperCut( const float  BetaInverse, unsigned nx ) const ; 
0920 
0921 
0922 /**
0923 NP* QCerenkovIntegral::getS2Integral_UpperCut
0924 -----------------------------------------
0925 
0926 BetaInverse beyond which get no photons could be used to determine 
0927 the bis range, avoiding a large chunk of the output array containing just zeros.
0928 
0929 **/
0930 
0931 template <typename T>
0932 NP* QCerenkovIntegral::getS2Integral_UpperCut( const NP* bis, unsigned nx ) const 
0933 {
0934     unsigned ni = bis->shape[0] ; 
0935     const T* bb = bis->cvalues<T>(); 
0936     unsigned nj = nx ; 
0937 
0938     NP* s2c = NP::Make<T>(ni, nj, 3) ; 
0939     LOG(info) << "[ creating s2c " << s2c->sstr() ; 
0940 
0941     for(unsigned i=0 ; i < ni ; i++)
0942     {
0943         const T BetaInverse = bb[i] ; 
0944         NP* s2c_one = getS2Integral_UpperCut<T>(BetaInverse, nx ); 
0945  
0946         if(s2c_one == nullptr)
0947         {
0948             LOG(info) 
0949                 << " s2c_one NULL " 
0950                 << " i " << i 
0951                 << " ni " << ni 
0952                 << " BetaInverse " << std::setw(10) << std::fixed << std::setprecision(4) << BetaInverse 
0953                 ; 
0954             continue; 
0955         }
0956         unsigned s2c_one_bytes = s2c_one->arr_bytes() ;  
0957         memcpy( s2c->bytes() + i*s2c_one_bytes, s2c_one->bytes(), s2c_one_bytes ); 
0958         s2c_one->clear(); 
0959     }
0960     LOG(info) << "] creating s2c " << s2c->sstr() ; 
0961     return s2c ; 
0962 }
0963 
0964 template NP* QCerenkovIntegral::getS2Integral_UpperCut<double>( const NP* , unsigned ) const ; 
0965 template NP* QCerenkovIntegral::getS2Integral_UpperCut<float>(  const NP* , unsigned ) const ; 
0966 
0967 
0968 /**
0969 QCK QCerenkovIntegral::makeICDF_UpperCut
0970 ------------------------------------
0971 
0972 ny 
0973     number BetaInverse values "height"
0974 nx
0975     number of energy domain values "width"
0976 
0977 **/
0978 
0979 template <typename T>
0980 QCK<T> QCerenkovIntegral::makeICDF_UpperCut( unsigned ny, unsigned nx, bool dump) const 
0981 {
0982     NP* bis = NP::Linspace<T>( 1. , rmx,  ny ) ;  
0983     NP* avph = getAverageNumberOfPhotons_s2<T>(bis ); 
0984     NP* s2c = getS2Integral_UpperCut<T>( bis, nx ); 
0985     NP* s2cn = s2c->copy(); 
0986     s2cn->divide_by_last<T>(); 
0987 
0988     unsigned nu = 1000 ; 
0989     unsigned hd_factor = 10 ;
0990  
0991     NP* icdf = NP::MakeICDF<T>(s2cn, nu, hd_factor, dump); 
0992     icdf->set_meta<std::string>("creator", "QCerenkovIntegral::makeICDF_UpperCut") ;  
0993     icdf->set_meta<unsigned>("hd_factor", hd_factor );
0994 
0995     NP* icdf_prop = NP::MakeProperty<T>( icdf, hd_factor ) ;
0996     icdf_prop->set_meta<std::string>("creator", "QCerenkovIntegral::makeICDF_UpperCut") ;  
0997     icdf_prop->set_meta<unsigned>("hd_factor", hd_factor );
0998 
0999 
1000     QCK<T> qck ; 
1001 
1002     qck.rindex = dsrc ; 
1003     qck.bis = bis ;  
1004     qck.avph = avph ; 
1005     qck.s2c = s2c ;
1006     qck.s2cn = s2cn ; 
1007     qck.icdf = icdf ; 
1008     qck.icdf_prop = icdf_prop ; 
1009 
1010     return qck ; 
1011 }
1012 template QCK<double> QCerenkovIntegral::makeICDF_UpperCut<double>( unsigned , unsigned, bool ) const ; 
1013 template QCK<float>  QCerenkovIntegral::makeICDF_UpperCut<float>(  unsigned , unsigned, bool ) const ; 
1014 
1015 
1016 
1017 
1018 template <typename T>
1019 QCK<T> QCerenkovIntegral::makeICDF_SplitBin( unsigned ny, unsigned mul, bool dump) const 
1020 {
1021     NP* bis = NP::Linspace<T>( 1. , rmx,  ny ) ;  
1022     std::stringstream ss ; 
1023     ss << "name:makeICDF_SplitBin,mul:" << mul ; 
1024     bis->meta = ss.str(); 
1025 
1026     NP* avph = getAverageNumberOfPhotons_s2<T>(bis ); 
1027     NP* s2c = getS2Integral_SplitBin<T>( bis, mul, dump ); 
1028     NP* s2cn = s2c->copy(); 
1029     s2cn->divide_by_last<T>(); 
1030 
1031     unsigned nu = 1000 ; 
1032     unsigned hd_factor = 10 ; 
1033     NP* icdf = NP::MakeICDF<T>(s2cn, nu, hd_factor, dump); 
1034     icdf->set_meta<std::string>("creator", "QCerenkovIntegral::makeICDF_SplitBin") ;  
1035     icdf->set_meta<unsigned>("hd_factor", hd_factor );
1036 
1037     NP* icdf_prop = NP::MakeProperty<T>( icdf, hd_factor ) ;
1038     icdf_prop->set_meta<std::string>("creator", "QCerenkovIntegral::makeICDF_UpperCut") ;  
1039     icdf_prop->set_meta<unsigned>("hd_factor", hd_factor );
1040 
1041 
1042     QCK<T> qck ; 
1043 
1044     qck.rindex = dsrc ; 
1045     qck.bis = bis ;  
1046     qck.avph = avph ; 
1047     qck.s2c = s2c ;
1048     qck.s2cn = s2cn ; 
1049     qck.icdf = icdf ; 
1050     qck.icdf_prop = icdf_prop ; 
1051 
1052     return qck ; 
1053 }
1054 template QCK<double> QCerenkovIntegral::makeICDF_SplitBin<double>( unsigned , unsigned, bool ) const ; 
1055 template QCK<float>  QCerenkovIntegral::makeICDF_SplitBin<float>(  unsigned , unsigned, bool ) const ; 
1056 
1057
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!