g4cx/tests/G4CXSimtraceTest.py

0001 #!/usr/bin/env python
0002 """
0003 G4CXSimtraceTest.py
0004 =====================
0005
0006 TODO: further reduce duplication with cx/tests/CSGOptiXSimtraceTest.py
0007
0008 Used from gxt.sh, eg::
0009
0010    SIMPLE=1 ./gxt.sh ana
0011    ./gxt.sh ana
0012
0013    MASK=non ./gxt.sh ana
0014        when the genstep grid is too small the view with default MASK=pos
0015        will be very zoomed in, use MASK=non to not mask the positions plotted
0016        by the genstep grid points
0017
0018    MASK=t ./gxt.sh ana
0019        MASK=t requires t>0 which excludes miss, but does not restrict to
0020        the genstep grid positions so can see intersects from outside the grid
0021
0022
0023 """
0024 import os, numpy as np, logging
0025 log = logging.getLogger(__name__)
0026
0027 from opticks.ana.fold import Fold
0028 from opticks.ana.p import *
0029 # including cf loaded from CFBASE
0030 # HMM: kinda bad form to instanciate cf from the module load
0031
0032 from opticks.ana.feature import SimtraceFeatures
0033 from opticks.ana.framegensteps import FrameGensteps
0034 from opticks.ana.simtrace_positions import SimtracePositions
0035 from opticks.ana.simtrace_plot import SimtracePlot, pv, mp
0036 from opticks.ana.pvplt import *
0037 from opticks.ana.eget import efloatlist_, elookce_, elook_epsilon_, eint_
0038
0039 from opticks.npy.mortonlib.morton2d import morton2d
0040
0041
0042
0043
0044 SPURIOUS = "SPURIOUS" in os.environ
0045 RERUN = "RERUN" in os.environ
0046 SELECTION = "SELECTION" in os.environ
0047 SIMPLE = "SIMPLE" in os.environ
0048 MASK = os.environ.get("MASK", "pos")
0049 FEAT = os.environ.get("FEAT", "pid" )
0050 PIDX = eint_("PIDX", "0")
0051
0052
0053 def simple(pl, gs, pos):
0054     """
0055     :param pl: pvplt_plotter instance
0056     :param gs: FrameGensteps instance
0057     :param pos: SimtracePositions instance
0058     """
0059     pvplt_simple(pl, gs.centers_local[:,:3], "gs.centers_local[:,:3]" )
0060     pvplt_simple(pl, pos.gpos[:,:3], "pos.gpos[:,:3]" )
0061     pvplt_simple(pl, pos.lpos[:,:3], "pos.lpos[:,:3]" )
0062
0063
0064
0065
0066 def spurious_2d_outliers(bbox, upos):
0067     """
0068     :param bbox: array of shape (2,3) eg from t.sframe.bbox
0069     :param upos: local frame positions expected to be within bbox, eg from t_pos.upos
0070
0071     """
0072     xbox = bbox[1] - bbox[0]
0073     fpos = np.zeros( (len(upos), 3), dtype=np.float32 )
0074     dim = []
0075     for d in [X,Y,Z]:
0076         if xbox[d] > 0.:
0077             dim.append(d)
0078             fpos[:,d] = (upos[:,d] - bbox[0,d])/xbox[d]
0079         pass
0080     pass
0081     assert len(dim) == 2, "expecting intersects to be in 2D plane"
0082     assert len(np.where(fpos > 1)[0]) == 0, "SPURIOUS running needs all isect within bbox, eg use MASK=t "
0083     assert len(np.where(fpos < 0)[0]) == 0, "SPURIOUS running needs all isect within bbox, eg use MASK=t "
0084     pass
0085     ipos = np.array( fpos*0xffffffff , dtype=np.uint64 )   ## 32bit range integer coordinates stored in 64 bit
0086     kpos = morton2d.Key(ipos[:,dim[0]], ipos[:,dim[1]])    ## morton interleave the two coordinates into one 64 bit code
0087
0088
0089     SPURIOUS_CUT = 1 + int(os.environ.get("SPURIOUS", "1"))
0090
0091
0092     ## scrub low bits and apply uniquing as data reduction : how many bits determines the coarseness
0093     u_kpos_0, i_kpos_0, c_kpos_0 = np.unique( kpos & ( 0xfff << 52 ), return_index=True, return_counts=True)
0094     sel = c_kpos_0 < SPURIOUS_CUT
0095
0096     ## finding outliers in 2d is reduced to finding outliers in sorted list of uint
0097     ## and can also use the counts : outliers are expected to be low count
0098
0099     u_kpos = u_kpos_0[sel]
0100     i_kpos = i_kpos_0[sel]  # original upos index of the unique
0101
0102     if len(i_kpos) < 10:
0103         log.info("spurious_2d_outliers SPURIOUS_CUT %d ", SPURIOUS_CUT )
0104         log.info("i_kpos\n%s" % str(i_kpos))
0105         log.info("upos[i_kpos]\n%s" % str(upos[i_kpos]) )
0106     pass
0107
0108     ## decode the selected isolated morton codes to find the probable spurious intersects
0109     d_kpos = morton2d.Decode(u_kpos)  ## HMM returns tuple of 2 which is bit inconvenient
0110     t_spos = np.zeros( (len(u_kpos),3), dtype=np.float32 )
0111     for idim,d in enumerate(dim):
0112         t_spos[:,d] = d_kpos[idim].astype(np.float32)/np.float32(0xffffffff)
0113         t_spos[:,d] *= xbox[d]
0114         t_spos[:,d] += bbox[0,d]
0115     pass
0116     return u_kpos, c_kpos_0, i_kpos, t_spos
0117
0118
0119
0120
0121 if __name__ == '__main__':
0122
0123     #fmt = '[%(asctime)s] p%(process)s {%(pathname)s:%(lineno)d} %(levelname)s - %(message)s'
0124     fmt = '{%(pathname)s:%(lineno)d} %(levelname)s - %(message)s'
0125     logging.basicConfig(level=logging.INFO, format=fmt)
0126
0127     t = Fold.Load(symbol="t")
0128     a = Fold.Load("$A_FOLD", symbol="a")
0129     b = Fold.Load("$B_FOLD", symbol="b")
0130     print("cf.cfbase : %s " % cf.cfbase)
0131
0132     print("---------Fold.Load.done")
0133     x = a
0134
0135     print(repr(t))
0136     print(repr(a))
0137     print(repr(b))
0138
0139     print("---------print.done")
0140
0141
0142     if not a is None and not a.seq is None:
0143         a_nib = seqnib_(a.seq[:,0])                  # valid steppoint records from seqhis count_nibbles
0144         a_gpos_ = a.record[PIDX,:a_nib[PIDX],0,:3]   # global frame photon step record positions of single PIDX photon
0145         a_gpos  = np.ones( (len(a_gpos_), 4 ) )
0146         a_gpos[:,:3] = a_gpos_
0147         a_lpos = np.dot( a_gpos, t.sframe.w2m )      # a global positions into gxt target frame
0148     else:
0149         a_lpos = None
0150     pass
0151
0152     if not b is None and not b.seq is None:
0153         b_nib = seqnib_(b.seq[:,0])                  # valid steppoint records from seqhis count_nibbles
0154         b_gpos_ = b.record[PIDX,:b_nib[PIDX],0,:3]  # global frame photon step record positions of single PIDX photon
0155         b_gpos  = np.ones( (len(b_gpos_), 4 ) )
0156         b_gpos[:,:3] = b_gpos_
0157         b_lpos = np.dot( b_gpos, t.sframe.w2m )
0158     else:
0159         b_lpos = None
0160     pass
0161
0162     x_lpos = a_lpos
0163
0164
0165
0166
0167     local = True
0168
0169     t_gs = FrameGensteps(t.genstep, t.sframe, local=local, symbol="gs" )  ## get gs positions in target frame
0170     print(t_gs)
0171
0172
0173     if RERUN:
0174         log.info("RERUN envvar switched on use of simtrace_rerun from CSG/SimtraceRerunTest.sh " )
0175         simtrace = t.simtrace_rerun
0176     else:
0177         simtrace = t.simtrace
0178     pass
0179
0180     t_pos = SimtracePositions(simtrace, t_gs, t.sframe, local=local, mask=MASK, symbol="t_pos" )
0181     print(t_pos)
0182
0183
0184     if SPURIOUS:
0185         log.info("SPURIOUS envvar switches on search for morton outliers using spurious_2d_outliers ")
0186         u_kpos, c_kpos, i_kpos, t_spos = spurious_2d_outliers( t.sframe.bbox, t_pos.upos )
0187         j_kpos = t_pos.upos2simtrace[i_kpos]
0188         log.info("j_kpos = t_pos.upos2simtrace[i_kpos]\n%s" % str(t_pos.upos2simtrace[i_kpos]) )
0189         log.info("simtrace[j_kpos]\n%s" % str(simtrace[j_kpos]) )
0190         simtrace_spurious = j_kpos
0191     else:
0192         t_spos = None
0193         simtrace_spurious = []
0194     pass
0195
0196     if SIMPLE:
0197         pl = pvplt_plotter()
0198         simple(pl, t_gs, t_pos)
0199         pl.show()
0200         raise Exception("SIMPLE done")
0201     pass
0202
0203     pf = SimtraceFeatures(t_pos, cf, featname=FEAT, symbol="pf" )
0204
0205     pl = SimtracePlot.MakePVPlotter()
0206
0207     plt = SimtracePlot(pl, pf.feat, t_gs, t.sframe, t_pos, outdir=os.path.join(t.base, "figs") )
0208
0209
0210     if not x_lpos is None:
0211         plt.x_lpos = x_lpos
0212     pass
0213     if not t_spos is None:
0214         plt.t_spos = t_spos     ## spurious intersects identified by morton2d isolation
0215     pass
0216
0217     ## created by CSG/SimtraceRerunTest.sh with SELECTION envvar picking simtrace indices to highlight
0218     ## but the SELECTION envvar used here just has to exist to trigger selection plotting
0219     if hasattr(t, "simtrace_selection") and SELECTION:
0220         plt.simtrace_selection = t.simtrace_selection
0221     elif len(simtrace_spurious) > 0:
0222         plt.simtrace_selection = simtrace[simtrace_spurious]
0223     else:
0224         pass
0225     pass
0226
0227     if not mp is None:
0228         plt.positions_mpplt()
0229         ax = plt.ax
0230     pass
0231
0232     if not pv is None:
0233         plt.positions_pvplt()
0234     pass
0235 pass
0236
0237