u4/tests/U4SimulateTest_fs.py

0001 #!/usr/bin/env python
0002 """
0003 U4SimulateTest_fs.py  : FastSim ModelTrigger etc.. plotting
0004 ===============================================================
0005
0006
0007 """
0008 import os, numpy as np
0009 from opticks.ana.fold import Fold
0010 from opticks.ana.p import *
0011
0012 MODE = int(os.environ.get("MODE", 0))
0013 if MODE > 0:
0014     from opticks.ana.pvplt import *
0015 pass
0016
0017
0018
0019
0020 if __name__ == '__main__':
0021     t = Fold.Load(symbol="t")
0022     print(repr(t))
0023
0024     uu = t.U4Recorder_G4StateRerun_726
0025
0026     fs = t.SFastSim_Debug
0027     kInfinity = 9.000e+99
0028
0029     pos = fs[:,0,:3]
0030     tim = fs[:,0,3]
0031
0032     mom = fs[:,1,:3]
0033     ds1 = fs[:,1,3]
0034
0035     pol = fs[:,2,:3]
0036     ds2 = fs[:,2,3]
0037
0038     trg = fs[:,3,0].astype(np.int64)  ## NB not ".view", as here wasting bits
0039     wai = fs[:,3,1].astype(np.int64)
0040     c   = fs[:,3,2]
0041     pid = fs[:,3,3].astype(np.int64)
0042
0043     ms1 = np.where(ds1 == kInfinity )[0]
0044     ht1 = np.where(ds1 != kInfinity )[0]
0045     ms2 = np.where(ds2 == kInfinity )[0]
0046     ht2 = np.where(ds2 != kInfinity )[0]
0047
0048     #  enum EWhereAmI { OutOfRegion, kInGlass, kInVacuum };
0049     glass      = wai == 1
0050     vacuum     = wai == 2
0051     trg_no     = trg == 0
0052     trg_yes    = trg == 1
0053     glass_no   = np.logical_and( wai == 1, trg == 0 )
0054     glass_yes  = np.logical_and( wai == 1, trg == 1 )
0055     vacuum_no  = np.logical_and( wai == 2, trg == 0 )
0056     vacuum_yes = np.logical_and( wai == 2, trg == 1 )
0057
0058
0059     PID = int(os.environ.get("PID", -1))
0060     if PID == -1: PID = int(os.environ.get("OPTICKS_G4STATE_RERUN", -1))
0061
0062     pidsel = np.where( pid == PID )[0]
0063     #pidsel = np.where( pid == PID )[0]
0064     num_pidsel = len(pidsel)
0065     print(" PID %d pidsel %s num_pidsel %d " % (PID, str(pidsel), num_pidsel))
0066
0067     if PID > -1 and num_pidsel > 0:
0068         extra = np.zeros( [num_pidsel,4,4], dtype=np.float64 )
0069         extra[:] = fs[pidsel, :4]                   # copy, not reference, fs
0070         extra[np.where( extra == kInfinity )] = -1  # replace obnoxious 9e99
0071
0072         pidsel_path = os.path.expandvars("$FOLD/py_pidsel_%d.npy" % PID)
0073         dirpath = os.path.dirname(pidsel_path)
0074         if not os.path.isdir(dirpath):
0075             os.makedirs(dirpath)
0076         pass
0077
0078         print("Save PID %d extra %s to pidsel_path %s (used from xxv.sh) " % (PID, str(extra.shape), pidsel_path))
0079         np.save(pidsel_path, extra )
0080
0081         ModelTrigger = extra[np.where(extra[:,3,0].astype(np.int64) == 1)]
0082     pass
0083
0084
0085     selnames = "glass vacuum trg_no trg_yes glass_no glass_yes vacuum_no vacuum_yes".split()
0086     for selname in selnames:
0087         sel = locals()[selname]
0088         print(" %20s : %s " % (selname, np.count_nonzero(sel) ))
0089     pass
0090
0091     #specs = "blue:trg_no:lines cyan:trg_yes:lines"
0092     #specs = "blue:glass_no:lines cyan:glass_yes:lines"
0093     specs = "blue:glass_no:lines cyan:glass_yes:lines red:vacuum_yes:lines yellow:vacuum_no:lines"
0094     label = "U4SimulateTest_fs.py:SFastSim_Debug " + specs
0095
0096     if MODE == 0:
0097         print("not plotting as MODE 0  in environ")
0098     else:
0099         pl = pvplt_plotter(label)
0100         os.environ["EYE"] = "0,100,165"
0101         os.environ["LOOK"] = "0,0,165"
0102         pvplt_viewpoint(pl)
0103
0104         for spec in specs.split():
0105             color,selname,stylename = spec.split(":")
0106             sel = locals()[selname]
0107             numsel = np.count_nonzero(sel)
0108             if numsel == 0: continue
0109
0110             if stylename == "lines":
0111                 pvplt_lines(pl, pos[sel], mom[sel], color=color, factor=10 )
0112             elif stylename == "points":
0113                 pl.add_points( pos[sel], color=color )
0114             else:
0115                 pass
0116             pass
0117         pass
0118         pl.show()
0119     pass
0120 pass