eic-opticks/ana/rotation.py

0001 #!/usr/bin/env python
0002 """
0003 rotate.py
0004 ===========
0005
0006 Make sense of GDML physvol/rotation and global to local transforms
0007
0008 ::
0009
0010      71423       <physvol copynumber="11336" name="pLPMT_Hamamatsu_R128600x353fc90">
0011      71424         <volumeref ref="HamamatsuR12860lMaskVirtual0x3290b70"/>
0012      71425         <position name="pLPMT_Hamamatsu_R128600x353fc90_pos" unit="mm" x="-7148.9484" y="17311.741" z="-5184.2567"/>
0013      71426         <rotation name="pLPMT_Hamamatsu_R128600x353fc90_rot" unit="deg" x="-73.3288783033161" y="-21.5835981926051" z="-96.2863976680901"/>
0014      71427       </physvol>
0015
0016
0017 ::
0018
0019     epsilon:src blyth$ grep \"rotation\" *.cc
0020     G4GDMLReadDefine.cc:      if (tag=="rotation") { RotationRead(child); } else
0021     G4GDMLReadParamvol.cc:      if (tag=="rotation") { VectorRead(child,rotation); } else
0022     G4GDMLReadSolids.cc:      if (tag=="rotation") { VectorRead(child,rotation); } else
0023     G4GDMLReadSolids.cc:      if (tag=="rotation") { VectorRead(child,rotation); } else
0024     G4GDMLReadStructure.cc:     else if (tag=="rotation")
0025     G4GDMLReadStructure.cc:      if (tag=="rotation")
0026     epsilon:src blyth$ pwd
0027     /usr/local/opticks_externals/g4_1042.build/geant4.10.04.p02/source/persistency/gdml/src
0028
0029
0030 ::
0031
0032     g4-cls Rotation
0033
0034     336   HepRotation & rotateX(double delta);
0035     337   // Rotation around the x-axis; equivalent to R = RotationX(delta) * R
0036     338
0037     339   HepRotation & rotateY(double delta);
0038     340   // Rotation around the y-axis; equivalent to R = RotationY(delta) * R
0039     341
0040     342   HepRotation & rotateZ(double delta);
0041     343   // Rotation around the z-axis; equivalent to R = RotationZ(delta) * R
0042     344
0043
0044
0045 """
0046
0047 import os, sympy as sp, numpy as np
0048 from sympy import pprint as pp
0049
0050
0051 def tr_inverse(tt):
0052     """
0053     decompose translate rotate matrix by inspection,
0054     negate the translation and transpose the rotation
0055
0056         In [167]: tr
0057         Out[167]:
0058         array([[   -0.1018,    -0.9243,     0.3679,     0.    ],
0059                [    0.2466,    -0.3817,    -0.8908,     0.    ],
0060                [    0.9638,     0.    ,     0.2668,     0.    ],
0061                [   -0.003 ,     0.0127, 19433.9994,     1.    ]])
0062
0063         In [168]: np.dot( hit, tr )
0064         Out[168]: array([-112.6704,  165.9216,  109.6381,    1.    ])
0065
0066
0067     """
0068     it = np.eye(4)
0069     it[3,:3] = -tt[3,:3]
0070
0071     ir = np.eye(4)
0072     ir[:3,:3] = tt[:3,:3].T
0073
0074     tr = np.dot(it,ir)
0075     return tr
0076
0077
0078 class Instance(object):
0079     def __init__(self, ridx):
0080         """
0081         epsilon:1 blyth$ inp GMergedMesh/?/iidentity.npy GMergedMesh/?/itransforms.npy
0082         a :                                  GMergedMesh/1/iidentity.npy :        (25600, 5, 4) : a4a7deb934cae243b9181c80ddc1066b : 20200719-2129
0083         b :                                GMergedMesh/1/itransforms.npy :        (25600, 4, 4) : 29a7bf21dabfd4a6f9228fadb7edabca : 20200719-2129
0084         c :                                  GMergedMesh/2/iidentity.npy :        (12612, 6, 4) : 4423ba6434c39aff488e6784df468ae1 : 20200719-2129
0085         d :                                GMergedMesh/2/itransforms.npy :        (12612, 4, 4) : 766b1e274449b0d9f2ecc35d58b52a71 : 20200719-2129
0086         e :                                  GMergedMesh/3/iidentity.npy :         (5000, 6, 4) : 52c59e1bb3179c404722c2df4c26ac81 : 20200719-2129
0087         f :                                GMergedMesh/3/itransforms.npy :         (5000, 4, 4) : 1ff4e96acee67137c4740b05e6684c93 : 20200719-2129
0088         g :                                  GMergedMesh/4/iidentity.npy :         (2400, 6, 4) : 08846aa446e53c50c1a7cea89674a398 : 20200719-2129
0089         h :                                GMergedMesh/4/itransforms.npy :         (2400, 4, 4) : aafe0245a283080c130d8575b7a83e67 : 20200719-2129
0090
0091         #. iidentity is now reshaped shortly after creation to have same item count as itransforms
0092         """
0093         tt = np.load(os.path.expandvars("$GC/GMergedMesh/%d/itransforms.npy" % ridx))
0094         ii = np.load(os.path.expandvars("$GC/GMergedMesh/%d/iidentity.npy" % ridx))
0095
0096         assert tt.shape[1:] == (4,4)
0097         assert len(ii.shape) == 3 and ii.shape[2] == 4
0098         assert len(tt) == len(ii)
0099         nvol = ii.shape(1)   # physvol per instance
0100
0101         self.ii = ii
0102         self.tt = tt
0103         self.nvol = nvol
0104
0105     def find_instance_index(self, pmtid):
0106         """
0107         Using vol 0 corresponds to the outer volume of the instance :
0108         but the copyNo is duplicated for all the volumes of the instance
0109         so all 0:nvol are the same.
0110         """
0111         return np.where( self.ii[:,0,3] == pmtid )[0][0]
0112     def find_instance_transform(self, pmtid):
0113         ix = self.find_instance_index(pmtid)
0114         return self.tt[ix]
0115     def find_instance_transform_inverse(self, pmtid):
0116         tt = self.find_instance_transform(pmtid)
0117         return tr_inverse(tt)
0118
0119     def find_local_pos(self, pmtid, global_pos):
0120         assert global_pos.shape == (4,)
0121         tr = self.find_instance_transform_inverse(pmtid)
0122         local_pos = np.dot( global_pos, tr )
0123         return local_pos
0124
0125
0126
0127 def three_to_four(M3):
0128     assert M3.shape == (3,3)
0129     M4 = sp.zeros(4)
0130     for i in range(3):
0131         for j in range(3):
0132             M4[i*4+j] = M3.row(i)[j]
0133         pass
0134     M4[15] = 1
0135     return M4
0136 pass
0137
0138
0139 alpha,beta,gamma = sp.symbols("alpha beta gamma")
0140
0141 row0 = rxx,ryx,rzx,rwx = sp.symbols("rxx,ryx,rzx,rwx")
0142 row1 = rxy,ryy,rzy,rwy = sp.symbols("rxy,ryy,rzy,rwy")
0143 row2 = rxz,ryz,rzz,rwz = sp.symbols("rxz,ryz,rzz,rwz")
0144 row3 = tx,ty,tz,tw     = sp.symbols("tx,ty,tz,tw")
0145 RTxyz = sp.Matrix([row0,row1,row2,row3])
0146
0147 v_rid = [
0148    (rxx,1),(ryx,0),(rzx,0),
0149    (rxy,0),(ryy,1),(rzy,0),
0150    (rxz,0),(ryz,0),(rzz,1) ]    # identity rotation
0151
0152 v_rw = [(rwx,0),(rwy,0),(rwz,0)]
0153 v_t0 = [(tx,0),(ty,0),(tz,0),(tw,1)] # identity translation
0154 v_tw = [(tw,1),]
0155
0156 RT = RTxyz.subs(v_rw+v_tw)
0157 R = RTxyz.subs(v_rw+v_t0)
0158 T = RTxyz.subs(v_rid+v_rw+v_tw)
0159
0160
0161 x,y,z,w = sp.symbols("x,y,z,w")
0162 P = sp.Matrix([[x,y,z,w]])
0163
0164 assert P.shape == (1,4)
0165 P1 = P.subs([(w,1)])    # position
0166 P0 = P.subs([(w,0)])    # direction vector
0167
0168
0169
0170 deg = np.pi/180.
0171 v_rot = [(alpha,-73.3288783033161*deg),(beta,-21.5835981926051*deg),(gamma, -96.2863976680901*deg)]
0172 v_pos = [(tx, -7148.9484),(ty,17311.741), (tz,-5184.2567)]
0173
0174 lhit0 = np.array([-112.67072395684227,165.92175413608675,109.63878699927591,1])  # from debug session in ProcessHits
0175 v_lhit0 = [(x,lhit0[0]), (y,lhit0[1]), (z,lhit0[2])]
0176
0177 hit = np.array([-7250.504552589168,17122.963751776308,-5263.596996014085, 1])  # global hit position
0178 v_hit = [(x,hit[0]),(y,hit[1]),(z,hit[2]),(w,hit[3])]
0179
0180
0181
0182 pmtid = 11336     # "BP=junoSD_PMT_v2::ProcessHits tds" debug session, see jnu/opticks-junoenv-hit.rst
0183 it = Instance(3)  # ridx 3 is Hamamatsu PMTs
0184 tt = it.find_instance_transform(pmtid)
0185 tr = it.find_instance_transform_inverse(pmtid)
0186 lhit = it.find_local_pos(pmtid, hit )
0187
0188 assert np.allclose(lhit0, lhit), lhit0-lhit
0189
0190
0191 zhit = lhit[2]
0192 rhit = np.sqrt(lhit[0]*lhit[0]+lhit[1]*lhit[1])
0193 # 249 185  axes of cathode ellipsoid for Hamamatsu (see ana/gpmt.py)
0194 one = ((zhit*zhit)/(185.*185.))+((rhit*rhit)/(249.*249.))   ## hmm generic way to get major axes of cathode ellipsoid ?
0195
0196
0197
0198
0199 Rx = three_to_four(sp.rot_axis1(alpha))
0200 IRx = three_to_four(sp.rot_axis1(-alpha))
0201
0202 Ry = three_to_four(sp.rot_axis2(beta))
0203 IRy = three_to_four(sp.rot_axis2(-beta))
0204
0205 Rz = three_to_four(sp.rot_axis3(gamma))
0206 IRz = three_to_four(sp.rot_axis3(-gamma))
0207
0208 R0 = Rz*Ry*Rx
0209 IR0 = IRx*IRy*IRz    # NB : reversed order
0210 assert IR0.transpose() == R0       # the inverse of a rotation matrix is its transpose
0211 assert R0.transpose() == IR0
0212
0213 R1 = Rx*Ry*Rz
0214 IR1 = IRz*IRy*IRx
0215 assert IR1.transpose() == R1
0216 assert R1.transpose() == IR1
0217
0218
0219
0220
0221
0222
0223
0224 # row3 as translation is used for simpler matching with glm/OpenGL standard practice when serializing transforms
0225 # (col3 as translation is the other possibility)
0226 R1T = sp.Matrix([R1.row(0), R1.row(1), R1.row(2), T.row(3)])
0227 R1T_check = R1*T
0228
0229 assert R1T == R1T_check
0230 assert P*R1*T == P*R1T
0231 # so that is rotate then translate
0232 #    as rotations are around the origin that is appropriate for orienting
0233 #    a PMT before translating it into place
0234
0235 """
0236 IR1T
0237 Matrix([
0238 [cos(beta)*cos(gamma), sin(alpha)*sin(beta)*cos(gamma) - sin(gamma)*cos(alpha),  sin(alpha)*sin(gamma) + sin(beta)*cos(alpha)*cos(gamma), 0],
0239 [sin(gamma)*cos(beta), sin(alpha)*sin(beta)*sin(gamma) + cos(alpha)*cos(gamma), -sin(alpha)*cos(gamma) + sin(beta)*sin(gamma)*cos(alpha), 0],
0240 [          -sin(beta),                                    sin(alpha)*cos(beta),                                     cos(alpha)*cos(beta), 0],
0241 [                  tx,                                                      ty,                                                       tz, 1]])
0242
0243 """
0244
0245
0246
0247 IR1T = sp.Matrix([IR1.row(0), IR1.row(1), IR1.row(2), T.row(3)])  # matrix in terms of alpha,beta,gamma,tx,ty,tz
0248
0249 M = IR1T.subs(v_rot+v_pos)
0250
0251
0252
0253 match = np.allclose( np.array(M, dtype=np.float32 ), tt )
0254 assert match
0255
0256
0257 MI = np.array(M.inv(), dtype=np.float64)
0258
0259 amx = np.abs(MI-tr).max()    # not close enough for np.allclose
0260 assert amx < 1e-3
0261
0262
0263
0264
0265 T2 = sp.Matrix([[1,0,0,0],
0266                [0,1,0,0],
0267                [0,0,1,0],
0268                [tx,ty,tz,1]])
0269
0270 assert T == T2, (T, T2)
0271
0272
0273 IT = sp.Matrix([[1,0,0,0],
0274                [0,1,0,0],
0275                [0,0,1,0],
0276                [-tx,-ty,-tz,1]])
0277
0278
0279
0280
0281
0282
0283
0284
0285 """
0286      71423       <physvol copynumber="11336" name="pLPMT_Hamamatsu_R128600x353fc90">
0287      71424         <volumeref ref="HamamatsuR12860lMaskVirtual0x3290b70"/>
0288      71425         <position name="pLPMT_Hamamatsu_R128600x353fc90_pos" unit="mm" x="-7148.9484" y="17311.741" z="-5184.2567"/>
0289      71426         <rotation name="pLPMT_Hamamatsu_R128600x353fc90_rot" unit="deg" x="-73.3288783033161" y="-21.5835981926051" z="-96.2863976680901"/>
0290      71427       </physvol>
0291 """
0292
0293
0294 if 0:
0295     print("\nRx")
0296     pp(Rx)
0297     print("\nRy")
0298     pp(Ry)
0299     print("\nRz")
0300     pp(Rz)
0301
0302
0303 def rotateX():
0304     """
0305     This demonstrates that HepRotation::rotateX is multiplying to the rhs::
0306
0307         Rxyz*Rx
0308
0309     66 HepRotation & HepRotation::rotateX(double a) {    // g4-cls Rotation
0310     67   double c1 = std::cos(a);
0311     68   double s1 = std::sin(a);
0312     69   double x1 = ryx, y1 = ryy, z1 = ryz;
0313     70   ryx = c1*x1 - s1*rzx;
0314     71   ryy = c1*y1 - s1*rzy;
0315     72   ryz = c1*z1 - s1*rzz;
0316     73   rzx = s1*x1 + c1*rzx;
0317     74   rzy = s1*y1 + c1*rzy;
0318     75   rzz = s1*z1 + c1*rzz;
0319     76   return *this;
0320     77 }
0321     """
0322     pass
0323
0324     x1,y1,z1 = sp.symbols("x1,y1,z1")
0325     v_ry = [(ryx,x1),(ryy,y1),(ryz,z1)]
0326
0327     Xlhs = (Rx*R).subs(v_ry)
0328     Xrhs = (R*Rx).subs(v_ry)   # HepRotation::rotateX is multiplying on the rhs
0329
0330
0331     print(rotateX.__doc__)
0332
0333     print("\nR")
0334     pp(R)
0335
0336     print("\nv_ry")
0337     pp(v_ry)
0338
0339     print("\nXrhs = (R*Rx).subs(v_ry) ")
0340     pp(Xrhs)
0341
0342     #print("\nXlhs = (Rx*R).subs(v_ry) ")  clearly rotateX not doing this
0343     #pp(Xlhs)
0344
0345 def rotateY():
0346     """
0347     This demonstrates that HepRotation::rotateY is multiplying to the rhs::
0348
0349         Rxyz*Ry
0350
0351     079 HepRotation & HepRotation::rotateY(double a){
0352      80   double c1 = std::cos(a);
0353      81   double s1 = std::sin(a);
0354      82   double x1 = rzx, y1 = rzy, z1 = rzz;
0355      83   rzx = c1*x1 - s1*rxx;
0356      84   rzy = c1*y1 - s1*rxy;
0357      85   rzz = c1*z1 - s1*rxz;
0358      86   rxx = s1*x1 + c1*rxx;
0359      87   rxy = s1*y1 + c1*rxy;
0360      88   rxz = s1*z1 + c1*rxz;
0361      89   return *this;
0362      90 }
0363     """
0364     x1,y1,z1 = sp.symbols("x1,y1,z1")
0365     v_rz = [(rzx,x1),(rzy,y1),(rzz,z1)]
0366
0367     Ylhs = (Ry*R).subs(v_rz)
0368     Yrhs = (R*Ry).subs(v_rz)
0369
0370     print(rotateY.__doc__)
0371
0372     print("\nR")
0373     pp(R)
0374
0375     print("\nv_rz")
0376     pp(v_rz)
0377
0378     print("\nYrhs = (R*Ry).subs(v_rz) ")
0379     pp(Yrhs)
0380
0381
0382
0383 def rotateZ():
0384     """
0385     This demonstrates that HepRotation::rotateZ is multiplying to the rhs::
0386
0387         Rxyz*Rz
0388
0389
0390     092 HepRotation & HepRotation::rotateZ(double a) {
0391      93   double c1 = std::cos(a);
0392      94   double s1 = std::sin(a);
0393      95   double x1 = rxx, y1 = rxy, z1 = rxz;
0394      96   rxx = c1*x1 - s1*ryx;
0395      97   rxy = c1*y1 - s1*ryy;
0396      98   rxz = c1*z1 - s1*ryz;
0397      99   ryx = s1*x1 + c1*ryx;
0398     100   ryy = s1*y1 + c1*ryy;
0399     101   ryz = s1*z1 + c1*ryz;
0400     102   return *this;
0401     103 }
0402     """
0403
0404     x1,y1,z1 = sp.symbols("x1,y1,z1")
0405     v_rx = [(rxx,x1),(rxy,y1),(rxz,z1)]
0406
0407     Zlhs = (Rz*R).subs(v_rx)
0408     Zrhs = (R*Rz).subs(v_rx)
0409
0410     print(rotateZ.__doc__)
0411
0412     print("\nR")
0413     pp(R)
0414
0415     print("\nv_rx")
0416     pp(v_rx)
0417
0418     print("\nZrhs = (R*Rz).subs(v_rx) ")
0419     pp(Zrhs)
0420
0421
0422
0423
0424 def G4GDMLReadStructure():
0425     """
0426
0427     289 void G4GDMLReadStructure::
0428     290 PhysvolRead(const xercesc::DOMElement* const physvolElement,
0429     291             G4AssemblyVolume* pAssembly)
0430     292 {
0431     ...
0432     372    G4Transform3D transform(GetRotationMatrix(rotation).inverse(),position);
0433     373    transform = transform*G4Scale3D(scale.x(),scale.y(),scale.z());
0434
0435
0436     132 G4RotationMatrix
0437     133 G4GDMLReadDefine::GetRotationMatrix(const G4ThreeVector& angles)
0438     134 {
0439     135    G4RotationMatrix rot;
0440     136
0441     137    rot.rotateX(angles.x());
0442     138    rot.rotateY(angles.y());
0443     139    rot.rotateZ(angles.z());
0444     140    rot.rectify();  // Rectify matrix from possible roundoff errors
0445     141
0446     142    return rot;
0447     143 }
0448
0449
0450     g4-cls Transform3D (icc)
0451
0452     029 inline
0453      30 Transform3D::Transform3D(const CLHEP::HepRotation & m, const CLHEP::Hep3Vector & v) {
0454      31   xx_= m.xx(); xy_= m.xy(); xz_= m.xz();
0455      32   yx_= m.yx(); yy_= m.yy(); yz_= m.yz();
0456      33   zx_= m.zx(); zy_= m.zy(); zz_= m.zz();
0457      34   dx_= v.x();  dy_= v.y();  dz_= v.z();
0458      35 }
0459
0460
0461     NB the order  (rotateX, rotateY, rotateZ).inverse()
0462
0463
0464     In [18]: t[3218]        # use the instance index to give the instance transform : rot matrix and tlate look familiar
0465     Out[18]:
0466     array([[   -0.1018,     0.2466,     0.9638,     0.    ],
0467            [   -0.9243,    -0.3817,     0.    ,     0.    ],
0468            [    0.3679,    -0.8908,     0.2668,     0.    ],
0469            [-7148.948 , 17311.74  , -5184.257 ,     1.    ]], dtype=float32)
0470
0471
0472     In [70]: pp((Rx*Ry*Rz).transpose().subs(v_rot)*T.subs(v_pos))
0473      -0.101820513179743   0.24656591428434    0.963762332221457    0
0474      -0.924290430171623  -0.381689927419044  8.32667268468867e-17  0
0475      0.367858374634817   -0.890796300632177   0.266762379264878    0
0476          -7148.9484          17311.741            -5184.2567       1
0477
0478
0479     In [52]: (P*IT*R1).subs(v_hit+v_pos+v_rot)    ### << MATCHES <<
0480     Out[52]: Matrix([[-112.670723956843, 165.921754136086, 109.638786999275, 1.0]])
0481
0482     In [100]: (P*IT*Rx*Ry*Rz).subs(v_hit+v_pos+v_rot)
0483     Out[100]: Matrix([[-112.670723956843, 165.921754136086, 109.638786999275, 1]])
0484
0485
0486
0487     In [101]: (P*Rx*Ry*Rz*IT).subs(v_hit+v_pos+v_rot)
0488     Out[101]: Matrix([[7036.28119031864, -17145.8318197405, -14140.1045440872, 1]])
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498     (gdb) p local_pos
0499     $7 = {dx = -112.67072395684227, dy = 165.92175413608675, dz = 109.63878699927591, static tolerance = 2.22045e-14}
0500     (gdb) p trans
0501
0502
0503
0504     ## 4th row matches the GDML
0505
0506      71423       <physvol copynumber="11336" name="pLPMT_Hamamatsu_R128600x353fc90">
0507      71424         <volumeref ref="HamamatsuR12860lMaskVirtual0x3290b70"/>
0508      71425         <position name="pLPMT_Hamamatsu_R128600x353fc90_pos" unit="mm" x="-7148.9484" y="17311.741" z="-5184.2567"/>
0509      71426         <rotation name="pLPMT_Hamamatsu_R128600x353fc90_rot" unit="deg" x="-73.3288783033161" y="-21.5835981926051" z="-96.2863976680901"/>
0510      71427       </physvol>
0511
0512
0513              { rxx = -0.10182051317974285,   rxy = -0.92429043017162327,   rxz =  0.36785837463481702,
0514                ryx =  0.24656591428433955,   ryy = -0.38168992741904467,   ryz = -0.89079630063217707,
0515                rzx =  0.96376233222145669,   rzy =  0,                     rzz =  0.26676237926487772,
0516                tx =  -0.0035142754759363015, ty =   0.012573876562782971,  tz  =  19434.000031086449}
0517
0518              THIS IS THE INVERSE TRANSFORM
0519
0520
0521     From examples/UseGeant4/UseGeant.cc
0522
0523     dbg_affine_trans   Transformation:
0524     rx/x,y,z: -0.101821 -0.92429 0.367858
0525     ry/x,y,z: 0.246566 -0.38169 -0.890796
0526     rz/x,y,z: 0.963762 0 0.266762
0527     tr/x,y,z: -0.00351428 0.0125739 19434
0528
0529
0530     """
0531     R0 = Rx*Ry*Rz
0532     R0T = (Rx*Ry*Rz).T
0533
0534     R1 = Rz*Ry*Rx
0535     R1T = (Rz*Ry*Rx).T
0536
0537
0538     print(G4GDMLReadStructure.__doc__)
0539
0540     if 0:
0541         print("\nR0 = Rx*Ry*Rz \n")
0542         pp(R0)
0543         pp(R)
0544         pp(R0.subs(va))
0545     pass
0546
0547     print("\nR0T = (Rx*Ry*Rz).T   THIS MATCHES THE ROTATION PART OF THE ITRANSFORM \n")
0548     pp(R0T)
0549     pp(R)
0550     pp(R0T.subs(v_rot))
0551
0552
0553     if 0:
0554         print("\nR1 = Rz*Ry*Rx\n")
0555         pp(R1)
0556         pp(R)
0557         pp(R1.subs(va))
0558     pass
0559
0560
0561
0562 def translate():
0563     """
0564     Using col3 for the translation as opposed to glm/OpenGL approach
0565     of row3 for the translation is a transpose of the translation matrix,
0566     which means need to transpose the point for shape consistency
0567     and multiply in other order.
0568     """
0569
0570     print("P1")
0571     pp(P1)
0572     print("T")
0573     pp(T)
0574     P1_T = P1*T
0575     print("P1*T")
0576     pp(P1_T)
0577     print("P1*T*T")
0578     pp(P1*T*T)
0579     P1_T_reverse_transpose_check = (T.T*P1.T).T
0580     print("(T.T*P1.T).T")
0581     pp(P1_T_reverse_transpose_check)
0582     assert P1_T_reverse_transpose_check == P1_T
0583
0584
0585 def translate_rotate():
0586     """
0587     P1
0588      x  y  z  1
0589     TR
0590      rxx  ryx  rzx  0
0591
0592      rxy  ryy  rzy  0
0593
0594      rxz  ryz  rzz  0
0595
0596      tx   ty   tz   tw
0597     P*TR
0598      rxx x + rxy y + rxz z + tx w  ryx x + ryy y + ryz z + ty w  rzx x + rzy y + rzz z + tz w  tw w
0599     P*TR.subs(v_rid)
0600      tx w + x  ty w + y  tz w + z  tw w
0601
0602     """
0603
0604     print("R")
0605     pp(R)
0606     print("T")
0607     pp(T)
0608     print("T*R : row3 has translation and rotation mixed up : ie translation first and then rotation which depends  ")
0609     pp(T*R)
0610     print("R*T : familiar row3 as translation : that means rotate then translate ")
0611     pp(R*T)
0612
0613     print("RT")
0614     pp(RT)
0615     assert RT == R*T
0616
0617
0618     print("P1")
0619     pp(P1)
0620     print("P*RT : notice that the translation just gets added to rotated coordinates : ie rotation first and then translation")
0621     pp(P*RT)
0622
0623     P_RT = P*RT
0624     print("P*RT.subs(v_rid) : setting rotation to identity ")
0625     pp(P_RT.subs(v_rid))
0626
0627
0628
0629
0630 if __name__ == '__main__':
0631
0632     #rotateX()
0633     #rotateY()
0634     #rotateZ()
0635
0636     #translate()
0637     #translate_rotate()
0638
0639     G4GDMLReadStructure()
0640
0641
0642