Opticks + JUNO : PMT Geometry + Optical Model Progress

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

Many Opticks + JUNO Geometry Bugs fixed : NO KNOWN ISSUES REMAIN
Layr.h Multi-Layer TMM CPU/GPU header-only impl : VALIDATED
junoPMTOpticalModel standalone tests reveal issues, all from FastSim use:
- over complex PMT geometry and Geant4 step histories
- LIVE BUGS: polarization and propagation time within PMT
PLAN FIX : replace junoPMTOpticalModel FastSim with CustomG4OpBoundaryProcess:
- can fix above issues using a more standard approach
- simpler geometry, simpler simulation history, less code
- facilitates bringing PMT Optical Model to GPU within Opticks
OTHER WORK : more direct Opticks geometry translation (AIM: massive code reduction)
- substantial progress U4Material/U4Surface/U4Solid/stree/snode/snd

Simon C Blyth, IHEP, CAS — 6 Feb 2023

JUNO + Opticks : Timeline, see Refs for more detail

2022 Aug-Sep : Fixed geometry issues:

JUNO : PMT, PMTMask overlaps
Opticks : ellipsoid apex hole, thin cylinder, cone apex

2022 Oct-Nov : MultiFilmSimSvc -> GPU/CPU Layr.h

17 November 2022, JOC Meeting

Show PMT geometry issues and fixes already merged
Describe MultiFilmSimSvc -> GPU/CPU Layr.h
opticks_20221117_mask_debug_and_tmm.html

2022 Dec : Standalone few-PMT junoPMTOpticalModel test:

overcomplex : PMT geometry and Geant4 histories (inappropriate FastSim)
LIVE BUGS : incorrect PMT reflect/refract polarization, incorrect propagation speed in PMT

20 December 2022, AFG Sim Meeting

Demonstrated CustomG4OpBoundaryProcess can fix issues (using standalone few-PMT tests)
my impression: General agreement to proceed with fix, Guofu convinced
opticks_20221220_junoPMTOpticalModel_FastSim_issues_and_CustomG4OpBoundaryProcess_fix.html

2023 Jan : Start adapting above for use within junosw "monolith" by decoupling PMT data access

low dependency PMT data access : PMTSimParamData (MERGE STALLED: Jan 17->Feb 3)
OTHER WORK : more direct geometry translation : U4Material/U4Surface/U4Solid/stree/snode/snd

Details in above presentations + notes : https://bitbucket.org/simoncblyth/opticks/src/master/notes/progress.rst

JUNO + Opticks : Summary of Geometry Bug Fixes

OPTICKS PRIM ISSUES

nmskSolidMask : ellipsoid hole at "apex" issue

few per million "apex" rays find hole in ellipsoid
manifests as spurious "waist" intersects
FIXED : add zcut safety margin to CSG_ZSPHERE

nmskSolidMaskTail : very thin cylinder "lip" issues

requires uncoincidence shift of subtractee
BUT : even after uncoincide some spurious remain !
BUG 1 : mis-translated as CSG_DISC, not CSG_CYLINDER
BUG 2 : intersect precision loss for near-axial rays
FIXED : reimplement CSG_CYLINDER with less FLOPS

nmskSolidMaskVirtual : cone precision + near-apex issues

spurious sprinkle remains after uncoincide
BUG 1 : CSG_CONE precision loss in one quadratic root
BUG 2 : rays close to infinite cone apex yield spurious
FIXED : CSG_CONE robust roots + fix bad assumption

https://simoncblyth.bitbucket.io/env/presentation/opticks_20221117_mask_debug_and_tmm.html

`nmsk_nnvt_solids_STUVWXY_nnvt_virtualMask_Mask_overlap.png`

NNVT : ONE overlap issue visible, MaskTail impinges MaskVirtual

(using ct.sh : Opticks CSG on CPU)

`hmsk_hama_solids_STUVWXY_BodySolid_x_MaskTail.png`

HAMA : ONE overlap issue, BodySolid impinges MaskTail

(using ct.sh : Opticks CSG on CPU)

`hmsk_hama_solids_STUVWXY_zoom_impinge.png`

HAMA : BodySolid impinges MaskTail (mct.sh)

`hmsk_hama_solids_STUVWXY_after_bug_33_fix.png`

HAMA : after bug 33 fix

Layr.h : simple GPU/CPU header-only re-impl. of MultiFilmSimSvc

