Opticks : Optical Photon Simulation for Particle Physics with NVIDIA OptiX

Opticks Benefits

Outline

• Context : Photomultiplier Tube based detectors
  • Photomultiplier tubes (PMTs)
  • Supernova Neutrino Burst Observation
  • Super-Kamiokande
  • Daya Bay, JUNO
    • Liquid Scintillator measurements, modelling
• Current simulation approach, Geant4
• How GPU acceleration can help
• Opticks : optical photon simulation with NVIDIA OptiX
  • geometry migration
  • porting optical physics
• Visualization
• Validation, performance comparisons
• Summary
• Short video

Photomultiplier Tube based Detectors : Primary Tool of Neutrino Physics

Understanding optical photon generation and propagation : critical for design, operation, analysis

Radical simulation speedup -> short development cycle -> improved understanding

Photomultiplier Tubes (PMTs)

Photomultiplier Tube Operation

Single photoelectron amplified by ~10 dynodes to yield measurable signal

Kamiokande II 1

Artists impression of SN1987A ejecta : based on ESO's VLT observation

© ESO/L. Calçada - Artwork based on ESO's VLT observations - http://www.eso.org/public/news/eso1032/

Kamiokande II 3

© ESO/L. Calçada - Artwork based on ESO's VLT observations - http://www.eso.org/public/news/eso1032/

Super-Kamiokande PMTs 1

© Kamioka Observatory, ICRR(Institute for Cosmic Ray Research), The University of Tokyo.

Super-Kamiokande PMTs 4

Super-Kamiokande Expt

• 50kTon Water Cherenkov Detector
• depth 2700 mwe, 42m high, 39.2 m diameter
• 13,000 20 inch PMTs

Solar/Atmospheric/Acclerator(T2K) neutrinos

e/mu Neutrinos from detected e/mu

• cone -> ring on PMTs -> l/Neutrino direction
• e-like : e+/e- EM shower smears ring
• mu-like : sharper ring

© Kamioka Observatory, ICRR(Institute for Cosmic Ray Research), The University of Tokyo.

Super-Kamiokande PMTs 3

© Kamioka Observatory, ICRR(Institute for Cosmic Ray Research), The University of Tokyo.

Super-Kamiokande PMTs 2

© Kamioka Observatory, ICRR(Institute for Cosmic Ray Research), The University of Tokyo.

Super-Kamiokande

Neutrino zenith angle : cos(th) < 0 upward

• (e-like) upward ~ downward
• (mu-like) upward ~ 0.5 downward

1998: Muon-neutrinos from opposite side of Earth disappearing : 1st evidence neutrinos have mass

2015 Nobel : T.Kajita (SuperK), A.McDonald (SNO)

Dayabay Reactor Neutrino Expt, Far Site

Daya Bay Far Site 2

Daya Bay Far Site 3

3 Sites, 8 Identical ADs

• Far site baseline ~ 1.9km
• Relative Far/Near measurement : Antineutrino Disappearance

Daya Bay PMT Wall Photo 1

© 2010 The Regents of the University of California, through the Lawrence Berkeley National Laboratory

Daya Bay Antineutrino Detection via Inverse Beta Decay 1

Daya Bay Antineutrino Detection via Inverse Beta Decay 2

Daya Bay nGd Analysis : Most Precise Theta13

Daya Bay : E(Prompt) -> Detector Unfolded Reactor Antineutrino Spectrum

• 217 days, 6 AD data : IBDs near(far) 296,721(41,580)
• Discrepancy : Huber-Mueller(ILL-Vogel) 2.6(2.4) sigma

Jiangmen Underground Neutrino Observatory (JUNO)

JUNO : Neutrino Mass Heirarchy (MH) Determination

Need Energy Resolution : 3% (at 1 MeV)

As many photons as possible:

• high light yield scintillator ~10k photons/MeV
• long attenuation length
• high photocathode coverage ~75% (large 20 inch PMTs)
• high QE PMTs ~35%

Uniformity, Low noise

• spherical design
• clean materials, quiet PMTs, non-scintillating buffer

Baseline Numbers: arXiv:1507.05613v2 (Oct 2015) Neutrino Physics with JUNO

JUNO : A Multi-purpose Neutrino Observatory

