JUNO+Opticks Photons : Validation and Deployment Plan

Open source, https://bitbucket.org/simoncblyth/opticks

Opticks replaces Geant4 optical photon simulation with an equivalent implementation that benefits from state-of-the-art GPU ray tracing from NVIDIA OptiX. All optically relevant aspects of Geant4 context must be translated+copied to GPU:

geometry : solids, structure, material+surface properties

generation : Cerenkov + Scintillation (using Gensteps from Geant4)

propagation : Rayleigh scattering, absorption, reemission, boundary

Achieving+maintaining equivalence is time consuming, BUT:

transformative performance benefits >1000x Geant4

Simon C Blyth, IHEP, CAS — 18 Jan 2022, JUNO Collaboration Meeting

Geant4OpticksWorkflow

Geant4OpticksWorkflow 2

Geometry/Translation Issues Summary

Fastener : interfering subtraction-of-subtraction issue

--additionacrylic-simplify-csg

TODO: further checks before making it default
could add "cavity" volume to avoid geo change

PMTs do not CSG subtract cavity in Pyrex for Vacuum...
favor volume heirarchy over complicated CSG

Cutdown PMT breaks Opticks translation

profligate G4IntersectionSolid cutdowns => overlarge CSG tree
FIXED with ZSolid : actually cuts the CSG tree
expect significant Opticks + Geant4 speedups

solidXJfixture : ~10/64 overlaps with fasteners

WIP : fastener/fixture positions need to work together

solidXJfixture : extreme coincidence, many spurious intersects

extreme coincident face issue, difficult to fix at Geant4 level
TODO: re-implement simpler, fewer constituents

Other Progress Summary

GPU Cerenkov float/double investigations

better precision numerical integration using S2 "sin^2cone_angle"
developed inverse transform lookup based Cerenkov wavelength generation
- faster than rejection sampling
- copes better with float precision than rejection sampling

Outline

2D Opticks Geometry Slicing
PMT Mask Modelling Fix
Fastener : interfering subtraction-of-subtraction issue
2D Geant4 Geometry Slicing
Cutdown PMT breaks Opticks translation
Render Speed Check
History Matching Check
Investigate XJfixtureConstruction Solid
Investigate XJfixtureConstruction Positions
Investigate XJanchorConstruction Solid
Investigate SJReceiverConstruction Solid

2D Opticks Geometry Slicing

Planar 2D ray tracing to give slice renders for translated geometry debug

/env/presentation/QEventTest/planar_gensteps.png

New transform enabled TORCH genstep (num_gs,6,4,4)

1*int4 : config enums etc.., number of photons to generate
1*float4 : local frame position
4*float4 : 4x4 local to global transform

opticks/sysrap/SEvent::MakeCenterExtentGensteps

termed "center-extent gensteps" (CEGS)
precise local targetting : used for rendering "slices"

Transform local frame photon source positions (eg XZ plane, XYZ grid) into global frame using transforms of any piece of geometry.

used GPU side with : quadarap/qsim.h/qsim::generate_photon_torch

671 template <typename T> inline QSIM_METHOD void qsim<T>::generate_photon_torch(quad4& p, curandStateXORWOW& rng, const quad6& gs)
672 {
673     p.q0.f = gs.q1.f ;  // start with local frame position, eg (0,0,0)
674
675     float u = curand_uniform(&rng);
676     float sinPhi, cosPhi;
677     sincosf(2.f*M_PIf*u,&sinPhi,&cosPhi);
678
679     //  local frame XZ plane directions
680     p.q1.f.x = cosPhi ;  p.q1.f.y = 0.f    ;  p.q1.f.z = sinPhi ;  p.q1.f.w = 0.f ;
681
682     qat4 qt(gs) ; // copy 4x4 transform from last 4 quads of genstep
683     qt.right_multiply_inplace( p.q0.f, 1.f );   // position            transform local->global
684     qt.right_multiply_inplace( p.q1.f, 0.f );   // direction
685 }

PMT Mask Modelling Fix

Hama_1_figs_positions_mpplt_pid.png

mask tail "bowl" cutting across PMT bulb ?

Hama_1_figs_positions_mpplt_pid.png 2

PMT Mask CSG Modelling Fix

coincident face CSG subtraction issue

easily fixed by CSG modelling change
enlarge subtracted shape
- avoid coincident face
- does not change shape, just modelling
Some such issues are not so easy to fix
Right : OpenGL(pyvista) render of Geant4 polygonized mesh
Below : before/after Opticks GPU ray trace renders

Fastener

Fastener : interfering subtraction-of-subtraction issue

quicktime 2 lAddition_uni_acrylic3

StickMPL_all

two horizontal green lines just below (0,0) are spurious
- => wrong Water/Acrylic properties used
- Opticks "SI BT BT BT BT AB" excess
- CAUSE : poor G4 geometry modelling
- => interfering subtraction of subtraction

`Daughter Cavity`

=> pointless
=> not "the Geant4 way"
=> causes sub-sub issue

Fastener : interfering subtraction-of-subtraction issue

Geant4: Each volume is created by describing its shape and its physical characteristics, and then placing it inside a containing volume.

daughter material replaces that of parent
favor volume hierarchy over complex CSG
CSG inherently limited, can be expensive/buggy
simple CSG avoids problems:

CSG Subtracting Daughter Cavities

not necessary, complicates CSG
easily doubles expense for no benefit
causes this issue

Vacuum cavity not subtracted from PMT Pyrex ...: So why do that here ?

If gap between acrylic cone and steel rods is not important

