eic-opticks/qudarap/QCK.hh

0001 #pragma once
0002
0003 #include <functional>
0004 #include "QUDARAP_API_EXPORT.hh"
0005
0006 /**
0007 QCK : carrier struct holding Cerenkov lookup arrays, created by QCerenkov
0008 ===========================================================================
0009
0010 Instances of QCK are created by *QCerenkov::makeICDF*
0011
0012 Contains NP arrays:
0013
0014 rindex
0015     (ni, 2)  refractive index values together with energy domain values (eV)
0016     for a single material
0017
0018 bis
0019     (ny,) : BetaInverse domain values in range from from 1. to RINDEX_max
0020
0021 s2c
0022     (ny,nx,2) : *ny* cumulative Cerenkov sin^2 theta "s2" energy integrals up to *nx* energy cuts.
0023     Last "payload" dimension is 2 as the energy domain values are kept together with the
0024     cumulative integral values. This is necessary as the energies are not common,
0025     with each BetaInverse having different energy ranges over which Cerenkov photons
0026     can be produced.
0027
0028 s2cn
0029     (ny,nx,2) : normalized version of *s2c* with all values divided by the maximum energy
0030     integral values giving the inverse-CDF ICDF. This is usable via NP::pdomain
0031     to perform domain lookups to generate Cerenkov wavelengths given input BetaInverse
0032     and a random number.
0033
0034 **/
0035
0036 struct NP ;
0037
0038 template <typename T>
0039 struct QUDARAP_API QCK
0040 {
0041     static const double DT_SCALE ;
0042
0043     std::string lfold ;
0044     NP* rindex ;
0045     NP* bis ;  // BetaInverse
0046     NP* avph ; // for each BetaInverse payload : (BetaInverse, emin, emax, averageNumberOfPhotons)
0047     NP* s2c ;  // cumulative s2 integral
0048     NP* s2cn ; // normalized *s2c*, aka *icdf*
0049     NP* icdf ;
0050     NP* icdf_prop ;
0051
0052     T emn ;
0053     T emx ;
0054     T rmn ;
0055     T rmx ;
0056
0057     void init() ;
0058
0059     void save(const char* base, const char* reldir=nullptr) const ;
0060     static QCK* Load(const char* base, const char* reldir=nullptr);
0061
0062
0063     int find_bisIdx( const T BetaInverse ) const ;
0064     void energy_range_s2cn( T& emin, T& emax, const T BetaInverse, bool dump ) const ;
0065     void energy_range_avph( T& emin, T& emax, const T BetaInverse, bool dump ) const ;
0066
0067
0068     // lookup from sets of ICDF, normalized s2 energy integrals
0069     bool is_permissable( const T BetaInverse) const ;
0070
0071     T   energy_lookup_( const T BetaInverse, const T u, double& dt, bool icdf_ ) const ;
0072     NP* energy_lookup(  const T BetaInverse, const NP* uu,  NP* tt, bool icdf_ ) const ;
0073
0074     // traditional s2 rejection sampling using rindex as function of energy
0075     T   energy_sample_( const T emin, const T emax, const T BetaInverse,  const std::function<T()>& rng, double& dt, unsigned& count ) const ;
0076     NP* energy_sample(  const T BetaInverse,  const std::function<T()>& rng, unsigned ni, NP* tt ) const ;
0077
0078 };
0079