Opticks + JUNO : More junoPMTOpticalModel Issues + Validation of Custom4 C4OpBoundaryProcess Fix

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

junoPMTOpticalModel issues:
- FORCED OVERLY COMPLEX: PMT geometry and Geant4 step histories
- BUGS: polarization and propagation time within PMT are wrong
- NEW: MORE BUGS: unphysical mid-vacuum reflect, refract, absorb, detect, "tunneling"
Custom Boundary POM (PMT Optical Model) Alternative, Custom4 mini-package
- C4OpBoundaryProcess, C4CustomART.h
- POM/PMT quadrant controls
Standalone PMT tests (Custom4 + j/PMTSim + Opticks)
- NEW: PMT=0,1 Chi2 comparision reveals ModelTrigger whereAmI bug
- Unphysical photon examples, mid-vacuum A,R,T,D,"Tunneling"
- ModelTrigger_Debug, mis-placed triggers, HAMA immunity
- finding and fixing junoPMTOpticalModel::ModelTrigger whereAmI bug
- PMT=0,1 Chi2 comparisons of photon histories
"Insitu" tests within junosw, with Opticks instrumentation
- input photons, photon history recorded into Opticks/SEvt
- PMT=0,1 Chi2 comparisons with input photons

Simon C Blyth, IHEP, CAS — 28 April 2023

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 : C4CustomART.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 propagate at Pyrex(not Vac.) speed
- FastSim->SlowSim transition misses speed fixup

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

ModelTrigger whereAmI BUG => Unphysical photons
- mid-vacuum A,R,T,D inside PMT
- "tunneling" through dynode/MCP

hamaLogicalPMT_fake.png

Biggest problem with FastSim based POM : forces using unnatural geometry

body hidden under inner1 inner2

body is FastSim envelope volume

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

hamaLogicalPMT_natural.png

Biggest advantage with Custom Boundary Process POM : enables natural geometry

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

FIX : Pivot to Custom Boundary Process For PMT Optical Model ?

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

OpSurfaceName[0] == '@'

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

Custom Boundary Process : Advantages

natural geometry, no fakes
standard Geant4 polarization, propagation, time
less code, simpler code
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 Custom4 C4OpBoundaryProcess
updating Geant4 may require care (depending on changes)

Advantages far outweigh disadvantages

C4OpBoundaryProcess::PostStepDoIt : 3-way (A,R,T) Customization

if( m_custom_status == 'Y' )
{
    G4double rand = G4UniformRand();
    if ( rand < theAbsorption )
    {
        DoAbsorption();     // A
    }
    else
    {
        DielectricDielectric();  // R or T
    }
}
else if (type == dielectric_metal)
{
      DielectricMetal();
}
...

Standard DoAbsorption DielectricDielectric reused

G4OpBoundaryProcess unless OpticalSurfaceNam[0] == '@'/'#'

if( OpticalSurfaceName0 == '@' || OpticalSurfaceName0 == '#' )
{
    if( m_custom_art->local_z(aTrack) < 0. ) // Lower hemi : Standard
    {
        m_custom_status = 'Z' ;
    }
    else if( OpticalSurfaceName0 == '@') // MultiFilm ART POM
    {
        m_custom_status = 'Y' ;
        m_custom_art->doIt(aTrack, aStep) ;

        type = dielectric_dielectric ;
        theModel = glisur ;
        theFinish = polished ;
    }
    else if( OpticalSurfaceName0 == '#' ) // Traditional POM
    {
        m_custom_status = '-' ;

        type = dielectric_metal ;
        theModel = glisur ;
        theReflectivity = 0. ;
        theTransmittance = 0. ;
        theEfficiency = 1. ;
    }
}

DielectricDielectric expects 2-way (theReflectivity,theTransmittance) => so C4CustomART.h rescales 3-way (A,R,T)

(theAbsorption,theReflectivity,theTransmittance) <= (A, R/(1-A), T/(1-A))

custom4/C4CustomART.h : include into C4OpBoundaryProcess.cc

#include "C4IPMTAccessor.h"
#include "C4MultiLayrStack.h"
#include "C4Touchable.h"

struct C4CustomART {
    const C4IPMTAccessor* accessor ; 
    G4double& theAbsorption ;   // doIt sets these
    G4double& theReflectivity ;
    G4double& theTransmittance ;
    G4double& theEfficiency ;

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

    C4CustomART(
        const C4IPMTAccessor* accessor,
        G4double& theAbsorption,
        G4double& theReflectivity,
        G4double& theTransmittance,
        G4double& theEfficiency,
        const G4ThreeVector& theGlobalPoint,
        const G4ThreeVector& OldMomentum,
        const G4ThreeVector& OldPolarization,
        const G4ThreeVector& theRecoveredNormal,
        const G4double& thePhotonMomentum
    );
    double local_z( const G4Track& aTrack );
    void doIt(const G4Track& aTrack, const G4Step& aStep ); 
};

https://github.com/simoncblyth/customgeant4/blob/main/C4CustomART.h

Custom4 : mini-package with Geant4 customizations

Avoid junosw-opticks circular dependency:
- by splitting off Custom4 mini-package
Custom4 depends only on Geant4
junosw depends on Custom4
opticks uses Custom4 for detailed testing

Added to : junoenv/junoenv-external-libs.sh

Auto-installed like other externals (eg Geant4), or:

cd $JUNOTOP/junoenv
bash junoenv libs get custom4
bash junoenv libs conf custom4
bash junoenv libs make custom4
bash junoenv libs install custom4

bash junoenv libs all custom4

https://github.com/simoncblyth/customgeant4/

PMT/POM quadrants : choice of PMT geometry and Optical Model impl.

4 POM*PMT quadrants

4 POM*PMT quadrants	Unnatural PMT Impl `--pmt-unnatural-geometry` m_UsePMTNaturalGeometry:false	Natural PMT Impl `--pmt-natural-geometry` m_UsePMTNaturalGeometry:true
Traditional POM (photons stop at photocathode, none enter PMT) PMT innards modelled with only (QE,CE) `--no-pmt-optical-model` m_UsePMTOpticalModel:false	OpticalSurfaceName without prefix	OpticalSurfaceName starting '#' C4OpBoundaryProcess active
MultiFilm POM (photons refract into PMT, complex rindex layers) `--pmt-optical-model` m_UsePMTOpticalModel:true	junoPMTOpticalModel FastSim in control (boundary process not run) OpticalSurfaceName without prefix	C4MultiLayrStack.h C4CustomART.h OpticalSurfaceName starting '@' C4OpBoundaryProcess Active DEFAULT QUADRANT

Unnatural PMT Impl

--pmt-unnatural-geometry m_UsePMTNaturalGeometry:false

Natural PMT Impl

--pmt-natural-geometry m_UsePMTNaturalGeometry:true

Traditional POM (photons stop at photocathode, none enter PMT)

PMT innards modelled with only (QE,CE)

--no-pmt-optical-model m_UsePMTOpticalModel:false

OpticalSurfaceName without prefix

OpticalSurfaceName starting '#'

C4OpBoundaryProcess active

MultiFilm POM (photons refract into PMT, complex rindex layers)

--pmt-optical-model m_UsePMTOpticalModel:true

junoPMTOpticalModel FastSim in control (boundary process not run)

OpticalSurfaceName without prefix

C4MultiLayrStack.h C4CustomART.h

OpticalSurfaceName starting '@'

C4OpBoundaryProcess Active DEFAULT QUADRANT

C4OpBoundaryProcess exactly same as G4OpBoundaryProcess, unless:

OpticalSurfaceName starts '@' or '#'
intersect onto solid "upper" (local_z > 0)

Standalone PMT tests : junosw extracts + Custom4 + Geant4

j:PMTSim Standalone PMT geometry and Optical Model

https://github.com/simoncblyth/j/tree/main/PMTSim
"virtual pkg" : reference junosw/blyth-88 branch sources
Depends on Custom4 mini-package, Geant4, opticks/SysRap

Quadrant Controls:

N=0/1 PMT geometry impl switch
POM=0/1 PMT Optical Model impl switch

Optical only Geant4 simulation, with full step point recording

opticks:u4/tests/U4SimulateTest.sh

N=0/1 POM=0/1 ./U4SimulateTest.sh

Geant4 Simtrace intersection : 2D geometry plotting

opticks:u4/tests/U4SimtraceTest.sh

N=0/1 ./U4SimtraceTest.sh

FewPMT_demo.png

two_pmt layout with HAMA, NNVT (Natural Geometry N=1)

FewPMT_demo0.png

