eic-opticks/ana/G4OKTest.py

0001 #!/usr/bin/env python
0002 """
0003 G4OKTest.py
0004 ============
0005
0006 Transform global positions of photons hitting sensors
0007 into the local frame and plot them all together.
0008
0009 TODO:
0010
0011 * show photon directions
0012 * plot step records
0013 * SDF distances of hits
0014
0015 """
0016 import os, logging, sys, numpy as np
0017 log = logging.getLogger(__name__)
0018 from opticks.ana.hismask import HisMask
0019 from opticks.ana.OpticksIdentity import OpticksIdentity
0020 from opticks.ana.ggeo import GGeo
0021 hismask = HisMask()
0022 ggeo = GGeo()
0023
0024 try:
0025     import pyvista as pv
0026 except ImportError:
0027     pv = None
0028
0029
0030 if __name__ == '__main__':
0031     logging.basicConfig(level=logging.INFO)
0032     if len(sys.argv) > 1 and os.path.isdir(sys.argv[1]):
0033         os.chdir(sys.argv[1])
0034         log.info("chdir %s " % os.getcwd())
0035     pass
0036     np.set_printoptions(suppress=True, linewidth=200)
0037
0038
0039     ## hmm need a standard place for detector level stuff like this
0040     sentid = np.load(os.path.expandvars("$TMP/G4OKTest/sensorData.npy")).view(np.uint32)[:,3]
0041
0042     #name = "ht"
0043     name = "ox"
0044
0045     ox = np.load("%s.npy" % name)
0046     ox_land = ox[ox.view(np.int32)[:,3,1] != -1]   # photons that land on sensors, only some become hits
0047
0048     lpos = np.zeros( (len(ox_land),4), dtype=np.float32 )
0049     ldir = np.zeros( (len(ox_land),4), dtype=np.float32 )
0050     lpol = np.zeros( (len(ox_land),4), dtype=np.float32 )
0051
0052     for i,oxr in enumerate(ox_land):
0053         oxf = oxr[3].view(np.int32)
0054
0055         ## use stomped on weight, for the "always there, not just sensors" node index of last intersected volume
0056         nidx = oxr[1,3].view(np.uint32)
0057         nrpo = ggeo.get_triplet_index(nidx)
0058         assert nrpo[0] == nidx
0059
0060
0061         bnd,sidx,idx,pflg = oxf
0062
0063         tid = sentid[sidx]  # sensor index to triplet id
0064         ## all volumes have a tid, so dont like this as only works for sensors
0065         ##  and depends on the non-standard loading of sensorData.npy
0066
0067         ridx,pidx,oidx = OpticksIdentity.Decode(tid)   # triplet from id
0068         tr = ggeo.get_transform(ridx,pidx,oidx)
0069         it = ggeo.get_inverse_transform(ridx,pidx,oidx)
0070         #it2 = np.linalg.inv(tr)   ## better to do this once within ggeo
0071         #assert np.allclose( it, it2 )
0072
0073         gpos = oxr[0]
0074         gdir = oxr[1]
0075         gpol = oxr[2]
0076
0077         gpos[3] = 1.  # replace time with 1. for transforming a position
0078         gdir[3] = 0.  # replace weight with 0. for transforming a direction
0079         gpol[3] = 0.  # replace wavelength with 0. for transforming a direction
0080
0081         lpos[i] = np.dot( gpos, it )
0082         ldir[i] = np.dot( gdir, it )
0083         lpol[i] = np.dot( gpol, it )
0084
0085         print("bnd/sidx/idx/pflg  %20s %15s  tid %8x  (ridx/pidx/oidx %d %4d %4d)  %s " % (oxf, hismask.label(pflg), tid,ridx,pidx,oidx, lpos[i] ))
0086         #print(tr)
0087     pass
0088
0089
0090     if pv:
0091         pl = pv.Plotter(window_size=2*np.array([1024,768], dtype=np.int32))
0092         pl.add_points(lpos[:,:3])
0093         pl.show_grid()
0094         log.info("Showing the VTK/pyvista plotter window, it may be hidden behind other windows. Enter q to quit.")
0095         cpos = pl.show()
0096         log.info(cpos)
0097
0098
0099
0100
0101