Opticks replaces Geant4 optical photon simulation with an equivalent implementation that benefits from state-of-the-art GPU ray tracing from NVIDIA OptiX, which can yield performance >1000x Geant4.
All optically relevant aspects of Geant4 context are translated+copied to GPU:
- geometry : solids, structure, material+surface properties
- generation : Cerenkov + Scintillation (using Gensteps from Geant4)
- propagation : Rayleigh scattering, absorption, reemission, boundary, +multi-layer TMM (WIP)
Tests Reveal : More Opticks+JUNO Geom. Issues ALL KNOWN, FIXED
- WIP : Bringing multi-layer TMM "FastSim" into Opticks
JUNO ISSUES
NNVT : MaskTail impinges MaskVirtual
HAMA : BodySolid impinges MaskTail
BirksConstant1 : 1,000,000x TOO BIG
OPTICKS PRIM ISSUES
nmskSolidMask : ellipsoid hole at "apex" issue
nmskSolidMaskTail : very thin cylinder "lip" issues
nmskSolidMaskVirtual : cone precision + near-apex issues
rare close to ellipsoid "apex" rays missed : fixed with zcut+safety
without uncoincidence the subtraction has coincident edge
much better after uncoinciding, but upper lip still has small issue
some spurious intersects in center and at the edge
illuminate problem areas with more gensteps, issues all in plane Z=-39 mm
extg4/x4t.sh : Geant4 reference intersects
extg4/cf_x4t.sh : blue:nmskTailOuter orange:nmskTailInner
CSG/cf_ct.sh : Opticks CSG running on CPU
Missing thin cylinder edge intersects
CAUSED : by mis-translation of very thin cylinder as "disc" not "cylinder"
CSG/ct.sh : PMT mask "lip" reveals issue with thin cylinder (hz 0.15 mm)
CSG/ct.sh : hz 0.15 mm thin cylinder precision loss issue
cd ~/opticks/CSG ; UNEXPECTED=1 NOLEGEND=1 ./ct.sh ana
near axial rays suffer precision loss -> spurious
Loss of precision issue : only apparent with very thin cylinders
CSG_OLDCYLINDER : CSG/csg_intersect_leaf_oldcylinder.h::intersect_leaf_oldcylinder
CSG_CYLINDER : CSG/csg_intersect_leaf.h::intersect_leaf_cylinder
CSG/tests/CSGIntersectComparisonTest.sh
Comparison of cylinder implementations
FOCUS=-257,-39,7 ./ct.sh ana ## using CSG/tests/CSGSimtraceTest.{cc/py}
new cylinder impl : avoids spurious + improves precision + simpler
Even with uncoincidence : still left with sprinkle between cylinder and cone
Uncoincidence : expanded cylinder up, in this case changing geometry
Without Uncoincidence : clear spurious lines cyl-cyl-cone
c2 t^2 + 2 c1 t + c0 = 0
disc = c1*c1 - c0*c2, for small c0 or c2 get precision loss in one root
Avoid precision loss subtraction
#ifdef NAIVE_QUADRATIC
t1 = (-b - sdisc)/d ;
t2 = (-b + sdisc)/d ; // t2 >= t1
#else
// pick root depending on sign of b
float q = b > 0.f ? -(b + sdisc) : -(b - sdisc) ;
float root1 = q/d ;
float root2 = c/q ;
t1 = fminf( root1, root2 );
t2 = fmaxf( root1, root2 );
#endif
Solution : pick t OR 1/t quadratic roots depending on coefficients : avoiding precision loss
Normal Quadratic Alternative quadratic in 1/t
----------------- --------------------------------
d t^2 + 2b t + c = 0 c (1/t)^2 + 2b (1/t) + d = 0
-2b +- sqrt((2b)^2 - 4dc ) -2b +- sqrt( (2b)^2 - 4dc )
t = -------------------------- 1/t = ----------------------------
2d 2c
-b +- sqrt( b^2 - d c ) 1/t = -b +- sqrt( b^2 - d c )
t = ----------------------- -------------------------
d c
c
t = ---------------------------
-b +- sqrt( b^2 - d c )
c # cd ~/opticks/CSG
./nmskSolidMaskVirtual.sh withnudge # GeoChain translate
./nmskSolidMaskVirtual.sh ct # intersect using ct.sh, CSGSimtraceTest.cc
./nmskSolidMaskVirtual.sh ana # plot with CSGSimtraceTest.py
./nmskSolidMaskVirtual.sh unx # select unexpected, save to array
./nmskSolidMaskVirtual.sh sample # rerun intersect for selected, using CSGSimtraceSample.cc
//intersect_leaf_oldcone t_near_alt 477.8 t_far_alt 1.225e+09 t_near_alt-t_near 17 t_far_alt-t_far 0
//intersect_leaf_oldcone r1 264.0500 z1 97.0000 r2 132.0250 z2 194.0500 : z0 291.1000
//intersect_leaf_oldcone c2 -0.0000 c1 365.0782 c0 -348871.4688 disc 133281.9219 disc > 0.f 1 : tth -1.3604
//intersect_leaf_oldcone c0 -3.489e+05 c1 365.1 c2 -5.96e-07 t_near 460.8 t_far 1.225e+09 sdisc 365.0780 (-c1-sdisc) -730.2 (-c1+sdisc) -0.0002747
//intersect_leaf_oldcone t_near_alt 477.8 t_far_alt 1.225e+09 t_near_alt-t_near 17 t_far_alt-t_far 0
- HIT
q0 norm t ( -0.5457 0.0000 0.8380 460.8000)
q1 ipos sd ( -212.8504 0.0000 114.4940 0.0000)- sd < SD_CUT : -0.0010
q2 ray_ori t_min ( 158.4300 0.0000 -158.4300)
q3 ray_dir gsid ( -0.8057 0.0000 0.5923 C4U ( 0 0 0 255 ) )
o.x 158.4300 v.x -0.8057 t0(-o.x/v.x) 196.6291 z0 -41.9699
2022-09-17 16:05:30.725 INFO [9037364] [CSGSimtraceSample::intersect@89] CSGSimtraceSample::desc
18 rays with unexpected cyl-cone intersects : all extend close to cone apex
CSG_OLDCONE
makes incorrect assumption that all rays intersect infinite cone
NOT TRUE FOR RAYS CLOSE TO APEX
https://bitbucket.org/simoncblyth/opticks/src/master/CSG/csg_intersect_leaf_oldcone.h
CSG_CONE
cap intersects indep of cone intersects + adopt robust_quadratic_roots
https://bitbucket.org/simoncblyth/opticks/src/master/CSG/csg_intersect_leaf_newcone.h
New CSG_CONE avoids apex and precision loss issues in nmskSolidMaskVirtual
NNVT : TWO overlap issues visible, one was fixed by updating to latest junosw, see next page
NNVT : ONE overlap issue visible, MaskTail impinges MaskVirtual
(using ct.sh : Opticks CSG on CPU)
NNVT : x4t.sh Geant4 spurious intersects visible (cf prior)
(using x4t.sh : Geant4 intersects)
HAMA : ONE overlap issue, BodySolid impinges MaskTail
(using ct.sh : Opticks CSG on CPU)
HAMA : following fix of bug 33 using ellipsoid shift
(using ct.sh : Opticks CSG on CPU)
HAMA : BodySolid impinges MaskTail (mct.sh)
HAMA : after bug 33 fix
NNVT : MaskTail impinges MaskVirtual (mct.sh)
J004_Hama:0:1000 g4cx/gxt.sh
J004_Hama:0:1000 g4cx/gxt.sh
J004_NNVT:0:1000 g4cx/gxt.sh
JPMT.h : header-only prop. access
Added property file (and directory tree) loading with:
NPFold.h and NP.hh are header-only impls
export NP_PROP_BASE=$JUNOTOP/data/Simulation/DetSim
TMM "ART" calc : complex math
Access to thickness and energy dependent : rindex, kindex
CPU : using NP interpolation : STATUS : works fine
GPU : using qprop.h interpolation : STATUS : works in isolation
"FastSim" integration of ART calc ? : STATUS : Initial development of CPU test machinery
Standalone test of single PMT with junoPMTOpticalModel : integrating Opticks photon history recording
junosw/Simulation/SimSvc/MultiFilmSimSvc
epsilon:MultiFilmSimSvc blyth$ find . -type f ./CMakeLists.txt ./python/MultiFilmSimSvc/__init__.py ./MultiFilmSimSvc/MultiFilmModel.h ./src/OpticalSystem.h ./src/Layer.h ./src/Material.h ./src/Layer.cc ## Layer, ThickLayer, ThinLayer ./src/Matrix.h ./src/OpticalSystem.cc ./src/MultiFilmModel.cc ./src/Material.cc ./src/Matrix.cc
template<typename T> struct Layr
{
T d ; // thickness : zero means incoherent
T pad ;
#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);
// all calculations in ctor
};
https://github.com/simoncblyth/j/blob/main/Layr/Layr.h Implemented simply : compiles with nvcc/gcc for GPU/CPU
Complex math on GPU
template<typename T, int N>
Stack<T,N>::Stack(
T wl,
T minus_cos_theta,
const StackSpec<T>& ss )
{
// minus_cos_theta, aka dot(mom,normal)
#ifdef WITH_THRUST
using thrust::complex ;
using thrust::norm ;
using thrust::conj ;
using thrust::exp ;
using thrust::sqrt ;
using thrust::sin ;
using thrust::cos ;
#else
using std::complex ;
using std::norm ;
using std::conj ;
using std::exp ;
using std::sqrt ;
using std::sin ;
using std::cos ;
#endif
// same complex TMM math works on both GPU and CPU
}
LayrMinimal output
$ name=LayrMinimal $ gcc $name.cc -std=c++11 -lstdc++ -o /tmp/$name $ /tmp/$name StackSpecL0 ( 1.4820 0.0000 ; 0.0000) L1 ( 1.9200 0.0000 ; 36.4900) L2 ( 2.4290 1.3660 ; 21.1300) L3 ( 1.0000 0.0000 ; 0.0000) ... comp Layr n:( 0.0000 0.0000)s d: 0.0000 st:( 0.0000 0.0000)s ct:( 0.0000 0.0000)s rs:( -0.0047 -0.2114)s rp:( 0.0047 0.2114)s ts:( 0.0067 0.6972)s tp:( 0.0067 0.6972)s S | ( 0.0138 -1.4342)s ( 0.1681 -0.8562)s | | ( -0.3032 0.0038)s ( -0.1772 0.4389)s | P | ( 0.0138 -1.4342)s ( -0.1681 0.8562)s | | ( 0.3032 -0.0038)s ( -0.1772 0.4389)s | ART R_s 0.0447 R_p 0.0447 T_s 0.3280 T_p 0.3280 A_s 0.6273 A_p 0.6273 R 0.0447 T 0.3280 A 0.6273 A_R_T 1.0000 wl 440.0000 mct -1.0000
#include "Layr.h"
//typedef double T ;
typedef float T ;
int main(int argc, char** argv)
{
T mct = argc > 1 ? std::atof(argv[1]) : -1.f ;
// minus_cos_theta
T wl = argc > 2 ? std::atof(argv[2]) : 440.f ;
// wavelength_nm
StackSpec<T> spec(StackSpec<T>::EGet());
// sensitive to L0, L1, L2, L3 envvars
Stack<T,4> stack(wl, mct, spec );
// ART calc done in ctor
std::cout << spec << stack ;
return 0 ;
}
Impl: arrays of simple struct
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)
All LayrTest use CPU interpolated rindex, uploaded to device (OK for testing, not for production)
https://bitbucket.org/simoncblyth/opticks/src/master/qudarap/QPMT.hh
https://bitbucket.org/simoncblyth/opticks/src/master/qudarap/tests/QPMTTest.sh
U4PMTFastSimTest.cc
21 U4PMTFastSimTest::U4PMTFastSimTest()
22 :
23 phy((G4VUserPhysicsList*)new U4Physics),
24 run(InitRunManager(phy)),
25 rec(new U4RecorderTest(run))
26 {
27 run->BeamOn(U::GetEnvInt("BeamOn",1));
28 }
29
30 int main(int argc, char** argv)
31 {
32 OPTICKS_LOG(argc, argv);
33
34 SEvt evt ;
35 SEvt::AddTorchGenstep();
36
37 U4PMTFastSimTest t ;
38
39 return 0 ;
40 }
j/PMTFastSim
u4/tests/U4PMTFastSimTest.sh
TODO: