u4/tests/U4SimulateTest_ph.py

0001 #!/usr/bin/env python
0002 """
0003 U4SimulateTest_ph.py
0004 ========================
0005
0006 ::
0007
0008     u4t
0009     ./U4SimulateTest.sh ph
0010     ./U4SimulateTest.sh nph
0011
0012 """
0013 import os, numpy as np
0014 from opticks.ana.fold import Fold
0015 from opticks.ana.p import *
0016
0017 LABEL = os.environ.get("LABEL", "U4SimulateTest_ph.py" )
0018 N =  int(os.environ.get("VERSION", "-1"))
0019 VERSION = N
0020 CMDLINE = "N=%d ./U4SimulateTest.sh ph" % N
0021
0022 MODE =  int(os.environ.get("MODE", "2"))
0023 assert MODE in [0,2,3]
0024 PID = int(os.environ.get("PID", -1))
0025 if PID == -1: PID = int(os.environ.get("OPTICKS_G4STATE_RERUN", -1))
0026
0027 if MODE > 0:
0028     from opticks.ana.pvplt import *
0029 pass
0030
0031 if __name__ == '__main__':
0032
0033     t = Fold.Load(symbol="t")
0034     tp = Fold.Load(symbol="tp", parent=True)
0035
0036     if not tp is None:
0037         rr = np.array([tp.run_meta.T_BeginOfRun, tp.run_meta.T_EndOfRun], dtype=np.uint64 )
0038         rre_ = "np.c_[rr.view(\"datetime64[us]\")]"
0039         print(rre_)
0040         print(eval(rre_))
0041     pass
0042
0043     print(repr(t))
0044     print( "MODE:%d" % (MODE) )
0045     print( "N:%d" % (N) )
0046
0047     if 'SPECS' in os.environ:
0048         SPECS = np.array(t.U4R_names.lines)
0049         st_ = t.aux[:,:,2,3].view(np.int32)
0050         st = SPECS[st_]   # get out of bounds, when using st_ for other checks
0051     else:
0052         SPEC, st_, st = None, None, None
0053     pass
0054
0055
0056     axes = 0, 2  # X,Z
0057     H,V = axes
0058     label = LABEL
0059
0060     #pos = t.photon[:,0,:3]
0061     pos = t.record[:,:,0,:3].reshape(-1,3)   # xyz : all photons, all steps
0062     q_ = t.seq[:,0]    #  t.seq shape eg (1000, 2, 2)
0063     q = ht.seqhis(q_)    # history label eg b'TO BT BT SA ... lots of blankspace...'
0064     qq = ht.Convert(q_)  # (n,32) int8 : for easy access to nibbles
0065     n = np.sum( seqnib_(q_), axis=1 )
0066
0067
0068     ReplicaNumber=False
0069     if ReplicaNumber:
0070         ## ReplicaNumber : but when not more than one of each type of volume this is -1
0071         rp = t.record[...,1,3].view(np.int32)
0072         np.set_printoptions(edgeitems=50)
0073
0074         u_rp, i_rp, v_rp, n_rp = np.unique(rp, axis=0, return_index=True, return_inverse=True, return_counts=True )
0075         print("\nnp.c_[np.arange(len(i_rp)),i_rp,n_rp,u_rp] ## unique ReplicaNumber sequences ")
0076         print(np.c_[np.arange(len(i_rp)),i_rp,n_rp,u_rp])
0077         print("\nlen(v_rp) : %d ## v_rp : unique array indices that reproduce original array  " % len(v_rp))
0078         assert len(rp) == len(v_rp)
0079     pass
0080
0081     ## uniqing the "|S96" label because NumPy lacks uint128
0082     qu, qi, qn = np.unique(q, return_index=True, return_counts=True)
0083     quo = np.argsort(qn)[::-1]
0084     expr_ = "np.c_[qn,qi,qu][quo]"
0085     expr = eval(expr_)
0086
0087
0088     print("\n%s  ## unique histories qu in descending count qn order, qi first index " % expr_ )
0089     print(expr)
0090
0091     numpystr_ = lambda _:"%5s" % _.strip().decode("utf-8")
0092     int_ = lambda _:"%5d" % _
0093     formatter = {'numpystr':numpystr_, 'int':int_}
0094     ## HMM: when concatenating all the types get converted to numpystr, so not much control
0095     # idea for more control of array presentation
0096     print("\n",np.array2string(expr,formatter=formatter))
0097
0098
0099
0100     ws_ = 0
0101     ws = np.where( q[:,0] == qu[quo][ws_] )   # select photons with the ws_ most prolific history
0102
0103
0104     #print("\nq[v_rp == 0]  ## history flag sequence for unique ReplicaNumber sequence 0"  )
0105     #print(repr(q[v_rp == 0]))
0106
0107     print("\nnp.unique(n, return_counts=True) ## occupied nibbles  ")
0108     print(repr(np.unique(n, return_counts=True)))
0109
0110     print("\nq[n > 16]  ## flag sequence of big bouncers  ")
0111     print(repr(q[n>16]))
0112
0113     n_cut = 10
0114
0115     expr = "q[n > %d]" % n_cut
0116     print("\n%s  ## flag sequence of big bouncers  " % expr )
0117     print(repr(eval(expr)))
0118
0119     expr = "np.c_[n,q][n>%d]" % (n_cut)
0120     print("\n%s  ## nibble count with flag sequence of big bouncers  " % expr )
0121     print(eval(expr))
0122
0123     print("\nnp.where(n > 28)  ## find index of big bouncer " )
0124     print(np.where(n > 28))
0125
0126     expr = "np.c_[np.where(n>%d)[0],q[n > %d]]" % (n_cut,n_cut)
0127     print("\n%s  ## show indices of multiple big bouncers together with history " % expr)
0128     print(eval(expr))
0129
0130     expr = " t.record[%d,:n[%d],0] " % (PID,PID)
0131     print("\n%s  ## show step record points for PID %d  " % (expr, PID))
0132     print(eval(expr))
0133
0134     expr = " np.where( t.record[:500,:,0,2] < 0 ) "
0135     print("\n%s ## look for records with -Z positions in the first half, that all start in +Z " % expr )
0136     print(eval(expr))
0137
0138     expr = " np.where( t.record[500:,:,0,2] > 0 ) "
0139     print("\n%s ## look for records with +Z positions in the second half, that all start in -Z " % expr )
0140     print(eval(expr))
0141
0142     print("\nt.base : %s  VERSION: %d " % (t.base, VERSION))
0143
0144
0145     pos = t.photon[:,0,:3]
0146     flagmask = t.photon[:,3,3].view(np.int32)
0147
0148     expr = " np.where( np.logical_and( pos[:,0] < -150, pos[:,2] > 200 )  ) "
0149     print("\n%s ## common break out line, use t.photon end position to get indices " % expr )
0150     print(eval(expr))
0151
0152     expr = " q[np.where( np.logical_and( pos[:,0] < -150, pos[:,2] > 200 ))][:3] "
0153     print("\n%s ## common break out line, use t.photon end position to get indices wline and see what the q histories are " % expr )
0154     print(eval(expr))
0155
0156
0157     SURFACE_DETECT = 0x1 << 6
0158     sd = flagmask & SURFACE_DETECT
0159     w_sd = np.where( sd )[0]
0160
0161
0162     x_midline = np.logical_and( pos[:,0] > -251, pos[:,0] < -249 )
0163     z_midline = np.logical_and( pos[:,2] > -250, pos[:,2] <  250 )
0164     xz_midline = np.logical_and( x_midline, z_midline )
0165     w_midline = np.where(xz_midline)[0]
0166
0167
0168     yy = t.record[:,:,0,1]
0169     myy = np.max( np.abs(yy), axis=1 )     ## max absolute y of all step points for each photon
0170
0171     wyy0 = np.where( myy == 0. )[0]        ## photon indices that always stay in the plane
0172     wyy1 = np.where( myy > 0 )[0]           ## photon indices with deviations out of plane
0173     sd0 = np.logical_and( sd, myy == 0. )   ## mask SD photons staying in plane
0174
0175
0176     #ppos0_ = "None"
0177     #ppos0_ = "pos #   "
0178     ppos0_ = "t.record[:,0,0,:3] # 0-position   "
0179
0180     #ppos1_ = "None"
0181     ppos1_ = "t.record[:,1,0,:3] # 1-position   "
0182
0183     #ppos2_ = "None"
0184     ppos2_ = "t.record[:,2,0,:3] # 2-position   "
0185
0186
0187     #ppos3_ = "None"
0188     ppos3_ = "t.photon[:,0,:3]"
0189
0190     ppos0  = eval(ppos0_)
0191     ppos1  = eval(ppos1_)
0192     ppos2  = eval(ppos2_)
0193     ppos3  = eval(ppos3_)
0194
0195     ppos = t.record[:,:,0,:3].reshape(-1,3)
0196
0197
0198
0199     elem = []
0200     elem.append(CMDLINE)
0201     if not ppos0 is None: elem.append("b:%s" % ppos0_)
0202     if not ppos1 is None: elem.append("r:%s" % ppos1_)
0203     if not ppos2 is None: elem.append("g:%s" % ppos2_)
0204     if not ppos3 is None: elem.append("c:%s" % ppos3_)
0205     label = "\n".join(elem)
0206
0207
0208     if MODE == 0:
0209         print("not plotting as MODE 0  in environ")
0210     elif MODE == 2:
0211         fig, axs = mpplt_plotter(label=label)
0212         assert len(axs) == 1
0213         ax = axs[0]
0214
0215         ax.set_ylim(-250,250)
0216         ax.set_xlim(-500,500)
0217
0218         if not ppos0 is None: ax.scatter( ppos0[:,H], ppos0[:,V], s=1, c="b" )
0219         if not ppos1 is None: ax.scatter( ppos1[:,H], ppos1[:,V], s=1, c="r" )
0220         if not ppos2 is None: ax.scatter( ppos2[:,H], ppos2[:,V], s=1, c="g" )
0221         if not ppos3 is None: ax.scatter( ppos3[:,H], ppos3[:,V], s=1, c="c" )
0222
0223         fig.show()
0224     elif MODE == 3:
0225         pl = pvplt_plotter(label)
0226         os.environ["EYE"] = "0,100,165"
0227         os.environ["LOOK"] = "0,0,165"
0228         pvplt_viewpoint(pl)
0229         pl.add_points(ppos )
0230         pl.show()
0231     pass
0232 pass