eic-opticks/ana/tpmt.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.py : PmtInBox Opticks vs G4 History comparisons
0023 ==========================================================
0024
0025 Loads test events from Opticks and Geant4 and
0026 compares their bounce histories.
0027
0028 Create the events by running tpmt- bash functions.
0029
0030 The convention is adopted of using positive tags for Opticks
0031 and negative ones of the same magnitude for the corresponding
0032 Geant4 simulated event.
0033
0034
0035 See Also
0036 -----------
0037
0038 :doc:`tpmt_debug`
0039       simulation debugging notes to acheive Opticks Geant4 match
0040
0041 :doc:`tpmt_distrib`
0042       comparison of distributions
0043
0044
0045 Expected Output
0046 ------------------
0047
0048 The expected output from the test is shown below.
0049 History step abbreviation:
0050
0051 * *TO* torch step
0052 * *BT* boundary transmit
0053 * *BR* boundary reflect
0054 * *SA* surface absorb
0055 * *SD* surface detect
0056 * *AB* bulk absorb
0057 * *SC* bulk scatter
0058
0059 Material abbreviations:
0060
0061 * *MO* Mineral Oil
0062 * *Py* Pyrex
0063 * *Vm* Vacuum
0064 * *OV* Opaque Vacuum
0065
0066
0067 .. code-block:: py
0068
0069     delta:ana blyth$ i
0070     SQLITE3_DATABASE=/usr/local/env/nuwa/mocknuwa.db
0071     Python 2.7.11 (default, Dec  5 2015, 23:51:51)
0072     Type "copyright", "credits" or "license" for more information.
0073
0074     IPython 1.2.1 -- An enhanced Interactive Python.
0075     ?         -> Introduction and overview of IPython's features.
0076     %quickref -> Quick reference.
0077     help      -> Python's own help system.
0078     object?   -> Details about 'object', use 'object??' for extra details.
0079
0080     IPython profile: g4opticks
0081
0082     In [1]: run pmt_test.py
0083      1.175 100.000
0084      0.377 100.002
0085                           4:PmtInBox   -4:PmtInBox           c2
0086                      8cd         67948        68252             0.68  [3 ] TO BT SA
0087                      7cd         21648        21369             1.81  [3 ] TO BT SD
0088                     8ccd          4581         4539             0.19  [4 ] TO BT BT SA
0089                       4d          3794         3864             0.64  [2 ] TO AB
0090                      86d           640          617             0.42  [3 ] TO SC SA
0091                      4cd           444          427             0.33  [3 ] TO BT AB
0092                     4ccd           350          362             0.20  [4 ] TO BT BT AB
0093                      8bd           283          259             1.06  [3 ] TO BR SA
0094                     8c6d            81           84             0.05  [4 ] TO SC BT SA
0095                    86ccd            51           57             0.33  [5 ] TO BT BT SC SA
0096                   8cbbcd            36           53             3.25  [6 ] TO BT BR BR BT SA
0097                      46d            40           30             1.43  [3 ] TO SC AB
0098                     7c6d            20           28             1.33  [4 ] TO SC BT SD
0099                      4bd            28           21             1.00  [3 ] TO BR AB
0100                 8cbc6ccd             9            3             0.00  [8 ] TO BT BT SC BT BR BT SA
0101                     866d             8            4             0.00  [4 ] TO SC SC SA
0102                    8cc6d             7            7             0.00  [5 ] TO SC BT BT SA
0103                     86bd             6            4             0.00  [4 ] TO BR SC SA
0104                     8b6d             3            6             0.00  [4 ] TO SC BR SA
0105               cbccbbbbcd             4            0             0.00  [10] TO BT BR BR BR BR BT BT BR BT
0106                               100000       100000         0.91
0107                           4:PmtInBox   -4:PmtInBox           c2
0108                      ee4         90040        90048             0.00  [3 ] MO Py Py
0109                     44e4          4931         4901             0.09  [4 ] MO Py MO MO
0110                       44          3794         3864             0.64  [2 ] MO MO
0111                      444           991          927             2.14  [3 ] MO MO MO
0112                     ee44           101          113             0.67  [4 ] MO MO Py Py
0113                    444e4            52           58             0.33  [5 ] MO Py MO MO MO
0114                   44eee4            40           54             2.09  [6 ] MO Py Py Py MO MO
0115                     4444            17           14             0.29  [4 ] MO MO MO MO
0116                    44e44             8            7             0.00  [5 ] MO MO Py MO MO
0117                 44ee44e4             6            3             0.00  [8 ] MO Py MO MO Py Py MO MO
0118                 444e44e4             5            0             0.00  [8 ] MO Py MO MO Py MO MO MO
0119               44e4eeeee4             4            0             0.00  [10] MO Py Py Py Py Py MO Py MO MO
0120                   ee44e4             0            4             0.00  [6 ] MO Py MO MO Py Py
0121                    ee444             2            0             0.00  [5 ] MO MO MO Py Py
0122               44edbe44e4             2            0             0.00  [10] MO Py MO MO Py OV Vm Py MO MO
0123                   4444e4             0            2             0.00  [6 ] MO Py MO MO MO MO
0124               4ebdbe44e4             0            1             0.00  [10] MO Py MO MO Py OV Vm OV Py MO
0125               4e5dbe44e4             0            1             0.00  [10] MO Py MO MO Py OV Vm Bk Py MO
0126               eebdbe44e4             1            0             0.00  [10] MO Py MO MO Py OV Vm OV Py Py
0127                  44ee444             1            0             0.00  [7 ] MO MO MO Py Py MO MO
0128                               100000       100000         0.78
0129
0130
0131
0132
0133
0134 """
0135 from __future__ import print_function
0136 import os, sys, logging, argparse
0137 log = logging.getLogger(__name__)
0138 try:
0139     import numpy as np
0140 except ImportError:
0141     np = None
0142 pass
0143
0144 if np == None:
0145     print("no numpy early abort : so this does not cause SSysTest to fail ")
0146     sys.exit(0)
0147 pass
0148
0149
0150 from opticks.ana.base import opticks_main
0151 from opticks.ana.nbase import vnorm
0152 from opticks.ana.evt  import Evt
0153
0154
0155
0156 if __name__ == '__main__':
0157
0158
0159
0160
0161     np.set_printoptions(precision=4, linewidth=200)
0162
0163     args = opticks_main(doc=__doc__, tag="10", src="torch", det="PmtInBox", c2max=2.0, tagoffset=0)
0164
0165     log.info("tag %s src %s det %s c2max %s  " % (args.utag,args.src,args.det, args.c2max))
0166
0167
0168     #seqs = ["TO BT BR BT BT BT BT SA"]
0169     #seqs = ["TO BT BR BR BT SA"]
0170     seqs=[]
0171
0172     try:
0173         a = Evt(tag="%s" % args.utag, src=args.src, det=args.det, seqs=seqs, args=args)
0174         b = Evt(tag="-%s" % args.utag , src=args.src, det=args.det, seqs=seqs, args=args)
0175     except IOError as err:
0176         log.fatal(err)
0177         #sys.exit(args.mrc)  this causes a sysrap-t test fail from lack of a tmp file
0178         sys.exit(0)
0179
0180
0181
0182     log.info( " a : %s " % a.brief)
0183     log.info( " b : %s " % b.brief )
0184
0185 if 0:
0186     if a.valid:
0187         a0 = a.rpost_(0)
0188         #a0r = np.linalg.norm(a0[:,:2],2,1)
0189         a0r = vnorm(a0[:,:2])
0190         if len(a0r)>0:
0191             print(" ".join(map(lambda _:"%6.3f" % _, (a0r.min(),a0r.max()))))
0192
0193     if b.valid:
0194         b0 = b.rpost_(0)
0195         #b0r = np.linalg.norm(b0[:,:2],2,1)
0196         b0r = vnorm(b0[:,:2])
0197         if len(b0r)>0:
0198             print(" ".join(map(lambda _:"%6.3f" % _, (b0r.min(),b0r.max()))))
0199
0200 if 1:
0201
0202     c2max = args.c2max
0203     c2max = None
0204
0205     Evt.compare_table(a,b, "seqhis_ana seqmat_ana".split(), lmx=20, c2max=c2max, cf=False)
0206     Evt.compare_table(a,b, "seqhis_ana seqhis_ana_2 seqhis_ana_3 seqhis_ana_4 seqmat_ana".split(), lmx=20, c2max=c2max, cf=True)
0207     Evt.compare_table(a,b, "pflags_ana hflags_ana".split(), lmx=20, c2max=c2max )
0208
0209     #a.history_table()
0210     #b.history_table()
0211
0212