eic-opticks/u4/U4Scint.h

0001 #pragma once
0002 /**
0003 U4Scint
0004 ================
0005
0006 Before 1100::
0007
0008     typedef G4PhysicsOrderedFreeVector G4MaterialPropertyVector;
0009
0010 From 1100::
0011
0012     typedef G4PhysicsFreeVector G4MaterialPropertyVector;
0013
0014 And from 1100 the "Ordered" methods have been consolidated into G4PhysicsFreeVector
0015 and the class G4PhysicsOrderedFreeVector is dropped.
0016 Try to cope with this without version barnching using edit::
0017
0018    :%s/PhysicsOrderedFree/MaterialProperty/gc
0019
0020 Maybe will need to add some casts too.
0021
0022 **/
0023
0024 #include <string>
0025 #include <iomanip>
0026
0027 #include "Randomize.hh"
0028 #include "G4MaterialPropertyVector.hh"
0029 #include "G4SystemOfUnits.hh"
0030 #include "G4PhysicalConstants.hh"
0031
0032 #include "ssys.h"
0033 #include "NPFold.h"
0034 #include "U4MaterialPropertyVector.h"
0035
0036 struct U4Scint
0037 {
0038     static constexpr const bool VERBOSE = false ;
0039     static constexpr const char* PROPS = "SLOWCOMPONENT,FASTCOMPONENT,REEMISSIONPROB" ;
0040     static U4Scint* Create(const NPFold* materials );
0041
0042     const NPFold* scint ;
0043     const char* name ;
0044
0045     const NP* fast ;
0046     const NP* slow ;
0047     const NP* reem ;
0048
0049     const double epsilon ;
0050     int mismatch_0 ;
0051     int mismatch_1 ;
0052     int mismatch ;
0053
0054     const G4MaterialPropertyVector* theFastLightVector ;
0055     const G4MaterialPropertyVector* theSlowLightVector ;
0056     const G4MaterialPropertyVector* ScintillationIntegral ;
0057
0058     const NP* icdf ;
0059
0060     const int num_wlsamp ;
0061     const NP* wlsamp ;
0062
0063
0064     U4Scint(const NPFold* fold, const char* name);
0065     void init();
0066
0067     std::string desc() const ;
0068     NPFold* make_fold() const ;
0069     void save(const char* base, const char* rel=nullptr ) const ;
0070
0071     NP* createWavelengthSamples( int num_samples );
0072     NP* createGeant4InterpolatedInverseCDF(
0073         int num_bins=4096,
0074         int hd_factor=20,
0075         const char* material_name="LS",
0076         bool energy_not_wavelength=false );
0077
0078     static G4MaterialPropertyVector* Integral( const G4MaterialPropertyVector* theFastLightVector ) ;
0079     static NP* CreateWavelengthSamples(
0080         const G4MaterialPropertyVector* ScintillatorIntegral,
0081         int num_samples
0082         );
0083     static NP* CreateGeant4InterpolatedInverseCDF(
0084         const G4MaterialPropertyVector* ScintillatorIntegral,
0085         int num_bins,
0086         int hd_factor,
0087         const char* name,
0088         bool energy_not_wavelength
0089         );
0090
0091 };
0092
0093 inline U4Scint* U4Scint::Create(const NPFold* materials ) // static
0094 {
0095     std::vector<const NPFold*> subs ;
0096     std::vector<std::string> names ;
0097     materials->find_subfold_with_all_keys( subs, names, PROPS );
0098
0099     int num_subs = subs.size();
0100     int num_names = names.size() ;
0101     assert( num_subs == num_names );
0102
0103     const char* name = num_names > 0 ? names[0].c_str() : nullptr ;
0104     const NPFold* sub = num_subs > 0 ? subs[0] : nullptr ;
0105     bool with_scint = name && sub ;
0106
0107     return with_scint ? new U4Scint(sub, name) : nullptr ;
0108 }
0109
0110
0111 inline U4Scint::U4Scint(const NPFold* scint_, const char* name_)
0112     :
0113     scint(scint_),
0114     name(strdup(name_)),
0115     fast(scint->get("FASTCOMPONENT")),
0116     slow(scint->get("SLOWCOMPONENT")),
0117     reem(scint->get("REEMISSIONPROB")),
0118     epsilon(0.),
0119     mismatch_0(NP::DumpCompare<double>(fast, slow, 0, 0, epsilon)),
0120     mismatch_1(NP::DumpCompare<double>(fast, slow, 1, 1, epsilon)),
0121     mismatch(mismatch_0+mismatch_1),
0122     theFastLightVector(U4MaterialPropertyVector::FromArray(fast)),
0123     theSlowLightVector(U4MaterialPropertyVector::FromArray(slow)),
0124     ScintillationIntegral(Integral(theFastLightVector)),
0125     icdf(nullptr),
0126     num_wlsamp(ssys::getenvint("U4Scint__num_wlsamp", 0)),
0127     wlsamp(nullptr)
0128 {
0129     init();
0130 }
0131
0132 inline void U4Scint::init()
0133 {
0134     if(mismatch > 0 ) std::cerr
0135         << " mismatch_0 " << mismatch_0
0136         << " mismatch_1 " << mismatch_1
0137         << " mismatch " << mismatch
0138         << std::endl
0139         ;
0140     assert( mismatch == 0 );
0141
0142     int num_bins = 4096 ;
0143     int hd_factor = 20 ;
0144     icdf = createGeant4InterpolatedInverseCDF(num_bins, hd_factor, name) ;
0145     wlsamp = createWavelengthSamples(num_wlsamp) ;
0146 }
0147
0148
0149
0150 inline std::string U4Scint::desc() const
0151 {
0152     std::stringstream ss ;
0153     ss << "U4Scint::desc" << std::endl
0154        << " name " << name
0155        << " fast " << ( fast ? fast->sstr() : "-" )
0156        << " slow " << ( slow ? slow->sstr() : "-" )
0157        << " reem " << ( reem ? reem->sstr() : "-" )
0158        << " icdf " << ( icdf ? icdf->sstr() : "-" )
0159        << " wlsamp " << ( wlsamp ? wlsamp->sstr() : "-" )
0160        << std::endl
0161        ;
0162
0163     std::string str = ss.str();
0164     return str ;
0165 }
0166
0167 inline NPFold* U4Scint::make_fold() const
0168 {
0169     NPFold* fold = new NPFold ;
0170     fold->add("fast", fast) ;
0171     fold->add("slow", slow) ;
0172     fold->add("reem", reem) ;
0173     fold->add("icdf", icdf) ;
0174     if(wlsamp) fold->add("wlsamp", wlsamp) ;
0175     return fold ;
0176 }
0177
0178 inline void U4Scint::save(const char* base, const char* rel ) const
0179 {
0180     NPFold* fold = make_fold();
0181     fold->save(base, rel);
0182 }
0183
0184
0185 inline NP* U4Scint::createWavelengthSamples( int num_samples )
0186 {
0187     return CreateWavelengthSamples(ScintillationIntegral, num_samples );
0188 }
0189
0190 inline NP* U4Scint::createGeant4InterpolatedInverseCDF(
0191     int num_bins,
0192     int hd_factor,
0193     const char* material_name,
0194     bool energy_not_wavelength
0195     )
0196 {
0197     return CreateGeant4InterpolatedInverseCDF(
0198                ScintillationIntegral,
0199                num_bins,
0200                hd_factor,
0201                material_name,
0202                energy_not_wavelength  );
0203 }
0204
0205
0206 /**
0207 U4Scint::Integral
0208 ---------------------------
0209
0210 Returns cumulative sum of the input property on the same energy domain,
0211 with values starting at 0. and increasing monotonically.
0212
0213 The is using trapezoidal numerical integration.
0214
0215
0216 **/
0217
0218 inline G4MaterialPropertyVector* U4Scint::Integral( const G4MaterialPropertyVector* theFastLightVector )
0219 {
0220      G4MaterialPropertyVector* aMaterialPropertyVector = new G4MaterialPropertyVector();
0221
0222           if (theFastLightVector) {
0223
0224                G4double currentIN = (*theFastLightVector)[0];
0225
0226                 if (currentIN >= 0.0) {
0227
0228                     // Create first (photon energy, Scintillation
0229                     // Integral pair
0230
0231                     G4double currentPM = theFastLightVector->
0232                         Energy(0);
0233
0234                     G4double currentCII = 0.0;
0235
0236                     aMaterialPropertyVector->
0237                         InsertValues(currentPM , currentCII);
0238
0239                     // Set previous values to current ones prior to loop
0240
0241                     G4double prevPM  = currentPM;
0242                     G4double prevCII = currentCII;
0243                     G4double prevIN  = currentIN;
0244
0245                     // loop over all (photon energy, intensity)
0246                     // pairs stored for this material
0247
0248                     for(size_t ii = 1;
0249                               ii < theFastLightVector->GetVectorLength();
0250                               ++ii)
0251                     {
0252                         currentPM = theFastLightVector->Energy(ii);
0253
0254                         currentIN= (*theFastLightVector)[ii];
0255
0256                         currentCII = 0.5 * (prevIN + currentIN);
0257
0258                         currentCII = prevCII +
0259                             (currentPM - prevPM) * currentCII;
0260
0261                         aMaterialPropertyVector->
0262                             InsertValues(currentPM, currentCII);
0263
0264                         prevPM  = currentPM;
0265                         prevCII = currentCII;
0266                         prevIN  = currentIN;
0267                     }
0268                }
0269             }
0270
0271     return aMaterialPropertyVector ;
0272 }
0273
0274
0275 inline NP* U4Scint::CreateWavelengthSamples(
0276     const G4MaterialPropertyVector* ScintillatorIntegral_,
0277     int num_samples )
0278 {
0279     if(num_samples == 0) return nullptr ;
0280
0281     G4MaterialPropertyVector* ScintillatorIntegral = const_cast<G4MaterialPropertyVector*>(ScintillatorIntegral_) ;
0282
0283     double mx = ScintillatorIntegral->GetMaxValue() ;
0284     std::cerr
0285         << "U4Scint::CreateWavelengthSamples"
0286         << " ScintillatorIntegral.max*1e9 "
0287         << std::fixed << std::setw(10) << std::setprecision(4) << mx*1e9
0288         ;
0289
0290     NP* wl = NP::Make<double>(num_samples);
0291     wl->fill<double>(0.);
0292     double* wl_v = wl->values<double>() ;
0293
0294     for(int i=0 ; i < num_samples ; i++)
0295     {
0296         G4double u = G4UniformRand() ;
0297         G4double CIIvalue = u*mx;
0298         G4double sampledEnergy = ScintillatorIntegral->GetEnergy(CIIvalue);  // from value to domain
0299
0300         G4double sampledWavelength_nm = h_Planck*c_light/sampledEnergy/nm ;
0301
0302         wl_v[i] = sampledWavelength_nm ;
0303
0304         if( i < 10 ) std::cout
0305             << " sampledEnergy/eV "
0306             << std::fixed << std::setw(10) << std::setprecision(4) << sampledEnergy/eV
0307             << " sampledWavelength_nm "
0308             <<  std::fixed << std::setw(10) << std::setprecision(4) << sampledWavelength_nm
0309             << std::endl
0310             ;
0311     }
0312     return wl ;
0313 }
0314
0315
0316 /**
0317 U4Scint::CreateGeant4InterpolatedInverseCDF
0318 -----------------------------------------------------
0319
0320 Reproducing the results of Geant4 dynamic bin finding interpolation
0321 using GPU texture lookups demands very high resolution textures for some
0322 ICDF shapes. This function prepares a three item buffer that can be used
0323 to create a 2D texture that effectively mimmicks variable bin sizing even
0324 though GPU hardware does not support that without paying the cost of
0325 high resolution across the entire range.
0326
0327 * item 0 : full range "standard" resolution
0328 * item 1: left hand side high resolution
0329 * item 2: right hand side high resolution
0330
0331 ::
0332
0333     hd_factor                LHS            RHS
0334     10          10x bins:    0.00->0.10     0.90->1.00
0335     20          20x bins:    0.00->0.05     0.95->1.00
0336
0337
0338 The ICDF is formed using Geant4s "domain lookup from value" functionality
0339 in the form of G4MaterialPropertyVector::GetEnergy
0340
0341 ::
0342
0343     g4-cls G4MaterialPropertyVector
0344
0345     096 G4double G4PhysicsOrderedFreeVector::GetEnergy(G4double aValue)
0346      97 {
0347      98         G4double e;
0348      99         if (aValue <= GetMinValue()) {
0349     100           e = edgeMin;
0350     101         } else if (aValue >= GetMaxValue()) {
0351     102           e = edgeMax;
0352     103         } else {
0353     104           size_t closestBin = FindValueBinLocation(aValue);
0354     105           e = LinearInterpolationOfEnergy(aValue, closestBin);
0355     106     }
0356     107         return e;
0357     108 }
0358
0359     118 G4double G4PhysicsOrderedFreeVector::LinearInterpolationOfEnergy(G4double aValue,
0360     119                                  size_t bin)
0361     120 {
0362     121         G4double res = binVector[bin];
0363     122         G4double del = dataVector[bin+1] - dataVector[bin];
0364     123         if(del > 0.0) {
0365     124           res += (aValue - dataVector[bin])*(binVector[bin+1] - res)/del;
0366     125         }
0367     126         return res;
0368     127 }
0369
0370
0371
0372
0373                                                         1  (x1,y1)     (  binVector[bin+1], dataVector[bin+1] )
0374                                                        /
0375                                                       /
0376                                                      *  ( xv,yv )       ( res, aValue )
0377                                                     /
0378                                                    /
0379                                                   0  (x0,y0)          (  binVector[bin], dataVector[bin] )
0380
0381
0382               Similar triangles::
0383
0384                  xv - x0       x1 - x0
0385                ---------- =   -----------
0386                  yv - y0       y1 - y0
0387
0388
0389
0390
0391                   res - binVector[bin]             binVector[bin+1] - binVector[bin]
0392                ----------------------------  =     -----------------------------------
0393                  aValue - dataVector[bin]          dataVector[bin+1] - dataVector[bin]
0394
0395
0396                                                                               binVector[bin+1] - binVector[bin]
0397                    res  = binVector[bin] +  ( aValue - dataVector[bin] ) *  -------------------------------------
0398                                                                               dataVector[bin+1] - dataVector[bin]
0399
0400                                                    x1 - x0
0401                    xv  =    x0  +   (yv - y0) *  --------------
0402                                                    y1 - y0
0403
0404 **/
0405
0406 inline NP* U4Scint::CreateGeant4InterpolatedInverseCDF(
0407        const G4MaterialPropertyVector* ScintillatorIntegral_,
0408        int num_bins,
0409        int hd_factor,
0410        const char* material_name,
0411        bool energy_not_wavelength
0412 )   // static
0413 {
0414
0415     assert( material_name );
0416
0417     G4MaterialPropertyVector* ScintillatorIntegral = const_cast<G4MaterialPropertyVector*>(ScintillatorIntegral_) ;  // tut tut : G4 GetMaxValue() GetEnergy() non-const
0418     double mx = ScintillatorIntegral->GetMaxValue() ;   // dataVector.back(); because its **ORDERED** to be increasing on Insert
0419
0420
0421     // hmm more extensible (eg for Cerenkov [BetaInverse,u,payload] icdf)
0422     // with the 3 for the different resolutions to be in the payload rather than as separate items ?
0423     // would of course use 4 to map to float4 after narrowing
0424
0425
0426     NP* icdf = NP::Make<double>(3, num_bins, 1);
0427     icdf->fill<double>(0.);
0428
0429     int ni = icdf->shape[0];
0430     int nj = icdf->shape[1];
0431     int nk = icdf->shape[2];
0432
0433     assert( ni == 3 );
0434     assert( nk == 1 );
0435     int k = 0 ;
0436
0437     assert( hd_factor == 10 || hd_factor == 20 );
0438     double edge = 1./double(hd_factor) ;
0439
0440     icdf->names.push_back(material_name) ;  // match X4/GGeo
0441     icdf->set_meta<std::string>("name", material_name );
0442
0443     icdf->set_meta<std::string>("creator", "U4Scint::CreateGeant4InterpolatedInverseCDF" );
0444     icdf->set_meta<int>("hd_factor", hd_factor );
0445     icdf->set_meta<int>("num_bins", num_bins );
0446     icdf->set_meta<double>("edge", edge );
0447
0448
0449     if(VERBOSE) std::cerr
0450         << "U4Scint::CreateGeant4InterpolatedInverseCDF"
0451         << " num_bins " << num_bins
0452         << " hd_factor " << hd_factor
0453         << " mx " << std::fixed << std::setw(10) << std::setprecision(4) << mx
0454         << " mx*1e9 " << std::fixed << std::setw(10) << std::setprecision(4) << mx*1e9
0455         << " edge " << std::fixed << std::setw(10) << std::setprecision(4) << edge
0456         << " icdf " << icdf->sstr()
0457         << std::endl
0458         ;
0459
0460     for(int j=0 ; j < nj ; j++)
0461     {
0462         double u_all = double(j)/double(nj) ;                         // 0 -> (nj-1)/nj = 1-1/nj
0463         double u_lhs = double(j)/double(hd_factor*nj) ;               // hd_factor=10(20) u_lhs=0.0->0.1  (0.0  ->~0.05)
0464         double u_rhs = 1. - edge + double(j)/double(hd_factor*nj) ;   // hd_factor=10(20) u_rhs=0.9->~1.0 (0.95 ->~1.0)
0465         // u_lhs and u_rhs mean there are hd_factor more lookups at the extremes
0466
0467         double energy_all = ScintillatorIntegral->GetEnergy( u_all*mx );
0468         double energy_lhs = ScintillatorIntegral->GetEnergy( u_lhs*mx );
0469         double energy_rhs = ScintillatorIntegral->GetEnergy( u_rhs*mx );
0470
0471         double wavelength_all = h_Planck*c_light/energy_all/nm ;
0472         double wavelength_lhs = h_Planck*c_light/energy_lhs/nm ;
0473         double wavelength_rhs = h_Planck*c_light/energy_rhs/nm ;
0474
0475         double v_all = energy_not_wavelength ? energy_all :  wavelength_all ;
0476         double v_lhs = energy_not_wavelength ? energy_lhs :  wavelength_lhs ;
0477         double v_rhs = energy_not_wavelength ? energy_rhs :  wavelength_rhs ;
0478
0479         icdf->set<double>(v_all, 0, j, k );
0480         icdf->set<double>(v_lhs, 1, j, k );
0481         icdf->set<double>(v_rhs, 2, j, k );
0482     }
0483     return icdf ;
0484 }
0485