eic-opticks/ana/ht.py

0001 #!/usr/bin/env python
0002 """
0003 ht.py : quick checks on hits
0004 ===================================
0005
0006 ::
0007
0008     In [2]: ox.view(np.int32)[:,3]
0009     Out[2]:
0010     array([[ -24,   -1,    0, 6400],
0011            [  41,   -1,    1, 6152],
0012            [  21,   -1,    2, 6656],
0013            ...,
0014            [  20,   -1, 9997, 6272],
0015            [  20,   -1, 9998, 6272],
0016            [  41,   -1, 9999, 6152]], dtype=int32)
0017
0018
0019     In [5]: np.all( ox.view(np.int32)[:,3,2] == np.arange(10000) )
0020     Out[5]: True
0021
0022
0023     In [10]: ox_flags[np.where(ox_flags[:,3] & hismask.code("SD"))]
0024     Out[10]:
0025     array([[ -30,  130,   19, 6208],
0026            [ -30,  139,   74, 6208],
0027            [ -30,  188,  217, 6208],
0028            [ -30,   76,  406, 6224],
0029            [ -30,  189,  546, 6208],
0030            [ -30,  181,  586, 7232],
0031            [ -30,  167,  690, 6208],
0032            [ -30,  185,  767, 6208],
0033              bnd sen-idx ox-idx hismask(OR of step flags)
0034
0035     In [29]: ht_flags[:,1].min()
0036     Out[29]: 8
0037
0038     In [30]: ht_flags[:,1].max()
0039     Out[30]: 189
0040
0041
0042 ::
0043
0044     In [36]: ht_flags.shape
0045     Out[36]: (87, 4)
0046
0047     In [35]: ox_flags[ox_flags[:,1] != -1].shape
0048     Out[35]: (431, 4)
0049     ## many photons land on sensors but are not classified as hits
0050
0051
0052 ::
0053
0054     In [43]: ox_lander[:,1].min()
0055     Out[43]: 0
0056
0057     In [44]: ox_lander[:,1].max()
0058     Out[44]: 190
0059
0060
0061
0062
0063 """
0064
0065 import os, sys, logging, numpy as np
0066 log = logging.getLogger(__name__)
0067
0068 from opticks.ana.histype import HisType
0069 from opticks.ana.mattype import MatType
0070 from opticks.ana.hismask import HisMask
0071 from opticks.ana.blib import BLib
0072
0073 histype = HisType()
0074 mattype = MatType()
0075 hismask = HisMask()
0076 blib = BLib()
0077
0078
0079 if __name__ == '__main__':
0080     logging.basicConfig(level=logging.INFO)
0081     if len(sys.argv) > 1 and os.path.isdir(sys.argv[1]):
0082         os.chdir(sys.argv[1])
0083         log.info("chdir %s " % os.getcwd())
0084     pass
0085     np.set_printoptions(suppress=True, linewidth=200)
0086
0087     ox = np.load("ox.npy")
0088     ht = np.load("ht.npy") # ht are a selection of the ox:photons with SD:SURFACE_DETECT flag set
0089
0090     ph = np.load("ph.npy") # seqhis, seqmat sequence histories for all photons
0091     seqhis = ph[:,0,0]
0092     seqmat = ph[:,0,1]
0093
0094     ## hmm need a standard place for detector level stuff like this
0095     sd = np.load(os.path.expandvars("$TMP/G4OKTest/sensorData.npy"))
0096     triplet_id = sd.view(np.uint32)[:,3]
0097     placement_id = ( triplet_id & 0x00ffff00 ) >> 8
0098
0099     dump = True
0100     if dump:
0101         for _ in seqhis[:10]:print("%16x : %s " %(_,histype.label(_)))
0102         for _ in seqmat[:10]:print("%16x : %s " %(_,mattype.label(_)))
0103         for htf in ht[:10,3].view(np.int32):print(htf, hismask.label(htf[3]))
0104     pass
0105
0106     ht_flags = ht.view(np.int32)[:,3]
0107     ox_flags = ox.view(np.int32)[:,3]
0108     ht_flags2 = ox_flags[np.where(ox_flags[:,3] & hismask.code("SD"))]
0109     # duplicate hit flags from the photons, using the .w flag mask
0110     # and knowledge that hit selection requires SD:SURFACE_DETECT
0111     assert np.all( ht_flags == ht_flags2 )
0112
0113     ht_idx = ht_flags[:,2]
0114     ox_lander = ox_flags[ox_flags[:,1] != -1]
0115
0116     print("ox_flags : %s " % repr(ox_flags.shape) )
0117     print("ht_flags : %s " % repr(ht_flags.shape) )
0118     print("ox_lander : %s : photons landing on sensor volumes  " % repr(ox_lander.shape))
0119
0120     bndidx = np.unique(np.abs(ox_lander[:,0])) - 1    # subtract 1 to get index as signed boundaries are 1-based
0121     print("bndidx : %s " % repr(bndidx))
0122     for _ in bndidx:print(blib.bname(_))
0123
0124
0125     tot = 0
0126     for oxf in ox_lander:
0127         sidx = oxf[1] # sensor index
0128         tid = triplet_id[sidx]
0129         pid = placement_id[sidx]
0130         idx = oxf[2]  # photon index
0131         sqh = seqhis[idx]
0132         sqm = seqmat[idx]
0133         is_hit = idx in ht_idx
0134         tot += int(is_hit)
0135         stat = "HIT" if is_hit else "-"
0136         print("%20s %15s %10s    %16x : %30s    %16s : %30s  %10x %3x %3d  " % (oxf, hismask.label(oxf[3]), stat, sqh, histype.label(sqh), sqm, mattype.label(sqm), tid, pid, pid ))
0137     pass
0138     assert tot == len(ht_idx)
0139