template<typename T> struct Layr
{
    T  d, pad   ; // d:thickness : zero means incoherent
#ifdef WITH_THRUST
    thrust::complex  n, st, ct, rs, rp, ts, tp ;
#else
    std::complex     n, st, ct, rs, rp, ts, tp ;
#endif
    Matx S, P ;
};

template <typename T, int N> struct Stack
{
    Layr<T> ll[N] ; // stack of N layers   
    Layr<T> comp ;  // composite for the N layers  
    ART_<T>  art ;  // results eg A,R,T

    Stack(T wl, T minus_cos_theta, const StackSpec<T>& ss);
};

https://github.com/simoncblyth/j/blob/main/Layr/Layr.h

simple impl : compiles with nvcc/gcc for GPU/CPU
using thrust::complex/std::complex
YES : struct, arrays, templates, enum, constexpr
NO : std::vector std::map std::string std::cout virtual meth.

Layr.h reproduces MultiFilmSimSvc ART results
- works on both CPU and GPU

ALT: multi-GB "ART" textures (Yuxiang Hu) TODO: compare

LayrTest.{sh,cc,h,cu,py} : AOI scan with float/double on CPU/GPU ...

 epsilon:Layr blyth$ ./LayrTest.sh ana
 CFLayrTest
  a :            EGet : scan__EGet__cpu_thr_double
  b :            EGet : scan__EGet__cpu_thr_float
  c :            EGet : scan__EGet__gpu_thr_double
  d :            EGet : scan__EGet__gpu_thr_float
  ...
  m :          R12860 : scan__R12860__cpu_pom_double
  n :          R12860 : scan__R12860__cpu_thr_double
  o :          R12860 : scan__R12860__cpu_thr_float
  p :          R12860 : scan__R12860__gpu_thr_double
  q :          R12860 : scan__R12860__gpu_thr_float

 In [1]: CF(m,q,0.05)
 Out[1]:
 CF(m,q,0.05) : scan__R12860__cpu_pom_double vs scan__R12860__gpu_thr_float
        lls :   0.000442 :   0.000442 :  -0.000414
      comps :   0.000341 :   0.000341 :  -6.17e-05
       arts :    6.2e-05 :    6.2e-05 :   -6.2e-05

 pmtcat:R12860 tt:5 lt:q : j/Layr/LayrTest scan__R12860__gpu_thr_float ni 900 wl 440
 +------------------------------+----------+----------+----------+----------+----------+
 |        R12860 arts\comps 0.05|     m:cpd|     n:ctd|     o:ctf|     p:gtd|     q:gtf|
 +==============================+==========+==========+==========+==========+==========+
 |                         m:cpd|         0| 0.0003325|  0.000301| 0.0003325| 0.0003407|
 +------------------------------+----------+----------+----------+----------+----------+   max difference of all param between scans 
 |                         n:ctd| 6.064e-05|         0| 4.829e-05| 7.445e-14| 4.829e-05|
 +------------------------------+----------+----------+----------+----------+----------+
 |                         o:ctf| 5.454e-05| 6.101e-06|         0| 4.829e-05| 3.977e-05|
 +------------------------------+----------+----------+----------+----------+----------+
 |                         p:gtd| 6.064e-05| 1.321e-14| 6.101e-06|         0| 4.829e-05|
 +------------------------------+----------+----------+----------+----------+----------+
 |                         q:gtf| 6.199e-05| 1.523e-06| 7.451e-06| 1.523e-06|         0|
 +------------------------------+----------+----------+----------+----------+----------+

 In [1]: ARTPlot(m, 0.05)

LayrTest_4_ARTQspx_R12860_4layer.png

S/P-Pol : perp./within plane

S/P Polarization : has huge effect on A,R,T in both directions : AOI > 90 is reversed stack

LayrTest_4_AQspx_R12860_4layer.png

4 : A

Absorption : large S/P-polarization difference in both directions

LayrTest_4_Rspx_R12860_4layer.png

4 : R

Only reflection at glancing incidence (90 deg) : large S/P-pol. difference for Vacuum->Pyrex

LayrTest_4_Tspx_R12860_4layer.png

4 : T

Pyrex->Vacuum : Above critical angle : no transmission, only reflect or absorb

Standalone few PMT test of Geometry and Optical Model

j:PMTFastSim Standalone PMT geometry and Optical Model