arXiv:1507.05613v2 (Oct 2015) Neutrino Physics with JUNO

Neutrino Physics with JUNO

~220 page "Yellow Book" compendium of JUNO Physics http://arxiv.org/abs/1507.05613

• Identifying the Neutrino Mass Hierarchy
• Precision Measurements of Neutrinos
• Supernova Burst Neutrinos
• Diffuse Supernova Neutrino Background (DSNB)
• Solar Neutrinos
• Atmospheric Neutrinos
• Geoneutrinos
• Sterile Neutrinos
• Nucleon Decays
• Neutrinos from Dark Matter
• Exotic Searches with Neutrinos

Liquid Scintillator (LS) Particle Dependent Non-Linear Light Yield

Electron Energy Response

Compton scattering benchtop measurement. Electron energy response of Daya Bay LS.

Chinese Physics C Vol. 39, No. 1 (2015) 016003 ; Zhang Fei-Hong et al.

Intrinsic scintillator quenching, described by Birks' Law

Less light at high dE/dx ; kB Birks constant

Energy model considerations:

• Cherenkov light emission
• absorption and reemission of photons
• energy transfers between solvent and solutes
• readout electronics non-linearity
• calibration from multiple locations
  • multiple radioactive gamma sources
  • beta spectrum of cosmogenic boron-12

Liquid Scintillator LS ; Ternary (LAB/PPO/bis-MSB) Optimization

Emission Spectra

Dis-entangle ternary scintillation stages:

• front face measurements : no self-absorption
• thin LS sample : only non-radiative
• PPO peak 360 nm -> excitation of bis-MSB
• Dilute PPO = 3g/L PPO : no dimerization
• bis-MSB shape not in LS : non-radiative negligible

Chinese Physics C 2010 34 (11): 1724-1728 : Xiao H.L. et al

Requirements:

• high light yield, long attenuation length
• emission spectrum matched to QE of PMTs
• stable, safe, low-radioactivity, low-cost liquids

Solvent

• LAB (Linear Alkylbenzene, 10-13 carbons + benzene)
• absorbs most ionization

Primary Solute, 3000 mg/L (Fluor)

• PPO ( 2,5-diphenyloxaxazole)
• primary is Fluor to maximize scintillation yield
• excitation LAB -> PPO, non-radiatively
• excitation quenches at high density, Birks' Law (empirical)

Secondary Solute, 15 mg/L (Wavelength Shifter)

• bis-MSB (p-bis(-o-methylstyryl)-benzene)
• wavelength shift : reduce self absorption, match PMT
• excitation PPO -> bis-MSB, radiatively:absorption+reemission

Fluor Scintillation + Cerenkov light -> Absorption + Reemission

Geant4 : Monte Carlo Simulation Toolkit

Geant4 : Monte Carlo Simulation Toolkit Generality

But, there is a problem...

JPMT Before Contact 2

Ray Traced Realistic Image Synthesis and Optical Photon Simulation

OptiX Glass Sample App

http://on-demand.gputechconf.com/siggraph/2013/presentation/SG3106-Building-Ray-Tracing-Applications-OptiX.pdf

• parallel problems, just different focus
• same limit : geometry intersection

Fast ray tracing can revolutionize optical photon simulation

NVIDIA OptiX Ray Tracing Engine:

• state-of-the-art GPU accelerated intersection
• regular releases: improvements, new GPU tuning
• shared C++/CUDA context eases development
• NVIDIA expertise on efficient GPU/multi-GPU usage
  • ~linear scaling with cores across multiple GPUs

NVIDIA OptiX makes massively parallel intersection accessible

NVIDIA OptiX 1

NVIDIA OptiX 2

Brief History of GPU Optical Photon Simulation Development

2013(Aug-) [liberate geometry from Geant4]

• Develop G4DAE geometry exporter that writes tesselated COLLADA 3D files, including all material and surface properties.

2014 [integrate geometry/event data with Chroma [1]]

• Integrate G4DAE geometries with Chroma
• Develop runtime event data bridge
• GPU Cerenkov/Scintillation photon generation

2015 [replace Chroma with Opticks]

• Develop Opticks with NVIDIA OptiX ray tracing
• Opticks/Geant4 matched for simple geometries
• Speedup factor of 200x with a mobile GPU

