eic-opticks/ana/wavelength.py

0001 #!/usr/bin/env python
0002 """
0003 wavelength.py
0004 ===============
0005
0006 ::
0007
0008    ARG=6 ipython -i wavelength.py
0009    ARG=7 ipython -i wavelength.py
0010    ARG=8 ipython -i wavelength.py
0011    ARG=11 ipython -i wavelength.py
0012    ARG=12 ipython -i wavelength.py
0013    ARG=13 ipython -i wavelength.py
0014    ARG=15 ipython -i wavelength.py
0015
0016 """
0017
0018 import os, numpy as np, logging
0019 log = logging.getLogger(__name__)
0020 from opticks.ana.key import keydir
0021 from opticks.ana.nbase import chi2
0022
0023 class Wavelength(object):
0024     """
0025     Comparing wavelength distribs between many different samples
0026     """
0027     FOLD = "/tmp/wavelength"
0028     def get_key(self, label):
0029         key = None
0030         for k,v in self.l.items():
0031             if v == label:
0032                 key = k
0033             pass
0034         pass
0035         return key
0036
0037     def get_keys(self, a_label, b_label):
0038         a = self.get_key(a_label)
0039         b = self.get_key(b_label)
0040         return a, b
0041
0042     def __call__(self, label):
0043         return self.get_key(label)
0044
0045
0046     def format(self, i):
0047         return " %2d : %d : %50s : %s " % (i, os.path.exists(self.p[i]), self.l[i], self.p[i])
0048
0049     def __init__(self, kd):
0050
0051
0052         p = {}
0053         l = {}
0054
0055         l[0] = "DsG4Scintillator_G4OpticksAnaMgr"    ## horses mouth
0056         p[0] = "/tmp/G4OpticksAnaMgr/WavelengthSamples.npy"
0057
0058         l[1] = "Opticks_QCtxTest_hd20"
0059         p[1] = os.path.join("/tmp/QCtxTest", "wavelength_20.npy")
0060
0061         l[2] = "Opticks_QCtxTest_hd0"
0062         p[2] = os.path.join("/tmp/QCtxTest", "wavelength_0.npy")
0063
0064         l[3] = "Opticks_QCtxTest_hd20_cudaFilterModePoint"
0065         p[3] = os.path.join("/tmp/QCtxTest", "wavelength_20_cudaFilterModePoint.npy")
0066
0067         l[4] = "Opticks_QCtxTest_hd0_cudaFilterModePoint"
0068         p[4] = os.path.join("/tmp/QCtxTest", "wavelength_0_cudaFilterModePoint.npy")
0069
0070         l[5] = "X4"
0071         p[5] = "/tmp/X4ScintillationTest/g4localSamples.npy"
0072
0073         l[6] = "GScintillatorLib_np_interp"
0074         p[6] = os.path.join(kd,"GScintillatorLib/GScintillatorLib.npy")
0075
0076         l[7] = "ck_photon_1M"
0077         p[7] = os.path.join("/tmp/QCtxTest", "cerenkov_photon_1000000.npy")
0078
0079         l[8] = "G4Cerenkov_modified_SKIP_CONTINUE"
0080         p[8] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_step_length_100000.000_SKIP_CONTINUE", "GenWavelength.npy")
0081
0082         l[9] = "G4Cerenkov_modified_ASIS"
0083         p[9] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_step_length_100000.000_ASIS", "GenWavelength.npy")
0084
0085         l[10] = "G4Cerenkov_modified_SKIP_CONTINUE_10k"
0086         p[10] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_override_fNumPhotons_10000_SKIP_CONTINUE", "GenWavelength.npy")
0087
0088         l[11] = "G4Cerenkov_modified_SKIP_CONTINUE_1M"
0089         p[11] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_override_fNumPhotons_1000000_SKIP_CONTINUE", "GenWavelength.npy")
0090
0091         l[12] = "G4Cerenkov_modified_SKIP_CONTINUE_1M_seed1"
0092         p[12] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_override_fNumPhotons_1000000_SKIP_CONTINUEseed_1_", "GenWavelength.npy")
0093
0094         l[13] = "G4Cerenkov_modified_SKIP_CONTINUE_1M_seed2"
0095         p[13] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_override_fNumPhotons_1000000_SKIP_CONTINUEseed_2_", "GenWavelength.npy")
0096
0097         l[14] = "G4Cerenkov_modified_SKIP_CONTINUE_1M_FLOAT_TEST"
0098         p[14] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_override_fNumPhotons_1000000_SKIP_CONTINUE_FLOAT_TEST", "GenWavelength.npy")
0099
0100         l[15] = "ck_photon_1M_FLIP_RANDOM"
0101         p[15] = os.path.join("/tmp/QCtxTest", "cerenkov_photon_FLIP_RANDOM_1000000.npy")
0102
0103         l[16] = "G4Cerenkov_modified_SKIP_CONTINUE_1M_seed1f"
0104         p[16] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_override_fNumPhotons_1000000_SKIP_CONTINUE_FLOAT_TEST_seed_1_", "GenWavelength.npy")
0105
0106         l[17] = "G4Cerenkov_modified_SKIP_CONTINUE_1M_seed2f"
0107         p[17] = os.path.join("/tmp/G4Cerenkov_modifiedTest", "BetaInverse_1.500_override_fNumPhotons_1000000_SKIP_CONTINUE_FLOAT_TEST_seed_2_", "GenWavelength.npy")
0108
0109         l[18] = "ana_ck_1M"
0110         p[18] = "/tmp/ck/ck_1000000.npy"
0111
0112         l[19] = "G4Cerenkov_modified_SKIP_CONTINUE_1M_PRECOOKED"
0113         p[19] = "/tmp/G4Cerenkov_modifiedTest/BetaInverse_1.500_override_fNumPhotons_1000000_SKIP_CONTINUE_PRECOOKED/GenWavelength.npy"
0114
0115         l[20] = "ck_photon_enprop_1M"
0116         p[20] = os.path.join("/tmp/QCtxTest", "cerenkov_photon_enprop_1000000.npy")
0117
0118         l[21] = "ck_photon_expt_1M"
0119         p[21] = os.path.join("/tmp/QCtxTest", "cerenkov_photon_expt_1000000.npy")
0120
0121         l[22] = "rindex_en_integrated_1M"
0122         p[22] = "/tmp/rindex/en_integrated_lookup_1M.npy"
0123
0124         self.p = p
0125         self.l = l
0126
0127         dom = np.arange(80, 400, 4)
0128         #dom = np.arange(300, 600, 1)
0129         #dom = np.arange(385, 475, 1)
0130
0131
0132         #dom = np.arange(350, 550, 1)
0133
0134
0135         a = {}
0136         e = {}
0137         w = {}
0138         h = {}
0139         r = {}
0140
0141         for i in range(len(l)):
0142             p_exists = os.path.exists(p[i])
0143             print(self.format(i))
0144             if not p_exists:
0145                 a[i] = None
0146                 w[i] = None
0147                 h[i] = None
0148             else:
0149                 a[i] = np.load(p[i])
0150                 if l[i].startswith("ck_photon"):
0151                     e[i] = a[i][:,0,0]
0152                     w[i] = a[i][:,0,1]
0153                     r[i] = a[i][:,0,2]
0154                 elif l[i].startswith("rindex_en_integrated"):
0155                     e[i] = a[i]
0156                     w[i] = 1240./e[i]
0157                     r[i] = np.zeros( len(a[i]) )
0158                 elif l[i].startswith("ana_ck"):
0159                     w[i] = a[i][:,0,1]
0160                 elif l[i].startswith("G4Cerenkov_modified"):
0161                     e[i] = a[i][:,0,0]
0162                     w[i] = a[i][:,0,1]
0163                     r[i] = a[i][:,0,2]
0164                 elif l[i] == "GScintillatorLib_np_interp":
0165                     aa = a[i]
0166                     self.aa = aa
0167                     aa0 = aa[0,:,0]
0168                     bb0 = np.linspace(0,1,len(aa0))
0169                     u = np.random.rand(1000000)
0170                     w[i] = np.interp(u, bb0, aa0 )
0171                 else:
0172                     w[i] = a[i]
0173                 pass
0174                 h[i] = np.histogram( w[i] , dom )
0175             pass
0176         pass
0177         self.w = w
0178         self.e = e
0179         self.r = r
0180         self.h = h
0181         self.a = a
0182         self.dom = dom
0183
0184     def interp(self, u):
0185         a = self.aa[0,:,0]
0186         b = np.linspace(0,1,len(a))
0187         return np.interp( u, b, a )
0188
0189
0190     def cf(self, arg):
0191         if arg == 0:
0192             a, b = self.get_keys('DsG4Scintillator_G4OpticksAnaMgr', "Opticks_QCtxTest_hd20")
0193         elif arg == 1:
0194             a, b = self.get_keys('DsG4Scintillator_G4OpticksAnaMgr', "Opticks_QCtxTest_hd0")
0195         elif arg == 2:
0196             a, b = self.get_keys('DsG4Scintillator_G4OpticksAnaMgr', 'Opticks_QCtxTest_hd20_cudaFilterModePoint')
0197         elif arg == 3:
0198             a, b = self.get_keys('DsG4Scintillator_G4OpticksAnaMgr', 'Opticks_QCtxTest_hd0_cudaFilterModePoint')
0199         elif arg == 4:
0200             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE', 'ck_photon' )
0201         elif arg == 5:
0202             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_10k', 'ck_photon_10k' )
0203         elif arg == 6:
0204             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M', 'ck_photon_1M' )
0205         elif arg == 7:
0206             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M_seed1', 'G4Cerenkov_modified_SKIP_CONTINUE_1M_seed2' )
0207         elif arg == 8:
0208             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M_FLOAT_TEST', 'ck_photon_1M' )
0209         elif arg == 9:
0210             a, b = self.get_keys('ck_photon_1M', 'ck_photon_1M_FLIP_RANDOM' )
0211         elif arg == 10:
0212             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M_seed1f', 'G4Cerenkov_modified_SKIP_CONTINUE_1M_seed2f' )
0213         elif arg == 11:
0214             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M_PRECOOKED', 'ana_ck_1M' )
0215         elif arg == 12:
0216             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M_PRECOOKED', 'ck_photon_enprop_1M' )
0217         elif arg == 13:
0218             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M', 'ck_photon_enprop_1M' )
0219         elif arg == 14:
0220             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M_FLOAT_TEST', 'ck_photon_enprop_1M' )
0221         elif arg == 15:
0222             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M', 'ck_photon_expt_1M' )
0223         elif arg == 16:
0224             a, b = self.get_keys('G4Cerenkov_modified_SKIP_CONTINUE_1M', 'rindex_en_integrated_1M' )
0225         else:
0226             assert 0
0227         pass
0228
0229         print("arg:%d -> a:%d b:%d " % (arg,a,b))
0230         print(self.format(a))
0231         print(self.format(b))
0232         return a, b
0233
0234
0235
0236 if __name__ == '__main__':
0237     kd = keydir(os.environ["OPTICKS_KEY"])
0238     ri = np.load(os.path.join(kd, "GScintillatorLib/LS_ori/RINDEX.npy"))
0239     ri[:,0] *= 1e6
0240     ri_ = lambda e:np.interp(e, ri[:,0], ri[:,1] )
0241
0242
0243     wl = Wavelength(kd)
0244     arg = int(os.environ.get("ARG","0"))
0245     ia,ib = wl.cf(arg)
0246
0247     pa = wl.p[ia]
0248     pb = wl.p[ib]
0249
0250     a = wl.a[ia]
0251     b = wl.a[ib]
0252
0253     wa = wl.w[ia]
0254     wb = wl.w[ib]
0255
0256     la = wl.l[ia]
0257     lb = wl.l[ib]
0258
0259     ea = wl.e[ia]
0260     eb = wl.e[ib]
0261
0262     ra = wl.r[ia]
0263     rb = wl.r[ib]
0264
0265
0266     print("la:%s" % la)
0267     print("pa:%s" % pa)
0268     print("lb:%s" % lb)
0269     print("pb:%s" % pb)
0270
0271
0272     if 0:
0273         wab = np.abs(wa-wb)
0274
0275         dev = wab > 1e-4
0276         num_dev = np.count_nonzero(dev)
0277         print("num_dev:%d " % num_dev )
0278
0279         np.set_printoptions(edgeitems=16)
0280
0281         print("a[dev]\n",a[dev].reshape(-1,8))
0282         print("b[dev]\n",b[dev].reshape(-1,16))
0283         b_ri = b[dev][:,0,2]
0284
0285     # are deviants mostly where sampled rindex bin values
0286     # TODO: enable dumping of deviants, to understand the reason
0287
0288     if 0:
0289         fig, ax = plt.subplots()
0290         ax.hist( b_ri, bins=50 )
0291         fig.show()
0292
0293         #mask = np.where(dev)[0]
0294         #np.save("/tmp/wavelength_deviant_mask.npy", mask)
0295
0296     if 0:
0297         # 2d
0298         n = 1000000
0299         fig, ax = plt.subplots()
0300         ax.scatter( ea[:n], ra[:n], color="r", label="a", s=0.1 )
0301         ax.scatter( eb[:n], rb[:n], color="b", label="b", s=0.1  )
0302         ax.legend()
0303         fig.show()
0304     pass
0305
0306     if 1:
0307         # energy histogram in a lot of bins
0308
0309         dom = 3,10,100
0310         #dom = 3,10,1001
0311         edom = np.linspace(*dom)
0312         ea_h = np.histogram( ea, edom )
0313         eb_h = np.histogram( eb, edom )
0314
0315         prefix = "energy_chi2_dom%d_arg%d" % (dom[2], arg)
0316         stem = "%s_%s_%s" % (prefix,la,lb)
0317         figpath = os.path.join(wl.FOLD, "%s.png" % stem )
0318
0319         c2 = chi2(ea_h[0],eb_h[0],cut=10)
0320         c2ndf = c2[0].sum()/c2[1]
0321         c2_smry = " c2/ndf = %6.2f/%3d = %4.2f " % ( c2[0].sum(), c2[1], c2ndf )
0322         dom_smry = " edom %s %s nbin:%s " % (dom[0], dom[1], dom[2] - 1)
0323         title = "\n".join([figpath,c2_smry,dom_smry])
0324
0325         print(title)
0326
0327         figsize = [12.8, 7.2]
0328         fig, ax = plt.subplots(1, figsize=figsize)
0329         fig.suptitle( title )
0330
0331         ax.plot( edom[:-1], ea_h[0], label=la, drawstyle="steps-post" )
0332         ax.plot( edom[:-1], eb_h[0], label=lb, drawstyle="steps-post" )
0333         ax.plot( edom[:-1], -2000*c2[0]/c2[0].max(), drawstyle="steps-post", label="chi2")
0334
0335         eri = ri[:,0]
0336         eris = eri[np.logical_and(eri>edom[0], eri<edom[-1])]
0337
0338         ylim = ax.get_ylim()
0339         for e in eris:
0340             ax.plot( [e,e], ylim, linestyle="dotted", color="b" )
0341         pass
0342
0343         ax.legend()
0344
0345         #x0,x1 = 3,7
0346         #x0,x1 = 7,10
0347         #ax.set_xlim(x0,x1)
0348         pass
0349         fig.show()
0350
0351
0352         figfold = os.path.dirname(figpath)
0353         if not os.path.exists(figfold):
0354             os.makedirs(figfold)
0355         pass
0356         print("figpath : %s " % figpath)
0357         fig.savefig(figpath)
0358
0359     pass
0360
0361
0362