g4cx/tests/G4CXTest_GEOM.py

0001 #!/usr/bin/env python
0002 """
0003 G4CXTest_GEOM.py
0004 ======================
0005
0006 This script is intended for simple plotting and
0007 statistical SEvt comparison.
0008
0009 The alternate script G4CXTest.py is used for
0010 more detailed plotting such as single photon path
0011 illustration.
0012
0013 ::
0014
0015      PICK=AB HSEL="TO BT BT SA" ~/o/G4CXTest_GEOM.sh ana
0016
0017
0018 """
0019 import os, logging, textwrap, numpy as np
0020 from opticks.ana.fold import Fold
0021 from opticks.sysrap.sevt import SEvt, SAB
0022 from opticks.ana.nbase import chi2_pvalue
0023
0024 from opticks.ana.p import *  # including cf boundary___
0025
0026
0027 MODE = int(os.environ.get("MODE","3"))
0028 SEL = int(os.environ.get("SEL","0"))
0029 HSEL = os.environ.get("HSEL","")    # History selection
0030 PICK = os.environ.get("PICK","A")
0031 SAB_ = "SAB" in os.environ
0032
0033 H,V = 0,2  # X, Z
0034
0035
0036
0037 if MODE in [2,3]:
0038     try:
0039         import pyvista as pv
0040         from opticks.ana.pvplt import pvplt_plotter, pvplt_viewpoint, mpplt_plotter, mpplt_focus_aspect
0041     except ImportError:
0042         pv = None
0043     pass
0044 else:
0045     pv = None
0046 pass
0047
0048
0049 def eprint(expr, _locals, _globals ):
0050     print(expr)
0051     try:
0052        val = eval(expr,_locals,_globals)
0053     except AttributeError:
0054        val = "eprint:AttributeError"
0055     pass
0056     print(val)
0057 pass
0058
0059
0060 if __name__ == '__main__':
0061     logging.basicConfig(level=logging.INFO)
0062     print("MODE:%s SEL:%s PICK:%s SAB_:%s " % (MODE, SEL, PICK, SAB_) )
0063
0064     if PICK == "A" or PICK == "AB" or SAB_:
0065         a = SEvt.Load("$AFOLD", symbol="a")
0066         print(repr(a))
0067     else:
0068         a = None
0069     pass
0070
0071     if PICK == "B" or PICK == "AB" or SAB_:
0072         b = SEvt.Load("$BFOLD", symbol="b")
0073         print(repr(b))
0074     else:
0075         b = None
0076     pass
0077
0078
0079     if not HSEL == "":
0080         print("doing HSEL [%s] history selection " % HSEL )
0081
0082         wa = a.q_startswith_or(HSEL) if not a is None else None
0083         wa_shape = "-" if wa is None else repr(wa.shape)
0084         print("wa.shape %s " % wa_shape)
0085
0086         wb = b.q_startswith_or(HSEL) if not b is None else None
0087         wb_shape = "-" if wb is None else repr(wb.shape)
0088         print("wb.shape %s " % wb_shape)
0089
0090         ra = a.f.record[wa] if not a is None else None
0091         ra_shape = "-" if ra is None else repr(ra.shape)
0092         print("ra.shape %s " % ra_shape)
0093
0094         rb = b.f.record[wb] if not b is None else None
0095         rb_shape = "-" if rb is None else repr(rb.shape)
0096         print("rb.shape %s " % rb_shape)
0097
0098         ns = len(HSEL.split(" "))
0099         st = slice(0,ns)
0100         bnd_ra = ra[:,st,3,0].view(np.uint32) >> 16 if not ra is None else None
0101         bnd_rb = rb[:,st,3,0].view(np.uint32) >> 16 if not rb is None else None
0102         assert np.all( bnd_rb == 0 ) # unfortunately no boundary info for B, YET
0103
0104         lbnd_ra = ra[:,st.stop-1,3,0].view(np.uint32) >> 16 if not ra is None else None
0105         u_lbnd_ra, n_lbnd_ra = np.unique(lbnd_ra, return_counts=True)
0106
0107         for i in range(len(u_lbnd_ra)):
0108             print(" u_lbnd_ra[%2d] %3d   n_lbnd_ra[%2d] %7d   cf.sim.bndnamedict.get(%3d) : %s " %
0109                      ( i, u_lbnd_ra[i], i, n_lbnd_ra[i], u_lbnd_ra[i], cf.sim.bndnamedict.get(u_lbnd_ra[i], "err-bnd") ))
0110             pass
0111         pass
0112
0113
0114         rra = np.sqrt(np.sum(ra[:,st,0,:3]*ra[:,st,0,:3],axis=2)) if not ra is None else None
0115         rra_shape = "-" if rra is None else repr(rra.shape)
0116         print("rra.shape %s " % rra_shape )
0117
0118         rrb = np.sqrt(np.sum(rb[:,st,0,:3]*rb[:,st,0,:3],axis=2)) if not rb is None else None
0119         rrb_shape = "-" if rrb is None else repr(rrb.shape)
0120         print("rrb.shape %s " % rrb_shape )
0121
0122     else:
0123         print("set HSEL to make history selection ")
0124         pass
0125     pass
0126
0127     ee = {'A':a, 'B':b}
0128
0129     if SAB_:
0130         print("[--- ab = SAB(a,b) ----")
0131         ab = None if a is None or b is None else SAB(a,b)
0132         print("]--- ab = SAB(a,b) ----")
0133
0134         print("[----- repr(ab) ")
0135         print(repr(ab))
0136         print("]----- repr(ab) ")
0137     else:
0138         print("ab = SAB(a,b) ## skipped as SAB not in os.environ")
0139     pass
0140
0141
0142     # EXPR_ = r"""
0143     # np.c_[np.unique(a.q, return_counts=True)]
0144     # np.c_[np.unique(b.q, return_counts=True)]
0145     # """
0146     # EXPR = list(filter(None,textwrap.dedent(EXPR_).split("\n")))
0147     # for expr in EXPR:eprint(expr, locals(), globals() )
0148
0149     context = "PICK=%s MODE=%d  ~/opticks/g4cx/tests/G4CXTest_GEOM.sh " % (PICK, MODE )
0150     print(context)
0151
0152     HITONLY = "HITONLY" in os.environ
0153
0154     for Q in PICK:
0155         e = ee.get(Q,None)
0156         if e is None:continue
0157
0158         elabel = "%s : %s " % ( e.symbol.upper(), e.f.base )
0159         label = context + " ## " + elabel
0160
0161
0162         if hasattr(e.f, 'hit'):
0163             hit = e.f.hit[:,0,:3]
0164         else:
0165             hit = None
0166         pass
0167
0168         if hasattr(e.f, 'photon') and not HITONLY:
0169             pos = e.f.photon[:,0,:3]
0170         else:
0171             pos = None
0172         pass
0173
0174         if hasattr(e.f, 'record') and not HITONLY:
0175             sel = np.where(e.f.record[:,:,2,3] > 0) # select on wavelength to avoid unfilled zeros
0176             poi = e.f.record[:,:,0,:3][sel]
0177         else:
0178             poi = None
0179         pass
0180
0181         if MODE == 2:
0182             pl = mpplt_plotter(label=label)
0183             fig, axs = pl
0184             assert len(axs) == 1
0185             ax = axs[0]
0186
0187             xlim, ylim = mpplt_focus_aspect()
0188             if not xlim is None:
0189                 ax.set_xlim(xlim)
0190                 ax.set_ylim(ylim)
0191             else:
0192                 log.info("mpplt_focus_aspect not enabled, use eg FOCUS=0,0,100 to enable ")
0193             pass
0194
0195             if SEL == 1:
0196                 #sel = np.logical_and( np.abs(u_pos[:,H]) < 500, np.abs(u_pos[:,V]) < 500 )
0197                 sel = pos[:,V] > 28000
0198
0199                 u_pos = pos[sel]
0200             else:
0201                 u_pos = pos
0202             pass
0203
0204             ax.scatter( u_pos[:,H], u_pos[:,V], s=0.1 )
0205
0206             fig.show()
0207
0208         elif MODE == 3 and not pv is None:
0209             pl = pvplt_plotter(label)
0210             pvplt_viewpoint(pl) # sensitive EYE, LOOK, UP, ZOOM envvars eg EYE=0,-3,0
0211
0212             if not poi is None:
0213                 pl.add_points( poi, color="green", point_size=3.0 )
0214             pass
0215             if not pos is None:
0216                 pl.add_points( pos, color="red", point_size=3.0 )
0217             pass
0218             if not hit is None:
0219                 pl.add_points( hit, color="cyan", point_size=3.0 )
0220             pass
0221
0222             if not "NOGRID" in os.environ:
0223                 pl.show_grid()
0224             pass
0225             cp = pl.show()
0226         pass
0227     pass
0228
0229