2016 [Opticks validation against Geant4]

• Achieved match for single PMT

[1] http://chroma.bitbucket.org Developed by Stan Seibert, University of Pennsylvania.

Migrating Geometry from Geant4 to Opticks

Translate Geant4 CSG tree solids into Opticks boundaries:

• inner material : self
• outer material : parent
• inner surface : outwards going photons (self to parent)
• outer surface : inwards going photons (parent to self)

Properties as function of wavelength interleaved into float4 texture.

Texture lookup simplifies OptiX programs

// material property lookup from OptiX programs
float nmi = (nm - wavelength_domain.x)/wavelength_domain.z + 0.5f ;
float4 material1 = tex2D(boundary_texture, nmi, line + 0.5f );

float refractive_index = material1.x ;
float absorption_length = material1.y ;
float scattering_length = material1.z ;
float reemission_probability = material1.w ;  // >0 for scintillators

Large Geometry/Event Techniques

Geometry analysed to find instances

JUNO: ~90M --> 0.1M triangles

• 19,000 20" PMTs
• 34,000 3" PMTs

OpenGL/OptiX instancing

Multiple Renderers

• global, full instances, bbox instances
• rapid switching controls

Split simulation from visualization

• compute mode
  • pure OptiX buffers ~7X interop
• visualize by separate event loading

Geometry Modelling : Tesselated vs Analytic Photomultiplier Tubes

Analytic : more realistic, faster, less memory, much more effort

For Dayabay PMT:

• partition CSG solids into 12 single primitive parts (instead of 2928 triangles)
• splitting at geometrical intersections avoids implementing general CSG boolean handling
• geometry provided to OptiX in form of ray intersection code

Tesselated Photomultiplier : unrealistic disco ball effect

Analytic PMTs together with triangulated geometry

Starting from convenient tesselated geometry, implement analytic geometry for critical parts in primary optical path, to allow very close Geant4/Opticks match.

Visualizing An Optical Photon Simulation

Compressed Steps, OpenGL Buffer, written by OptiX

• space for 16 steps per photon
  • position, time, polarization, wavelength

OpenGL Geometry Shader: lines -> points/line_strip

• event time is uniform input
• pair validity from step times
• pair selection from step times relative to input time
• interpolation between steps gives position

Interactive: time scrubbing + history/material selection

