eic-opticks/ana/wavelength_plt.py

0001 #!/usr/bin/env python
0002 """
0003 wavelength.py
0004 ===============
0005
0006 w0
0007     qudarap/tests/QCtxTest wavelength samples from GPU texture
0008
0009 w1
0010     /tmp/G4OpticksAnaMgr/wavelength.npy from the horses mouth DsG4Scintillation
0011     collected by G4OpticksAnaMgr at BeginOfRun
0012
0013 w2
0014     np.interp results from the icdf obtained from GScintillatorLib, this
0015     closely matches w0 showing that the GPU texture lookups is working fine.
0016
0017 But w0 and w2 clearly exhibit some artifact of 20nm binning
0018 in the preparation of the icdf.
0019
0020 TODO: look into GScintillatorLib is should be using all the information
0021 provided from material properties, not using coarse bins.
0022
0023 """
0024 import os, numpy as np, logging
0025 log = logging.getLogger(__name__)
0026
0027 from opticks.ana.main import opticks_main
0028 from opticks.ana.key import keydir
0029 from opticks.ana.material import Material
0030 from opticks.ana.wavelength import Wavelength
0031
0032 from matplotlib import pyplot as plt
0033
0034 if __name__ == '__main__':
0035     ok = opticks_main()
0036     kd = keydir(os.environ["OPTICKS_KEY"])
0037     wl = Wavelength(kd)
0038
0039     h = wl.h
0040     l = wl.l
0041     dom = wl.dom[:-1]
0042
0043
0044     plt.ion()
0045     fig, ax = plt.subplots(figsize=ok.figsize)
0046
0047     fig.suptitle("LS : Absorption Length/Scattering Length/Reemission Probability/Spectrum Distribution")
0048
0049     ax.set_xlabel("wavelength (nm) [2 nm bins]" )
0050
0051     ax.set_ylabel("1M wavelength sample, 2nm bin counts")
0052     ax.plot( dom, h[0], drawstyle="steps", label=l[0] )
0053     ax.plot( dom, h[1], drawstyle="steps", label=l[1] )
0054     #ax.plot( dom, h[2], drawstyle="steps", label=l[2] )
0055     ax.legend(loc="center left")
0056
0057
0058     # https://stackoverflow.com/questions/9103166/multiple-axis-in-matplotlib-with-different-scales
0059     axr = ax.twinx()
0060     axr2 = ax.twinx()
0061
0062     ls = Material("LS")
0063     reempr = ls.reemission_prob(dom)
0064     abslen = ls.absorption_length(dom)
0065     scatlen = ls.scattering_length(dom)
0066
0067     axr.set_ylabel("length (m)")
0068     p1, = axr.plot( dom, abslen/1000., drawstyle="steps", label="abslen (m)", color="r" )
0069     p2, = axr.plot( dom, scatlen/1000., drawstyle="steps", label="scatlen (m)", color="g" )
0070
0071     axr2.set_ylabel("Re-emission Probability")
0072     axr2.spines['right'].set_position(('outward', 60))
0073     p3, = axr2.plot( dom, reempr, drawstyle="steps", label="reemprob", color="b" )
0074
0075     axr.legend(handles=[p1,p2,p3], loc="best")
0076
0077     lines = False
0078     if lines:
0079         ylim = np.array(ax.get_ylim())
0080         for w in [410,430]:
0081             ax.plot( [w,w], ylim*1.1 )
0082         pass
0083     pass
0084
0085     ax.axvspan(410, 430, color='grey', alpha=0.3, lw=0)
0086
0087     ax.text( 420, 100, "410-430 nm", ha="center" )
0088
0089     fig.tight_layout()
0090     plt.show()
0091
0092     path="/tmp/ana/LS_wavelength.png"
0093     fold=os.path.dirname(path)
0094     if not os.path.isdir(fold):
0095        os.makedirs(fold)
0096     pass
0097     log.info("save to %s " % path)
0098     fig.savefig(path)
0099
0100
0101