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