File indexing completed on 2026-04-09 07:48:50
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 """
0022 tbox.py : BoxInBox Opticks vs Geant4 comparisons
0023 =======================================================
0024
0025 Without and with cfg4 runs::
0026
0027 tbox-;tbox-t
0028
0029 Visualize the cg4 created evt in interop mode viewer::
0030
0031 ggv-;ggv-box-test --tcfg4 --load
0032
0033
0034 .. code-block:: py
0035
0036 0.000 100.004 0.129 100.002
0037 1:BoxInBox -1:BoxInBox c2
0038 8d 432190 431510 0.54 [2 ] TO SA
0039 4d 56135 56764 3.50 [2 ] TO AB
0040 86d 10836 10828 0.00 [3 ] TO SC SA
0041 46d 695 723 0.55 [3 ] TO SC AB
0042 866d 140 161 1.47 [4 ] TO SC SC SA
0043 466d 3 10 0.00 [4 ] TO SC SC AB
0044 8666d 0 3 0.00 [5 ] TO SC SC SC SA
0045 4666d 0 1 0.00 [5 ] TO SC SC SC AB
0046 3d 1 0 0.00 [2 ] TO MI
0047 500000 500000 1.21
0048 1:BoxInBox -1:BoxInBox c2
0049 44 488325 488274 0.00 [2 ] MO MO
0050 444 11531 11551 0.02 [3 ] MO MO MO
0051 4444 143 171 2.50 [4 ] MO MO MO MO
0052 44444 0 4 0.00 [5 ] MO MO MO MO MO
0053 4 1 0 0.00 [1 ] MO
0054 500000 500000 0.84
0055
0056
0057 Issues
0058 -------
0059
0060 Refractive index of MO mismatch, or somehow accessing wrong material prop::
0061
0062 INFO:opticks.ana.evt:Evt seqs ['TO SA']
0063 sa(Op)
0064 0 z 300.000 300.000 300.000 r 0.000 100.004 50.002 t 0.098 0.098 0.098 smry m1/m2 4/ 12 MO/Rk 124 (123) 13:TO
0065 1 z -300.000 -300.000 -300.000 r 0.000 100.004 50.002 t 3.070 3.070 3.070 smry m1/m2 4/ 12 MO/Rk 124 (123) 8:SA
0066 sb(G4)
0067 0 z 300.000 300.000 300.000 r 0.129 100.002 50.066 t 0.098 0.098 0.098 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 13:TO
0068 1 z -300.000 -300.000 -300.000 r 0.129 100.002 50.066 t 3.265 3.265 3.265 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 8:SA
0069
0070
0071 Running this in debugger::
0072
0073 ggv-;ggv-box-test --cfg4 --dbg
0074
0075 Shows that steps exiting the world are not well handled by G4OpBoundaryProcess, as
0076 a step status of fWorldBoundary results in cop out with NotAtBoundary status::
0077
0078
0079 183 G4bool isOnBoundary =
0080 184 (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
0081 185
0082 186 if (isOnBoundary) {
0083 187 Material1 = pStep->GetPreStepPoint()->GetMaterial();
0084 188 Material2 = pStep->GetPostStepPoint()->GetMaterial();
0085 189 } else {
0086 190 theStatus = NotAtBoundary;
0087 191 if ( verboseLevel > 0) BoundaryProcessVerbose();
0088 192 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0089 193 }
0090
0091
0092 So add a Pyrex cube to the geometry::
0093
0094 182 local test_config=(
0095 183 mode=BoxInBox
0096 184 analytic=1
0097 185 shape=box
0098 186 boundary=Rock//perfectAbsorbSurface/MineralOil
0099 187 parameters=0,0,0,300
0100 188
0101 189 shape=box
0102 190 boundary=MineralOil///Pyrex
0103 191 parameters=0,0,0,100
0104 192 )
0105
0106
0107
0108 Checking refractive index with pmt-test gives the expected value by comparison with ggv --mat 3::
0109
0110 ggv-;ggv-pmt-test --cfg4 --dbg
0111
0112
0113 (lldb) p thePhotonMomentum*1e6
0114 (double) $6 = 3.2627417774210459
0115
0116 (lldb) p Rindex1
0117 (G4double) $7 = 1.4826403856277466
0118
0119
0120 ::
0121
0122 simon:env blyth$ ggv --mat 3
0123 /Users/blyth/env/bin/ggv.sh dumping cmdline arguments
0124 --mat
0125 3
0126 [2016-03-04 13:14:14.853950] [0x000007fff7057a31] [info] Opticks::preconfigure argv[0] /usr/local/env/optix/ggeo/bin/GMaterialLibTest mode Interop detector dayabay
0127 [2016-Mar-04 13:14:14.858301]:info: GMaterialLib::dump NumMaterials 38
0128 [2016-Mar-04 13:14:14.858686]:info: GPropertyMap<T>:: 3 material m:MineralOil k:refractive_index absorption_length scattering_length reemission_prob MineralOil
0129 domain refractive_index absorption_length scattering_length reemission_prob
0130 60 1.434 11.1 850 0
0131 80 1.434 11.1 850 0
0132 100 1.434 11.1 850 0
0133 120 1.434 11.1 850 0
0134 140 1.64207 11.1 850 0
0135 160 1.75844 11.1 850 0
0136 180 1.50693 11.1 850 0
0137 200 1.59558 10.7949 851.716 0
0138 220 1.57716 10.6971 2201.56 0
0139 240 1.55875 11.5 3551.41 0
0140 260 1.54033 11.3937 4901.25 0
0141 280 1.52192 10.9 6251.1 0
0142 300 1.5035 39.6393 7602.84 0
0143 320 1.49829 117.679 11675 0
0144 340 1.49307 490.025 15747.2 0
0145 360 1.48786 1078.9 19819.4 0
0146 380 1.48264 4941.76 23891.6 0
0147 400 1.47743 11655.2 27963.7 0
0148 420 1.47458 24706.1 36028.8 0
0149 440 1.47251 25254.7 45367.7 0
0150 460 1.47063 24925.3 52039 0
0151 480 1.46876 24277 58710.2 0
0152 500 1.46734 23628.8 68425 0
0153 520 1.4661 22980.5 81100.8 0
0154 540 1.46487 22332.2 93776.7 0
0155 560 1.46369 21277.4 117807 0
0156 580 1.46252 18523.2 152790 0
0157 600 1.46158 14966.4 181999 0
0158 620 1.46081 7061.42 205618 0
0159 640 1.46004 4159.07 229236 0
0160 660 1.45928 5311.87 252855 0
0161 680 1.45851 5615.17 276473 0
0162 700 1.45796 4603.84 300155 0
0163 720 1.45764 3697.27 340165 0
0164 740 1.45733 1365.95 380175 0
0165 760 1.45702 837.71 420184 0
0166 780 1.45671 2274.95 460194 0
0167 800 1.4564 2672.76 500000 0
0168 820 1.4564 1614.62 500000 0
0169
0170
0171 Timing still off despite a quick check suggesting the refractive indices are matching::
0172
0173 INFO:opticks.ana.evt:Evt seqs ['TO BT BT SA']
0174 sa(Op)
0175 0 z 300.000 300.000 300.000 r 0.000 100.004 50.002 t 0.098 0.098 0.098 smry m1/m2 4/ 14 MO/Py -28 ( 27) 13:TO
0176 1 z 99.997 99.997 99.997 r 0.000 100.004 50.002 t 1.086 1.086 1.086 smry m1/m2 14/ 4 Py/MO 28 ( 27) 12:BT
0177 2 z -99.997 -99.997 -99.997 r 0.000 100.004 50.002 t 2.063 2.063 2.063 smry m1/m2 4/ 12 MO/Rk 124 (123) 12:BT
0178 3 z -300.000 -300.000 -300.000 r 0.000 100.004 50.002 t 3.052 3.052 3.052 smry m1/m2 4/ 12 MO/Rk 124 (123) 8:SA
0179 sb(G4)
0180 0 z 300.000 300.000 300.000 r 0.208 100.003 50.105 t 0.098 0.098 0.098 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 13:TO
0181 1 z 99.997 99.997 99.997 r 0.208 100.003 50.105 t 1.154 1.154 1.154 smry m1/m2 14/ 0 Py/?0? 0 ( -1) 12:BT
0182 2 z -99.997 -99.997 -99.997 r 0.208 100.003 50.105 t 2.130 2.130 2.130 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 12:BT
0183 3 z -300.000 -300.000 -300.000 r 0.208 100.003 50.105 t 3.180 3.180 3.180 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 8:SA
0184
0185 Ahha, default time domain of 0:200 ns leads to excess imprecision over 0:5, so set timemax to 10::
0186
0187 INFO:opticks.ana.evt:Evt seqs ['TO BT BT SA']
0188 sa(Op)
0189 0 z 300.000 300.000 300.000 r 0.000 100.004 50.002 t 0.100 0.100 0.100 smry m1/m2 4/ 14 MO/Py -28 ( 27) 13:TO
0190 1 z 99.997 99.997 99.997 r 0.000 100.004 50.002 t 1.089 1.089 1.089 smry m1/m2 14/ 4 Py/MO 28 ( 27) 12:BT
0191 2 z -99.997 -99.997 -99.997 r 0.000 100.004 50.002 t 2.062 2.062 2.062 smry m1/m2 4/ 12 MO/Rk 124 (123) 12:BT
0192 3 z -300.000 -300.000 -300.000 r 0.000 100.004 50.002 t 3.051 3.051 3.051 smry m1/m2 4/ 12 MO/Rk 124 (123) 8:SA
0193 sb(G4)
0194 0 z 300.000 300.000 300.000 r 0.208 100.003 50.105 t 0.100 0.100 0.100 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 13:TO
0195 1 z 99.997 99.997 99.997 r 0.208 100.003 50.105 t 1.155 1.155 1.155 smry m1/m2 14/ 0 Py/?0? 0 ( -1) 12:BT
0196 2 z -99.997 -99.997 -99.997 r 0.208 100.003 50.105 t 2.128 2.128 2.128 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 12:BT
0197 3 z -300.000 -300.000 -300.000 r 0.208 100.003 50.105 t 3.183 3.183 3.183 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 8:SA
0198
0199 In [3]: 1.155 - 1.089
0200 Out[3]: 0.06600000000000006
0201
0202 In [4]: 2.128 - 2.062
0203 Out[4]: 0.06600000000000028
0204
0205 In [5]: 3.183 - 3.051
0206 Out[5]: 0.13199999999999967
0207
0208 In [6]: 0.066*2
0209 Out[6]: 0.132
0210
0211
0212 Chased where times and velocities are set in G4 in g4op- conclusion
0213 is that a sneaky GROUPVEL property is added to the MaterialPropertyTable
0214 for materials with RINDEX. This GROUPVEL is used for optical photon
0215 velocity and resulting times. Subverting this by imposing a GROUPVEL which is
0216 actually the phase velocity (c_light/RINDEX) brings Opticks and CFG4 timings
0217 into alignment.
0218
0219 Running with CPropLib::m_groupvel_kludge = true gets times to agree with Opticks::
0220
0221 INFO:opticks.ana.evt:Evt seqs ['TO BT BT SA']
0222 sa(Op)
0223 0 z 300.000 300.000 300.000 r 0.000 100.004 50.002 t 0.100 0.100 0.100 smry m1/m2 4/ 14 MO/Py -28 ( 27) 13:TO
0224 1 z 99.997 99.997 99.997 r 0.000 100.004 50.002 t 1.089 1.089 1.089 smry m1/m2 14/ 4 Py/MO 28 ( 27) 12:BT
0225 2 z -99.997 -99.997 -99.997 r 0.000 100.004 50.002 t 2.062 2.062 2.062 smry m1/m2 4/ 12 MO/Rk 124 (123) 12:BT
0226 3 z -300.000 -300.000 -300.000 r 0.000 100.004 50.002 t 3.051 3.051 3.051 smry m1/m2 4/ 12 MO/Rk 124 (123) 8:SA
0227 sb(G4)
0228 0 z 300.000 300.000 300.000 r 0.208 100.003 50.105 t 0.100 0.100 0.100 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 13:TO
0229 1 z 99.997 99.997 99.997 r 0.208 100.003 50.105 t 1.089 1.089 1.089 smry m1/m2 14/ 0 Py/?0? 0 ( -1) 12:BT
0230 2 z -99.997 -99.997 -99.997 r 0.208 100.003 50.105 t 2.062 2.062 2.062 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 12:BT
0231 3 z -300.000 -300.000 -300.000 r 0.208 100.003 50.105 t 3.051 3.051 3.051 smry m1/m2 4/ 0 MO/?0? 0 ( -1) 8:SA
0232
0233
0234 TODO: examine the GROUPVEL calc and see how best to do in Opticks ?
0235
0236
0237 """
0238 import os, sys, logging, numpy as np
0239 log = logging.getLogger(__name__)
0240
0241 import matplotlib.pyplot as plt
0242
0243 from opticks.ana.base import opticks_main
0244 from opticks.ana.nbase import vnorm
0245 from opticks.ana.evt import Evt
0246
0247 X,Y,Z,W = 0,1,2,3
0248
0249
0250 if __name__ == '__main__':
0251 args = opticks_main(tag="1",src="torch",det="BoxInBox")
0252
0253 np.set_printoptions(precision=4, linewidth=200)
0254 np.set_printoptions(formatter={'int':hex})
0255
0256 plt.ion()
0257 plt.close()
0258
0259 try:
0260 a = Evt(tag="%s" % args.tag, src=args.src, det=args.det, args=args)
0261 b = Evt(tag="-%s" % args.tag , src=args.src, det=args.det, args=args )
0262 except IOError as err:
0263 log.fatal(err)
0264 sys.exit(args.mrc)
0265
0266 log.info(" A : %s " % a.brief )
0267 log.info(" B : %s " % b.brief )
0268
0269 if 1:
0270 c2max = args.c2max
0271 c2max = None
0272
0273 log.info("\n\nEvt.compare_table.0\n")
0274 Evt.compare_table(a,b, "seqhis_ana seqmat_ana".split(), lmx=20, c2max=c2max, cf=False)
0275 log.info("\n\nEvt.compare_table.1\n")
0276 Evt.compare_table(a,b, "seqhis_ana seqhis_ana_2 seqhis_ana_3 seqhis_ana_4 seqmat_ana".split(), lmx=20, c2max=c2max, cf=True)
0277 log.info("\n\nEvt.compare_table.2\n")
0278 Evt.compare_table(a,b, "pflags_ana hflags_ana".split(), lmx=20, c2max=c2max )
0279
0280
0281
0282 if 0:
0283 if a.valid:
0284 a0 = a.rpost_(0)
0285
0286 a0r = vnorm(a0[:,:2])
0287 print " ".join(map(lambda _:"%6.3f" % _, (a0r.min(),a0r.max())))
0288
0289 if b.valid:
0290 b0 = b.rpost_(0)
0291
0292 b0r = vnorm(b0[:,:2])
0293 print " ".join(map(lambda _:"%6.3f" % _, (b0r.min(),b0r.max())))
0294
0295
0296 if a.valid:
0297 ahm = np.dstack([a.seqhis,a.seqmat])[0]
0298 print "ahm (Op)\n", ahm
0299
0300 if b.valid:
0301 bhm = np.dstack([b.seqhis,b.seqmat])[0]
0302 print "bhm (G4)\n", bhm
0303
0304
0305 if 0:
0306 sel = "TO BT BT SA"
0307 nst = len(sel.split())
0308
0309 sa = Evt(tag="%s" % args.tag, src=args.src, det=args.det, seqs=[sel], args=args)
0310 sb = Evt(tag="-%s" % args.tag, src=args.src, det=args.det, seqs=[sel], args=args )
0311
0312 if sa.valid:
0313 print "sa(Op)"
0314 sa_zrt = sa.zrt_profile(nst)
0315
0316 if sb.valid:
0317 print "sb(G4)"
0318 sb_zrt = sb.zrt_profile(nst)
0319
0320
0321
0322