https://github.com/simoncblyth/j/
j:PMTFastSim/junoPMTOpticalModel.cc
j:PMTFastSim/HamamatsuR12860PMTManager.cc
N=0/1 natural geometry switch

Optical only Geant4 simulation, with full step point recording

opticks:u4/tests/U4SimulateTest.sh

N=0/1 ./U4SimulateTest.sh

Geant4 Simtrace intersection : 2D geometry plotting

opticks:u4/tests/U4SimtraceTest.sh

N=0/1 ./U4SimtraceTest.sh

junoPMTOpticalModel Issues : All from using FastSim

simple 2 volume (Pyrex+Vacuum) geometry
standard G4OpBoundaryProcess reflect/refract/absorb/TIR
much less code : reuse of standard Geant4
easier Opticks on GPU : simpler, more standard
- no near-degenerate 1e-3 mm boundaries
- allows same geometry on CPU and GPU
- avoid history rewriting for comparisons

BUT : need to change very nasty Geant4 impl:

G4OpBoundaryProcess::PostStepDoIt: ugly monolith, rats nest of nested ifs

Unpleasant to get familiar with
Change palatable, using : CustomART.h
Keep change separated, detailed validation vital

Four volume PMT due to FastSim limitation
- Two volumes (Pyrex + Vacuum) would be natural

Fake volumes yield many fake intersects
- Pyrex/Pyrex + Vacuum/Vacuum same material fakes
- complicates Geant4 step point history
- complicates comparison with Opticks

Non-FastSim in PMT propagates at Pyrex (not Vacuum) speed
- FastSim->SlowSim transition misses speed fixup

reflected + refracted polarization doesnt follow Geant4
- should follow G4OpBoundaryProcess, depending on S/P

All issues stem from using FastSim (designed for regions)
- not good fit for 1 boundary (Pyrex/ARC/PHC/Vacuum)

Polarization on FresnelReflection/FresnelRefraction/TotalInternalReflection

_E2_t(make_float2( 2.f*n1c1/(n1c1+n2c2) ,
                   2.f*n1c1/(n2c1+n1c2) )),// (ts,tp)
_E2_r(make_float2( _E2_t.x - 1.f       ,
                   n2*_E2_t.y/n1 - 1.f  )),// (rs,rp)
A_transverse(cross(p.mom, orient*normal)),
E1_perp(dot(p.pol, A_transverse)),
E1(make_float2(E1_perp,sqrtf(1.f-E1_perp*E1_perp)))),
E2_t(_E2_t*E1),RR(normalize(E2_r)),
E2_r(_E2_r*E1),TT(normalize(E2_t)),
{
p.mom = reflect
      ?
         p.mom + 2.0f*c1*orient*normal
      :
         eta*(p.mom) + (eta*c1 - c2)*orient*normal
      ;
A_parallel = normalize(cross(p.mom, A_transverse));

p.pol =
        ( reflect ?
             ( tir ?
                   -p.pol + 2.f*EdotN*orient*normal
             :
                   RR.x*A_transverse + RR.y*A_parallel
             )
          :
             TT.x*A_transverse + TT.y*A_parallel
        )
      ;
}  // IMPL BASED ON + VALIDATED AGAINST GEANT4

Reflect/Refract : INCORRECT POLARIZATION

 void junoPMTOpticalModel::Reflect()
 {
     dir -= 2.*(dir*norm)*norm;
     pol -= 2.*(pol*norm)*norm;
 }
 void junoPMTOpticalModel::Refract()
 {
     dir = (real(_cos_theta4) -
            _cos_theta1*_n1/_n4)*norm + (_n1/_n4)*dir;
     pol = (pol-(pol*dir)*dir).unit();
 }

G4OpBoundaryProcess polarization from:

S/P Fresnel Coeff
S/P Directions : out/in plane of incidence

Opticks : sysrap/sboundary.h qudarap/qsim.h follows G4, ===>

opticks/src/master/sysrap/sboundary.h

Derive Fresnel from Maxwell Boundary Conditions

Fresnel Equations, Alexander I. Lvovsky: Encylopedia of Optical Engineering (fresnel-eoe.pdf)

Compare Reflected Polarization Impls for Brewster Angle Incidence

opticks/src/master/sysrap/tests/sboundary_test.sh : build, run, plot

sboundary_test/figs/sboundary_test/pvcap/AOI_BREWSTER_linear_polarization_by_reflection.png

