CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
makeMuonMisalignmentScenario.py
Go to the documentation of this file.
1 from optparse import OptionParser
2 from random import gauss
3 from math import sqrt
4 import os
5 execfile("Alignment/MuonAlignment/data/idealTransformation.py")
6 
7 ### Get variances and covariances from the commandline
8 
9 parser = OptionParser(usage="Usage: python %prog outputName [options] (default is unit matrix times 1e-15)")
10 
11 parser.add_option("--xx", dest="xx", help="variance of x (cm*cm)", default="1e-15")
12 parser.add_option("--xy", dest="xy", help="covariance of x and y (cm*cm)", default="0")
13 parser.add_option("--xz", dest="xz", help="covariance of x and z (cm*cm)", default="0")
14 parser.add_option("--xphix", dest="xphix", help="covariance of x and phix (cm*rad)", default="0")
15 parser.add_option("--xphiy", dest="xphiy", help="covariance of x and phiy (cm*rad)", default="0")
16 parser.add_option("--xphiz", dest="xphiz", help="covariance of x and phiz (cm*rad)", default="0")
17 
18 parser.add_option("--yy", dest="yy", help="variance of y (cm*cm)", default="1e-15")
19 parser.add_option("--yz", dest="yz", help="covariance of y and z (cm*cm)", default="0")
20 parser.add_option("--yphix", dest="yphix", help="covariance of y and phix (cm*rad)", default="0")
21 parser.add_option("--yphiy", dest="yphiy", help="covariance of y and phiy (cm*rad)", default="0")
22 parser.add_option("--yphiz", dest="yphiz", help="covariance of y and phiz (cm*rad)", default="0")
23 
24 parser.add_option("--zz", dest="zz", help="variance of z (cm*cm)", default="1e-15")
25 parser.add_option("--zphix", dest="zphix", help="covariance of z and phix (cm*rad)", default="0")
26 parser.add_option("--zphiy", dest="zphiy", help="covariance of z and phiy (cm*rad)", default="0")
27 parser.add_option("--zphiz", dest="zphiz", help="covariance of z and phiz (cm*rad)", default="0")
28 
29 parser.add_option("--phixphix", dest="phixphix", help="variance of phix (rad*rad)", default="1e-15")
30 parser.add_option("--phixphiy", dest="phixphiy", help="covariance of phix and phiy (rad*rad)", default="0")
31 parser.add_option("--phixphiz", dest="phixphiz", help="covariance of phix and phiz (rad*rad)", default="0")
32 
33 parser.add_option("--phiyphiy", dest="phiyphiy", help="variance of phiy (rad*rad)", default="1e-15")
34 parser.add_option("--phiyphiz", dest="phiyphiz", help="covariance of phiy and phiz (rad*rad)", default="0")
35 
36 parser.add_option("--phizphiz", dest="phizphiz", help="variance of phiz (rad*rad)", default="1e-15")
37 
38 parser.add_option("-f", dest="force", help="force overwrite of output files", action="store_true")
39 
40 options, args = parser.parse_args()
41 if args is None or len(args) != 1:
42  parser.print_help()
43  exit(-1)
44 outputName = args[0]
45 
46 if not options.force:
47  if os.path.exists(outputName + ".xml"):
48  raise Exception(outputName + ".xml exists!")
49  if os.path.exists(outputName + "_convert_cfg.py"):
50  raise Exception(outputName + "_convert_cfg.py exists!")
51  if os.path.exists(outputName + ".db"):
52  raise Exception(outputName + ".db exists!")
53  if os.path.exists(outputName + "_correlations.txt"):
54  raise Exception(outputName + "_correlations.txt exists!")
55 
56 components = "xx", "xy", "xz", "xphix", "xphiy", "xphiz", "yy", "yz", "yphix", "yphiy", "yphiz", "zz", "zphix", "zphiy", "zphiz", "phixphix", "phixphiy", "phixphiz", "phiyphiy", "phiyphiz", "phizphiz"
57 for component in components:
58  exec("%s = float(options.%s)" % (component, component))
59 
60 ### Print out user's choices as diagnostics
61 
62 print "Spread in each parameter: x %g mm" % (sqrt(xx)*10.)
63 print " y %g mm" % (sqrt(yy)*10.)
64 print " z %g mm" % (sqrt(zz)*10.)
65 print " phix %g mrad" % (sqrt(phixphix)*1000.)
66 print " phiy %g mrad" % (sqrt(phiyphiy)*1000.)
67 print " phiz %g mrad" % (sqrt(phizphiz)*1000.)
68 print
69 
70 print "Covariance matrix (x, y, z, phix, phiy, phiz):"
71 print "%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xx, xy, xz, xphix, xphiy, xphiz)
72 print "%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xy, yy, yz, yphix, yphiy, yphiz)
73 print "%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xz, yz, zz, zphix, zphiy, zphiz)
74 print "%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphix, yphix, zphix, phixphix, phixphiy, phixphiz)
75 print "%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphiy, yphiy, zphiy, phixphiy, phiyphiy, phiyphiz)
76 print "%11.8f %11.8f %11.8f %11.8f %11.8f %11.8f" % (xphiz, yphiz, zphiz, phixphiz, phiyphiz, phizphiz)
77 print
78 
79 print "Correlation (x, y, z, phix, phiy, phiz):"
80 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))
81 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))
82 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))
83 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))
84 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))
85 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))
86 print
87 
88 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)), \
89  abs(yz/sqrt(yy)/sqrt(zz)), abs(yphix/sqrt(yy)/sqrt(phixphix)), abs(yphiy/sqrt(yy)/sqrt(phiyphiy)), abs(yphiz/sqrt(yy)/sqrt(phizphiz)),
90  abs(zphix/sqrt(zz)/sqrt(phixphix)), abs(zphiy/sqrt(zz)/sqrt(phiyphiy)), abs(zphiz/sqrt(zz)/sqrt(phizphiz)),
91  abs(phixphiy/sqrt(phixphix)/sqrt(phiyphiy)), abs(phixphiz/sqrt(phixphix)/sqrt(phizphiz)),
92  abs(phiyphiz/sqrt(phiyphiy)/sqrt(phizphiz))]:
93  if correlation_coefficient > 1.:
94  raise Exception("Correlations must not be larger than one!")
95 
96 ### Some useful mathematical transformations (why don't we have access to numpy?)
97 
98 def mmult(a, b):
99  """Matrix multiplication: mmult([[11, 12], [21, 22]], [[-1, 0], [0, 1]]) returns [[-11, 12], [-21, 22]]"""
100  return [[sum([i*j for i, j in zip(row, col)]) for col in zip(*b)] for row in a]
101 
102 def mvdot(m, v):
103  """Applies matrix m to vector v: mvdot([[-1, 0], [0, 1]], [12, 55]) returns [-12, 55]"""
104  return [i[0] for i in mmult(m, [[vi] for vi in v])]
105 
106 def mtrans(a):
107  """Matrix transposition: mtrans([[11, 12], [21, 22]]) returns [[11, 21], [12, 22]]"""
108  return [[a[j][i] for j in range(len(a[i]))] for i in range(len(a))]
109 
110 def cholesky(A):
111  """Cholesky decomposition of the correlation matrix to properly normalize the transformed random deviates"""
112 
113  # 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
114  # find L and D from A using recurrence relations
115  L = {}
116  D = {}
117  for j in range(len(A[0])):
118  D[j] = A[j][j] - sum([L[j,k]**2 * D[k] for k in range(j)])
119  for i in range(len(A)):
120  if i > j:
121  L[i,j] = (A[i][j] - sum([L[i,k] * L[j,k] * D[k] for k in range(j)])) / D[j]
122 
123  L = [[ 1., 0., 0., 0., 0., 0.],
124  [L[1,0], 1., 0., 0., 0., 0.],
125  [L[2,0], L[2,1], 1., 0., 0., 0.],
126  [L[3,0], L[3,1], L[3,2], 1., 0., 0.],
127  [L[4,0], L[4,1], L[4,2], L[4,1], 1., 0.],
128  [L[5,0], L[5,1], L[5,2], L[5,1], L[5,0], 1.]]
129 
130  Dsqrt = [[sqrt(D[0]), 0., 0., 0., 0., 0.],
131  [ 0., sqrt(D[1]), 0., 0., 0., 0.],
132  [ 0., 0., sqrt(D[2]), 0., 0., 0.],
133  [ 0., 0., 0., sqrt(D[3]), 0., 0.],
134  [ 0., 0., 0., 0., sqrt(D[4]), 0.],
135  [ 0., 0., 0., 0., 0., sqrt(D[5])]]
136 
137  return mmult(L, Dsqrt)
138 
139 matrix = [[ xx, xy, xz, xphix, xphiy, xphiz],
140  [ xy, yy, yz, yphix, yphiy, yphiz],
141  [ xz, yz, zz, zphix, zphiy, zphiz],
142  [xphix, yphix, zphix, phixphix, phixphiy, phixphiz],
143  [xphiy, yphiy, zphiy, phixphiy, phiyphiy, phiyphiz],
144  [xphiz, yphiz, zphiz, phixphiz, phiyphiz, phizphiz]]
145 
146 chomat = cholesky(matrix)
147 
148 ### Generate correlated random misalignments for all chambers
149 
150 def random6dof():
151  randomunit = [gauss(0., 1.), gauss(0., 1.), gauss(0., 1.), gauss(0., 1.), gauss(0., 1.), gauss(0., 1.)]
152  return mvdot(chomat, randomunit)
153 
154 misal = {}
155 
156 for wheel in -2, -1, 0, 1, 2:
157  for station in 1, 2, 3, 4:
158  for sector in range(1, 14+1):
159  if station != 4 and sector > 12: continue
160 
161  misal["DT", wheel, station, 0, sector] = random6dof()
162 
163 for endcap in 1, 2:
164  for station in 1, 2, 3, 4:
165  for ring in 1, 2, 3:
166  if station > 1 and ring == 3: continue
167  for sector in range(1, 36+1):
168  if station > 1 and ring == 1 and sector > 18: continue
169 
170  misal["CSC", endcap, station, ring, sector] = random6dof()
171 
172 ### More diagnostics
173 
174 sum_x = 0.
175 sum_y = 0.
176 sum_z = 0.
177 sum_phix = 0.
178 sum_phiy = 0.
179 sum_phiz = 0.
180 
181 sum_xx = 0.
182 sum_xy = 0.
183 sum_xz = 0.
184 sum_xphix = 0.
185 sum_xphiy = 0.
186 sum_xphiz = 0.
187 sum_yy = 0.
188 sum_yz = 0.
189 sum_yphix = 0.
190 sum_yphiy = 0.
191 sum_yphiz = 0.
192 sum_zz = 0.
193 sum_zphix = 0.
194 sum_zphiy = 0.
195 sum_zphiz = 0.
196 sum_phixphix = 0.
197 sum_phixphiy = 0.
198 sum_phixphiz = 0.
199 sum_phiyphiy = 0.
200 sum_phiyphiz = 0.
201 sum_phizphiz = 0.
202 
203 for xi, yi, zi, phixi, phiyi, phizi in misal.values():
204  sum_x += xi
205  sum_y += yi
206  sum_z += zi
207  sum_phix += phixi
208  sum_phiy += phiyi
209  sum_phiz += phizi
210 
211  sum_xx += xi*xi
212  sum_xy += xi*yi
213  sum_xz += xi*zi
214  sum_xphix += xi*phixi
215  sum_xphiy += xi*phiyi
216  sum_xphiz += xi*phizi
217  sum_yy += yi*yi
218  sum_yz += yi*zi
219  sum_yphix += yi*phixi
220  sum_yphiy += yi*phiyi
221  sum_yphiz += yi*phizi
222  sum_zz += zi*zi
223  sum_zphix += zi*phixi
224  sum_zphiy += zi*phiyi
225  sum_zphiz += zi*phizi
226  sum_phixphix += phixi*phixi
227  sum_phixphiy += phixi*phiyi
228  sum_phixphiz += phixi*phizi
229  sum_phiyphiy += phiyi*phiyi
230  sum_phiyphiz += phiyi*phizi
231  sum_phizphiz += phizi*phizi
232 
233 ave_x = sum_x/float(len(misal))
234 ave_y = sum_y/float(len(misal))
235 ave_z = sum_z/float(len(misal))
236 ave_phix = sum_phix/float(len(misal))
237 ave_phiy = sum_phiy/float(len(misal))
238 ave_phiz = sum_phiz/float(len(misal))
239 
240 ave_xx = sum_xx/float(len(misal))
241 ave_xy = sum_xy/float(len(misal))
242 ave_xz = sum_xz/float(len(misal))
243 ave_xphix = sum_xphix/float(len(misal))
244 ave_xphiy = sum_xphiy/float(len(misal))
245 ave_xphiz = sum_xphiz/float(len(misal))
246 ave_yy = sum_yy/float(len(misal))
247 ave_yz = sum_yz/float(len(misal))
248 ave_yphix = sum_yphix/float(len(misal))
249 ave_yphiy = sum_yphiy/float(len(misal))
250 ave_yphiz = sum_yphiz/float(len(misal))
251 ave_zz = sum_zz/float(len(misal))
252 ave_zphix = sum_zphix/float(len(misal))
253 ave_zphiy = sum_zphiy/float(len(misal))
254 ave_zphiz = sum_zphiz/float(len(misal))
255 ave_phixphix = sum_phixphix/float(len(misal))
256 ave_phixphiy = sum_phixphiy/float(len(misal))
257 ave_phixphiz = sum_phixphiz/float(len(misal))
258 ave_phiyphiy = sum_phiyphiy/float(len(misal))
259 ave_phiyphiz = sum_phiyphiz/float(len(misal))
260 ave_phizphiz = sum_phizphiz/float(len(misal))
261 
262 print "Estimated covariance matrix from %d chambers (x, y, z, phix, phiy, phiz):" % len(misal)
263 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)
264 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)
265 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)
266 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)
267 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)
268 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)
269 print
270 
271 print "Estimated correlation from %d chambers (x, y, z, phix, phiy, phiz):" % len(misal)
272 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))
273 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))
274 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))
275 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))
276 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))
277 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))
278 print
279 
280 ### Delete all three files at once to make sure the user never sees
281 ### stale data (e.g. from a stopped process due to failed conversion)
282 
283 if os.path.exists(outputName + ".xml"):
284  os.unlink(outputName + ".xml")
285 if os.path.exists(outputName + "_convert_cfg.py"):
286  os.unlink(outputName + "_convert_cfg.py")
287 if os.path.exists(outputName + ".db"):
288  os.unlink(outputName + ".db")
289 if os.path.exists(outputName + "_correlations.txt"):
290  os.unlink(outputName + "_correlations.txt")
291 
292 ### Print out the list of correlations
293 
294 txtfile = file(outputName + "_correlations.txt", "w")
295 for wheel in -2, -1, 0, 1, 2:
296  for station in 1, 2, 3, 4:
297  for sector in range(1, 14+1):
298  if station != 4 and sector > 12: continue
299  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())
300 
301 for endcap in 1, 2:
302  for station in 1, 2, 3, 4:
303  for ring in 1, 2, 3:
304  if station > 1 and ring == 3: continue
305  for sector in range(1, 36+1):
306  if station > 1 and ring == 1 and sector > 18: continue
307  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())
308 
309 ### Make an XML representation of the misalignment
310 
311 xmlfile = file(outputName + ".xml", "w")
312 xmlfile.write("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
313 xmlfile.write("<?xml-stylesheet type=\"text/xml\" href=\"MuonAlignment.xsl\"?>\n")
314 xmlfile.write("<MuonAlignment>\n\n")
315 
316 for (system, whendcap, station, ring, sector), (xi, yi, zi, phixi, phiyi, phizi) in misal.items():
317  if system == "DT": wheel = whendcap
318  if system == "CSC": endcap = whendcap
319 
320  rot = rotation[system, whendcap, station, ring, sector]
321  localape = [[xx, xy, xz], [xy, yy, yz], [xz, yz, zz]]
322  globalape = mmult(rot, mmult(localape, mtrans(rot)))
323  globalxx = globalape[0][0]
324  globalxy = globalape[0][1]
325  globalxz = globalape[0][2]
326  globalyy = globalape[1][1]
327  globalyz = globalape[1][2]
328  globalzz = globalape[2][2]
329 
330  xmlfile.write("<operation>\n")
331 
332  if system == "DT":
333  xmlfile.write(" <DTChamber wheel=\"%(wheel)d\" station=\"%(station)d\" sector=\"%(sector)d\" />\n" % vars())
334  if system == "CSC":
335  xmlfile.write(" <CSCChamber endcap=\"%(endcap)d\" station=\"%(station)d\" ring=\"%(ring)d\" chamber=\"%(sector)d\" />\n" % vars())
336 
337  # ME1/1a is called "ring 4", but it should always get exactly the same alignment constants as the corresponding ME1/1b ("ring 1")
338  if (station, ring) == (1, 1):
339  xmlfile.write(" <CSCChamber endcap=\"%(endcap)d\" station=\"%(station)d\" ring=\"4\" chamber=\"%(sector)d\" />\n" % vars())
340 
341  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())
342  xmlfile.write(" <setape xx=\"%(globalxx)g\" xy=\"%(globalxy)g\" xz=\"%(globalxz)g\" yy=\"%(globalyy)g\" yz=\"%(globalyz)g\" zz=\"%(globalzz)g\" />\n" % vars())
343  xmlfile.write("</operation>\n\n")
344 
345 xmlfile.write("</MuonAlignment>\n")
346 xmlfile.close()
347 
348 ### Convert it to an SQLite file with CMSSW
349 
350 cfgfile = file(outputName + "_convert_cfg.py", "w")
351 
352 cfgfile.write("""import FWCore.ParameterSet.Config as cms
353 
354 process = cms.Process("CONVERT")
355 process.source = cms.Source("EmptySource")
356 process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1))
357 
358 process.load("Configuration.StandardSequences.GeometryDB_cff")
359 process.load("Geometry.MuonNumbering.muonNumberingInitialization_cfi")
360 
361 process.MuonGeometryDBConverter = cms.EDAnalyzer("MuonGeometryDBConverter",
362  input = cms.string("xml"),
363  fileName = cms.string("%(outputName)s.xml"),
364  shiftErr = cms.double(1000.),
365  angleErr = cms.double(6.28),
366 
367  output = cms.string("db")
368  )
369 
370 process.load("CondCore.DBCommon.CondDBSetup_cfi")
371 process.PoolDBOutputService = cms.Service("PoolDBOutputService",
372  process.CondDBSetup,
373  connect = cms.string("sqlite_file:%(outputName)s.db"),
374  toPut = cms.VPSet(cms.PSet(record = cms.string("DTAlignmentRcd"), tag = cms.string("DTAlignmentRcd")),
375  cms.PSet(record = cms.string("DTAlignmentErrorExtendedRcd"), tag = cms.string("DTAlignmentErrorExtendedRcd")),
376  cms.PSet(record = cms.string("CSCAlignmentRcd"), tag = cms.string("CSCAlignmentRcd")),
377  cms.PSet(record = cms.string("CSCAlignmentErrorExtendedRcd"), tag = cms.string("CSCAlignmentErrorExtendedRcd"))))
378 
379 process.Path = cms.Path(process.MuonGeometryDBConverter)
380 """ % vars())
381 
382 print "To create an SQLite file for this geometry (%(outputName)s.db), run the following:" % vars()
383 print
384 os.system("echo cmsRun %s_convert_cfg.py" % outputName)
385 
makeMuonMisalignmentScenario.mmult
def mmult
Print out user&#39;s choices as diagnostics.
Definition: makeMuonMisalignmentScenario.py:98
makeMuonMisalignmentScenario.random6dof
def random6dof
Generate correlated random misalignments for all chambers.
Definition: makeMuonMisalignmentScenario.py:150
ComparisonHelper::zip
OutputIterator zip(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2, OutputIterator result, Compare comp)
Definition: L1TStage2CaloLayer1.h:38
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Generated for CMSSW Reference Manual by doxygen 1.8.5