eic-opticks/ana/hit.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 hit.py
0024 ========
0025
0026 Checking relationship betweet hit counts and photon counts
0027 from a photon scan.  The "hits" are not real ones, just some
0028 photon flag mask chosen to give some statistics for machinery testing.
0029
0030 From 1M up to 60M get a very uniform gradient, from 70M continue
0031 with the same gradient but with a great big offset, as if suddenly
0032 lost around 67M photons.
0033
0034 Sawtooth plot::
0035
0036     In [9]: np.unique(n, axis=0)
0037     Out[9]:
0038     array([[  1000000,       516],
0039            [ 10000000,      5212],
0040            [ 20000000,     10406],
0041            [ 30000000,     15675],
0042            [ 40000000,     20802],
0043            [ 50000000,     26031],
0044            [ 60000000,     31126],
0045            [ 70000000,      1562],    ## glitch down to a hit count would expect from around 3M photons
0046            [ 80000000,      6663],
0047            [ 90000000,     11963],
0048            [100000000,     17122]], dtype=int32)
0049
0050
0051     In [18]:  u[:,0]/u[:,1]
0052     Out[18]:
0053     array([ 1937,  1918,  1921,  1913,  1922,  1920,  1927, 44814, 12006,
0054             7523,  5840], dtype=int32)
0055
0056     In [20]:  u[:7,0]/u[:7,1]
0057     Out[20]: array([1937, 1918, 1921, 1913, 1922, 1920, 1927], dtype=int32)
0058
0059     In [21]: np.average( u[:7,0]/u[:7,1])
0060     Out[21]: 1922.5714285714287
0061
0062     ## gradient same beyond the glitch
0063
0064     In [37]: (u[8:,0]-u[7,0])/(u[8:,1]-u[7,1])
0065     Out[37]: array([1960, 1922, 1928], dtype=int32)
0066
0067
0068 """
0069 import os, sys, re, logging, numpy as np
0070 from collections import OrderedDict as odict
0071
0072 from opticks.ana.num import Num
0073 from opticks.ana.base import findfile
0074
0075 import matplotlib.pyplot as plt
0076
0077 log = logging.getLogger(__name__)
0078
0079
0080 class Arr(object):
0081
0082     @classmethod
0083     def Find(cls, pfx, name, qcat=None):
0084         base = os.path.expandvars(os.path.join("$TMP", pfx ))
0085         rpaths = findfile(base, name )
0086         aa = []
0087         for rp in rpaths:
0088             path = os.path.join(base, rp)
0089             elem = path.split("/")
0090             cat = elem[elem.index("evt")+1]
0091             if not qcat is None and not cat.startswith(qcat): continue
0092             snpho = cat.split("_")[-1]
0093             npho = Num.Int(snpho)
0094             aa.append(Arr(path,cat, npho))
0095         pass
0096         return aa
0097
0098     def __init__(self, path, cat, npho):
0099         self.path = path
0100         self.cat = cat
0101         self.npho = npho
0102         self.a = np.load(path)
0103         self.items = len(self.a)
0104
0105
0106 if __name__ == '__main__':
0107
0108
0109     np.set_printoptions(suppress=True, precision=3)
0110
0111     plt.ion()
0112
0113     aa = Arr.Find("scan-ph", "ht.npy", "cvd_1_rtx_1_" )
0114     print(" aa %s " % len(aa))
0115
0116     n = np.zeros( [len(aa),2] , dtype=np.int32 )
0117     for i,a in enumerate(sorted(aa, key=lambda a:a.npho)):
0118         n[i,0] = a.npho
0119         n[i,1] = a.items
0120     pass
0121     print(n)
0122
0123     fig = plt.figure()
0124     ax = fig.add_subplot(111)
0125
0126     ax.plot( n[:,0], n[:,1], "o--" )
0127
0128     plt.show()
0129
0130
0131
0132
0133
0134