Opticks+JUNO : junoPMTOpticalModel FastSim issues and CustomG4OpBoundaryProcess fix
Opticks + JUNO : junoPMTOpticalModel FastSim issues and proposed fix using a CustomG4OpBoundaryProcess
Open source, https://bitbucket.org/simoncblyth/opticks
junoPMTOpticalModel CPU tests reveal issues, all from FastSim use:
- over complex PMT geometry and Geant4 step histories
- inconsistent polarization and propagation time within PMT
Replacing FastSim use with a CustomG4OpBoundaryProcess:
- fixes all issues using a more standard approach
- simpler geometry, simpler simulation history, less code
- would facilitate bringing PMT Optical Model to GPU within Opticks
- as more standard : mostly done already
Simon C Blyth, IHEP, CAS — 20 Dec 2022
junoPMTOpticalModel Issues : All from using FastSim
CustomG4OpBoundaryProcess ?
- 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
inline sboundary::sboundary():
_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
)
;
}
Reflect/Refract Polarization
void junoPMTOpticalModel::Reflect()
{
dir -= 2.*(dir*norm)*norm;
pol -= 2.*(pol*norm)*norm; // similar to G4 TIR pol, NOT Fresnel
}
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, ===>
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
| 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
AOI_BREWSTER_linear_polarization_by_reflection.png
AOI_BREWSTER_reflect_alt_pol.png
Compare Refracted Polarization Impls for Brewster Angle Incidence
- opticks/src/master/sysrap/tests/sboundary_test.sh : build, run, plot
| 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 )
AOI_BREWSTER_brewster_refract.png
AOI_BREWSTER_brewster_refract_alt_pol.png
LayrTest_4_ARTQspx_R12860_4layer.png
LayrTest_4_AQspx_R12860_4layer.png
LayrTest_4_Rspx_R12860_4layer.png
LayrTest_4_Tspx_R12860_4layer.png
LayrTest_2_ARTQspx_R12860_2layer.png
LayrTest_2_AQspx_R12860_2layer.png
LayrTest_2_Rspx_R12860_2layer.png
LayrTest_2_Tspx_R12860_2layer.png
Standalone few PMT test : N=0/1 junoPMTOpticalModel/CustomART.h
Minor changes for standalone use #ifdef PMTFASTSIM_STANDALONE
- JPMT.h parameter access (not PMTSimParamSvc, PMTParamSvc)
- https://github.com/simoncblyth/j/blob/main/PMTFastSim/junoPMTOpticalModel.cc
- https://github.com/simoncblyth/j/blob/main/PMTFastSim/HamamatsuR12860PMTManager.cc
- N=0/1 bool m_natural_geometry = false/true
Optical only Geant4 simulation, with full step point recording
N=0/1 ./U4SimulateTest.sh
Geant4 Simtrace intersection : for 2D geometry plotting
N=0/1 ./U4SimtraceTest.sh
CustomART integrated with InstrumentedG4OpBoundaryProcess
hamaLogicalPMTWrapLV_old.png
hama_UsePMTOpticalModel=0 REVERSE=1 NOGS=1 ~/opticks/extg4/xxv.sh
body hidden under pmt : WHY?
Old PMT Optical Model : G4 Optical Surfaces photocathode, mirror
Old Surface POM
+---------------pmt-Pyrex----------------+ | +-------------body-Pyrex-------------+ | | | | | | | | | | | +------------------------+ | | | | | | | | | | | | | | | | | inner1-Vacuum | |-| | | | | |1e-3 | | | | | | | | +~~coincident~face~~~~~~~+ | | | | | | | | | | | | | | | | | inner2-Vacuum | | | | | | | | | | | | | | | | | +------------------------+ | | | | | | | | | | | +------------------------------------+ | +----------------------------------------+
photons DO NOT enter PMT vacuum
| m_Photocathode_opsurf | inner1<->body | pv<->pv |
| m_mirror_opsurf | inner2<->body | pv<->pv |
Why pmt+body Pyrex (1e-3 mm separation) ? GUESS:
- avoid G4LogicalBorderSurface for every PMT, as 1 body pv
- keeps PMT manager self-contained
Why split vacuum inner1/inner2 ?
- convenient for photocathode and mirror properties
- BUT: fake Vacuum/Vacuum coincident face inside PMT
- no problem : as no photon enter PMT
Alternatively : with CustomG4OpBoundaryProcess
- Two Volume PMT : Pyrex + Vacuum (Instead of Four)
- mirror_opsurf G4LogicalSkinSurface
- Skin not Border: so only need LV
- custom override to Photocathode_opsurf for local_z > 0.
- Effectively : -Z:mirror +Z:photocathode
- NOT YET IMPLEMENTED : SEEMS DOABLE
hamaLogicalPMT_natural.png
hamaLogicalPMT_fake.png
body hidden under inner1 inner2
body is FastSim envelope volume
junoPMTOpticalModel FastSim : Fake Volumes -> Contorted ModelTrigger
Current MultiLayer POM
+---------------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 }
FastSim Big Bouncer : Collect ModelTrigger positions, decisions
SFastSim_Debug.hh
struct SYSRAP_API SFastSim_Debug
{
static std::vector<SFastSim_Debug> record ;
static constexpr const char* NAME = "SFastSim_Debug.npy" ;
static constexpr int LIMIT = 100000 ;
static void Save(const char* dir);
void add();
double posx ; // fs[:,0,:3]
double posy ;
double posz ;
double time ; // fs[:,0,3]
double dirx ; // fs[:,1,:3]
double diry ;
double dirz ;
double dist1 ; // fs[:,1,3]
double polx ; // fs[:,2,:3]
double poly ;
double polz ;
double dist2 ; // fs[:,2,3]
double ModelTrigger; // fs[:,3,0].astype(np.int64)
double whereAmI ; // fs[:,3,1].astype(np.int64)
double c ; // fs[:,3,2]
double PhotonId ; // fs[:,3,3].astype(np.int64)
};
epsilon:tests blyth$ PID=726 ./U4SimulateTest.sh nana
In [1]: extra.shape
Out[1]: (14, 4, 4)
In [2]: extra
Out[2]:
array([[[-112.83 , 0. , 164.918, 0.164],
[ 0.032, 0. , -0.999, 0.001], // dist1:1e-3 mm
[ 0. , -1. , 0. , 165.005],
[ 1. , 1. , 0. , 726. ]], ## ModelTrigger:1
[[-112.83 , 0. , 164.917, 0.164],
[ -0.138, 0. , -0.99 , 166.512], // dist1 and dist2 equal
[ 0. , -1. , 0. , 166.512],
[ 0. , 2. , 0. , 726. ]], ## ModelTrigger:0
[[-156.577, 0. , -148.846, 1.778],
[ 0.81 , 0. , 0.587, 253.614],
[ 0. , 1. , 0. , 0. ],
[ 0. , 1. , 0. , 726. ]],
[[ -95. , 0. , -104.211, 2.166],
[ -0.81 , 0. , 0.587, 177.561],
[ 0. , -1. , 0. , 169.042],
[ 0. , 1. , 0. , 726. ]],
[[-238.764, 0. , -0. , 3.071],
[ -0.81 , 0. , 0.587, 12.404],
[ 0. , -1. , 0. , -1. ],
[ 1. , 2. , 0. , 726. ]],
...
ModelTrigger "pre-triggering" dist1 ahead as DoIt will shift by dist1
- Complication arising from unnatural geometry
- Workaround G4 FastSim needing Region, not boundary
hamaLogicalPMTWrapLV_ModelTriggerYES.png
Starts at upper left, here -->
- FastSim Workaround -> DoIt dist1 shifts -> Effectively "Pre-Triggering" : slide red dots forwards
hamaLogicalPMTWrapLV_ModelTriggerNO.png
- blue dots : ModelTrigger:NO : within FastSim region, but handled by SlowSim
junoPMTOpticalModelSimple : Natural Geometry -> Simple ModelTrigger
Simplified MultiLayer POM
+---------------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
Simplified MultiLayer POM
+---------------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 ?
Custom Boundary POM
+---------------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
- same TMM calc, less code, simpler code
- 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 |
u4/CustomART.h : include into CustomG4OpBoundaryProcess.cc
CustomART<J>::maybe_doIt
template<typename J>
inline char CustomART<J>::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 "JPMT.h"
#include "Layr.h"
#include "U4Touchable.h"
template<typename J> struct CustomART
{
J* parameter_accessor ;
G4double& theTransmittance ;
G4double& theReflectivity ;
G4double& theEfficiency ;
const G4ThreeVector& theGlobalPoint ;
const G4ThreeVector& OldMomentum ;
const G4ThreeVector& OldPolarization ;
const G4ThreeVector& theRecoveredNormal ;
const G4double& thePhotonMomentum ;
CustomART(
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 );
};
hamaLogicalPMTWrapLV_full_history.png
TO BT BT BT BT SR SR BT BR BR BT SR SR SR BT BR BT SR BT SA 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 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 |
Rerunning Single Photons : Very Useful Debug Technique
Save/Load g4state into NP array
Efficiently save/restore thousands of states
MixMaxRng (G4 1042 random engine)
State of engine : 38*uint32
U4Engine::SaveState( NP* states, int idx ); U4Engine::RestoreState( NP* states, int idx );
Rerun "big bouncer" photon
vi U4SimulateTest.sh # SRM_G4STATE_SAVE N=0 ./U4SimulateTest.sh # save g4state into ALL0 vi U4SimulateTest.sh # SRM_G4STATE_RERUN N=1 ./U4SimulateTest.sh # load ALL0, save SEL1
- save g4state before sim each photon
- enables re-simulation of single photons
- compare between geometry versions
- add random "burns" matching consumption
- manual fixup, only for debug
hamaLogicalPMT_two_pmt_geom.png
hamaLogicalPMT_two_pmt_cf.png
junoPMTOpticalModel Issues : All from using FastSim
CustomG4OpBoundaryProcess ?
- 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)
Next Steps
CPU
- Extend JPMT.h to access QE + pmtcat using pmtid
- Extend to NNVT (currently only HAMA)
- Compare absorption, efficiency between implementations
GPU
- Implement CUDA equivalent of CustomART.h
- Integrate wih Opticks qudarap/qsim.h
More standard approach and simpler geometry makes it easier