eic-opticks/ana/ana.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 ana.py geometrical and plotting utils
0023 ========================================
0024
0025 TODO: reposition these into more appropriate locations
0026
0027 """
0028 import os, logging, sys
0029 import numpy as np
0030
0031 from opticks.ana.base import opticks_main
0032 from opticks.ana.nbase import count_unique, vnorm
0033 from opticks.ana.evt import Evt, costheta_
0034
0035 deg = np.pi/180.
0036
0037 log = logging.getLogger(__name__)
0038
0039 X,Y,Z,W = 0,1,2,3
0040
0041 def theta(xyz):
0042     """
0043     :param xyz: array of cartesian coordinates
0044     :return: array of Spherical Coordinare polar angle, theta in degrees
0045
0046     First subtract off any needed translations to align
0047     the coordinate system with focal points.
0048
0049     Spherical Coordinates (where theta is polar angle 0:pi, phi is azimuthal 0:2pi)
0050
0051      x = r sin(th) cos(ph)  = r st cp
0052      y = r sin(th) sin(ph)  = r st sp
0053      z = r cos(th)          = r ct
0054
0055
0056      sqrt(x*x + y*y) = r sin(th)
0057                   z  = r cos(th)
0058
0059      atan( sqrt(x*x+y*y) / z ) = th
0060
0061     """
0062     #r = np.linalg.norm(xyz, ord=2, axis=1)
0063     r = vnorm(xyz)
0064     z = xyz[:,2]
0065     th = np.arccos(z/r)*180./np.pi
0066     return th
0067
0068
0069 def scatter3d(fig,  xyz):
0070     ax = fig.add_subplot(111, projection='3d')
0071     ax.scatter(xyz[:,0], xyz[:,1], xyz[:,2])
0072
0073 def histo(fig,  vals):
0074     ax = fig.add_subplot(111)
0075     ax.hist(vals, bins=91,range=[0,90])
0076
0077 def xyz3d(fig, path):
0078     xyz = np.load(path).reshape(-1,3)
0079     ax = fig.add_subplot(111, projection='3d')
0080     ax.scatter(xyz[:,0], xyz[:,1], xyz[:,2])
0081
0082
0083
0084
0085 class Rat(object):
0086     def __init__(self, n, d, label=""):
0087         n = len(n)
0088         d = len(d)
0089         r = float(n)/float(d)
0090         self.n = n
0091         self.d = d
0092         self.r = r
0093         self.label = label
0094     def __repr__(self):
0095         return "Rat %s %s/%s %5.3f " % (self.label,self.n, self.d, self.r)
0096
0097
0098 def recpos_plot(fig, evts, irec=0, nb=100, origin=[0,0,0] ):
0099
0100     origin = np.asarray(origin)
0101
0102     nr = len(evts)
0103     nc = 3
0104     clab = ["X","Y","Z"]
0105
0106     for ir,evt in enumerate(evts):
0107         pos = evt.rpost_(irec)[:,:3] - origin
0108
0109         for ic, lab in enumerate(clab):
0110             ax = fig.add_subplot(nr,nc,1+ic+nc*ir)
0111             ax.hist(pos[:,ic],bins=nb)
0112             ax.set_xlabel(lab)
0113
0114
0115 def angle_plot(fig, evts, irec=0, axis=[0,0,1], origin=[0,0,-200], nb=100):
0116
0117     origin = np.asarray(origin)
0118     nc = len(evts)
0119     nr = 1
0120
0121     for ic,evt in enumerate(evts):
0122         pos = evt.rpost_(irec)[:,:3] - origin
0123
0124         axis_ = np.tile(axis, len(pos)).reshape(-1, len(axis))
0125
0126         ct = costheta_(pos, axis_)
0127         th = np.arccos(ct)/deg
0128
0129         ax = fig.add_subplot(nr,nc,1+ic)
0130         ax.hist(th, bins=nb)
0131         ax.set_xlabel("angle to axis %s " % str(axis))
0132
0133
0134
0135
0136 if __name__ == '__main__':
0137     args = opticks_main(tag="5", det="rainbow", src="torch")
0138     try:
0139         evt = Evt(tag=args.tag, det=args.det, src=args.src, args=args)
0140     except IOError as err:
0141         log.fatal(err)
0142         sys.exit(args.mrc)
0143
0144
0145