Compare Unnatural/Natural N=0/1 geometry simulations by skipping N=0 Pyrex/Pyrex + Vacuum/Vacuum fake points

Statistical Comparison of Simulation Histories

	Geometry
A	N=0	Unnatural 4-volume Pyrex+Pyrex+Vacuum+Vacuum PMT Geometry (FastSim kludge)
B	N=1	Natural 2-volume Pyrex+Vacuum PMT Geometry (Using Custom4 C4OpBoundaryProcess)

As comparing simulations from different geometries : random aligned comparison is impractical, so:

skip recording fake points (Vacuum/Vacuum + Pyrex/Pyrex) -> make histories comparable
validate using statistical comparison of photon simulation histories : chi2 history matching

Initial comparisons :

HUGE Chi2
- very different history counts between simulations
find : ~2.5% of "whereAmI:Pyrex" junoPMTOpticalModel::ModelTrigger actually in Vacuum
- -> Unphysical reflect/refract/absorb/detect in middle of PMT vacuum
- unphysical N=0 behaviour -> huge chi2 in A/B comparison

Unphysical junoPMTOpticalModel::ModelTrigger FastSim Examples

FastSim junoPMTOpticalModel::ModelTrigger should happen:

on Vac/Vac inner1/inner2 border in middle of PMT
on Pyrex/Pyrex pmt/body border around upper hemi-ellipsoid

ModelTrigger elsewhere is a BUG -> unphysical results, including "tunneling" through geometry

Commandlines to create plots showing : NNVT Unphysical mid-Vacuum Reflect/Refract/Detect/Absorb

caused by "Undefined" behaviour of junoPMTOpticalModel::ModelTrigger

APID=2181 AOPT=idx N=0 SUBTITLE="9:reflect at vac/vac (but obscured)" ./U4SimtraceTest.sh ana

APID=2563 AOPT=idx N=0 SUBTITLE="7:SD at vac/vac" ./U4SimtraceTest.sh ana

APID=2912 AOPT=idx N=0 SUBTITLE="7:vac/vac refract 11:SD detect in vac" ./U4SimtraceTest.sh ana

APID=2942 AOPT=idx N=0 SUBTITLE="6->7:cross-geometry 7:vac/vac-reflect" ./U4SimtraceTest.sh ana

APID=2982 AOPT=idx N=0 SUBTITLE="7:unphysical refract" ./U4SimtraceTest.sh ana

APID=5866 AOPT=idx N=0 SUBTITLE="7:reflect at vac/vac" ./U4SimtraceTest.sh ana

APID=7625 AOPT=idx N=0 SUBTITLE="13:detect mid vacuum" ./U4SimtraceTest.sh ana

APID=8499 AOPT=idx N=0 SUBTITLE="13:detect mid vacuum" ./U4SimtraceTest.sh ana

FewPMT_2563_Unphysical_SD_in_vacuum.png

FewPMT_2912_Unphysical_refract_and_detect.png

FewPMT_2942_Unphysical_cross_geometry_vac_reflect.png

FewPMT_2982_Unphysical_vac_refract.png

FewPMT_5866_vac_reflect.png

FewPMT_7625_Unphysical_detect_mid_vacuum.png

FewPMT_8499_Unphysical_detect_mid_vacuum.png

Debug junoPMTOpticalModel::ModelTrigger with simple struct

 #include "NPX.h"

 struct ModelTrigger_Debug
 {
    static constexpr const char* NAME = "ModelTrigger_Debug.npy" ;
    static std::vector<ModelTrigger_Debug> RECORD ;
    static const int NUM_QUAD = 7 ;

    void add(){ RECORD.push_back(*this); }

    static NP* Array(){ return
    NPX::ArrayFromVec<double,ModelTrigger_Debug>(RECORD,NUM_QUAD,4);
    }

    static void Save(const char* savedir);

    double pos_x ;       // 00
    double pos_y ;       // 01
    double pos_z ;       // 02
    double time ;        // 03
    ...

    double   dist1 ;     // 20
    double   dist2 ;     // 21
    uint64_t mlv   ;     // 22
    uc8      etrig ;     // 23

    uint64_t index ;     // 30
    uint64_t pv ;        // 31
    uint64_t whereAmI ;  // 32
    uint64_t trig ;      // 33
    ...
 };

Simulation/DetSimV2/PMTSim/include/ModelTrigger_Debug.h

FewPMT_NNVT_ModelTriggerBuggy_inner2_Pyrex_trig.png

Reflection off NNVT innards => misplaced "Pyrex" ModelTrigger => FastSim too soon => Unphysical photons

FewPMT_HAMA_ModelTriggerBuggy_no_inner2_Pyrex_trig.png

Why no similar ModelTriggerBuggy problem for HAMA ? Undefined behaviour of G4UnionSolid::DistanceToIn from Inside

junoPMTOpticalModel::ModelTrigger : NNVT : ~2.5% UNPHYSICAL TRIG

u4/tests/U4SimulateTest.sh
u4/tests/U4SimulateTest_mt.py

Check ModelTrigger position, whereAmI : 271/11214 ~2.5% unphysical

mtd_trig = mtd.trig == 1
mtd_upper = mtd.pos[:,2] > 1e-4
mtd_mid   = np.abs( mtd.pos[:,2]) < 1e-4
mtd_lower = mtd.pos[:,2] < -1e-4

mtd_trig_vacuum = np.logical_and(mtd.trig == 1, mtd.whereAmI_ == 2 )
mtd_trig_vacuum_upper = np.logical_and(mtd_trig_vacuum, mtd_upper )
mtd_trig_vacuum_mid   = np.logical_and(mtd_trig_vacuum, mtd_mid )
mtd_trig_vacuum_lower = np.logical_and(mtd_trig_vacuum, mtd_lower )

mtd_trig_pyrex  = np.logical_and(mtd.trig == 1, mtd.whereAmI_ == 1 )
mtd_trig_pyrex_upper = np.logical_and(mtd_trig_pyrex, mtd_upper )
mtd_trig_pyrex_mid   = np.logical_and(mtd_trig_pyrex, mtd_mid )
mtd_trig_pyrex_lower = np.logical_and(mtd_trig_pyrex, mtd_lower )

NumPy mask expressions used for ModelTrigger_Debug analysis

Implications of whereAmI bug for NNVT ModelTrigger results

~/opticks/u4/tests/G4VSolid_Test.{sh,cc}

G4UnionSolid::DistanceToIn(DistanceToOut) from Inside(Outside) is "Undefined" BUT GEANT4 DOES NOT CHECK

> u4t ; ./G4VSolid_Test.sh
                                                             NNVT INNER2: Just Lower-Hemi-Ellipsoid
        pos    INNER1  DistToIn DistToOut      Dist        INNER2  DistToIn DistToOut      Dist     trig
   (0,0,20)   kInside  170.0000  170.0000  170.0000      kOutside kInfinity    0.0000 kInfinity      YES
   (0,0,15)   kInside  175.0000  175.0000  175.0000      kOutside kInfinity    0.0000 kInfinity      YES
   (0,0,10)   kInside  180.0000  180.0000  180.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,5)   kInside  185.0000  185.0000  185.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,4)   kInside  186.0000  186.0000  186.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,3)   kInside  187.0000  187.0000  187.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,2)   kInside  188.0000  188.0000  188.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,1)   kInside  189.0000  189.0000  189.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,0)  kSurface    0.0000  190.0000  190.0000      kSurface kInfinity    0.0000    0.0000      YES
   (0,0,-1)  kOutside    1.0000  191.0000    1.0000       kInside kInfinity    1.0000    1.0000      YES
   (0,0,-2)  kOutside    2.0000  192.0000    2.0000       kInside kInfinity    2.0000    2.0000      YES
   (0,0,-3)  kOutside    3.0000  193.0000    3.0000       kInside kInfinity    3.0000    3.0000      YES
   (0,0,-4)  kOutside    4.0000  194.0000    4.0000       kInside kInfinity    4.0000    4.0000      YES
   (0,0,-5)  kOutside    5.0000  195.0000    5.0000       kInside kInfinity    5.0000    5.0000      YES
   (0,0,-6)  kOutside    6.0000  196.0000    6.0000       kInside kInfinity    6.0000    6.0000      YES
   (0,0,-7)  kOutside    7.0000  197.0000    7.0000       kInside kInfinity    7.0000    7.0000      YES
   (0,0,-8)  kOutside    8.0000  198.0000    8.0000       kInside kInfinity    8.0000    8.0000      YES
   (0,0,-9)  kOutside    9.0000  199.0000    9.0000       kInside kInfinity    9.0000    9.0000      YES
  (0,0,-10)  kOutside   10.0000  200.0000   10.0000       kInside kInfinity   10.0000   10.0000      YES

