eic-opticks/ana/genstep.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 genstep.py: Fit genstep xyz vs time, used for viewpoint tracking
0023 =================================================================
0024
0025 See okc/TrackView.cc for the usage of this
0026
0027
0028 Fit the genstep xyz vs time to obtain parametric eqn of genstep position
0029 with time parameter.::
0030
0031     In [3]: tk = A.load_("cerenkov","1_track","dayabay")
0032
0033     In [4]: tk
0034     Out[4]: A(cerenkov,1_track,dayabay)
0035
0036     In [5]: print tk
0037     [[ -16390.518 -802295.938   -7059.101]
0038      [   -162.573     251.993       0.172]]
0039
0040 ::
0041
0042     In [1]: run genstep.py
0043     [[[   0.177   -1.583    4.94     1.   ]
0044       [-252.339  -45.677 -155.278    0.   ]
0045       [   0.      82.83     0.       0.   ]]]
0046     INFO:opticks.ana.nload:saving derivative of A(cerenkov,1,juno) to /usr/local/env/opticks/juno/cerenkov/1_track.npy
0047
0048 ::
0049
0050     simon:npy blyth$ /usr/local/opticks.ana/bin/NumpyEvtTest
0051     [2016-Mar-25 12:08:16.965734]:info: NumpyEvt::loadGenstepDerivativeFromFile   typ cerenkov tag 1_track det dayabay
0052     [2016-Mar-25 12:08:16.966326]:info: NumpyEvt::loadGenstepDerivativeFromFile (3,4)
0053
0054     (  0)  -16390.518  -802295.938   -7059.101       1.000
0055     (  1)    -162.573     251.993       0.172       0.000
0056     (  2)       0.844      27.423       0.000       0.000
0057
0058 """
0059 import os, logging
0060 import numpy as np
0061 import matplotlib.pyplot as plt
0062
0063 from opticks.ana.base import opticks_main
0064 from opticks.ana.nload import A, I, II, path_
0065
0066 log = logging.getLogger(__name__)
0067
0068 X,Y,Z,W,T = 0,1,2,3,3
0069
0070
0071 if __name__ == '__main__':
0072
0073     args = opticks_main(det="juno", src="cerenkov", tag="1")
0074
0075     try:
0076         a = A.load_("gensteps",args.src,args.tag,args.det)
0077     except IOError as err:
0078         log.fatal(err)
0079         sys.exit(args.mrc)
0080
0081
0082     log.info("loaded gensteps %s %s %s " % (a.path, a.stamp, repr(a.shape)))
0083
0084     #path = os.path.expandvars("$LOCAL_BASE/opticks/opticksdata/gensteps/dayabay/cerenkov/1.npy")
0085     #path = os.path.expandvars("$LOCAL_BASE/opticks/opticksdata/gensteps/juno/cerenkov/1.npy")
0086     #a = np.load(path)
0087
0088
0089     xyzt = a[:,1]
0090     #print xyzt
0091     x,y,z,t = xyzt[:,X], xyzt[:,Y], xyzt[:,Z], xyzt[:,T]
0092
0093     tr = [t.min(), t.max()]
0094     xr = [x.min(), x.max()]
0095     yr = [y.min(), y.max()]
0096     zr = [z.min(), z.max()]
0097
0098
0099     plt.close()
0100     plt.ion()
0101     ny,nx = 1,3
0102
0103     fig = plt.figure()
0104
0105     ax = fig.add_subplot(ny,nx,1)
0106     ax.scatter(t, x)
0107     xf = np.polyfit(t,x,1,full=True)
0108     xm, xc = xf[0]        ## xm is "x**1" coefficient, xc is "x**0" (constant term)
0109     xl = [xm*tt + xc for tt in tr]
0110     ax.plot(tr, xl, '-r')
0111
0112     ax = fig.add_subplot(ny,nx,2)
0113     ax.scatter(t, y)
0114     yf = np.polyfit(t,y,1,full=True)
0115     ym, yc = yf[0]
0116     yl = [ym*tt + yc for tt in tr]
0117     ax.plot(tr, yl, '-r')
0118
0119     ax = fig.add_subplot(ny,nx,3)
0120     ax.scatter(t, z)
0121     zf = np.polyfit(t,z,1,full=True)
0122     zm, zc = zf[0]
0123     zl = [zm*tt + zc for tt in tr]
0124     ax.plot(tr, zl, '-r')
0125
0126
0127     ## origin, direction, range
0128     track = np.array(
0129                 [[
0130                   [xc,yc,zc,1.0],
0131                   [xm,ym,zm,0.0],
0132                   [tr[0],tr[1],0.0,0.0]
0133                 ]], dtype=np.float32)
0134     print track
0135
0136
0137     #a.derivative_save(track, "track")
0138
0139
0140
0141