eic-opticks/CSGOptiX/cxt_precision.py

0001 #!/usr/bin/env python
0002 """
0003 cxt_precision.py
0004 ==================
0005
0006
0007 """
0008 import os, sys, re, logging, textwrap, numpy as np
0009 np.set_printoptions(linewidth=200,edgeitems=20)
0010
0011 log = logging.getLogger(__name__)
0012
0013 from opticks.ana.fold import Fold
0014 from opticks.sysrap.sevt import SEvt
0015 import matplotlib.pyplot as mp
0016 SIZE=np.array([1280, 720])
0017
0018
0019 if __name__ == '__main__':
0020     r = SEvt.Load("$RFOLD", symbol="r")
0021     a = SEvt.Load("$AFOLD", symbol="a")
0022
0023     RFOLD = os.environ["RFOLD"]
0024     AFOLD = os.environ["AFOLD"]
0025
0026     FIGPATH = os.environ.get("FIGPATH", None)
0027
0028     iprs = os.environ["OPTICKS_INPUT_PHOTON_RECORD_SLICE"]
0029     iprt = os.environ["OPTICKS_INPUT_PHOTON_RECORD_TIME"]
0030     opr = int(os.environ.get("OPTICKS_PROPAGATE_REFINE", "0"))
0031     oprd = float(os.environ.get("OPTICKS_PROPAGATE_REFINE_DISTANCE", "-1"))
0032
0033     _title = [
0034                "TEST=%s ~/opticks/CSGOptiX/cxt_precision.sh" % os.environ["TEST"],
0035                "OPTICKS_INPUT_PHOTON_RECORD_SLICE %s" % iprs,
0036                "OPTICKS_INPUT_PHOTON_RECORD_TIME %s" % iprt ,
0037                "OPTICKS_PROPAGATE_REFINE %s OPTICKS_PROPAGATE_REFINE_DISTANCE %s" % (opr,oprd)
0038              ]
0039     title = "\n".join(_title)
0040
0041     rw = r.q__(iprs)
0042     tt = np.arange(*list(map(float,iprt[1:-1].split(":"))))
0043
0044
0045     _nt = "a.f.simtrace.shape[0]//rw.shape[0]"
0046     nt = eval(_nt)
0047
0048     _s = "a.f.simtrace.reshape(-1,nt,4,4)"
0049     s = eval(_s)
0050
0051     _e = "np.linalg.norm((s[:,:,1,:3] - s[:,-1,1,:3][:,np.newaxis]).reshape(-1,3),axis=1).reshape(-1,nt)[:,:-1]"
0052     e = eval(_e)
0053
0054     _d = "s[:,:,0,3][:,:-1]"
0055     d = eval(_d)
0056
0057     _dav = "np.average(d, axis=0)"
0058     dav = eval(_dav)
0059
0060     _eav = "np.average(e, axis=0)"
0061     eav = eval(_eav)
0062
0063     _emx = "np.max(e, axis=0 )"
0064     emx = eval(_emx)
0065
0066     _emn = "np.min(e, axis=0 )"
0067     emn = eval(_emn)
0068
0069     _emd = "np.median(e, axis=0 )"
0070     emd = eval(_emd)
0071
0072
0073     EXPR = r"""
0074     title
0075
0076     RFOLD
0077     r
0078
0079     AFOLD  # special simtrace precision run
0080     a
0081
0082     iprs   # OPTICKS_INPUT_PHOTON_RECORD_SLICE
0083
0084     iprt   # OPTICKS_INPUT_PHOTON_RECORD_TIME
0085     tt.shape
0086     tt
0087
0088     a.f.simtrace.shape
0089     rw.shape
0090     a.f.simtrace.shape[0]//rw.shape[0]     # nt
0091     a.f.simtrace.shape[0] % rw.shape[0]    # expect zero
0092
0093     _s
0094     #s
0095     s.shape
0096
0097     _nt
0098     nt
0099
0100     len(tt) == nt
0101
0102     # distance from intersect positions to last(and most precise) intersect position
0103     _e
0104     e
0105     e.shape
0106
0107     # distance from origin to intersect
0108     _d
0109     d
0110     d.shape
0111
0112     _dav   # average over clumps
0113     dav
0114
0115     _eav   # average over clumps
0116     eav
0117
0118     _emx
0119     emx
0120
0121     _emn
0122     emn
0123
0124     _emd
0125     emd
0126
0127     """
0128     for expr in list(map(str.strip, textwrap.dedent(EXPR).split("\n"))):
0129         print(expr)
0130         if expr == "" or expr[0] == "#": continue
0131         print(repr(eval(expr)))
0132     pass
0133
0134
0135
0136
0137     fig, ax = mp.subplots(figsize=SIZE/100.)
0138     fig.suptitle(title)
0139
0140
0141     ax.scatter( d.ravel(), e.ravel(), s=0.2 )
0142
0143     #ax.plot( dav, eav, label="avg" )
0144     #ax.plot( dav, emn, label="min" )
0145
0146     ax.plot( dav, emx, label="max" )
0147     ax.plot( dav, emd, label="median" )
0148
0149     ax.axhline(0.05, linestyle='--', label="EPSILON (0.05 mm)")
0150
0151     ax.set_ylabel("Intersect Error (mm)")
0152     ax.set_xlabel("Ray trace distance (mm)")
0153
0154     #ax.set_xscale('log')
0155     ax.set_yscale('log')
0156
0157     ax.set_ylim( 4e-5, 1.1 )
0158
0159     ax.legend()
0160
0161     fig.show()
0162
0163     if not FIGPATH is None:
0164         print("savefig to FIGPATH:%s" % FIGPATH)
0165         mp.savefig(FIGPATH)
0166     else:
0167         print("no-savefig as on FIGPATH")
0168     pass
0169
0170
0171
0172