qudarap/tests/QCKTest.py

0001 #!/usr/bin/env python
0002 """
0003 QCKTest.py
0004 ==================
0005
0006 ::
0007
0008     QCerenkovIntegralTest         # create and persist ICDF
0009     QCKTest                       # load ICDF and create energy samples, persisting into ICDF dirs
0010     ipython -i tests/QCKTest.py   # analysis/plotting comparing rejection and lookup sampling
0011
0012
0013 Hmm largest c2 for different BetaInverse always at 7.6eV
0014 just to left of rindex peak.
0015 Possibly there is a one bin shifted issue, that is showing up
0016 the most in the region where rindex is changing fastest.
0017
0018 Perhaps could check this with an artifical rindex pattern,
0019 such as a step function.
0020
0021 Actually just setting BetaInverse to 1.792 just less than rmx 1.793
0022 is informative as then there is only a very small range
0023 of possible energies.
0024
0025 Hmm: generating millions of photons just there is a kinda
0026 extreme test, as in reality will be less than 1.
0027
0028 Hmm maybe should exclude BetaInverse where the average number
0029 of photons is less than 1
0030
0031 ::
0032
0033     In [18]: np.c_[t.s2c[-120:-100,-1],t.bis[-120:-100]]
0034     Out[18]:
0035     array([[8.45 , 1.113, 1.746],
0036            [8.448, 1.1  , 1.746],
0037            [8.446, 1.087, 1.747],
0038
0039
0040 See also::
0041
0042     ana/rindex.py
0043     ana/ckn.py
0044
0045 """
0046 import os, logging, numpy as np
0047
0048 from opticks.ana.nbase import chi2
0049 from opticks.ana.edges import divide_bins
0050 from opticks.ana.rsttable import RSTTable
0051
0052 log = logging.getLogger(__name__)
0053
0054 def isfloat(value):
0055     try:
0056         float(value)
0057         return True
0058     except ValueError:
0059         return False
0060     pass
0061
0062 class QCKTest(object):
0063     FOLD = os.path.expandvars("/tmp/$USER/opticks/QCerenkovIntegralTest")
0064
0065     def __init__(self, approach="UpperCut", use_icdf=False):
0066         """
0067         :param approach: choosing the ICDF folder
0068         :param use_icdf: choosing between energy samples ending _s2cn.npy and _icdf.npy
0069         """
0070         assert approach in ["UpperCut", "SplitBin"]
0071         self.approach = approach
0072         self.use_icdf = use_icdf
0073
0074         base_path = os.path.join(self.FOLD, "test_makeICDF_%s" % approach)
0075         if not os.path.exists(base_path):
0076             log.fatal("base_path %s does not exist" % base_path)
0077             assert 0
0078         pass
0079         names = os.listdir(base_path)
0080         log.info("loading from base_path %s " % base_path)
0081         for name in filter(lambda _:_.endswith(".npy"), names):
0082             path = os.path.join(base_path, name)
0083             stem = name[:-4]
0084             a = np.load(path)
0085             print( " t.%5s  %s " % (stem, str(a.shape)))
0086             setattr(self, stem, a )
0087         pass
0088         sample_base = os.path.join(base_path, "QCKTest")
0089         self.sample_base = sample_base
0090     pass
0091
0092     def bislist(self):
0093         names = sorted(os.listdir(os.path.expandvars(self.sample_base)))
0094         names = list(filter(lambda n:isfloat(n), names))
0095         print(names)
0096         bis = list(map(float, names))
0097         return bis
0098
0099     def s2cn_plot(self, istep):
0100         """
0101         :param ii: list of first dimension indices, corresponding to BetaInverse values
0102         """
0103         s2cn = self.s2cn
0104         bis = self.bis
0105         sample_base = self.sample_base
0106
0107         assert len(s2cn) == len(bis)
0108
0109         ii = np.arange( 0,len(s2cn), istep )
0110
0111         log.info("s2cn %s istep %s ii %s " % (str(s2cn.shape),istep, str(ii) ))
0112
0113         sbis = "%6.4f ... %6.4f " % (bis[0], bis[-1])
0114
0115         title_ = "QCKTest.py : s2cn_plot : s2cn.shape %s bis.shape %s bis %s istep %d " % (str(s2cn.shape), str(bis.shape), sbis, istep)
0116         desc_ = "JUNO LS : Cerenkov S2 integral CDF for sample of BetaInverse values"
0117
0118         title = "\n".join([title_, desc_])
0119         fig, ax = plt.subplots(figsize=[12.8, 7.2])
0120         fig.suptitle(title)
0121         for i in ii:
0122             label = "%6.4f" % bis[i]
0123             #ulabel = label if i % 100 == 0 else None
0124             ulabel = label
0125             ax.plot( s2cn[i,:,0], s2cn[i,:,-1] , label=ulabel )
0126         pass
0127         ax.legend()
0128         fig.show()
0129         figpath = os.path.join(sample_base, "s2cn_plot.png")
0130         fig.savefig(figpath)
0131         log.info(figpath)
0132
0133
0134     def one_s2cn_plot(self, BetaInverse ):
0135         s2cn = self.s2cn
0136         ibi = self.getBetaInverseIndex(BetaInverse)
0137         title_ = "QCKTest.py : one_s2cn_plot BetaInverse %6.4f  ibi %d s2cn[ibi] %s " % (BetaInverse, ibi, str(s2cn[ibi].shape))
0138         desc_ = " cdf (normalized s2 integral) for single BetaInverse "
0139
0140         title = "\n".join([title_, desc_]) ;
0141         fig, ax = plt.subplots(figsize=[12.8, 7.2])
0142         fig.suptitle(title)
0143
0144         ax.plot( s2cn[ibi,:,0], s2cn[ibi,:,-1] , label="s2cn[%d]" % ibi )
0145         ax.legend()
0146         fig.show()
0147
0148
0149     def getBetaInverseIndex(self, BetaInverse):
0150         bis = self.bis
0151         ibi = np.abs(bis - BetaInverse).argmin()
0152         return ibi
0153
0154     def rindex_plot(self):
0155         ri = self.rindex
0156         c2 = self.c2
0157         c2poppy = self.c2poppy
0158         bi = self.bi
0159         edges = self.edges
0160         c2riscale = self.c2riscale
0161
0162         title = "\n".join(["QCKTest.py : rindex_plot"]) ;
0163         fig, ax = plt.subplots(figsize=[12.8, 7.2])
0164         fig.suptitle(title)
0165
0166         ax.scatter( ri[:,0], ri[:,1], label="ri" )
0167         ax.plot(    ri[:,0], ri[:,1], label="ri" )
0168         ax.plot(    edges[:-1], c2*c2riscale,    label="c2", drawstyle="steps-post" )
0169         ax.plot( [ri[0,0], ri[-1,0]], [bi, bi], label="bi %6.4f " % bi )
0170
0171         bi0 = 1.75
0172         ax.plot( [ri[0,0], ri[-1,0]], [bi0, bi0], label="bi0 %6.4f " % bi0 )
0173
0174         ylim = ax.get_ylim()
0175
0176         for i in c2poppy:
0177             ax.plot( [edges[i], edges[i]], ylim , label="edge %d " % i, linestyle="dotted" )
0178             ax.plot( [edges[i+1], edges[i+1]], ylim , label="edge+1 %d " % (i+1), linestyle="dotted" )
0179         pass
0180
0181         ax.legend()
0182         fig.show()
0183
0184     def en_load(self, bi):
0185         bi_base = os.path.expandvars("%s/%6.4f" % (self.sample_base, bi) )
0186
0187         use_icdf = self.use_icdf
0188         ext = ["s2cn","icdf"][int(use_icdf)] ;
0189         log.info("load from bi_base %s ext %s " % (bi_base, ext))
0190
0191         el = np.load(os.path.join(bi_base,"test_energy_lookup_many_%s.npy" % ext ))
0192         es = np.load(os.path.join(bi_base,"test_energy_sample_many.npy"))
0193
0194         tl = np.load(os.path.join(bi_base,"test_energy_lookup_many_tt.npy"))
0195         ts = np.load(os.path.join(bi_base,"test_energy_sample_many_tt.npy"))
0196
0197         self.bi_base = bi_base
0198         self.el = el
0199         self.es = es
0200         self.tl = tl
0201         self.ts = ts
0202
0203
0204     def check_s2c_monotonic(self):
0205         s2c = self.s2c
0206         for i in range(len(s2c)):
0207             w = np.where( np.diff(s2c[i,:,2]) < 0 )[0]
0208             print(" %5d : %s " % (i, str(w)))
0209         pass
0210
0211     def en_compare(self, bi, num_edges=101):
0212         """
0213         Compare the energy samples created by QCKTest for a single BetaInverse
0214         """
0215         ri = self.rindex
0216         el = self.el
0217         es = self.es
0218         s2cn = self.s2cn
0219         avph = self.avph
0220         s2c  = self.s2c
0221
0222         ibi = self.getBetaInverseIndex(bi)
0223
0224         approach = self.approach
0225         if approach == "UpperCut":  # see QCerenkov::getS2Integral_UpperCut
0226             en_slot = 0
0227             s2_slot = 1
0228             cdf_slot = 2
0229             emn = s2cn[ibi, 0,en_slot]
0230             emx = s2cn[ibi,-1,en_slot]
0231             avp = s2c[ibi, -1,cdf_slot]
0232         elif approach == "SplitBin":  # see QCerenkov::getS2Integral_SplitBin
0233             en_slot = 0   # en_b
0234             s2_slot = 5   # s2_b
0235             cdf_slot = 7  # s2integral
0236             emn = avph[ibi, 1]
0237             emx = avph[ibi, 2]
0238             avp = avph[ibi, 3]
0239         else:
0240             assert 0, "unknown approach %s " % approach
0241         pass
0242
0243         self.en_slot = en_slot
0244         self.s2_slot = s2_slot
0245         self.cdf_slot = cdf_slot
0246         self.emn = emn
0247         self.emx = emx
0248         self.avp = avp
0249
0250         edom = emx - emn
0251         edif = edom/(num_edges-1)
0252         edges0 = np.linspace( emn, emx, num_edges )    # across Cerenkov permissable range
0253         edges = np.linspace( emn-edif, emx+edif, num_edges + 2 )   # push out with extra bins either side
0254
0255         #edges = np.linspace(1.55,15.5,100)  # including rightmost
0256         #edges = np.linspace(1.55,15.5,200)  # including rightmost
0257         #edges = divide_bins( ri[:,0], mul=4 )
0258
0259         hl = np.histogram( el, bins=edges )
0260         hs = np.histogram( es, bins=edges )
0261
0262         c2, c2n, c2c = chi2( hl[0], hs[0] )
0263         ndf = max(c2n - 1, 1)
0264         c2sum = c2.sum()
0265         c2p = c2sum/ndf
0266         c2label = "chi2/ndf %4.2f [%d] %.2f " % (c2p, ndf, c2sum)
0267
0268         c2amx = c2.argmax()
0269         rimax = ri[:,1].max()
0270         c2max = c2.max()
0271         c2riscale = rimax/c2max
0272         c2poppy = np.where( c2 > c2max/3. )[0]
0273
0274         hmax = max(hl[0].max(), hs[0].max())
0275         c2hscale = hmax/c2max
0276
0277         cf = " c2max:%4.2f c2amx:%d c2[c2amx] %4.2f edges[c2amx] %5.3f edges[c2amx+1] %5.3f " % (c2max, c2amx, c2[c2amx], edges[c2amx], edges[c2amx+1] )
0278
0279         print("cf", cf)
0280
0281         #print("c2", c2)
0282         print("c2n", c2n)
0283         print("c2c", c2c)
0284
0285         qq = "hl hs c2 c2label c2n c2c c2riscale c2hscale hmax edges c2max c2poppy cf bi ibi"
0286         for q in qq.split():
0287             globals()[q] = locals()[q]
0288             setattr(self, q, locals()[q] )
0289         pass
0290
0291         t = self
0292         print("np.c_[t.c2, t.hs[0], t.hl[0]][t.c2 > 0]")
0293         print(np.c_[t.c2, t.hs[0], t.hl[0]][t.c2 > 0] )
0294
0295         return [bi, c2sum, ndf, c2p, emn, emx, avp ]
0296
0297     LABELS = "bi c2sum ndf c2p emn emx avp".split()
0298
0299     def en_plot(self, c2overlay=0., c2poppy_=True):
0300         """
0301
0302         Using divide_edges is good for chi2 checking as it prevents
0303         bin migrations or "edge" effects.  But it leads to rather
0304         differently sized bins resulting in a strange histogram shape.
0305
0306         """
0307         ri = self.rindex
0308         s2c = self.s2c
0309         s2cn = self.s2cn
0310         bi = self.bi
0311         ibi = self.ibi
0312         c2 = self.c2
0313         c2poppy = self.c2poppy
0314         c2hscale = self.c2hscale
0315         hmax = self.hmax
0316         hl = self.hl
0317         hs = self.hs
0318         edges = self.edges
0319         en_slot = self.en_slot
0320         s2_slot = self.s2_slot
0321         cdf_slot = self.cdf_slot
0322
0323         emn = self.emn
0324         emx = self.emx
0325         avp = self.avp
0326
0327
0328         icdf_shape = str(s2cn.shape)
0329
0330         title_ = ["QCKTest.py : en_plot : lookup cf sampled : icdf_shape %s : %s " % ( icdf_shape, self.bi_base ),
0331                   "%s : %s " % ( self.c2label, self.cf),
0332                   "approach:%s use_icdf:%s avp %6.2f " % (self.approach, self.use_icdf, avp )
0333                  ]
0334
0335         title = "\n".join(title_)
0336         print(title)
0337         fig, ax = plt.subplots(figsize=[12.8, 7.2])
0338         fig.suptitle(title)
0339
0340         ax.plot( edges[:-1], hl[0], drawstyle="steps-post", label="lookup" )
0341         ax.plot( edges[:-1], hs[0], drawstyle="steps-post", label="sampled" )
0342
0343         if c2overlay != 0.:
0344             ax.plot( edges[:-1], c2*c2hscale*c2overlay , label="c2", drawstyle="steps-post" )
0345         pass
0346
0347         ylim = ax.get_ylim()
0348         xlim = ax.get_xlim()
0349
0350         s2max = s2cn[ibi,:,s2_slot].max()
0351         ax.plot( s2cn[ibi,:,en_slot], s2cn[ibi,:,s2_slot]*hmax/s2max , label="s2cn[%d,:,%d]*hmax/s2max (s2)" % (ibi,s2_slot) )
0352         ax.plot( s2cn[ibi,:,en_slot], s2cn[ibi,:,cdf_slot]*hmax ,       label="s2cn[%d,:,%d]*hmax      (cdf)" % (ibi, cdf_slot) )
0353
0354         ax.set_xlim(xlim)
0355         ax.set_ylim(ylim)
0356
0357         ax.plot( [emx, emx], ylim, linestyle="dotted" )
0358         ax.plot( [emn, emn], ylim, linestyle="dotted" )
0359
0360         if c2poppy_:
0361             for i in c2poppy:
0362                 ax.plot( [edges[i], edges[i]], ylim ,     label="c2poppy edge   %d " % i    , linestyle="dotted" )
0363                 ax.plot( [edges[i+1], edges[i+1]], ylim , label="c2poppy edge+1 %d " % (i+1), linestyle="dotted" )
0364             pass
0365         pass
0366
0367         ax.legend()
0368         fig.show()
0369         figpath = os.path.join(self.bi_base, "en_plot.png")
0370         log.info("savefig %s " % figpath)
0371         fig.savefig(figpath)
0372
0373     def compare(t, bis):
0374         res = np.zeros( (len(bis), len(t.LABELS)) )
0375         for i, bi in enumerate(bis):
0376             t.en_load(bi)
0377             res[i] = t.en_compare(bi)
0378             t.en_plot(c2overlay=0.5, c2poppy_=False)
0379             #t.rindex_plot()
0380         pass
0381         t.res = res
0382
0383     def __str__(self):
0384         title = "%s use_icdf:%s" % ( self.sample_base, self.use_icdf )
0385         underline = "=" * len(title)
0386         rst = RSTTable.Rdr(self.res, self.LABELS )
0387         return "\n".join(["", title, underline, "", rst])
0388
0389
0390 if __name__ == '__main__':
0391     logging.basicConfig(level=logging.INFO)
0392
0393     arg = sys.argv[1] if len(sys.argv) > 1 else None
0394
0395     #approach = "UpperCut"
0396     approach = "SplitBin"
0397     use_icdf = False
0398
0399     t = QCKTest(approach=approach, use_icdf=use_icdf)
0400
0401     if arg == "I":
0402         t.s2cn_plot(istep=20)
0403     else:
0404         bis = t.bislist()
0405         #bis = bis[-2:-1]
0406         #bis = [1.45,]
0407         #bis = [1.6,]
0408         t.compare(bis)
0409         print(t)
0410     pass
0411
0412
0413