eic-opticks/sysrap/sevt.py

0001 #!/usr/bin/env python
0002 """
0003 sevt.py (formerly opticks.u4.tests.U4SimulateTest)
0004 ====================================================
0005
0006
0007 """
0008 import os, re, logging, textwrap, numpy as np
0009 log = logging.getLogger(__name__)
0010
0011 from opticks.ana.fold import Fold
0012
0013 print("[from opticks.ana.p import * ")
0014 from opticks.ana.p import *
0015 print("]from opticks.ana.p import * ")
0016
0017 from opticks.ana.eget import eslice_
0018 from opticks.ana.base import PhotonCodeFlags
0019 from opticks.ana.qcf import QU,QCF,QCFZero
0020 from opticks.ana.npmeta import NPMeta
0021 from opticks.ana.hismask import HisMask
0022
0023 hm = HisMask()
0024
0025 pcf = PhotonCodeFlags()
0026 fln = pcf.fln
0027 fla = pcf.fla
0028
0029
0030 MODE = int(os.environ.get("MODE","2"))
0031 if MODE > 0:
0032     from opticks.ana.pvplt import *
0033 else:
0034     pass
0035 pass
0036
0037 QLIM=eslice_("QLIM", "0:10")
0038
0039
0040 axes = 0, 2  # X,Z
0041 H,V = axes
0042
0043
0044 class SEvt(object):
0045     """
0046     Higher level wrapper for an Opticks SEvt folder of arrays
0047
0048     WIP: Concatenation of multiple SEvt, starting from Fold concatenation
0049     """
0050     @classmethod
0051     def Load(cls, *args, **kwa):
0052         NEVT = int(kwa.get("NEVT",0))
0053         log.info("SEvt.Load NEVT:%d " % NEVT)
0054         if NEVT > 0:
0055             f = cls.LoadConcat(*args,**kwa)
0056         else:
0057             f = Fold.Load(*args, **kwa )
0058         pass
0059         return None if f is None else cls(f, **kwa)
0060
0061     @classmethod
0062     def LoadConcat(cls, *args, **kwa):
0063         """
0064         HMM maybe should load the separate SEvt and then concat those,
0065         rather than trying to concat at the lower Fold level
0066         """
0067         NEVT = int(kwa.get("NEVT",0))
0068         assert NEVT > 0
0069         f = Fold.LoadConcat(*args, **kwa)
0070         assert hasattr(f, 'ff')
0071         return f
0072
0073     def __init__(self, f, **kwa):
0074         """
0075         :param f: Fold instance
0076         :param kwa: override param, only W2M accepted
0077         """
0078
0079         self.f = f
0080         self.symbol = f.symbol
0081         self._r = None
0082         self.r = None
0083         self._a = None
0084         self.a = None
0085         self._pid = -1
0086         self.W2M = kwa.get("W2M", None)
0087         print("sevt.init W2M\n", self.W2M )
0088
0089         self.init_run(f)
0090         self.init_meta(f)
0091         self.init_fold_meta_timestamp(f)
0092         self.init_fold_meta(f)
0093         self.init_U4R_names(f)
0094         self.init_photon(f)
0095         self.init_record(f)
0096         self.init_record_sframe(f)
0097         self.init_SEventConfig_meta(f)
0098         self.init_seq(f)
0099         self.init_aux(f)
0100         self.init_ee(f)    ## handling of concatenated SEvt is done currently only for ee
0101         self.init_aux_t(f)
0102         self.init_sup(f)
0103         self.init_junoSD_PMT_v2_SProfile(f)
0104         self.init_env(f)  # setting pid needs frames so has to be late
0105
0106     @classmethod
0107     def CommonRunBase(cls, f):
0108         """
0109         :param f: Fold instance
0110         :return urun_base: str
0111         """
0112         run_bases = []
0113         if type(f.base) is str:
0114             run_base = os.path.dirname(f.base)
0115             run_bases.append(run_base)
0116         elif type(f.base) is list:
0117             for base in f.base:
0118                 run_base = os.path.dirname(base)
0119                 if not run_base in run_bases:
0120                     run_bases.append(run_base)
0121                 pass
0122         else:
0123             pass
0124         pass
0125         assert len(run_bases) == 1
0126         urun_base = run_bases[0]
0127         return urun_base
0128
0129     def init_run(self, f):
0130         """
0131         HMM the run meta should be same for all of concatenated SEvts
0132         """
0133         urun_base = self.CommonRunBase(f)
0134         urun_path = os.path.join( urun_base, "run.npy" )
0135
0136         # THIS WAS RECURSING AND LOADING THE n001 and p001 folders
0137         #print("[ sevt.init_run Fold.Load urun_path %s loading urun_base %s " % (urun_path,urun_base) )
0138         #fp = Fold.Load(urun_base) if os.path.exists(urun_path) else None
0139         #print("] sevt.init_run Fold.Load ")
0140         #run_meta = not fp is None and getattr(fp,'run_meta', None) != None
0141         #self.fp = fp
0142
0143         urun_meta = os.path.join( urun_base, "run_meta.txt" )
0144         run_meta = NPMeta.Load(urun_meta) if os.path.exists(urun_meta) else None
0145         has_run_meta = not run_meta is None
0146
0147         prof2stamp_ = lambda prof:prof.split(",")[0]
0148
0149         rr_ = [prof2stamp_(run_meta.SEvt__BeginOfRun),
0150                prof2stamp_(run_meta.SEvt__EndOfRun  )] if has_run_meta else [0,0]
0151
0152         rr = np.array(rr_, dtype=np.uint64 )
0153
0154         self.rr = rr
0155         self.run_meta = run_meta
0156
0157         qwns = ["FAKES_SKIP",]
0158         for qwn in qwns:
0159             mval = list(map(int,getattr(run_meta, qwn, ["-1"] ) if has_run_meta else ["-2"]))
0160             setattr(self, qwn, mval[0] )
0161         pass
0162
0163
0164     def init_env(self, f):
0165         """
0166         """
0167         pid = int(os.environ.get("%sPID" % f.symbol.upper(), "-1"))
0168         opt = os.environ.get("%sOPT" % f.symbol.upper(), "")
0169         off_ = os.environ.get("%sOFF" % f.symbol.upper(), "0,0,0")
0170         off = np.array(list(map(float,off_.split(","))))
0171         log.info("SEvt.__init__  symbol %s pid %d opt %s off %s " % (f.symbol, pid, opt, str(off)) )
0172
0173         self.pid = pid
0174         self.opt = opt
0175         self.off = off
0176
0177     def init_meta(self, f):
0178         """
0179         """
0180         metakey = os.environ.get("METAKEY", "junoSD_PMT_v2_Opticks_meta" )
0181         meta = getattr(f, metakey, None)
0182         self.metakey = metakey
0183         self.meta = meta
0184
0185     def init_fold_meta_timestamp(self, f):
0186         """
0187         """
0188         if getattr(f, "NPFold_meta", None) == None: return
0189         msrc = f.NPFold_meta
0190         if msrc is None: return
0191
0192         def mts_(key):
0193             val = msrc.get_value(key)
0194             if val == '': val = "0"
0195             return np.uint64(val)
0196
0197         bor = mts_("T_BeginOfRun")
0198         boe = mts_("t_BeginOfEvent")
0199         gs0 = mts_("t_setGenstep_0")
0200         gs1 = mts_("t_setGenstep_1")
0201         gs2 = mts_("t_setGenstep_2")
0202         gs3 = mts_("t_setGenstep_3")
0203         gs4 = mts_("t_setGenstep_4")
0204         gs5 = mts_("t_setGenstep_5")
0205         gs6 = mts_("t_setGenstep_6")
0206         gs7 = mts_("t_setGenstep_7")
0207         gs8 = mts_("t_setGenstep_8")
0208         pre = mts_("t_PreLaunch")
0209         pos = mts_("t_PostLaunch")
0210         eoe = mts_("t_EndOfEvent")
0211
0212         b2e = eoe - boe
0213         ee = np.array([boe, eoe, b2e], dtype=np.uint64 )
0214         ett = np.array([bor, boe, gs0, gs1, gs2, gs3, gs4, gs5, gs6, gs7, gs8, pre, pos, eoe], dtype=np.uint64 )
0215         ettl = "bor boe gs0 gs1 gs2 gs3 gs4 gs5 gs6 gs7 gs8 pre pos eoe".split()
0216         ettd = np.zeros( len(ett), dtype=np.uint64 )
0217         ettd[1:] = np.diff(ett)
0218
0219         ettv = ett.view("datetime64[us]")
0220         ettc = np.c_[ettl, ett, ettv.astype("|S26"), ettd/1e6, ettd]
0221
0222         self.ee = ee
0223         self.ett = ett
0224         self.ettd = ettd
0225         self.ettv = ettv
0226         self.ettl = ettl
0227         self.ettc = ettc
0228
0229
0230     def init_fold_meta(self, f):
0231         """
0232         """
0233         if getattr(f, "NPFold_meta", None) == None: return
0234         msrc = f.NPFold_meta
0235         if msrc is None: return
0236
0237         GEOM = msrc.get_value("GEOM")
0238         CHECK = msrc.get_value("CHECK")
0239         LAYOUT = msrc.get_value("LAYOUT")
0240         VERSION = msrc.get_value("VERSION")
0241         if LAYOUT.find(" ") > -1: LAYOUT = ""
0242
0243         SCRIPT = os.environ.get("SCRIPT", "")
0244
0245         GEOMList = msrc.get_value("${GEOM}_GEOMList", "")
0246         if GEOMList.endswith("GEOMList"): GEOMList = ""
0247
0248         TITLE = "N=%s %s # %s/%s " % (VERSION, SCRIPT, LAYOUT, CHECK )
0249         IDE = list(filter(None,[f.symbol.upper(), GEOM, GEOMList, "N%s"%VERSION, LAYOUT, CHECK]))
0250         ID = "_".join(IDE)
0251
0252         self.CHECK = CHECK
0253         self.LAYOUT = LAYOUT
0254         self.VERSION = VERSION
0255         self.SCRIPT = SCRIPT
0256         self.GEOM = GEOM
0257         self.GEOMList = GEOMList
0258         self.TITLE = TITLE
0259         self.IDE = IDE
0260         self.ID = ID
0261
0262     def init_U4R_names(self, f):
0263         """
0264         Now its an array with UNSET
0265         """
0266         U4R_names = getattr(f, "U4R_names", None)
0267         #SPECS = np.array(U4R_names.lines) if not U4R_names is None else None
0268         SPECS = U4R_names if not U4R_names is None else None
0269         self.SPECS = SPECS
0270
0271     @classmethod
0272     def Make_flagmask_table(cls, flagmask):
0273         """
0274
0275             In [30]: np.c_[ label, fmtab[:,1] ][:20]
0276             Out[30]:
0277             array([['TO|BT|SA', '326185'],
0278                    ['EC|TO|BT|SD', '312177'],
0279                    ['TO|BT|SR|SA', '56805'],
0280                    ['TO|BT|BR|SR|SA', '53334'],
0281                    ['TO|BT|BR|SA', '42888'],
0282                    ['TO|BT|BR|AB', '37716'],
0283                    ['EC|TO|BT|BR|SD', '21709'],
0284                    ['TO|BT|BR|SA|SC', '20683'],
0285                    ['TO|BT|BR|SC|AB', '15776'],
0286                    ['EC|TO|BT|BR|SD|SC', '13679'],
0287                    ['TO|BT|AB', '12490'],
0288                    ['TO|BT|BR|SR|SA|SC', '8668'],
0289
0290         """
0291         _fmtabraw = np.c_[np.unique(flagmask,return_counts=True)]
0292         fmtabraw = _fmtabraw[np.argsort(_fmtabraw[:,1])][::-1]
0293         label = hm.label(fmtabraw[:,0])
0294         fmtab = np.c_[ label, fmtabraw[:,1] ]
0295         return fmtab
0296
0297     def init_photon(self, f):
0298         if hasattr(f,'photon') and not f.photon is None:
0299             iix = f.photon[:,1,3].view(np.int32)
0300             flagmask = f.photon[:,3,3].view(np.int32)
0301             fmtab = self.Make_flagmask_table(flagmask)
0302         else:
0303             iix = None
0304             flagmask = None
0305             fmtab = None
0306         pass
0307         self.iix = iix
0308         self.flagmask = flagmask
0309         self.fmtab = fmtab
0310
0311
0312     def init_record(self, f):
0313         """
0314         :param f: fold instance
0315         """
0316         if hasattr(f,'record') and not f.record is None:
0317
0318             if f.record.dtype.name == 'float32':
0319                 iir = f.record[:,:,1,3].view(np.int32)
0320                 idr = f.record[:,:,1,3].view(np.int32)
0321             elif f.record.dtype.name == 'float64':
0322                 iir = f.record[:,:,1,3].view(np.int64)
0323                 idr = f.record[:,:,1,3].view(np.int64)
0324             else:
0325                 assert False
0326             pass
0327         else:
0328             iir = None
0329             idr = None
0330         pass
0331         self.iir = iir
0332         self.idr = idr
0333
0334     def init_record_sframe(self, f):
0335
0336         with_sframe = not getattr(f, "sframe", None) is None
0337         with_record = not getattr(f, "record", None) is None
0338
0339         W2M = self.W2M
0340
0341         if with_sframe and with_record:
0342             gpos = np.ones( f.record.shape[:-1] )  ## trim last dimension reducing eg (10000,32,4,4) to (10000, 32, 4)
0343             gpos[:,:,:3] = f.record[:,:,0,:3]      ## point positions of all photons
0344             w2m = W2M if not W2M is None else f.sframe.w2m
0345             lpos = np.dot( gpos, w2m )  ## transform all points from global to local frame
0346         else:
0347             w2m = W2M if not W2M is None else None
0348             gpos = None
0349             lpos = None
0350         pass
0351         self.w2m = w2m
0352         self.gpos = gpos
0353         self.lpos = lpos
0354
0355
0356
0357     def init_SEventConfig_meta(self, f):
0358         meta = getattr(f, 'SEventConfig_meta', {})
0359
0360         ipl = getattr(meta, "InputPhoton", [])
0361         ip = ipl[0] if len(ipl)==1 else None
0362
0363         ipfl = getattr(meta, "InputPhotonFrame", [])
0364         ipf = ipfl[0] if len(ipfl)==1 else None
0365
0366         ipcl = getattr(meta, "InputPhotonCheck", [])
0367         ipc = ipc[0] if len(ipcl)==1 else None
0368
0369         self.ip = ip
0370         self.ipf = ipf
0371         self.ipc = ipc
0372
0373
0374     def q_find(self, txt="EC", encoding="utf-8"):
0375         """
0376         :param txt:
0377         :return a: array of indices of self.q array elements that contain *txt*
0378
0379         np.char.find returns -1 when the txt is not found in array elements
0380         so adding one makes that zero which is excluded with np.flatnonzero
0381
0382         Note that because the self.q array elements typically end with blanks
0383         it is not useful to implement q_endswith, instead use this q_find
0384         to find array indices with elements containing the txt.
0385         """
0386         return np.flatnonzero(1+np.char.find(self.q, txt.encode(encoding) ))
0387
0388     def q_find_or(self, txt1="EC", txt2="EX", encoding="utf-8"):
0389         return np.flatnonzero(np.logical_or(
0390                   1+np.char.find(self.q,txt1.encode(encoding)),
0391                   1+np.char.find(self.q,txt2.encode(encoding))
0392                ))
0393
0394     def q_find_or_(self, txts_="EC,EX", delim=",", encoding="utf-8"):
0395         txts= txts_.split(delim)
0396         aa = []
0397         for txt in txts:
0398             aa.append( 1+np.char.find(self.q, txt.encode(encoding)))
0399         pass
0400         return np.flatnonzero(np.logical_or.reduce(aa))
0401
0402     def q__(self, prefix="TO BT SD", encoding="utf-8"):
0403         return self.q_startswith_or_(prefix) if "," in prefix else self.q_startswith(prefix=prefix)
0404
0405     def q_startswith(self, prefix="TO BT SD", encoding="utf-8"):
0406         return np.flatnonzero(np.char.startswith(self.q, prefix.encode(encoding) ))
0407
0408     def q_startswith_or(self, pfx1="TO BT BT BT SA", pfx2="TO BT BT BT SD", encoding="utf-8"):
0409         return np.flatnonzero(np.logical_or(
0410                    np.char.startswith(self.q,pfx1.encode(encoding)),
0411                    np.char.startswith(self.q,pfx2.encode(encoding))
0412                ))
0413
0414     def q_startswith_or_(self, pfxs_="TO BT BT BT SA,TO BT BT BT SD", delim=",", encoding="utf-8"):
0415         """
0416         :param pfxs_: delimited string list with potentially multiple history strings
0417         :param delim:
0418
0419         Another way to implement this is with np.any
0420         """
0421         pfxs = pfxs_.split(delim)
0422         aa = []
0423         for pfx in pfxs:
0424             aa.append( np.char.startswith(self.q, pfx.encode(encoding)))
0425         pass
0426         return np.flatnonzero(np.logical_or.reduce(aa))
0427
0428     def init_seq(self, f):
0429         """
0430
0431         +---------------------+----------------------------------------------------------+
0432         | a.f.seq.shape       | (1000000, 2, 2)                                          |
0433         +---------------------+----------------------------------------------------------+
0434         | a.f.seq[:,0].shape  | (1000000, 2)                                             |
0435         +---------------------+----------------------------------------------------------+
0436         | a.f.seq[0,0]        | array([14748335283153718477, 37762157772], dtype=uint64) |
0437         +---------------------+----------------------------------------------------------+
0438
0439         """
0440         symbol = f.symbol
0441         qlim = QLIM
0442         qtab_ = "np.c_[qn,qi,qu][quo][qlim]"
0443
0444         if hasattr(f,'seq') and not f.seq is None:
0445             q_ = f.seq[:,0]
0446             q  = ht.seqhis(q_)  # ht from opticks.ana.p
0447             qq = ht.Convert(q_)  # (n,32) int8 : for easy access to nibbles
0448             n = np.sum( seqnib_(q_), axis=1 )  # occupied nibbles, ie number of history points
0449
0450             qu, qi, qn = np.unique(q, return_index=True, return_counts=True)
0451             quo = np.argsort(qn)[::-1]
0452
0453             qtab = eval(qtab_)
0454             qtab_ = qtab_.replace("q","%s.q" % symbol)
0455
0456             nosc = np.ones(len(qq), dtype=np.bool_ )       # start all true
0457             nosc[np.where(qq == pcf.SC)[0]] = 0     # knock out photons with scatter in their histories
0458
0459             nore = np.ones(len(qq), dtype=np.bool_ )       # start all true
0460             nore[np.where(qq == pcf.RE)[0]] = 0     # knock out photons with reemission in their histories
0461
0462             noscab = np.ones(len(qq), dtype=np.bool_ )     # start all true
0463             noscab[np.where(qq == pcf.SC)[0]] = 0   # knock out photons with scatter in their histories
0464             noscab[np.where(qq == pcf.AB)[0]] = 0   # knock out photons with bulk absorb  in their histories
0465         else:
0466             q_ = None
0467             q = None
0468             qq = None
0469             n = None
0470             qu, qi, qn = None, None, None
0471             quo = None
0472             qtab = None
0473
0474             nosc = None
0475             nore = None
0476             noscab = None
0477         pass
0478
0479         self.q_ = q_
0480         self.q = q
0481         self.qq = qq
0482         self.n = n
0483
0484         self.qu = qu
0485         self.qi = qi
0486         self.qn = qn
0487
0488         self.quo = quo
0489         self.qtab = qtab
0490
0491         self.qlim = qlim
0492         self.qtab_ = qtab_
0493
0494         self.nosc = nosc   # mask of photons without SC in their histories
0495         self.nore = nore   # mask of photons without RE in their histories
0496         self.noscab = noscab   # mask of photons without SC or AB in their histories
0497
0498
0499
0500     def minimal_qtab(self, sli="[:10]", dump=False):
0501         """
0502         :return qtab: history table in descending count order
0503         """
0504         e = self
0505         uq,iq,nq = np.unique(e.q, return_index=True, return_counts=True)
0506         oq = np.argsort(nq)[::-1]
0507         expr ="np.c_[nq,iq,uq][oq]%(sli)s" % locals()
0508         qtab = eval(expr)
0509         if dump:
0510             log.info("minimal_qtab : %s " % expr)
0511             print(qtab)
0512         pass
0513         return qtab
0514
0515     def init_aux(self, f):
0516         """
0517         Note aux is currently CPU only
0518         """
0519         if hasattr(f, 'aux') and not f.aux is None:
0520             fk = f.aux[:,:,2,2].view(np.uint32)    ## fakemask : for investigating fakes when FAKES_SKIP is disabled
0521             spec_ = f.aux[:,:,2,3].view(np.int32)   ## step spec
0522             max_point = f.aux.shape[1]   # instead of hardcoding typical values like 32 or 10, use array shape
0523             uc4 = f.aux[:,:,2,2].copy().view(np.uint8).reshape(-1,max_point,4) ## see sysrap/spho.h c4/C4Pho.h
0524             eph = uc4[:,:,1]          # .y    ProcessHits enum at step point level
0525             ep = np.max(eph, axis=1 ) #       ProcessHits enum at photon level
0526         else:
0527             fk = None
0528             spec_ = None
0529             uc4 = None
0530             eph = None
0531             ep = None
0532             t = None
0533         pass
0534         self.fk = fk
0535         self.spec_ = spec_
0536         self.uc4 = uc4
0537         self.eph = eph   # ProcessHits EPH enum at step point level
0538         self.ep  = ep    # ProcessHits EPH enum at photon level
0539
0540     def init_ee(self, f):
0541         """
0542         microsecond[us] timestamps (UTC epoch) [BeginOfEvent,EndOfEvent,Begin2End]
0543
0544         For concatenated SEvt the total of all the EndOfEvent-BeginOfEvent
0545         is placed into ee[-1]
0546
0547         TIMING METADATA HAS BEEN MOVED FROM THE PHOTON TO THE FOLD
0548         with_ff = not getattr(f, 'ff', None) is None
0549         log.info("init_ee with_ff:%d" % (with_ff))
0550         if with_ff:
0551             kk = f.ff.keys()
0552             boe = np.uint64(0)
0553             eoe = np.uint64(0)
0554             b2e = np.uint64(0)
0555             for k in kk:
0556                 fk = f.ff[k]
0557                 k_boe = np.uint64(fk.photon_meta.t_BeginOfEvent[0])
0558                 k_eoe = np.uint64(fk.photon_meta.t_EndOfEvent[0])
0559                 k_b2e = np.uint64(k_eoe-k_boe)
0560                 b2e += k_b2e
0561             pass
0562         else:
0563             boe, eoe, b2e = 0,0,0
0564         pass
0565         self.ee = np.array([boe, eoe, b2e], dtype=np.uint64 )
0566         """
0567
0568     def init_junoSD_PMT_v2_SProfile(self, f):
0569         """
0570         The timestamps come from sysrap/sstamp.h and are datetime64[us] (UTC) compliant
0571
0572         pf
0573             uint64_t microsecond timestamps collected by SProfile.h
0574         pfr
0575             uint64_t (last - first) timestamp difference for each SProfile (currently ProcessHits call)
0576
0577
0578         More than 10% of time spent in ProcessHits::
0579
0580             In [16]: np.sum(a.pfr)/a.ee[-1]
0581             Out[16]: 0.10450881266649384
0582
0583             In [17]: np.sum(b.pfr)/b.ee[-1]
0584             Out[17]: 0.11134881257006096
0585
0586         """
0587         if hasattr(f, 'junoSD_PMT_v2_SProfile') and not f.junoSD_PMT_v2_SProfile is None:
0588             pf = f.junoSD_PMT_v2_SProfile
0589             pfmx = np.max(pf[:,1:], axis=1 )
0590             pfmi = pf[:,1]
0591             pfr = pfmx - pfmi
0592         else:
0593             pf = None
0594             pfmx = None
0595             pfmi = None
0596             pfr = None
0597         pass
0598         self.pf = pf  ## CAUTION: multiple ProcessHits calls per photon, so not in photon index order
0599         self.pfmx = pfmx
0600         self.pfmi = pfmi
0601         self.pfr  = pfr
0602
0603
0604
0605     def init_aux_t(self, f):
0606         if hasattr(f, 'aux') and not f.aux is None:
0607             t = f.aux[:,:,3,:2].copy().view(np.uint64).squeeze()   # step point timestamps
0608         else:
0609             t = None
0610         pass
0611         self.t = t      # array of photon step point time stamps (UTC epoch)
0612
0613     def init_sup(self, f):
0614         """
0615         sup is CPU only
0616
0617         photon level beginPhoton endPhoton time stamps from the sup quad4
0618         """
0619         if hasattr(f,'sup') and not f.sup is None:
0620             s0 = f.sup[:,0,:2].copy().view(np.uint64).squeeze()  # SEvt::beginPhoton (xsup.q0.w.x)
0621             s1 = f.sup[:,0,2:].copy().view(np.uint64).squeeze()  # SEvt::finalPhoton (xsup.q0.w.y)
0622             ss = s1 - s0      # endPhoton - beginPhoton
0623
0624             f0 = f.sup[:,1,:2].copy().view(np.uint64).squeeze()  # SEvt::finalPhoton (xsup.q1.w.x) t_PenultimatePoint
0625             f1 = f.sup[:,1,2:].copy().view(np.uint64).squeeze()  # SEvt::finalPhoton (xsup.q1.w.y) t_LastPoint
0626             ff = f1 - f0      # LastPoint - PenultimatePoint
0627
0628             h0 = f.sup[:,2,:2].copy().view(np.uint64).squeeze()  # SEvt::addProcessHitsStamp(0) (xsup.q2.w.x)
0629             h1 = f.sup[:,2,2:].copy().view(np.uint64).squeeze()  # SEvt::addProcessHitsStamp(0) (xsup.q2.w.y)
0630             hh = h1 - h0                      # timestamp range of SEvt::AddProcessHitsStamp(0) calls
0631             hc = f.sup[:,3,0].view(np.int32)
0632
0633             i0 = f.sup[:,4,:2].copy().view(np.uint64).squeeze()  # SEvt::addProcessHitsStamp(1) (xsup.q4.w.x)
0634             i1 = f.sup[:,4,2:].copy().view(np.uint64).squeeze()  # SEvt::addProcessHitsStamp(1) (xsup.q4.w.y)
0635             ii = i1 - i0                      # timestamp range of SEvt::AddProcessHitsStamp(1) calls
0636             ic = f.sup[:,5,0].view(np.int32)
0637
0638             hi0 = i0 - h0
0639             hi1 = i1 - h1
0640         else:
0641             s0 = None
0642             s1 = None
0643             ss = None
0644
0645             f0 = None
0646             f1 = None
0647             ff = None
0648
0649             h0 = None
0650             h1 = None
0651             hh = None
0652             hc = None
0653
0654             i0 = None
0655             i1 = None
0656             ii = None
0657             ic = None
0658
0659             hi0 = None
0660             hi1 = None
0661         pass
0662         if not getattr(self, 't', None) is None and not s0 is None:
0663             t0 = s0.min()
0664             tt = self.t.copy()
0665             tt[np.where( tt != 0 )] -= t0   # subtract earliest time "pedestal", but leave the zeros
0666         else:
0667             t0 = None
0668             tt = None
0669         pass
0670
0671         self.t0 = t0    # scalar : event minimum time stamp (which is minimum s0:beginPhoton timestamp)
0672         self.tt = tt    # array of photon step point time stamps relative to t0
0673
0674         self.s0 = s0    # array of beginPhoton time stamps (UTC epoch)
0675         self.s1 = s1    # array of endPhoton time stamps (UTC epoch)
0676         self.ss = ss    # array of endPhoton-beginPhoton timestamp differences
0677
0678         self.f0 = f0    # array of PenultimatePoint timestamps (UTC epoch)
0679         self.f1 = f1    # array of LastPoint timestamps (UTC epoch)
0680         self.ff = ff    # array of LastPoint-PenultimatePoint timestamp differences
0681
0682         self.h0 = h0    # array of SEvt::AddProcessHitsStamp(0) range begin (UTC epoch)
0683         self.h1 = h1    # array of SEvt::AddProcessHitsStamp(0) range end (UTC epoch)
0684         self.hh = hh    # array of SEvt::AddProcessHitsStamp(0) ranges (microseconds [us])
0685         self.hc = hc    # array of SEvt::AddProcessHitsStamp(0) call counts for each photon
0686
0687         self.i0 = i0    # array of SEvt::AddProcessHitsStamp(1) range begin (UTC epoch)
0688         self.i1 = i1    # array of SEvt::AddProcessHitsStamp(1) range end (UTC epoch)
0689         self.ii = ii    # array of SEvt::AddProcessHitsStamp(1) ranges (microseconds [us])
0690         self.ic = ic    # array of SEvt::AddProcessHitsStamp(1) call counts for each photon
0691
0692         self.hi0 = hi0  # array of range begin SEvt::AddProcessHitsStamp(1)-SEvt::AddProcessHitsStamp(0)
0693         self.hi1 = hi1  # array of range end   SEvt::AddProcessHitsStamp(1)-SEvt::AddProcessHitsStamp(0)
0694
0695     def __repr__(self):
0696         fmt = "SEvt symbol %s pid %s opt %s off %s %s.f.base %s "
0697         return fmt % ( self.symbol, self.pid, self.opt, str(self.off), self.symbol, self.f.base )
0698
0699     def _get_pid(self):
0700         return self._pid
0701     def _set_pid(self, pid):
0702         """
0703         """
0704         f = self.f
0705         q = self.q
0706         W2M = self.W2M
0707         symbol = self.f.symbol
0708         if pid > -1 and hasattr(f,'record') and pid < len(f.record):
0709             _r = f.record[pid]
0710             wl = _r[:,2,3]
0711             r = _r[wl > 0]    ## select wl>0 to avoid lots of record buffer zeros, example shape (13, 4, 4)
0712
0713             g = np.ones([len(r),4])
0714             g[:,:3] = r[:,0,:3]
0715
0716             w2m = W2M if not W2M is None else f.sframe.w2m
0717             l = np.dot( g, w2m )
0718             ld = np.diff(l,axis=0)
0719
0720             pass
0721             _aux = getattr(f, 'aux', None)
0722             _a = None if _aux is None else _aux[pid]
0723             if _a is None:
0724                 a = None
0725                 ast = "no-ast"
0726             else:
0727                 a = _a[wl > 0]
0728                 ast = a[:len(a),1,3].copy().view(np.int8)[::4]
0729             pass
0730             qpid = q[pid][0].decode("utf-8").strip()
0731
0732             npoi = (len(qpid)+1)//3
0733             qkey = " ".join(map(lambda _:"%0.2d"%_, range(npoi)))
0734
0735             #label = "%s %sPID:%d\n%s " % (qpid, symbol.upper(), pid, qkey)
0736             label = "%s\n%s" % (qpid, qkey)
0737         else:
0738             _r = None
0739             r = None
0740             _a = None
0741             a = None
0742             ast = None
0743             qpid = None
0744             label = ""
0745
0746             g = None
0747             l = None
0748             ld = None
0749         pass
0750         self._pid = pid
0751         self._r = _r
0752         self.r = r
0753         self._a = _a
0754         self.a = a
0755         self.ast = ast
0756         self.qpid = qpid
0757         self.label = label
0758         self.g = g
0759         self.l = l
0760         self.ld = ld
0761
0762     pid = property(_get_pid, _set_pid)
0763
0764
0765     #def _get_w2m(self):
0766     #    return self._w2m
0767     #
0768     #def _set_w2m(self, w2m=None):
0769     #    self._w2m = w2m if not w2m is None else f.sframe.w2m
0770
0771
0772
0773
0774     def spec_(self, i):
0775         """
0776         ## need np.abs as detected fakes that are not skipped are negated
0777         In [10]: np.c_[a.spec[5,:a.n[5]],a.SPECS[np.abs(a.spec[5,:a.n[5]])]]
0778         Out[10]:
0779         array([['0', 'UNSET'],
0780                ['1', 'Water/Water:pInnerWater/pLPMT_Hamamatsu_R12860'],
0781                ['2', 'Water/AcrylicMask:pLPMT_Hamamatsu_R12860/HamamatsuR12860pMask'],
0782                ['3', 'AcrylicMask/Water:HamamatsuR12860pMask/pLPMT_Hamamatsu_R12860'],
0783                ['4', 'Water/Pyrex:pLPMT_Hamamatsu_R12860/HamamatsuR12860_PMT_20inch_log_phys'],
0784                ['-5', 'Pyrex/Pyrex:HamamatsuR12860_PMT_20inch_log_phys/HamamatsuR12860_PMT_20inch_body_phys'],
0785                ['6', 'Pyrex/Pyrex:HamamatsuR12860_PMT_20inch_body_phys/HamamatsuR12860_PMT_20inch_body_phys']], dtype='<U94')
0786
0787         In [1]: a.spec_(5)
0788         Out[1]:
0789         array([['0', 'UNSET'],
0790                ['1', 'Water/Water:pInnerWater/pLPMT_Hamamatsu_R12860'],
0791                ['2', 'Water/AcrylicMask:pLPMT_Hamamatsu_R12860/HamamatsuR12860pMask'],
0792                ['3', 'AcrylicMask/Water:HamamatsuR12860pMask/pLPMT_Hamamatsu_R12860'],
0793                ['4', 'Water/Pyrex:pLPMT_Hamamatsu_R12860/HamamatsuR12860_PMT_20inch_log_phys'],
0794                ['-5', 'Pyrex/Pyrex:HamamatsuR12860_PMT_20inch_log_phys/HamamatsuR12860_PMT_20inch_body_phys'],
0795                ['6', 'Pyrex/Pyrex:HamamatsuR12860_PMT_20inch_body_phys/HamamatsuR12860_PMT_20inch_body_phys']], dtype='<U94')
0796
0797         In [2]: b.spec_(5)
0798         Out[2]:
0799         array([['0', 'UNSET'],
0800                ['1', 'Water/Water:pInnerWater/pLPMT_Hamamatsu_R12860'],
0801                ['2', 'Water/AcrylicMask:pLPMT_Hamamatsu_R12860/HamamatsuR12860pMask'],
0802                ['3', 'AcrylicMask/Water:HamamatsuR12860pMask/pLPMT_Hamamatsu_R12860'],
0803                ['4', 'Water/Pyrex:pLPMT_Hamamatsu_R12860/HamamatsuR12860_PMT_20inch_log_phys'],
0804                ['5', 'Pyrex/Vacuum:HamamatsuR12860_PMT_20inch_log_phys/HamamatsuR12860_PMT_20inch_inner_phys']], dtype='<U93')
0805
0806         """
0807         t = self
0808         n = t.n[i]
0809         spec = t.spec[i,:n]
0810         return np.c_[spec,t.SPECS[np.abs(spec)]]
0811
0812
0813 class SABHit(object):
0814     def __init__(self, a, b):
0815         self.a = a
0816         self.b = b
0817         assert not a is None and not b is None
0818
0819     def __repr__(self):
0820         a = len(self.a.f.hit) if hasattr(self.a.f, "hit") else 0
0821         a_std = np.sqrt(a)
0822         a_rel = a_std/a  if a != 0 else 0 # estimate of rel error from sqrt of count, as Poisson
0823
0824         b = len(self.b.f.hit) if hasattr(self.b.f, "hit") else 0
0825         b_std = np.sqrt(b)
0826         b_rel = b_std/b if b != 0 else 0  # estimate of rel error from sqrt of count, as Poisson
0827
0828         ab_rel = a_rel + b_rel # assume uncorrelated : so add relative errors
0829         ba_rel = a_rel + b_rel # assume uncorrelated : so add relative errors
0830
0831         ab = float(a)/float(b) if b != 0 else 0
0832         ab_std = (ab_rel)*ab
0833
0834         ba = float(b)/float(a) if a != 0 else 0
0835         ba_std = (ba_rel)*ba
0836
0837         fmt = "\n\nSABHit  a_nhit : ( %7d +- %4d )  b_nhit:(  %7d +- %4d )  a/b: ( %6.3f +- %6.3f )  b/a: ( %6.3f +- %6.3f)  [assume uncorr.]\n\n"
0838         return fmt % ( a, int(a_std), b, int(b_std),  ab, ab_std, ba, ba_std )
0839
0840
0841
0842 class SAB(object):
0843     """
0844     Comparison of pairs of SEvt
0845
0846     This is based on the seq array histories, requiring the SEvt to have been created with
0847     OPTICKS_EVENT_MODE such as::
0848
0849          HitPhotonSeq
0850          HitSeq
0851
0852     See ~/o/cxs_min.sh
0853     """
0854     def __getitem__(self, sli):
0855         self.sli = sli
0856         return self
0857
0858     def __init__(self, a, b):
0859
0860         self.sli = slice(None)
0861         log.info("[")
0862         if a.q is None or b.q is None:
0863             qcf = None
0864             qcf0 = None
0865         else:
0866             qcf = QCF( a.q, b.q, symbol="qcf")
0867             qcf0 = QCFZero(qcf) if "ZERO" in os.environ else None
0868         pass
0869         if a.meta is None or b.meta is None:
0870             meta = None
0871         else:
0872             meta = NPMeta.Compare( a.meta, b.meta  )
0873         pass
0874
0875         self.a = a
0876         self.b = b
0877         self.qcf = qcf
0878         self.qcf0 = qcf0
0879         self.meta = meta
0880         log.info("]")
0881
0882
0883     def q_and(self, a_prefix, b_prefix):
0884         """
0885         :param a_prefix: eg TO BT BT BT BT BT SA
0886         :param b_prefix: eg TO BT BT BT BT SA
0887         :return wab: indices where self.a.q startswith a_prefix and self.b.q startswith b_prefix
0888
0889         For example with input photons which start exactly the same this
0890         provides a way to look at history migration, eg an extra BT on one side
0891         compared to the other
0892         """
0893         return np.flatnonzero( np.logical_and( np.char.startswith(self.a.q, a_prefix.encode("utf-8")), np.char.startswith(self.b.q, b_prefix.encode("utf-8")) ))
0894
0895
0896     def q_startswith(self, prefix="TO BT BT SA"):
0897         """
0898         This is a cheating way to find accidentally random aligned photons
0899         for debugging purposes.
0900         It is especially useful with eg rain_point storch.h where all photons
0901         start with same position and direction.
0902         """
0903         return np.where( np.logical_and( self.a.q == self.b.q, np.char.startswith(self.a.q, prefix.encode("utf-8")) ))
0904
0905     def __repr__(self):
0906         ab = self
0907         a = self.a
0908         b = self.b
0909         qcf = self.qcf
0910         qcf0 = self.qcf0
0911         meta = self.meta
0912
0913         lines = []
0914         lines.append("SAB")
0915
0916         if not "BRIEF" in os.environ:
0917             lines.append(str(a))
0918             lines.append(repr(a.f))
0919             lines.append(str(b))
0920             lines.append(repr(b.f))
0921             if not qcf is None:
0922                 lines.append("ab.qcf.aqu")
0923                 lines.append(repr(qcf.aqu))
0924                 lines.append("ab.qcf.bqu")
0925                 lines.append(repr(qcf.bqu))
0926             pass
0927         pass
0928         lines.append("a.CHECK : %s " % a.CHECK)
0929         lines.append("b.CHECK : %s " % b.CHECK)
0930         if not ab.qcf is None:
0931             lines.append("ab.qcf[:40]")
0932             lines.append(repr(ab.qcf[:40]))
0933         pass
0934         if not ab.qcf0 is None:
0935             lines.append("ab.qcf0[:30])")
0936             lines.append(repr(ab.qcf0[:30]))
0937         pass
0938         if not meta is None:
0939             lines.append(str(meta))
0940         pass
0941         return "\n".join(lines)
0942
0943
0944
0945
0946 class KeyIdx(object):
0947     def __init__(self, kk, symbol="msab.ik"):
0948         self.kk = kk
0949         self.symbol = symbol
0950     def __getattr__(self, k):
0951         wk = np.where( self.kk == k)[0]
0952         return wk[0] if len(wk) == 1 else -1
0953     def __repr__(self):
0954         kk = self.kk
0955         symbol = self.symbol
0956         ii = np.arange(len(kk))
0957         lines = []
0958         lines.append("KeyIdx %s" % symbol)
0959         for i, k in enumerate(kk):
0960             lines.append("%s.%s = %d " % (symbol, k, i))
0961         pass
0962         return "\n".join(lines)
0963
0964 class MSAB(object):
0965     """
0966     Comparison of metadata across NEVT pairs of SEvt
0967     Typically NEVT is small, eg 10
0968     """
0969
0970     EXPRS = r"""
0971
0972     %(sym)s.c_itab[%(sym)s.ik.YSAVE,1].sum()/%(sym)s.c_itab[%(sym)s.ik.YSAVE,0].sum()  # N=1/N=0 FMT:%%10.4f
0973
0974     %(sym)s.c_itab[%(sym)s.ik.YSAVE,0].sum()/%(sym)s.c_itab[%(sym)s.ik.YSAVE,1].sum()  # N=0/N=1 FMT:%%10.4f
0975
0976     %(sym)s.c_ftab[0,1].sum()/%(sym)s.c_ftab[0,0].sum()  # ratio of duration totals N=1/N=0 FMT:%%10.4f
0977
0978     np.diff(%(sym)s.c_ttab)/1e6   # seconds between event starts
0979
0980     np.diff(%(sym)s.c_ttab)[0,1].sum()/np.diff(%(sym)s.c_ttab)[0,0].sum() # ratio N=1/N=0 FMT:%%10.4f
0981
0982     np.c_[%(sym)s.c_ttab[0].T].view('datetime64[us]') # SEvt start times (UTC)
0983
0984     %(sym)s.c2tab  # c2sum, c2n, c2per for each event
0985
0986     %(sym)s.c2tab[0,:].sum()/%(sym)s.c2tab[1,:].sum() # c2per_tot FMT:%%10.4f
0987
0988     """
0989
0990     def __init__(self, NEVT, AFOLD, BFOLD, symbol="%(sym)s"):
0991
0992         efmt = "%0.3d"
0993         assert efmt in AFOLD and efmt in BFOLD
0994         self.symbol = symbol
0995         self.exprs = list(filter(None,textwrap.dedent(self.EXPRS).split("\n")))
0996
0997         itabs = []
0998         ftabs = []
0999         ttabs = []
1000
1001         first = True
1002         kk0 = None
1003         skk0 = None
1004         ffield0 = None
1005         ifield0 = None
1006         tfield0 = None
1007         c2tab = np.zeros( (3, NEVT)  )
1008
1009         mab = {}
1010
1011         print("NEVT:%d" % NEVT)
1012         assert( NEVT > 0 )
1013
1014         for i in range(NEVT):
1015
1016             afold = AFOLD % i
1017             bfold = BFOLD % i
1018             a = SEvt.Load(afold,symbol="a", quiet=True)
1019             b = SEvt.Load(bfold,symbol="b", quiet=True)
1020             if a is None:
1021                 log.fatal("FAILED TO LOAD SEVT A FROM AFOLD: %s " % afold )
1022             pass
1023             if b is None:
1024                 log.fatal("FAILED TO LOAD SEVT B FROM BFOLD: %s " % bfold )
1025             pass
1026
1027             ab = SAB(a,b)
1028             mab[i] = ab
1029
1030
1031             kk = ab.meta.kk
1032             #print("ab.meta.kk (A,B common keys):%s" % str(kk))
1033
1034             skk = ab.meta.skk
1035             tab = ab.meta.tab
1036
1037             # metadata fields that look like floats, integers, timestamps
1038             ffield = np.unique(np.where(np.char.find(tab,'.') > -1 )[0])
1039             ifield = np.unique(np.where( np.logical_and( np.char.str_len( tab ) < 15, np.char.find(tab, '.') == -1 ))[0])
1040             tfield = np.unique(np.where( np.logical_and( np.char.str_len( tab ) > 15, np.char.find(tab, '.') == -1 ))[0])
1041
1042             if first:
1043                 first = False
1044                 skk0 = skk
1045                 kk0 = kk
1046                 tab0 = tab
1047                 ffield0 = ffield
1048                 ifield0 = ifield
1049                 tfield0 = tfield
1050             else:
1051                 assert np.all( skk == skk0 )
1052                 assert kk == kk0
1053                 assert tab.shape == tab0.shape
1054                 assert np.all( ffield == ffield0 )
1055                 assert np.all( ifield == ifield0 )
1056                 assert np.all( tfield == tfield0 )
1057             pass
1058
1059             itab = tab[ifield]
1060             ftab = tab[ffield]
1061             ttab = tab[tfield]
1062
1063             itabs.append(itab)
1064             ftabs.append(ftab)
1065             ttabs.append(ttab)
1066
1067             if not ab.qcf is None:
1068                 c2tab[0,i] = ab.qcf.c2sum
1069                 c2tab[1,i] = ab.qcf.c2n
1070                 c2tab[2,i] = ab.qcf.c2per
1071             pass
1072
1073             # hmm how to present together with c2sum, c2n, c2per ?
1074             if "c2desc" in os.environ:
1075                 fmt = " %s : %%s " % efmt
1076                 print( fmt % ( i, ab.qcf.c2desc))
1077             else:
1078                 if "DUMP" in os.environ:print(repr(ab))
1079             pass
1080         pass
1081
1082
1083         self.c2tab = c2tab
1084         self.c2tab_zero = np.all( c2tab == 0. )
1085
1086         self.mab = mab
1087         akk = np.array( list(kk0))
1088         ikk = akk[ifield]
1089         fkk = akk[ffield]
1090         tkk = akk[tfield]
1091
1092         ik = KeyIdx(ikk, symbol="%s.ik" % self.symbol )
1093         fk = KeyIdx(fkk, symbol="%s.fk" % self.symbol )
1094         tk = KeyIdx(tkk, symbol="%s.tk" % self.symbol )
1095
1096         c_itab = np.zeros( (itab.shape[0], itab.shape[1], NEVT ), dtype=np.uint64 )
1097         for i in range(len(itabs)):
1098             c_itab[:,:,i] = itabs[i]
1099         pass
1100         c_ftab = np.zeros( (ftab.shape[0], ftab.shape[1], NEVT ), dtype=np.float64 )
1101         for i in range(len(ftabs)):
1102             c_ftab[:,:,i] = ftabs[i]
1103         pass
1104         c_ttab = np.zeros( (ttab.shape[0], ttab.shape[1], NEVT ), dtype=np.uint64 )
1105         for i in range(len(ttabs)):
1106             c_ttab[:,:,i] = ttabs[i]
1107         pass
1108
1109         self.ikk = ikk
1110         self.fkk = fkk
1111         self.tkk = tkk
1112
1113         self.ik = ik
1114         self.fk = fk
1115         self.tk = tk
1116
1117         self.c_itab = c_itab
1118         self.c_ftab = c_ftab
1119         self.c_ttab = c_ttab
1120
1121
1122
1123     def annotab(self, tabsym, keys, nline=3):
1124         """
1125         :param tabsym: attribute symbol of table
1126         :param keys: array of names for each table item
1127         :param nline: number of lines for each table item
1128         :return str: annoted table string
1129
1130         Annotate a numpy array repr with keys
1131         TODO: split off into reusable AnnoTab? object
1132         """
1133         expr = "%%(sym)s.%s" % tabsym
1134         lines = []
1135         lines.append("\n%s\n" % ( expr % {'sym':self.symbol} ))
1136         rawlines = repr(eval(expr % {'sym':"self" })).split("\n")
1137         for i, line in enumerate(rawlines):
1138             j = i//nline
1139             anno = "%2d:%s" % (j,keys[j]) if i % nline == 0 else ""
1140             lines.append("%s       %s " % (line, anno ))
1141         pass
1142         return "\n".join(lines)
1143
1144     def __repr__(self):
1145         lines = []
1146         lines.append(self.annotab("c_itab", self.ikk, 3))
1147         lines.append(self.annotab("c_ftab", self.fkk, 3))
1148         pass
1149         fmt_ptn = re.compile("FMT:(.*)\\s*")
1150
1151         for expr in self.exprs:
1152             fmt_match = fmt_ptn.search(expr)
1153             fmt = fmt_match.groups()[0].replace("%%","%") if not fmt_match is None else None
1154             if "c2tab" in expr and self.c2tab_zero: continue
1155             lines.append("\n%s\n" % ( expr % {'sym':self.symbol} ))
1156             value = eval( expr % {'sym':"self"} )
1157             urep = repr(value) if fmt is None else fmt % value
1158             lines.append(urep)
1159         pass
1160         return "\n".join(lines)
1161
1162
1163
1164 if __name__ == '__main__':
1165     logging.basicConfig(level=logging.INFO)
1166
1167     FOLD = os.environ.get("FOLD", None)
1168     log.info(" -- SEvt.Load FOLD" )
1169     a = SEvt.Load(FOLD, symbol="a")   # optional photon histories
1170     print(a)
1171
1172     beg_ = "a.f.record[:,0,0,:3]"
1173     beg = eval(beg_)
1174
1175     end_ = "a.f.photon[:,0,:3]"
1176     end = eval(end_)
1177
1178
1179     #label0, ppos0 = None, None
1180     label0, ppos0 = "b:%s" % beg_ , beg
1181
1182     #label0, ppos0 = None, None
1183     label1, ppos1 = "r:%s" % end_ , end
1184
1185
1186     HEADLINE = "%s %s" % ( a.LAYOUT, a.CHECK )
1187     label = "\n".join( filter(None, [HEADLINE, label0, label1]))
1188     print(label)
1189
1190     if MODE == 0:
1191         print("not plotting as MODE 0  in environ")
1192     elif MODE == 2:
1193         fig, ax = mpplt_plotter(label=label)
1194
1195         ax.set_ylim(-250,250)
1196         ax.set_xlim(-500,500)
1197
1198         if not ppos0 is None: ax.scatter( ppos0[:,H], ppos0[:,V], s=1, c="b" )
1199         if not ppos1 is None: ax.scatter( ppos1[:,H], ppos1[:,V], s=1, c="r" )
1200
1201         fig.show()
1202     elif MODE == 3:
1203         pl = pvplt_plotter(label)
1204         os.environ["EYE"] = "0,100,165"
1205         os.environ["LOOK"] = "0,0,165"
1206         pvplt_viewpoint(pl)
1207         pl.add_points(ppos0)
1208         pl.show()
1209     pass
1210 pass
1211