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::

    201 guint4 GVolume::getIdentity()
    202 {
    203     unsigned node_index = m_index ;
    204 
    205     //unsigned identity_index = getSensorSurfaceIndex() ;   
    206     unsigned identity_index = m_copyNumber  ;
    207 
    208     return guint4(
    209                    node_index,
    210                    getMeshIndex(),
    211                    m_boundary,
    212                    identity_index
    213                  );
    214 }



* 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