eic-opticks/u4/U4SurfaceArray.h

0001 #pragma once
0002 /**
0003 U4SurfaceArray.h
0004 ==================
0005
0006 Formerly did this in U4Surface::MakeStandardArray
0007
0008 * this is a reimplementation that follows the intent of GSurfaceLib::createStandardSurface
0009 * for surfaces only payload group zero is filled with four payload
0010   probabilities that sum to one::
0011
0012     (detect, absorb, reflect_specular, reflect_diffuse)
0013
0014 * EFFICIENCY, REFLECTIVITY and specular/diffuse optical surface nature
0015   are inputs to setting the four probabilities that must sum to 1.
0016
0017 * notice that implicit surfaces are like perfects in that there
0018   is no wavelength variation of properties : the implicits
0019   always just perfectly absorb corresponding to Geant4 fStopAndKill
0020   behaviour when going from RINDEX to non-RINDEX material.
0021
0022
0023 Only simple surfaces are covered by this
0024 ------------------------------------------
0025
0026 U4SurfaceArray currently only handles simple surfaces where each surface
0027 is either:
0028
0029 * sensor
0030 * absorber
0031 * specular refector
0032 * diffuse reflector
0033
0034 As implemented in qsim::propagate_at_surface
0035
0036 Other types of surfaces need to use a different approach and
0037 steered to the appropriate qsim method by qsim::propagate
0038 based on the "ems" or smatsur.h enum which is planted
0039 into the optical buffer by sstandard::make_optical
0040
0041
0042 **/
0043
0044 #include <vector>
0045 #include <string>
0046
0047 #include "G4SystemOfUnits.hh"
0048 #include "G4LogicalSurface.hh"
0049 #include "G4OpticalSurface.hh"
0050 #include "G4MaterialPropertiesTable.hh"
0051
0052 #include "NP.hh"
0053 #include "sdomain.h"
0054 #include "sprop.h"
0055
0056 #include "U4SurfacePerfect.h"
0057 #include "U4OpticalSurfaceFinish.h"
0058
0059 struct U4SurfaceArray
0060 {
0061     static constexpr const bool VERBOSE = false ;
0062     static constexpr const double UNSET = -1. ;
0063     struct D4 { double x,y,z,w ; } ;
0064
0065     sdomain dom ;
0066     int num_surface ;
0067     int num_implicit ;
0068     int num_perfect ;
0069     int ni ;
0070     int nj ;
0071     int nk ;
0072     int nl ;
0073     int j  ;  // 0:payload group
0074
0075     NP* sur ;
0076     double* sur_v ;
0077
0078     U4SurfaceArray(
0079         const std::vector<const G4LogicalSurface*>& surface,
0080         const std::vector<std::string>& implicit,
0081         const std::vector<U4SurfacePerfect>& perfect
0082         );
0083
0084     void addSurface( int i, const G4LogicalSurface* ls);
0085     void addImplicit(int i, const char* name);
0086     void addPerfect( int i, const U4SurfacePerfect& perfect );
0087 };
0088
0089
0090 inline U4SurfaceArray::U4SurfaceArray(
0091         const std::vector<const G4LogicalSurface*>& surface,
0092         const std::vector<std::string>& implicit,
0093         const std::vector<U4SurfacePerfect>& perfect  )
0094     :
0095     num_surface(surface.size()),
0096     num_implicit(implicit.size()),
0097     num_perfect(perfect.size()),
0098     ni(num_surface + num_implicit + num_perfect),
0099     nj(sprop::NUM_PAYLOAD_GRP),
0100     nk(dom.length),
0101     nl(sprop::NUM_PAYLOAD_VAL),
0102     j(0),   // payload group
0103     sur(NP::Make<double>(ni, nj, nk, nl )),
0104     sur_v(sur->values<double>())
0105 {
0106     sur->fill<double>(UNSET); // matching X4/GGeo
0107
0108     for(int i=0 ; i < ni ; i++)
0109     {
0110         if( i < num_surface )
0111         {
0112             const G4LogicalSurface* ls = surface[i] ;
0113             addSurface(i, ls);
0114         }
0115         else if( i < num_surface + num_implicit )
0116         {
0117             const std::string& impl = implicit[i-num_surface] ;
0118             addImplicit(i, impl.c_str());
0119         }
0120         else if( i < num_surface + num_implicit + num_perfect )
0121         {
0122             const U4SurfacePerfect& perf = perfect[i-num_surface-num_implicit] ;
0123             addPerfect(i, perf);
0124         }
0125     }
0126     assert( int(sur->names.size()) == ni );
0127 }
0128
0129
0130
0131
0132 /**
0133 U4SurfaceArray::addSurface
0134 ----------------------------
0135
0136 Payload values::
0137
0138    detect
0139    absorb
0140    reflect_specular
0141    reflect_diffuse
0142
0143 depend on REFLECTIVITY and EFFICIENCY properties.
0144 BUT, there is not a one-to-one relationship between
0145 the properties and the content of the sur array.
0146
0147 Rather there is dependency on the surface being a sensor
0148 (having any efficiency greater than zero) and also
0149 the U4OpticalSurfaceFinish::IsPolished specular nature.
0150 Those inputs determine what is filled into the payload
0151 slots that all add to one and are treated as probabilities.
0152
0153 Note assumption that sensors do not reflect, this is old traditional POM
0154 for more involved modelling need to use special surfaces, see smatsur.h.
0155
0156
0157 **/
0158
0159 inline void U4SurfaceArray::addSurface(int i, const G4LogicalSurface* ls)
0160 {
0161     const G4String& name = ls->GetName() ;
0162     sur->names.push_back(name.c_str()) ;
0163
0164     G4OpticalSurface* os = dynamic_cast<G4OpticalSurface*>(ls->GetSurfaceProperty());
0165     unsigned finish = os->GetFinish() ;
0166     bool is_specular = U4OpticalSurfaceFinish::IsPolished( finish );
0167
0168     G4MaterialPropertiesTable* mpt = os->GetMaterialPropertiesTable() ;
0169     assert( mpt );
0170     //if( mpt == nullptr ) std::cerr << "U4Surface::MakeStandardArray NO MPT " << name << std::endl ;
0171
0172     G4MaterialPropertyVector* EFFICIENCY = mpt ? mpt->GetProperty("EFFICIENCY") : nullptr ;
0173     G4MaterialPropertyVector* REFLECTIVITY = mpt ? mpt->GetProperty("REFLECTIVITY") : nullptr ;
0174
0175     double max_effi = 0. ;
0176     double max_refl = 0. ;
0177
0178     for(int k=0 ; k < nk ; k++)  // energy/wavelength domain
0179     {
0180         double energy = dom.energy_eV[k] * eV ;
0181         double effi = EFFICIENCY ? EFFICIENCY->Value(energy) : 0. ;
0182         double refl = REFLECTIVITY ? REFLECTIVITY->Value(energy) : 0. ;
0183         if( effi > max_effi ) max_effi = effi ;
0184         if( refl > max_refl ) max_refl = refl ;
0185     }
0186
0187     bool is_sensor = max_effi > 0. ;
0188
0189     if(VERBOSE) std::cout
0190         << "U4SurfaceArray::addSurface"
0191         << " name " << std::setw(30) << name
0192         << " max_effi " << std::setw(10) << std::fixed << std::setprecision(4) << max_effi
0193         << " max_refl " << std::setw(10) << std::fixed << std::setprecision(4) << max_refl
0194         << " is_sensor " << ( is_sensor ? "YES" : "NO " )
0195         << " is_specular " << ( is_specular ? "YES" : "NO " )
0196         << " finish " << U4OpticalSurfaceFinish::Name(finish)
0197         << std::endl
0198         ;
0199
0200     for(int k=0 ; k < nk ; k++)  // energy/wavelength domain
0201     {
0202         double energy = dom.energy_eV[k] * eV ;
0203         double effi = EFFICIENCY ? EFFICIENCY->Value(energy) : 0. ;
0204         double refl = REFLECTIVITY ? REFLECTIVITY->Value(energy) : 0. ;
0205
0206        D4 d4 ;
0207
0208         if( is_sensor )
0209         {
0210             d4.x = effi ;       // detect
0211             d4.y = 1. - effi ;  // absorb
0212             d4.z = 0. ;         // reflect_specular
0213             d4.w = 0. ;         // reflect_diffuse
0214         }
0215         else
0216         {
0217             if( is_specular )
0218             {
0219                 d4.x = 0. ;          // detect
0220                 d4.y = 1. - refl  ;  // absorb
0221                 d4.z = refl ;        // reflect_specular
0222                 d4.w = 0. ;          // reflect_diffuse
0223             }
0224             else
0225             {
0226                 d4.x = 0. ;          // detect
0227                 d4.y = 1. - refl  ;  // absorb
0228                 d4.z = 0. ;          // reflect_specular
0229                 d4.w = refl ;        // reflect_diffuse
0230             }
0231         }
0232         double d4_sum = d4.x + d4.y + d4.z + d4.w ;
0233         assert( std::abs( d4_sum - 1. ) < 1e-9 );
0234         int index = i*nj*nk*nl + j*nk*nl + k*nl ;
0235         sur_v[index+0] = d4.x ;
0236         sur_v[index+1] = d4.y ;
0237         sur_v[index+2] = d4.z ;
0238         sur_v[index+3] = d4.w ;
0239     }
0240 }
0241
0242 /**
0243 U4SurfaceArray::addImplicit
0244 ----------------------------
0245
0246 Implicits are perfect absorbers that mimic within Opticks on GPU
0247 the implicit fStopAndKill absorption behaviour of Geant4 when
0248 photons reach a boundary from a material with RINDEX property
0249 to a material without RINDEX property.
0250 Typically this is photons hitting non-transparent materials.
0251
0252 **/
0253
0254 inline void U4SurfaceArray::addImplicit(int i, const char* name)
0255 {
0256     U4SurfacePerfect impl = { name, 0., 1., 0., 0. } ;
0257     addPerfect(i, impl );
0258 }
0259
0260 /**
0261 U4SurfaceArray::addPerfect
0262 ----------------------------
0263
0264 Perfect surfaces are used for debugging and unrealistic tests.
0265 They have constant properties across the energy/wavelength domain,
0266 implemented by simply copying the values.
0267
0268 **/
0269
0270 inline void U4SurfaceArray::addPerfect(int i, const U4SurfacePerfect& perf )
0271 {
0272     sur->names.push_back(perf.name.c_str()) ;
0273     for(int k=0 ; k < nk ; k++)  // energy/wavelength domain
0274     {
0275         assert( std::abs( perf.sum() - 1. ) < 1e-9 );
0276         int index = i*nj*nk*nl + j*nk*nl + k*nl ;
0277         sur_v[index+0] = perf.detect ;
0278         sur_v[index+1] = perf.absorb ;
0279         sur_v[index+2] = perf.reflect_specular ;
0280         sur_v[index+3] = perf.reflect_diffuse ;
0281     }
0282 }
0283