eic-opticks/ana/ScintillationIntegral.py

0001 #!/usr/bin/env python
0002 """
0003
0004 See also::
0005
0006    ana/wavelength.py
0007    ana/wavelength_plt.py
0008    ana/wavelength_cfplot.py
0009    extg4/tests/X4ScintillationIntegralTest.py
0010
0011
0012
0013 """
0014 import os, numpy as np
0015 import matplotlib.pyplot as plt
0016 from opticks.ana.key import keydir
0017
0018 def wavelengthSample(icdf, num_sample=1000000):
0019     u = np.random.rand(num_sample)
0020     w = np.interp( u, icdf.si, icdf.nm )
0021     return w
0022
0023
0024 class Scint(object):
0025     def __repr__(self):
0026         return "%s : nm %s si %s " % (self.LABEL, str(self.nm.shape), str(self.si.shape))
0027
0028 class ScintillationIntegral(Scint):
0029     """
0030     aa,fc,sc
0031         obtained from horses mouth DsG4Scintillation by G4OpticksAnaMgr
0032         at BeginOfRun and grabbed with ana/G4OpticksAnaMgr.sh
0033     bb
0034         integrated from fc by extg4/tests/X4ScintillationIntegralTest.cc
0035
0036
0037     1239.84198/(si.fc[:,0]*1e6)
0038
0039     """
0040     LABEL = "G4.ScintillationIntegral"
0041     FOLD = "/tmp/G4OpticksAnaMgr"
0042
0043     def __init__(self):
0044         aa = np.load(os.path.join(self.FOLD,"ScintillationIntegral.npy"))
0045         fc = np.load(os.path.join(self.FOLD,"FASTCOMPONENT.npy"))
0046         sc = np.load(os.path.join(self.FOLD,"SLOWCOMPONENT.npy"))
0047         assert np.all(fc == sc)
0048
0049         bb = np.load(os.path.join(self.FOLD,"X4ScintillationIntegralTest.npy"))
0050         assert np.all(aa == bb)
0051
0052         eV = aa[:,0]*1e6
0053         si = aa[:,1]/aa[:,1].max()
0054
0055
0056         hc_eVnm = 1239.84198
0057         #hc_eVnm = 1240.
0058
0059         nm = hc_eVnm/eV    ## 330,331,...,599,600  very close to integer nm with 1239.84198
0060
0061         self.aa = aa
0062         self.bb = bb
0063         self.fc = fc
0064         self.sc = sc
0065         self.eV = eV
0066         self.nm = nm
0067         self.si = si
0068         self.wl = wavelengthSample(self)
0069
0070
0071 class GScintillator(Scint):
0072     LABEL = "OK.GScint"
0073     def __init__(self, kd):
0074         aa = np.load(os.path.join(kd,"GScintillatorLib/GScintillatorLib.npy"))
0075         fc = np.load(os.path.join(kd,"GScintillatorLib/LS/FASTCOMPONENT.npy"))
0076         sc = np.load(os.path.join(kd,"GScintillatorLib/LS/SLOWCOMPONENT.npy"))
0077
0078         a = aa[0,:,0]
0079         b = np.linspace(0,1,len(a))
0080
0081         self.nm = a
0082         self.si = b
0083         self.wl = wavelengthSample(self)
0084
0085
0086 def compare_icdf(gs, si):
0087     fig = plt.figure()
0088     plt.title("Inverted Cumulative Distribution Function : for Scintillator Reemission " )
0089     ax = fig.add_subplot(1,1,1)
0090
0091     ax.plot( gs.si, gs.nm, label=gs.LABEL)
0092     ax.plot( si.si, si.nm, label=si.LABEL)
0093
0094     ax.set_ylabel("Wavelength (nm)")
0095     ax.set_xlabel("Probability")
0096     ax.legend()
0097
0098     #ax.set_xscale('log')
0099
0100     fig.show()
0101
0102 def compare_samples(gs, si):
0103
0104     bins = np.arange(300, 600, 2 )
0105     gsh = np.histogram( gs.wl, bins )
0106     sih = np.histogram( si.wl, bins )
0107
0108     fig, ax = plt.subplots()
0109     ax.plot( bins[:-1], gsh[0], drawstyle="steps-post", label=gs.LABEL )
0110     ax.plot( bins[:-1], sih[0], drawstyle="steps-post", label=si.LABEL )
0111     ax.legend()
0112
0113     fig.show()
0114
0115
0116
0117 def integral(sc):
0118     out = np.zeros((len(sc),2))
0119
0120     currentIN = sc[0,1]
0121     assert currentIN >= 0.0
0122
0123     currentPM = sc[0,0]
0124     currentCII = 0.0;
0125
0126     out[0] = [currentPM, currentCII]
0127
0128     prevPM  = currentPM
0129     prevCII = currentCII
0130     prevIN  = currentIN
0131
0132     for ii in range(1,len(sc)):
0133         currentPM = sc[ii,0]
0134         currentIN = sc[ii,1]
0135         currentCII = 0.5 * (prevIN + currentIN)
0136         currentCII = prevCII + (currentPM - prevPM) * currentCII
0137         out[ii] = [currentPM, currentCII]
0138         prevPM  = currentPM
0139         prevCII = currentCII
0140     pass
0141     return out
0142
0143
0144
0145
0146 if __name__ == '__main__':
0147
0148     ok = os.environ["OPTICKS_KEY"]
0149     kd = keydir(ok)
0150
0151     gs = GScintillator(kd)
0152     si = ScintillationIntegral()
0153
0154     print("gs:%r " % gs )
0155     print("si:%r " % si )
0156
0157
0158     compare_icdf(gs, si)
0159     #compare_samples(gs, si)
0160
0161     ff = integral(si.sc)
0162
0163