sboundary_test/figs/sboundary_test/pvcap/AOI_BREWSTER_reflect_alt_pol.png

G4OpBoundaryProcess/qsim.h/sboundary.h : Only S-polarized survives

junoPMTOpticalModel::Reflect : very different

incident from left (-X), surface normal vertically upwards (+Z), intersection point in middle,
Colored lines represent polarization directions of 128 photons before and after Reflection
Reflected using sboundary.h (validated against G4OpBoundaryProcess)
Compared with sboundary.h:alt_pol that duplicates junoPMTOpticalModel::Reflect

Brewster (or polarizing) incident angle th1 : tan(th1) = n2/n1 ; th1 + th2 = pi/2

Compare Refracted Polarization Impls for Brewster Angle Incidence

opticks/src/master/sysrap/tests/sboundary_test.sh : build, run, plot

$sboundary_test/figs/sboundary_test/pvcap/AOI_BREWSTER_brewster_refract.png$ $sboundary_test/figs/sboundary_test/pvcap/AOI_BREWSTER_brewster_refract_alt_pol.png$

G4OpBoundaryProcess/qsim.h/sboundary.h : partial pol

junoPMTOpticalModel::Refract : not partially pol

incident from left (-X), surface normal vertically upwards (+Z), intersection point in middle,
Colored lines represent polarization directions of 128 photons before and after Refraction
Refracted using sboundary.h (compared with alt_pol that duplicates junoPMTOpticalModel::Refract )

hamaLogicalPMT_natural.png

Natural 2-volume PMT (Pyrex + Vacuum) : NOT COMPATIBLE WITH FastSim

hamaLogicalPMT_fake.png

body hidden under inner1 inner2

body is FastSim envelope volume

CURRENT : Unnatural 4-volume PMT (Pyrex+Pyrex+Vacuum+Vacuum) => Half Fake

junoPMTOpticalModel FastSim : Fake Volumes -> Contorted ModelTrigger

  +---------------pmt-Pyrex----------------+
  |                                        |
  |                                        |
  |     +----body-Pyrex-(FSim-env)---+     |
  |     | +------------------------+ |     |
  |     | |                        | |     |
  |     | |                        | |     |
  |     | |        inner1-Vacuum   |-|     |
  |     | |                        |1e-3   |
  |     | |                        | |     |
  |     | +~~coincident~face~~~~~~~+ |     |
  |     | |                        | |     |
  |     | |                        | |     |
  |     | |        inner2-Vacuum   | |     |
  |     | |                        | |     |
  |     | |                        | |     |
  |     | +------------------------+ |     |
  |     +----------------------------+     |
  |                                        |
  |                                        |
  +----------------------------------------+

ModelTrigger : unnatural trig @ Vac/Vac bnd.
DoIt : shift to edge, for ART calc
only m_mirror_opsurf (inner2<->body) "active"

   80 G4bool junoPMTOpticalModel::ModelTrigger(const G4FastTrack &fastTrack)
   81 {   //  Contorted as "pre-trigger" on where next 
   83     if(fastTrack.GetPrimaryTrack()->GetVolume() == _inner2_phys){
   84         return false;
   85     }
   ...
   89     pos     = fastTrack.GetPrimaryTrackLocalPosition();
   90     dir     = fastTrack.GetPrimaryTrackLocalDirection();
   94
   95     if(fastTrack.GetPrimaryTrack()->GetVolume() == _inner1_phys){
   96         whereAmI = kInVacuum;
   97     }else{
   98         whereAmI = kInGlass;
   99     }
  100
  101     if(whereAmI == kInGlass){
  102         dist1 = _inner1_solid->DistanceToIn(pos, dir);
  103         dist2 = _inner2_solid->DistanceToIn(pos, dir);
  104
  105         if(dist1 == kInfinity){
  106             return false;
  107         }else if(dist1>dist2){
  108             return false;
  109         }else{
  110             return true;
  111         }
  112     }else{
  113         dist1 = _inner1_solid->DistanceToOut(pos, dir);
  114         dist2 = _inner2_solid->DistanceToIn(pos, dir);
  115
  116         if(dist2 == kInfinity){
  117             return true;
  118         }
  119     }
  120     return false;
  121 }

