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