eic-opticks/ana/ckmplot.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 ::
0023
0024     In [14]: evt.rpost_(slice(0,5)).shape
0025     Out[14]: (500000, 5, 4)
0026
0027
0028 """
0029 import os, sys, re, logging, numpy as np
0030 from collections import OrderedDict as odict
0031
0032 import matplotlib as mpl
0033 from mpl_toolkits.mplot3d import Axes3D
0034 import matplotlib.pyplot as plt
0035
0036 log = logging.getLogger(__name__)
0037
0038 from opticks.ana.base import opticks_main, gcp_
0039 from opticks.ana.evt import Evt
0040
0041
0042
0043 def plot(evt, fig, rposti=False, rpostn=True, ox=False, mesh=True):
0044     """
0045     This will not work for millions of photons like Opticks geometry shaders do
0046     but handy anyway for an alternative interface to the record points
0047     of a few hundred photons
0048     """
0049     ax = fig.gca(projection='3d')
0050
0051     if mesh:
0052         ptn = re.compile("\d+")
0053         meshidxs = map(int, filter(lambda n:ptn.match(n), os.listdir(gcp_("GMeshLib"))))
0054         for idx in meshidxs:
0055             vtx = np.load(gcp_("GMeshLib/%d/vertices.npy" % idx))
0056             ax.scatter( vtx[:,0], vtx[:,1], vtx[:,2])
0057         pass
0058     pass
0059
0060     if rposti:
0061         for i in range(min(300, len(evt.seqhis))):
0062             xyzt = evt.rposti(i)
0063             x = xyzt[:,0]
0064             y = xyzt[:,1]
0065             z = xyzt[:,2]
0066             ax.plot(x, y, z )
0067         pass
0068     pass
0069
0070     if rpostn:
0071         for n in range(16):
0072             rpn = evt.rpostn(n)
0073             if len(rpn) == 0:continue
0074             for i in range(n-1):
0075                 x = rpn[:,i,0]
0076                 y = rpn[:,i,1]
0077                 z = rpn[:,i,2]
0078                 u = rpn[:,i+1,0] - rpn[:,i,0]
0079                 v = rpn[:,i+1,1] - rpn[:,i,1]
0080                 w = rpn[:,i+1,2] - rpn[:,i,2]
0081                 ax.quiver(x, y, z, u, v, w )
0082             pass
0083         pass
0084
0085     if ox:
0086         xs = evt.ox[:,0,0]
0087         ys = evt.ox[:,0,1]
0088         zs = evt.ox[:,0,2]
0089         ax.scatter( xs, ys, zs)
0090     pass
0091
0092
0093
0094 if __name__ == '__main__':
0095     args = opticks_main(tag="1",src="natural", det="g4live", doc=__doc__)
0096     np.set_printoptions(suppress=True, precision=3)
0097
0098     evt = Evt(tag=args.tag, src=args.src, det=args.det, seqs=[], args=args)
0099
0100     log.debug("evt")
0101     print(evt)
0102
0103     log.debug("evt.history_table")
0104     evt.history_table(slice(0,20))
0105     log.debug("evt.history_table DONE")
0106
0107
0108     plt.ion()
0109     fig = plt.figure()
0110
0111     plot(evt, fig)
0112
0113     plt.show()
0114
0115
0116
0117
0118
0119