junoPMTOpticalModelSimple : Natural Geometry -> Simple ModelTrigger

  +---------------pmt-Pyrex----------------+
  |                                        |
  |                                        |
  |                                        |
  |       +~inner~Vacuum~(FSim~env)+       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       +                        +       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       +------------------------+       |
  |                                        |
  |                                        |
  |                                        |
  +----------------------------------------+

ModelTrigger : z>0 : no unnatural trig
override mirror_opsurf (inner<->pmt) @ z > 0
simpler DoIt : no geom. queries, no shifting

 G4bool junoPMTOpticalModelSimple::ModelTrigger(
         const G4FastTrack &fastTrack)
 {
      //  Simple + Cheap due to natural geometry 
     return fastTrack.GetPrimaryTrackLocalPosition().z() > 0. ;
 }

https://github.com/simoncblyth/j/blob/main/PMTFastSim/junoPMTOpticalModelSimple.cc

junoPMTOpticalModelSimple : Natural Geometry -> Simple DoIt

  +---------------pmt-Pyrex----------------+
  |                                        |
  |                                        |
  |                                        |
  |       +~inner~Vacuum~(FSim~env)+       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       +                        +       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       +------------------------+       |
  |                                        |
  |                                        |
  |                                        |
  +----------------------------------------+

ModelTrigger : z>0 : no unnatural trig
override mirror_opsurf (inner<->pmt) @ z > 0
simpler DoIt : no geom. queries, no shifting

 void junoPMTOpticalModelSimple::DoIt(const G4FastTrack& fastTrack,
       G4FastStep &fastStep)
 {
     G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
     wavelength_nm = twopi*hbarc/energy/nm ;

     position = fastTrack.GetPrimaryTrackLocalPosition();
     direction = fastTrack.GetPrimaryTrackLocalDirection();
     polarization = fastTrack.GetPrimaryTrackLocalPolarization();

     G4VSolid* envelope_solid = fastTrack.GetEnvelopeSolid();
     surface_normal = envelope_solid->SurfaceNormal(position);

     minus_cos_theta = direction*surface_normal ;
     whereAmI = minus_cos_theta < 0. ? kInGlass : kInVacuum ;

     StackSpec spec ;
     // ... skip : collect refractive indices, thickness into spec ...

     Stack stack(      wavelength_nm, minus_cos_theta, spec );
     Stack stackNormal(wavelength_nm, -1.            , spec );
     // ...
}

Simple : BUT NOT WORKABLE
Cannot stop: G4OpBoundaryProcess
- REAL REASON FOR THE FAKE VOLUMES
- Workaround limitations of Geant4 FastSim

Pivot to CustomG4OpBoundaryProcess For PMT Optical Model ?

  +---------------pmt-Pyrex----------------+
  |                                        |
  |                                        |
  |                                        |
  |       +~inner~Vacuum~~~~~~~~~~~+       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       !                        !       |
  |       +                        +       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       |                        |       |
  |       +------------------------+       |
  |                                        |
  |                                        |
  |                                        |
  +----------------------------------------+

OpSurfaceName[0] == '@'

local_z>0 : does MultiLayer ART calc
!(local_z>0) : standard mirror_opsurf

Cannot use natural geometry with FastSim
- SO: CustomG4OpBoundaryProcess

Custom Boundary Process : Advantages

standard Geant4 polarization, propagation
less code, simpler code
natural geometry, avoiding fakes
simpler Geant4 step history (no fakes)
same geometry on GPU+CPU, easier Opticks validation
half the geometry objects to model PMT (4->2)

Current FastSim POM	4 Solid, 4 LV, 4 PV
Custom Boundary POM	2 Solid, 2 LV, 2 PV

Disadvantages

maintain custom G4OpBoundaryProcess
updating Geant4 may require care (depending on changes)

Advantages far outweigh disadvantages

u4/CustomART.h : include into CustomG4OpBoundaryProcess.cc

inline char CustomART::maybe_doIt(
     const char* OpticalSurfaceName,
     const G4Track& aTrack, const G4Step& aStep )
{
    if( OpticalSurfaceName == nullptr ||
        OpticalSurfaceName[0] != '@') return 'N' ;

    const G4AffineTransform& transform =
      aTrack.GetTouchable()
      ->GetHistory()->GetTopTransform();

    G4ThreeVector localPoint =
      transform.TransformPoint(theGlobalPoint);

    if(localPoint.z() <= 0) return 'Z' ;

    return doIt(aTrack, aStep) ;
}

doIt TMM calc, runs only for:

