g4cx/tests/G4CXTest_raindrop_simtrace.py

0001 #!/usr/bin/env python
0002 """
0003 G4CXTest_raindrop_simtrace.py
0004 ===============================
0005
0006 Simplify CSGOptiX/cxt_min.py
0007
0008
0009 """
0010 import os, re, logging, numpy as np
0011 log = logging.getLogger(__name__)
0012
0013 from opticks.sysrap.sevt import SEvt
0014 from opticks.ana.p import cf    # load CSGFoundry geometry info
0015 from opticks.CSG.CSGFoundry import KeyNameConfig
0016
0017 GLOBAL = 1 == int(os.environ.get("GLOBAL","0"))
0018 GSGRID = 1 == int(os.environ.get("GSGRID","0"))
0019 FRAME = 1 == int(os.environ.get("FRAME","0"))
0020 MODE = int(os.environ.get("MODE","2"))
0021
0022 X,Y,Z = 0,1,2
0023 H,V = X,Z
0024
0025
0026 if MODE in [2,3]:
0027     import opticks.ana.pvplt as pvp
0028 pass
0029
0030 if __name__ == '__main__':
0031     print("GLOBAL:%d MODE:%d" % (GLOBAL,MODE))
0032
0033     t = SEvt.Load("$TFOLD", symbol="t")
0034     print(repr(t))
0035
0036     label = t.f.base
0037     sf = t.f.sframe
0038     w2m = sf.w2m
0039     gs = t.f.genstep
0040     st = t.f.simtrace
0041
0042     st_t = st[:,0,3]
0043     st_gp_bn = st[:,2,3].view(np.int32)
0044     st_gp = st_gp_bn >> 16  # globalPrimIdx
0045     st_bn = st_gp_bn & 0xffff   ## boundary
0046
0047
0048     gs_pos = gs[:,1]
0049     all_one = np.all( gs_pos[:,3] == 1. )
0050     all_zero = np.all( gs_pos[:,:3]  == 0 )
0051     assert all_one  # SHOULD ALWAYS BE 1.
0052     #assert all_zero # NOT ZERO WHEN USING CE_OFFSET=CE CE_SCALE=1
0053     gs_tra = gs[:,2:]
0054     assert gs_tra.shape[1:] == (4,4)
0055
0056     ggrid = np.zeros( (len(gs), 4 ), dtype=np.float32 )
0057     for i in range(len(ggrid)): ggrid[i] = np.dot(gs_pos[i], gs_tra[i])
0058     ## np.matmul/np.tensordot/np.einsum can probably do this without the loop
0059     lgrid = np.dot( ggrid, w2m )
0060     ugrid = ggrid if GLOBAL else lgrid
0061
0062
0063
0064     #presel = slice(None)
0065     #presel = st_t > 0.
0066     presel = st_gp > -1
0067     #presel = np.logical_and( st_gp > -1, st_t > 0. )
0068
0069
0070
0071     ust = st[presel]
0072
0073     gp_bn = ust[:,2,3].view(np.int32)    ## simtrace intersect boundary indices
0074     gp = gp_bn >> 16      ## globalPrimIdx
0075     bn = gp_bn & 0xffff   ## boundary
0076
0077
0078     u_gp, n_gp = np.unique(gp, return_counts=True )
0079     o_gp = np.argsort(n_gp)[::-1]
0080     gp_tab = np.c_[u_gp,n_gp][o_gp]
0081     print(repr(gp_tab))
0082     u_gp = gp_tab[:,0]
0083     n_gp = gp_tab[:,1]
0084
0085
0086     ## simtrace layout see sysrap/sevent.h
0087     gpos = ust[:,1].copy()
0088     gpos[...,3] = 1.   ## transform as position
0089     lpos = np.dot( gpos, w2m )
0090     upos = gpos if GLOBAL else lpos
0091
0092
0093
0094     if MODE == 3:
0095         pl = pvp.pvplt_plotter(label)
0096         pvp.pvplt_viewpoint(pl)   # sensitive to EYE, LOOK, UP envvars
0097         if FRAME: pvp.pvplt_frame(pl, sf, local=not GLOBAL, pickable=False )
0098
0099
0100         colors = ["red", "green", "blue", "yellow", "magenta", "grey", "white" ]
0101
0102         for i, g in enumerate(u_gp):
0103             w = np.where( gp == g )[0]
0104             color = colors[i % len(colors)]
0105
0106             prn = cf.primname[g] if g > -1 and g < len(cf.primname) else "???"
0107
0108             label = " gp:%d len(w):%d prn:%s col:%s " % ( g, len(w), prn, color )
0109             pvp.pvplt_add_points(pl, upos[w,:3], color=color, label=label )
0110         pass
0111
0112         if GSGRID: pl.add_points(ugrid[:,:3], color="r", pickable=False )
0113         cp = pvp.pvplt_show(pl, incpoi=-5, legend=True, title=None )
0114     pass
0115
0116
0117