u4/tests/U4SimulateTest_mt.py

0001 #!/usr/bin/env python
0002 """
0003 U4SimulateTest_mt.py
0004 ========================
0005
0006 ::
0007
0008     u4t
0009     ./U4SimulateTest.sh mt
0010
0011 """
0012 import os, textwrap, numpy as np
0013 from opticks.ana.fold import Fold, AttrBase
0014 from opticks.ana.p import *
0015
0016 # TODO: these should be coming from standard place, not duplicate
0017 BULK_ABSORB = 0x1 <<  3
0018 SURFACE_DETECT = 0x1 << 6
0019 SURFACE_ABSORB = 0x1 << 7
0020 SURFACE_DREFLECT  = 0x1 << 8
0021 SURFACE_SREFLECT = 0x1 <<  9
0022 BOUNDARY_REFLECT  = 0x1 << 10
0023 BOUNDARY_TRANSMIT = 0x1 << 11
0024 TORCH = 0x1 << 12
0025
0026 AB = BULK_ABSORB.bit_length()
0027 SD = SURFACE_DETECT.bit_length()
0028 SA = SURFACE_ABSORB.bit_length()
0029 DR = SURFACE_DREFLECT.bit_length()
0030 SR = SURFACE_SREFLECT.bit_length()
0031 BR = BOUNDARY_REFLECT.bit_length()
0032 BT = BOUNDARY_TRANSMIT.bit_length()
0033 TO = TORCH.bit_length()
0034
0035 N = int(os.environ.get("VERSION", "-1"))
0036 CMDLINE = "N=%d ./U4SimulateTest.sh mt" % N
0037
0038 TEST = os.environ.get("TEST", "Manual")
0039 GEOM = os.environ.get("GEOM", "DummyGEOM")
0040 GEOMList = os.environ.get("%s_GEOMList" % GEOM, "DummyGEOMList")
0041
0042 VERSION = N
0043 MODE =  int(os.environ.get("MODE", "2"))
0044 assert MODE in [0,2,3]
0045 PIDX = int(os.environ.get("PIDX", "123"))
0046
0047
0048 if MODE > 0:
0049     from opticks.ana.pvplt import *
0050 pass
0051
0052 from opticks.u4.tests.ModelTrigger_Debug import ModelTrigger_Debug
0053
0054 axes = 0, 2  # X,Z
0055 H,V = axes
0056
0057 if __name__ == '__main__':
0058     t = Fold.Load(symbol="t")
0059     print(repr(t))
0060
0061     end = t.photon[:,0,:3]
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_)               # array of step point flags
0065     n = np.sum( seqnib_(q_), axis=1 ) # nibble count
0066
0067     w_sr = np.where( qq == SR )       # selecting step points with SR
0068     w_to = np.where( qq == TO )
0069     w_br = np.where( qq == BR )
0070
0071     sr_pos = t.record[tuple(w_sr[0]), tuple(w_sr[1]),0,:3]   # step point positions of SR
0072     to_pos = t.record[tuple(w_to[0]), tuple(w_to[1]),0,:3]
0073     br_pos = t.record[tuple(w_br[0]), tuple(w_br[1]),0,:3]
0074
0075
0076     #LAYOUT = os.environ.get("LAYOUT", "")
0077     ## ana env may not match run env, making it fragile to assume so
0078     LAYOUT = t.photon_meta.LAYOUT[0]
0079     CHECK = t.photon_meta.CHECK[0]
0080
0081     print("TEST:%s" % (TEST) )
0082     print("CMDLINE:%s" % (CMDLINE) )
0083     print("LAYOUT:%s" % (LAYOUT) )
0084     print("CHECK:%s" % (CHECK) )
0085     print("MODE:%d" % (MODE) )
0086     print("PIDX:%d" % (PIDX) )
0087     print("N:%d" % (N) )
0088
0089     mtd = ModelTrigger_Debug(t, symbol="mtd", publish=False)  # publish:True crashing
0090     print(mtd)
0091
0092     if 'SPECS' in os.environ:
0093         SPECS = np.array(t.U4R_names.lines)     # step specification, for skip fake debugging
0094         st_ = t.aux[:,:,2,3].view(np.int32)
0095         st = SPECS[st_]
0096         st_dump = True
0097     else:
0098         SPEC, st_, st = None, None, None
0099         st_dump = False
0100     pass
0101
0102     if st_dump:
0103         u_st, n_st = np.unique(st, return_counts=True)
0104         expr = "np.c_[n_st,u_st][np.argsort(n_st)[::-1]]"
0105         print(expr)
0106         print(eval(expr))
0107     pass
0108
0109     exprs = r"""
0110     q[PIDX]
0111     t.record[PIDX,:n[PIDX],0]
0112     mtd.pv[mtd.index==PIDX]
0113     np.unique(mtd.whereAmI[mtd.trig==1],return_counts=True)
0114     """
0115     exprs_ = list(filter(None,textwrap.dedent(exprs).split("\n")))
0116     for expr in exprs_:
0117         print("\n%s ## " % expr)
0118         print(eval(expr))
0119     pass
0120
0121
0122     mtd_outside = np.logical_and(mtd.trig == 1, mtd.EInside1 == 0 )
0123
0124     mtd_trig = mtd.trig == 1
0125     mtd_upper = mtd.pos[:,2] > 1e-4
0126     mtd_mid   = np.abs( mtd.pos[:,2]) < 1e-4
0127     mtd_lower = mtd.pos[:,2] < -1e-4
0128     mtd_pyrex  = mtd.whereAmI_ == 1
0129     mtd_vacuum = mtd.whereAmI_ == 2
0130
0131     mtd_trig_vacuum = np.logical_and(mtd_trig, mtd_vacuum)
0132     mtd_trig_vacuum_upper = np.logical_and(mtd_trig_vacuum, mtd_upper )
0133     mtd_trig_vacuum_mid   = np.logical_and(mtd_trig_vacuum, mtd_mid )
0134     mtd_trig_vacuum_lower = np.logical_and(mtd_trig_vacuum, mtd_lower )
0135
0136     mtd_trig_pyrex  = np.logical_and(mtd_trig, mtd_pyrex )
0137     mtd_trig_pyrex_upper = np.logical_and(mtd_trig_pyrex, mtd_upper )
0138     mtd_trig_pyrex_mid   = np.logical_and(mtd_trig_pyrex, mtd_mid )
0139     mtd_trig_pyrex_lower = np.logical_and(mtd_trig_pyrex, mtd_lower )
0140
0141
0142
0143     rqwns = textwrap.dedent("""
0144     GEOM
0145     GEOMList
0146     LAYOUT
0147     CHECK
0148     TEST
0149     """)
0150
0151     lqwns = textwrap.dedent("""
0152     mtd.IMPL
0153     t.photon.shape
0154     mtd.pos.shape
0155
0156     mtd_trig
0157     mtd_upper
0158     mtd_mid
0159     mtd_lower
0160
0161     mtd_pyrex
0162     mtd_vacuum
0163
0164
0165
0166     mtd_trig_pyrex
0167     mtd_trig_pyrex_upper
0168     mtd_trig_pyrex_mid
0169     mtd_trig_pyrex_lower
0170
0171     mtd_trig_vacuum
0172     mtd_trig_vacuum_upper
0173     mtd_trig_vacuum_mid
0174     mtd_trig_vacuum_lower
0175     """)
0176
0177
0178     llines = []
0179     for qwn in lqwns.split("\n"):
0180         if len(qwn) == 0:
0181             llines.append("")
0182         elif qwn.find(".") > -1:
0183             llines.append("%30s : %s" % ( qwn, eval(qwn)))
0184         else:
0185             num = np.count_nonzero(eval(qwn))
0186             llines.append("%30s : %d " % ( qwn, num ))
0187         pass
0188     pass
0189     lanno = "\n".join(llines)
0190     print(lanno)
0191
0192     if not "NANNO" in os.environ:
0193         os.environ["LHSANNO"] = lanno
0194     pass
0195
0196     rlines = []
0197     for qwn in rqwns.split("\n"):
0198         line = "" if len(qwn) == 0 else "%s : %s" % ( qwn, eval(qwn))
0199         rlines.append(line)
0200     pass
0201     ranno = "\n".join(rlines)
0202     print(ranno)
0203     if not "NANNO" in os.environ:
0204         os.environ["RHSANNO"] = ranno
0205     pass
0206
0207
0208
0209
0210
0211
0212     idxs = np.unique(mtd.index[mtd_trig_pyrex_lower])   # photon indices
0213
0214     flagmask = t.photon[:,3,3].view(np.int32)
0215     sd = flagmask & SURFACE_DETECT != 0
0216     sa = flagmask & SURFACE_ABSORB != 0
0217
0218     ## branch on layout as coordinates are specific to each (eg midline of the left hand PMT in two_pmt)
0219     if LAYOUT == "two_pmt":
0220         x_midline = np.logical_and( end[:,0] > -251, end[:,0] < -249 )
0221         z_midline = np.logical_and( end[:,2] > -250, end[:,2] <  250 )
0222         xz_midline = np.logical_and( x_midline, z_midline )
0223     elif LAYOUT == "one_pmt":
0224         x_midline = np.logical_and( end[:,0] > -245, end[:,0] < 245 )
0225         z_midline = np.logical_and( end[:,2] > -1,   end[:,2] <  1 )
0226         xz_midline = np.logical_and( x_midline, z_midline )
0227     else:
0228         xz_midline = None
0229     pass
0230     w_midline = np.where(xz_midline)[0]
0231
0232
0233     nn = n[w_midline]
0234
0235     end2 = t.record[w_midline, tuple(nn-1), 0, :3]   # recreate photon endpoint pos from record array
0236     assert(np.all( end[w_midline] == end2 ))
0237
0238     penultimate =  t.record[w_midline, tuple(nn-2), 0, :3]
0239     prior =  t.record[w_midline, tuple(nn-3), 0, :3]
0240
0241
0242     #ppos0_ = "None"
0243     #ppos0_ = "end"
0244     #ppos0_ = "end[sd] # photon SD endpoints around the upper hemi"
0245     #ppos0_ = "end[sa] # photon SA endpoints around the upper hemi and elsewhere"
0246     #ppos0_ = "end[w_midline]  # photons ending on midline "
0247     #ppos0_ = "mtd.pos[mtd_outside] # just around upper hemi "
0248     ppos0_  = "mtd.pos[mtd_trig]"
0249     #ppos0_ = "mtd.pos[mtd_trig_pyrex]  # Simple:just around upper hemi, Buggy:also dynode/MCP sprinkle "
0250     #ppos0_ = "mtd.pos[mtd_trig_pyrex_lower] # "
0251     #ppos0_ = "mtd.pos[mtd_trig_vacuum] # mostly on midline, sprinkle of obliques around upper hemi "
0252     #ppos0_ = "mtd.next_pos[mtd_trig] # just around upper hemi"
0253     #ppos0_ = "sr_pos # SR positions  "
0254     #ppos0_ = "to_pos # TO positions  "
0255     #ppos0_ = "br_pos # BR positions  "
0256     #ppos0_ = "mtd.pos[mtd_pyrex] #   "
0257     #ppos0_ = "mtd.pos[mtd_trig_pyrex] #   "
0258
0259     #ppos1_ = "None"
0260     #ppos1_ = "mtd.pos[mtd_trig_vacuum_upper]"
0261     #ppos1_ = "mtd.next_pos[mtd_trig_vacuum_upper]"
0262     ppos1_ = "mtd.pos[mtd_trig_pyrex_lower] # BUG : Pyrex triggers inside inner2 : UNPHYSICAL  "
0263     #ppos1_  = "end[xz_midline]"
0264     #ppos1_  = "penultimate  # photon position prior to terminal one"
0265     #ppos1_  = "prior  # two positions before last"
0266
0267     ppos2_ = "None"
0268
0269
0270     if TEST == "quiver":
0271         qsel = "mtd_trig_vacuum_upper"
0272         #qsel = "np.where(mtd_trig_vacuum_upper)[0][:10]"
0273         #qsel = "mtd_trig"
0274         ppos0_ = "mtd.pos[%(qsel)s] #QUIVER "  % locals()
0275         ppos1_ = "mtd.dir[%(qsel)s] #QUIVER "  % locals()
0276         ppos2_ = "mtd.next_pos[%(qsel)s] #QUIVER " % locals()
0277     else:
0278         print("using default ppos0_ ppos1_ ppos2_ exprs ")
0279     pass
0280
0281
0282     ppos0  = eval(ppos0_)
0283     ppos1  = eval(ppos1_)
0284     ppos2  = eval(ppos2_)
0285
0286     elem = []
0287     elem.append(CMDLINE)
0288     if not ppos0 is None: elem.append("blue:%s" % ppos0_)
0289     if not ppos1 is None: elem.append("red:%s" % ppos1_)
0290     if not ppos2 is None: elem.append("green:%s" % ppos2_)
0291     label = "\n".join(elem)
0292
0293
0294     exprs = textwrap.dedent("""
0295     np.count_nonzero(np.logical_and(np.logical_and(mtd.trig==1,mtd.whereAmI_==2),mtd.EInside1==0))
0296     np.count_nonzero(np.logical_and(np.logical_and(mtd.trig==1,mtd.whereAmI_==2),mtd.EInside1==1))
0297     np.count_nonzero(np.logical_and(np.logical_and(mtd.trig==1,mtd.whereAmI_==2),mtd.EInside1==2))
0298     np.count_nonzero(np.logical_and(np.logical_and(mtd.trig==1,mtd.whereAmI_==1),mtd.EInside1==0))
0299     np.count_nonzero(np.logical_and(np.logical_and(mtd.trig==1,mtd.whereAmI_==1),mtd.EInside1==1))
0300     np.count_nonzero(np.logical_and(np.logical_and(mtd.trig==1,mtd.whereAmI_==1),mtd.EInside1==2))
0301     """)
0302
0303     for expr in exprs.split():
0304         print(expr)
0305         print(eval(expr))
0306     pass
0307
0308
0309     if MODE == 0:
0310         print("not plotting as MODE 0  in environ")
0311     elif MODE == 2:
0312         fig, axs = mpplt_plotter(label=label)
0313         ax = axs[0]
0314         ax.set_ylim(-250,250)
0315         ax.set_xlim(-500,500)
0316
0317         if "#QUIVER" in ppos0_ and "#QUIVER" in ppos1_ and "#QUIVER" in ppos2_:
0318             assert not ppos0 is None
0319             assert not ppos1 is None
0320             assert not ppos2 is None
0321             ax.quiver( ppos0[:,H], ppos0[:,V],  ppos1[:,H], ppos1[:,V], units="width", width=0.0002, scale=10.0 )
0322             ax.scatter( ppos0[:,H], ppos0[:,V], s=1, c="r")
0323             ax.scatter( ppos2[:,H], ppos2[:,V], s=1, c="g")
0324         else:
0325             if not ppos0 is None: ax.scatter( ppos0[:,H], ppos0[:,V], s=1 )
0326             if not ppos1 is None: ax.scatter( ppos1[:,H], ppos1[:,V], s=1, c="r" )
0327         pass
0328         fig.show()
0329
0330     elif MODE == 3:
0331         pl = pvplt_plotter(label)
0332         os.environ["EYE"] = "0,100,165"
0333         os.environ["LOOK"] = "0,0,165"
0334         pvplt_viewpoint(pl)
0335         pl.add_points(ppos0)
0336         pl.show()
0337     pass
0338 pass