Model with hierarchy of 2 volumes:

Outer Volume : Acrylic polycone (without cavity for rods)
Inner Volume : Steel Rods

This is the approach taken with option:

--additionacrylic-simplify-csg

problem is avoided, but some geometry change
currently not the default

If gap between acrylic cone and steel rods is important

Model with hierarchy of 3 volumes:

Outer Volume : Acrylic polycone
Middle Volume : Water Cylinders
Inner Volume : Steel Rods

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/geometryDef.html

How to handle Opticks CSG colocated sub-sub limitation A-(B-C)

Multiple non-interfering CSG subtractions work fine

problem arises from >1 subtractions "on top of each other",
- colocation breaks un-coinciding by subtractee expansion

CSG Sub. "pipe" cylinder => spurious inner intersects

"pipe" implemented as CSG subtraction
implementing pipe directly as primitive may fix issue
- significant work, but doable, not a general fix

Problem Arises from inherent CSG fragility regarding coincident faces

Possibility : Volume-centric CSG Implementation more like Geant4 ?

more robust (presumably)
slower, keep asking "is this point inside this constituent subtree ?"
uncertain development work to implement on GPU, eg:
- CUDA code generation of SDF (signed distance functions) for all CSG subtrees of all solids

Opticks CSG is surface-centric

more fragile but inherently faster that Geant4 implementation
identical results to Geant4 so long as coincident faces are avoided

2D Geant4 Geometry Slicing

X4IntersectSolidTest : Geant4 2D cross-sections of single G4VSolid

get G4VSolid from j/PMTSim with PMTSim::GetSolid
SEvent::GenerateCenterExtentGenstepsPhotons creates grid of 2D planar "rays"
X4Intersect::scan uses X4Intersect::Distance and collects intersect positions into NP array
X4IntersectSolidTest.py presents the intersect positions with scatter plots

G4double X4Intersect::Distance(const G4VSolid* solid,
          const G4ThreeVector& pos, const G4ThreeVector& dir)
{
    EInside in =  solid->Inside(pos) ; G4double t = kInfinity ;
    switch( in ) {
        case kInside:  t = solid->DistanceToOut( pos, dir ) ; break ;
        case kSurface: t = solid->DistanceToOut( pos, dir ) ; break ;
        case kOutside: t = solid->DistanceToIn(  pos, dir ) ; break ; }
    return t ;
}

body_solid_nurs

2D scatter plot of Geant4 intersects from planar genstep grid
spurious Geant4 cylinder-torus neck intersects visible

body_solid_nurs_pcnk

G4Polycone neck avoids the spurious intersects

body_solid

after real surface "cutdown" still have spurious G4 intersects

body_solid_pcnk

polycone neck (now default) avoids issue

X4IntersectVolumeTest : Geant4 2D cross-sections of G4VPhysicalVolume

G4VPhysicalVolume has no convenient "Distance" methods ... so scan solids individually and present together after applying structure transforms.

get G4VPhysicalVolume from j/PMTSim with PMTSim::GetPV
collect and save structure transforms PMTSim::SaveTransforms for each solid
X4Intersect::Scan each solid and persist in NP arrays
present together using X4IntersectVolumeTest.py

Gives rapid 2D cross-sections of G4PhysicalVolume (very little Opticks code used)

Usage of xxv.sh script which runs executable and ipython:

cd ~/opticks/extg4
./xxv.sh
GEOM=body_phys ./xxv.sh

body_phys_nurs_pdyn

with dynode geometry and G4Tubs-G4Torus neck

body_phys_nurs_pcnk_pdyn

with dynode geometry included and G4Polycone neck

body_phys_pdyn_pcnk

Cutdown PMT

Cutdown PMT breaks Opticks translation

Hamamatsu PMT Solid breaking Opticks : CSG tree height 8 : TOO DEEP

ZSolid::Draw [-1] nameprefix _body_solid_  NODE:19 PRIM:10 UNDEFINED:19 EXCLUDE:0 INCLUDE:0 MIXED:0 Order:IN                            Int
                                                                                                                                        U
                                                                                                                                        17

                                                                                                                        Uni                     Pol
                                                                                                                        U                       U
                                                                                                                        15                      18

                                                                                                        Uni                     Pol
                                                                                                        U                       U
                                                                                                        13                      16

                                                                                        Uni                     Pol
                                                                                        U                       U
                                                                                        11                      14

                                                                        Uni                     Ell
                                                                        U                       U
                                                                        9                       12

                                        Uni                                     Pol
                                        U                                       U
                                        5                                       10

                        Uni                             Sub
                        U                               U
                        3                               7

        Uni                     Ell             Pol             Tor
        U                       U               U               U
        1                       4               6               8

Ell             Pol
U               U
0               2

0.0             -2.5            -5.0            -179.2          -210.0          -242.5          -275.0          -385.0          -420.0          3.4     zdelta

190.0           0.0             -5.0            -148.4          -130.0          -210.0          -275.0          -350.0          -420.0          190.0   az1
0.0             -5.0            -195.0          -210.0          -290.0          -275.0          -365.0          -420.0          -450.0          -183.2  az0

I       1_2     II      1_3     III     1_4     IV_tub  IV      IV_tor  1_5     V       1_6     VI      1_8     VIII    1_9     IX      cut     intubs

ZSolid::Draw Uni:Union Ell:Ellipsoid Pol:Polycone Tor:Torus Int:Intersection Sub:Subtraction
19 Nodes, 10 prims, depth 8 Use of G4IntersectionSolid to make z-cut prevents tree balancing

