eic-opticks/sysrap/sevt_tt.py

0001 #!/usr/bin/env python
0002 """
0003 sevt_tt.py (formerly u4/tests/U4SimulateTest_tt.py)
0004 ======================================================
0005
0006 Python script (not module) usage from two bash scripts
0007 --------------------------------------------------------
0008
0009 ::
0010
0011     PLOT=STAMP STAMP_ANNO=1 STAMP_LEGEND=1 ~/opticks/u4/tests/tt.sh
0012
0013     PLOT=STAMP STAMP_TT=90000,5000 STAMP_ANNO=1 ~/j/ntds/ntds.sh tt
0014
0015     PLOT=PHO_AVG ~/j/ntds/ntds.sh tt
0016
0017
0018 PLOT envvar selects what to plot
0019 ----------------------------------
0020
0021 PLOT=STAMP
0022    time stamp illustration
0023
0024    STAMP_TT=200000,5000
0025    STAMP_TT=200k,5k
0026       control time window of STAMP plot in microseconds [us]
0027
0028    STAMP_ANNO=1
0029       enable photon index and history annotation
0030
0031 PLOT=PHO_AVG
0032    average beginPhoton->endPhoton CPU time vs photon point count
0033
0034 PLOT=PHO_HIS
0035    histogram of beginPhoton->endPhoton CPU times
0036
0037    PLOT=PHO_HIS ~/j/ntds/ntds.sh tt
0038
0039 PLOT=PHO_N
0040    A, B photon step point counts (fakes must be skipped in A to get a match)
0041
0042    PLOT=PHO_N ~/j/ntds/ntds.sh tt
0043
0044
0045 Dev Notes
0046 ------------
0047
0048 HMM: for reusability, its tempting to
0049 just relocate this into sysrap/tests/tt.py
0050 as a runnable python script which is used from
0051 wherever with bash level envvars controlling the events
0052 to be loaded rather than doing python dev to
0053 bring in functionality at python level.
0054
0055 When is this bash in control approach appropriate
0056 as opposed to moving into sevt.py ?
0057
0058 The level of generality determines what is appropriate
0059 If the functionality really is general and likely
0060 to be usable from all over the place then it belongs
0061 into sevt.py with python control and flexibility of use.
0062
0063 However if only likely to use functionality from a few places
0064 pointing at differnt events that doing
0065 at bash level seems appropriate and avoids python
0066 module development.
0067
0068 ::
0069
0070     u4t
0071     ./U4SimulateTest.sh tt
0072
0073     PLOT=STAMP   ~/opticks/u4/tests/tt.sh
0074     PLOT=PHO_HIS ~/opticks/u4/tests/tt.sh
0075     PLOT=PHO_AVG ~/opticks/u4/tests/tt.sh
0076
0077
0078 """
0079 import os, numpy as np, textwrap, logging
0080 log = logging.getLogger(__name__)
0081 from opticks.ana.fold import Fold
0082 from opticks.ana.p import *
0083 from opticks.ana.eget import efloatarray_
0084 from opticks.sysrap.sevt import SEvt
0085
0086
0087 SCRIPT = "./U4SimulateTest.sh tt"
0088 os.environ["SCRIPT"] = SCRIPT
0089 LABEL = os.environ.get("LABEL", "U4SimulateTest_tt.py" )
0090
0091
0092 N = int(os.environ.get("N", "-1"))
0093 MODE =  int(os.environ.get("MODE", "2"))
0094 assert MODE in [0,2,-2,3,-3]
0095 FLIP = int(os.environ.get("FLIP", "0")) == 1
0096 TIGHT = int(os.environ.get("TIGHT", "0")) == 1
0097 STAMP_TT = efloatarray_("STAMP_TT", "162307,3000") # "t0,dt"  microseconds 1M=1s
0098
0099 PLOT = os.environ.get("PLOT", None)
0100 NEVT = int(os.environ.get("NEVT", 0))  # when NEVT>0 SEvt.LoadConcat loads and concatenates the SEvt
0101
0102 ENVOUT = os.environ.get("ENVOUT", None)
0103 CAP_STEM = PLOT
0104 CAP_BASE = None # set below to a.f.base or b.f.base
0105
0106 if MODE != 0:
0107     from opticks.ana.pvplt import *
0108     WM = mp.pyplot.get_current_fig_manager()
0109 else:
0110     WM = None
0111 pass
0112
0113
0114
0115 def GetAnnotation(title, EXPRS_, syms):
0116     EXPRS = list(filter(None, textwrap.dedent(EXPRS_).split("\n")))
0117
0118     rlines = []
0119     rlines.append(title)
0120
0121     for i,sym in enumerate(syms):
0122         for expr_ in EXPRS:
0123             label = expr_[expr_.find("#"):] if expr_.find("#") > -1 else ""
0124             expr_ = expr_.split("#")[0].strip()
0125             expr = expr_ % locals()
0126             print(expr)
0127             try:
0128                 val = eval(expr)
0129             except ValueError:
0130                 log.fatal("FAILED EVAL : %s " % expr )
0131                 val = "?"
0132             pass
0133             fmt_ = {}
0134             fmt_["str"] = "%30s : %s : %s "
0135             fmt_["int"] = "%30s : %8d : %s "
0136             fmt_["float"] = "%30s : %8.3f : %s "
0137             fmt_[""] = fmt_["float"]
0138
0139             if type(val) is list: val = " ".join(val)
0140             fmt = fmt_.get(type(val).__name__, fmt_["float"] )
0141             rlines.append(fmt % ( expr, val, label  ))
0142         pass
0143         rlines.append("")
0144     pass
0145     ranno = "\n".join(rlines)
0146     return ranno
0147
0148
0149
0150 if __name__ == '__main__':
0151     logging.basicConfig(level=logging.INFO)
0152
0153     print("AFOLD:%s" % os.environ.get("AFOLD", "-"))
0154     print("BFOLD:%s" % os.environ.get("BFOLD", "-"))
0155     print("N:%d " % N )
0156
0157     if N == -1:
0158         a = SEvt.Load("$AFOLD",symbol="a", NEVT=NEVT)
0159         b = SEvt.Load("$BFOLD",symbol="b", NEVT=NEVT)
0160         syms = ['a','b']
0161         evts = [a,b]
0162     elif N == 0:
0163         a = SEvt.Load("$AFOLD",symbol="a", NEVT=NEVT)
0164         b = None
0165         syms = ['a']
0166         evts = [a,]
0167     elif N == 1:
0168         a = None
0169         b = SEvt.Load("$BFOLD",symbol="b", NEVT=NEVT)
0170         syms = ['b']
0171         evts = [b,]
0172     else:
0173         assert(0)
0174     pass
0175
0176     ## CAP_BASE is passed via ENVOUT to the invoking bash script
0177     ## to control where the figs folder with screen captures should be
0178     if not a is None:
0179         CAP_BASE = a.f.base
0180     elif not b is None:
0181         CAP_BASE = b.f.base
0182     pass
0183
0184     # concatenated SEvt have list of bases
0185     if type(CAP_BASE) is list: CAP_BASE = CAP_BASE[0]
0186
0187     if not a is None:print(repr(a))
0188     if not b is None:print(repr(b))
0189
0190
0191     if PLOT == "EXPRS":
0192         EXPRS = r"""
0193         %(sym)s
0194         %(sym)s.symbol
0195         w
0196         %(sym)s.q[w]
0197         %(sym)s.n[w]
0198         np.diff(%(sym)s.tt[w])
0199         np.c_[%(sym)s.t[w].view("datetime64[us]")]
0200         """
0201
0202         for sym in syms:
0203
0204             n = eval("%(sym)s.n" % locals() )
0205             ww = np.where( n > 24)[0]
0206
0207             for w in ww:
0208                 for e_ in textwrap.dedent(EXPRS).split("\n"):
0209                     e = e_ % locals()
0210                     print(e)
0211                     if len(e) == 0 or e[0] == "#": continue
0212                     print(eval(e))
0213                 pass
0214             pass
0215         pass
0216     pass
0217
0218
0219     # needs: syms, a, b
0220
0221     if PLOT == "PHO_HIS":
0222         msg = "histogram of beginPhoton->endPhoton durations"
0223         label = "%s : %s" % (PLOT, msg)
0224         fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False)
0225         ax = axs[0]
0226         e_ = "%(sym)s.s1 - %(sym)s.s0  # PHO_HIS "
0227
0228         #bins = np.logspace( np.log10(0.1),np.log10(500.0), 25 )
0229         bins = np.linspace( 0, 500, 50 )
0230
0231         h = {}
0232         for sym in syms:
0233             e = e_ % locals()
0234             h[sym] = np.histogram( eval(e) , bins=bins )
0235             ax.plot( h[sym][1][:-1], h[sym][0], label=e  )
0236         pass
0237         ax.legend()
0238         fig.show()
0239     pass
0240
0241
0242     if PLOT == "TABLE":
0243
0244         HEAD = "np.c_["
0245         TAIL = "]"
0246         FIELDS = list(filter(None,textwrap.dedent(r"""
0247         %(sym)s.s1[w] -  %(sym)s.s0[w]                          # beginPhoton->endPhoton
0248         %(sym)s.t[w,0] - %(sym)s.s0[w]                          # beginPhoton->firstPoint
0249         %(sym)s.t[w,%(n_1)s] - %(sym)s.t[w,0]                   # firstPoint->lastPoint
0250         %(sym)s.s1[w] - %(sym)s.t[w,%(n_1)s]                    # lastPoint->endPhoton
0251         """).split("\n")))
0252
0253         LABELS = list(map(lambda _:_[_.find("#")+1:], FIELDS))
0254         FIELDS = list(map(lambda _:_[:_.find("#")].strip(), FIELDS))
0255
0256         print("compare first and last point stamp range with beginPhoton endPhoton range")
0257
0258         nn = np.arange(2,33)
0259         for n in nn:
0260             n_1 = int(n - 1)
0261             for sym in syms:
0262                 w_ = "np.where( %(sym)s.n == %(n)s )[0]" % locals()
0263                 w = eval(w_)
0264                 expr__ = "".join([HEAD, ",".join(FIELDS), TAIL ])
0265                 expr_ = expr__ % locals()
0266                 expr = eval(expr_)
0267                 label = " ".join(LABELS)
0268                 if len(expr) == 0: continue
0269
0270                 print(expr_)
0271                 print(w_)
0272                 print(label)
0273                 print(expr)
0274             pass
0275         pass
0276     pass
0277
0278
0279     print(" average photonBegin->photonEnd for different step counts ")
0280     ssa = np.zeros((33,len(syms),4), dtype=np.float64 )
0281     for n in range(33):
0282         for i,sym in enumerate(syms):
0283             w_ = "np.where( %(sym)s.n == %(n)s )[0]" % locals()
0284             w = eval(w_)
0285             nw = len(w)
0286             ss_ = "%(sym)s.ss[%(sym)s.n == %(n)s ]" % locals()
0287             ss = eval(ss_)
0288
0289             ssa[n,i,0] = n
0290             ssa[n,i,1] = nw
0291
0292             if len(ss) == 0: continue
0293
0294             ssa[n,i,2] = np.average(ss)
0295             ssa[n,i,3] = np.std(ss, ddof=1)/np.sqrt(np.size(ss))
0296         pass
0297     pass
0298
0299
0300     expr_ = "np.c_[ssa.reshape(-1,8)]"
0301     print(expr_)
0302     print(eval(expr_))
0303
0304     if PLOT == "PHO_AVG":
0305         EXPRS_ = r"""
0306         %(sym)s.f.baselabel                        #
0307         %(sym)s.FAKES_SKIP                         # from run_meta
0308         np.diff(%(sym)s.rr)[0]/1e6                 # Run
0309         %(sym)s.ee[-1]/1e6                         # Evt
0310         np.sum(%(sym)s.ss)/1e6                     # Pho
0311         np.sum(%(sym)s.ss)/%(sym)s.ee[-1]          # Pho/Evt
0312         """
0313         title = "sevt_tt.PHO_AVG.table NEVT:%(NEVT)d ## " % locals()
0314         ranno = GetAnnotation(title, EXPRS_, syms )
0315         os.environ["RHSANNO"] = "\n".join(filter(lambda line:line.find("##")==-1, ranno.split("\n")))
0316         os.environ["RHSANNO_POS"]="0.45,-0.02"
0317
0318         cut = 10
0319         labelx = "stat_cut:%(cut)d NEVT:%(NEVT)d " % locals()
0320         label = "PHO_AVG : average beginPhoton->endPhoton duration [us] vs (step point count - 2)   # " + labelx
0321         print(label)
0322
0323         fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False)
0324         ax = axs[0]
0325
0326         #sl = slice(2,16)
0327         sl = slice(None)
0328         style = {'a':"ro-", 'b':"bo-" }
0329         fitstyle = {'a':"r:", 'b':"b:" }
0330
0331
0332         mc = np.zeros( (len(syms), 2), dtype=np.float64 )
0333
0334
0335         for i,sym in enumerate(syms):
0336             x_ = ssa[sl,i,0] - 2   # step point count - 2 (so starts from 0)
0337             n_ = ssa[sl,i,1]
0338             y_ = ssa[sl,i,2]
0339             ye_ = ssa[sl,i,3]
0340
0341             w = np.where( n_ > cut)[0]
0342
0343             x = x_[w]
0344             y = y_[w]
0345             ye = ye_[w]
0346
0347             label = "%(sym)s : avg(ss) vs n-2  " % locals()
0348
0349             #ax.plot( x, y , style[sym], label=label )
0350
0351             ax.errorbar( x, y,yerr=ye, fmt=style[sym], label=label  )
0352
0353
0354
0355             wf = np.where( n_ > 10)[0]  ## only fit points with some stats
0356             x = x_[wf]
0357             y = y_[wf]
0358
0359             A = np.vstack([x, np.ones(len(x))]).T
0360             m, c = np.linalg.lstsq(A, y, rcond=None)[0]
0361             mc[i,0] = m
0362             mc[i,1] = c
0363             ax.plot(x, m*x + c, fitstyle[sym], label='Fit: m %10.3f c %10.3f ' % (m,c) )
0364         pass
0365         ax.legend()
0366         fig.show()
0367
0368         print(label)
0369         expr = "np.c_[mc]"
0370         print(expr)
0371         print(eval(expr))
0372
0373         print("ranno")
0374         print(ranno)
0375     pass
0376
0377
0378     if PLOT == "PHO_N":
0379
0380         EXPRS_ = r"""
0381         %(sym)s.f.baselabel                        #
0382         %(sym)s.FAKES_SKIP                         # from run_meta
0383         """
0384         title = "sevt_tt.PHO_N.table NEVT:%(NEVT)d ## " % locals()
0385         ranno0 = GetAnnotation(title, EXPRS_, syms )
0386         ranno = "\n".join(filter(lambda line:line.find("##")==-1, ranno0.split("\n")))
0387         os.environ["RHSANNO"] = ranno
0388         os.environ["RHSANNO_POS"]="0.1,-0.02"
0389
0390         a_n, a_nc = np.unique(a.n, return_counts=True)
0391         b_n, b_nc = np.unique(b.n, return_counts=True)
0392
0393         #msg = "consistent when fakes are skipped"
0394         msg = "much more in A when fakes not skipped : NEVT:%(NEVT)d" % locals()
0395         label = "PHO_N : A(N=0),B(N=1) photon step counts : %s " % msg
0396         fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False)
0397         ax = axs[0]
0398         ax.set_yscale("log")
0399         #ax.plot( a_n, a_nc, "ro", label="A"  )
0400         #ax.plot( b_n, b_nc, "bo", label="B" )
0401
0402         ax.errorbar( a_n, a_nc,yerr=np.sqrt(a_nc.astype(np.float64)),fmt="ro-", label="A"  )
0403         ax.errorbar( b_n, b_nc,yerr=np.sqrt(b_nc.astype(np.float64)),fmt="bo-", label="B" )
0404
0405         ax.legend()
0406         fig.show()
0407     pass
0408
0409
0410
0411
0412     if PLOT == "PHO_SCAT":
0413         msg = "scatter plot of beginPhoton->endPhoton CPU time vs step point count"
0414         label = "%s : %s" % (PLOT, msg)
0415
0416         fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False)
0417         ax = axs[0]
0418
0419         ax.set_ylim(0,1000)
0420
0421         for i,sym in enumerate(syms):
0422             x_ = "%(sym)s.n.astype(np.float64)" % locals()
0423             y_ = "%(sym)s.ss" % locals()
0424             x = eval(x_)
0425             y = eval(y_)
0426             if i == 1: x += 0.2
0427
0428             ax.scatter( x[x<16], y[x<16], s=2 )
0429         pass
0430         ax.legend()
0431         fig.show()
0432     pass
0433
0434
0435     """
0436     STAMP NOTES
0437
0438     HMM: how to annotate the stamp plot with photon indices ?
0439
0440     For comparibility between A and B it makes more
0441     sense to select based on time stamps rather than photon indices.
0442
0443     To find time range with big bouncers::
0444
0445          a.s0[np.where(a.n > 20)]-a.ee[0]
0446
0447     """
0448     if PLOT == "STAMP":
0449
0450
0451         t0 = STAMP_TT[0]
0452         t1 = STAMP_TT[0]+STAMP_TT[1]
0453         stt = os.environ.get("STAMP_TT","-")
0454
0455         label = "A:lhs, B:rhs :  STAMP_TT=%s #  (t0,t1):(%d,%d)  " % (stt,t0,t1)
0456
0457         #os.environ["SUBTITLE"] = subtitle
0458         #os.environ["THIRDLINE"] = subtitle
0459
0460         fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False)
0461         ax = axs[0]
0462
0463         #sl = slice(-100,None,1)  # last 100 photons, which are actually processed first
0464         #sl = slice(1040,1550,1)
0465         sl = slice(None)
0466
0467         s0_ = {}
0468         s1_ = {}
0469         tt_ = {}
0470         wt_ = {}
0471         h0_ = {}
0472         h1_ = {}
0473         i0_ = {}
0474         i1_ = {}
0475
0476         s0 = {}
0477         s1 = {}
0478         tt = {}
0479         wt = {}
0480         h0 = {}
0481         h1 = {}
0482         i0 = {}
0483         i1 = {}
0484
0485         sz = 0.01
0486         xx = {'a':[-0.5, -sz], 'b':[sz,0.5] }
0487         zz = {'a':[-0.2, -sz], 'b':[sz,0.2] }
0488         yy = {'a':-0.5, 'b':0.20 }
0489
0490         qwns = "s0 s1 t h0 h1 i0 i1".split()
0491         q_expr = "%(sym)s.%(q)s[sl] - %(sym)s.ee[0]"
0492
0493         print(" (t0,t1)  (%(t0)d,%(t1)d) " % locals() )
0494
0495         for i, sym in enumerate(syms):
0496             print("sym:%s" % sym)
0497             for j in range(len(qwns)):
0498                 q = qwns[j]
0499                 expr = q_expr % locals()
0500                 t = eval(expr)
0501                 w = eval("np.where(np.logical_and(t > t0, t < t1 ))[0]")
0502                 ws = str(w.shape)
0503                 qmn = t.min()
0504                 qmx = t.max()
0505                 fmt = " %(q)2s  %(expr)22s : (%(qmn)8d,%(qmx)8d) : ws:%(ws)s "
0506                 print( fmt % locals() )
0507             pass
0508         pass
0509
0510         print(" when all the ws are 0: adjust the time range to find some stamps")
0511
0512
0513         for i, sym in enumerate(syms):
0514             r0 = eval("%(sym)s.rr[0]" % locals())
0515             #if i == 1:continue
0516             print( "sym:%s r0:%10.3f " % (sym, r0))
0517
0518             ## where selecting on times greater than zero messes up the indices,
0519             ## so instead just rely on exclusion of crazy times for unfilled cases
0520
0521             cols = "r b g c m c m".split()
0522             assert len(cols) == len(qwns)
0523
0524             for j in range(len(qwns)):
0525                 q = qwns[j]
0526                 expr = q_expr % locals()
0527                 t = eval(expr)
0528                 w = eval("np.where(np.logical_and(t > t0, t < t1 ))[0]")
0529                 # photon indices of times within time window
0530                 kk = list(range(len(w)))
0531                 xmin,xmax = zz[sym] if q == "t" else xx[sym]
0532
0533                 ax.hlines( t[np.logical_and(t > t0, t < t1 )], xmin, xmax, cols[j], label=expr )
0534
0535                 if "STAMP_ANNO" in os.environ:
0536                     for k in kk:
0537                         idx = w[k]
0538                         anno = "(%s) %s : %d : " % (q, sym.upper(), idx)
0539                         if q == "s0":
0540                             his = eval("%(sym)s.q[%(idx)s][0].decode(\"utf-8\").strip()" % locals())
0541                             anno +=  his
0542                         elif q == "h0":
0543                             hc = eval("%(sym)s.hc[%(idx)s]" % locals())
0544                             anno += " hc:%d " % (hc)
0545                         elif q == "i0":
0546                             ic = eval("%(sym)s.ic[%(idx)s]" % locals())
0547                             hi0 = eval("%(sym)s.hi0[%(idx)s]" % locals())
0548                             anno += "ic:%d hi0:%d "% (ic, hi0)
0549                         else:
0550                             anno = None
0551                         pass
0552                         if not anno is None:
0553                             ax.text( yy[sym], t[idx], anno )
0554                         pass
0555                     pass
0556                 pass
0557             pass
0558         pass
0559
0560         if "STAMP_LEGEND" in os.environ:
0561             ax.legend()
0562         pass
0563         fig.show()
0564     pass
0565     if not ENVOUT is None and not CAP_STEM is None and not CAP_BASE is None:
0566         envout = "\n".join([
0567                        "export ENVOUT_PATH=%s" % ENVOUT,
0568                        "export ENVOUT_CAP_STEM=%s" % CAP_STEM,
0569                        "export ENVOUT_CAP_BASE=%s" % CAP_BASE,
0570                        ""
0571                        ])
0572         open(ENVOUT, "w").write(envout)
0573         print(envout)
0574     pass
0575
0576