eic-opticks/ana/wavelength_cfplot.py

0001 #!/usr/bin/env python
0002 """
0003 wavelength_cfplot.py
0004 =======================
0005
0006 ::
0007
0008     an ;
0009
0010     ARG=0 ipython -i wavelength_cfplot.py
0011     ARG=1 ipython -i wavelength_cfplot.py
0012     ARG=2 ipython -i wavelength_cfplot.py
0013     ARG=3 ipython -i wavelength_cfplot.py
0014
0015     ARG=4 ipython -i wavelength_cfplot.py
0016     ARG=5 ipython -i wavelength_cfplot.py
0017     ARG=6 ipython -i wavelength_cfplot.py
0018     ARG=7 ipython -i wavelength_cfplot.py
0019     ARG=8 ipython -i wavelength_cfplot.py
0020     ARG=9 ipython -i wavelength_cfplot.py
0021     ARG=10 ipython -i wavelength_cfplot.py
0022     ARG=11 ipython -i wavelength_cfplot.py
0023     ARG=15 ipython -i wavelength_cfplot.py
0024
0025
0026     mkdir -p ~/simoncblyth.bitbucket.io/env/presentation/ana/wavelength_cfplot
0027     cp /tmp/ana/wavelength_cfplot/*.png ~/simoncblyth.bitbucket.io/env/presentation/ana/wavelength_cfplot/
0028
0029
0030 """
0031
0032 import os, numpy as np, logging
0033 log = logging.getLogger(__name__)
0034
0035 from opticks.ana.main import opticks_main
0036 from opticks.ana.key import keydir
0037 from opticks.ana.material import Material
0038 from opticks.ana.wavelength import Wavelength
0039 from opticks.ana.cfplot import one_cfplot
0040 from opticks.ana.cfh import CFH
0041
0042 from matplotlib import pyplot as plt
0043
0044 if __name__ == '__main__':
0045     ok = opticks_main()
0046     kd = keydir(os.environ["OPTICKS_KEY"])
0047     wl = Wavelength(kd)
0048     arg = int(os.environ.get("ARG","0"))
0049     a,b = wl.cf(arg)
0050
0051     wa = wl.w[a]
0052     wb = wl.w[b]
0053
0054     la = wl.l[a]
0055     lb = wl.l[b]
0056
0057     energy = False
0058     hc_eVnm = 1240.
0059
0060     dom = wl.dom[:-1]
0061
0062     if energy:
0063         wa = hc_eVnm/wa
0064         wb = hc_eVnm/wb
0065         dom = hc_eVnm/dom[::-1]
0066     pass
0067
0068     h = CFH()
0069
0070     if arg < 4:
0071         h.log = True
0072         h.ylim = (50., 5e4 )
0073         title = "Compare two 1M sample LS scintillation wavelength distributions in 1nm bins"
0074         xline = [wl.interp(0.05), wl.interp(0.95)]
0075     else:
0076         h.log = True
0077         title = "Compare two 2.82M samples of Cerenkov wavelength distributions in 1nm bins : poor chi2, interpolation artifact ?  "
0078         xline = []
0079     pass
0080
0081     h.suptitle = title
0082     c2cut = 10
0083     h(dom, wa, wb, [la, lb], c2cut=c2cut )
0084
0085     fig, axs = one_cfplot(ok, h, xline=xline )
0086
0087     if 1:
0088
0089         rindex = np.load(os.path.join(kd, "GScintillatorLib/LS_ori/RINDEX.npy"))
0090         rindex[:,0] *= 1e6
0091         if not energy:
0092             rindex[:,0] = hc_eVnm/rindex[:,0]
0093             rindex = rindex[::-1]
0094         pass
0095
0096         ax = axs[0]
0097
0098         ylim = ax.get_ylim()
0099         for i in range(len(rindex)):
0100             ax.plot( [rindex[i,0], rindex[i,0]], ylim , linestyle="dotted", color="b" )
0101         pass
0102
0103         axr = ax.twinx()
0104         axr.set_ylabel("rindex")
0105         axr.spines['right'].set_position(('outward', 0))
0106
0107         p3, = axr.plot( rindex[:,0] , rindex[:,1], drawstyle="steps-post", label="rindex", color="b" )
0108
0109         #wmin, wmax = 80., 800.
0110         wmin, wmax = 100., 400.
0111         emin, emax = hc_eVnm/wmax, hc_eVnm/wmin
0112
0113         if energy:
0114             xlim = [emin, emax]
0115         else:
0116             xlim = [wmin, wmax]
0117         pass
0118
0119         axs[0].set_xlim(*xlim)
0120         axs[1].set_xlim(*xlim)
0121
0122         axs[1].set_ylim(0,h.chi2.max()*1.1)
0123     pass
0124
0125     plt.ion()
0126     plt.show()
0127
0128     print(h.lhabc)
0129
0130     path="/tmp/ana/wavelength_cfplot/%s_%s.png" % ( wl.l[a], wl.l[b] )
0131     fold=os.path.dirname(path)
0132     if not os.path.isdir(fold):
0133        os.makedirs(fold)
0134     pass
0135     log.info("save to %s " % path)
0136     fig.savefig(path)
0137
0138