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