eic-opticks/ana/tpmt_skimmer.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 tpmt_skimmer.py: Following positions of PMT skimmers
0023 ======================================================
0024
0025 Creates plot showing step by step average positions of
0026 all photons with a specific history, namely: "TO BT BR BR BT SA"
0027 and tabulates the min/max/mid positions.
0028
0029 Expected Output
0030 -----------------
0031
0032 .. code-block:: py
0033
0034     In [1]: run tpmt_skimmer.py
0035     WARNING:opticks.ana.evt:init_index PmtInBox/torch/-5 : TO BT BR BR BT SA finds too few (ps)phosel uniques : 1
0036     WARNING:opticks.ana.evt:init_index PmtInBox/torch/-5 : TO BT BR BR BT SA finds too few (rs)recsel uniques : 1
0037     WARNING:opticks.ana.evt:init_index PmtInBox/torch/-5 : TO BT BR BR BT SA finds too few (rsr)reshaped-recsel uniques : 1
0038     A(Op) PmtInBox/torch/5 : TO BT BR BR BT SA
0039       0 z:    300.000    300.000    300.000   r:     98.999     98.999     98.999   t:      0.098      0.098      0.098   smry m1/m2   4/ 14 MO/Py  -28 ( 27)  13:TO
0040       1 z:     67.559     67.559     67.559   r:     98.999     98.999     98.999   t:      1.251      1.251      1.251   smry m1/m2  14/  4 Py/MO   28 ( 27)  12:BT
0041       2 z:     50.832     50.832     50.832   r:    100.372    100.372    100.372   t:      1.331      1.331      1.331   smry m1/m2  14/ 11 Py/OV -125 (124)  11:BR
0042       3 z:     35.551     35.551     35.551   r:     93.176     93.176     93.176   t:      1.416      1.416      1.416   smry m1/m2  14/  4 Py/MO   28 ( 27)  11:BR
0043       4 z:     19.181     19.181     19.181   r:     89.001     89.001     89.001   t:      1.495      1.495      1.495   smry m1/m2   4/ 12 MO/Rk  124 (123)  12:BT
0044       5 z:   -300.000   -300.000   -300.000   r:     26.569     26.569     26.569   t:      3.107      3.107      3.107   smry m1/m2   4/ 12 MO/Rk  124 (123)   8:SA
0045     B(G4) PmtInBox/torch/-5 : TO BT BR BR BT SA
0046       0 z:    300.000    300.000    300.000   r:     98.999     98.999     98.999   t:      0.098      0.098      0.098   smry m1/m2   4/  0 MO/?0?    0 ( -1)  13:TO
0047       1 z:     67.559     67.559     67.559   r:     98.999     98.999     98.999   t:      1.251      1.251      1.251   smry m1/m2  14/  0 Py/?0?    0 ( -1)  12:BT
0048       2 z:     50.832     50.832     50.832   r:    100.372    100.372    100.372   t:      1.331      1.331      1.331   smry m1/m2  14/  0 Py/?0?    0 ( -1)  11:BR
0049       3 z:     35.551     35.551     35.551   r:     93.176     93.176     93.176   t:      1.416      1.416      1.416   smry m1/m2  14/  0 Py/?0?    0 ( -1)  11:BR
0050       4 z:     19.181     19.181     19.181   r:     89.001     89.001     89.001   t:      1.495      1.495      1.495   smry m1/m2   4/  0 MO/?0?    0 ( -1)  12:BT
0051       5 z:   -300.000   -300.000   -300.000   r:     26.569     26.569     26.569   t:      3.107      3.107      3.107   smry m1/m2   4/  0 MO/?0?    0 ( -1)   8:SA
0052
0053
0054 See Also
0055 ------------
0056
0057 :doc:`pmt_skimmer_debug`
0058        Debugging Opticks TIR with pmt_skimmer.py
0059
0060
0061 """
0062
0063 import os, sys, logging, numpy as np
0064 log = logging.getLogger(__name__)
0065
0066 import matplotlib.pyplot as plt
0067
0068 from opticks.ana.base import opticks_main
0069 from opticks.ana.evt import Evt
0070 from opticks.ana.pmt.plot import Pmt, one_plot
0071
0072
0073 X,Y,Z,W = 0,1,2,3
0074
0075 ZX = [Z,X]
0076 ZY = [Z,Y]
0077 XY = [X,Y]
0078
0079 ZMIN,ZMAX,ZAVG,RMIN,RMAX,RAVG,TMIN,TMAX,TAVG = 0,1,2,3,4,5,6,7,8
0080
0081
0082 def zr_plot(data, neg=False):
0083     xx = data[:,ZAVG]
0084     yy = data[:,RAVG]
0085     if neg:yy=-yy
0086
0087     labels = map(str, range(len(data)))
0088     plt.plot(xx, yy, "o")
0089     for label, x, y in zip(labels, xx,yy):
0090         plt.annotate(label, xy = (x, y), xytext = (-3,10), textcoords = 'offset points')
0091
0092
0093
0094 if __name__ == '__main__':
0095     args = opticks_main(tag="10", src="torch", det="PmtInBox")
0096     np.set_printoptions(precision=4, linewidth=200)
0097
0098     plt.ion()
0099     plt.close()
0100
0101     fig = plt.figure()
0102
0103     pmt = Pmt()
0104     ALL, PYREX, VACUUM, CATHODE, BOTTOM, DYNODE = None,0,1,2,3,4
0105     pts = pmt.parts(ALL)
0106     one_plot(fig, pmt, pts, axes=ZX, clip=True)
0107
0108     pol = False
0109
0110     # aseqs=["TO BT BR BT BT BT BT BT BT SA"]   before fixing the TIR bug this was what was happening
0111     aseqs=["TO BT BR BR BT SA"]
0112     bseqs=["TO BT BR BR BT SA"]
0113     na = len(aseqs[0].split())
0114     nb = len(bseqs[0].split())
0115
0116     try:
0117         a = Evt(tag="%s" % args.tag, src=args.src, det=args.det, seqs=aseqs, args=args)
0118         b = Evt(tag="-%s" % args.tag, src=args.src, det=args.det, seqs=bseqs, args=args)
0119     except IOError as err:
0120         log.fatal(err)
0121         sys.exit(args.mrc)
0122
0123     if a.valid:
0124         log.info("A : plot zrt_profile %s " % a.label)
0125         print "A(Op) %s " % a.label
0126         a_zrt = a.zrt_profile(na, pol=pol)
0127         zr_plot(a_zrt)
0128     else:
0129         log.warning("failed to load A")
0130
0131     if b.valid:
0132         log.info("B : plot zrt_profile %s " % b.label)
0133         print "B(G4) %s " % b.label
0134         b_zrt = b.zrt_profile(nb, pol=pol)
0135         zr_plot(b_zrt, neg=True)
0136     else:
0137         log.warning("failed to load B")
0138
0139
0140
0141
0142