OpticalSurfaceName[0]=='@' && zlocal > 0
-> Transmittance, Reflectivity, Efficiency
TMM calc using Layr.h multi-layer Stack

#include "IPMTAccessor.h"
#include "Layr.h"
#include "U4Touchable.h"

struct CustomART
{
    const IPMTAccessor* accessor ; // JPMT.h OR PMTAccessor.h 

    G4double& theTransmittance ;   // doIt sets these
    G4double& theReflectivity ;
    G4double& theEfficiency ;

    const G4ThreeVector& theGlobalPoint ;
    const G4ThreeVector& OldMomentum ;
    const G4ThreeVector& OldPolarization ;
    const G4ThreeVector& theRecoveredNormal ;
    const G4double& thePhotonMomentum ;

    CustomART(
        const IPMTAccessor* accessor,
        G4double& theTransmittance,
        G4double& theReflectivity,
        G4double& theEfficiency,
        const G4ThreeVector& theGlobalPoint,
        const G4ThreeVector& OldMomentum,
        const G4ThreeVector& OldPolarization,
        const G4ThreeVector& theRecoveredNormal,
        const G4double& thePhotonMomentum
    );
    char maybe_doIt(const char* OpticalSurfaceName,
              const G4Track& aTrack, const G4Step& aStep );
    char doIt(const G4Track& aTrack, const G4Step& aStep );
};

opticks/src/master/u4/CustomART.h

hamaLogicalPMTWrapLV_full_history.png

TO  BT  BT  BT  BT  SR  SR  BT  BR  BR  BT  SR  SR  SR  BT  BR  BT  SR  BT  SA    Lots of fakes 
00  01 [02] 03 [04] 05  06 [07] 08  09 [10] 11  12  13 [14] 15 [16] 17 [18] 19   (7/20 Fake)

hamaLogicalPMTWrapLV_natural.png

TO  BT  BT  SR  SR  BR  BR  SR  SR  SR  BR  SR  BR  SR  SA
00  01  02  03  04  05  06  07  08  09  10  11  12  13  14       Simpler history, no fakes 
00  01  03  05  06  08  09  11  12  13  15  17  19

Compare "big bouncer" position, time between N=0,1 (1)

u4t ; ./U4SimulateTest.sh cf

