u4/tests/U4RecorderTest.py

0001 #!/usr/bin/env python
0002 """
0003 U4RecorderTest.py
0004 ==================
0005
0006 ::
0007
0008      18 struct spho
0009      19 {
0010      20     int gs ; // 0-based genstep index within the event
0011      21     int ix ; // 0-based photon index within the genstep
0012      22     int id ; // 0-based photon identity index within the event
0013      23     int gn ; // 0-based reemission index incremented at each reemission
0014
0015
0016 """
0017 import numpy as np
0018 from opticks.ana.fold import Fold
0019 from opticks.ana.p import *
0020
0021 from opticks.sysrap.stag import stag
0022 from opticks.u4.U4Stack import U4Stack
0023 stack = U4Stack()
0024
0025
0026 def check_pho_labels(l):
0027     """
0028     :param l: spho labels
0029
0030     When reemission is enabled this would fail for pho0 (push_back labels)
0031     but should pass for pho (labels slotted in by event photon id)
0032
0033     1. not expecting gaps in list of unique genstep index gs_u
0034        as there should always be at least one photon per genstep
0035
0036     2. expecting the photon identity index to be unique within event,
0037        so id_c should all be 1, otherwise likely a labelling issue
0038
0039     """
0040     gs, ix, id_, gn = l[:,0], l[:,1], l[:,2], l[:,3]
0041
0042     gs_u, gs_c = np.unique(gs, return_counts=True )
0043     np.all( np.arange( len(gs_u) ) == gs_u )
0044
0045     id_u, id_c = np.unique( id_, return_counts=True  )
0046     all_one =  np.all( id_c == 1 )
0047
0048     if not all_one: print("check_pho_labels all_one FAIL")
0049     #assert all_one
0050
0051     ix_u, ix_c = np.unique( ix, return_counts=True )
0052
0053     gn_u, gn_c = np.unique( gn, return_counts=True )
0054     print(gn_u)
0055     print(gn_c)
0056
0057
0058 def parse_slice(ekey):
0059     if not ekey in os.environ: return slice(None)
0060     elem = os.environ[ekey].split(":")
0061     start = None if elem[0] == "" else int(elem[0])
0062     stop = None if elem[1] == "" else int(elem[1])
0063     step = None if elem[2] == "" else int(elem[2])
0064     return slice(start, stop, step)
0065
0066
0067
0068 if __name__ == '__main__':
0069
0070     t = Fold.Load()
0071     PIDX = int(os.environ.get("PIDX","-1"))
0072     print(t)
0073
0074     if "RECS_PLOT" in os.environ:
0075         """
0076         RECS_PLOT=1 ./U4RecorderTest.sh ana
0077
0078
0079         RECS_PLOT=1 SEQHIS="TO BT BT SA" RSLI="0::4" ./U4RecorderTest.sh ana
0080             select photons with one history and then just plot the first record step
0081
0082         RECS_PLOT=1 SEQHIS="TO BT BT SA" RSLI="1::4" ./U4RecorderTest.sh ana
0083             expected to show positions of first BT
0084
0085             np.all(np.tile(np.arange(4),100)[slice(0,None,4)]==0)
0086             np.all(np.tile(np.arange(4),100)[slice(1,None,4)]==1)
0087             np.all(np.tile(np.arange(4),100)[slice(2,None,4)]==2)
0088             np.all(np.tile(np.arange(4),100)[slice(3,None,4)]==3)
0089
0090         RECS_PLOT=1 SEQHIS="TO BT BT SA" RSLI="2::4" ./U4RecorderTest.sh ana
0091
0092         RECS_PLOT=1 SEQHIS="TO BT BT SA" RSLI="3::4" ./U4RecorderTest.sh ana
0093
0094         RECS_PLOT=1 SEQHIS="TO BT BT SA" IREC=2 ./U4RecorderTest.sh ana
0095
0096
0097         RECS_PLOT=1 SEQHIS="TO BT BT BT BT BT SA" IREC=1 ./U4RecorderTest.sh ana
0098
0099         RECS_PLOT=1 SEQHIS="TO BT BT BT BT BT SA" ./U4RecorderTest.sh ana
0100              plot all step points
0101
0102         """
0103         from opticks.ana.pvplt import *
0104         pl = pvplt_plotter(label="dump all record point positions")
0105         if "SEQHIS" in os.environ:
0106             #w_sel = np.where( t.seq[:,0] == cseqhis_("TO BT BT SA") )[0]
0107             w_sel = np.where( t.seq[:,0] == cseqhis_(os.environ["SEQHIS"]) )[0]
0108         else:
0109             w_sel = slice(None)
0110         pass
0111
0112         if "RSLI" in os.environ:
0113             ## This way of selecting record points is overly complicated, much better
0114             ## to use the more direct IREC approach.
0115             w_rec = np.where( t.record[w_sel,:,3,3].view(np.int32) != 0 )
0116             #recs = t.record[w_rec]   ## flattens step points to give shape eg (35152, 4, 4) where 8788*4  = 35152
0117             ## the above would be an error as it is using w_sel sub-selection indices as if they were full array indices.
0118             recs = t.record[w_sel][w_rec]
0119             ## To correctly use a where selection the base selection must be the same as that
0120             ## used to create the where selection.
0121             rsli = parse_slice("RSLI")
0122             rpos = recs[rsli,0,:3]
0123         elif "IREC" in os.environ:
0124             ## using IREC is a simpler more direct way of selecting step points
0125             irec = int(os.environ["IREC"])
0126             rpos = t.record[w_sel, irec, 0,:3]
0127         else:
0128             w_rec = np.where( t.record[w_sel,:,3,3].view(np.int32) != 0 )
0129             ## Above w_rec skips unfilled record points.
0130             ## Notice that for w_sel slice(None) for example w_rec is selecting
0131             ## different numbers of step points for each photon : so the below rpos contains
0132             ## all positions from all step points within the w_sel selection. This includes
0133             ## bulk scatter and absorption points.
0134             ## TODO: find way to do similar but select based on the flag of the step points
0135             ## so can color plots according to the flag.
0136             rpos = t.record[w_sel][w_rec][:,0,:3]
0137         pass
0138
0139         pl.add_points(rpos)
0140         pl.show()
0141     pass
0142
0143
0144
0145
0146     # pho: labels are collected within U4Recorder::PreUserTrackingAction
0147     l = t.pho if hasattr(t, "pho") else None      # labels slotted in using spho::id
0148     check_pho_labels(l)
0149
0150     gs, ix, id_, gn = l[:,0], l[:,1], l[:,2], l[:,3]
0151
0152     st = stag.Unpack(t.tag) if hasattr(t,"tag") else None
0153
0154     p = t.photon if hasattr(t, "photon") else None
0155     r = t.record if hasattr(t, "record") else None
0156     seq = t.seq if hasattr(t, "seq") else None
0157     nib = seqnib_(seq[:,0])  if not seq is None else None
0158
0159     lim = min(len(p), 3)
0160     for i in range(lim):
0161         if not (PIDX == -1 or PIDX == i): continue
0162         if PIDX > -1: print("PIDX %d " % PIDX)
0163         print("r[%d,:,:3]" % i)
0164         print(r[i,:nib[i],:3])
0165         print("\n\nbflagdesc_(r[%d,j])" % i)
0166         for j in range(nib[i]):
0167             print(bflagdesc_(r[i,j]))
0168         pass
0169
0170         #print("ridiff_(r[%d])*1000." % i)
0171         #print(ridiff_(r[i])*1000.)
0172
0173         print("\n")
0174         print("p[%d]" % i)
0175         print(p[i])
0176         print("\n")
0177         print("bflagdesc_(p[%d])" % i)
0178         print(bflagdesc_(p[i]))
0179         print("\n")
0180         if not seq is None:
0181             print("seqhis_(seq[%d,0]) nib[%d]  " % (i,i) )
0182             print(" %s : %s "% ( seqhis_(seq[i,0]), nib[i] ))
0183             print("\n")
0184         pass
0185         print("\n\n")
0186     pass
0187     idx = p.view(np.uint32)[:,3,2]
0188     check_idx = np.all( np.arange( len(p) ) == idx )
0189     if not check_idx: print("check_idx FAIL")
0190     #assert check_idx
0191
0192     flagmask_u, flagmask_c = np.unique(p.view(np.uint32)[:,3,3], return_counts=True)
0193     print("flagmask_u:%s " % str(flagmask_u))
0194     print("flagmask_c:%s " % str(flagmask_c))
0195
0196     #print("\n".join(seqhis_( t.seq[:,0] )))
0197     for i in range(lim):
0198         print("%4d : %s " % (i, seqhis_(t.seq[i,0])))
0199     pass
0200
0201     seq_u, seq_c = np.unique(seq[:,0], return_counts=True)
0202     for i in range(len(seq_u)):
0203         print("%4d : %4d : %20s  " % (i, seq_c[i], seqhis_(seq_u[i]) ))
0204     pass
0205     print("t.base : %s " %  t.base)
0206 pass
0207