Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:20

0001 #!/usr/bin/env python
0002 
0003 from __future__ import print_function
0004 from builtins import range
0005 from Alignment.MuonAlignment.geometryXMLparser import MuonGeometry, dtorder, cscorder
0006 import sys, getopt
0007 
0008 usage = "Usage: geometryDiff.py [-h|--help] [-e|--epsilon epsilon] geometry1.xml geometry2.xml"
0009 
0010 try:
0011   opts, args = getopt.getopt(sys.argv[1:], "he:", ["help", "epsilon="])
0012 except getopt.GetoptError as msg:
0013   print(usage, file=sys.stderr)
0014   sys.exit(2)
0015 
0016 if len(args) != 2:
0017   print(usage, file=sys.stderr)
0018   sys.exit(2)
0019 
0020 opts = dict(opts)
0021 
0022 if "-h" in opts or "--help" in opts:
0023   print(usage)
0024   sys.exit(0)
0025 
0026 epsilon = 1e-6
0027 if "-e" in opts: epsilon = float(opts["-e"])
0028 if "--epsilon" in opts: epsilon = float(opts["--epsilon"])
0029 
0030 geom1 = MuonGeometry(file(args[0]))
0031 geom2 = MuonGeometry(file(args[1]))
0032 
0033 from math import sin, cos, sqrt
0034 sqrtepsilon = sqrt(epsilon)
0035 
0036 def matrixmult(a, b):
0037   return [[sum([i*j for i, j in zip(row, col)]) for col in zip(*b)] for row in a]
0038 
0039 def transpose(a):
0040   return [[a[j][i] for j in range(len(a[i]))] for i in range(len(a))]
0041 
0042 def rotFromPhi(g):
0043   phix, phiy, phiz = g.phix, g.phiy, g.phiz
0044   rotX = [[1.,         0.,         0.,       ],
0045           [0.,         cos(phix),  sin(phix),],
0046           [0.,        -sin(phix),  cos(phix),]]
0047   rotY = [[cos(phiy),  0.,        -sin(phiy),],
0048           [0.,         1.,         0.,       ],
0049           [sin(phiy),  0.,         cos(phiy),]]
0050   rotZ = [[cos(phiz),  sin(phiz),  0.,       ],
0051           [-sin(phiz), cos(phiz),  0.,       ],
0052           [0.,         0.,         1.,       ]]
0053   return matrixmult(rotX, matrixmult(rotY, rotZ))
0054 
0055 def rotFromEuler(g):
0056   s1, s2, s3 = sin(g.alpha), sin(g.beta), sin(g.gamma)
0057   c1, c2, c3 = cos(g.alpha), cos(g.beta), cos(g.gamma)
0058   return [[c2 * c3,    c1 * s3 + s1 * s2 * c3,   s1 * s3 - c1 * s2 * c3,],
0059           [-c2 * s3,   c1 * c3 - s1 * s2 * s3,   s1 * c3 + c1 * s2 * s3,],
0060           [s2,        -s1 * c2,                  c1 * c2,               ]]
0061 
0062 def loopover(which):
0063   if which == "DT":
0064     keys = geom1.dt.keys()
0065     keys.sort(dtorder)
0066 
0067   elif which == "CSC":
0068     keys = geom1.csc.keys()
0069     keys.sort(cscorder)
0070 
0071   else: raise Exception
0072 
0073   for key in keys:
0074     if which == "DT":
0075       g1 = geom1.dt[key]
0076       g2 = geom2.dt[key]
0077     else:
0078       g1 = geom1.csc[key]
0079       g2 = geom2.csc[key]
0080 
0081     if g1.relativeto != g2.relativeto:
0082       print("%s %s relativeto=\"%s\" versus relativeto=\"%s\"" % (which, str(key), g1.relativeto, g2.relativeto))
0083 
0084     if abs(g1.x - g2.x) > epsilon or abs(g1.y - g2.y) > epsilon or abs(g1.z - g2.z) > epsilon:
0085       print("%s %s position difference: (%g, %g, %g) - (%g, %g, %g) = (%g, %g, %g)" % \
0086             (which, str(key), g1.x, g1.y, g1.z, g2.x, g2.y, g2.z, g1.x - g2.x, g1.y - g2.y, g1.z - g2.z))
0087 
0088     if "phix" in g1.__dict__:
0089       g1type = "phi"
0090       g1a, g1b, g1c = g1.phix, g1.phiy, g1.phiz
0091       g1rot = rotFromPhi(g1)
0092     else:
0093       g1type = "euler"
0094       g1a, g1b, g1c = g1.alpha, g1.beta, g1.gamma
0095       g1rot = rotFromEuler(g1)
0096 
0097     if "phix" in g2.__dict__:
0098       g2type = "phi"
0099       g2a, g2b, g2c = g2.phix, g2.phiy, g2.phiz
0100       g2rot = rotFromPhi(g2)
0101     else:
0102       g2type = "euler"
0103       g2a, g2b, g2c = g2.alpha, g2.beta, g2.gamma
0104       g2rot = rotFromEuler(g2)
0105     
0106     diff = matrixmult(g1rot, transpose(g2rot))
0107     if abs(diff[0][0] - 1.) > sqrtepsilon or abs(diff[1][1] - 1.) > sqrtepsilon or abs(diff[2][2] - 1.) > sqrtepsilon or \
0108        abs(diff[0][1]) > epsilon or abs(diff[0][2]) > epsilon or abs(diff[1][2]) > epsilon:
0109       print("%s %s rotation difference: %s(%g, %g, %g) - %s(%g, %g, %g) = %s" % \
0110             (which, str(key), g1type, g1a, g1b, g1c, g2type, g2a, g2b, g2c, str(diff)))
0111 
0112 loopover("DT")
0113 loopover("CSC")
0114 
0115