In [12]: np.c_[ar[:,0],np.arange(len(ar)),br[:,0]]
Out[12]:
array([[-113.   ,    0.   ,  200.   ,    0.   ,    0.   , -113.   ,    0.   ,  200.   ,    0.   ],
       [-113.   ,    0.   ,  170.163,    0.137,    1.   , -113.   ,    0.   ,  170.163,    0.137],
       [-112.83 ,    0.   ,  164.918,    0.164,    2.   , -112.83 ,    0.   ,  164.917,    0.164],
       [-112.83 ,    0.   ,  164.917,    0.164,    3.   , -156.577,    0.   , -148.846,    1.22 ],
       [-135.824,    0.   ,    0.   ,    1.012,    4.   ,  -95.   ,    0.   , -104.211,    1.474],
       [-156.577,    0.   , -148.846,    1.778,    5.   , -248.807,    0.   ,    7.28 ,    2.108],
       [ -95.   ,    0.   , -104.211,    2.166,    6.   ,   53.206,    0.   ,  180.727,    3.269],
       [-238.764,    0.   ,   -0.   ,    3.071,    7.   ,  245.605,    0.   ,  -35.443,    4.235],
       [-248.807,    0.   ,    7.28 ,    3.112,    8.   ,   95.   ,    0.   ,  -99.428,    4.781],
       [  53.205,    0.   ,  180.727,    4.274,    9.   ,  177.724,    0.   , -134.574,    5.08 ],
       [ 214.06 ,    0.   ,    0.   ,    5.507,   10.   ,  141.059,    0.   ,  152.451,    6.046],
       [ 245.605,    0.   ,  -35.443,    5.749,   11.   , -239.66 ,    0.   ,  -55.195,    7.492],
       [  95.   ,    0.   ,  -99.428,    6.583,   12.   ,  237.91 ,    0.   ,   54.597,    9.127],
       [ 177.724,    0.   , -134.574,    7.041,   13.   ,   50.   ,    0.   ,  -63.74 ,    9.867],
       [ 160.533,    0.   ,    0.   ,    7.732,   14.   ,   58.352,    0.   ,  -69.   ,    9.9  ],
       [ 141.059,    0.   ,  152.451,    8.245,   15.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [-138.46 ,    0.   ,    0.   ,    9.867,   16.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [-239.66 ,    0.   ,  -55.195,   10.455,   17.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [   0.427,    0.   ,    0.   ,   11.71 ,   18.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [ 237.91 ,    0.   ,   54.596,   12.523,   19.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [   0.   ,    0.   ,    0.   ,    0.   ,   20.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [   0.   ,    0.   ,    0.   ,    0.   ,   21.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [   0.   ,    0.   ,    0.   ,    0.   ,   22.   ,    0.   ,    0.   ,    0.   ,    0.   ],
       [   0.   ,    0.   ,    0.   ,    0.   ,   23.   ,    0.   ,    0.   ,    0.   ,    0.   ],

ar	a.record[PID]	N=0 current geom	degenerate and fake intersect points
bb	b.record[PID]	N=1 natural geom	less points, simpler history

Need point-to-point mapping to compare

Compare "big bouncer" position, time between N=0,1 (2)

 u4t ; ./U4SimulateTest.sh cf

 In [2]: b2a  ## point-to-point mapping to skip a fakes
 Out[2]: array([ 0,  1,  3,  5,  6,  8,  9, 11, 12, 13, 15, 17, 19])

 In [4]: abr = np.c_[ar[b2a,0],br[:len(b2a),0],ar[b2a,0]-br[:len(b2a),0]] ; abr
 Out[4]:
 array([[-113.   ,    0.   ,  200.   ,    0.   , -113.   ,    0.   ,  200.   ,    0.   ,    0.   ,    0.   ,    0.   ,    0.   ],
        [-113.   ,    0.   ,  170.163,    0.137, -113.   ,    0.   ,  170.163,    0.137,    0.   ,    0.   ,    0.   ,    0.   ],
        [-112.83 ,    0.   ,  164.917,    0.164, -112.83 ,    0.   ,  164.917,    0.164,    0.   ,    0.   ,    0.   ,   -0.   ],
        [-156.577,    0.   , -148.846,    1.778, -156.577,    0.   , -148.846,    1.22 ,   -0.   ,    0.   ,    0.   ,    0.558],
        [ -95.   ,    0.   , -104.211,    2.166,  -95.   ,    0.   , -104.211,    1.474,    0.   ,    0.   ,    0.   ,    0.692],
        [-248.807,    0.   ,    7.28 ,    3.112, -248.807,    0.   ,    7.28 ,    2.108,    0.   ,    0.   ,   -0.   ,    1.004],
        [  53.205,    0.   ,  180.727,    4.274,   53.206,    0.   ,  180.727,    3.269,   -0.   ,    0.   ,    0.   ,    1.004],
        [ 245.605,    0.   ,  -35.443,    5.749,  245.605,    0.   ,  -35.443,    4.235,    0.   ,    0.   ,    0.   ,    1.514],
        [  95.   ,    0.   ,  -99.428,    6.583,   95.   ,    0.   ,  -99.428,    4.781,    0.   ,    0.   ,    0.   ,    1.802],
        [ 177.724,    0.   , -134.574,    7.041,  177.724,    0.   , -134.574,    5.08 ,    0.   ,    0.   ,    0.   ,    1.96 ],
        [ 141.059,    0.   ,  152.451,    8.245,  141.059,    0.   ,  152.451,    6.046,   -0.   ,    0.   ,    0.   ,    2.199],
        [-239.66 ,    0.   ,  -55.195,   10.455, -239.66 ,    0.   ,  -55.195,    7.492,    0.   ,    0.   ,    0.   ,    2.963],
        [ 237.91 ,    0.   ,   54.596,   12.523,  237.91 ,    0.   ,   54.597,    9.127,    0.   ,    0.   ,   -0.   ,    3.397]], dtype=float32)

Positions match closely, some times are way off

Compare "big bouncer" position, time, dist, speed between N=0,1

rvtime_ = lambda r:np.diff(r[:,0,3])
rvstep_ = lambda r:np.diff(r[:,0,:3],axis=0 )
rvdist_ = lambda r:np.sqrt(np.sum(rvstep_(r)*rvstep_(r),axis=1))
rvspeed_ = lambda r:rvdist_(r)/rvtime_(r)

In [14]: np.c_[rvtime_(ar[b2a]), rvtime_(br[:len(b2a)]),
               rvdist_(ar[b2a]), rvdist_(br[:len(b2a)]),
              rvspeed_(ar[b2a]),rvspeed_(br[:len(b2a)])]

Out[14]:
array([[  0.137,   0.137,  29.837,  29.837, 218.038, 218.038],  ## Water
       [  0.027,   0.027,   5.249,   5.249, 196.216, 196.215],  ## Pyrex
       [  1.615,   1.057, 316.798, 316.798, 196.215, 299.792],  ## Vacuum
       [  0.388,   0.254,  76.053,  76.053, 196.215, 299.792],
       [  0.946,   0.634, 189.965, 189.965, 200.744, 299.792],  ## comb. of Vacuum and Pyrex speeds, split at Fake 
       [  1.162,   1.162, 348.275, 348.275, 299.792, 299.792],
       [  1.475,   0.965, 289.392, 289.392, 196.215, 299.792],
       [  0.834,   0.546, 163.634, 163.634, 196.215, 299.792],
       [  0.458,   0.3  ,  89.881,  89.881, 196.215, 299.792],
       [  1.204,   0.965, 289.357, 289.357, 240.315, 299.792],  ## comb. of Vacuum and Pyrex speeds, split at Fake
       [  2.21 ,   1.447, 433.663, 433.663, 196.215, 299.792],
       [  2.068,   1.635, 490.027, 490.027, 236.919, 299.792]], dtype=float32)

N=0	Pyrex speed within PMT Vacuum	FastSim->SlowSim transitions miss speed setup ?
N=1	Always Vacuum speed in Vacuum	All standard Geant4 with customized G4OpBoundaryProcess

hamaLogicalPMT_two_pmt_geom.png

hamaLogicalPMT_two_pmt_cf.png

Low dependency access to PMT data (26x faster)

CODING BEST PRACTICE
- less deps -> more useful code
test separate from the "monolith"
- development cycle < 1 s (vs minutes)
serialization with NP.hh,NPFold.h
- https://github.com/simoncblyth/np/
- 0-deps (hdr-only) -> PMT data anywhere

QE scan over 100 energy points, all PMTs

get_pmtid_qe( pmtid, en );

IPMTSimParamSvc	7.70 s
PMTSimParamData	0.30 s
IPMTSimParamSvc/PMTSimParamData	26.10

exactly same results, 26x faster

junosw/-/merge_requests/126 MERGED AFTER ~3 WEEKS

junosw/-/issues/66

junosw/-/tree/blyth-66-low-dependency-PMT-data-access

PMTSimParamSvc and PMTParamSvc: => too many deps for CustomG4OpBoundaryProcess

hdrs for collecting data are irrelevant to access
-> separate data struct from Svc
- header-only impl : PMTSimParamData PMTParamData
- _PMTSimParamData : serialization save/load
reference data struct from Svc
- minimal changes to Svc, API
- consistency guaranteed
- no change with diff src : file, DB, Frontier
speedup not the motivation, BUT 26x faster by:
- calc qe_scale ONCE, not for every photon
- large pmtid std::map -> std::vector

U4Tree.h/stree.h : Direct Geometry Translation

Geant4 -> U4Tree.h/stree.h -> CSG

Minimal intermediate stree.h model

stree.h : n-ary tree of snode volumes
scsg.hh : n-ary trees of snd constituent solids
(much less code)

U4Tree: G4VPhysicalVolume -> stree/snode
U4Solid: G4VSolid -> scsg/snd
U4Material: G4Material -> NPFold/NP
U4Surface: G4LogicalSurface -> NPFold/NP

stree/snode/scsg/snd -> CSGFoundry/CSGSolid/CSGPrim/CSGNode

Next Steps : Custom Boundary + Natural Geometry Branch

junosw : Natural PMT Geometry Branch

almost all standalone tested already, need to bring into junosw monolith
need switch : FastSim/standard-Sim, natural geometry, custom boundary

Primary	Support
CustomG4OpBoundaryProcess CustomART.h Layr.h	PMTSimParamData (merged in MR 126) PMTAccessor/IPMTAccessor (j/Layr) HamamatsuR12860PMTManager NNVTMCPPMTManager

Opticks : Custom boundary equivalent

PMTSimParamData PMTAccessor -> GPU
natural geometry => easier to validate (no fakes)
standard propagation + polarization => easier (already implemented qudarap/qsim.h )