DistanceToIn, DistanceToOut, Distance_ in +Z direction from points along Z axis
CSG constituent structure of NNVT inner1/inner2 yield mis-placed Pyrex triggers for Z > -185 mm

Distance uses DistanceToOut/In depending on G4VSolid::Inside

Returns distance to first intersect from position pos in direction dir, just like a ray tracer

inline G4double Distance(
     const G4VSolid* solid,
     const G4ThreeVector& pos, const G4ThreeVector& dir, EInside& in )
{
    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 ;
}

Avoid "Undefined" result by using appropriate DistanceToOut or DistanceToIn method.

Implications of whereAmI bug for HAMA ModelTrigger results

DistanceToIn, DistanceToOut, Distance in +Z direction from points along Z axis

> u4t ; ./G4VSolid_Test.sh
                                                             HAMA INNER2: Union of Polycone and Lower-Hemi-Ellipsoid
        pos    INNER1  DistToIn DistToOut      Dist        INNER2  DistToIn DistToOut      Dist     trig
   (0,0,20)   kInside  170.0000  170.0000  170.0000      kOutside kInfinity    0.0000 kInfinity      YES
   (0,0,15)   kInside  175.0000  175.0000  175.0000      kOutside kInfinity    0.0000 kInfinity      YES
   (0,0,10)   kInside  180.0000  180.0000  180.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,5)   kInside  185.0000  185.0000  185.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,4)   kInside  186.0000  186.0000  186.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,3)   kInside  187.0000  187.0000  187.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,2)   kInside  188.0000  188.0000  188.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,1)   kInside  189.0000  189.0000  189.0000      kOutside kInfinity    0.0000 kInfinity      YES
    (0,0,0)  kSurface    0.0000  190.0000  190.0000      kSurface kInfinity    0.0000    0.0000      YES
   (0,0,-1)  kOutside    1.0000  191.0000    1.0000       kInside kInfinity    1.0000    1.0000      YES
   (0,0,-2)  kOutside    2.0000  192.0000    2.0000       kInside kInfinity    2.0000    2.0000      YES
   (0,0,-3)  kOutside    3.0000  193.0000    3.0000       kInside kInfinity    3.0000    3.0000      YES
   (0,0,-4)  kOutside    4.0000  194.0000    4.0000       kInside kInfinity    4.0000    4.0000      YES
   (0,0,-5)  kOutside    5.0000  195.0000    5.0000       kInside    0.0000    5.0000    5.0000      NO
   (0,0,-6)  kOutside    6.0000  196.0000    6.0000       kInside    1.0000    6.0000    6.0000      NO
   (0,0,-7)  kOutside    7.0000  197.0000    7.0000       kInside    2.0000    7.0000    7.0000      NO
   (0,0,-8)  kOutside    8.0000  198.0000    8.0000       kInside    3.0000    8.0000    8.0000      NO
   (0,0,-9)  kOutside    9.0000  199.0000    9.0000       kInside    4.0000    9.0000    9.0000      NO
  (0,0,-10)  kOutside   10.0000  200.0000   10.0000       kInside    5.0000   10.0000   10.0000      NO

HAMA INNER2 has equatorial cylinder from Z=-5 mm to Z=0 mm (NNVT has no such cylinder)
The equatorial cylinder changes the "Undefined" behavior preventing triggers below Z=-5 mm
So ModelTrigger bug is avoided for HAMA by "Undefined" behaviour (NOT GOOD)

junoPMTOpticalModel::ModelTrigger BUG => Unphysical Photons

+----------*----pmt-Pyrex----------------+
|                                        |
|                                        |
|     +------*-------body-Pyrex----+     |
|     | +-----*-------inner1-Vac-+ |     |
|     | |      \                 | |     |
|     | |       \                | |     |
|     | |        \               |-|     |
|     | |         \              |1e-3   |
|     | |          \             | |     |
|     | +~~vac/vac~~\~~~~~!~~~~~~+ |     |
|     | |            \   / \     | |     |
|     | |             \ /   \    | |     |
|     | |     +--------*---+ \   | |     |
|     | |     | dynode/MCP |  \  | |     |
|     | |     +------------+   \ | |     |
|     | +-inner2-Vac------------*+ |     |
|     +----------------------------+     |
|                                        |
|                                        |
+----------------------------------------+

dyn/MCP reflect => whereAmI:kInGlass
-> unphysical mid Vacuum A/R/T/D
"Undefined" G4VSolid::DistanceToIn

  G4bool junoPMTOpticalModel::ModelTrigger(const G4FastTrack &fastTrack)
  {   // GetVolume : CAUTION NEEDED AT BORDERS 
    ...
    if(fastTrack.GetPrimaryTrack()->GetVolume() == _inner2_phys){
        return false;
    }   // SIMPLISTIC VIEW OF GEOMETRY 
    if(fastTrack.GetPrimaryTrack()->GetVolume() == _inner1_phys){
        whereAmI = kInVacuum;
    }else{
        whereAmI = kInGlass;
        // BUG: REFLECTION FROM dynode/MCP -> kInGlass 
    }
    if(whereAmI == kInGlass){
        dist1 = _inner1_solid->DistanceToIn(pos, dir);
        dist2 = _inner2_solid->DistanceToIn(pos, dir);
        // BUG: dist2 "Undefined" at reflect inside inner2
        if(dist1 == kInfinity){ return false;
        }else if(dist1>dist2){  return false;
        }else{                  return true;}
        // BUG: depends on "Undefined" behaviour (+CSG constituents) 
    }else{
        dist1 = _inner1_solid->DistanceToOut(pos, dir);
        dist2 = _inner2_solid->DistanceToIn(pos, dir);
        if(dist2 == kInfinity){ return true;}
         // reliance on edge handling : not a good approach
    }
    return false;
  }

junoPMTOpticalModel::ModelTrigger FIXED

G4bool ModelTrigger(const G4FastTrack &fastTrack)
{
  bool trig = false ;
  dist1 = Distance( _inner1_solid, pos, dir, in1 );
  if( dist1 != kInfinity && dist1 > EPSILON )
  {
    next_pos = pos + dir*dist1 ;
    trig = next_pos.z() > EPSILON ;
    // disqualify flat face intersects

    next_norm = _inner1_solid->SurfaceNormal(next_pos);
    next_mct  = next_norm * dir ;
    whereAmI  = next_mct < 0. ? kInGlass : kInVacuum ;
    // against/with normal is Glass/Vacuum
  }
  return trig ;
}
const G4double EPSILON = 2e-4 ;

Use: Distance, SurfaceNormal -- NOT Volumes

 G4bool junoPMTOpticalModel::ModelTrigger(const G4FastTrack &fastTrack)
 {   // GetVolume : NOT RELIABLE AT BORDERS 
   ...
   if(fastTrack.GetPrimaryTrack()->GetVolume() == _inner2_phys){
       return false;
   }   // SIMPLISTIC VIEW OF GEOMETRY 
   if(fastTrack.GetPrimaryTrack()->GetVolume() == _inner1_phys){
       whereAmI = kInVacuum;
   }else{
       whereAmI = kInGlass;
       // BUG: REFLECTION FROM dynode/MCP -> kInGlass 
   }
   if(whereAmI == kInGlass){
       dist1 = _inner1_solid->DistanceToIn(pos, dir);
       dist2 = _inner2_solid->DistanceToIn(pos, dir);
       // BUG: dist2 "Undefined" at reflect inside inner2
       if(dist1 == kInfinity){ return false;
       }else if(dist1>dist2){  return false;
       }else{                  return true;}
       // BUG: depends on "Undefined" behaviour (+CSG constituents) 
   }else{
       dist1 = _inner1_solid->DistanceToOut(pos, dir);
       dist2 = _inner2_solid->DistanceToIn(pos, dir);
       if(dist2 == kInfinity){ return true;}
        // reliance on edge handling : not a good approach
   }
   return false;
 }

Need to fix FastSim approach : to have something to validate new natural geometry approach with

FewPMT_NNVT_ModelTriggerSimple_avoids_issue.png

ModelTriggerSimple : Avoids Undefined trig behavior + Unphysical mid-Vacuum Reflect/Refract/Detect/Absorb, "Tunneling"

FewPMT_check_mtd_trig_vacuum_upper_due_to_obliques.png

mtd_trig_vacuum_upper initially seems surprising : but actually it is expected from "oblique" py->vac paths

