qudarap/tests/QCerenkovIntegralTest.py

0001 #!/usr/bin/env python
0002 """
0003 QCerenkovIntegralTest.py
0004 ============================
0005
0006 ::
0007
0008     QCerenkovIntegralTest
0009     ipython -i tests/QCerenkovIntegralTest.py
0010
0011
0012 Getting very close match with symbolic integral results::
0013
0014
0015     In [10]: 1e12*(t.p_s2c[0,:,-1] - t.b_s2c[0,:,-1])
0016     Out[10]: array([ 0.   ,  0.006,  0.009,  0.007,  0.007,  0.014,  0.018,  0.014,  0.028,  0.014,  0.   , -0.014,  0.028,  0.028,  0.028,  0.028,  0.057,  0.057])
0017
0018
0019 """
0020 import os, logging, numpy as np
0021 log = logging.getLogger(__name__)
0022
0023 class QCerenkovIntegralTest(object):
0024     FOLD = os.path.expandvars("/tmp/$USER/opticks/QCerenkovIntegralTest")
0025     def __init__(self):
0026         pass
0027
0028     def load(self, qwn, pfx, nam):
0029         name = os.path.join(pfx, nam)
0030         path = name if name[0] == "/" else os.path.join(self.FOLD, name)
0031         log.debug("load %s " % path )
0032         a = np.load(path) if os.path.exists(path) else None
0033         if a is None:
0034             log.fatal("failed to load %s " % path )
0035         else:
0036             os.system("ls -l %s" % path )
0037         pass
0038         setattr(self, qwn, a )
0039
0040     def load_getS2Integral_UpperCut(self):
0041         pfx = "test_getS2Integral_UpperCut"
0042         self.load("a_s2c",  pfx, "s2c.npy")
0043         self.load("a_s2cn", pfx, "s2cn.npy")
0044
0045     def load_getS2Integral_SplitBin(self):
0046         pfx = "test_getS2Integral_SplitBin"
0047         self.load("b_s2c",  pfx, "s2c.npy")
0048         self.load("b_s2cn", pfx, "s2cn.npy")
0049         self.load("b_bis",  pfx, "bis.npy")
0050
0051     def load_piecewise(self):
0052         pfx = "/tmp/ana/piecewise"
0053         self.load("p_s2c", pfx, "scan.npy")
0054         self.load("p_bis", pfx, "bis.npy")
0055
0056     def plot(self, qwns, ii):
0057         fig, ax = plt.subplots(figsize=[12.8, 7.2])
0058         for qwn in qwns:
0059             a = getattr(self, qwn)
0060             #ken = 1 if qwn[0] == "b" else 0
0061             ken = 0   # reajjanged to put en_b into payload slot 0
0062             log.debug("qwn %s shape %s ken %d " % (qwn, str(a.shape), ken))
0063             assert len(a.shape) == 3
0064             for i in ii:
0065                 if qwn[0] == "b":
0066                     ax.scatter( a[i,:,ken], a[i,:,-1] , label="%s  %d" % (qwn, i) )
0067                 else:
0068                     ax.plot( a[i,:,ken], a[i,:,-1] , label="%s  %d" % (qwn, i) )
0069                 pass
0070             pass
0071             ax.legend()
0072         pass
0073         fig.show()
0074
0075
0076     def compare(self, ii, sa, sb, ):
0077         """
0078         * with mul=1 this is giving excellent agreement less than 1e-12
0079         * huh, after adding SUB handling discrep up to almost 0.1 photon
0080         """
0081         t = self
0082         a = getattr(t, sa)
0083         b = getattr(t, sb)
0084         log.info(" sa:%s a:%s sb:%s b:%s " % (sa, str(a.shape), sb, str(b.shape)))
0085         for i in ii:
0086             BetaInverse = t.b_bis[i]
0087             df = np.abs(a[i,:,-1] - b[i,:,-1])
0088             dfmax = df.max()
0089
0090             #print("BetaInverse : %10.4f  dfmax %10.4g  df %s " % (BetaInverse, dfmax, str(df))  )
0091             print("BetaInverse : %10.4f  dfmax %10.4g " % (BetaInverse, dfmax )  )
0092         pass
0093
0094
0095
0096 def plot_s2(t):
0097
0098     en = t.s2c[0][:,0]
0099     s2 = t.s2c[0][:,1]
0100
0101     fig, ax = plt.subplots(figsize=[12.8, 7.2])
0102     ax.plot( en, s2, label="s2" )
0103     ax.scatter( en, s2, label="s2" )
0104     ax.legend()
0105     fig.show()
0106
0107 if __name__ == '__main__':
0108     logging.basicConfig(level=logging.INFO)
0109     t = QCerenkovIntegralTest()
0110
0111     #ii = np.arange( 0, 1000, 100 )
0112     ii = np.arange(9)
0113     #ii = [0]
0114
0115     #t.load_getS2Integral_UpperCut()
0116     t.load_getS2Integral_SplitBin()
0117     t.load_piecewise()
0118
0119     #t.plot(["a_s2cn",], ii)
0120     #t.plot(["b_s2cn",], ii)
0121
0122     t.plot(["b_s2c","p_s2c"], ii)
0123
0124     #plot_s2(t)
0125
0126     assert np.all( t.p_bis == t.b_bis )
0127     t.compare(ii, "p_s2c", "b_s2c")
0128
0129
0130
0131