eic-opticks/ana/ox.py

0001 #!/usr/bin/env python
0002 """
0003 ox.py : quick checks on photons
0004 ===================================
0005
0006 ::
0007
0008     In [36]: ox.view(np.int32)[:,3]
0009     Out[36]:
0010     array([[-1507329,     3159,        0,     6400],
0011            [ 2752511,     4430,        1,     6152],
0012            [ 1441791,     4425,        2,     6656],
0013            ...,
0014            [ 1376255,     3155,     4997,     6272],
0015            [-1376257,     3157,     4998,     6416],
0016            [ 1376255,     3155,     4999,     6272]], dtype=int32)
0017
0018            bnd_sidx
0019            two int16       nidx     phidx     flags
0020
0021     In [37]: np.all( ox.view(np.int32)[:,3,2] == np.arange(5000) )
0022     Out[37]: True
0023
0024     In [38]: ox_flags[np.where(ox_flags[:,3] & hismask.code("SD"))]
0025     Out[38]:
0026     array([[-1965950,     3981,       19,     6208],
0027            [-1965941,     4035,       74,     6208],
0028            [-1965892,     4329,      217,     6208],
0029            [-1966004,     3657,      406,     6224],
0030            [-1965891,     4335,      546,     6208],
0031            [-1965899,     4287,      586,     7232],
0032            [-1965913,     4203,      690,     6208],
0033
0034     In [41]: ox_flags[np.where(ox_flags[:,3] & hismask.code("SD"))][:,0] & 0xffff
0035     Out[41]:
0036     array([130, 139, 188,  76, 189, 181, 167, 185, 152, 150,  29,  89,  97, 160, 183,  37, 132,  50,  13, 169, 141,  84,  73,  85, 144, 128,  87,  19, 187, 174,  76, 180, 101,  82, 116,  66,  29,  63,
0037             88, 165,  36, 169,   9, 121,  23, 129, 143], dtype=int32)
0038
0039     In [42]: ox_flags[np.where(ox_flags[:,3] & hismask.code("SD"))][:,0] >> 16
0040     Out[42]:
0041     array([-30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30, -30,  30,  30,  30, -30, -30, -30, -30, -30, -30, -30, -30,  30, -30, -30,
0042            -30, -30, -30, -30,  30, -30, -30, -30, -30], dtype=int32)
0043
0044
0045 """
0046
0047 import os, sys, logging, numpy as np
0048 log = logging.getLogger(__name__)
0049
0050 from opticks.ana.histype import HisType
0051 from opticks.ana.mattype import MatType
0052 from opticks.ana.hismask import HisMask
0053 from opticks.ana.blib import BLib
0054 from opticks.ana.ggeo import GGeo
0055
0056 histype = HisType()
0057 mattype = MatType()
0058 hismask = HisMask()
0059 blib = BLib()
0060 ggeo = GGeo()
0061
0062
0063 def dump_boundaries(ox):
0064     bndidx = (ox[:,3,0].view(np.uint32) >> 16).view(np.int16)[0::2]
0065     u_bndidx, u_bndidx_counts = np.unique(bndidx, return_counts=True)
0066     tot = 0
0067     print("dump_boundaries")
0068     for bnd,bnd_count in sorted(zip(u_bndidx,u_bndidx_counts), key=lambda _:_[1], reverse=True):
0069         name = blib.bname(np.abs(bnd)-1)  # subtract 1 to get index as signed boundaries are 1-based
0070         print("%4d : %7d  : %s " % (bnd, bnd_count, name))
0071         tot += bnd_count
0072     pass
0073     print("%4s : %7d " % ("TOT",tot))
0074
0075 def dump_sensorIndex(ox):
0076     sidx = (ox[:,3,0].view(np.uint32) & 0xffff).view(np.int16)[0::2]
0077     u_sidx, u_sidx_counts = np.unique(sidx, return_counts=True)
0078     tot = 0
0079     print("dump_sensorIndex")
0080     for sid,sid_count in sorted(zip(u_sidx,u_sidx_counts), key=lambda _:_[1], reverse=True):
0081         print("%4d : %7d  : %s " % (sid, sid_count, ""))
0082         tot += sid_count
0083     pass
0084     print("%4s : %7d " % ("TOT",tot))
0085
0086
0087 if __name__ == '__main__':
0088     logging.basicConfig(level=logging.INFO)
0089     if len(sys.argv) > 1 and os.path.isdir(sys.argv[1]):
0090         os.chdir(sys.argv[1])
0091         log.info("chdir %s " % os.getcwd())
0092     pass
0093     np.set_printoptions(suppress=True, linewidth=200)
0094
0095     ox = np.load("ox.npy")
0096     ph = np.load("ph.npy") # seqhis, seqmat sequence histories for all photons
0097     seqhis = ph[:,0,0]
0098     seqmat = ph[:,0,1]
0099
0100     ox_flags = ox.view(np.int32)[:,3]
0101     ox_lander = ox_flags[ox_flags[:,1] != -1]
0102
0103     print("ox_flags : %s " % repr(ox_flags.shape) )
0104     print("ox_lander : %s : photons landing on sensor volumes  " % repr(ox_lander.shape))
0105
0106     dump_boundaries(ox)
0107     #dump_sensorIndex(ox)
0108
0109     for i, oxr in enumerate(ox):
0110         oxf = oxr[3].view(np.int32)
0111
0112         # see okc/OpticksPhotonFlags optixrap/cu/generate.cu
0113         bnd_sidx,nidx,idx,pflg  = oxf   ## nidx3 will soon become "the one"
0114
0115         nrpo = ggeo.get_triplet_index(nidx)
0116         nidx2,ridx,pidx,oidx = nrpo
0117         assert nidx2 == nidx
0118         #if ridx > 0: continue   # skip photons with last intersect on instanced geometry
0119         if ridx == 0: continue   # skip photons with last intersect on remainder geometry
0120
0121         bnd = np.int16(bnd_sidx >> 16)
0122         sidx = np.int16(bnd_sidx & 0xffff)
0123
0124         sqh = seqhis[idx]  # photon index
0125         sqm = seqmat[idx]
0126         msk = " %15s " % hismask.label(pflg)
0127         his = "( %16x : %30s ) " % (sqh, histype.label(sqh))
0128         mat = "( %16x : %30s ) " % (sqm, mattype.label(sqm))
0129         print(" %5d : %6s %6s : %15s :  %s %s %s : %s " % (i, bnd, sidx, oxf[1:], msk,his,mat, nrpo) )
0130     pass
0131
0132     dump_boundaries(ox)
0133     #dump_sensorIndex(ox)
0134
0135
0136