Opticks Autumn Progress

• Cerenkov energy sampling via ICDF lookups ?
  • workaround poor float/double matching with rejection sampling
  • ICDF lookup working on CPU, good enough chi2 match to rejection
  • TODO: bring over to 2d GPU float texture
• Opticks updates for Geant4 11 API changes, G4 bugs
  • Unclear how many more API changes coming, ongoing work
• 2D planar ray tracing to create cross-section views from intersects
  • Useful technique to investigate geometry issues
• possible CSG bug found : A-(B-C) yielding spurious intersects
  • TODO: confirm with simpler solids (elim. tree balancing)
  • specific case fixable by simpler modelling, general fix unclear

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

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);

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

G4Cerenkov_modified::GetAverageNumberOfPhotons_s2

G4double G4Cerenkov_modified::GetAverageNumberOfPhotons_s2(
    const G4double charge, const G4double beta, const G4Material*, G4MaterialPropertyVector* Rindex) const
{
    G4double BetaInverse = 1./beta;
    G4double s2integral(0.) ;

    for(unsigned i=0 ; i < Rindex->GetVectorLength()-1 ; i++)
    {
        G4double en_0 = Rindex->Energy(i)   ; G4double en_1 = Rindex->Energy(i+1) ;
        G4double ri_0 = (*Rindex)[i]        ; G4double ri_1 = (*Rindex)[i+1] ;
        G4double ct_0 = BetaInverse/ri_0    ; G4double ct_1 = BetaInverse/ri_1 ;
        G4double s2_0 = (1.-ct_0)*(1.+ct_0) ; G4double s2_1 = (1.-ct_1)*(1.+ct_1) ;

        G4bool cross = s2_0*s2_1 < 0. ;
        G4double en_cross =  cross ? en_0 + (BetaInverse - ri_0)*(en_1 - en_0)/(ri_1 - ri_0) : -1. ;
        // linear crossing more precision than s2 zeros

        if( s2_0 <= 0. and s2_1 <= 0. )  // no CK
        {
            // noop
        }
        else if( s2_0 < 0. and s2_1 > 0. )  // s2 becomes +ve within the bin  left edge triangle
        {
            s2integral +=  (en_1 - en_cross)*s2_1*0.5 ;
        }
        else if( s2_0 >= 0. and s2_1 >= 0. )   // s2 +ve across full bin    trapezoid
        {
            s2integral += (en_1 - en_0)*(s2_0 + s2_1)*0.5 ;
        }
        else if( s2_0 > 0. and s2_1 < 0. )  // s2 becomes -ve within the bin   right edge triangle
        {
            s2integral +=  (en_cross - en_0)*s2_0*0.5 ;
        }
    }
    const G4double Rfact = 369.81/(eV * cm);
    return Rfact * charge/eplus * charge/eplus * s2integral ;
}

test_GetAverageNumberOfPhotons_plot

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.0000_en_plot.png

QCKTest_1.0500_en_plot.png

QCKTest_1.1000_en_plot.png

QCKTest_1.1500_en_plot.png

QCKTest_1.2000_en_plot.png

QCKTest_1.2500_en_plot.png

QCKTest_1.3000_en_plot.png

QCKTest_1.3500_en_plot.png

QCKTest_1.4000_en_plot.png

QCKTest_1.4500_en_plot.png

QCKTest_1.5000_en_plot.png

QCKTest_1.5500_en_plot.png

QCKTest_1.6000_en_plot.png

QCKTest_1.6500_en_plot.png

QCKTest_1.7000_en_plot.png

QCKTest_1.7500_en_plot.png

QCKTest_1.7920_en_plot.png

Cerenkov Energy Sample : chi2 compare Rejection vs Lookup