quiver plot (mtd_pos, mtd_dir) with scatter plot of mtd_pos and mtd_next_pos

N0_one_pmt_FewPMT_hamaLogicalPMT_rain_line_012last.png

N=0 Unnatural : HAMA "rain_line" photon positions : 0th:blue 1st:red last:cyan

N1_one_pmt_FewPMT_hamaLogicalPMT_rain_line_012last.png

N=1 Natural : HAMA "rain_line" photon positions : 0th:blue 1st:red last:cyan

NNVT : "rain_line" Example of Statistical History Comparison

 ./U4SimulateTest.sh cf ## PMT Geometry : A(N=0) Unnatural+FastSim, B(N=1) Natural+CustomBoundary
 GEOM/GEOMList/IMPL/LAYOUT/CHECK : FewPMT/nnvtLogicalPMT/ModelTriggerSimple/one_pmt/rain_line
 c2sum :    60.2333 c2n :    67.0000 c2per:     0.8990  C2CUT:   30 

 np.c_[siq,_quo,siq,sabo2,sc2,sabo1][:25]  ## A-B history frequency chi2 comparison
 [[' 0' TO BT SD                             ' 37975  37525' ' 2.6821' '  1268   1268']
  [' 1' TO BT SA                             ' 32161  32735' ' 5.0770' '  1174   1200']
  [' 2' TO BT BT SR SA                       '  4980   4998' ' 0.0325' ' 11979  12125']
  [' 3' TO BT BT SA                          '  4392   4362' ' 0.1028' ' 10275  10285']
  [' 4' TO BT BT SR BT BT SA                 '  4366   4366' ' 0.0000' ' 21394  21384']
  [' 5' TO BT BR BT SA                       '  3225   3153' ' 0.8128' '  1286   1269']
  [' 6' TO SA                                '  2281   2283' ' 0.0009' '     0      0']
  [' 7' TO BT BT SR BR SR SA                 '  1185   1280' ' 3.6613' ' 22175  22191']
  [' 8' TO BT BT SR BR SA                    '  1151   1150' ' 0.0004' ' 21484  21399']
  [' 9' TO BR SA                             '  1051   1091' ' 0.7470' '  1154   1154']
  ['10' TO BT BT SR SR SA                    '   742    730' ' 0.0978' ' 11888  11906']
  ['11' TO BT BT SR BR SR BT BT SA           '   714    694' ' 0.2841' ' 22163  22168']
  ['12' TO BT BT SR BR SR SR SA              '   396    324' ' 7.2000' ' 22136  23759']
  ['13' TO BT BT SR SR BT BT SA              '   368    357' ' 0.1669' ' 18370  18360']
  ['14' TO BT AB                             '   341    367' ' 0.9548' '  1483   1417']
  ['15' TO BT BT SR BR SR SR BT BT SA        '   363    338' ' 0.8916' ' 22134  22135']
  ['16' TO BT BT SR SR SR SA                 '   285    276' ' 0.1444' ' 11903  12085']
  ['17' TO BT BT SR SR SR BT BT SA           '   239    238' ' 0.0021' ' 13690  13690']
  ['18' TO BT BT SR BR SR BR SR BT BT SA     '   225    236' ' 0.2625' ' 22345  22388']

                                                   A      B    (A-B)^2   APID[0] BPID[0]
                                                             -------
                                                               (A+B)

Compare simulations by chi2 comparison of photon history counts

TO:Torch, BT:BoundaryTransmit, SD:SurfDetect, SA:SurfAbsorb, SR:SpecularReflect, BR:FresnelReflect, AB:BulkAbsorb
chi2 value indicates the two simulations have consistent photon histories

HAMA : "rain_line" Example of Statistical History Comparison

 ./U4SimulateTest.sh cf ## PMT Geometry : A(N=0) Unnatural+FastSim, B(N=1) Natural+CustomBoundary
 GEOM/GEOMList/IMPL/LAYOUT/CHECK : FewPMT/hamaLogicalPMT/ModelTriggerSimple/one_pmt/rain_line
 c2sum :    41.2803 c2n :    43.0000 c2per:     0.9600  C2CUT:   30

 np.c_[siq,_quo,siq,sabo2,sc2,sabo1][:25]  ## A-B history frequency chi2 comparison
 [[' 0' 'TO BT SD                        '' 37923  37958' ' 0.0161' '  1262   1274']
  [' 1' 'TO BT SA                        '' 32806  32562' ' 0.9108' '  1185   1207']
  [' 2' 'TO BT BT SA                     '' 12071  12233' ' 1.0798' ' 10532  10430']
  [' 3' 'TO BT BT SR SA                  ''  3757   3741' ' 0.0341' ' 12754  12407']
  [' 4' 'TO BT BR BT SA                  ''  3204   3250' ' 0.3279' '  1465   1464']
  [' 5' 'TO SA                           ''  2277   2274' ' 0.0020' '     1      0']
  [' 6' 'TO BT BT SR SR SA               ''  1857   1859' ' 0.0011' ' 12391  12393']
  [' 7' 'TO BR SA                        ''  1129   1215' ' 3.1553' '  1154   1154']
  [' 8' 'TO BT BT SR SR SR BT BT SA      ''   613    580' ' 0.9128' ' 16112  16120']
  [' 9' 'TO BT BT SR SR SR SA            ''   584    576' ' 0.0552' ' 13493  16172']
  ['10' 'TO BT BT SR BT BT SA            ''   357    375' ' 0.4426' ' 34234  34225']
  ['11' 'TO BT AB                        ''   356    325' ' 1.4112' '  1367   1317']
  ['12' 'TO BT BT SR BR SA               ''   287    305' ' 0.5473' ' 34230  34221']
  ['13' 'TO BT BT SR SR BT BT SA         ''   243    226' ' 0.6162' ' 27497  27530']
  ['14' 'TO BT BT SR SR SR BR SA         ''   232    206' ' 1.5434' ' 16755  16130']
  ['15' 'TO BT BT SR SR SR BR BT BT SA   ''   200    206' ' 0.0887' ' 16715  16760']
  ['16' 'TO BT BT BR SR SA               ''   181    179' ' 0.0111' ' 10909  10984']
  ['17' 'TO BT BT SR SR SR BR SR SA      ''   146    180' ' 3.5460' ' 16223  16205']
  ['18' 'TO BT BT SR SR SR BR BR SR SA   ''   148    131' ' 1.0358' ' 17957  16775']
  ['19' 'TO AB                           ''   147    126' ' 1.6154' '    86     68']
  ['20' 'TO BT BT SR SR SR BR BR SA      ''   113    111' ' 0.0179' ' 16724  16934']

                                                A      B    (A-B)^2   APID[0] BPID[0]
                                                           -------
                                                            (A+B)

Same for NNVT with this check
chi2 value indicates the two simulations have consistent photon histories

Checks with different initial photon pos/dir/shapes

opticks/u4/tests/U4SimulateTest.sh checks

FewPMT/(hama,nnvt)LogicalPMT/ModelTriggerSimple/one_pmt/... using 100k/50k photons

check		c2sum/c2n:c2per(C2CUT)
check		HAMA \| NNVT
rain_disc	3D beam directly onto PMT front face	32.56/41 0.794 (30)	41.91/46:0.911 (30)
rain_line	2D line directly onto PMT midline	41.28/43 0.960 (30)	28.85/45:0.641 (30)
up_rain_line	2D line shooting PMT from behind	10.80/24 0.450 (30)	6.41/12:0.534 (30)
rain_dynode	2D line from vacuuum directed onto dynode	12.08/11 1.098 (30)	23.76/35:0.679 (30)
rain_dynode_diag	2D line diagonnally directed from mid-vacuum	8.31/26 0.320 (30)	42.82/56:0.765 (30)
lhs_window_line	2D line from left directed onto PMT window	37.53/34 1.104 (30)	26.74/43:0.622 (30)
lhs_reflector_line	2D line from left directed onto PMT reflector	14.35/18 0.797 (30)	15.75/19:0.829 (30)

Obtained using opticks/u4/tests/cf.sh (runs both simulations and compares)

Reasonable Chi2 for all checks => consistent simulation histories with Unnatural/Natural geometry

N0_one_pmt_FewPMT_hamaLogicalPMT_lhs_window_line_inwin.png

(N=0) Shoot from side to check PMT window : many photons enter PMT

N1_one_pmt_FewPMT_hamaLogicalPMT_lhs_window_line_inwin.png

(N=1) Shoot from side to check PMT window : many photons enter PMT

N0_one_pmt_FewPMT_hamaLogicalPMT_lhs_reflector_line_reflected.png

(N=0) Shoot from side to check PMT reflector : no photons get into PMT

N1_one_pmt_FewPMT_hamaLogicalPMT_lhs_reflector_line_reflected.png

(N=1) Shoot from side to check PMT reflector : no photons get into PMT

Standalone then Insitu Testing

Standalone Test Architecture

small "FewPMT" geometries
Opticks executable such as U4SimulateTest

j/PMTSim "virtual" package : PMT Geometry+Model

github.com/simoncblyth/j/blob/main/PMTSim/CMakeLists.txt
CMakeLists.txt references junosw sources

No Known Issues In Standalone Tests

=> move to insitu tests with junosw tut_detsim.py

Insitu Testing with: --opticks-mode 2

junosw Geant4 simulation
integrated with Opticks Instrumentation (U4Recorder, SEvt)
no GPU usage in this mode

junosw + Opticks : Input Photon Running : ntds bash function

 ntds(){
     local opts=""
     opts="$opts --opticks-mode 2"  # G4 with Opticks Instrumentation

     case $POM in     ## passed into UsePMTOpticalModel
       0) opts="$opts --no-pmt-optical-model"  ;;
       1) opts="$opts --pmt-optical-model"     ;;
     esac

     case $PMT in  ## passed into UsePMTNaturalGeometry
       0) opts="$opts --pmt-unnatural-geometry" ;;
       1) opts="$opts --pmt-natural-geometry"   ;;
     esac

     opts="$opts --opticks-anamgr"  ## Full Photon Point recording
     export OPTICKS_EVENT_MODE=StandardFullDebug
     export OPTICKS_MAX_BOUNCE=31
     [ "$PMT" == "0" ] && export U4Recorder__FAKES_SKIP=1

     export OPTICKS_INPUT_PHOTON=RainXZ_Z230_10k_f8.npy
     export OPTICKS_INPUT_PHOTON_FRAME=Hama:0:1000 

     python tut_detsim.py $opts opticks  ## GtOpticksTool
}

