qudarap/tests/QCtxTest.py

0001 #!/usr/bin/env python
0002 """
0003 QCtxTest.py
0004 =============
0005
0006 ::
0007
0008     qudarap
0009     TEST=E ipython -i tests/QCtxTest.py
0010
0011
0012
0013 """
0014 import os, sys, numpy as np
0015 from opticks.ana.nload import np_load
0016 from opticks.ana.key import keydir
0017 from matplotlib import pyplot as plt
0018
0019
0020
0021 class QCtxTest(object):
0022     FOLD = "/tmp/QCtxTest"
0023     hc_eVnm = 1240.
0024     colors = "rgbcmyk"
0025     figsize = [12.8,7.2]
0026     num = int(os.environ.get("NUM", "1000000"))
0027
0028     def globals(self, *args):
0029         assert len(args) % 2 == 0
0030         for i in range(int(len(args)//2)):
0031             k = args[2*i+0]
0032             v = args[2*i+1]
0033             print(" %10s : %s " % (k, str(v.shape)))
0034             globals()[k] = v
0035         pass
0036
0037
0038     def name(self, stem, suffix=None, num=None, flip_random=False):
0039         ss = ""
0040         ss += stem
0041         if not suffix is None:
0042            ss += suffix
0043         pass
0044         if flip_random:
0045            ss += "_FLIP_RANDOM"
0046         pass
0047         if not num is None:
0048            ss += "_%d" % num
0049         pass
0050         ss += ".npy"
0051         return ss
0052
0053     def path(self, stem, suffix=None, num=None, flip_random=False):
0054         name = self.name(stem, suffix=suffix, num=num, flip_random=flip_random)
0055         return os.path.join(self.FOLD, name)
0056
0057     def load(self, stem, suffix=None, num=None, flip_random=False):
0058         path = self.path(stem, suffix=suffix, num=num, flip_random=flip_random)
0059         a = np.load(path)
0060         print("QCtxTest.load %s : %s " % (str(a.shape), path))
0061         return a
0062
0063     def np_load(self, reldir):
0064         dirpath = os.path.join(self.FOLD, reldir)
0065         a, paths = np_load(dirpath)
0066         print("QCtxTest.np_load dirpath %s loaded %d paths " % (dirpath, len(paths)))
0067         print("\n".join(paths))
0068         return a
0069
0070
0071 class PropLookup(QCtxTest):
0072     def __init__(self):
0073         pp = self.load("prop_lookup_pp")
0074         x = self.load("prop_lookup_x")
0075         yy = self.load("prop_lookup_yy")
0076
0077         colors = self.colors
0078         self.pp = pp
0079         self.x = x
0080         self.yy = yy
0081
0082         fig, ax = plt.subplots(figsize=self.figsize)
0083         for i,y in enumerate(yy):
0084             ax.plot( x, y, color=colors[i] )
0085         pass
0086         for i,p in enumerate(pp):
0087             ni = p.view(np.uint32)[-1,-1]
0088             ax.scatter( p[:ni,0], p[:ni,1], color=colors[i] )
0089         pass
0090         fig.show()
0091
0092         globals()["pp"] = pp
0093         globals()["x"] = x
0094         globals()["yy"] = yy
0095
0096
0097
0098
0099 class CerenkovPhotonBase(QCtxTest):
0100     def __init__(self, suffix, num):
0101         p = self.load("cerenkov_photon", suffix=suffix, num=num)
0102
0103         en = p[:,0,0]
0104         wl = p[:,0,1]
0105         ri = p[:,0,2]
0106         ct = p[:,0,3]
0107
0108         s2 = p[:,1,0]
0109         bi = p[:,1,3]
0110
0111         w0 = p[:,2,0]
0112         w1 = p[:,2,1]
0113         u0 = p[:,2,2]
0114         u1 = p[:,2,3]
0115
0116         li = p[:,3,0].view(np.int32)
0117         lo = p[:,3,1].view(np.int32)
0118
0119         self.globals("p",p,"en",en,"wl",wl)
0120
0121
0122 class CerenkovPhoton(CerenkovPhotonBase):
0123     def __init__(self):
0124         CerenkovPhotonBase.__init__(self, suffix="", num=1000000)
0125
0126 class CerenkovPhotonEnprop(CerenkovPhotonBase):
0127     def __init__(self):
0128         CerenkovPhotonBase.__init__(self, suffix="_enprop", num=1000000)
0129
0130 class CerenkovPhotonExpt(CerenkovPhotonBase):
0131     def __init__(self):
0132         CerenkovPhotonBase.__init__(self, suffix="_expt", num=1000000)
0133
0134
0135
0136
0137 class RngSequence(QCtxTest):
0138     """
0139     Note slow histogramming + plotting as 256M randoms, shape (1M,16,16)
0140     """
0141     def __init__(self, reldir):
0142         r = self.np_load(reldir)
0143
0144         fig, axs = plt.subplots()
0145         fig.suptitle(reldir)
0146
0147         r_dom = np.linspace(0,1,256)
0148         h_r = np.histogram(r, r_dom )
0149
0150         ax = axs
0151         ax.plot( h_r[1][:-1], h_r[0], label="h_r", drawstyle="steps" )
0152
0153         ax.set_ylim( 0, h_r[0].max()*2. )
0154         ax.legend()
0155         fig.show()
0156
0157         self.globals("r",r)
0158
0159
0160
0161 class OldQCtxTest(object):
0162     FOLD = "/tmp/QCtxTest"
0163     hc_eVnm = 1240.
0164
0165     @classmethod
0166     def LoadCK(cls, num=10000):
0167         path = os.path.join( cls.FOLD, "cerenkov_photon_%d.npy" % num )
0168         p = np.load(path)
0169         return p
0170
0171     def scint_wavelength(self):
0172         """
0173         See::
0174
0175              ana/wavelength.py
0176              ana/wavelength_cfplot.py
0177
0178         """
0179         w0 = np.load(os.path.join(self.FOLD, "wavelength_scint_hd20.npy"))
0180
0181         path1 = "/tmp/G4OpticksAnaMgr/wavelength.npy"
0182         w1 = np.load(path1) if os.path.exists(path1) else None
0183
0184         kd = keydir(os.environ["OPTICKS_KEY"])
0185         aa = np.load(os.path.join(kd,"GScintillatorLib/GScintillatorLib.npy"))
0186         a = aa[0,:,0]
0187         b = np.linspace(0,1,len(a))
0188         u = np.random.rand(1000000)
0189         w2 = np.interp(u, b, a )
0190
0191         #bins = np.arange(80, 800, 4)
0192         bins = np.arange(300, 600, 4)
0193
0194         h0 = np.histogram( w0 , bins )
0195         h1 = np.histogram( w1 , bins )
0196         h2 = np.histogram( w2 , bins )
0197
0198         fig, ax = plt.subplots()
0199
0200         ax.plot( bins[:-1], h0[0], drawstyle="steps-post", label="OK.QCtxTest" )
0201         ax.plot( bins[:-1], h1[0], drawstyle="steps-post", label="G4" )
0202         ax.plot( bins[:-1], h2[0], drawstyle="steps-post", label="OK.GScint.interp" )
0203
0204         ylim = ax.get_ylim()
0205
0206         for w in [320,340,360,380,400,420,440,460,480,500,520,540]:
0207             ax.plot( [w,w], ylim )
0208         pass
0209
0210         ax.legend()
0211
0212         plt.show()
0213
0214         self.w0 = w0
0215         self.w1 = w1
0216         self.w2 = w2
0217
0218     def boundary_lookup_all(self):
0219         l = np.load(os.path.join(self.FOLD, "boundary_lookup_all.npy"))
0220         s_ = np.load(os.path.join(self.FOLD, "boundary_lookup_all_src.npy"))
0221         s = s_.reshape(l.shape)
0222         assert np.allclose(s, l)
0223
0224         self.s_ = s_
0225         self.s = s
0226         self.l = l
0227
0228     def boundary_lookup_line(self):
0229         p = np.load(os.path.join(self.FOLD, "boundary_lookup_line_props.npy"))
0230         w = np.load(os.path.join(self.FOLD, "boundary_lookup_line_wavelength.npy"))
0231
0232         path = "/tmp/RINDEXTest/g4_line_lookup.npy"
0233         g = np.load(path) if os.path.exists(path) else None
0234
0235         fig, ax = plt.subplots()
0236         ax.plot( w, p[:,0], drawstyle="steps", label="ri.qctx" )
0237
0238         if not g is None:
0239             assert np.all( w  == g[:,0] )
0240             ax.plot( w, g[:,1], drawstyle="steps", label="ri.g4" )
0241             ri_diff = p[:,0] - g[:,1]
0242             self.ri_diff = ri_diff
0243             print("ri_diff.min %s ri_diff.max %s " % (ri_diff.min(), ri_diff.max()))
0244         pass
0245
0246         #ax.plot( w, p[:,1], drawstyle="steps", label="abslen" )
0247         #ax.plot( w, p[:,2], drawstyle="steps", label="scatlen" )
0248         #ax.plot( w, p[:,3], drawstyle="steps", label="reemprob" )
0249
0250         ax.legend()
0251         fig.show()
0252
0253         self.p = p
0254         self.w = w
0255         self.g = g
0256
0257
0258
0259     def cerenkov_wavelength(self):
0260
0261         nm_dom = [80, 800, 1]
0262
0263         name = "wavelength_cerenkov"
0264         w0 = np.load(os.path.join(self.FOLD, "%s.npy" % name))
0265         e0 = self.hc_eVnm/w0
0266         wdom = np.arange(*nm_dom)
0267
0268         edom0 = self.hc_eVnm/wdom[::-1]
0269         # Convert to energy and reverse order for ascending energy domain
0270         # but thats a funny variable energy bin domain (with fixed 1 nm wavelength bin)
0271         #
0272         # Instead use equal energy bin domain across same range as wavelenth one
0273         # with same number of bins.
0274         edom = np.linspace( self.hc_eVnm/nm_dom[1], self.hc_eVnm/nm_dom[0], len(wdom) )
0275
0276
0277         h_w0 = np.histogram(w0, wdom )
0278         h_e0 = np.histogram(e0, edom )
0279
0280         h_w100 = np.histogram(w0, 100 )
0281         h_e100 = np.histogram(e0, edom )
0282
0283
0284
0285         fig, axs = plt.subplots(2,2)
0286         fig.suptitle(name)
0287
0288         ax = axs[0,0]
0289         ax.plot( h_w0[1][:-1], h_w0[0], label="h_w0", drawstyle="steps" )
0290         ax.legend()
0291
0292         ax = axs[0,1]
0293         ax.plot( h_e0[1][:-1], h_e0[0], label="h_e0", drawstyle="steps" )
0294         ax.legend()
0295
0296         ax = axs[1,0]
0297         ax.plot( h_w100[1][:-1], h_w100[0], label="h_w100", drawstyle="steps" )
0298         ax.legend()
0299
0300         ax = axs[1,1]
0301         ax.plot( h_e100[1][:-1], h_e100[0], label="h_e100", drawstyle="steps" )
0302         ax.legend()
0303
0304
0305
0306         fig.show()
0307
0308         self.w0 = w0
0309         self.e0 = e0
0310         self.h_w0 = h_w0
0311         self.h_e0 = h_e0
0312
0313         self.edom = edom
0314
0315
0316 if __name__ == '__main__':
0317
0318     test = os.environ.get("TEST", "X")
0319
0320     print("test [%s] " % (test))
0321     if test == 'Y':
0322         t = PropLookup()
0323     elif test == 'K':
0324         t = CerenkovPhoton()
0325     elif test == 'E':
0326         t = CerenkovPhotonEnprop()
0327     elif test == 'X':
0328         t = CerenkovPhotonExpt()
0329     elif test == 'F':
0330         reldir = "rng_sequence_f_ni1000000_nj16_nk16_tranche100000"
0331         t = RngSequence(reldir)
0332     else:
0333         print("test [%s] is not implemented" % test )
0334         pass
0335     pass
0336
0337