g4cx/tests/G4CXTest.py

0001 #!/usr/bin/env python
0002 """
0003 G4CXTest.py
0004 ============
0005
0006 PICK=AB HSEL="TO BT BT SA" ~/o/G4CXTest_GEOM.sh dna
0007
0008
0009 """
0010 import os, logging, numpy as np
0011 from opticks.ana.fold import Fold, IsRemoteSession
0012 from opticks.sysrap.sevt import SEvt, SAB
0013 from opticks.ana.p import *  # including cf boundary___
0014
0015 log = logging.getLogger(__name__)
0016
0017 GLOBAL = int(os.environ.get("GLOBAL","0")) == 1
0018 D = int(os.environ.get("D","3"))
0019 MODE = int(os.environ.get("MODE",D))
0020
0021 print(" D:%d MODE:%d " % (D, MODE))
0022
0023
0024 SEL = os.environ.get("SEL","")
0025 POI = int(os.environ.get("POI","-1"))
0026
0027 PICK = os.environ.get("PICK","AB")
0028 AIDX = int(os.environ.get("AIDX","0"))
0029 BIDX = int(os.environ.get("BIDX","0"))
0030
0031 H,V = 0,2  # X, Z
0032
0033
0034 if IsRemoteSession():  # HMM: maybe do this inside pvplt ?
0035     MODE = 0
0036     print("detect fold.IsRemoteSession forcing MODE:%d" % MODE)
0037 elif MODE in [2,3]:
0038     from opticks.ana.pvplt import *  # HMM this import overrides MODE, so need to keep defaults the same
0039 pass
0040
0041
0042
0043
0044 def onephotonplot(pl, e ):
0045     """
0046     :param pl: MODE 2/2 plotter objects
0047     :param e: SEvt instance
0048     """
0049     if e is None: return
0050     if e.pid < 0: return
0051
0052     print("onephotonplot e.pid %d " % e.pid )
0053
0054     if not hasattr(e,'l'): return
0055     if e.l is None: return
0056     print("onephotonplot e.pid %d PROCEED " % e.pid )
0057
0058     rpos =  e.l[:,:3] + e.off
0059
0060     if MODE == 2:
0061         fig, axs = pl
0062         assert len(axs) == 1
0063         ax = axs[0]
0064         if True:
0065             mpplt_add_contiguous_line_segments(ax, rpos, axes=(H,V), label=None )
0066             #self.mp_plab(ax, f)
0067         pass
0068         #if "nrm" in f.opt:
0069         #    self.mp_a_normal(ax, f)
0070         #pass
0071     elif MODE == 3:
0072         pass
0073         pvplt_add_contiguous_line_segments(pl, rpos )
0074     pass
0075
0076
0077
0078
0079
0080
0081
0082 if __name__ == '__main__':
0083     logging.basicConfig(level=logging.INFO)
0084     print("GLOBAL:%d MODE:%d SEL:%s PICK:%s " % (GLOBAL,MODE, SEL, PICK))
0085
0086     a = SEvt.Load("$AFOLD", symbol="a")
0087     print(repr(a))
0088     b = SEvt.Load("$BFOLD", symbol="b")
0089     print(repr(b))
0090
0091
0092     print(cf)
0093
0094
0095     #sli="[:]"  allowing everything makes for big tables of low stat histories
0096     sli="[:15]"
0097     at = None if a is None else a.minimal_qtab(sli=sli)
0098     bt = None if b is None else b.minimal_qtab(sli=sli)
0099
0100     print("at\n",at)
0101     print("bt\n",bt)
0102
0103     if "SAB" in os.environ:
0104         ab = SAB(a,b) if not a is None and not b is None else None
0105         print(repr(ab))
0106     else:
0107         print(" ab = SAB(a,b) ## SKIPPED as SAB not on os.environ ")
0108     pass
0109
0110
0111     assert PICK in ["A","B","AB","BA", "CF"]
0112     if PICK == "A":
0113         ee = [a,]
0114     elif PICK == "B":
0115         ee = [b,]
0116     elif PICK == "AB":
0117         ee = [a,b,]
0118     elif PICK == "BA":
0119         ee = [b,a,]
0120     elif PICK == "CF":
0121         ee = []
0122     pass
0123
0124     if len(ee) == 0:
0125         log.fatal("set PICK envvar to one of : A,B,AB,BA in order to load SEvt from AFOLD BFOLD " )
0126         assert(0)
0127     pass
0128
0129
0130     context = "PICK=%s MODE=%d SEL=%s POI=%d ~/opticks/g4cx/tests/G4CXTest.sh ana " % (PICK, MODE, SEL, POI )
0131     print(context)
0132
0133     hsel = None
0134     #hsel = "TO BT BT BT SD,TO BT BT BT SA"
0135     HSEL = os.environ.get("HSEL", hsel)
0136
0137
0138     for e in ee:
0139         ew = e.q_startswith_or_(HSEL) if not HSEL is None else None
0140
0141         elabel = "%s : %s " % ( e.symbol.upper(), e.f.base )
0142         if not ew is None:
0143             elabel += " HSEL=%s " % HSEL
0144         pass
0145         label = context + " ## " + elabel
0146
0147         qtab = e.minimal_qtab()
0148
0149
0150         if MODE in [0,1]:
0151             print("not plotting as MODE %d in environ" % MODE )
0152         elif MODE == 2:
0153             pl = mpplt_plotter(label=label)
0154             fig, axs = pl
0155             assert len(axs) == 1
0156             ax = axs[0]
0157
0158             xlim, ylim = mpplt_focus_aspect()
0159             if not xlim is None:
0160                 ax.set_xlim(xlim)
0161                 ax.set_ylim(ylim)
0162             else:
0163                 log.info("mpplt_focus_aspect not enabled, use eg FOCUS=0,0,100 to enable ")
0164             pass
0165
0166             #ax.set_xlim(-356,356)
0167             #ax.set_ylim(-201,201)
0168
0169         elif MODE == 3:
0170             pl = pvplt_plotter(label)
0171             pvplt_viewpoint(pl)   # sensitive to EYE, LOOK, UP envvars
0172             pvplt_frame(pl, e.f.sframe, local=not GLOBAL )
0173         pass
0174
0175         #pp = e.f.inphoton[:,0,:3]
0176         #pp = e.f.photon[:,0,:3]
0177         #pp = e.f.record[:,:,0,:3]
0178
0179         if POI == -1:
0180             pp = e.f.record[:,:,0,:3]
0181         elif POI > -1:
0182             pp = e.f.record[:,POI,0,:3]
0183         else:
0184             pp = None
0185         pass
0186
0187         ww = pp[ew] if not ew is None else None
0188
0189         ppf = pp.reshape(-1,3)
0190         print("ppf.shape : %s " % repr(ppf.shape) )
0191
0192         g_pos = np.ones( [len(ppf), 4 ] )
0193         g_pos[:,:3] = ppf
0194         l_pos = np.dot( g_pos, e.f.sframe.w2m )
0195         u_pos = g_pos if GLOBAL else l_pos
0196
0197
0198
0199         wwf = ww.reshape(-1,3) if not ww is None else None
0200         print("wwf.shape : %s " % repr(wwf.shape) )
0201
0202         if not wwf is None:
0203             h_pos = np.ones( [len(wwf), 4 ] )
0204             h_pos[:,:3] = wwf
0205             i_pos = np.dot( h_pos, e.f.sframe.w2m )
0206             v_pos = h_pos if GLOBAL else i_pos
0207         else:
0208             h_pos = None
0209             i_pos = None
0210             v_pos = None
0211         pass
0212
0213
0214
0215         if SEL == "1":
0216             sel = np.logical_and( np.abs(u_pos[:,H]) < 500, np.abs(u_pos[:,V]) < 500 )
0217             s_pos = u_pos[sel]
0218         else:
0219             s_pos = None
0220         pass
0221
0222
0223         if s_pos is None:
0224             print("s_pos:None skip plotting : use SEL=1 to plot ")
0225         else:
0226             if MODE == 2:
0227                 ax.scatter( s_pos[:,H], s_pos[:,V], s=0.1 )
0228             elif MODE == 3:
0229                 pl.add_points(s_pos[:,:3])
0230             else:
0231                 pass
0232             pass
0233         pass
0234
0235         if not v_pos is None:
0236             print("v_pos:NOT-None proceed with selected point plotting")
0237             if MODE == 2:
0238                 ax.scatter( v_pos[:,H], v_pos[:,V], s=0.1, c="r" )
0239             elif MODE == 3:
0240                 if len(v_pos) > 0:
0241                     pl.add_points(v_pos[:,:3], color="red")
0242                 else:
0243                     print(" SKIP HSEL PLOTTING : pl.add_points(v_pos[:,:3] : AS v_pos EMPTY " )
0244                 pass
0245             else:
0246                 pass
0247             pass
0248         else:
0249             print("v_pos:None skip selected point plotting")
0250             pass
0251         pass
0252
0253
0254         if e.pid > -1:
0255             onephotonplot(pl, e)
0256         pass
0257
0258
0259
0260         if MODE == 2:
0261             fig.show()
0262         elif MODE == 3:
0263             pl.show()
0264         pass
0265     pass    ## ee loop
0266 pass