Actual ntds bash function : https://github.com/simoncblyth/j/blob/main/jx.bash

Geant4 + Opticks + NVIDIA OptiX 7 : Hybrid Workflow

Geant4 + Opticks Instrumentation (no GPU usage)

U4Recorder : uses CPU Opticks to persist SEvt photons/histories/.. into NumPy .npy files for analysis

Opticks SEvt loaded into ipython : arrays detailing every photon point

In [3]: b
 Out[3]:
 b.base:/tmp/blyth/opticks/GEOM/V1J008/ntds2/ALL1/000
   : b.genstep                                          :            (1, 6, 4) : 14:01:49.017610
   : b.seq                                              :       (100000, 2, 2) : 13:59:54.715319
   : b.pho0                                             :          (100000, 4) : 14:01:40.913446
   : b.record                                           :   (100000, 32, 4, 4) : 13:59:55.174871
   : b.junoSD_PMT_v2                                    :                 (1,) : 14:01:44.544258
   : b.sframe                                           :            (4, 4, 4) : 13:59:54.714881
   : b.inphoton                                         :       (100000, 4, 4) : 14:01:44.544879
   : b.pho                                              :          (100000, 4) : 14:01:42.596733
   : b.flat                                             :         (100000, 64) : 14:01:49.018186
   : b.prd                                              :   (100000, 32, 2, 4) : 14:01:26.149970
   : b.photon                                           :       (100000, 4, 4) : 14:01:26.466269
   : b.gs                                               :               (1, 4) : 14:01:49.017187
   : b.aux                                              :   (100000, 32, 4, 4) : 14:01:49.191892
   : b.tag                                              :          (100000, 4) : 13:59:54.705722

array	shape	notes
inphoton	(n,4,4)	input photons
photon	(n,4,4)	final photon positions (at AB, SA, SD)
record	(n,32,4,4)	photon histories (up to 32 points)
seq	(n,2,2)	uint64 packed histories eg "TO BT BT SD"
aux	(n,32,4,4)	extra step point info eg A,R,T
sframe	(4,4,4)	target frame with M2W W2M transforms

B_V1J008_N1_ip_MOI_Hama:0:1000_yy_target_frame.png

cd ~/j/ntds ; N=1 GLOBAL=1 MODE=3 ./ntds.sh ana ## interactive 3D pyvista

Target Frame Hama:0:1000

solidName:ordinalIdx:instanceIdx

B_V1J008_N1_ip_MOI_Hama:0:1000_yy_frame_close.png

Green : start position (100k input photons)

Red : end position, Cyan : other position

Transform photon history points into target frame

NumPy transform all photon positions to target frame (eg Hama:0:1000)

target frame m2w w2m transforms persisted within opticks/SEvt folder

w2m = evt.f.sframe.w2m          # (4,4)         : world-to-model
gpos = b.f.record[:,:,0].copy() # (100000,32,4) : global (pos, time)
gpos[...,3] = 1.                # use 1. : to transform as position
lpos = np.dot( gpos, w2m )      # (100000,32,4) : local (pos, 1)

Suitable : input photon shape, position, direction, target frame

-> action mostly in plane -> 2D plotting (often clearer than 3D)

cd ~/j/ntds ; MODE=2 ./ntds.sh ana

https://github.com/simoncblyth/j/blob/main/ntds/ntds.sh

A_V0J008_N0_ip_MOI_Hama:0:1000_a_no_skip_fake.png

cd ~/j/ntds ; N=0 ./ntds.sh ana

10k : Unnatural Pyrex+Pyrex+Vac+Vac PMT geometry, without : export U4Recorder__FAKES_SKIP=1

A_V0J008_N0_ip_MOI_Hama:0:1000_a.png

cd ~/j/ntds ; N=0 ./ntds.sh ana

10k : Unnatural Pyrex+Pyrex+Vac+Vac HAMA PMT geometry, with : export U4Recorder__FAKES_SKIP=1

B_V1J008_N1_ip_MOI_Hama:0:1000_b.png

cd ~/j/ntds ; N=1 ./ntds.sh ana

10k : Natural Pyrex+Vac HAMA PMT geometry, no fakes, no skipping needed

100k Input Photons targetting Hama:0:1000 Chi2 History Comparison