bi	c2sum	ndf	c2p	emn	emx	avp
1.0000	82.0388	99.0000	0.8287	1.5500	15.5000	293.2454
1.0500	91.5019	99.0000	0.9243	1.5500	15.5000	270.4207
1.1000	92.9688	99.0000	0.9391	1.5500	15.5000	246.4823
1.1500	79.5928	99.0000	0.8040	1.5500	15.5000	221.4304
1.2000	100.4855	99.0000	1.0150	1.5500	15.5000	195.2648
1.2500	116.2938	99.0000	1.1747	1.5500	15.5000	167.9857
1.3000	110.0716	99.0000	1.1118	1.5500	15.5000	139.5929
1.3500	84.7895	99.0000	0.8565	1.5500	15.5000	110.0865
1.4000	90.6946	99.0000	0.9161	1.5500	15.5000	79.4666
1.4500	374.1477	99.0000	3.7793	1.5500	15.5000	47.7330
1.5000	147.4192	99.0000	1.4891	3.1136	9.9651	28.7508
1.5500	222.8197	92.0000	2.4220	4.6647	9.5725	16.5258
1.6000	254.6680	78.0000	3.2650	5.7872	9.2549	9.1571
1.6500	156.4206	99.0000	1.5800	7.4768	8.9444	5.3764
1.7000	209.5865	100.0000	2.0959	7.5724	8.6780	2.8470
1.7500	401.4795	101.0000	3.9750	7.6680	8.4288	0.9762
1.7920	54036.1523	101.0000	535.0114	7.7485	7.7895	0.0007

• /tmp/blyth/opticks/QCerenkovIntegralTest/test_makeICDF_SplitBin/QCKTest use_icdf:False

Opticks Updates for Geant4 11.beta (1100)

Aiming to include Opticks example in Geant4 11 Distribution

Geant4 Fermilab Group (Hans Wentzel, Soon Yung Jun) working on this

My role : Opticks updates to cope with API changes

• keep > 500 unittests passing
• frequent snapshots bumping OPTICKS_VERSION_NUMBER (OpticksVersionNumber.hh)
  • https://github.com/simoncblyth/opticks/tags

Geant4 11 API+behavior changes causing problems

• Geant4 11 forces c++17 : Opticks follows (c++14->c++17)
  • requires: om-cleaninstall
• G4PhysicsOrderedFreeVector removed, consolidated to G4PhysicsFreeVector
• G4MaterialPropertiesTable changes caused several Opticks unittest FAILs
  • G4GDMLReadMaterials unable to read GDML with custom properties
  • GetProperty returning bad pointer fMP[-1] when expect nullptr
  • unintended(?) removal of introspection
  • attempt to access a non-existing property/const-property throws fatal exceptions
  • ConstPropertyExists method never able to return false, as fatal exception

From last tarball tried : all tests passing BUT notable slowdown, suspect G4VSolid::CalculateExtent

opticks-t : results with Geant4 1042 and 1100 (geant4.master100721.tar)

1042 : all pass, no slow

1100, Soon 2nd tarball geant4.master100721.tar:

SLOW: tests taking longer that 15 seconds
  3  /45  Test #3  : CFG4Test.CTestDetectorTest   Passed   35.47
  5  /45  Test #5  : CFG4Test.CGDMLDetectorTest   Passed   35.54
  7  /45  Test #7  : CFG4Test.CGeometryTest       Passed   35.34
  27 /45  Test #27 : CFG4Test.CInterpolationTest  Passed   37.76

FAILS:  0   / 501   :  Sat Oct  9 04:57:42 2021
(base) [simon@localhost opticks]$

CTraverser::AncestorTraverse finding extents of 300k solids

• 1042 : 3s, 1100 : 30s
• Suspected culprit G4VSolid::CalculateExtent
• No longer used in mainline Opticks, just slows the 4 tests.

quicktime 2 lAddition_uni_acrylic3

Planar 2d ray tracing to create CSG geometry cross-sections

qudarap/QEvent.cc/QEvent::MakeCenterExtentGensteps

• genstep arrays : (num_gs,6,4) 6*float4 per genstep
• NEW: transform enabled TORCH gensteps
  • 1*int4 : config enums etc.., number of photons to generate
  • 1*float4 : local center position
  • 4*float4 : 4x4 transform

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

• local targetting global 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 }

Planar 2d ray tracing to create CSG geometry cross-sections : code

CSGOptiX/cxs.sh

• control via envvars

CSGOptiX/tests/CSGOptiXSimulateTest.cc

