eic-opticks/ana/trainbow.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 trainbow.py : Rainbow deviation angle comparison between Opticks and Geant4
0023 ====================================================================================
0024
0025 To simulate the rainbow events::
0026
0027    trainbow-
0028
0029    trainbow-- --spol
0030    trainbow-- --ppol
0031
0032    trainbow-- --spol --tcfg4
0033    trainbow-- --ppol --tcfg4
0034
0035
0036
0037 Expected output, is a scattering angle plot and history comparison table:
0038
0039 .. code-block:: py
0040
0041     In [28]: run rainbow_cfg4.py
0042     WARNING:opticks.ana.evt:init_index S-Pol G4 finds too few (ps)phosel uniques : 1
0043     WARNING:opticks.ana.evt:init_index S-Pol G4 finds too few (rs)recsel uniques : 1
0044     WARNING:opticks.ana.evt:init_index S-Pol G4 finds too few (rsr)reshaped-recsel uniques : 1
0045                            5:rainbow   -5:rainbow           c2
0046                     8ccd        819160       819654             0.15  [4 ] TO BT BT SA
0047                      8bd        102089       101615             1.10  [3 ] TO BR SA
0048                    8cbcd         61869        61890             0.00  [5 ] TO BT BR BT SA
0049                   8cbbcd          9617         9577             0.08  [6 ] TO BT BR BR BT SA
0050                  8cbbbcd          2604         2687             1.30  [7 ] TO BT BR BR BR BT SA
0051                 8cbbbbcd          1056         1030             0.32  [8 ] TO BT BR BR BR BR BT SA
0052                    86ccd          1014         1000             0.10  [5 ] TO BT BT SC SA
0053                8cbbbbbcd           472          516             1.96  [9 ] TO BT BR BR BR BR BR BT SA
0054                      86d           498          473             0.64  [3 ] TO SC SA
0055               bbbbbbbbcd           304          294             0.17  [10] TO BT BR BR BR BR BR BR BR BR
0056               8cbbbbbbcd           272          247             1.20  [10] TO BT BR BR BR BR BR BR BT SA
0057               cbbbbbbbcd           183          161             1.41  [10] TO BT BR BR BR BR BR BR BR BT
0058                      4cd           161          139             1.61  [3 ] TO BT AB
0059                     86bd           138          142             0.06  [4 ] TO BR SC SA
0060                    8c6cd           126          106             1.72  [5 ] TO BT SC BT SA
0061                     4ccd           100          117             1.33  [4 ] TO BT BT AB
0062                   86cbcd            88          110             2.44  [6 ] TO BT BR BT SC SA
0063                       4d            51           54             0.09  [2 ] TO AB
0064                    8cc6d            38           40             0.05  [5 ] TO SC BT BT SA
0065                  8cc6ccd            19           33             3.77  [7 ] TO BT BT SC BT BT SA
0066                              1000000      1000000         1.09
0067
0068
0069
0070 * see also: :doc:`rainbow_cfg4_notes`
0071
0072
0073
0074 """
0075 import os, sys, logging, numpy as np
0076 log = logging.getLogger(__name__)
0077
0078 import matplotlib.pyplot as plt
0079 import matplotlib.gridspec as gridspec
0080 from mpl_toolkits.mplot3d import Axes3D
0081 from matplotlib.patches import Rectangle
0082
0083
0084 from opticks.ana.base import opticks_main
0085 from opticks.ana.evt import Evt, costheta_, cross_
0086 from opticks.ana.boundary import Boundary
0087 from opticks.ana.droplet import Droplet
0088 from opticks.ana.nbase import chi2
0089
0090 X,Y,Z,W = 0,1,2,3
0091
0092 deg = np.pi/180.
0093 n2ref = 1.33257
0094
0095
0096 def a_scatter_plot_cf(ax, a_evt, b_evt, log_=False):
0097     """
0098     :param ax: mpl axis
0099     :param a_evt:
0100     :param b_evt:
0101     :param log_: log scale
0102
0103     :return: cnt, bns dicts keyed with 0,1
0104
0105     Scattering angle in degrees 0 to 360
0106
0107     Histograms result of a_deviation_angle for each evt,
0108     storing counts, bins and patches in dicts
0109     with keys 0, 1
0110
0111     """
0112     db = np.arange(0,360,1)
0113
0114     incident = np.array([0,0,-1])
0115     cnt = {}
0116     bns = {}
0117     ptc = {}
0118     j = -1
0119     for i,evt in enumerate([a_evt, b_evt]):
0120         dv = evt.a_deviation_angle(axis=X, incident=incident)/deg
0121         ax.set_xlim(0,360)
0122         if len(dv) > 0:
0123             cnt[i], bns[i], ptc[i] = ax.hist(dv, bins=db,  log=log_, histtype='step', label=evt.label)
0124             j = i
0125     pass
0126     if len(bns) == 2:
0127         assert np.all( bns[0] == bns[1] )
0128
0129     if j == -1:
0130         bns = None
0131     else:
0132         bns = bns[j]
0133
0134     return cnt, bns
0135
0136
0137 def cf_plot(evt_a, evt_b, label="", log_=False, ylim=[1,1e5], ylim2=[0,10], sli=slice(0,10)):
0138
0139     tim_a = " ".join(map(lambda f:"%5.2f" % f, map(float, filter(None, evt_a.tdii['propagate']) )[sli]))
0140     tim_b = " ".join(map(lambda f:"%5.2f" % f, map(float, filter(None, evt_b.tdii['propagate']) )[sli]))
0141
0142     fig = plt.figure()
0143     suptitle = "Rainbow cfg4 " + label + "[" + tim_a + "] [" + tim_b + "]"
0144     log.info("plotting %s " % suptitle)
0145     fig.suptitle(suptitle  )
0146
0147     gs = gridspec.GridSpec(2, 1, height_ratios=[3,1])
0148
0149     ax = fig.add_subplot(gs[0])
0150
0151     c, bns = a_scatter_plot_cf(ax, evt_a, evt_b, log_=log_)
0152     droplet.bow_angle_rectangles()
0153     ax.set_ylim(ylim)
0154     ax.legend()
0155
0156     xlim = ax.get_xlim()
0157
0158     ax = fig.add_subplot(gs[1])
0159
0160     if len(c) == 2:
0161         a,b = c[0],c[1]
0162
0163         c2, c2n, c2nn = chi2(a, b, cut=30)
0164         c2p = c2.sum()/c2n
0165
0166         plt.plot( bns[:-1], c2, drawstyle='steps', label="chi2/ndf %4.2f" % c2p )
0167         ax.set_xlim(xlim)
0168         ax.legend()
0169
0170         ax.set_ylim(ylim2)
0171
0172         droplet.bow_angle_rectangles()
0173
0174
0175
0176 if __name__ == '__main__':
0177     args = opticks_main(tag="5",src="torch",det="rainbow",doc=__doc__)
0178
0179
0180     boundary = Boundary("Vacuum///MainH2OHale")
0181     droplet = Droplet(boundary)
0182
0183     plt.ion()
0184     plt.close()
0185
0186     tag = args.tag
0187     src = args.src
0188     det = args.det
0189     rec = True
0190     log_ = True
0191     not_ = False
0192
0193     if det == "rainbow":
0194        if tag == "5":
0195            label = "S-Pol"
0196        elif tag == "6":
0197            label = "P-Pol"
0198        else:
0199            label = "no label"
0200
0201
0202     seqs = []
0203
0204     try:
0205         a = Evt(tag=tag, src=src, det=det, label="%s Op" % label, seqs=seqs, not_=not_, rec=rec, args=args)
0206         b = Evt(tag="-%s" % tag, src=src, det=det, label="%s G4" % label, seqs=seqs, not_=not_, rec=rec, args=args)
0207     except IOError as err:
0208         log.fatal(err)
0209         sys.exit(args.mrc)
0210
0211     print a.brief
0212     print b.brief
0213
0214     if not (a.valid and b.valid):
0215         log.fatal("need two valid events to compare ")
0216         sys.exit(1)
0217
0218
0219     lmx = 20
0220     #hcf = a.history.table.compare(b.history.table)
0221     hcf = a.his.compare(b.his)
0222     if len(hcf.lines) > lmx:
0223         hcf.sli = slice(0,lmx)
0224     print hcf
0225
0226     cf_plot(a, b, label=label, log_=log_, ylim=[0.8,4e4],ylim2=None)
0227
0228
0229 ## EXERCISE : ENABLE THE BELOW PLOTS AND INTERPRET WHAT YOU GET
0230
0231 if 0:
0232     #seqs = a.history.table.labels[:5]
0233     seqs = a.his.labels[:5]
0234     for seq in seqs:
0235         qa =  Evt(tag=tag, src=src, det=det, label="%s Op" % label, seqs=[seq], not_=not_, rec=rec, args=args)
0236         qb =  Evt(tag="-%s" % tag, src=src, det=det, label="%s G4" % label, seqs=[seq], not_=not_, rec=rec, args=args)
0237
0238         cf_plot(qa, qb, label=label + " " + seq, log_=log_, ylim=[0.8,4e4],ylim2=None)
0239
0240
0241