cd ~/j/ntds ; ./ntds.sh cf
QCF qcf c2sum/c2n:c2per(C2CUT)  103.55/103:1.005 (30)
np.c_[siq,_quo,siq,sabo2,sc2,sabo1][:25]  ## A-B history frequency chi2 comparison
[[' 0' 'TO BT BT BT BT SA                                  ' ' 0' ' 37664  37569'   0.1200 '     2      1']
 [' 1' 'TO BT BT BT BT SD                                  ' ' 1' ' 30711  30904'   0.6045 '     4      2']
 [' 2' 'TO BT BT BT BT BT SA                               ' ' 2' ' 12432  12470'   0.0580 '  9199   8859']
 [' 3' 'TO BT BT BT BT BT SR SA                            ' ' 3' '  3824   3761'   0.5233 ' 11002  11014']
 [' 4' 'TO BT BT BT BT BT SR SR SA                         ' ' 4' '  1984   1967'   0.0731 ' 10886  10878']
 [' 5' 'TO BT BT AB                                        ' ' 5' '   853    918'   2.3857 '     8     27']
 [' 6' 'TO BT BT BT BT BT SR SR SR SA                      ' ' 6' '   583    609'   0.5671 ' 14779  14727']
 [' 7' 'TO BT BT BT BT BR BT BT BT BT BT SA                ' ' 7' '   444    420'   0.6667 '  1042   1016']
 [' 8' 'TO BT BT BT BT BR BT BT BT BT BT BT AB             ' ' 8' '   435    430'   0.0289 ' 10280   7615']
 [' 9' 'TO BT BT BT BT BR BT BT BT BT BT SD                ' ' 9' '   354    328'   0.9912 '  5253   5256']
 ['10' 'TO BT BT BT BT AB                                  ' '10' '   331    341'   0.1488 '    84    143']
 ['11' 'TO BT BT BT BT BT SR BR SA                         ' '11' '   258    314'   5.4825 ' 33580  33573']
 ['12' 'TO BT BT BR BT BT BT SA                            ' '12' '   295    250'   3.7156 '     0      3']
 ['13' 'TO BT BT BT BR BT BT BT BT SA                      ' '13' '   291    259'   1.8618 '   213    192']
 ['14' 'TO BT BT BT BT BR BT BT BT BT AB                   ' '14' '   289    275'   0.3475 '  8511   9213']
 ['15' 'TO BT BT BT BT BR BT BT BT BT BT BT BT BT SA       ' '15' '   245    276'   1.8445 '  9353   9360']
 ['16' 'TO BT BT BT BT BR BT BT BT BT BT BT BT BT SD       ' '16' '   261    246'   0.4438 '  9350   9416']
 ['17' 'TO BT BT BT BT BT SR SR SR BR BT BT BT BT BT BT SA ' '17' '   247    199'   5.1659 ' 15391  15388']
 ['18' 'TO BT BT BT BT BT SR SR SR BR SA                   ' '18' '   246    227'   0.7632 ' 14754  15177']
 ['19' 'TO BT BT BT BT BT SR SR BT BT BT BT BT BT SA       ' '19' '   215    194'   1.0782 ' 26564  26588']
                                                                       A,B counts    Chi2   A,B 1st photon idx

https://github.com/simoncblyth/j/blob/main/ntds/ntds.sh

Chi2 comparison of photon histories for A:Unnatural(fakes skipped) B:Natural geometry

Any bug changing simulation histories -> huge Chi2 -> USEFUL METRIC

BUT: does not catch all issues (eg time offsets)

Instrument junoSD_PMT_v2::ProcessHits : Compare N=0,1 enum counts

359 #ifdef WITH_G4CXOPTICKS
360 G4bool junoSD_PMT_v2::ProcessHits(G4Step* step,...)
361 {
362     G4bool is_hit = ProcessHits_(step, nullptr) ;
363     m_jpmt_dbg->add( m_eph, is_hit );
373     return is_hit ;
374 }
382 #endif
384 #ifdef WITH_G4CXOPTICKS
385 G4bool junoSD_PMT_v2::ProcessHits_(G4Step* step,...)
386 #else
387 G4bool junoSD_PMT_v2::ProcessHits(G4Step* step,...)
388 #endif
389 {
390     if (m_disable) {
391 #ifdef WITH_G4CXOPTICKS
392         m_eph = EPH::NDIS ;
393 #endif
394         return false;
395     }
...

Counts stored into SEvt metadata

cd ~/j/ntds ; NEVT=3 ./ntds.sh cfmeta
array([[[31967, 31903, 32017],  N=0   0:ProcessHits_count
        [33305, 33464, 33416]], N=1

       [[ 8732,  8693,  8745],        1:ProcessHits_true
        [ 8696,  8772,  8717]],

       [[23235, 23210, 23272],        2:ProcessHits_false
        [24609, 24692, 24699]],

       [[ 8732,  8693,  8745],        3:SaveNormHit_count
        [ 8696,  8772,  8717]],
      ...

       [[ 2155,  2149,  2086],        8:NEDEP
        [ 3443,  3577,  3552]],

       [[  134,   114,   130],        9:NBOUND
        [  134,   131,   129]],
       ...

       [[20946, 20947, 21056],       12:NDECULL
        [21032, 20984, 21018]],
       ...

       [[ 8732,  8693,  8745],       14:YSAVE
        [ 8696,  8772,  8717]],

       [[    2,     2,     2],       15:opticksMode
        [    2,     2,     2]],

       [[    0,     1,     2],       16:eventID
        [    0,     1,     2]],

       [[    0,     0,     0],       17:VERSION
        [    1,     1,     1]]], dtype=uint64)

HAMA InputPhoton Chi2 (10 evt) `~/j/ntds/ntds.sh cfmeta`

msab.c_itab
array([[[15617, 15661, 15635, 15890, 15447, 15871, 15500, 15824, 15598, 15627],        0:ProcessHits_count  N=0 
        [18490, 18210, 18490, 18174, 18467, 18069, 18608, 18198, 18291, 18284]],                            N=1 

       [[ 3083,  3070,  3035,  2978,  3020,  2966,  2976,  2944,  3044,  3029],        1:ProcessHits_true
        [ 3042,  3001,  2904,  2958,  2925,  2996,  3023,  3007,  3068,  2991]],

       [[12534, 12591, 12600, 12912, 12427, 12905, 12524, 12880, 12554, 12598],        2:ProcessHits_false
        [15448, 15209, 15586, 15216, 15542, 15073, 15585, 15191, 15223, 15293]],
       ...

msab.c_itab[msab.ik.YSAVE,0].sum()/msab.c_itab[msab.ik.YSAVE,1].sum()     1.0077  # ratio of YSAVE counts N=0/N=1  
msab.c_ftab[0,1].sum()/msab.c_ftab[0,0].sum()                             0.8472  # ratio of durations N=1/N=0  
np.diff(msab.c_ttab)[0,1].sum()/np.diff(msab.c_ttab)[0,0].sum()           0.8582  # ratio of evt times N=1/N=0  

np.diff(msab.c_ttab)/1e6                                                  # seconds between event starts 
array([[[2.589, 2.983, 3.004, 3.01 , 2.994, 2.665, 2.619, 2.641, 2.61 ],
        [2.304, 2.526, 2.453, 2.443, 2.495, 2.144, 2.179, 2.478, 2.531]]])      Little difference in durations, evt times

np.c_[msab.c_ttab[0].T].view('datetime64[us]') # SEvt start times (UTC)  opticks/sysrap/stimer.h (std::chrono)
array([['2023-04-15T17:35:25.726622', '2023-04-15T17:41:01.654181'],
       ['2023-04-15T17:35:28.315248', '2023-04-15T17:41:03.957860'],
       ...
       ['2023-04-15T17:35:50.840385', '2023-04-15T17:41:23.207810']], dtype='datetime64[us]')

msab.c2tab  # c2sum, c2n, c2per for each event
array([[28.08 , 21.765, 24.036, 14.987, 12.353, 29.081, 23.99 , 18.305, 17.485, 20.135],
       [25.   , 25.   , 26.   , 25.   , 25.   , 27.   , 25.   , 25.   , 27.   , 28.   ],
       [ 1.123,  0.871,  0.924,  0.599,  0.494,  1.077,  0.96 ,  0.732,  0.648,  0.719]])  # c2per for 10 evts 

msab.c2tab[0,:].sum()/msab.c2tab[1,:].sum()      0.8148             # c2per_tot

N=0,1 History Chi2 (target HAMA) : Consistent across 10 evt Overall Chi2 : 0.8148

A_V0J008_N0_ip_MOI_NNVT:0:1000_10k.png

cd ~/j/ntds ; N=0 ./ntds.sh ana

10k : Unnatural Pyrex+Pyrex+Vac+Vac NNVT PMT geometry, with : export U4Recorder__FAKES_SKIP=1

B_V1J008_N1_ip_MOI_NNVT:0:1000_10k.png

cd ~/j/ntds ; N=1 ./ntds.sh ana

10k : Natural Pyrex+Vacuum NNVT PMT geometry, no fakes, no skipping needed

100k Input Photons targetting NNVT:0:1000 Chi2 History Comparison

cd ~/j/ntds ; ./ntds.sh cf
QCF qcf : c2sum/c2n:c2per(C2CUT) 197.53/194:1.018 (30)
np.c_[siq,_quo,siq,sabo2,sc2,sabo1][:25]  ## A-B history frequency chi2 comparison
[[' 0' 'TO BT BT BT BT SD                                          ' ' 0' ' 33025  33400' ' 2.1170' '     0      3']
 [' 1' 'TO BT BT BT BT SA                                          ' ' 1' ' 28317  27980' ' 2.0173' '     3      4']
 [' 2' 'TO BT BT BT BT BT SR SA                                    ' ' 2' '  6205   6328' ' 1.2071' ' 10427  10543']
 [' 3' 'TO BT BT BT BT BT SA                                       ' ' 3' '  4684   4693' ' 0.0086' '  8523   8481']
 [' 4' 'TO BT BT BT BT BT SR BR SR SA                              ' ' 4' '  1188   1174' ' 0.0830' ' 21107  21015']
 [' 5' 'TO BT BT BT BT BR BT BT BT BT BT BT AB                     ' ' 5' '  1008    961' ' 1.1219' '  7713   9959']
 [' 6' 'TO BT BT BT BT BT SR BR SA                                 ' ' 6' '   960    944' ' 0.1345' ' 20304  20311']
 [' 7' 'TO BT BT BT BT BT SR SR SA                                 ' ' 7' '   951    900' ' 1.4052' ' 10387  10424']
 [' 8' 'TO BT BT AB                                                ' ' 8' '   904    828' ' 3.3349' '    14    126']
 [' 9' 'TO BT BT BT BT BT SR BT BT BT BT BT BT BT AB               ' ' 9' '   629    587' ' 1.4507' ' 20503  20814']
 ['10' 'TO BT BT BT BT BR BT BT BT BT AB                           ' '10' '   590    585' ' 0.0213' '    86   8872']
 ['11' 'TO BT BT BT BT BR BT BT BT BT BT BT BT BT SD               ' '11' '   445    462' ' 0.3186' '  9137   9101']
 ['12' 'TO BT BT BT BT BR BT BT BT BT BT BT BT BT BT BT BT BT SD   ' '12' '   445    435' ' 0.1136' ' 11465   7184']
 ['13' 'TO BT BT BT BT BR BT BT BT BT BT SA                        ' '13' '   445    437' ' 0.0726' '   666    675']
 ['14' 'TO BT BT BT BT BR BT BT BT BT BT BT BT BT SA               ' '14' '   423    426' ' 0.0106' '  9090   9118']
 ['15' 'TO BT BT BT BT BR BT BT BT BT BT BT BT BT BT BT BT BT SA   ' '15' '   387    414' ' 0.9101' ' 11478   7091']
 ['16' 'TO BT BT BT BT BR BT BT BT BT BT BT SC BT BT BT BT BT BT SD' '16' '   381    404' ' 0.6739' ' 16538  16757']
 ['17' 'TO BT BT BT BT BT SR BR SR SR SA                           ' '17' '   335    383' ' 3.2089' ' 22766  20996']
 ['18' 'TO BT BT BT BT BR BT BT BT BT BT BT SC BT BT BT BT BT BT SA' '18' '   356    327' ' 1.2313' ' 16336  16142']
 ['19' 'TO BT BT BT BT BT SR SR SR SA                              ' '19' '   354    348' ' 0.0513' ' 10746  10369']

https://github.com/simoncblyth/j/blob/main/ntds/ntds.sh

Chi2 comparison of photon histories for:

Unnatural geometry with fakes skipped
Natural geometry

NNVT InputPhoton Chi2 (10 evt) `~/j/ntds/ntds.sh cfmeta`

msab.c_itab
array([[[22063, 22278, 21789, 21762, 21833, 22246, 21838, 21914, 22052, 21912],        0:ProcessHits_count
        [25986, 25956, 25693, 25834, 26009, 25528, 25944, 25642, 25787, 25522]],

       [[ 3512,  3517,  3541,  3496,  3576,  3434,  3607,  3484,  3542,  3552],        1:ProcessHits_true
        [ 3541,  3529,  3510,  3535,  3491,  3534,  3601,  3506,  3562,  3529]],

       [[18551, 18761, 18248, 18266, 18257, 18812, 18231, 18430, 18510, 18360],        2:ProcessHits_false
        [22445, 22427, 22183, 22299, 22518, 21994, 22343, 22136, 22225, 21993]],
...

msab.c_itab[msab.ik.YSAVE,0].sum()/msab.c_itab[msab.ik.YSAVE,1].sum()     0.9978       # N=0/N=1
msab.c_ftab[0,1].sum()/msab.c_ftab[0,0].sum()                             0.8766       # ratio of durations N=1/N=0
np.diff(msab.c_ttab)[0,1].sum()/np.diff(msab.c_ttab)[0,0].sum()           0.9107       # ratio of evt times N=1/N=0

np.diff(msab.c_ttab)/1e6   # seconds between event starts
array([[[3.086, 3.643, 3.474, 3.063, 3.012, 3.149, 3.164, 3.192, 3.177],
        [2.942, 2.859, 2.8  , 2.798, 2.788, 2.827, 2.795, 2.8  , 3.767]]])

np.c_[msab.c_ttab[0].T].view('datetime64[us]') # SEvt start times (UTC)

array([['2023-04-15T19:18:04.210151', '2023-04-15T19:22:38.213798'],
       ['2023-04-15T19:18:07.296470', '2023-04-15T19:22:41.155448'],
       ...
      ['2023-04-15T19:18:33.170123', '2023-04-15T19:23:04.588021']], dtype='datetime64[us]')

msab.c2tab  # c2sum, c2n, c2per for each event
array([[21.007, 30.807, 31.403, 33.485, 14.285, 42.094, 40.497, 21.573, 31.372, 22.241],
       [31.   , 33.   , 33.   , 32.   , 32.   , 31.   , 30.   , 31.   , 34.   , 31.   ],
       [ 0.678,  0.934,  0.952,  1.046,  0.446,  1.358,  1.35 ,  0.696,  0.923,  0.717]])    # >> 1. would indicate bug

msab.c2tab[0,:].sum()/msab.c2tab[1,:].sum()                              0.9081         # c2per_tot

N=0,1 History Chi2 (target NNVT) : Consistent across 10 evt : Overall Chi2 0.9081

Gun Running N=0,1 Chi2 comparison across 10 evts (opticksMode:2)

cd ~/j/ntds ; ./ntds.sh cfmeta
c_itab
array([[[11274, 11012, 10851, 11422, 11372, 11151, 11535, 11250, 11123, 11001],  N=0  ProcessHits_count
        [13990, 13709, 13649, 13709, 14365, 14609, 13600, 13869, 14686, 14014]], N=1

       [[ 1820,  1724,  1696,  1740,  1806,  1721,  1748,  1745,  1763,  1713],       ProcessHits_true
        [ 1722,  1654,  1702,  1757,  1791,  1761,  1690,  1769,  1763,  1820]],
       ...
       [[ 9214,  9051,  8905,  9427,  9326,  9206,  9531,  9242,  9131,  9055],       NEDEP
        [11990, 11762, 11674, 11658, 12290, 12558, 11621, 11829, 12641, 11923]],

       [[    2,     1,     1,     0,     1,     2,     3,     1,     1,     0],       NBOUND
        [   21,    31,    29,    28,    33,    36,    26,    34,    32,    32]], ##  More sensitive Pyrex with N=1 
       ...
       [[ 1820,  1724,  1696,  1740,  1806,  1721,  1748,  1745,  1763,  1713],       YSAVE
        [ 1722,  1654,  1702,  1757,  1791,  1761,  1690,  1769,  1763,  1820]],

np.diff(c_ttab)/1e6   # seconds between event starts
[[[2.76  2.586 3.106 3.033 3.231 2.62  2.772 2.69  3.112]  ##  Little difference, timings dominated by recording overheads however 
  [3.026 2.71  2.858 2.842 3.094 2.874 2.836 2.841 2.708]]]

np.c_[c_ttab[0,0].view('datetime64[us]'),c_ttab[0,1].view('datetime64[us]')] # event start times (UTC)
[['2023-04-14T13:16:14.802047' '2023-04-14T13:24:48.696600']
 ['2023-04-14T13:16:17.561901' '2023-04-14T13:24:51.722786']
 ...

c2tab  # c2sum, c2n, c2per for each event : Chi2 comparing simulation histories
array([[67.968, 54.336, 57.879, 50.578, 68.72 , 64.923, 60.431, 79.85 , 68.574, 60.395],
       [63.   , 62.   , 63.   , 69.   , 60.   , 63.   , 62.   , 66.   , 64.   , 64.   ],
       [ 1.079,  0.876,  0.919,  0.733,  1.145,  1.031,  0.975,  1.21 ,  1.071,  0.944]]) ##  Consistent histories between N=0, N=1 

c2per_tot:  0.9963   ##  Combined Chi2 over all 10 events : close to 1.

Default Gun Check : again no surprises
N=1 Larger "NBOUND" due to more sensitive Pyrex => more BULK_ABSORB not-at-boundary

B_V1J008_N1_OIPF_NNVT:0:1000_gridxy.png

export OPTICKS_INPUT_PHOTON=GridXY_X1000_Z1000_40k_f8.npy
export OPTICKS_INPUT_PHOTON_FRAME=NNVT:0:1000

MODE=3 EDL=1 N=0 EYE=500,0,2300 CHECK=not_first ~/j/ntds/ntds.sh ana

Photon step points from grid of input photons target NNVT:0:1000 (POM:1)

B_V1J008_N1_OIPF_NNVT:0:1000_POM0.png

POM=0 ntds2_cf ## Check No PMT Optical model

Photon step points from grid of input photons target NNVT:0:1000 (POM:0)

Check With PMT Optical Model disabled : N=0,1 => Consistent Simulations

EVT=001 ~/j/ntds/ntds.sh cf    ## Chi2 history check for 1 evt
c2sum/c2n:c2per(C2CUT)  64.20/69:0.930 (30)

np.c_[siq,_quo,siq,sabo2,sc2,sabo1][:25]  ## A-B history frequency chi2 comparison
[[' 0' 'TO BT BT BT BT SD                    ' ' 0' ' 27052  27168' ' 0.2482' '     0      0']
 [' 1' 'TO SA                                ' ' 1' '  3099   3078' ' 0.0714' '    78     77']
 [' 2' 'TO AB                                ' ' 2' '   948    908' ' 0.8621' '    33      2']
 [' 3' 'TO DR BT SA                          ' ' 3' '   872    929' ' 1.8040' '   132    138']
 [' 4' 'TO BT SD                             ' ' 4' '   889    881' ' 0.0362' '  8241   8241']
 [' 5' 'TO BT BT AB                          ' ' 5' '   343    314' ' 1.2801' '    74    223']
 [' 6' 'TO BT BT SA                          ' ' 6' '   278    274' ' 0.0290' '   479   1032']
 [' 7' 'TO DR BT DR BT SA                    ' ' 7' '   239    226' ' 0.3634' '   154    134']
 [' 8' 'TO BT BT BR BT BT SA                 ' ' 8' '   200    229' ' 1.9604' '   128    176']
 [' 9' 'TO BT SA                             ' ' 9' '   181    177' ' 0.0447' '  8439   8240']
 ['10' 'TO DR BT BT BT BT SD                 ' '10' '   162    176' ' 0.5799' '   162    537']

NEVT=3 ~/j/ntds/ntds.sh cfmeta   ## Compare counts, chi2 for 3 events

array([[[31967, 31903, 32017],  N=0   0:ProcessHits_count        array([[[6.679, 6.614, 6.632],        0:SEvt__TimerDone
        [33305, 33464, 33416]], N=1                                      [6.192, 6.146, 6.107]]])

       [[ 8732,  8693,  8745],        1:ProcessHits_true         msab.c2tab  # c2sum, c2n, c2per for each event
        [ 8696,  8772,  8717]],
                                                                 array([[42.234, 64.198, 51.987],
       [[ 2155,  2149,  2086],        8:NEDEP                           [66.   , 69.   , 67.   ],
        [ 3443,  3577,  3552]],                                         [ 0.64 ,  0.93 ,  0.776]])

       [[  134,   114,   130],        9:NBOUND
        [  134,   131,   129]],

POM:0 simpler histories as photons do not enter PMT
More ProcessHits_false for N=1 (from more Sensitive volume)
Consistent ProcessHits_true
N=0,1 simulation history Chi2 matching

Issue 88, Branch, Merge Request

junosw/-/issues/88

Geant4 FastSim based junoPMTOpticalModel implementation has incorrect polarization and propagation time within PMT and also overcomplex 4 volume PMT geometry

junosw/-/tree/blyth-88-pivot-PMT-optical-model-from-FastSim-to-CustomG4OpBoundaryProcess

Pivot Branch

junosw/-/merge_requests/180

Merge Request

Suggested Further Sanity Checks from MR Reviewer

My checks are mostly between the N=0,N=1 implementations in the branch.

Some sanity checks (eg hit counts, time distribution of hits) between:

current main with the bugs
branch N=0
branch N=1

Good to check both with and without the PMT Optical Model

Aim is to confirm that any differences are small.

Small effect expected as:

proportion of photons effected is small
QE for "backwards" photons gets set to zero
- so need multi-PMT photon history for the bug to change hits

Although expect small effect on hits, needs to be fixed for Opticks full photon validation to match.

Opticks + JUNO : More junoPMTOpticalModel Issues + Validation of Custom4 C4OpBoundaryProcess Fix

Opticks + JUNO : More junoPMTOpticalModel Issues + Validation of Custom4 C4OpBoundaryProcess Fix

Simon C Blyth, IHEP, CAS — 28 April 2023

junoPMTOpticalModel Issues : All from using FastSim

hamaLogicalPMT_fake.png

hamaLogicalPMT_natural.png

FIX : Pivot to Custom Boundary Process For PMT Optical Model ?

C4OpBoundaryProcess::PostStepDoIt : 3-way (A,R,T) Customization

custom4/C4CustomART.h : include into C4OpBoundaryProcess.cc

Custom4 : mini-package with Geant4 customizations

PMT/POM quadrants : choice of PMT geometry and Optical Model impl.

Standalone PMT tests : junosw extracts + Custom4 + Geant4

FewPMT_demo.png

FewPMT_demo0.png

Statistical Comparison of Simulation Histories

Unphysical junoPMTOpticalModel::ModelTrigger FastSim Examples

FewPMT_2563_Unphysical_SD_in_vacuum.png

FewPMT_2912_Unphysical_refract_and_detect.png

FewPMT_2942_Unphysical_cross_geometry_vac_reflect.png

FewPMT_2982_Unphysical_vac_refract.png

FewPMT_5866_vac_reflect.png

FewPMT_7625_Unphysical_detect_mid_vacuum.png

FewPMT_8499_Unphysical_detect_mid_vacuum.png

Debug junoPMTOpticalModel::ModelTrigger with simple struct

FewPMT_NNVT_ModelTriggerBuggy_inner2_Pyrex_trig.png

FewPMT_HAMA_ModelTriggerBuggy_no_inner2_Pyrex_trig.png

junoPMTOpticalModel::ModelTrigger : NNVT : ~2.5% UNPHYSICAL TRIG

Implications of whereAmI bug for NNVT ModelTrigger results

Distance uses DistanceToOut/In depending on G4VSolid::Inside

Implications of whereAmI bug for HAMA ModelTrigger results

junoPMTOpticalModel::ModelTrigger BUG => Unphysical Photons

junoPMTOpticalModel::ModelTrigger FIXED

FewPMT_NNVT_ModelTriggerSimple_avoids_issue.png

FewPMT_check_mtd_trig_vacuum_upper_due_to_obliques.png

N0_one_pmt_FewPMT_hamaLogicalPMT_rain_line_012last.png

N1_one_pmt_FewPMT_hamaLogicalPMT_rain_line_012last.png

NNVT : "rain_line" Example of Statistical History Comparison

HAMA : "rain_line" Example of Statistical History Comparison

Checks with different initial photon pos/dir/shapes

N0_one_pmt_FewPMT_hamaLogicalPMT_lhs_window_line_inwin.png

N1_one_pmt_FewPMT_hamaLogicalPMT_lhs_window_line_inwin.png

N0_one_pmt_FewPMT_hamaLogicalPMT_lhs_reflector_line_reflected.png

N1_one_pmt_FewPMT_hamaLogicalPMT_lhs_reflector_line_reflected.png

Standalone then Insitu Testing

junosw + Opticks : Input Photon Running : ntds bash function

Geant4 + Opticks + NVIDIA OptiX 7 : Hybrid Workflow

Geant4 + Opticks Instrumentation (no GPU usage)

Opticks SEvt loaded into ipython : arrays detailing every photon point

B_V1J008_N1_ip_MOI_Hama:0:1000_yy_target_frame.png

B_V1J008_N1_ip_MOI_Hama:0:1000_yy_frame_close.png

Transform photon history points into target frame

A_V0J008_N0_ip_MOI_Hama:0:1000_a_no_skip_fake.png

A_V0J008_N0_ip_MOI_Hama:0:1000_a.png

B_V1J008_N1_ip_MOI_Hama:0:1000_b.png

100k Input Photons targetting Hama:0:1000 Chi2 History Comparison

Instrument junoSD_PMT_v2::ProcessHits : Compare N=0,1 enum counts

HAMA InputPhoton Chi2 (10 evt) ~/j/ntds/ntds.sh cfmeta

A_V0J008_N0_ip_MOI_NNVT:0:1000_10k.png

B_V1J008_N1_ip_MOI_NNVT:0:1000_10k.png

100k Input Photons targetting NNVT:0:1000 Chi2 History Comparison

NNVT InputPhoton Chi2 (10 evt) ~/j/ntds/ntds.sh cfmeta

Gun Running N=0,1 Chi2 comparison across 10 evts (opticksMode:2)

B_V1J008_N1_OIPF_NNVT:0:1000_gridxy.png

B_V1J008_N1_OIPF_NNVT:0:1000_POM0.png

Check With PMT Optical Model disabled : N=0,1 => Consistent Simulations

Issue 88, Branch, Merge Request

Suggested Further Sanity Checks from MR Reviewer

HAMA InputPhoton Chi2 (10 evt) `~/j/ntds/ntds.sh cfmeta`

NNVT InputPhoton Chi2 (10 evt) `~/j/ntds/ntds.sh cfmeta`