• generates TORCH photons from gensteps
• collects quad4 : 4*union(float4,uint4,int4) at first intersect
• saves array of quad4 intersects with relevant transforms

CSGOptiX/tests/CSGOptiXSimulateTest.py

• loads photons and transforms
• presents intersects using matplotlib (2d) OR pyvista (3d)
• identity info (eg boundary material pairs) used to color intersects

CSGOptiX/cxs.sh : Planar 2d ray tracing, CSG geometry cross-sections

##!/bin/bash -l
cxs=${CXS:-20}   # collect sets of config underneath CXS

if [ "$cxs" == "1" ]; then
    moi=Hama
    ##cegs=16:0:9:1000:18700:0:0:100
    cegs=16:0:9:500
    gridscale=0.05
elif [ "$cxs" == "2" ]; then
    moi=uni_acrylic3
    cegs=16:0:9:100
    ##cegs=0:0:0:1000
    ##cegs=16:4:9:100
    gridscale=0.05
elif [ "$cxs" == "20" ]; then
    note="very tight grid to get into close corners"
    moi=uni_acrylic3
    cegs=16:0:9:100
    gridscale=0.025
fi
export MOI=${MOI:-$moi}
export CEGS=${CEGS:-$cegs}
export GRIDSCALE=${GRIDSCALE:-$gridscale}
export CXS=${CXS:-$cxs}
export TOPLINE="cxs.sh CSGOptiXSimulateTest CXS $CXS MOI $MOI CEGS $CEGS GRIDSCALE $GRIDSCALE ISEL $ISEL ZZ $ZZ"
export BOTLINE="ZOOM $ZOOM LOOK $LOOK"

if [ "$1" == "py" ]; then
    ipython --pdb -i tests/CSGOptiXSimulateTest.py
else
    CSGOptiXSimulateTest
fi

StickMPL_all

• CSG geometry 2D cross-section
• composed of ray-geometry intersection positions
• global intersects transformed into local frame of target

gplt lAddition_uni_acrylic3

• 2D plot created with GDML parser/plotter ana/gplt.py (adhoc "manual" 3d->2d r-z)
• LV : lAddition, Solid : uni_acrylic3, Class : AdditionAcrylicConstruction

• PREV: Suspected : subtracting sph. segment, sagitta 5.68mm, was performance killer.
  • this was causing huge bounding box, now fixed in OptiX 7 geometry
  • TODO: fix in OptiX < 7 geom, and measure performance for both
    • overlarge bbox is expected to cause performance problem, so maybe...

StickMPL_skip_steel

• cylinder with inner removed + cylinder rods are subtracted from cone
• two horizontal green lines just below (0,0) should not be there

• Possible explanation of Opticks "SI BT BT BT BT AB" excess
• Spurious boundaries may causes wrong Water/Acrylic properties used

AdditionAcrylicConstruction::makeAdditionLogical

