Opticks Event Data¶
Array components of Opticks Events¶
Opticks Events (optickscore/OpticksEvent.hh) comprise multiple arrays.
- nopstep
- non-optical steps
- genstep
- scintillation or cerenkov or artificial “torch” gensteps
- records
- photon step records
- photons
- last photon step at absorption, detection, truncation
- sequence
- photon level material/flag sequence histories (seqmat/seqhis), 64 bit ints holding up to 16 steps of 4-bit flags
- phosel
- photon history sequence indices, obtained by indexing sequence
- recsel
- obtained by repeating phosel by maxrec
- hits
- same structure as photons array, the subset of the photons array that passes hitmask selection
Photons¶
Optical photon generation and propagation through to absorption or detection or truncation after a configurable maximum number of bounces are performed on the GPU using the NVIDIA OptiX ray generation program, optixrap/cu/generate.cu.
The Photon struct p is serialized into the photon_buffer with psave.
optixrap/cu/generate.cu:
654 psave(p, photon_buffer, photon_offset );
optixrap/cu/photon.h:
41 struct Photon
42 {
43 float3 position ;
44 float time ;
45
46 float3 direction ;
47 float weight ;
48
49 float3 polarization ;
50 float wavelength ;
51
52 quad flags ;
53
54 // flags.i.x : 1-based signed boundary index, 0 means no intersect
55 // flags.u.y : sensor index (TODO: revive)
56 // flags.u.z : 4 debug bytes
57 // flags.u.w : history mask (bitwise OR of all step flags)
58
59 };
..
90 __device__ void psave( Photon& p, optix::buffer<float4>& pbuffer, unsigned int photon_offset)
91 {
92 pbuffer[photon_offset+0] = make_float4( p.position.x, p.position.y, p.position.z, p.time );
93 pbuffer[photon_offset+1] = make_float4( p.direction.x, p.direction.y, p.direction.z, p.weight );
94 pbuffer[photon_offset+2] = make_float4( p.polarization.x,p.polarization.y,p.polarization.z, p.wavelength );
95 pbuffer[photon_offset+3] = make_float4( p.flags.f.x, p.flags.f.y, p.flags.f.z, p.flags.f.w);
96 }
The quad type is a union of float4, uint4 and int4 types allowing integers or unsigned integers to be conveniently carried and persisted within the float buffer/array.
Cerenkov and scintillation photons and artificial “torch” photons are generated by methods from the below headers. The 4 flags values are initialized to zero at generation.
- optixrap/cu/cerenkovstep.h
- optixrap/cu/scintillationstep.h
- optixrap/cu/torchstep.h
Subsequent code that changes p.flags:
- optixrap/cu/generate.cu
- optixrap/cu/photon.h
- optixrap/cu/propagate.h
optixrap/cu/generate.cu:
201 /**
202 FLAGS
203 -------
204
205 Set the photon flags p.flags using values from state s and per-ray-data prd
206
207 **/
208
209 #define FLAGS(p, s, prd) \
210 { \
211 p.flags.i.x = prd.boundary ; \
212 p.flags.u.y = s.identity.w ; \
213 p.flags.u.w |= s.flag ; \
214 } \
...
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 ;
...
574 break ;
575 }
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 );
580
581 s.distance_to_boundary = prd.distance_to_boundary ;
582 s.surface_normal = prd.surface_normal ;
583 s.cos_theta = prd.cos_theta ;
584
585 // setting p.flags for things like boundary, history flags
586 FLAGS(p, s, prd);
The closest hit program closest_hit_propagate populates prd with information from the geometry. Including:
- prd.distance_to_boundary (float)
- t distance in ray direction from ray origin to closest hit intersection point, this is closest solution of the polynomial in t that fulfils the minimum t requirement
- prd.identity (uint4)
- geometry identity information including GVolume/GNode index, mesh/solid index, boundary index, sensor index
ggeo/GVolume.cc:
259 glm::uvec4 GVolume::getIdentity() const
260 {
261 glm::uvec4 id(getIndex(), getTripletIdentity(), getShapeIdentity(), getSensorIndex()) ;
262 return id ;
263 }
- Opticks GVolume/GNode corresponds to Geant4 placed volumes (ie Physical Volume, G4PVPlacement). There are often many tens of thousands of these in the detector geometry
- There are typically much fewer distinct GMesh in a geometry, these correspond to Geant4 solids that are associated with a logical volume. Despite the “mesh” name GMesh instances contain both triangulated and analytic geometry information.
The NVIDIA OptiX ray tracing pipeline uses its acceleration structure to determine which geometry intersects to call. Subsequently the closest hit program is called. This populates the PerRayData_propagate prd struct in order to communicate geometry and material/surface information back to the ray generation program.
optixrap/cu/material1_propagate.cu:
20 #include <optix.h>
21 #include "PerRayData_propagate.h"
22 #include "wavelength_lookup.h"
23
//
// attributes set in the geometry specific intersect programs
//
26 rtDeclareVariable(float3, geometricNormal, attribute geometric_normal, );
27 rtDeclareVariable(uint4, instanceIdentity, attribute instance_identity, );
28
29 rtDeclareVariable(PerRayData_propagate, prd, rtPayload, );
30 rtDeclareVariable(optix::Ray, ray, rtCurrentRay, );
31 rtDeclareVariable(float, t, rtIntersectionDistance, );
32
33
34 RT_PROGRAM void closest_hit_propagate()
35 {
36 const float3 n = normalize(rtTransformNormal(RT_OBJECT_TO_WORLD, geometricNormal)) ;
37 float cos_theta = dot(n,ray.direction);
38
39 prd.cos_theta = cos_theta ;
40 prd.distance_to_boundary = t ; // huh: there is an standard attrib for this
41 unsigned int boundaryIndex = instanceIdentity.z ;
42 prd.boundary = cos_theta < 0.f ? -(boundaryIndex + 1) : boundaryIndex + 1 ;
43 prd.identity = instanceIdentity ;
44 prd.surface_normal = cos_theta > 0.f ? -n : n ;
45 }
optixrap/cu/intersect_analytic.cu:
067 rtDeclareVariable(unsigned int, instance_index, ,);
68 // optix::GeometryInstance instance_index into the identity buffer,
69 // set by oxrap/OGeo.cc, 0 for non-instanced
...
74 rtBuffer<Part> partBuffer;
75 rtBuffer<Matrix4x4> tranBuffer;
76
77 rtBuffer<Prim> primBuffer;
78 rtBuffer<uint4> identityBuffer; // from GMergedMesh::getAnalyticInstanceIdentityBuffer()
..
84 // attributes communicate between intersect (which is geometry specific) and the general closest hit program,
85 // the attributes must be set inbetween rtPotentialIntersection and rtReportIntersection in the intersect program
86
87 rtDeclareVariable(uint4, instanceIdentity, attribute instance_identity,);
88 rtDeclareVariable(float3, geometric_normal, attribute geometric_normal, );
89 rtDeclareVariable(float3, shading_normal, attribute shading_normal, );
90
...
287 RT_PROGRAM void intersect(int primIdx)
288 {
289 const Prim& prim = primBuffer[primIdx];
290
291 unsigned partOffset = prim.partOffset() ;
292 unsigned numParts = prim.numParts() ;
293 unsigned primFlag = prim.primFlag() ;
294
295 uint4 identity = identityBuffer[instance_index] ;
296
297 if(primFlag == CSG_FLAGNODETREE)
298 {
299 Part pt0 = partBuffer[partOffset + 0] ;
300
301 identity.z = pt0.boundary() ; // replace placeholder zero with test analytic geometry root node boundary
302
303 evaluative_csg( prim, identity );
306 }
optixrap/cu/csg_intersect_boolean.h:
563 static __device__
564 void evaluative_csg( const Prim& prim, const uint4& identity )
565 {
...
... CSG boolean intersection
...
882 if(rtPotentialIntersection( fabsf(ret.w) ))
883 {
884 shading_normal = geometric_normal = make_float3(ret.x, ret.y, ret.z) ;
885 instanceIdentity = identity ;
...
892 rtReportIntersection(0);
893 }
894 }
...
928 }
- identity (uint4) for a specific geometry instance (eg a particular PMT) is obtained from the AnalyticInstanceIdentityBuffer with the instance_index
- NVIDIA OptiX attributes are used to communicate between the intersect program (which is geometry specific) and the general closest hit program
The FLAGS macro sets p.flags values obtained from from the s and prd structs.
- prd : PerRayData_propagate struct
- used to communicate between the ray trace intersect program and ray generation program
- s : State struct
- filled with geometry and material/surface information by optixrap/cu/state.h:fill_state
Photon flags¶
- p.flags.i.x
- 1-based signed integer boundary index of last intersect, 0 means no intersect, sign encodes in/out direction relative to the surface normal of the geometry at the intersect
- p.flags.u.y
- (s.identity.w) sensor index (TODO: revive this)
- p.flags.u.z
- 4 debugging bytes : 3 of which are currently spare
- p.flags.u.w
- history mask created from the bitwise-or of s.flag (state flags) for all steps of the propagation
Note that in addition to the history mask Opticks also records the history sequence into the sequence array. The so called seqhis and seqmat record 4-bit flags from up to 16 steps of the photon being into a 64 bit integer providing the step by step flag and material history of the photon propagation.
- optickscore/OpticksPhoton.h
enum of possible flags:
22 enum 23 { 24 CERENKOV = 0x1 << 0, 25 SCINTILLATION = 0x1 << 1, 26 MISS = 0x1 << 2, 27 BULK_ABSORB = 0x1 << 3, 28 BULK_REEMIT = 0x1 << 4, 29 BULK_SCATTER = 0x1 << 5, 30 SURFACE_DETECT = 0x1 << 6, 31 SURFACE_ABSORB = 0x1 << 7, 32 SURFACE_DREFLECT = 0x1 << 8, 33 SURFACE_SREFLECT = 0x1 << 9, 34 BOUNDARY_REFLECT = 0x1 << 10, 35 BOUNDARY_TRANSMIT = 0x1 << 11, 36 TORCH = 0x1 << 12, ...
- optickscore/OpticksFlags.{hh,cc}
- flag manipulation/presentation methods
- ana/base.py ana/hismask.py
- python machinery for flag manipulation
In [1]: from opticks.ana.base import PhotonMaskFlags
In [2]: flags = PhotonMaskFlags()
[2020-06-17 12:38:14,608] p69482 {__init__ :enum.py :37} INFO - parsing $OPTICKS_PREFIX/include/OpticksCore/OpticksPhoton.h
[2020-06-17 12:38:14,608] p69482 {__init__ :enum.py :39} INFO - path expands to /usr/local/opticks/include/OpticksCore/OpticksPhoton.h
In [8]: flags.code2name
Out[8]:
{1: 'CERENKOV',
2: 'SCINTILLATION',
4: 'MISS',
8: 'BULK_ABSORB',
16: 'BULK_REEMIT',
32: 'BULK_SCATTER',
64: 'SURFACE_DETECT',
128: 'SURFACE_ABSORB',
256: 'SURFACE_DREFLECT',
512: 'SURFACE_SREFLECT',
1024: 'BOUNDARY_REFLECT',
2048: 'BOUNDARY_TRANSMIT',
4096: 'TORCH',
8192: 'NAN_ABORT',
16384: 'G4GUN',
32768: 'FABRICATED',
65536: 'NATURAL',
131072: 'MACHINERY',
262144: 'EMITSOURCE',
524288: 'PRIMARYSOURCE',
1048576: 'GENSTEPSOURCE'}
Higher level flag and abbreviation manipulations are provided by ana/hismask.py:
epsilon:ana blyth$ ipython
In [1]: from opticks.ana.hismask import HisMask
In [2]: hm = HisMask()
[2020-06-17 12:53:07,046] p69692 {__init__ :enum.py :37} INFO - parsing $OPTICKS_PREFIX/include/OpticksCore/OpticksPhoton.h
[2020-06-17 12:53:07,046] p69692 {__init__ :enum.py :39} INFO - path expands to /usr/local/opticks/include/OpticksCore/OpticksPhoton.h
In [3]: hm.label(65)
Out[3]: 'SD|CK'
In [4]: hm.label(129)
Out[4]: 'SA|CK'
In [5]: hm.abbrev.abbr2name
Out[5]:
{'/usr/local/opticks/include/OpticksCore/OpticksFlags_Abbrev.json': 'jsonLoadPath',
'AB': 'BULK_ABSORB',
'BR': 'BOUNDARY_REFLECT',
'BT': 'BOUNDARY_TRANSMIT',
'CK': 'CERENKOV',
'DR': 'SURFACE_DREFLECT',
'FD': 'FABRICATED',
'GN': 'G4GUN',
'GS': 'GENSTEPSOURCE',
'MI': 'MISS',
'MY': 'MACHINERY',
'NA': 'NAN_ABORT',
'NL': 'NATURAL',
'PS': 'PRIMARYSOURCE',
'RE': 'BULK_REEMIT',
'SA': 'SURFACE_ABSORB',
'SC': 'BULK_SCATTER',
'SD': 'SURFACE_DETECT',
'SI': 'SCINTILLATION',
'SO': 'EMITSOURCE',
'SR': 'SURFACE_SREFLECT',
'TO': 'TORCH',
'XX': 'BAD_FLAG'}
Opticks analysis and debugging is based on python and NumPy (plus matplotlib for plotting). You may need to set PYTHONPATH to load the opticks module, eg with:
export PYTHONHOME=$HOME