eic-opticks/sysrap/SABTest.py

0001 #!/usr/bin/env python
0002 """
0003 SABTest.py : analysis script used from SABTest.sh
0004 ==================================================
0005
0006 Started from ~/j/InputPhotonsCheck/InputPhotonsCheck.py
0007
0008 """
0009
0010 import os, logging, numpy as np
0011 log = logging.getLogger(__name__)
0012 os.environ["MODE"] = "3"
0013
0014 from opticks.sysrap.sevt import SEvt, SAB, SABHit
0015 from opticks.ana.p import cf
0016
0017 MODE = int(os.environ.get("MODE","3"))
0018 MODE0 = MODE
0019 print("SABTest.py initial MODE:%d " % MODE)
0020
0021
0022 GLOBAL = int(os.environ.get("GLOBAL","0")) == 1
0023 PICK = os.environ.get("PICK","A")
0024 INCPOI = float(os.environ.get("INCPOI","0"))   ## increment pyvista presentation point and line size, often need -5
0025
0026
0027 if MODE in [2,3]:
0028     import opticks.ana.pvplt as pvp
0029     print(pvp.__doc__)
0030     print("after import MODE:%d " % MODE)
0031     MODE = MODE0
0032     # HMM this import overrides MODE, so need to keep defaults the same
0033 pass
0034
0035
0036 if __name__ == '__main__':
0037     logging.basicConfig(level=logging.INFO)
0038
0039     INST_FRAME = int(os.environ.get("INST_FRAME","-1"))
0040     if INST_FRAME == -1:
0041         print("INST_FRAME %d : using default sframe saved with the sevt" % INST_FRAME )
0042         M2W_OVERRIDE = None
0043         W2M_OVERRIDE = None
0044     else:
0045         print("INST_FRAME %d : W2M_OVERRIDE obtained from cf.inst controlled by envvar INST_FRAME " % INST_FRAME )
0046         M2W_OVERRIDE = np.eye(4)
0047         M2W_OVERRIDE[:,:3] = cf.inst[INST_FRAME][:,:3]
0048         W2M_OVERRIDE = np.linalg.inv(M2W_OVERRIDE)
0049     pass
0050
0051
0052     print(__doc__)
0053
0054     a = SEvt.Load("$AFOLD", symbol="a", W2M=W2M_OVERRIDE)
0055     b = SEvt.Load("$BFOLD", symbol="b", W2M=W2M_OVERRIDE)
0056     print(repr(a))
0057     print(repr(b))
0058
0059
0060     if not a is None and b is None:
0061         PICK = "A"
0062     elif not b is None and a is None:
0063         PICK = "B"
0064     else:
0065         pass
0066     pass
0067
0068     if not a is None and not b is None:
0069         if "SAB" in os.environ:
0070             print("[--- ab = SAB(a,b) ----")
0071             ab = SAB(a,b)
0072             print("]--- ab = SAB(a,b) ----")
0073
0074             print("[----- repr(ab) ")
0075             print(repr(ab))
0076             print("]----- repr(ab) ")
0077         else:
0078             print(" NOT doing ab = SAB(a,b) # do that with : export SAB=1" )
0079         pass
0080
0081         print("[--- abh = SABHit(a,b) ----")
0082         abh = SABHit(a,b)
0083         print("]--- abh = SABHit(a,b) ----")
0084         print("[----- repr(abh) ")
0085         print(repr(abh))
0086         print("]----- repr(abh) ")
0087     pass
0088
0089     assert PICK in ["A","B","AB","BA", "CF"]
0090     if PICK == "A":
0091         ee = [a,]
0092     elif PICK == "B":
0093         ee = [b,]
0094     elif PICK == "AB":
0095         ee = [a,b,]
0096     elif PICK == "BA":
0097         ee = [b,a,]
0098     elif PICK == "CF":
0099         ee = []
0100     pass
0101
0102     context = "PICK=%s MODE=%d  ~/j/jtds/jtds.sh " % (PICK, MODE )
0103     print(context)
0104
0105     pl = None
0106     for e in ee:
0107         if e is None:continue
0108         elabel = "%s : %s " % ( e.symbol.upper(), e.f.base )
0109
0110         if hasattr(e.f, 'photon'):
0111             pos = e.f.photon[:,0,:3]
0112         elif hasattr(e.f, 'hit'):
0113             pos = e.f.hit[:,0,:3]
0114         else:
0115             log.info("%s:sevt lacks photon or hit" % e.symbol)
0116             pos = None
0117         pass
0118
0119         if hasattr(e.f, 'record'):
0120             sel = np.where(e.f.record[:,:,2,3] > 0) # select on wavelength to avoid unfilled zeros
0121             poi = e.f.record[:,:,0,:3][sel]
0122         else:
0123             log.info("%s:sevt lacks record" % e.symbol)
0124             poi = None
0125         pass
0126
0127         if hasattr(e.f, 'sframe'):
0128             _W2M = e.f.sframe.w2m
0129             log.info("%s:sevt has sframe giving _W2M " % e.symbol)
0130         else:
0131             _W2M = np.eye(4)
0132             log.info("%s:sevt lacks sframe" % e.symbol)
0133         pass
0134
0135         W2M = _W2M if W2M_OVERRIDE is None else W2M_OVERRIDE
0136
0137
0138         print("W2M_OVERRIDE (from manual INST_FRAME, no longer typical)\n",W2M_OVERRIDE)
0139         print("_W2M(from %s.f.sframe.w2m, now standard)\n" % e.symbol,_W2M)
0140         print("W2M\n",W2M)
0141
0142         if not pos is None:
0143             print("transform gpos (from photon or hit array) to lpos using W2M : ie into target frame : GLOBAL %d " % GLOBAL)
0144             gpos = np.ones( [len(pos), 4 ] )
0145             gpos[:,:3] = pos
0146             lpos = np.dot( gpos, W2M )   # hmm unfilled global zeros getting transformed somewhere
0147             upos = gpos if GLOBAL else lpos
0148         else:
0149             upos = None
0150         pass
0151
0152         if not poi is None:
0153             print("transform gpoi (from record array) to lpoi using W2M : ie into target frame : GLOBAL %d " % GLOBAL)
0154             gpoi = np.ones( [len(poi), 4 ] )
0155             gpoi[:,:3] = poi
0156             lpoi = np.dot( gpoi, W2M )
0157             upoi = gpoi if GLOBAL else lpoi
0158         else:
0159             upoi = None
0160         pass
0161
0162
0163         label = context + " ## " + elabel
0164
0165         if MODE == 3 and not "SAB" in os.environ:
0166             if pvp.pv is None:
0167                 print("MODE is 3 but pvp.pv is None, get into ana environment first eg with hookup_conda_ok")
0168             else:
0169                 pl = pvp.pvplt_plotter(label)
0170                 pvp.pvplt_viewpoint(pl) # sensitive EYE, LOOK, UP, ZOOM envvars eg EYE=0,-3,0
0171
0172                 if not upoi is None:
0173                     print("add_points green for the upoi : record points ")
0174                     pl.add_points( upoi[:,:3], color="green", point_size=3.0 )
0175                 pass
0176                 if not upos is None:
0177                     print("add_points red  for the upos : photon/hit points ")
0178                     pl.add_points( upos[:,:3], color="red", point_size=3.0 )
0179                 pass
0180                 pl.show_grid()
0181                 pl.increment_point_size_and_line_width(INCPOI)
0182                 cp = pl.show()
0183             pass
0184         pass
0185     pass
0186
0187
0188
0189 pass
0190
0191