eic-opticks/ana/treflect.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 treflect.py
0023 =======================
0024
0025 See :doc:`../tests/treflect`
0026
0027 Loads S and P polarized events, selects photons that
0028 have histories : "TO BR SA" or "TO BR AB",
0029 ie photons that reflect at the first interface and
0030 are either bulk absorbed (AB) or surface absorbed (SA).
0031
0032 A ratio of histograms of the angle of incidence (theta)
0033 for all photons with reflecting photons is compared
0034 against expectations of the Fresnel formula.
0035
0036
0037
0038 """
0039 import os, sys, logging
0040 import numpy as np
0041 log = logging.getLogger(__name__)
0042
0043 import matplotlib.pyplot as plt
0044
0045 from opticks.ana.base import opticks_main
0046 from opticks.ana.nbase import count_unique
0047 from opticks.ana.evt import Evt, costheta_
0048 from opticks.ana.ana import Rat, theta, recpos_plot, angle_plot
0049 from opticks.ana.boundary import Boundary
0050 from opticks.ana.fresnel import Fresnel
0051
0052 deg = np.pi/180.
0053
0054 np.set_printoptions(suppress=True, precision=3)
0055
0056 rat_ = lambda n,d:float(len(n))/float(len(d))
0057
0058 X,Y,Z,W = 0,1,2,3
0059
0060
0061 class Reflect(object):
0062
0063     def theta(self, vecs, normal):
0064         N = len(vecs)
0065         nrm = np.tile(normal, N).reshape(-1, len(normal))
0066         ct = costheta_(vecs, nrm)
0067         th = np.arccos(ct)*180./np.pi
0068         return th
0069
0070     def __init__(self, evt, focus, normal):
0071         """
0072         :param evt:
0073         :param focus: coordinates of reflection point
0074
0075
0076         p0a initial position with no selection applied
0077
0078         """
0079
0080         al = evt
0081         p0a = al.rpost_(0)[:,:W] - focus
0082         tha = self.theta( p0a, normal )
0083
0084         br = Evt.selection(evt, seqs=["TO BR SA","TO BR AB"], label="TO BR ..")
0085
0086         p0 = br.rpost_(0)[:,:W] - focus   # start
0087         p1 = br.rpost_(1)[:,:W] - focus   # refl point
0088         p2 = br.rpost_(2)[:,:W] - focus   # end
0089
0090         #miss = p1[:,2] != p1[0,2]   # this assumes a particular geometry
0091         # missers confirmed to be a triangulation crack effect
0092         # avoid by not targetting cracks with millions of photons whilst doing reflection tests
0093         # ...get **ZERO** missers when avoid cracks ...
0094
0095
0096         if 0:
0097             miss = np.tile(False, len(p1))
0098             self.miss = miss
0099             msk = ~miss
0100             p0u = p0[msk]
0101             p2u = p2[msk]
0102         else:
0103             p0u = p0
0104             p2u = p2
0105
0106
0107         th0 = self.theta(p0u, normal)
0108         th2 = self.theta(p2u, normal)
0109         th = th0
0110
0111         missrat = Rat(th0,p0,"th0/p0")
0112
0113         rats = []
0114         rats.append(Rat(p0,tha,  "p0/tha"))
0115         rats.append(Rat(th0,tha, "th0/tha"))
0116         rats.append(missrat)
0117
0118         msg = " ".join(map(repr,rats))
0119         log.info("%s : %s " % (repr(evt),msg))
0120
0121
0122         # th0 and th2 are incident and reflected angles
0123         # matching very well after exclude miss-ers
0124
0125         self.evt = evt
0126         self.focus = focus
0127         self.normal = normal
0128         self.al = al
0129         self.br = br
0130         self.p0a = p0a
0131         self.tha = tha
0132         self.p0 = p0
0133         self.p1 = p1
0134         self.p2 = p2
0135         self.th0 = th0
0136         self.th2 = th2
0137         self.th = th
0138         self.missfrac = 1. - missrat.r
0139         self.rats = rats
0140
0141     def plot_misser(self):
0142         """
0143         6 percent, 3d plot shows on sides of inner block
0144         hmm the triangulated face has a crack down the diagonal
0145         may getting thru the gap ?
0146         so their reflection not in intended place
0147
0148         distinct structure, with islands, somehow inversely related to the
0149         position of cracks in geometry
0150
0151         YEP: confirmed, easily avoided
0152         by shifting the focus point to avoid cracks
0153         """
0154         miss = self.miss
0155         p0m = self.p0[miss]  # destined to miss
0156         p1m = self.p1[miss]
0157
0158         fig = plt.figure()
0159         scatter3d(fig, p0m)
0160         fig.show()
0161
0162
0163 class ReflectionPlot(object):
0164     def __init__(self, s, p, fr):
0165         self.s = s
0166         self.p = p
0167         self.fr = fr
0168         self.dom = fr.dom
0169
0170     def ratio(self, xsel, xall):
0171         dom = self.dom
0172         fa, ba = np.histogram(xall, bins=dom)
0173         fs, bs = np.histogram(xsel, bins=dom)
0174         assert np.all( ba == dom )
0175         assert np.all( bs == dom )
0176         rat = fs.astype(np.float32)/fa.astype(np.float32)
0177         return rat
0178
0179     def abs_norm(self):
0180         """
0181         Absolute normalization using initial "all" distrib with no selection
0182         """
0183         dom = self.dom
0184         s = self.s
0185         p = self.p
0186
0187         sr = self.ratio( s.th, s.tha )
0188         pr = self.ratio( p.th, p.tha )
0189
0190         plt.plot( dom[:-1], sr, drawstyle="steps", label="s", c="r")
0191         plt.plot( dom[:-1], pr, drawstyle="steps", label="p", c="b")
0192
0193     def theory(self):
0194         fr = self.fr
0195         if fr is None:
0196             return
0197         fr.pl()
0198         fr.angles()
0199
0200     def title(self):
0201         msize = self.s.evt.msize()
0202         return "npy-/reflection.py %3.1fM: %s " % (msize,self.fr.title())
0203
0204     def legend(self, ax, log_=False):
0205         if log_:
0206             ax.set_yscale('log')
0207             loc = 'lower left'
0208         else:
0209             loc = 'upper left'
0210         legend = ax.legend(loc=loc, shadow=True)
0211
0212     def plot(self, ax, log_=False):
0213         self.abs_norm()
0214         self.theory()
0215         self.legend(ax, log_=log_)
0216
0217
0218 def oneplot(fr,rp, log_=False):
0219     fig = plt.figure()
0220     plt.title(rp.title())
0221     ax = fig.add_subplot(111)
0222     rp.plot(ax,log_=log_)
0223     fig.show()
0224
0225 def twoplot(fr,rp):
0226     fig = plt.figure()
0227     plt.title(rp.title())
0228     ax = fig.add_subplot(121)
0229     rp.plot(ax,log_=False)
0230     ax = fig.add_subplot(122)
0231     rp.plot(ax,log_=True)
0232     fig.show()
0233
0234 def attic():
0235     transform = "0.500,0.866,0.000,0.000,-0.866,0.500,0.000,0.000,0.000,0.000,1.000,0.000,-86.603,0.000,0.000,1.000"
0236     tx = np.fromstring(transform, sep=",").reshape(4,4)
0237     focus = np.dot([0,0,0,1],tx)[:3]
0238     normal = np.dot([0,1,0,0],tx)[:3]
0239
0240
0241
0242 if __name__ == '__main__':
0243
0244     args = opticks_main(det="reflect",stag="1", ptag="2")
0245
0246     try:
0247         es = Evt(tag=args.stag, label="S", det=args.det, args=args)
0248         ep = Evt(tag=args.ptag, label="P", det=args.det, args=args)
0249     except IOError as err:
0250         log.fatal(err)
0251         sys.exit(args.mrc)
0252
0253     log.info(" es : %s " % es.brief )
0254     log.info(" ep : %s " % ep.brief )
0255
0256
0257     if not (es.valid and ep.valid):
0258         log.fatal("both es and ep must be valid")
0259         sys.exit(1)
0260     pass
0261
0262
0263     normal = [0,0,-1]
0264     source = [0,0,-200]
0265
0266     plt.ion()
0267     fig = plt.figure()
0268     recpos_plot(fig, [es,ep], origin=source)
0269
0270
0271
0272     fig = plt.figure()
0273     angle_plot(fig, [es,ep], axis=normal, origin=source)
0274
0275     swl = es.unique_wavelength()
0276     pwl = ep.unique_wavelength()
0277     assert swl == pwl
0278
0279     #boundary = Boundary("Vacuum///GlassSchottF2")
0280     boundary = Boundary("Vacuum///MainH2OHale")
0281
0282     wl = swl
0283     n1 = boundary.omat.refractive_index(wl)
0284     n2 = boundary.imat.refractive_index(wl)
0285
0286     dom = np.linspace(0,90,250)
0287     fr = Fresnel(n1, n2, dom=dom)
0288
0289     s = Reflect(es, focus=source, normal=normal)
0290     p = Reflect(ep, focus=source, normal=normal)
0291
0292     rp = ReflectionPlot(s, p, fr)
0293
0294     #oneplot(fr,rp,log_=False)
0295     oneplot(fr,rp,log_=True)
0296     #twoplot(fr,rp)
0297
0298
0299
0300