092 AdditionAcrylicConstruction::makeAdditionLogical(){
108         double ZNodes3[3];
109         double RminNodes3[3];
110         double RmaxNodes3[3];
111         ZNodes3[0] = 5.7*mm; RminNodes3[0] = 0*mm; RmaxNodes3[0] = 450.*mm;
112         ZNodes3[1] = 0.0*mm; RminNodes3[1] = 0*mm; RmaxNodes3[1] = 450.*mm;
113         ZNodes3[2] = -140.0*mm; RminNodes3[2] = 0*mm; RmaxNodes3[2] = 200.*mm;
115         solidAddition_down = new G4Polycone("solidAddition_down",0.0*deg,360.0*deg,3,ZNodes3,RminNodes3,RmaxNodes3);
...
122     solidAddition_up = new G4Sphere("solidAddition_up",0*mm,17820*mm,0.0*deg,360.0*deg,0.0*deg,180.*deg); 
123
124     uni_acrylic1 = new G4SubtractionSolid("uni_acrylic1",solidAddition_down,solidAddition_up,0,G4ThreeVector(0*mm,0*mm,+17820.0*mm));
125
126     solidAddition_up1 = new G4Tubs("solidAddition_up1",120*mm,208*mm,15.2*mm,0.0*deg,360.0*deg);
127     uni_acrylic2 = new G4SubtractionSolid("uni_acrylic2",uni_acrylic1,solidAddition_up1,0,G4ThreeVector(0.*mm,0.*mm,-20*mm));
128     solidAddition_up2 = new G4Tubs("solidAddition_up2",0,14*mm,52.5*mm,0.0*deg,360.0*deg);
130     for(int i=0;i<8;i++)
131     {
132     uni_acrylic3 = new G4SubtractionSolid("uni_acrylic3",uni_acrylic2,solidAddition_up2,0,G4ThreeVector(164.*cos(i*pi/4)*mm,164.*sin(i*pi/4)*mm,-87.5));
133     uni_acrylic2 = uni_acrylic3;
135     }

Start from Polycone (thin cylinder+large cone) and subtract:

1. "sagitta" of huge acrylic sphere
2. cylinder with inner radius (tube/ring shape)
3. 8x cylinder rods

• 1 : subtracting sliver of Acrylic from Acrylic cylinder : no physical effect ?
• 2+3 : pointless subtracting holes for daughter volumes ?

JUST the Polycone should work fine : and give same results

CSG sub sub bug ? Subtracted subtraction yielding spurious intersects

solidAddition_down (GPolycone) Z: 5.7, 0.0, -140.0

solidAddition_up (G4Sphere)

uni_acrylic1 (G4SubtractionSolid)

solidAddition_down - solidAddition_up [ Z : +17820.0 ] subtract sagitta cap

solidAddition_up1 (G4Tubs) rmin/rmax 120/208 hz 15.2

uni_acrylic2 (G4SubtractionSolid)

uni_acrylic1 - solidAddition_up1 [ Z : -20 ]

• [15.2, -15.2] -20. = [-4.8, -35.2]
• [15.2, -15.2]*1.01 - 20.0 = [ -4.648, -35.352] 1% expand subtracted

solidAddition_up2 (G4Tubs) rmin/rmax 0/14 hz 52.5

8 x uni_acrylic3 (G4SubtractionSolid)

uni_acrylic2 - solidAddition_up2 [ Z : -87.5 ]

• [ 52.5, -52.5] - 87.5 = [ -35., -140.]

20_figs_positions_pvplt_0,1,2,3_csgsubsubbug.png

• .
• close look at intersect Z positions : see jump --->
• show enlarged C, Z positions   A - (B - C)

StickPV_rod

Mind the gap:

• is it important ?

Geant4 Geometry : is the cavity between acrylic and screw important ?

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

NB NO CSG subtraction for the water cylinders or steel rods

If gap between acrylic cone and steel rods is not important

Model with hierarchy of 2 volumes:

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

NB NO CSG subtraction for the steel rods

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

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

Multiple non-interfering CSG subtractions work fine

• problem arises from >1 subtractions "on top of each other",
  • colocation prevents coincident faces avoidance trick of "nudge expanding" the subtracted solid
  • TODO: try without the nudge expansion : expect will just shift the spurious intersects

CSG Subtraction of a "pipe" cylinder (with inner rmin>0) : gives spurious inner intersects

• TODO: test with simpler solids, to eliminate tree balancing
• problem due to "pipe" being implemented as CSG subtraction
• could avoid by implementing pipe directly as primitive
  • significant work, but doable, not a general fix

Problem Arises from inherent CSG fragility regarding coincident faces

• general solution unknown
  • need to seek inspiration from other peoples CSG implementations
  • not aware of any other GPU implementations however

CSG Deep Tree : Positivize tree using De Morgans laws

1st step to allow balancing : Positivize : remove CSG difference di operators

                                                     ...    ...

                                               un          cy

                                       un          cy

                               un          cy

                       un          cy

               un          cy

       di          cy

   cy      cy

                                                     ...    ...

                                               un          cy

                                       un          cy

                               un          cy

                       un          cy

               un          cy

       in          cy

   cy      !cy

CSG Examples

LHCb-RICH (Yunlong Li, Manchester University) : Opticks OpenGL Render (OKTest)