sysrap/tests/sreport_abcd.py

0001 #!/usr/bin/env python
0002
0003 import os, numpy as np
0004 from opticks.ana.fold import Fold
0005 from opticks.ana.npmeta import NPMeta
0006 from opticks.sysrap.tests.sreport import Substamp
0007 from opticks.sysrap.tests.sreport import RUN_META
0008
0009 MODE = 2
0010 PLOT =  os.environ.get("PLOT","?plot")
0011 COMMANDLINE = os.environ.get("COMMANDLINE", "")
0012
0013
0014 if MODE != 0:
0015     from opticks.ana.pvplt import *
0016 pass
0017
0018
0019 class ABCD_Substamp_ALL_Etime_vs_Photon(object):
0020     def __init__(self, *rr):
0021         title = RUN_META.ABCD_Title(*rr)
0022
0023         ratt = "abcd"
0024         assert len(rr) <= len(ratt)
0025         _qn0 = None
0026         for i,r in enumerate(rr):
0027             qf = "%sf" % ratt[i]
0028             qt = "%st" % ratt[i]
0029             qn = "%sn" % ratt[i]
0030             qh = "%sh" % ratt[i]
0031
0032             _qf = r.substamp.a
0033             _qt = Substamp.ETime(_qf)
0034             _qn = r.submeta_NumPhotonCollected.a
0035             _qh = Substamp.Subcount(_qf, "hit")
0036
0037             setattr(self, qf, _qf )
0038             setattr(self, qt, _qt )
0039             setattr(self, qn, _qn )
0040             setattr(self, qh, _qh )
0041
0042             if i == 0:
0043                 _qn0 = _qn
0044             else:
0045                 assert np.all( _qn0 == _qn )
0046             pass
0047         pass
0048         photon = _qn0[:,0]/1e6
0049
0050
0051         fontsize = 20
0052         YSCALE_ = ["log", "linear" ]
0053         YSCALE = os.environ.get("YSCALE", YSCALE_[1])
0054         assert YSCALE in YSCALE_
0055
0056
0057
0058         YMIN = float(os.environ.get("YMIN", "1e-2"))
0059         YMAX = float(os.environ.get("YMAX", "45"))
0060
0061         if MODE == 2:
0062             fig, axs = mpplt_plotter(nrows=1, ncols=1, label=title, equal=False)
0063             ax = axs[0]
0064
0065             #sli = slice(3,12)
0066             #sli = slice(None)
0067
0068             deg = 1  # 1:lin 2:par
0069
0070
0071             for i,r in enumerate(rr):
0072                 qt = "%st" % ratt[i]
0073                 _qt = getattr(self, qt)
0074
0075
0076                 # from A_SLI, B_SLI, if present otherwise from SLI, default "3:12"
0077
0078                 Q_SLI = "%s_SLI" % ratt[i].upper()
0079                 if Q_SLI in os.environ:
0080                     _SLI = os.environ.get(Q_SLI)
0081                 else:
0082                     SLI_DEF="3:12"
0083                     _SLI = os.environ.get("SLI", SLI_DEF)
0084                 pass
0085                 SLI = list(map(int, _SLI.split(":")))
0086                 sli = slice(*SLI)
0087
0088
0089
0090                 q_fit = "%s_fit" % ratt[i]
0091                 r_fit = np.poly1d(np.polyfit(photon[sli], _qt[sli], deg))
0092                 r_fit_label = "%s : linefit( slope %10.3f  intercept %10.3f )" % (r.symbol, r_fit.coef[0], r_fit.coef[1])
0093
0094                 setattr(self, q_fit, r_fit )
0095
0096                 r_GPU_label = RUN_META.GPUMeta(r)
0097                 r_RNG_label = RUN_META.QSim__RNGLabel(r)
0098                 r_label = "%s (%s)" % (r_GPU_label, r_RNG_label)
0099
0100                 ax.scatter( photon, _qt, label=r_label )
0101                 ax.plot( photon[sli], r_fit(photon[sli]), label=r_fit_label )
0102             pass
0103             ax.set_xlim( -5, 105 );
0104             ax.set_ylim( YMIN, YMAX );  # 50*200 = 1e4
0105             ax.set_yscale(YSCALE)
0106             ax.set_ylabel("Event time (seconds)", fontsize=fontsize )
0107             ax.set_xlabel("Number of Photons (Millions)", fontsize=fontsize )
0108             ax.legend()
0109             ax.legend()
0110             fig.show()
0111         pass
0112
0113
0114
0115
0116
0117 if __name__ == '__main__':
0118
0119     print("[sreport_abcd.py:PLOT[%s]" % PLOT )
0120
0121     REPS="ABCD"
0122     rr = []
0123     for i, Q in enumerate(REPS):
0124
0125         q = os.environ.get(Q, None)
0126         if q is None:
0127             continue
0128         pass
0129         QRF = "$%s_SREPORT_FOLD" % Q
0130
0131         print("[sreport_abcd.py:r = Fold.Load(%s, symbol=%s)" % (QRF, q) )
0132         r = Fold.Load(QRF, symbol=q )
0133         print("]sreport_abcd.py:r = Fold.Load(%s, symbol=%s)" % (QRF, q) )
0134
0135         print("[sreport_abcd.py:repr(r)" )
0136         print(repr(r))
0137         print("]sreport_abcd.py:repr(r)" )
0138         rr.append(r)
0139     pass
0140     print("]sreport_ab.py:PLOT[%s]" % PLOT )
0141
0142     if PLOT.startswith("ABCD_Substamp_ALL_Etime_vs_Photon"):
0143         abcd = ABCD_Substamp_ALL_Etime_vs_Photon(*rr)
0144     pass
0145
0146
0147
0148
0149