g4cx/tests/recplot.py

0001 #!/usr/bin/env python
0002 """
0003 recplot.py
0004 ============
0005
0006 Comparing photon record points between A and B in two different histories,
0007 under the assumption that the histories should really be the same.
0008 Also geometry intersect position comparison between A and B as a check of the
0009 translation.
0010
0011
0012 """
0013 import os, numpy as np
0014 from opticks.ana.fold import Fold
0015 from opticks.ana.p import *       # including cf loaded from CFBASE
0016 from opticks.ana.eget import efloatlist_, elookce_, elook_epsilon_, eint_
0017 from opticks.ana.pvplt import *
0018
0019 GEOM = os.environ.get("GEOM", None)
0020 GDMLPath = os.environ.get("%s_GDMLPath" % GEOM, None)
0021 GDMLSub  = os.environ.get("%s_GDMLSub"  % GEOM, None)
0022
0023
0024 BOFF = efloatlist_("BOFF", "0,0,0")
0025
0026
0027 def make_local(gpos_, w2m):
0028     gpos = np.zeros( (len(gpos_), 4 ), dtype=np.float32 )
0029     gpos[:,:3] = gpos_
0030     gpos[:,3] = 1.
0031     lpos_ = np.dot( gpos, w2m )
0032     lpos = lpos_[:,:3]
0033     return lpos
0034
0035
0036 if __name__ == '__main__':
0037     a = Fold.Load("$A_FOLD", symbol="a")
0038     b = Fold.Load("$B_FOLD", symbol="b")
0039
0040     if GEOM == "J000":
0041         if GDMLSub == None:
0042             qa = "TO BT BT BT SD"
0043             qb = "TO BT BT BT BT BT SD"
0044         else:
0045             qa = "TO BT BT BT BT SD"       # GDMLSub wrap does not skip the hatbox
0046             qb = "TO BT BT BT BT BT SD"
0047         pass
0048     elif GEOM == "hama_body_log":
0049         qa = "TO BT SD"
0050         qb = "TO BT SD"
0051     else:
0052         assert 0, GEOM
0053     pass
0054
0055     na = len(qa.split())
0056     nb = len(qb.split())
0057     wa = np.where(a.seq[:,0] == cseqhis_(qa))[0]
0058     wb = np.where(b.seq[:,0] == cseqhis_(qb))[0]
0059     wc = np.intersect1d(wa, wb)  # indices in common
0060
0061     print("qa : %20s : len(wa) %d " % (qa,len(wa)))
0062     print("qb : %20s : len(wb) %d " % (qb,len(wb)))
0063     print(" len(wc) %d " % len(wc) )
0064
0065     # pt: point index
0066     apos_ = lambda pt:a.record[wc,pt,0,:3]
0067     bpos_ = lambda pt:b.record[wc,pt,0,:3]
0068
0069
0070     assert np.all( a.sframe.m2w == b.sframe.m2w )
0071     assert np.all( a.sframe.w2m == b.sframe.w2m )
0072
0073     w2m = a.sframe.w2m
0074
0075
0076     pl = pvplt_plotter()
0077     pvplt_viewpoint(pl)
0078
0079     for pt in range(1,na):
0080         gpos_ = apos_(pt)
0081         lpos = make_local(gpos_, w2m )
0082         pl.add_points( lpos, color="red" )
0083     pass
0084     for pt in range(1,nb):
0085         gpos_ = bpos_(pt)
0086         lpos = make_local(gpos_, w2m )
0087         lpos += BOFF
0088         pl.add_points( lpos, color="blue" )
0089     pass
0090
0091     pl.show()
0092
0093