Hamamatsu PMT Solid : PROFLIGATE G4ItersectionSolid Z-Cut

Bad Practice : Using G4IntersectionSolid to apply Z-cut to PMT

flexible and convenient cut solid definition
BUT EXCEEDINGLY EXPENSIVE (PROFLIGATE) : and it breaks Opticks translation
every single intersect must traverse CSG tree checking intersects with multiple cutaway constituent primitives
-> pay the price of convenient solid definition easily millions of times per event


485     // Reduce the size when real surface is enabled.
487     if (m_useRealSurface ) {
...
546         const double body_height = m_pmt_h;
547         const double body_half_height = body_height / 2;
548         const G4ThreeVector cut_body_displacement(0., 0., m_z_equator-pmt_half_height);
549         G4VSolid* cut_body_solid = new G4Tubs( GetName() + "_body_solid_intubs",
550                                               0.,
551                                               helper_sep_tube_r+1E-9*mm,
552                                               body_half_height,
553                                               0., 360.*degree);
554         body_solid = new G4IntersectionSolid( GetName() + "_body_solid_cut",
555                                                 body_solid,
556                                                 cut_body_solid,
557                                                 NULL,
558                                                 cut_body_displacement);

SOLUTION ACTUALLY CUT THE CSG TREE using https://github.com/simoncblyth/j/blob/main/PMTSim/ZSolid.hh

G4VSolid* ZSolid::ApplyZCutTree( const G4VSolid* original, double zcut)

j/PMTSim : ZSolid::ApplyZCutTree

G4VSolid* ZSolid::ApplyZCutTree(const G4VSolid* orig, double zcut)

G4BooleanSolid has no SetRight SetLeft: so placement new replace at same address

trick to re-use ctor with different args

G4BooleanSolid* src = dynamic_cast<G4BooleanSolid*>(node) ;
G4String name = src->GetName() ;
G4VSolid* left = src->GetConstituentSolid(0) ;

G4SolidStore::GetInstance()->DeRegister(src);

src1 = new (src) G4UnionSolid(name, left, new_right)
assert( src1 == src );

Similar trick used to cut G4Polycone ZSolid::ApplyZCut_G4Polycone

https://github.com/simoncblyth/j/blob/main/PMTSim/ZSolid.hh

Applying Z-cut to G4VSolid CSG Tree

classify tree nodes INCLUDE/EXCLUDE against Z-cut
identify "crux" nodes dividing INCLUDE/EXCLUDE subtrees
prune at crux nodes by tree surgery reconnections
repeat until all tree nodes are INCLUDE

Pruning and Reconnection

pruning at non-root "crux" node
- identify surviving child
- change crux parent Left/Right to the survivor
pruning at root "crux" node
- identify surviving child
- change root to the survivor

Geant4 Booleans not intended to be modified: -> required "placement new" trickery

ZSolid::ApplyZCutTree classify tree nodes against Z-cut

ZSolid::apply_cut before prune [ 0] nameprefix maker_body_solid_zcut-183.2246_  NODE:15 PRIM:8 UNDEFINED:0 EXCLUDE:4 INCLUDE:7 MIXED:4 Order:IN
                                                                                                        Uni
             I : include                                                                                IE
                                                                                                        13
             E : exclude
                                                                                        Uni                     Pol
             IE: mixed include/exclude                                                  IE                      E
                                                                                        11                      14
             X : "crux" node
                                                                        Uni                     Pol
             S : survivor node                                          IE                      E
                                                                        9                       12

                                                        Uni                     Ell
                                                        IE                      E
                                                        7                       10
                                                        X
                                        Uni                     Pol
                                        I                       E
                                        5                       8
                                        S
                        Uni                     Pol
                        I                       I
                        3                       6

        Uni                     Ell
        I                       I
        1                       4

Ell             Pol
I               I
0               2

0.0             -2.5            -5.0            -179.2          -242.5          -275.0          -385.0          -420.0  zdelta

190.0           0.0             -5.0            -162.0          -210.0          -275.0          -350.0          -420.0  az1
0.0             -5.0            -183.2          -183.2          -275.0          -365.0          -420.0          -450.0  az0

I       1_2     II      1_3     III     1_4     IV      1_5     V       1_6     VI      1_8     VIII    1_9     IX

CSG constituents are classified INCLUDE/EXCLUDE/MIXED against the Z-cut

ZSolid::ApplyZCutTree : (2) prune + reconnect tree

CSG node tree			Stage
node	leaf	height	Stage
15	8	7	Original Solid
17	9	8	After "cutting" with G4IntersectionSolid
7	4	3	After actual tree cutting ZSolid::ApplyZCutTree

same shape using 7 CSG nodes (not 17)
still flexible
avoids breaking Opticks translation
=> expect substantial Geant4 speedup

tree::apply_cut after prune and re-classify [ 3] nameprefix maker_body_solid_zcut-183.2246_
NODE:7 PRIM:4 UNDEFINED:0 EXCLUDE:0 INCLUDE:7 MIXED:0 Order:IN

                                        Uni
                                        I
                                        5
                                        S

                        Uni                     Pol
                        I                       I
                        3                       6


        Uni                     Ell
        I                       I
        1                       4


Ell             Pol
I               I
0               2


0.0             -2.5            -5.0            -179.2  zdelta

190.0           0.0             -5.0            -162.0  az1
0.0             -5.0            -183.2          -183.2  az0
I       1_2     II      1_3     III     1_4     IV

Getting JUNO to work on GPU, improves it on CPU too

maker_body_solid_zcut-500.0

maker_body_solid_zcut-400.0

maker_body_solid_zcut-350.0

maker_body_solid_zcut-300.0

maker_body_solid_zcut-200.0

maker_body_solid_zcut-183.25

maker_body_solid_zcut-100.0

maker_body_solid_zcut-0.0

maker_body_solid_zcutp100.0

Render Speed Check

cxr_overview_emm_t0_moi_-1_ALL.jpg

cxr_overview_emm_t0,_moi_-1.jpg

-e t0, : NOT 0 : 3084:sWorld : exclude global remainder volumes

cxr_overview_emm_image_grid_overview

Comparison of ray traced render times of different geometry: simple way to find issues, eg over complex CSG, overlarge BBox

sWaterTube_image_grid_cxr_view

Same viewpoint inside JUNO Central Detector, vary included volumes: ray trace performance very sensitive to geometry and its modelling => BVH structure

[Dec 2021] JUNO : OptiX 7 Ray Trace Times ~2M-pix : TITAN RTX

-e : controls components : "t" means ~ (NOT)
time(s) : GPU ray trace CUDA launch time
relative : compares to "ONLY PMT" baseline

Recent Geometry Fixes

fix profligate PMT modelling
simpler Acrylic fastener
avoid overlarge Fastener BBox
fix Polycone inner overlap

Small-ish time range 1:15 (previously 1:600): suggests no bad bottlenecks remaining

Substantial speedup after fixing geometry issues

ALL PMTs : 0.0097 -> 0.0061 (x1.6 faster)

ALL : 0.6240 -> 0.0054 (x155 faster)

idx	-e	time(s)	relative	enabled geometry description 3dbec4dc
0	5,	0.0004	0.0643	ONLY: 1:sStrutBallhead
1	9,	0.0004	0.0658	ONLY: 130:sPanel
2	7,	0.0005	0.0782	ONLY: 1:base_steel
3	8,	0.0006	0.0966	ONLY: 1:uni_acrylic1
4	6,	0.0006	0.1009	ONLY: 1:uni1
5	1,	0.0009	0.1476	ONLY: 5:PMT_3inch_pmt_solid FAST cf 20in
6	4,	0.0015	0.2386	ONLY: 4:mask_PMT_20inch_vetosMask
7	3,	0.0033	0.5373	ONLY: 5:HamamatsuR12860sMask SLOW cf 3in
8	0,	0.0040	0.6556	ONLY: 3084:sWorld
9	2,	0.0040	0.6627	ONLY: 5:NNVTMCPPMTsMask SLOW cf 3in
10	t4,	0.0050	0.8307	EXCL: 4:mask_PMT_20inch_vetosMask
11	t2,	0.0051	0.8391	EXCL: 5:NNVTMCPPMTsMask
12	t3,	0.0052	0.8514	EXCL: 5:HamamatsuR12860sMask
13	t6,	0.0053	0.8799	EXCL: 1:uni1
14	t7,	0.0054	0.8809	EXCL: 1:base_steel
15	t0	0.0054	0.8843	ALL
16	t5,	0.0054	0.8843	EXCL: 1:sStrutBallhead
17	t9,	0.0054	0.8855	EXCL: 130:sPanel
18	t1,	0.0054	0.8860	EXCL: 5:PMT_3inch_pmt_solid
19	t8,	0.0055	0.9013	EXCL: 1:uni_acrylic1
20	t0,	0.0059	0.9753	EXCL: 3084:sWorld
21	1,2,3,4	0.0061	1.0000	ONLY PMT
22	t8,0	0.0062	1.0217	EXCL: 1:uni_acrylic1 3084:sWorld

Perhaps solid angle correction would even out render speeds

History Matching Check

tds3gun.sh history check

 epsilon:ana blyth$ tds3gun.sh 1 ab.ahis .    all_seqhis_ana .           TOTALS: n        iseq 0000    7cccccc2 0001          42 0002   7cccccc62 0003   7cccccc52 0004         452 0005         462 0006          41 0007  7cccccc662 0008  7cccccc652 0009  7cccccc552 0010        4552 0011        8cc2 0012        4652 0013  cccccc6662 0014  cccccc6652 0015        4662 .        TOTALS:    11300    11300

## seqhis: 64bit uint : 16x4bit step flags for each photon cfo:sum 1:g4live:tds3gun -1:g4live:tds3gun 11300 11300 94.69 94.69/62 = 1.53 pvalue:P[C2>]:1.000 1-pvalue:P[C2<]:0.000 a b a-b (a-b)^2/(a+b) a/b b/a [ns] label 1832 1795 37 0.38 1.021 +- 0.024 0.980 +- 0.023 [8 ] SI BT BT BT BT BT BT SD 1727 1683 44 0.57 1.026 +- 0.025 0.975 +- 0.024 [2 ] SI AB 766 713 53 1.90 1.074 +- 0.039 0.931 +- 0.035 [9 ] SI SC BT BT BT BT BT BT SD 522 518 4 0.02 1.008 +- 0.044 0.992 +- 0.044 [9 ] SI RE BT BT BT BT BT BT SD 515 505 10 0.10 1.020 +- 0.045 0.981 +- 0.044 [3 ] SI RE AB 379 377 2 0.01 1.005 +- 0.052 0.995 +- 0.051 [3 ] SI SC AB 256 262 -6 0.07 0.977 +- 0.061 1.023 +- 0.063 [2 ] CK AB 246 249 -3 0.02 0.988 +- 0.063 1.012 +- 0.064 [10] SI SC SC BT BT BT BT BT BT SD 212 230 -18 0.73 0.922 +- 0.063 1.085 +- 0.072 [10] SI RE SC BT BT BT BT BT BT SD 167 168 -1 0.00 0.994 +- 0.077 1.006 +- 0.078 [10] SI RE RE BT BT BT BT BT BT SD 173 162 11 0.36 1.068 +- 0.081 0.936 +- 0.074 [4 ] SI RE RE AB 118 173 -55 10.40 0.682 +- 0.063 1.466 +- 0.111 [4 ] SI BT BT SA 119 120 -1 0.00 0.992 +- 0.091 1.008 +- 0.092 [4 ] SI RE SC AB 124 104 20 1.75 1.192 +- 0.107 0.839 +- 0.082 [10] SI SC SC SC BT BT BT BT BT BT 111 114 -3 0.04 0.974 +- 0.092 1.027 +- 0.096 [10] SI RE SC SC BT BT BT BT BT BT 100 110 -10 0.48 0.909 +- 0.091 1.100 +- 0.105 [4 ] SI SC SC AB 94.69 94.69/62 = 1.53 pvalue:P[C2>]:1.000 1-pvalue:P[C2<]:0.000 class="b">prior check excluded sticks geometry due to Fastener, giving chi2/ndf ~ 1 : (now with full geometry : ~1.5 )

Largest discrepancy is an Opticks lack of "SI BT BT SA" compared to Geant4

fraction of discrepant photons 55/11300 ~0.5%

History Check : largest Chi2 contribs

Selecting the largest chi2 contribs:

In [4]: np.set_printoptions(linewidth=100)

In [5]: ab.his.ss[ab.his.c2 > 2]
Out[5]:
array(['0011             8cc2      118     173       -55           10.40     SI BT BT SA',
       '0018        7ccccccc2       70      100      -30            5.29     SI BT BT BT BT BT BT BT SD',
       '0019       cccccccc62       77       59       18            2.38     SI SC BT BT BT BT BT BT BT BT',
       '0031       ccccc66652       35       54      -19            4.06     SI RE SC SC SC BT BT BT BT BT',
       '0042            4ccc2       54        8       46           34.13     SI BT BT BT AB',
       '0045       ccccc65552       21       36      -15            3.95     SI RE RE RE SC BT BT BT BT BT',
       '0059       cccbcccc52       24       13       11            3.27     SI RE BT BT BT BT BR BT BT BT'],
      dtype='|S98')

In [6]:

"Counterpart" Opticks excess is "SI BT BT BT AB" which happens very little for Geant4.

Fixing that one pair of discrepancies would get chi2/ndf below 1.

Looks like only one issue big enough to chase at this sample size.

Delving in => PE_PA solid is discrepant => XJfixtureConstruction

In [8]: a.sel = "SI BT BT BT AB"   ## select Opticks photons with this history

In [9]: a.bn.shape
Out[9]: (54, 1, 4)

In [13]: a.bn.view(np.int8).shape
Out[13]: (54, 1, 16)

In [14]: a.bn.view(np.int8)   ## step-by-step boundary indices within selection, suspicious -24, 24 pairs
Out[14]:
A([[[ 18, -17,  17,  24,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]],
   [[ 18, -17, -24,  24,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]],
   [[ 18, -17, -24,  24,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]],
   [[ 18, -17, -24,  24,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]],
   ...

In [26]: a.blib.bname(17-1)
Out[26]: 'Water///Acrylic'

In [27]: a.blib.bname(24-1)
Out[27]: 'Water///PE_PA'

epsilon:issues blyth$ jgr PE_PA
...
./Simulation/DetSimV2/CentralDetector/include/XJfixtureConstruction.hh:    G4Material* PE_PA ;
./Simulation/DetSimV2/CentralDetector/src/XJfixtureConstruction.cc:    PE_PA = G4Material::GetMaterial("PE_PA");
./Simulation/DetSimV2/CentralDetector/src/XJfixtureConstruction.cc:     PE_PA,

Delving in shows the excess Opticks photons absorbed by "PE_PA" and with suspicions sign flip last boundaries.

Looks like an issue with XJfixtureConstruction

Investigate XJfixtureConstruction Solid

X4MeshTest_XJfixtureConstruction_png

cd ~/opticks/extg4 ; GEOM=XJfixtureConstruction ./X4MeshTest.sh: view of Geant4 polygonized solidXJfixture
holy geometry : celtic cross (tubs + FOUR box) atop an altar (TWO box): union of annular tubs and SIX boxes

cxs_custom_XJfixtureConstructionXY

COMP=custom_XY ./cxs_solidXJfixture.sh: cross section view of standalone solidXJfixture
All interior red intersects are spurious : not easy to fix.: probable cause of history discrepancy

So many spurious intersects : Easier to start over with fewer constituent solids

XJfixtureConstruction CSG Model

   Unexpected spurious vertical at 35          In Y box is at 50 +- 15
   (only appears with 4 box, not 2)                35    50     65
                                                    :
                                                    :     :     :                            altar frame      fixture frame
                                                    :
     -------------+                             +---+---+-+-----+        - - - - - - - - - - 18.5+13  = 31.5          6.5
     |            |                             |   :   |      13/2=6.5
     +            +                             +   :   + :     :         - - - - - - - - -  18.5+6.5 = 25            0.0
     |            |                             |   :   |       :
     +------------+----------------+-----25-----+---20--+-+-----+         - - - - - - -       8.5+10 =  18.5         -6.5
     |                                                          |
     +    up2                      +                            +       - - - - - - - - - -   8.5+5  = 13.5         -11.5
     |                                                          |
     +---------+^^^^^^^^^^^^^^^^^^^+^^^^^^^^^^^^^^^^^^+---------+       - - - - - - - - - -             8.5         -16.5
               |                                      |
               |                                     17/2=8.5
               +  up1                                 +                - - - - - - - - - -              0.0         -25.0
               |                                      |
               |                                      |
               +-------------------+-------40---------+            - - - - - - - - - - - - -           -8.5         -33.5

                                   |            |    :   |
        Z                          0            25   35  45
        |                                       |    :   |
        |                                       |    :   |
        +-- Y                                   |    :   outer tubs
       /                                        |    :
      X                                         |    spurious vertical from box edge (why? it is within the tubs ring)
                                                inner tubs

Due to construction technique and Geant4 limitations : simple anti-coincidence tricks cannot be used

Investigate XJfixtureConstruction Positions

`RTP (Radial-Theta-Phi) Tangential Frames`

RTP Tangential frame : target geometry without using instance frame: => can be used for any geometry, eg non-instanced XJfixtureConstruction => consistent control of viewpoints or genstep grids for all placements

Render 64 solidXJfixture from same viewpoint in RTP frame of each

cd ~/opticks/CSGOptiX
for n in $(seq 0 63) ; do
  EYE=2,-1,-1 LOOK=0,0,0 UP=1,0,0 MOI=solidXJfixture:$n:-3 ./cxr_view.sh
done

MOI identifier picks volume using solid names and repeat indices
Instanced Geometry	Global Geometry
`NNVT:0:1000`	`solidXJfixture:10` (default cartesian XYZ)
`Hama:0:2000`	`solidXJfixture:10:-3` (RTP tangential frame)

Components of MOI Mname:mOrdinal:InstanceIndex:

Mname: selects first volume with solid name starting with Mname
mOrdinal: usually 0 with instanced geometry, or selects between global repeats
InstanceIndex: selects instanced repeat or with global geometry selects between frame types

image_grid_cxr_solidXJfixture:XX:-3

Observations of the 64 solidXJfixture renders

First 8 (0-7) very different to the 56 others (8-63)
- 0-7 are at smaller radius than 8-63
- 0,1 positioned mid-chimney and overlap each other
- 2-7 located around sTarget radius, are obscured by uni_acrylic1
36,40,46 close to uni_acrylic1 of fastener
37,41,44,47,50,52,55,58,61 overlapped or cut by uni_acrylic1

probably inconsistency between the set of pos files ?

LSExpDetectorConstruction::setupCD_Sticks Implemented poorly

12 large copy/pasted blocks of code with small changes in each

obfuscates meaning by burying it in boilerplate (had to create table to understand it)
coding style that invites bugs, now and in future

Clarity requires higher level code (simply use static method to handle the repetitive parts)

-> less code, clearer code, understandable code

LSExpDetectorConstruction::setupCD_Sticks

member	class	pos_file	radius	#pos
m_strut_name	StrutAcrylicConstruction	m_strut_pos_file	strut_r	370
m_strut2_name	StrutBar2AcrylicConstruction	m_strut2_pos_file	strut2_r	220
m_strutballhead_name	StrutBallheadAcrylicConstruction	m_fastener_pos_file	strutballhead_r	590
m_fastener_name	FastenerAcrylicConstruction	m_fastener_pos_file	fastener_r
m_upper_name	UpperAcrylicConstruction	m_fastener_pos_file	upper_r
m_addition_name	AdditionAcrylicConstruction	m_fastener_pos_file	addition_r
m_xjanchor_name	XJanchorConstruction	m_xjanchor_pos_file	xjanchor_r	56
m_xjfixture_name	XJfixtureConstruction	m_xjanchor_pos_file	xjfixture_r	56
m_sjclsanchor_name	SJCLSanchorConstruction	m_sjclsanchor_pos_file	sjclsanchor_r	2
m_sjfixture_name	SJFixtureConstruction	m_sjfixture_pos_file	sjfixture_r	36
m_sjreceiver_name	SJReceiverConstruction	m_sjreceiver_pos_file	sjreceiver_r	8
m_sjreceiver_fastener_name	XJfixtureConstruction	m_sjreceiver_pos_file	sjreceiver_fastener_r	8

XJfixtureConstruction appears twice with different radius => two position files => 8+56 = 64 solidXJfixture
Position files loaded from $JUNOTOP/offline/Simulation/DetSimV2/DetSimOptions/data/Strut_Acrylic.csv,...

solidXJfixture64radii.png

cxr_view_solidXJfixture:55:-3_cam_1_eye_8,-4,-4_zoom_1_tmin_0.1

EYE=8,-4,-4 LOOK=0,0,0 MOI=solidXJfixture:55:-3 ./cxr_view.sh: view directly at the fixture, but not visible as uni_acrylic1 in front

cxr_view_solidXJfixture:55:-3_cam_1_eye_1,-0.5,-0.5_zoom_1_tmin_0.1

EYE=1,-0.5,-0.5 LOOK=0,0,0 MOI=solidXJfixture:55:-3 ./cxr_view.sh: closer again, at same angle of view : fixture now visible, coincidence speckle between spherically curved uni_acrylic1 base and sAcrylic

cxs_XJfixtureConstructionPR_55

COMP=PR_55 ./cxs_solidXJfixture.sh: PR cross section view of solidXJfixture:55:-3 from the +ve T axis : solidXJfixture is bumping into uni1

solidXJfixture:0:-3_cam_1_t0_cxr_view___eye_32,-48,0__look_32,0,0__zoom_1__tmin_48_solidXJfixture:0:-3.jpg

EYE=32,-48,0 LOOK=32,0,0 CAM=1 ./cxr_view.sh: wider view using CAM=1 parallel projection : fixture visible at base of chimney

solidXJfixture:0:-3_cam_1_t0_cxr_view___eye_32,-48,0__look_32,0,0__zoom_0.25__tmin_48_solidXJfixture:0:-3.jpg

EYE=32,-48,0 LOOK=32,0,0 CAM=1 ZOOM=0.25 ./cxr_view.sh: even wider view from same position, but with ZOOM=0.25 increasing the field-of-view : TT now visible

cxs_XJfixtureConstructionTR_0

COMP=TR_0 ./cxs_solidXJfixture.sh: TR cross section view of solidXJfixture:0:-3 which is mid-chimney showing overlap with solidSJReceiver

cxs_XJfixtureConstructionTP_0_Rshift

COMP=TP_0 ./cxs_solidXJfixture.sh: Same TP "tangential" cross section view of solidXJfixture:0:-3 but lower in radial direction showing crossed "celtic-crosses"

XJfixtureConstructionTR_2

COMP=TR_2 ./cxs_solidXJfixture.sh: TR cross section view of solidXJfixture:2:-3 showing positioning at sTarget radius (first 0-7 are at this lower radius)

XJfixtureConstructionTR_52

COMP=TR_52 ./cxs_solidXJfixture.sh: TR cross section view of solidXJfixture:52:-3 showing positioning at higher radius inside uni_acrylic1

solidXJfixture:41:-3_cam_1_t0_cxr_view___eye_2,-1,-1__zoom_1__tmin_0.1_solidXJfixture:41:-3.jpg

EYE=2,-1,-1 MOI=solidXJfixture:41:-3 ./cxr_view.sh: cxr render of solidXJfixture:41:-3 shows it poking out of uni_acrylic1

XJfixtureConstructionPR_41

COMP=PR_41 ./cxs_solidXJfixture.sh: PR cross section view of solidXJfixture:41:-3 showing uni_acrylic1 cutting across the fixture

XJfixtureConstructionCloserPR_41

COMP=PR_41 ./cxs_solidXJfixture.sh: closer PR cross section view of solidXJfixture:41:-3 showing uni_acrylic1 cutting across the fixture

Investigate XJanchorConstruction Solid

X4IntersectSolidTest_XJanchorConstruction_YZ

Spurious Geant4 Intersects Visible on line between cone and cylinder: not fixed by expanding cylinder upwards, goes away when dont subtract acrylic sphere

Investigate SJReceiverConstruction Solid

X4MeshTest_SJReceiverConstruction_png

GEOM=SJReceiverConstruction EYE=1,1,0.2 ZOOM=2.5 ./X4MeshTest.sh: union of two unconnected parts : inner cone face is spherically curved

Geant4 docs : component solids must not be disjoint

X4IntersectSolidTest_SJReceiverConstruction_XZ_png

cd ~/opticks/extg4 ; GEOM=SJReceiverConstruction_XZ ./xxs.sh

Geant4_docs_Boolean_Solids_png

https://geant4.web.cern.ch/sites/default/files/geant4/collaboration/working_groups/geometry/training/D2-Basics.pdf

Geant4_docs_Boolean_Solids_html

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Detector/Geometry/geomSolids.html?highlight=boolean#constructed-solid-geometry-csg-solids

Stick Geometry Issue Overview

class	solid	lv	mat	G4VSolid constituents
StrutAcrylicConstruction	sStrut	lSteel	StrutSteel	Tubs
StrutBar2AcrylicConstruction	sStrut	lSteel2	StrutSteel	Tubs
StrutBallheadAcrylicConstruction	sStrutBallhead	lSteel	Steel	Orb
FastenerAcrylicConstruction	uni1	lFasteners	Steel	Tubs, Union
UpperAcrylicConstruction (1)	base_steel	lUpper	Steel	Polycone
AdditionAcrylicConstruction (2)	uni_acrylic3/1	lAddition	Acrylic	Polycone, Sphere, Sub.
XJanchorConstruction (3)	solidXJanchor	lXJanchor	Acrylic	Tubs, Cons, Sphere, Sub., Union
XJfixtureConstruction (4)(5)	solidXJfixture	lXJfixture	PE_PA	Tubs, Box, Union
SJCLSanchorConstruction	solidSJCLSanchor	lSJCLSanchor	Acrylic	Box, Cons, Sphere, Sub., Union
SJFixtureConstruction	solidSJFixture	lSJFixture	Acrylic	Box, Cons, Sphere, Sub., Union
SJReceiverConstruction (6)	solidSJReceiver	lSJReceiver	Acrylic	Tubs, Cons, Box, Sphere, Sub., Union

DONE : fixed Opticks Polycone with multiple Rmin translation
needs --additionacrylic-simplify-csg to fix daughter volume cavity subtraction (TODO: make default)
spurious Geant4 intersects observed (subtracting big sphere implicated, TODO: try without)
overlaps uni_acylic1, uni1, self (WIP: Yuxiang assigned to update positions)
spurious coincidences (TODO: try simpler CSG modelling)
disjoint union + overlap with solidXJfixture

Overview

JUNO Simulation -> Opticks -> GPU
Benefits:	drastically faster GPU sim. in production improved CPU sim. speed + correctness
Status:	running out of geometry to have issues issues getting smaller, are converging

extremely sensitive to source geometry issues
- DONE: many fixes to PMTs, Fastener
- => improves CPU sim, enables GPU sim

Many Thanks to : Yuxiang Hu, Wang Zike

More pioneer users very welcome

https://bitbucket.org/simoncblyth/opticks	code repository
https://simoncblyth.bitbucket.io	presentations and videos

Extra Slides : G4Cerenkov Improvement, Cerenkov on GPU

Extra Slides : G4Cerenkov Improvement, CK on GPU

Frank-Tamm Formula : Cerenkov Photon Yield /mm at BetaInverse
- G4Cerenkov/G4Cerenkov_modified makes poor approximation => -ve photon yield close to RINDEX peak
Numerical Integration Imprecision close to RINDEX peak
Numerical Integration Error (number of CK photons)
Improved Precision "S2" Integration
Monte Carlo Sampling : "Rejection" vs "Inverse Transform"
Cerenkov ICDF for various BetaInverse
Chi2 Comparisons of the two sampling methods
- New Cerenkov ICDF sampling : faster, avoids double->float sensitivity

Frank-Tamm Formula : Cerenkov Photon Yield /mm at BetaInverse

                                                BetaInverse^2
          N_photon/369.81  =    Integral ( 1 - -----------------  )  dE         where   BetaInverse < ri[E]
                                                   ri(E)^2

                            =   Integral [ 1 ] dE -  BetaInverse^2 * Integral[ 1./(ri[E]*ri[E]) ] dE

G4Cerenkov::BuildThePhysicsTable -> CerenkovAngleIntegrals (misnomer)

Integral[ 1/ri^2 ] dE cumulative trapezoidal approx. integral over RINDEX[E] of material

Problems with G4Cerenkov::GetAverageNumberOfPhotons integration

assumes monotonic RINDEX : only one permitted energy range
putting together the split integral leads to -ve NumPhotons when close to RINDEX peak
G4PhysicsVector::GetValue applies linear interpolation to cumulative integral of 1/ri^2 <-- POOR APPROX

636   G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();
637
638   G4double dp, ge;
642   if (nMax < BetaInverse)        // ... no photons
649   else if (nMin > BetaInverse) {
650      dp = Pmax - Pmin;
651      ge = CAImax;
660   } else {
661      Pmin = Rindex->GetEnergy(BetaInverse);
662      dp = Pmax - Pmin;
667      G4double CAImin = CerenkovAngleIntegrals->GetValue(Pmin, isOutRange);
668      ge = CAImax - CAImin;
674   }
677   G4double NumPhotons = Rfact * charge/eplus * charge/eplus * (dp - ge * BetaInverse*BetaInverse);

scan_GetAverageNumberOfPhotons_plot_1.7000_1.8000

scan_GetAverageNumberOfPhotons_difference_plot

Alternative "s2" integral approach : more precise, simpler, faster


                 BetaInverse*BetaInverse
Integral [ 1. -  ----------------------- ] (for BetaInverse < RINDEX)
                 RINDEX * RINDEX

Integral [ 1. - cos^2 theta ]

Integral [ sin^2 theta ]

Integral [ s2 ]             ( s2 > 0 )

Do not split the integral, do "s2" integral in one pass. Advantages:

avoids one level of linear approximation, ca
cannot give -ve values
simple one pass code, no separate find_crossings
also s2 is faster other than when numPhotons is zero
- AVOID by not calling when BetaInverse > RINDEX_max

Monte Carlo "Rejection Sampling" VS "Inverse Transform Sampling"

Sampling : generate (x0,x1,x2,... ) that follow distribution f(x)

Rejection Sampling (used by current G4Cerenkov)

throw 2 random numbers (x,y) in (domain, value) range
reject x when y > f(x)
Pros : simple, no integration or lookup tables
Cons :
- poor efficiency for some f(x) : many candidates rejected
- observed float/double sensitivity : problem on GPU

Inverse Transform (ICDF) Sampling

compute CDF Cumulative Distribution Function, CDF(x) : integrals of f(x) across domain
normalize : CDF(xmin) = 0, CDF(xmax) = 1
domain x and CDF(x) monotonic -> invert (swap y<->x)
equate random u(0->1) to y enables to lookup sample x
Cons : requires computation of large ICDF lookup tables
Pros :
- simple sampling : one random lookup (eg GPU texture lookup)
- ICDF lookups tables can be created once for a geometry (often depending only on constant properties of materials)

Cerenkov radiation is electromagnetic equivalent of sonic boom in air or bow waves in water

test_makeICDF_SplitBin_QCKTest_s2cn_plot.png

Chi2 comparison of Rejection and ICDF Lookup Cerenkov energy samples

Created with:

cd ~/opticks/qudarap
QCerenkovIntegralTest   # create ICDF lookups for many BetaInverse
QCKTest                 # load ICDF, create lookup+rejection samples
ipython -i tests/QCKTest.py   # chi2 comparisons on energy distribs

Plots show:

histogram of energy distrib from ICDF lookup sample
histogram of energy distrib from rejection sample
"s2(E)" sin^2(cerenkov_cone_angle)[E] : distribution the samples seek to follow
ICDF : cumulative integral curve
relative chi2 contribution for each energy bin

Note: ICDF becomes flat in disallowed energy regions

QCKTest_1.5000_en_plot.png

QCKTest_1.6000_en_plot.png