MuonAlignment/python/makeMuonMisalignmentScenario.py

0001 from __future__ import print_function
0002 from builtins import range
0003 from optparse import OptionParser
0004 from random import gauss
0005 from math import sqrt
0006 import os
0007 execfile("Alignment/MuonAlignment/data/idealTransformation.py")
0008
0009 ### Get variances and covariances from the commandline
0010
0011 parser = OptionParser(usage="Usage: python %prog outputName [options] (default is unit matrix times 1e-15)")
0012
0013 parser.add_option("--xx", dest="xx", help="variance of x (cm*cm)", default="1e-15")
0014 parser.add_option("--xy", dest="xy", help="covariance of x and y (cm*cm)", default="0")
0015 parser.add_option("--xz", dest="xz", help="covariance of x and z (cm*cm)", default="0")
0016 parser.add_option("--xphix", dest="xphix", help="covariance of x and phix (cm*rad)", default="0")
0017 parser.add_option("--xphiy", dest="xphiy", help="covariance of x and phiy (cm*rad)", default="0")
0018 parser.add_option("--xphiz", dest="xphiz", help="covariance of x and phiz (cm*rad)", default="0")
0019
0020 parser.add_option("--yy", dest="yy", help="variance of y (cm*cm)", default="1e-15")
0021 parser.add_option("--yz", dest="yz", help="covariance of y and z (cm*cm)", default="0")
0022 parser.add_option("--yphix", dest="yphix", help="covariance of y and phix (cm*rad)", default="0")
0023 parser.add_option("--yphiy", dest="yphiy", help="covariance of y and phiy (cm*rad)", default="0")
0024 parser.add_option("--yphiz", dest="yphiz", help="covariance of y and phiz (cm*rad)", default="0")
0025
0026 parser.add_option("--zz", dest="zz", help="variance of z (cm*cm)", default="1e-15")
0027 parser.add_option("--zphix", dest="zphix", help="covariance of z and phix (cm*rad)", default="0")
0028 parser.add_option("--zphiy", dest="zphiy", help="covariance of z and phiy (cm*rad)", default="0")
0029 parser.add_option("--zphiz", dest="zphiz", help="covariance of z and phiz (cm*rad)", default="0")
0030
0031 parser.add_option("--phixphix", dest="phixphix", help="variance of phix (rad*rad)", default="1e-15")
0032 parser.add_option("--phixphiy", dest="phixphiy", help="covariance of phix and phiy (rad*rad)", default="0")
0033 parser.add_option("--phixphiz", dest="phixphiz", help="covariance of phix and phiz (rad*rad)", default="0")
0034
0035 parser.add_option("--phiyphiy", dest="phiyphiy", help="variance of phiy (rad*rad)", default="1e-15")
0036 parser.add_option("--phiyphiz", dest="phiyphiz", help="covariance of phiy and phiz (rad*rad)", default="0")
0037
0038 parser.add_option("--phizphiz", dest="phizphiz", help="variance of phiz (rad*rad)", default="1e-15")
0039
0040 parser.add_option("-f", dest="force", help="force overwrite of output files", action="store_true")
0041
0042 options, args = parser.parse_args()
0043 if args is None or len(args) != 1:
0044     parser.print_help()
0045     exit(-1)
0046 outputName = args[0]
0047
0048 if not options.force:
0049     if os.path.exists(outputName + ".xml"):
0050         raise Exception(outputName + ".xml exists!")
0051     if os.path.exists(outputName + "_convert_cfg.py"):
0052         raise Exception(outputName + "_convert_cfg.py exists!")
0053     if os.path.exists(outputName + ".db"):
0054         raise Exception(outputName + ".db exists!")
0055     if os.path.exists(outputName + "_correlations.txt"):
0056         raise Exception(outputName + "_correlations.txt exists!")
0057
0058 components = "xx", "xy", "xz", "xphix", "xphiy", "xphiz", "yy", "yz", "yphix", "yphiy", "yphiz", "zz", "zphix", "zphiy", "zphiz", "phixphix", "phixphiy", "phixphiz", "phiyphiy", "phiyphiz", "phizphiz"
0059 for component in components:
0060     exec("%s = float(options.%s)" % (component, component))
0061
0062 ### Print out user's choices as diagnostics
0063
0064 print("Spread in each parameter: x %g mm" % (sqrt(xx)*10.))
0065 print("                          y %g mm" % (sqrt(yy)*10.))
0066 print("                          z %g mm" % (sqrt(zz)*10.))
0067 print("                          phix %g mrad" % (sqrt(phixphix)*1000.))
0068 print("                          phiy %g mrad" % (sqrt(phiyphiy)*1000.))
0069 print("                          phiz %g mrad" % (sqrt(phizphiz)*1000.))
0070 print()
0071
0072 print("Covariance matrix (x, y, z, phix, phiy, phiz):")
0073 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xx, xy, xz, xphix, xphiy, xphiz))
0074 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xy, yy, yz, yphix, yphiy, yphiz))
0075 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xz, yz, zz, zphix, zphiy, zphiz))
0076 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphix, yphix, zphix, phixphix, phixphiy, phixphiz))
0077 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphiy, yphiy, zphiy, phixphiy, phiyphiy, phiyphiz))
0078 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphiz, yphiz, zphiz, phixphiz, phiyphiz, phizphiz))
0079 print()
0080
0081 print("Correlation (x, y, z, phix, phiy, phiz):")
0082 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xx/sqrt(xx)/sqrt(xx), xy/sqrt(xx)/sqrt(yy), xz/sqrt(xx)/sqrt(zz), xphix/sqrt(xx)/sqrt(phixphix), xphiy/sqrt(xx)/sqrt(phiyphiy), xphiz/sqrt(xx)/sqrt(phizphiz)))
0083 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xy/sqrt(yy)/sqrt(xx), yy/sqrt(yy)/sqrt(yy), yz/sqrt(yy)/sqrt(zz), yphix/sqrt(yy)/sqrt(phixphix), yphiy/sqrt(yy)/sqrt(phiyphiy), yphiz/sqrt(yy)/sqrt(phizphiz)))
0084 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xz/sqrt(zz)/sqrt(xx), yz/sqrt(zz)/sqrt(yy), zz/sqrt(zz)/sqrt(zz), zphix/sqrt(zz)/sqrt(phixphix), zphiy/sqrt(zz)/sqrt(phiyphiy), zphiz/sqrt(zz)/sqrt(phizphiz)))
0085 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphix/sqrt(phixphix)/sqrt(xx), yphix/sqrt(phixphix)/sqrt(yy), zphix/sqrt(phixphix)/sqrt(zz), phixphix/sqrt(phixphix)/sqrt(phixphix), phixphiy/sqrt(phixphix)/sqrt(phiyphiy), phixphiz/sqrt(phixphix)/sqrt(phizphiz)))
0086 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphiy/sqrt(phiyphiy)/sqrt(xx), yphiy/sqrt(phiyphiy)/sqrt(yy), zphiy/sqrt(phiyphiy)/sqrt(zz), phixphiy/sqrt(phiyphiy)/sqrt(phixphix), phiyphiy/sqrt(phiyphiy)/sqrt(phiyphiy), phiyphiz/sqrt(phiyphiy)/sqrt(phizphiz)))
0087 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphiz/sqrt(phizphiz)/sqrt(xx), yphiz/sqrt(phizphiz)/sqrt(yy), zphiz/sqrt(phizphiz)/sqrt(zz), phixphiz/sqrt(phizphiz)/sqrt(phixphix), phiyphiz/sqrt(phizphiz)/sqrt(phiyphiy), phizphiz/sqrt(phizphiz)/sqrt(phizphiz)))
0088 print()
0089
0090 for correlation_coefficient in [abs(xy/sqrt(xx)/sqrt(yy)), abs(xz/sqrt(xx)/sqrt(zz)), abs(xphix/sqrt(xx)/sqrt(phixphix)), abs(xphiy/sqrt(xx)/sqrt(phiyphiy)), abs(xphiz/sqrt(xx)/sqrt(phizphiz)), \
0091                                 abs(yz/sqrt(yy)/sqrt(zz)), abs(yphix/sqrt(yy)/sqrt(phixphix)), abs(yphiy/sqrt(yy)/sqrt(phiyphiy)), abs(yphiz/sqrt(yy)/sqrt(phizphiz)),
0092                                 abs(zphix/sqrt(zz)/sqrt(phixphix)), abs(zphiy/sqrt(zz)/sqrt(phiyphiy)), abs(zphiz/sqrt(zz)/sqrt(phizphiz)),
0093                                 abs(phixphiy/sqrt(phixphix)/sqrt(phiyphiy)), abs(phixphiz/sqrt(phixphix)/sqrt(phizphiz)),
0094                                 abs(phiyphiz/sqrt(phiyphiy)/sqrt(phizphiz))]:
0095     if correlation_coefficient > 1.:
0096         raise Exception("Correlations must not be larger than one!")
0097
0098 ### Some useful mathematical transformations (why don't we have access to numpy?)
0099
0100 def mmult(a, b):
0101     """Matrix multiplication: mmult([[11, 12], [21, 22]], [[-1, 0], [0, 1]]) returns [[-11, 12], [-21, 22]]"""
0102     return [[sum([i*j for i, j in zip(row, col)]) for col in zip(*b)] for row in a]
0103
0104 def mvdot(m, v):
0105     """Applies matrix m to vector v: mvdot([[-1, 0], [0, 1]], [12, 55]) returns [-12, 55]"""
0106     return [i[0] for i in mmult(m, [[vi] for vi in v])]
0107
0108 def mtrans(a):
0109     """Matrix transposition: mtrans([[11, 12], [21, 22]]) returns [[11, 21], [12, 22]]"""
0110     return [[a[j][i] for j in range(len(a[i]))] for i in range(len(a))]
0111
0112 def cholesky(A):
0113     """Cholesky decomposition of the correlation matrix to properly normalize the transformed random deviates"""
0114
0115     # A = L * D * L^T = (L * D^0.5) * (L * D^0.5)^T where we want (L * D^0.5), the "square root" of A
0116     # find L and D from A using recurrence relations
0117     L = {}
0118     D = {}
0119     for j in range(len(A[0])):
0120         D[j] = A[j][j] - sum([L[j,k]**2 * D[k] for k in range(j)])
0121         for i in range(len(A)):
0122             if i > j:
0123                 L[i,j] = (A[i][j] - sum([L[i,k] * L[j,k] * D[k] for k in range(j)])) / D[j]
0124
0125     L = [[    1.,     0.,     0.,     0.,     0., 0.],
0126          [L[1,0],     1.,     0.,     0.,     0., 0.],
0127          [L[2,0], L[2,1],     1.,     0.,     0., 0.],
0128          [L[3,0], L[3,1], L[3,2],     1.,     0., 0.],
0129          [L[4,0], L[4,1], L[4,2], L[4,1],     1., 0.],
0130          [L[5,0], L[5,1], L[5,2], L[5,1], L[5,0], 1.]]
0131
0132     Dsqrt = [[sqrt(D[0]),         0.,          0.,         0.,         0.,         0.],
0133              [        0., sqrt(D[1]),          0.,         0.,         0.,         0.],
0134              [        0.,         0., sqrt(D[2]),          0.,         0.,         0.],
0135              [        0.,         0.,          0., sqrt(D[3]),         0.,         0.],
0136              [        0.,         0.,          0.,         0., sqrt(D[4]),         0.],
0137              [        0.,         0.,          0.,         0.,         0., sqrt(D[5])]]
0138
0139     return mmult(L, Dsqrt)
0140
0141 matrix = [[   xx,    xy,    xz,    xphix,    xphiy,    xphiz],
0142           [   xy,    yy,    yz,    yphix,    yphiy,    yphiz],
0143           [   xz,    yz,    zz,    zphix,    zphiy,    zphiz],
0144           [xphix, yphix, zphix, phixphix, phixphiy, phixphiz],
0145           [xphiy, yphiy, zphiy, phixphiy, phiyphiy, phiyphiz],
0146           [xphiz, yphiz, zphiz, phixphiz, phiyphiz, phizphiz]]
0147
0148 chomat = cholesky(matrix)
0149
0150 ### Generate correlated random misalignments for all chambers
0151
0152 def random6dof():
0153     randomunit = [gauss(0., 1.), gauss(0., 1.), gauss(0., 1.), gauss(0., 1.), gauss(0., 1.), gauss(0., 1.)]
0154     return mvdot(chomat, randomunit)
0155
0156 misal = {}
0157
0158 for wheel in -2, -1, 0, 1, 2:
0159     for station in 1, 2, 3, 4:
0160         for sector in range(1, 14+1):
0161             if station != 4 and sector > 12: continue
0162
0163             misal["DT", wheel, station, 0, sector] = random6dof()
0164
0165 for endcap in 1, 2:
0166     for station in 1, 2, 3, 4:
0167         for ring in 1, 2, 3:
0168             if station > 1 and ring == 3: continue
0169             for sector in range(1, 36+1):
0170                 if station > 1 and ring == 1 and sector > 18: continue
0171
0172                 misal["CSC", endcap, station, ring, sector] = random6dof()
0173
0174 ### More diagnostics
0175
0176 sum_x = 0.
0177 sum_y = 0.
0178 sum_z = 0.
0179 sum_phix = 0.
0180 sum_phiy = 0.
0181 sum_phiz = 0.
0182
0183 sum_xx = 0.
0184 sum_xy = 0.
0185 sum_xz = 0.
0186 sum_xphix = 0.
0187 sum_xphiy = 0.
0188 sum_xphiz = 0.
0189 sum_yy = 0.
0190 sum_yz = 0.
0191 sum_yphix = 0.
0192 sum_yphiy = 0.
0193 sum_yphiz = 0.
0194 sum_zz = 0.
0195 sum_zphix = 0.
0196 sum_zphiy = 0.
0197 sum_zphiz = 0.
0198 sum_phixphix = 0.
0199 sum_phixphiy = 0.
0200 sum_phixphiz = 0.
0201 sum_phiyphiy = 0.
0202 sum_phiyphiz = 0.
0203 sum_phizphiz = 0.
0204
0205 for xi, yi, zi, phixi, phiyi, phizi in misal.values():
0206     sum_x += xi
0207     sum_y += yi
0208     sum_z += zi
0209     sum_phix += phixi
0210     sum_phiy += phiyi
0211     sum_phiz += phizi
0212
0213     sum_xx += xi*xi
0214     sum_xy += xi*yi
0215     sum_xz += xi*zi
0216     sum_xphix += xi*phixi
0217     sum_xphiy += xi*phiyi
0218     sum_xphiz += xi*phizi
0219     sum_yy += yi*yi
0220     sum_yz += yi*zi
0221     sum_yphix += yi*phixi
0222     sum_yphiy += yi*phiyi
0223     sum_yphiz += yi*phizi
0224     sum_zz += zi*zi
0225     sum_zphix += zi*phixi
0226     sum_zphiy += zi*phiyi
0227     sum_zphiz += zi*phizi
0228     sum_phixphix += phixi*phixi
0229     sum_phixphiy += phixi*phiyi
0230     sum_phixphiz += phixi*phizi
0231     sum_phiyphiy += phiyi*phiyi
0232     sum_phiyphiz += phiyi*phizi
0233     sum_phizphiz += phizi*phizi
0234
0235 ave_x = sum_x/float(len(misal))
0236 ave_y = sum_y/float(len(misal))
0237 ave_z = sum_z/float(len(misal))
0238 ave_phix = sum_phix/float(len(misal))
0239 ave_phiy = sum_phiy/float(len(misal))
0240 ave_phiz = sum_phiz/float(len(misal))
0241
0242 ave_xx = sum_xx/float(len(misal))
0243 ave_xy = sum_xy/float(len(misal))
0244 ave_xz = sum_xz/float(len(misal))
0245 ave_xphix = sum_xphix/float(len(misal))
0246 ave_xphiy = sum_xphiy/float(len(misal))
0247 ave_xphiz = sum_xphiz/float(len(misal))
0248 ave_yy = sum_yy/float(len(misal))
0249 ave_yz = sum_yz/float(len(misal))
0250 ave_yphix = sum_yphix/float(len(misal))
0251 ave_yphiy = sum_yphiy/float(len(misal))
0252 ave_yphiz = sum_yphiz/float(len(misal))
0253 ave_zz = sum_zz/float(len(misal))
0254 ave_zphix = sum_zphix/float(len(misal))
0255 ave_zphiy = sum_zphiy/float(len(misal))
0256 ave_zphiz = sum_zphiz/float(len(misal))
0257 ave_phixphix = sum_phixphix/float(len(misal))
0258 ave_phixphiy = sum_phixphiy/float(len(misal))
0259 ave_phixphiz = sum_phixphiz/float(len(misal))
0260 ave_phiyphiy = sum_phiyphiy/float(len(misal))
0261 ave_phiyphiz = sum_phiyphiz/float(len(misal))
0262 ave_phizphiz = sum_phizphiz/float(len(misal))
0263
0264 print("Estimated covariance matrix from %d chambers (x, y, z, phix, phiy, phiz):" % len(misal))
0265 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xx, ave_xy, ave_xz, ave_xphix, ave_xphiy, ave_xphiz))
0266 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xy, ave_yy, ave_yz, ave_yphix, ave_yphiy, ave_yphiz))
0267 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xz, ave_yz, ave_zz, ave_zphix, ave_zphiy, ave_zphiz))
0268 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xphix, ave_yphix, ave_zphix, ave_phixphix, ave_phixphiy, ave_phixphiz))
0269 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xphiy, ave_yphiy, ave_zphiy, ave_phixphiy, ave_phiyphiy, ave_phiyphiz))
0270 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xphiz, ave_yphiz, ave_zphiz, ave_phixphiz, ave_phiyphiz, ave_phizphiz))
0271 print()
0272
0273 print("Estimated correlation from %d chambers (x, y, z, phix, phiy, phiz):" % len(misal))
0274 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xx/sqrt(ave_xx)/sqrt(ave_xx), ave_xy/sqrt(ave_xx)/sqrt(ave_yy), ave_xz/sqrt(ave_xx)/sqrt(ave_zz), ave_xphix/sqrt(ave_xx)/sqrt(ave_phixphix), ave_xphiy/sqrt(ave_xx)/sqrt(ave_phiyphiy), ave_xphiz/sqrt(ave_xx)/sqrt(ave_phizphiz)))
0275 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xy/sqrt(ave_yy)/sqrt(ave_xx), ave_yy/sqrt(ave_yy)/sqrt(ave_yy), ave_yz/sqrt(ave_yy)/sqrt(ave_zz), ave_yphix/sqrt(ave_yy)/sqrt(ave_phixphix), ave_yphiy/sqrt(ave_yy)/sqrt(ave_phiyphiy), ave_yphiz/sqrt(ave_yy)/sqrt(ave_phizphiz)))
0276 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xz/sqrt(ave_zz)/sqrt(ave_xx), ave_yz/sqrt(ave_zz)/sqrt(ave_yy), ave_zz/sqrt(ave_zz)/sqrt(ave_zz), ave_zphix/sqrt(ave_zz)/sqrt(ave_phixphix), ave_zphiy/sqrt(ave_zz)/sqrt(ave_phiyphiy), ave_zphiz/sqrt(ave_zz)/sqrt(ave_phizphiz)))
0277 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xphix/sqrt(ave_phixphix)/sqrt(ave_xx), ave_yphix/sqrt(ave_phixphix)/sqrt(ave_yy), ave_zphix/sqrt(ave_phixphix)/sqrt(ave_zz), ave_phixphix/sqrt(ave_phixphix)/sqrt(ave_phixphix), ave_phixphiy/sqrt(ave_phixphix)/sqrt(ave_phiyphiy), ave_phixphiz/sqrt(ave_phixphix)/sqrt(ave_phizphiz)))
0278 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xphiy/sqrt(ave_phiyphiy)/sqrt(ave_xx), ave_yphiy/sqrt(ave_phiyphiy)/sqrt(ave_yy), ave_zphiy/sqrt(ave_phiyphiy)/sqrt(ave_zz), ave_phixphiy/sqrt(ave_phiyphiy)/sqrt(ave_phixphix), ave_phiyphiy/sqrt(ave_phiyphiy)/sqrt(ave_phiyphiy), ave_phiyphiz/sqrt(ave_phiyphiy)/sqrt(ave_phizphiz)))
0279 print("%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (ave_xphiz/sqrt(ave_phizphiz)/sqrt(ave_xx), ave_yphiz/sqrt(ave_phizphiz)/sqrt(ave_yy), ave_zphiz/sqrt(ave_phizphiz)/sqrt(ave_zz), ave_phixphiz/sqrt(ave_phizphiz)/sqrt(ave_phixphix), ave_phiyphiz/sqrt(ave_phizphiz)/sqrt(ave_phiyphiy), ave_phizphiz/sqrt(ave_phizphiz)/sqrt(ave_phizphiz)))
0280 print()
0281
0282 ### Delete all three files at once to make sure the user never sees
0283 ### stale data (e.g. from a stopped process due to failed conversion)
0284
0285 if os.path.exists(outputName + ".xml"):
0286     os.unlink(outputName + ".xml")
0287 if os.path.exists(outputName + "_convert_cfg.py"):
0288     os.unlink(outputName + "_convert_cfg.py")
0289 if os.path.exists(outputName + ".db"):
0290     os.unlink(outputName + ".db")
0291 if os.path.exists(outputName + "_correlations.txt"):
0292     os.unlink(outputName + "_correlations.txt")
0293
0294 ### Print out the list of correlations
0295
0296 txtfile = file(outputName + "_correlations.txt", "w")
0297 for wheel in -2, -1, 0, 1, 2:
0298     for station in 1, 2, 3, 4:
0299         for sector in range(1, 14+1):
0300             if station != 4 and sector > 12: continue
0301             txtfile.write("DT %(wheel)d %(station)d %(sector)d %(xx)g %(xy)g %(xz)g %(xphix)g %(xphiy)g %(xphiz)g %(yy)g %(yz)g %(yphix)g %(yphiy)g %(yphiz)g %(zz)g %(zphix)g %(zphiy)g %(zphiz)g %(phixphix)g %(phixphiy)g %(phixphiz)g %(phiyphiy)g %(phiyphiz)g %(phizphiz)g\n" % vars())
0302
0303 for endcap in 1, 2:
0304     for station in 1, 2, 3, 4:
0305         for ring in 1, 2, 3:
0306             if station > 1 and ring == 3: continue
0307             for sector in range(1, 36+1):
0308                 if station > 1 and ring == 1 and sector > 18: continue
0309                 txtfile.write("CSC %(endcap)d %(station)d %(ring)d %(sector)d %(xx)g %(xy)g %(xz)g %(xphix)g %(xphiy)g %(xphiz)g %(yy)g %(yz)g %(yphix)g %(yphiy)g %(yphiz)g %(zz)g %(zphix)g %(zphiy)g %(zphiz)g %(phixphix)g %(phixphiy)g %(phixphiz)g %(phiyphiy)g %(phiyphiz)g %(phizphiz)g\n" % vars())
0310
0311 ### Make an XML representation of the misalignment
0312
0313 xmlfile = file(outputName + ".xml", "w")
0314 xmlfile.write("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
0315 xmlfile.write("<?xml-stylesheet type=\"text/xml\" href=\"MuonAlignment.xsl\"?>\n")
0316 xmlfile.write("<MuonAlignment>\n\n")
0317
0318 for (system, whendcap, station, ring, sector), (xi, yi, zi, phixi, phiyi, phizi) in misal.items():
0319     if system == "DT": wheel = whendcap
0320     if system == "CSC": endcap = whendcap
0321
0322     rot = rotation[system, whendcap, station, ring, sector]
0323     localape = [[xx, xy, xz], [xy, yy, yz], [xz, yz, zz]]
0324     globalape = mmult(rot, mmult(localape, mtrans(rot)))
0325     globalxx = globalape[0][0]
0326     globalxy = globalape[0][1]
0327     globalxz = globalape[0][2]
0328     globalyy = globalape[1][1]
0329     globalyz = globalape[1][2]
0330     globalzz = globalape[2][2]
0331
0332     xmlfile.write("<operation>\n")
0333
0334     if system == "DT":
0335         xmlfile.write("    <DTChamber wheel=\"%(wheel)d\" station=\"%(station)d\" sector=\"%(sector)d\" />\n" % vars())
0336     if system == "CSC":
0337         xmlfile.write("    <CSCChamber endcap=\"%(endcap)d\" station=\"%(station)d\" ring=\"%(ring)d\" chamber=\"%(sector)d\" />\n" % vars())
0338
0339         # ME1/1a is called "ring 4", but it should always get exactly the same alignment constants as the corresponding ME1/1b ("ring 1")
0340         if (station, ring) == (1, 1):
0341             xmlfile.write("    <CSCChamber endcap=\"%(endcap)d\" station=\"%(station)d\" ring=\"4\" chamber=\"%(sector)d\" />\n" % vars())
0342
0343     xmlfile.write("    <setposition relativeto=\"ideal\" x=\"%(xi)g\" y=\"%(yi)g\" z=\"%(zi)g\" phix=\"%(phixi)g\" phiy=\"%(phiyi)g\" phiz=\"%(phizi)g\" />\n" % vars())
0344     xmlfile.write("    <setape xx=\"%(globalxx)g\" xy=\"%(globalxy)g\" xz=\"%(globalxz)g\" yy=\"%(globalyy)g\" yz=\"%(globalyz)g\" zz=\"%(globalzz)g\" />\n" % vars())
0345     xmlfile.write("</operation>\n\n")
0346
0347 xmlfile.write("</MuonAlignment>\n")
0348 xmlfile.close()
0349
0350 ### Convert it to an SQLite file with CMSSW
0351
0352 cfgfile = file(outputName + "_convert_cfg.py", "w")
0353
0354 cfgfile.write("""import FWCore.ParameterSet.Config as cms
0355
0356 process = cms.Process("CONVERT")
0357 process.source = cms.Source("EmptySource")
0358 process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1))
0359
0360 process.load("Configuration.StandardSequences.GeometryDB_cff")
0361 process.load("Geometry.MuonNumbering.muonNumberingInitialization_cfi")
0362
0363 process.MuonGeometryDBConverter = cms.EDAnalyzer("MuonGeometryDBConverter",
0364                                                  input = cms.string("xml"),
0365                                                  fileName = cms.string("%(outputName)s.xml"),
0366                                                  shiftErr = cms.double(1000.),
0367                                                  angleErr = cms.double(6.28),
0368
0369                                                  output = cms.string("db")
0370                                                  )
0371
0372 process.load("CondCore.DBCommon.CondDBSetup_cfi")
0373 process.PoolDBOutputService = cms.Service("PoolDBOutputService",
0374                                           process.CondDBSetup,
0375                                           connect = cms.string("sqlite_file:%(outputName)s.db"),
0376                                           toPut = cms.VPSet(cms.PSet(record = cms.string("DTAlignmentRcd"), tag = cms.string("DTAlignmentRcd")),
0377                                                             cms.PSet(record = cms.string("DTAlignmentErrorExtendedRcd"), tag = cms.string("DTAlignmentErrorExtendedRcd")),
0378                                                             cms.PSet(record = cms.string("CSCAlignmentRcd"), tag = cms.string("CSCAlignmentRcd")),
0379                                                             cms.PSet(record = cms.string("CSCAlignmentErrorExtendedRcd"), tag = cms.string("CSCAlignmentErrorExtendedRcd"))))
0380
0381 process.Path = cms.Path(process.MuonGeometryDBConverter)
0382 """ % vars())
0383
0384 print("To create an SQLite file for this geometry (%(outputName)s.db), run the following:" % vars())
0385 print()
0386 os.system("echo cmsRun %s_convert_cfg.py" % outputName)
0387