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