// geom.glsl
layout (lines) in;
layout (line_strip, max_vertices = 2) out;
void main ()
{
    vec4 p0 = gl_in[0].gl_Position  ;   // p0.w : step begin time
    vec4 p1 = gl_in[1].gl_Position  ;   // p1.w : step end time
    ...

Compare Opticks/Geant4 Simulations with Simple Lights/Geometries

1M Photons -> Water Sphere (S-Polarized)

0.5M Photons -> Dayabay PMT

• record steps of all photons : position, time, wavelength, polarization vector, material/history codes
• 128 bit per step (2*short4) : snorm/uchar compressed using known domains
• water sphere : yields multiple rainbows depending on reflections

Opticks/Geant4 Rainbow Step Sequence Comparison

Statistically consistent photon histories in the two simulations

• BT/BR: boundary transmit/reflect
• TO/SC/SA: torch/scatter/surface absorb

 64-bit uint  Opticks    Geant4    chi2                                      (tag:5,-5)

        8ccd   819160    819654    0.15  [4 ] TO BT BT SA                    (cross droplet)
         8bd   102087    101615    1.09  [3 ] TO BR SA                       (external reflect)
       8cbcd    61869     61890    0.00  [5 ] TO BT BR BT SA                 (bow 1)
      8cbbcd     9618      9577    0.09  [6 ] TO BT BR BR BT SA              (bow 2)
     8cbbbcd     2604      2687    1.30  [7 ] TO BT BR BR BR BT SA           (bow 3)
    8cbbbbcd     1056      1030    0.32  [8 ] TO BT BR BR BR BR BT SA        (bow 4)
       86ccd     1014      1000    0.10  [5 ] TO BT BT SC SA
   8cbbbbbcd      472       516    1.96  [9 ] TO BT BR BR BR BR BR BT SA     (bow 5)
         86d      498       473    0.64  [3 ] TO SC SA
  bbbbbbbbcd      304       294    0.17  [10] TO BT BR BR BR BR BR BR BR BR  (bow 8+ truncated)
  8cbbbbbbcd      272       247    1.20  [10] TO BT BR BR BR BR BR BR BT SA  (bow 6)
  cbbbbbbbcd      183       161    1.41  [10] TO BT BR BR BR BR BR BR BR BT  (bow 7 truncated)

1M 64 bit uint flag sequences indexed using CUDA Thrust, 0.040 s (sparse histogram sort)

Rainbow Spectrum for 1st six bows

Spectra obtained by selecting photons by internal reflection counts. Colors obtained from spectra of each bin using CIEXYZ weighting functions converted into sRGB/D65 colorspace. Exposures by normalizing to bin with maximum luminance (CIE-Y) of each bow. White lines indicate geometric optics prediction of deviation angle ranges of the visible range 380-780nm. 180-360 degrees signifies exit on same side of droplet as incidence.

1M Rainbow S-Polarized, Comparison Opticks/Geant4

Deviation angle(degrees) of 1M parallel monochromatic photons in disc shaped beam incident on water sphere. Numbered bands are visible range expectations of first 11 rainbows. S-Polarized intersection (E field perpendicular to plane of incidence) arranged by directing polarization radially.

PmtInBox after 1

PmtInBox after 2

PMT Opticks/Geant4 step distribution comparison TO BT [SD]

Good agreement reached, after several fixes: geometry, total internal reflection, group velocity

position(xyz), time(t)

polarization(abc), radius(r)

PMT Opticks/Geant4 step distribution comparison : chi2/ndf

Photon Propagation Times Geant4 cf Opticks

• Opticks > 200X Geant4 with only 384 core mobile GPU[1] (multi-GPU workstation up to 20x more cores)
• photon propagation time will become effectively zero

• Interop uses OpenGL buffers allowing visualization, Compute uses OptiX buffers
• Interop/Compute : perfectly identical results, monitored by digest

[1] MacBook Pro (2013), NVIDIA GeForce GT 750M, 2048 MB, 384 cores

Summary

Opticks enables particle physics to benefit from optical photon simulation taking effectively zero time and zero CPU memory, thanks to massive parallelism made accessible by NVIDIA OptiX.

• The more photons the bigger the overall speedup (99% -> 100x)
• Drastic speedup -> better detector understanding -> greater precision
• Large PMT based neutrino experiments, such as JUNO, can benefit the most

https://bitbucket.org/simoncblyth/opticks

YouTube Video

https://www.youtube.com/watch?v=QzH6y0pKXk4

YouKu Video

http://v.youku.com/v_show/id_XMTUxNjM3MzA4OA==.html

JPMT Wide

http://simoncblyth.bitbucket.io/env/graphics/ggeoview/jpmt-wide_half.png

JPMT After Contact

http://simoncblyth.bitbucket.io/env/graphics/ggeoview/jpmt-after-contact_half.png

JPMT Approach

http://simoncblyth.bitbucket.io/env/graphics/ggeoview/jpmt-approach_half.png

JPMT Arrival

http://simoncblyth.bitbucket.io/env/graphics/ggeoview/jpmt-arrival_half.png

JPMT Inside Outside

http://simoncblyth.bitbucket.io/env/graphics/ggeoview/jpmt-inside-outside_half.png

OptiX Experience

Geometry Issues Fixed/Avoidable

• triangulated geometry leaks when shoot millions of photons at cracks (more a problem for testing than usage)
• bounding boxes that touch geometry cause leaks, slightly enlarging the boxes avoids leak (<1 in 3M)

Use cuRAND in OptiX without large stack

Init in separate CUDA launch, persist cuRAND state, subsequently load into OptiX context : avoids large stack performance hit

Compositing OptiX Ray Trace and OpenGL Rasterized

• best of both worlds : rasterization speed, raytrace accuracy
• OptiX : replace uchar4 alpha with Z clip depth
• OpenGL : fragment shader, set gl_FragDepth

out vec4 frag_colour;
uniform sampler2D ColorTex ;
void main ()
{
   frag_colour = texture(ColorTex, texcoord);
   float depth = frag_colour.w ;
   frag_colour.w = 1.0 ;
   gl_FragDepth = depth  ;
   ...