Material and Surface Properties =================================== Opticks handling of material and surface properties is all about preparing GPU textures : which means they have to be standardized to fit into arrays with specific properties and common wavelength domain. Materials must have a non-zero REEMISSIONPROB for it to be regarded as a scintillator and collected by GScintillatorLib. The below guide to the relevant classes is intended to show how material/surface properties get translated into GPU textures and how those are accessed GPU side. GPU texture preparation -------------------------- ggeo/GScintillatorLib scintillator lib is populated by GGeo::prepareScintillatorLib (invoked by GGeo::prepare) from raw materials with the 3 properties : SLOWCOMPONENT,FASTCOMPONENT,REEMISSIONPROB For a material to be regarded as a scintillator it must have all these three properties. Also GScintillatorLib prepares the reemission inverse CDF array. optixrap/OScintillatorLib Converts the GScintillatorLib array into the reemission GPU texture. Lookups on the reemission texture with uniform randoms generates wavelengths according to the desired distribution. ggeo/GMaterialLib,GMaterial materials get standardized into 8 properties and a common wavelength domain, as they are destined to get interleaved into the boundary GPU texture. For a material to be a scintillator must have non-zero "reemission_prob:REEMISSIONPROB,” The 8 material properties (3 are spare). :: 63 const char* GMaterialLib::keyspec = 64 "refractive_index:RINDEX," 65 "absorption_length:ABSLENGTH," 66 "scattering_length:RAYLEIGH," 67 "reemission_prob:REEMISSIONPROB," 68 "group_velocity:GROUPVEL," 69 "extra_y:EXTRA_Y," 70 "extra_z:EXTRA_Z," 71 "extra_w:EXTRA_W," ggeo/GSurfaceLib,GSurface surface properties are standardized into 8 properties (4 are spare) :: 51 const char* GSurfaceLib::detect = "detect" ; 52 const char* GSurfaceLib::absorb = "absorb" ; 53 const char* GSurfaceLib::reflect_specular = "reflect_specular" ; 54 const char* GSurfaceLib::reflect_diffuse = "reflect_diffuse" ; 55 56 const char* GSurfaceLib::extra_x = "extra_x" ; 57 const char* GSurfaceLib::extra_y = "extra_y" ; 58 const char* GSurfaceLib::extra_z = "extra_z" ; 59 const char* GSurfaceLib::extra_w = "extra_w" ; 60 87 const char* GSurfaceLib::keyspec = 88 "detect:EFFICIENCY," 89 "absorb:DUMMY," 90 "reflect_specular:REFLECTIVITY," 91 "reflect_diffuse:REFLECTIVITY," 92 ; ggeo/GBnd A boundary is the 4 integers (omat,osur,isur,imat) representing materials and surfaces on either side of a transition been materials. :: omat:outer-material osur:outer-surface isur:inner-surface imat:inner-material ggeo/GBndLib collects all unique boundaries GBnd for the entire geometry, boundary indices for all bits of the geometry are stored with the geometry interleaves the standardized material and surface property arrays from GMaterialLib and GSurfaceLib into the boundary dynamic array optixrap/OBndLib converts GBndLib into the GPU boundary texture and uploads its content to GPU GPU texture usage -------------------------- optixrap/cu/boundary_lookup.h static __device__ __inline__ float4 boundary_lookup( float nm, unsigned int line, unsigned int k) optixrap/cu/state.h fill_state does the boundary lookups All pieces of geometry have a boundary index assigned to it, which means that on intersection with a ray you get the boundary index which via offsets into the boundary texture enable you to access the material and surface properties omat/osur/isur/imat relevant to that intersect. This approach copies material and surface properties multiple times, but the boundary texture is very small compared to what GPUs are designed to handle so its not a problem and it simplifies property lookup. :: 23 struct State 24 { 25 unsigned int flag ; 26 float4 material1 ; // refractive_index/absorption_length/scattering_length/reemission_prob 27 float4 m1group2 ; // group_velocity/spare1/spare2/spare3 28 float4 material2 ; 29 float4 surface ; // detect/absorb/reflect_specular/reflect_diffuse 30 float3 surface_normal ; 31 float cos_theta ; 32 float distance_to_boundary ; 33 uint4 optical ; // x/y/z/w index/type/finish/value 34 uint4 index ; // indices of m1/m2/surf/sensor 35 uint4 identity ; // node/mesh/boundary/sensor indices of last intersection 36 float ureflectcheat ; 37 }; 38 .. 48 __device__ void fill_state( State& s, int boundary, uint4 identity, float wavelength ) 49 { 50 // boundary : 1 based code, signed by cos_theta of photon direction to outward geometric normal 51 // >0 outward going photon 52 // <0 inward going photon 53 // 54 // NB the line is above the details of the payload (ie how many float4 per matsur) 55 // it is just 56 // boundaryIndex*4 + 0/1/2/3 for OMAT/OSUR/ISUR/IMAT 57 // 58 59 int line = boundary > 0 ? (boundary - 1)*BOUNDARY_NUM_MATSUR : (-boundary - 1)*BOUNDARY_NUM_MATSUR ; 60 61 // pick relevant lines depening on boundary sign, ie photon direction relative to normal 62 // 63 int m1_line = boundary > 0 ? line + IMAT : line + OMAT ; 64 int m2_line = boundary > 0 ? line + OMAT : line + IMAT ; 65 int su_line = boundary > 0 ? line + ISUR : line + OSUR ; 66 67 // consider photons arriving at PMT cathode surface 68 // geometry normals are expected to be out of the PMT 69 // 70 // boundary sign will be -ve : so line+3 outer-surface is the relevant one 71 72 s.material1 = boundary_lookup( wavelength, m1_line, 0); 73 s.m1group2 = boundary_lookup( wavelength, m1_line, 1); 74 75 s.material2 = boundary_lookup( wavelength, m2_line, 0); 76 s.surface = boundary_lookup( wavelength, su_line, 0); 77 78 s.optical = optical_buffer[su_line] ; // index/type/finish/value 79 optixrap/cu/generate.cu fill_state is the first thing down after an OptiX ray trace intersection :: 564 rtTrace(top_object, optix::make_Ray(p.position, p.direction, propagate_ray_type, propagate_epsilon, RT_DEFAULT_MAX), prd ); 565 566 if(prd.boundary == 0) 567 { 568 s.flag = MISS ; // overwrite CERENKOV/SCINTILLATION for the no hitters 569 // zero out no-hitters to avoid leftovers 570 s.index.x = 0 ; 571 s.index.y = 0 ; 572 s.index.z = 0 ; 573 s.index.w = 0 ; 574 break ; 575 } 576 // initial and CONTINUE-ing records 577 578 // use boundary index at intersection point to do optical constant + material/surface property lookups 579 fill_state(s, prd.boundary, prd.identity, p.wavelength ); optixrap/cu/wavelength_lookup.h reemission_lookup is used depending on reemission_prob and a random throw, note use of the state s.material1.w optixrap/cu/propagate.h :: 122 123 if (absorption_distance <= scattering_distance) 124 { 125 if (absorption_distance <= s.distance_to_boundary) 126 { 127 p.time += absorption_distance/speed ; 128 p.position += absorption_distance*p.direction; 129 130 const float& reemission_prob = s.material1.w ; 131 float u_reemit = reemission_prob == 0.f ? 2.f : curand_uniform(&rng); // avoid consumption at absorption when not scintillator 132 133 if (u_reemit < reemission_prob) 134 { 135 // no materialIndex input to reemission_lookup as both scintillators share same CDF 136 // non-scintillators have zero reemission_prob 137 p.wavelength = reemission_lookup(curand_uniform(&rng)); 138 p.direction = uniform_sphere(&rng); 139 p.polarization = normalize(cross(uniform_sphere(&rng), p.direction)); 140 p.flags.i.x = 0 ; // no-boundary-yet for new direction 141 142 s.flag = BULK_REEMIT ; 143 return CONTINUE; 144 } 145 else 146 { 147 s.flag = BULK_ABSORB ; 148 return BREAK; 149 } 150 } 151 // otherwise sail to boundary 152 } 153 else