MuonAlignment/python/makeMuonMisalignmentScenario.py

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