CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
LASBarrelAlgorithm.cc
Go to the documentation of this file.
1 
2 
3 #include "Alignment/LaserAlignment/interface/LASBarrelAlgorithm.h"
4 
5 // this is ugly but we need it for Minuit
6 // until a better solution is at hand
7 LASGlobalData<LASCoordinateSet>* aMeasuredCoordinates;
8 LASGlobalData<LASCoordinateSet>* aNominalCoordinates;
9 
13 LASBarrelAlgorithm::LASBarrelAlgorithm() { minuit = new TMinuit(52); }
14 
20 LASBarrelAlignmentParameterSet LASBarrelAlgorithm::CalculateParameters(
21  LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
22  std::cout << " [LASBarrelAlgorithm::CalculateParameters] -- Starting." << std::endl;
23 
25  // for testing..
26  //ReadMisalignmentFromFile( "misalign-var.txt", measuredCoordinates, nominalCoordinates );
28 
29  // statics for minuit
30  aMeasuredCoordinates = &measuredCoordinates;
31  aNominalCoordinates = &nominalCoordinates;
32 
33  // minimizer and variables for it
34  minuit->SetFCN(fcn);
35  double arglist[10];
36  int _ierflg = 0;
37 
38  // toggle minuit blabla
39  arglist[0] = -1;
40  minuit->mnexcm("SET PRI", arglist, 1, _ierflg);
41 
42  // set par errors
43  arglist[0] = 1;
44  minuit->mnexcm("SET ERR", arglist, 1, _ierflg);
45 
46  //
47  // define 52 parameters
48  //
49 
50  // start values: to be evacuated to cfg
51  static float _vstart[52] = {
52  0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TIB+
53  0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TIB-
54  0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TOB+
55  0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TOB-
56  0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TEC+
57  0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TEC-
58  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, // beams 0-3
59  0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 // beams 4-7
60  };
61 
63  // ReadStartParametersFromFile( "startParameters.txt", _vstart ); // debug
65 
66  // step sizes: to be tuned, to be evacuated to cfg
67  static float _vstep[52] = {
68  0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TIB+
69  0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TIB-
70  0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TOB+
71  0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TOB-
72  0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TEC+
73  0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TEC-
74  0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, // beams 0-3
75  0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001 // beams 4-7
76  };
77 
78  // subdetector parameters for TIB+:
79 
80  // rotation around z of first end face
81  minuit->mnparm(0, "subRot1TIB+", _vstart[0], _vstep[0], 0, 0, _ierflg);
82  // rotation around z of second end face
83  minuit->mnparm(1, "subRot2TIB+", _vstart[1], _vstep[1], 0, 0, _ierflg);
84  // translation along x of first end face
85  minuit->mnparm(2, "subTransX1TIB+", _vstart[2], _vstep[2], 0, 0, _ierflg);
86  // translation along x of second end face
87  minuit->mnparm(3, "subTransX2TIB+", _vstart[3], _vstep[3], 0, 0, _ierflg);
88  // translation along y of first end face
89  minuit->mnparm(4, "subTransY1TIB+", _vstart[4], _vstep[4], 0, 0, _ierflg);
90  // translation along y of second end face
91  minuit->mnparm(5, "subTransY2TIB+", _vstart[5], _vstep[5], 0, 0, _ierflg);
92 
93  // subdetector parameters for TIB-:
94 
95  // rotation around z of first end face
96  minuit->mnparm(6, "subRot1TIB-", _vstart[6], _vstep[6], 0, 0, _ierflg);
97  // rotation around z of second end face
98  minuit->mnparm(7, "subRot2TIB-", _vstart[7], _vstep[7], 0, 0, _ierflg);
99  // translation along x of first end face
100  minuit->mnparm(8, "subTransX1TIB-", _vstart[8], _vstep[8], 0, 0, _ierflg);
101  // translation along x of second end face
102  minuit->mnparm(9, "subTransX2TIB-", _vstart[9], _vstep[9], 0, 0, _ierflg);
103  // translation along y of first end face
104  minuit->mnparm(10, "subTransY1TIB-", _vstart[10], _vstep[10], 0, 0, _ierflg);
105  // translation along y of second end face
106  minuit->mnparm(11, "subTransY2TIB-", _vstart[11], _vstep[11], 0, 0, _ierflg);
107 
108  // subdetector parameters for TOB+:
109 
110  // rotation around z of first end face
111  minuit->mnparm(12, "subRot1TOB+", _vstart[12], _vstep[12], 0, 0, _ierflg);
112  // rotation around z of second end face
113  minuit->mnparm(13, "subRot2TOB+", _vstart[13], _vstep[13], 0, 0, _ierflg);
114  // translation along x of first end face
115  minuit->mnparm(14, "subTransX1TOB+", _vstart[14], _vstep[14], 0, 0, _ierflg);
116  // translation along x of second end face
117  minuit->mnparm(15, "subTransX2TOB+", _vstart[15], _vstep[15], 0, 0, _ierflg);
118  // translation along y of first end face
119  minuit->mnparm(16, "subTransY1TOB+", _vstart[16], _vstep[16], 0, 0, _ierflg);
120  // translation along y of second end face
121  minuit->mnparm(17, "subTransY2TOB+", _vstart[17], _vstep[17], 0, 0, _ierflg);
122 
123  // subdetector parameters for TOB-:
124 
125  // rotation around z of first end face
126  minuit->mnparm(18, "subRot1TOB-", _vstart[18], _vstep[18], 0, 0, _ierflg);
127  // rotation around z of second end face
128  minuit->mnparm(19, "subRot2TOB-", _vstart[19], _vstep[19], 0, 0, _ierflg);
129  // translation along x of first end face
130  minuit->mnparm(20, "subTransX1TOB-", _vstart[20], _vstep[20], 0, 0, _ierflg);
131  // translation along x of second end face
132  minuit->mnparm(21, "subTransX2TOB-", _vstart[21], _vstep[21], 0, 0, _ierflg);
133  // translation along y of first end face
134  minuit->mnparm(22, "subTransY1TOB-", _vstart[22], _vstep[22], 0, 0, _ierflg);
135  // translation along y of second end face
136  minuit->mnparm(23, "subTransY2TOB-", _vstart[23], _vstep[23], 0, 0, _ierflg);
137 
138  // subdetector parameters for TEC+:
139 
140  // rotation around z of first end face
141  minuit->mnparm(24, "subRot1TEC+", _vstart[24], _vstep[24], 0, 0, _ierflg);
142  // rotation around z of second end face
143  minuit->mnparm(25, "subRot2TEC+", _vstart[25], _vstep[25], 0, 0, _ierflg);
144  // translation along x of first end face
145  minuit->mnparm(26, "subTransX1TEC+", _vstart[26], _vstep[26], 0, 0, _ierflg);
146  // translation along x of second end face
147  minuit->mnparm(27, "subTransX2TEC+", _vstart[27], _vstep[27], 0, 0, _ierflg);
148  // translation along y of first end face
149  minuit->mnparm(28, "subTransY1TEC+", _vstart[28], _vstep[28], 0, 0, _ierflg);
150  // translation along y of second end face
151  minuit->mnparm(29, "subTransY2TEC+", _vstart[29], _vstep[29], 0, 0, _ierflg);
152 
153  // subdetector parameters for TEC-:
154 
155  // rotation around z of first end face
156  minuit->mnparm(30, "subRot1TEC-", _vstart[30], _vstep[30], 0, 0, _ierflg);
157  // rotation around z of second end face
158  minuit->mnparm(31, "subRot2TEC-", _vstart[31], _vstep[31], 0, 0, _ierflg);
159  // translation along x of first end face
160  minuit->mnparm(32, "subTransX1TEC-", _vstart[32], _vstep[32], 0, 0, _ierflg);
161  // translation along x of second end face
162  minuit->mnparm(33, "subTransX2TEC-", _vstart[33], _vstep[33], 0, 0, _ierflg);
163  // translation along y of first end face
164  minuit->mnparm(34, "subTransY1TEC-", _vstart[34], _vstep[34], 0, 0, _ierflg);
165  // translation along y of second end face
166  minuit->mnparm(35, "subTransY2TEC-", _vstart[35], _vstep[35], 0, 0, _ierflg);
167 
168  // beam parameters (+-z pairs, duplicated for beams 0-7):
169 
170  // rotation around z at zt1
171  minuit->mnparm(36, "beamRot1Beam0", _vstart[36], _vstep[36], 0, 0, _ierflg);
172  // rotation around z at zt2
173  minuit->mnparm(37, "beamRot2Beam0", _vstart[37], _vstep[37], 0, 0, _ierflg);
174 
175  // rotation around z at zt1
176  minuit->mnparm(38, "beamRot1Beam1", _vstart[38], _vstep[38], 0, 0, _ierflg);
177  // rotation around z at zt2
178  minuit->mnparm(39, "beamRot2Beam1", _vstart[39], _vstep[39], 0, 0, _ierflg);
179 
180  // rotation around z at zt1
181  minuit->mnparm(40, "beamRot1Beam2", _vstart[40], _vstep[40], 0, 0, _ierflg);
182  // rotation around z at zt2
183  minuit->mnparm(41, "beamRot2Beam2", _vstart[41], _vstep[41], 0, 0, _ierflg);
184 
185  // rotation around z at zt1
186  minuit->mnparm(42, "beamRot1Beam3", _vstart[42], _vstep[42], 0, 0, _ierflg);
187  // rotation around z at zt2
188  minuit->mnparm(43, "beamRot2Beam3", _vstart[43], _vstep[43], 0, 0, _ierflg);
189 
190  // rotation around z at zt1
191  minuit->mnparm(44, "beamRot1Beam4", _vstart[44], _vstep[44], 0, 0, _ierflg);
192  // rotation around z at zt2
193  minuit->mnparm(45, "beamRot2Beam4", _vstart[45], _vstep[45], 0, 0, _ierflg);
194 
195  // rotation around z at zt1
196  minuit->mnparm(46, "beamRot1Beam5", _vstart[46], _vstep[46], 0, 0, _ierflg);
197  // rotation around z at zt2
198  minuit->mnparm(47, "beamRot2Beam5", _vstart[47], _vstep[47], 0, 0, _ierflg);
199 
200  // rotation around z at zt1
201  minuit->mnparm(48, "beamRot1Beam6", _vstart[48], _vstep[48], 0, 0, _ierflg);
202  // rotation around z at zt2
203  minuit->mnparm(49, "beamRot2Beam6", _vstart[49], _vstep[49], 0, 0, _ierflg);
204 
205  // rotation around z at zt1
206  minuit->mnparm(50, "beamRot1Beam7", _vstart[50], _vstep[50], 0, 0, _ierflg);
207  // rotation around z at zt2
208  minuit->mnparm(51, "beamRot2Beam7", _vstart[51], _vstep[51], 0, 0, _ierflg);
209 
210  // we fix the respective outer disks 9 of each endcap
211  // as a reference system (pars 25,27,29 & 30,32,34)
212  // note: minuit numbering is fortran style...
213  arglist[0] = 26;
214  arglist[1] = 28;
215  arglist[2] = 30;
216  // minuit->mnexcm( "FIX", arglist ,3, _ierflg ); // TEC+
217  arglist[0] = 31;
218  arglist[1] = 33;
219  arglist[2] = 35;
220  // minuit->mnexcm( "FIX", arglist ,3, _ierflg ); // TEC-
221 
223  // DEBUG: FIX BEAM PARAMETERS /////////////////////////////////////////////////////////////////////
224  double parlist[16];
225  for (int par = 37; par <= 52; ++par)
226  parlist[par - 37] = par;
227  minuit->mnexcm("FIX", parlist, 16, _ierflg);
229 
231  // DEBUG: FIX ALGN PARAMETERS /////////////////////////////////////////////////////////////////////
232  // double parlist[36];
233  // for( int par = 1; par <= 36; ++par ) parlist[par-1] = par;
234  // minuit->mnexcm( "FIX", parlist ,36, _ierflg );
236 
237  // now ready for minimization step
238  arglist[0] = 10000;
239  arglist[1] = 0.1;
240  minuit->mnexcm("MIGRAD", arglist, 2, _ierflg); // minimizer
241  // minuit->mnexcm( "MINOS", arglist , 1, _ierflg ); // error recalculation
242 
243  // now fill the result vector.
244  // turned out that the parameter numbering is stupid, change this later..
245  LASBarrelAlignmentParameterSet theResult;
246  double par = 0., parError = 0.;
247 
248  // TEC+ rot
249  minuit->GetParameter(24, par, parError);
250  theResult.GetParameter(0, 0, 0).first = par;
251  theResult.GetParameter(0, 0, 0).second = parError;
252  minuit->GetParameter(25, par, parError);
253  theResult.GetParameter(0, 1, 0).first = par;
254  theResult.GetParameter(0, 1, 0).second = parError;
255  // TEC+ x
256  minuit->GetParameter(26, par, parError);
257  theResult.GetParameter(0, 0, 1).first = par;
258  theResult.GetParameter(0, 0, 1).second = parError;
259  minuit->GetParameter(27, par, parError);
260  theResult.GetParameter(0, 1, 1).first = par;
261  theResult.GetParameter(0, 1, 1).second = parError;
262  // TEC+ x
263  minuit->GetParameter(28, par, parError);
264  theResult.GetParameter(0, 0, 2).first = par;
265  theResult.GetParameter(0, 0, 2).second = parError;
266  minuit->GetParameter(29, par, parError);
267  theResult.GetParameter(0, 1, 2).first = par;
268  theResult.GetParameter(0, 1, 2).second = parError;
269 
270  // TEC- rot
271  minuit->GetParameter(30, par, parError);
272  theResult.GetParameter(1, 0, 0).first = par;
273  theResult.GetParameter(1, 0, 0).second = parError;
274  minuit->GetParameter(31, par, parError);
275  theResult.GetParameter(1, 1, 0).first = par;
276  theResult.GetParameter(1, 1, 0).second = parError;
277  // TEC- x
278  minuit->GetParameter(32, par, parError);
279  theResult.GetParameter(1, 0, 1).first = par;
280  theResult.GetParameter(1, 0, 1).second = parError;
281  minuit->GetParameter(33, par, parError);
282  theResult.GetParameter(1, 1, 1).first = par;
283  theResult.GetParameter(1, 1, 1).second = parError;
284  // TEC- x
285  minuit->GetParameter(34, par, parError);
286  theResult.GetParameter(1, 0, 2).first = par;
287  theResult.GetParameter(1, 0, 2).second = parError;
288  minuit->GetParameter(35, par, parError);
289  theResult.GetParameter(1, 1, 2).first = par;
290  theResult.GetParameter(1, 1, 2).second = parError;
291 
292  // TIB+ rot
293  minuit->GetParameter(0, par, parError);
294  theResult.GetParameter(2, 0, 0).first = par;
295  theResult.GetParameter(2, 0, 0).second = parError;
296  minuit->GetParameter(1, par, parError);
297  theResult.GetParameter(2, 1, 0).first = par;
298  theResult.GetParameter(2, 1, 0).second = parError;
299  // TIB+ x
300  minuit->GetParameter(2, par, parError);
301  theResult.GetParameter(2, 0, 1).first = par;
302  theResult.GetParameter(2, 0, 1).second = parError;
303  minuit->GetParameter(3, par, parError);
304  theResult.GetParameter(2, 1, 1).first = par;
305  theResult.GetParameter(2, 1, 1).second = parError;
306  // TIB+ x
307  minuit->GetParameter(4, par, parError);
308  theResult.GetParameter(2, 0, 2).first = par;
309  theResult.GetParameter(2, 0, 2).second = parError;
310  minuit->GetParameter(5, par, parError);
311  theResult.GetParameter(2, 1, 2).first = par;
312  theResult.GetParameter(2, 1, 2).second = parError;
313 
314  // TIB- rot
315  minuit->GetParameter(6, par, parError);
316  theResult.GetParameter(3, 0, 0).first = par;
317  theResult.GetParameter(3, 0, 0).second = parError;
318  minuit->GetParameter(7, par, parError);
319  theResult.GetParameter(3, 1, 0).first = par;
320  theResult.GetParameter(3, 1, 0).second = parError;
321  // TIB- x
322  minuit->GetParameter(8, par, parError);
323  theResult.GetParameter(3, 0, 1).first = par;
324  theResult.GetParameter(3, 0, 1).second = parError;
325  minuit->GetParameter(9, par, parError);
326  theResult.GetParameter(3, 1, 1).first = par;
327  theResult.GetParameter(3, 1, 1).second = parError;
328  // TIB- x
329  minuit->GetParameter(10, par, parError);
330  theResult.GetParameter(3, 0, 2).first = par;
331  theResult.GetParameter(3, 0, 2).second = parError;
332  minuit->GetParameter(11, par, parError);
333  theResult.GetParameter(3, 1, 2).first = par;
334  theResult.GetParameter(3, 1, 2).second = parError;
335 
336  // TOB+ rot
337  minuit->GetParameter(12, par, parError);
338  theResult.GetParameter(4, 0, 0).first = par;
339  theResult.GetParameter(4, 0, 0).second = parError;
340  minuit->GetParameter(13, par, parError);
341  theResult.GetParameter(4, 1, 0).first = par;
342  theResult.GetParameter(4, 1, 0).second = parError;
343  // TOB+ x
344  minuit->GetParameter(14, par, parError);
345  theResult.GetParameter(4, 0, 1).first = par;
346  theResult.GetParameter(4, 0, 1).second = parError;
347  minuit->GetParameter(15, par, parError);
348  theResult.GetParameter(4, 1, 1).first = par;
349  theResult.GetParameter(4, 1, 1).second = parError;
350  // TOB+ x
351  minuit->GetParameter(16, par, parError);
352  theResult.GetParameter(4, 0, 2).first = par;
353  theResult.GetParameter(4, 0, 2).second = parError;
354  minuit->GetParameter(17, par, parError);
355  theResult.GetParameter(4, 1, 2).first = par;
356  theResult.GetParameter(4, 1, 2).second = parError;
357 
358  // TOB- rot
359  minuit->GetParameter(18, par, parError);
360  theResult.GetParameter(5, 0, 0).first = par;
361  theResult.GetParameter(5, 0, 0).second = parError;
362  minuit->GetParameter(19, par, parError);
363  theResult.GetParameter(5, 1, 0).first = par;
364  theResult.GetParameter(5, 1, 0).second = parError;
365  // TOB- x
366  minuit->GetParameter(20, par, parError);
367  theResult.GetParameter(5, 0, 1).first = par;
368  theResult.GetParameter(5, 0, 1).second = parError;
369  minuit->GetParameter(21, par, parError);
370  theResult.GetParameter(5, 1, 1).first = par;
371  theResult.GetParameter(5, 1, 1).second = parError;
372  // TOB- x
373  minuit->GetParameter(22, par, parError);
374  theResult.GetParameter(5, 0, 2).first = par;
375  theResult.GetParameter(5, 0, 2).second = parError;
376  minuit->GetParameter(23, par, parError);
377  theResult.GetParameter(5, 1, 2).first = par;
378  theResult.GetParameter(5, 1, 2).second = parError;
379 
380  std::cout << " [LASBarrelAlgorithm::CalculateParameters] -- Done." << std::endl;
381 
382  return theResult;
383 }
384 
388 void fcn(int& npar, double* gin, double& f, double* par, int iflag) {
389  double chisquare = 0.;
390 
391  // the loop object and its variables
392  LASGlobalLoop moduleLoop;
393  int det, beam, pos, disk;
394 
396  // ADJUST THIS ALSO IN LASGeometryUpdater
398 
399  // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
400  // parameters are: det, side(0=+/1=-), z(lower/upper). TECs have no side (use side = 0)
401  std::vector<std::vector<std::vector<double> > > endFaceZPositions(
402  4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
403  endFaceZPositions.at(0).at(0).at(0) = 1322.5; // TEC+, *, disk1 ///
404  endFaceZPositions.at(0).at(0).at(1) = 2667.5; // TEC+, *, disk9 /// SIDE INFORMATION
405  endFaceZPositions.at(1).at(0).at(0) = -2667; // TEC-, *, disk1 /// MEANINGLESS FOR TEC -> USE .at(0)!
406  endFaceZPositions.at(1).at(0).at(1) = -1322.5; // TEC-, *, disk9 ///
407  endFaceZPositions.at(2).at(1).at(0) = -700.; // TIB, -, small z
408  endFaceZPositions.at(2).at(1).at(1) = -300.; // TIB, -, large z
409  endFaceZPositions.at(2).at(0).at(0) = 300.; // TIB, +, small z
410  endFaceZPositions.at(2).at(0).at(1) = 700.; // TIB, +, large z
411  endFaceZPositions.at(3).at(1).at(0) = -1090.; // TOB, -, small z
412  endFaceZPositions.at(3).at(1).at(1) = -300.; // TOB, -, large z
413  endFaceZPositions.at(3).at(0).at(0) = 300.; // TOB, +, small z
414  endFaceZPositions.at(3).at(0).at(1) = 1090.; // TOB, +, large z
415 
416  // the z positions of the virtual planes at which the beam parameters are measured
417  std::vector<double> disk9EndFaceZPositions(2, 0.);
418  disk9EndFaceZPositions.at(0) = -2667.5; // TEC- disk9
419  disk9EndFaceZPositions.at(1) = 2667.5; // TEC+ disk9
420 
421  // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
422  double detReducedZ[2] = {0., 0.};
423  // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
424  double beamReducedZ[2] = {0., 0.};
425 
426  // calculate residual for TIBTOB
427  det = 2;
428  beam = 0;
429  pos = 0;
430  do {
431  // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
432  const int theSide = pos < 3 ? 0 : 1;
433 
434  // this is the path the beam has to travel radially after being reflected
435  // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
436  const double radialOffset = det == 2 ? 50. : 36.;
437 
438  // reduced module's z position with respect to the subdetector endfaces
439  detReducedZ[0] =
440  aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ() - endFaceZPositions.at(det).at(theSide).at(0);
441  detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
442  detReducedZ[1] =
443  endFaceZPositions.at(det).at(theSide).at(1) - aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ();
444  detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
445 
446  // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
447  beamReducedZ[0] =
448  (aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset) - disk9EndFaceZPositions.at(0);
449  beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
450  beamReducedZ[1] =
451  disk9EndFaceZPositions.at(1) - (aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset);
452  beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
453 
454  // phi residual for this module as measured
455  const double measuredResidual = aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhi() - //&
456  aNominalCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhi();
457 
458  // shortcuts for speed
459  const double currentPhi = aNominalCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhi();
460  const double currentR = aNominalCoordinates->GetTIBTOBEntry(det, beam, pos).GetR();
461 
462  // phi residual for this module calculated from the parameters and nominal coordinates:
463  // this is the sum over the contributions from all parameters
464  double calculatedResidual = 0.;
465 
466  // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
467  // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
468 
469  // unfortunately, minuit keeps parameters in a 1-dim array,
470  // so we need to address the correct par[] for the 4 cases TIB+/TIB-/TOB+/TOB-
471  int indexBase = 0;
472  if (det == 2) { // TIB
473  if (theSide == 0)
474  indexBase = 0; // TIB+
475  if (theSide == 1)
476  indexBase = 6; // TIB-
477  }
478  if (det == 3) { // TOB
479  if (theSide == 0)
480  indexBase = 12; // TOB+
481  if (theSide == 1)
482  indexBase = 18; // TOB-
483  }
484 
485  // par[0] ("subRot1"): rotation around z of first end face
486  calculatedResidual += detReducedZ[1] * par[indexBase + 0];
487 
488  // par[1] ("subRot2"): rotation around z of second end face
489  calculatedResidual += detReducedZ[0] * par[indexBase + 1];
490 
491  // par[2] ("subTransX1"): translation along x of first end face
492  calculatedResidual += detReducedZ[1] * sin(currentPhi) / currentR * par[indexBase + 2];
493 
494  // par[3] ("subTransX2"): translation along x of second end face
495  calculatedResidual += detReducedZ[0] * sin(currentPhi) / currentR * par[indexBase + 3];
496 
497  // par[4] ("subTransY1"): translation along y of first end face
498  calculatedResidual += -1. * detReducedZ[1] * cos(currentPhi) / currentR * par[indexBase + 4];
499 
500  // par[5] ("subTransY2"): translation along y of second end face
501  calculatedResidual += -1. * detReducedZ[0] * cos(currentPhi) / currentR * par[indexBase + 5];
502 
503  // now come the 8*2 beam parameters, calculate the respective parameter index base first (-> which beam)
504  indexBase = 36 + beam * 2;
505 
506  // (there's no TIB/TOB/+/- distinction here for the beams)
507 
508  // ("beamRot1"): rotation around z at zt1
509  calculatedResidual += beamReducedZ[1] * par[indexBase];
510 
511  // ("beamRot2"): rotation around z at zt2
512  calculatedResidual += beamReducedZ[0] * par[indexBase + 1];
513 
514  // now calculate the chisquare
515  chisquare += pow(measuredResidual - calculatedResidual, 2) /
516  pow(aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhiError(), 2);
517 
518  } while (moduleLoop.TIBTOBLoop(det, beam, pos));
519 
520  // calculate residual for TEC AT
521  det = 0;
522  beam = 0;
523  disk = 0;
524  do {
525  // define the side: TECs sides already disentangled by the "det" index, so fix this to zero
526  const int theSide = 0;
527 
528  // reduced module's z position with respect to the subdetector endfaces
529  detReducedZ[0] =
530  aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ() - endFaceZPositions.at(det).at(theSide).at(0);
531  detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
532  detReducedZ[1] =
533  endFaceZPositions.at(det).at(theSide).at(1) - aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ();
534  detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
535 
536  // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
537  beamReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ() - disk9EndFaceZPositions.at(0);
538  beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
539  beamReducedZ[1] = disk9EndFaceZPositions.at(1) - aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ();
540  beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
541 
542  // phi residual for this module as measured
543  const double measuredResidual = aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhi() - //&
544  aNominalCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhi();
545 
546  // shortcuts for speed
547  const double currentPhi = aNominalCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhi();
548  const double currentR = aNominalCoordinates->GetTEC2TECEntry(det, beam, disk).GetR();
549 
550  // phi residual for this module calculated from the parameters and nominal coordinates:
551  // this is the sum over the contributions from all parameters
552  double calculatedResidual = 0.;
553 
554  // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
555  // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
556 
557  // there's also a distinction between TEC+/- parameters in situ (det==0 ? <TEC+> : <TEC->)
558 
559  // par[0] ("subRot1"): rotation around z of first end face
560  calculatedResidual += detReducedZ[1] * (det == 0 ? par[24] : par[30]);
561 
562  // par[1] ("subRot2"): rotation around z of second end face
563  calculatedResidual += detReducedZ[0] * (det == 0 ? par[25] : par[31]);
564 
565  // par[2] ("subTransX1"): translation along x of first end face
566  calculatedResidual += detReducedZ[1] * sin(currentPhi) * (det == 0 ? par[26] : par[32]) / currentR;
567 
568  // par[3] ("subTransX2"): translation along x of second end face
569  calculatedResidual += detReducedZ[0] * sin(currentPhi) * (det == 0 ? par[27] : par[33]) / currentR;
570 
571  // par[4] ("subTransY1"): translation along y of first end face
572  calculatedResidual += -1. * detReducedZ[1] * cos(currentPhi) * (det == 0 ? par[28] : par[34]) / currentR;
573 
574  // par[5] ("subTransY2"): translation along y of second end face
575  calculatedResidual += -1. * detReducedZ[0] * cos(currentPhi) * (det == 0 ? par[29] : par[35]) / currentR;
576 
577  // now come the 8*2 beam parameters; calculate the respective parameter index base first (-> which beam)
578  const unsigned int indexBase = 36 + beam * 2;
579 
580  // there's no TEC+/- distinction here
581 
582  // par[6] ("beamRot1"): rotation around z at zt1
583  calculatedResidual += beamReducedZ[1] * par[indexBase];
584 
585  // par[7] ("beamRot2"): rotation around z at zt2
586  calculatedResidual += beamReducedZ[0] * par[indexBase + 1];
587 
588  // now calculate the chisquare
589  // TODO: check for phi != 0 !!!
590  chisquare += pow(measuredResidual - calculatedResidual, 2) /
591  pow(aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhiError(), 2);
592 
593  } while (moduleLoop.TEC2TECLoop(det, beam, disk));
594 
595  // return the chisquare by ref
596  f = chisquare;
597 }
598 
603 void LASBarrelAlgorithm::Dump(void) {
604  if (!minuit) {
605  std::cerr << " [LASBarrelAlgorithm::Dump] ** WARNING: minimizer object uninitialized." << std::endl;
606  return;
607  }
608 
609  std::cout << std::endl << " [LASBarrelAlgorithm::Dump] -- Parameter dump: " << std::endl;
610 
611  const int subdetParMap[6] = {24, 30, 0, 6, 12, 18}; // map to one-dim array
612  const std::string subdetNames[6] = {" TEC+ ", " TEC- ", " TIB+ ", " TIB- ", " TOB+ ", " TOB- "};
613  double value, error;
614 
615  std::cout << " Detector parameters: " << std::endl;
616  std::cout << " -------------------" << std::endl;
617  std::cout << " Values: PHI1 X1 Y1 PHI2 X2 Y2 " << std::endl;
618  for (int subdet = 0; subdet < 6; ++subdet) {
619  std::cout << subdetNames[subdet];
620  for (int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2) {
621  minuit->GetParameter(par, value, error);
622  std::cout << std::setw(12) << std::setprecision(6) << std::fixed << value;
623  }
624  for (int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2) {
625  minuit->GetParameter(par, value, error);
626  std::cout << std::setw(12) << std::setprecision(6) << std::fixed << value;
627  }
628  std::cout << std::endl;
629  }
630 
631  std::cout << " Errors: PHI1 X1 Y1 PHI2 X2 Y2 " << std::endl;
632  for (int subdet = 0; subdet < 6; ++subdet) {
633  std::cout << subdetNames[subdet];
634  for (int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2) {
635  minuit->GetParameter(par, value, error);
636  std::cout << std::setw(12) << std::setprecision(6) << std::fixed << error;
637  }
638  for (int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2) {
639  minuit->GetParameter(par, value, error);
640  std::cout << std::setw(12) << std::setprecision(6) << std::fixed << error;
641  }
642  std::cout << std::endl;
643  }
644 
645  std::cout << std::endl;
646  std::cout << " Beam parameters: " << std::endl;
647  std::cout << " ---------------" << std::endl;
648  std::cout << " Values: PHI1 PHI2" << std::endl;
649  for (int beam = 0; beam < 8; ++beam) {
650  std::cout << " " << beam << " ";
651  for (int z = 0; z < 2; ++z) {
652  minuit->GetParameter(36 + 2 * beam + z, value, error);
653  std::cout << std::setw(12) << std::setprecision(6) << std::fixed << value;
654  }
655  std::cout << std::endl;
656  }
657 
658  std::cout << " Errors: PHI1 PHI2" << std::endl;
659  for (int beam = 0; beam < 8; ++beam) {
660  std::cout << " " << beam << " ";
661  for (int z = 0; z < 2; ++z) {
662  minuit->GetParameter(36 + 2 * beam + z, value, error);
663  std::cout << std::setw(12) << std::setprecision(6) << std::fixed << error;
664  }
665  std::cout << std::endl;
666  }
667 
668  // det parameters once again without leading column (for easy read-in), into a file
669  std::ofstream file("/afs/cern.ch/user/o/olzem/public/parameters_det.txt");
670  for (int subdet = 0; subdet < 6; ++subdet) {
671  for (int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2) {
672  minuit->GetParameter(par, value, error);
673  file << std::setw(12) << std::setprecision(6) << std::fixed << value;
674  }
675  for (int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2) {
676  minuit->GetParameter(par, value, error);
677  file << std::setw(12) << std::setprecision(6) << std::fixed << value;
678  }
679  file << std::endl;
680  }
681  file.close();
682 
683  // same for beam parameters
684  file.open("/afs/cern.ch/user/o/olzem/public/parameters_beam.txt");
685  for (int beam = 0; beam < 8; ++beam) {
686  for (int z = 0; z < 2; ++z) {
687  minuit->GetParameter(36 + 2 * beam + z, value, error);
688  file << std::setw(12) << std::setprecision(6) << std::fixed << value;
689  }
690  file << std::endl;
691  }
692  file.close();
693 
694  std::cout << " [LASBarrelAlgorithm::Dump] -- End parameter dump." << std::endl;
695  std::cout << std::endl;
696 }
697 
707 void LASBarrelAlgorithm::ReadMisalignmentFromFile(const char* filename,
708  LASGlobalData<LASCoordinateSet>& measuredCoordinates,
709  LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
710  std::ifstream file(filename);
711  if (file.bad()) {
712  std::cerr << " [LASBarrelAlgorithm::ReadMisalignmentFromFile] ** ERROR: cannot open file \"" << filename << "\"."
713  << std::endl;
714  return;
715  }
716 
717  std::cerr
718  << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
719  << std::endl;
720  std::cerr
721  << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
722  << std::endl;
723  std::cerr
724  << " [LASBarrelAlgorithm::ReadMisalignmentFromFile] ** WARNING: you are reading a fake measurement from a file!"
725  << std::endl;
726  std::cerr
727  << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
728  << std::endl;
729  std::cerr
730  << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
731  << std::endl;
732 
733  // the measured coordinates will finally be overwritten;
734  // first, set them to the nominal values
735  measuredCoordinates = nominalCoordinates;
736 
737  // and put large errors on all values;
738  {
739  LASGlobalLoop moduleLoop;
740  int det, ring, beam, disk, pos;
741 
742  det = 0;
743  ring = 0;
744  beam = 0;
745  disk = 0;
746  do {
747  measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhiError(1000.);
748  } while (moduleLoop.TECLoop(det, ring, beam, disk));
749 
750  det = 2;
751  beam = 0;
752  pos = 0;
753  do {
754  measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhiError(1000.);
755  } while (moduleLoop.TIBTOBLoop(det, beam, pos));
756 
757  det = 0;
758  beam = 0;
759  disk = 0;
760  do {
761  measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhiError(1000.);
762  } while (moduleLoop.TEC2TECLoop(det, beam, disk));
763  }
764 
765  // buffers for read-in
766  int det, beam, z, ring;
767  double phi, phiError;
768 
769  while (!file.eof()) {
770  file >> det;
771  if (file.eof())
772  break; // do not read the last line twice, do not fill trash if file empty
773 
774  file >> beam;
775  file >> z;
776  file >> ring;
777  file >> phi;
778  file >> phiError;
779 
780  if (det > 1) { // TIB/TOB
781  measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhi(phi);
782  measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhiError(phiError);
783  } else { // TEC or TEC(at)
784  if (ring > -1) { // TEC
785  measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhi(phi);
786  measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhiError(phiError);
787  } else { // TEC(at)
788  measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhi(phi);
789  measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhiError(phiError);
790  }
791  }
792  }
793 
794  file.close();
795 }
796 
805 void LASBarrelAlgorithm::ReadStartParametersFromFile(const char* filename, float values[52]) {
806  std::ifstream file(filename);
807  if (file.bad()) {
808  std::cerr << " [LASBarrelAlgorithm::ReadStartParametersFromFile] ** ERROR: cannot open file \"" << filename << "\"."
809  << std::endl;
810  return;
811  }
812 
813  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
814  "@@@@@@@@@@@@"
815  << std::endl;
816  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
817  "@@@@@@@@@@@@"
818  << std::endl;
819  std::cerr << " [LASBarrelAlgorithm::ReadStartParametersFrom File] ** WARNING: you are reading parameter start values "
820  "from a file!"
821  << std::endl;
822  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
823  "@@@@@@@@@@@@"
824  << std::endl;
825  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
826  "@@@@@@@@@@@@"
827  << std::endl;
828 
829  // map to the minuit par array
830  const int subdetParMap[6] = {24, 30, 0, 6, 12, 18};
831 
832  for (int det = 0; det < 6; ++det) {
833  file >> values[subdetParMap[det]]; // phi1
834  file >> values[subdetParMap[det] + 2]; // x1
835  file >> values[subdetParMap[det] + 4]; // y1
836  file >> values[subdetParMap[det] + 1]; // phi2
837  file >> values[subdetParMap[det] + 3]; // x2
838  file >> values[subdetParMap[det] + 5]; // y2
839  }
840 }
LASCoordinateSet::SetPhi
void SetPhi(double aPhi)
Definition: LASCoordinateSet.h:27
LASBarrelAlignmentParameterSet::GetParameter
std::pair< double, double > & GetParameter(int aSubdetector, int aDisk, int aParameter)
Definition: LASBarrelAlignmentParameterSet.cc:37
LASCoordinateSet::SetPhiError
void SetPhiError(double aPhiError)
Definition: LASCoordinateSet.h:32
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
aMeasuredCoordinates
LASGlobalData< LASCoordinateSet > * aMeasuredCoordinates
Definition: LASBarrelAlgorithm.cc:7
LASBarrelAlgorithm::ReadMisalignmentFromFile
void ReadMisalignmentFromFile(const char *, LASGlobalData< LASCoordinateSet > &, LASGlobalData< LASCoordinateSet > &)
Definition: LASBarrelAlgorithm.cc:707
LASCoordinateSet::GetPhi
double GetPhi(void) const
Definition: LASCoordinateSet.h:18
aNominalCoordinates
LASGlobalData< LASCoordinateSet > * aNominalCoordinates
Definition: LASBarrelAlgorithm.cc:8
EcalCondDBWriter_cfi.beam
string beam
Definition: EcalCondDBWriter_cfi.py:45
LASBarrelAlgorithm::CalculateParameters
LASBarrelAlignmentParameterSet CalculateParameters(LASGlobalData< LASCoordinateSet > &, LASGlobalData< LASCoordinateSet > &)
Definition: LASBarrelAlgorithm.cc:20
LASGlobalLoop::TEC2TECLoop
bool TEC2TECLoop(int &, int &, int &) const
Definition: LASGlobalLoop.cc:74
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
LASBarrelAlgorithm::ReadStartParametersFromFile
void ReadStartParametersFromFile(const char *, float[52])
Definition: LASBarrelAlgorithm.cc:805
LASGlobalData::GetTIBTOBEntry
T & GetTIBTOBEntry(int subdetector, int beam, int tibTobPosition)
Definition: LASGlobalData.h:107
LASGlobalData::GetTEC2TECEntry
T & GetTEC2TECEntry(int subdetector, int beam, int tecDisk)
Definition: LASGlobalData.h:129
fcn
void fcn(int &, double *, double &, double *, int)
Definition: LASBarrelAlgorithm.cc:388
LASGlobalData::GetTECEntry
T & GetTECEntry(int subdetector, int tecRing, int beam, int tecDisk)
Definition: LASGlobalData.h:84
LASCoordinateSet::GetZ
double GetZ(void) const
Definition: LASCoordinateSet.h:22
lut2db_cfg.filename
tuple filename
Definition: lut2db_cfg.py:20
LASGlobalLoop::TECLoop
bool TECLoop(int &, int &, int &, int &) const
Definition: LASGlobalLoop.cc:13
LASGlobalLoop::TIBTOBLoop
bool TIBTOBLoop(int &, int &, int &) const
Definition: LASGlobalLoop.cc:46
LASCoordinateSet::GetPhiError
double GetPhiError(void) const
Definition: LASCoordinateSet.h:19
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0.cerr
tuple cerr
Definition: EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0.py:8
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
LASCoordinateSet::GetR
double GetR(void) const
Definition: LASCoordinateSet.h:20
subdetNames
static const char * subdetNames[]
Definition: parseHcalDetId.cc:11
Generated for CMSSW Reference Manual by doxygen 1.8.5