CSG/tests/CSGFoundry_getFrame_Test.py

0001 #!/usr/bin/env python
0002 """
0003
0004 Majority have position differences at 4e-4 level::
0005
0006     In [23]: np.unique(np.where( d_tpos > 4e-4 )[0]).shape
0007     Out[23]: (4447,)
0008
0009     In [24]: np.unique(np.where( d_tpos > 5e-4 )[0]).shape
0010     Out[24]: (0,)
0011
0012 """
0013
0014 import os, numpy as np
0015 from opticks.ana.fold import Fold
0016
0017
0018 def py_transform_0(ip, m2w):
0019
0020     ipt = np.zeros( (len(ip),4,4) )
0021     ipt[:,0,:3] = ip[:,0,:3]
0022     ipt[:,0,3]  = 1.      # pos transformed as position
0023
0024     ipt[:,1,:3] = ip[:,1,:3]
0025     ipt[:,1,3]  = 0.      # mom tranformed as direction
0026
0027     ipt[:,2,:3] = ip[:,2,:3]
0028     ipt[:,2,3]  = 0.      # pol transformed as direction
0029
0030
0031     ipt[:,0] = np.dot( ipt[:,0], m2w )
0032     ipt[:,1] = np.dot( ipt[:,1], m2w )
0033     ipt[:,2] = np.dot( ipt[:,2], m2w )
0034
0035     return ipt
0036
0037 def py_transform_1(ip, m2w):
0038
0039     pos = np.zeros( (len(ip), 4) )
0040     mom = np.zeros( (len(ip), 4) )
0041     pol = np.zeros( (len(ip), 4) )
0042
0043     pos[:,:3] = ip[:,0,:3]
0044     pos[:,3] = 1.
0045
0046     mom[:,:3] = ip[:,1,:3]
0047     mom[:,3] = 0.
0048
0049     pol[:,:3] = ip[:,2,:3]
0050     pol[:,3] = 0.
0051
0052     tpos = np.dot( pos, m2w )
0053     tmom = np.dot( mom, m2w )
0054     tpol = np.dot( pol, m2w )
0055
0056     ipt = np.zeros( (len(ip),4,4) )
0057     ipt[:,0] = tpos
0058     ipt[:,1] = tmom
0059     ipt[:,2] = tpol
0060     return ipt
0061
0062
0063
0064 if __name__ == '__main__':
0065     t = Fold.Load()
0066     print("t.base:%s " % t.base)
0067     print(t)
0068     print(repr(t))
0069
0070     print("t.sframe")
0071     print(repr(t.sframe))
0072
0073
0074     m2w = t.m2w[0]
0075     ip = t.ip          # local frame input photons
0076
0077     pos = np.zeros( (len(ip), 4) )
0078     mom = np.zeros( (len(ip), 4) )
0079     pol = np.zeros( (len(ip), 4) )
0080
0081     pos[:,:3] = ip[:,0,:3] ; pos[:,3] = 1.    # position
0082     mom[:,:3] = ip[:,1,:3] ; mom[:,3] = 0.    # direction
0083     pol[:,:3] = ip[:,2,:3] ; pol[:,3] = 0.    # direction
0084
0085     tpos = np.dot( pos, m2w )
0086     tmom = np.dot( mom, m2w )
0087     tpol = np.dot( pol, m2w )
0088
0089     tpho = np.zeros( (len(ip),4,4) )
0090     tpho[:,0] = tpos
0091     tpho[:,1] = tmom
0092     tpho[:,2] = tpol
0093
0094
0095
0096
0097     ipt = {}
0098     ipt[0] = t.ipt0
0099     ipt[1] = t.ipt1
0100     ipt[2] = t.ipt2
0101     ipt[3] = tpho
0102
0103     lpt = {}
0104     lpt[0] = "ipt0 : C++ sframe transformed input photons with normalize:false "
0105     lpt[1] = "ipt1 : C++ sframe transformed input photons with normalize:true "
0106     lpt[2] = "ipt2 : C++ SEvt::getInputPhotons transformed input photons with normalize:true "
0107     lpt[3] = "ipt3 : python transformed input photons "
0108
0109
0110     for i in range(4):
0111         print(lpt[i])
0112     pass
0113
0114     print("cross comparison (q:3,i:4,j:4) cross comparison  ")
0115     print("HUH: ipt2 is standing out as worst position match " )
0116     print("EXPECTED ipt2 and ipt1 to be perfect match : but they are not ??? ")
0117
0118     qij = np.zeros( (3,4,4) )
0119     for q in range(3):
0120         for i in range(4):
0121             for j in range(4):
0122                 _qij = np.abs( ipt[i][:,q,:3] - ipt[j][:,q,:3] )
0123                 qij[q,i,j] = _qij.max()
0124             pass
0125         pass
0126     pass
0127     print("q:pos,mom,pol ")
0128     print("qij*1e6:\n",qij*1e6)
0129
0130
0131     normdiff_ = lambda v3:np.abs( 1. - np.sqrt(np.sum( np.power(v3, 2), axis=1)) ).max()
0132
0133     nd = np.zeros( (3,4) )
0134     for q in np.arange(1,3):
0135         for i in range(4):
0136             nd[q,i] = normdiff_( ipt[i][:,q,:3] )
0137         pass
0138     pass
0139     print("\nnormdiff_ difference from 1. (pos dummy zero)  ")
0140     print("nd*1e6:\n", nd*1e6)
0141
0142