qudarap/tests/QCerenkovIntegralTest.cc

0001 /**
0002 QCerenkovIntegralTest
0003 ======================
0004
0005 Performs many RINDEX s2 integrals for many BetaInverse values
0006 in order to construct an ICDF that is saved into $TMP/QCerenkovIntegralTest
0007
0008 **/
0009
0010 #include "spath.h"
0011 #include "NP.hh"
0012 #include "NPX.h"
0013 #include "QCerenkovIntegral.hh"
0014 #include "QCK.hh"
0015
0016 #include "scuda.h"
0017 #include "OPTICKS_LOG.hh"
0018
0019
0020 /**
0021 test_GetAverageNumberOfPhotons_s2
0022 -----------------------------------
0023
0024 See ana/ckn.py:compareOtherScans for check that these
0025 results match those from the python prototype and cks
0026
0027 **/
0028
0029 const char* BASE = "$TMP/QCerenkovIntegralTest" ;
0030
0031
0032 void test_GetAverageNumberOfPhotons_s2(const QCerenkovIntegral& ck)
0033 {
0034     LOG(info) ;
0035     const double charge = 1. ;
0036     double numPhotons, emin, emax ;
0037
0038     unsigned ni = 1001u ;
0039     NP* bis = NP::Linspace<double>(1., 2., ni );
0040     const double* bb = bis->cvalues<double>();
0041
0042     NP* scan = NP::Make<double>(ni, 4);
0043     double* ss = scan->values<double>();
0044
0045     for(unsigned i=0 ; i < unsigned(bis->shape[0]) ; i++ )
0046     {
0047         const double BetaInverse = bb[i] ;
0048         numPhotons = ck.GetAverageNumberOfPhotons_s2<double>(emin, emax, BetaInverse, charge );
0049
0050         ss[4*i+0] = BetaInverse ;
0051         ss[4*i+1] = numPhotons ;
0052         ss[4*i+2] = emin ;
0053         ss[4*i+3] = emax ;
0054
0055         if(numPhotons != 0.)
0056             std::cout
0057                 << " i " << std::setw(5) << i
0058                 << " BetaInverse " << std::fixed << std::setw(10) << std::setprecision(4) << BetaInverse
0059                 << " numPhotons " << std::fixed << std::setw(10) << std::setprecision(4) << numPhotons
0060                 << " emin " << std::fixed << std::setw(10) << std::setprecision(4) << emin
0061                 << " emax " << std::fixed << std::setw(10) << std::setprecision(4) << emax
0062                 << std::endl
0063                 ;
0064     }
0065
0066     const char* path = spath::Resolve(BASE, "test_GetAverageNumberOfPhotons_s2.npy" );
0067     LOG(info) << " save to " << path ;
0068     scan->save(path);
0069 }
0070
0071
0072 void test_getS2Integral_UpperCut_one( const QCerenkovIntegral& ck, double BetaInverse )
0073 {
0074     unsigned nx = 100 ;
0075     NP* s2c = ck.getS2Integral_UpperCut(BetaInverse, nx);
0076
0077     LOG(info)
0078         << " BetaInverse " <<  std::fixed << std::setw(10) << std::setprecision(4) << BetaInverse
0079         << " s2c " << s2c->desc()
0080         ;
0081
0082     const char* path = spath::Resolve(BASE, "test_getS2Integral_UpperCut_one.npy" );
0083     LOG(info) << " save to " << path ;
0084     s2c->save(path);
0085 }
0086
0087
0088
0089
0090 void test_getS2Integral_UpperCut(const QCerenkovIntegral& ck )
0091 {
0092     const NP* bis  = NP::Linspace<double>( 1., 2. , 1001 );     // BetaInverse
0093     unsigned nx = 100u ;
0094     NP* s2c = ck.getS2Integral_UpperCut<double>(bis, nx );
0095
0096
0097     LOG(info)
0098         << std::endl
0099         << " bis    " << bis->desc()
0100         << std::endl
0101         << " s2c " << s2c->desc()
0102         ;
0103
0104     const char* fold = spath::Resolve(BASE, "test_getS2Integral_UpperCut" );
0105     LOG(info) << " save to " << fold ;
0106     s2c->save(fold,"s2c.npy");
0107
0108     s2c->divide_by_last<double>();
0109     s2c->save(fold,"s2cn.npy");
0110 }
0111
0112
0113 void test_getS2Integral_SplitBin( const QCerenkovIntegral& ck, const char* bis_, unsigned mul, bool dump )
0114 {
0115     LOG(info) << "[" ;
0116     const NP* bis  = bis_ == nullptr ? NP::Linspace<double>( 1., 2. , 1001 ) : NPX::FromString<double>(bis_);
0117     NP* s2c = ck.getS2Integral_SplitBin<double>(bis, mul, dump);
0118
0119     const char* fold = spath::Resolve(BASE, "test_getS2Integral_SplitBin");
0120
0121     LOG(info)
0122         << " bis_ " << bis_
0123         << " mul " << mul
0124         << " dump " << dump
0125         << " s2c " << s2c->sstr()
0126         << std::endl
0127         << " fold " << fold
0128         ;
0129
0130     bis->save(fold, "bis.npy");
0131     s2c->save(fold, "s2c.npy");
0132     s2c->divide_by_last<double>();
0133     s2c->save(fold, "s2cn.npy");
0134     LOG(info) << "]" ;
0135 }
0136
0137
0138 void test_makeICDF_UpperCut(const QCerenkovIntegral& ck, unsigned ny, unsigned nx, bool dump )
0139 {
0140     QCK<double> qck = ck.makeICDF_UpperCut<double>( ny, nx, dump );
0141
0142     LOG(info)
0143         << std::endl
0144         << " qck.bis  " << qck.bis->desc()
0145         << std::endl
0146         << " qck.s2c  " << qck.s2c->desc()
0147         << std::endl
0148         << " qck.s2cn " << qck.s2cn->desc()
0149         << std::endl
0150         ;
0151
0152     const char* qck_path = spath::Resolve(BASE, "test_makeICDF_UpperCut");
0153     qck.save(qck_path);
0154 }
0155
0156
0157
0158 void test_makeICDF_SplitBin(const QCerenkovIntegral& ck, unsigned ny, unsigned mul, bool dump )
0159 {
0160     QCK<double> qck = ck.makeICDF_SplitBin<double>( ny, mul, dump );
0161
0162     LOG(info)
0163         << std::endl
0164         << " qck.bis  " << qck.bis->desc()
0165         << std::endl
0166         << " qck.s2c  " << qck.s2c->desc()
0167         << std::endl
0168         << " qck.s2cn " << qck.s2cn->desc()
0169         << std::endl
0170         ;
0171
0172     const char* qck_path = spath::Resolve(BASE, "test_makeICDF_SplitBin");
0173     qck.save(qck_path);
0174 }
0175
0176
0177
0178
0179 int main(int argc, char** argv)
0180 {
0181     OPTICKS_LOG(argc, argv);
0182
0183     const char* _path = "$HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/material/LS/RINDEX.npy" ;
0184     const char* path = spath::Resolve(_path) ;
0185     LOG(info)
0186         << " _path " << path
0187         << " path " << path
0188         ;
0189
0190     char d = 'N' ;
0191     char t = argc > 1 ? argv[1][0] : d  ;
0192     LOG(info) << " t " << t ;
0193
0194     QCerenkovIntegral ck(path) ;
0195     if(ck.dsrc == nullptr) return 0 ;
0196
0197
0198     if( t == 'N' )
0199     {
0200         test_GetAverageNumberOfPhotons_s2(ck);
0201     }
0202     else if( t == 'A' )
0203     {
0204         const char* bis = "1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.792" ;
0205         //const char* bis = "1.0" ;
0206         unsigned mul = 10 ;
0207         bool dump = true ;
0208         test_getS2Integral_SplitBin(ck, bis, mul, dump) ;
0209     }
0210     else if ( t == 'U' )
0211     {
0212         unsigned ny = 2000u ;
0213         unsigned nx = 2000u ;
0214         bool dump = true ;
0215         test_makeICDF_UpperCut(ck, ny, nx, dump );
0216     }
0217     else if ( t == 'S' )
0218     {
0219         unsigned ny = 1000u ;
0220         unsigned mul = 100u ;
0221         bool dump = true ;
0222         test_makeICDF_SplitBin(ck, ny, mul, dump );
0223     }
0224     return 0 ;
0225 }
0226