CMS 3D CMS Logo
LASAlignmentTubeAlgorithm.cc
Go to the documentation of this file.
1 
2 #include "Alignment/LaserAlignment/interface/LASAlignmentTubeAlgorithm.h"
3 
7 LASAlignmentTubeAlgorithm::LASAlignmentTubeAlgorithm() {}
8 
12 LASBarrelAlignmentParameterSet LASAlignmentTubeAlgorithm::CalculateParameters(
13  LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
14  std::cout << " [LASAlignmentTubeAlgorithm::CalculateParameters] -- Starting." << std::endl;
15 
16  // for debugging only
17  //######################################################################################
18  //ReadMisalignmentFromFile( "misalign-var.txt", measuredCoordinates, nominalCoordinates );
19  //######################################################################################
20 
21  // loop object
22  LASGlobalLoop globalLoop;
23  int det, beam, disk, pos;
24 
25  // phi positions of the AT beams in rad
26  const double phiPositions[8] = {0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};
27  std::vector<double> beamPhiPositions(8, 0.);
28  for (beam = 0; beam < 8; ++beam)
29  beamPhiPositions.at(beam) = phiPositions[beam];
30 
31  // the radii of the alignment tube beams for each halfbarrel.
32  // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
33  // in TIB/TOB modules these radii differ from the beam radius..
34  // ..due to the radial offsets (after the semitransparent mirrors)
35  const double radii[6] = {564., 564., 514., 514., 600., 600.};
36  std::vector<double> beamRadii(6, 0.);
37  for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel)
38  beamRadii.at(aHalfbarrel) = radii[aHalfbarrel];
39 
40  // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
41  // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
42  std::vector<std::vector<std::vector<double> > > endFaceZPositions(
43  4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
44  endFaceZPositions.at(0).at(0).at(0) = 1322.5; // TEC+, *, disk1 ///
45  endFaceZPositions.at(0).at(0).at(1) = 2667.5; // TEC+, *, disk9 /// SIDE INFORMATION
46  endFaceZPositions.at(1).at(0).at(0) = -2667.5; // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
47  endFaceZPositions.at(1).at(0).at(1) = -1322.5; // TEC-, *, disk1 ///
48  endFaceZPositions.at(2).at(1).at(0) = -700.; // TIB, -, outer
49  endFaceZPositions.at(2).at(1).at(1) = -300.; // TIB, -, inner
50  endFaceZPositions.at(2).at(0).at(0) = 300.; // TIB, +, inner
51  endFaceZPositions.at(2).at(0).at(1) = 700.; // TIB, +, outer
52  endFaceZPositions.at(3).at(1).at(0) = -1090.; // TOB, -, outer
53  endFaceZPositions.at(3).at(1).at(1) = -300.; // TOB, -, inner
54  endFaceZPositions.at(3).at(0).at(0) = 300.; // TOB, +, inner
55  endFaceZPositions.at(3).at(0).at(1) = 1090.; // TOB, +, outer
56 
57  // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
58  double detReducedZ[2] = {0., 0.};
59  // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
60  double beamReducedZ[2] = {0., 0.};
61 
62  // the z positions of the virtual planes at which the beam parameters are measured
63  std::vector<double> disk9EndFaceZPositions(2, 0.);
64  disk9EndFaceZPositions.at(0) = -2667.5; // TEC- disk9
65  disk9EndFaceZPositions.at(1) = 2667.5; // TEC+ disk9
66 
67  // define sums over measured values to "simplify" the beam parameter formulas
68 
69  // all these have 6 entries, one for each halfbarrel (TEC+,TEC-,TIB+,TIB-,TOB+,TOB-)
70  std::vector<double> sumOverPhiZPrime(6, 0.);
71  std::vector<double> sumOverPhiZPrimePrime(6, 0.);
72  std::vector<double> sumOverPhiZPrimeSinTheta(6, 0.);
73  std::vector<double> sumOverPhiZPrimePrimeSinTheta(6, 0.);
74  std::vector<double> sumOverPhiZPrimeCosTheta(6, 0.);
75  std::vector<double> sumOverPhiZPrimePrimeCosTheta(6, 0.);
76 
77  // these have 8 entries, one for each beam
78  std::vector<double> sumOverPhiZTPrime(8, 0.);
79  std::vector<double> sumOverPhiZTPrimePrime(8, 0.);
80 
81  // define sums over nominal values
82 
83  // all these have 6 entries, one for each halfbarrel (TEC+,TEC-,TIB+,TIB-,TOB+,TOB-)
84  std::vector<double> sumOverZPrimeSquared(6, 0.);
85  std::vector<double> sumOverZPrimePrimeSquared(6, 0.);
86  std::vector<double> sumOverZPrimeZPrimePrime(6, 0.);
87  std::vector<double> sumOverZPrimeZTPrime(6, 0.);
88  std::vector<double> sumOverZPrimeZTPrimePrime(6, 0.);
89  std::vector<double> sumOverZPrimePrimeZTPrime(6, 0.);
90  std::vector<double> sumOverZPrimePrimeZTPrimePrime(6, 0.);
91 
92  // all these are scalars
93  double sumOverZTPrimeSquared = 0.;
94  double sumOverZTPrimePrimeSquared = 0.;
95  double sumOverZTPrimeZTPrimePrime = 0.;
96 
97  // now calculate them for TIBTOB
98  det = 2;
99  beam = 0;
100  pos = 0;
101  do {
102  // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
103  const int theSide = pos < 3 ? 0 : 1;
104 
105  // define the halfbarrel number from det/side
106  const int halfbarrel = det == 2 ? det + theSide : det + 1 + theSide; // TIB:TOB
107 
108  // this is the path the beam has to travel radially after being reflected
109  // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
110  const double radialOffset = det == 2 ? 50. : 36.;
111 
112  // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
113  detReducedZ[0] = measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() -
114  endFaceZPositions.at(det).at(theSide).at(0); // = zPrime
115  detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
116  detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
117  measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ(); // = zPrimePrime
118  detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
119 
120  // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
121  beamReducedZ[0] = (measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset) -
122  disk9EndFaceZPositions.at(0); // = ZTPrime
123  beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
124  beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
125  (measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset); // ZTPrimePrime
126  beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
127 
128  // residual in phi (in the formulas this corresponds to y_ik/R)
129  const double phiResidual = measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetPhi() -
130  nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetPhi();
131 
132  // sums over measured values
133  sumOverPhiZPrime.at(halfbarrel) += phiResidual * detReducedZ[0];
134  sumOverPhiZPrimePrime.at(halfbarrel) += phiResidual * detReducedZ[1];
135  sumOverPhiZPrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * sin(beamPhiPositions.at(beam));
136  sumOverPhiZPrimePrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * sin(beamPhiPositions.at(beam));
137  sumOverPhiZPrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * cos(beamPhiPositions.at(beam));
138  sumOverPhiZPrimePrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * cos(beamPhiPositions.at(beam));
139 
140  sumOverPhiZTPrime.at(beam) += phiResidual * beamReducedZ[0]; // note the index change here..
141  sumOverPhiZTPrimePrime.at(beam) += phiResidual * beamReducedZ[1];
142 
143  // sums over nominal values
144  sumOverZPrimeSquared.at(halfbarrel) += pow(detReducedZ[0], 2) / 8.; // these are defined beam-wise, so: / 8.
145  sumOverZPrimePrimeSquared.at(halfbarrel) += pow(detReducedZ[1], 2) / 8.;
146  sumOverZPrimeZPrimePrime.at(halfbarrel) += detReducedZ[0] * detReducedZ[1] / 8.;
147  sumOverZPrimeZTPrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[0] / 8.;
148  sumOverZPrimeZTPrimePrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[1] / 8.;
149  sumOverZPrimePrimeZTPrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[0] / 8.;
150  sumOverZPrimePrimeZTPrimePrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[1] / 8.;
151 
152  sumOverZTPrimeSquared += pow(beamReducedZ[0], 2) / 8.;
153  sumOverZTPrimePrimeSquared += pow(beamReducedZ[1], 2) / 8.;
154  sumOverZTPrimeZTPrimePrime += beamReducedZ[0] * beamReducedZ[1] / 8.;
155 
156  } while (globalLoop.TIBTOBLoop(det, beam, pos));
157 
158  // now for TEC2TEC
159  det = 0;
160  beam = 0;
161  disk = 0;
162  do {
163  // for the tec, the halfbarrel numbers are equal to the det numbers...
164  const int halfbarrel = det;
165 
166  // ...so there's no side distinction for the TEC
167  const int theSide = 0;
168 
169  // also, there's no radial offset for the TEC
170  const double radialOffset = 0.;
171 
172  // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
173  detReducedZ[0] = measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() -
174  endFaceZPositions.at(det).at(theSide).at(0); // = zPrime
175  detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
176  detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
177  measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ(); // = zPrimePrime
178  detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
179 
180  // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
181  beamReducedZ[0] = (measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset) -
182  disk9EndFaceZPositions.at(0); // = ZTPrime
183  beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
184  beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
185  (measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset); // ZTPrimePrime
186  beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
187 
188  // residual in phi (in the formulas this corresponds to y_ik/R)
189  const double phiResidual = measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetPhi() -
190  nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetPhi();
191 
192  // sums over measured values
193  sumOverPhiZPrime.at(halfbarrel) += phiResidual * detReducedZ[0];
194  sumOverPhiZPrimePrime.at(halfbarrel) += phiResidual * detReducedZ[1];
195  sumOverPhiZPrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * sin(beamPhiPositions.at(beam));
196  sumOverPhiZPrimePrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * sin(beamPhiPositions.at(beam));
197  sumOverPhiZPrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * cos(beamPhiPositions.at(beam));
198  sumOverPhiZPrimePrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * cos(beamPhiPositions.at(beam));
199 
200  sumOverPhiZTPrime.at(beam) += phiResidual * beamReducedZ[0]; // note the index change here..
201  sumOverPhiZTPrimePrime.at(beam) += phiResidual * beamReducedZ[1];
202 
203  // sums over nominal values
204  sumOverZPrimeSquared.at(halfbarrel) += pow(detReducedZ[0], 2) / 8.; // these are defined beam-wise, so: / 8.
205  sumOverZPrimePrimeSquared.at(halfbarrel) += pow(detReducedZ[1], 2) / 8.;
206  sumOverZPrimeZPrimePrime.at(halfbarrel) += detReducedZ[0] * detReducedZ[1] / 8.;
207  sumOverZPrimeZTPrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[0] / 8.;
208  sumOverZPrimeZTPrimePrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[1] / 8.;
209  sumOverZPrimePrimeZTPrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[0] / 8.;
210  sumOverZPrimePrimeZTPrimePrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[1] / 8.;
211 
212  sumOverZTPrimeSquared += pow(beamReducedZ[0], 2) / 8.;
213  sumOverZTPrimePrimeSquared += pow(beamReducedZ[1], 2) / 8.;
214  sumOverZTPrimeZTPrimePrime += beamReducedZ[0] * beamReducedZ[1] / 8.;
215 
216  } while (globalLoop.TEC2TECLoop(det, beam, disk));
217 
218  // more "simplification" terms...
219  // these here are functions of theta and can be calculated directly
220  double sumOverSinTheta = 0.;
221  double sumOverCosTheta = 0.;
222  double sumOverSinThetaSquared = 0.;
223  double sumOverCosThetaSquared = 0.;
224  double sumOverCosThetaSinTheta = 0.;
225  double mixedTrigonometricTerm = 0.;
226 
227  for (beam = 0; beam < 8; ++beam) {
228  sumOverSinTheta += sin(beamPhiPositions.at(beam));
229  sumOverCosTheta += cos(beamPhiPositions.at(beam));
230  sumOverSinThetaSquared += pow(sin(beamPhiPositions.at(beam)), 2);
231  sumOverCosThetaSquared += pow(cos(beamPhiPositions.at(beam)), 2);
232  sumOverCosThetaSinTheta += cos(beamPhiPositions.at(beam)) * sin(beamPhiPositions.at(beam));
233  }
234 
235  mixedTrigonometricTerm = 8. * (sumOverCosThetaSquared * sumOverSinThetaSquared - pow(sumOverCosThetaSinTheta, 2)) -
236  pow(sumOverCosTheta, 2) * sumOverSinThetaSquared -
237  pow(sumOverSinTheta, 2) * sumOverCosThetaSquared +
238  2. * sumOverCosTheta * sumOverSinTheta * sumOverCosThetaSinTheta;
239 
240  // even more shortcuts before we can calculate the parameters
241  double beamDenominator =
242  (pow(sumOverZTPrimeZTPrimePrime, 2) - sumOverZTPrimeSquared * sumOverZTPrimePrimeSquared) * beamRadii.at(0);
243  std::vector<double> alignmentDenominator(6, 0.);
244  std::vector<double> termA(6, 0.), termB(6, 0.), termC(6, 0.), termD(6, 0.);
245  for (unsigned int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel) {
246  alignmentDenominator.at(aHalfbarrel) =
247  (pow(sumOverZPrimeZPrimePrime.at(aHalfbarrel), 2) -
248  sumOverZPrimeSquared.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel)) *
249  mixedTrigonometricTerm;
250  termA.at(aHalfbarrel) = sumOverZPrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
251  sumOverZPrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeSquared;
252  termB.at(aHalfbarrel) = sumOverZPrimePrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
253  sumOverZPrimePrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeSquared;
254  termC.at(aHalfbarrel) = sumOverZPrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
255  sumOverZPrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimePrimeSquared;
256  termD.at(aHalfbarrel) = sumOverZPrimePrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
257  sumOverZPrimePrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimePrimeSquared;
258  }
259 
260  // have eight alignment tube beams..
261  const int numberOfBeams = 8;
262 
263  // that's all for preparation, now it gets ugly:
264  // calculate the alignment parameters
265  LASBarrelAlignmentParameterSet theResult;
266 
267  // can do this in one go for all halfbarrels
268  for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel) {
269  // no errors yet
270 
271  // rotation angles of the lower z endfaces (in rad)
272  theResult.GetParameter(aHalfbarrel, 0, 0).first =
273  (sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
274  sumOverSinThetaSquared -
275  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
276  sumOverSinThetaSquared -
277  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
278  sumOverSinThetaSquared +
279  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
280  sumOverSinThetaSquared +
281  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
282  sumOverSinTheta -
283  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) *
284  sumOverCosThetaSinTheta * sumOverSinTheta -
285  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
286  sumOverSinTheta +
287  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
288  sumOverSinTheta -
289  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * pow(sumOverCosThetaSinTheta, 2) +
290  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) *
291  pow(sumOverCosThetaSinTheta, 2) +
292  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
293  sumOverCosThetaSinTheta -
294  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
295  sumOverCosThetaSinTheta) /
296  alignmentDenominator.at(aHalfbarrel);
297 
298  // rotation angles of the upper z endfaces (in rad)
299  theResult.GetParameter(aHalfbarrel, 1, 0).first =
300  (sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
301  sumOverSinThetaSquared -
302  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
303  sumOverSinThetaSquared -
304  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
305  sumOverSinThetaSquared +
306  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
307  sumOverSinThetaSquared +
308  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) *
309  sumOverCosThetaSinTheta * sumOverSinTheta -
310  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) *
311  sumOverCosThetaSinTheta * sumOverSinTheta -
312  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) *
313  sumOverCosThetaSquared * sumOverSinTheta +
314  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
315  sumOverSinTheta -
316  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
317  sumOverCosThetaSinTheta +
318  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSinTheta *
319  sumOverCosThetaSinTheta +
320  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
321  sumOverCosThetaSinTheta -
322  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
323  sumOverCosThetaSinTheta) /
324  alignmentDenominator.at(aHalfbarrel);
325 
326  // x deviations of the lower z endfaces (in mm)
327  theResult.GetParameter(aHalfbarrel, 0, 1).first =
328  -1. *
329  (sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
330  sumOverSinTheta -
331  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
332  sumOverSinTheta -
333  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
334  sumOverSinTheta +
335  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
336  sumOverSinTheta -
337  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
338  sumOverCosThetaSinTheta +
339  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
340  sumOverCosThetaSinTheta +
341  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
342  sumOverCosThetaSinTheta -
343  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
344  sumOverCosThetaSinTheta -
345  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
346  sumOverCosThetaSquared +
347  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
348  sumOverCosThetaSquared +
349  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
350  sumOverCosTheta -
351  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
352  sumOverCosTheta) /
353  alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
354 
355  // x deviations of the upper z endfaces (in mm)
356  theResult.GetParameter(aHalfbarrel, 1, 1).first =
357  -1. *
358  (sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
359  sumOverSinTheta -
360  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
361  sumOverSinTheta -
362  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
363  sumOverSinTheta +
364  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
365  sumOverSinTheta -
366  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
367  sumOverCosThetaSinTheta +
368  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
369  sumOverCosThetaSinTheta +
370  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
371  sumOverCosThetaSinTheta -
372  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
373  sumOverCosThetaSinTheta -
374  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
375  sumOverCosThetaSquared +
376  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
377  sumOverCosThetaSquared +
378  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
379  sumOverCosTheta -
380  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
381  sumOverCosTheta) /
382  alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
383 
384  // y deviations of the lower z endfaces (in mm)
385  theResult.GetParameter(aHalfbarrel, 0, 2).first =
386  (sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
387  sumOverSinThetaSquared -
388  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
389  sumOverSinThetaSquared -
390  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
391  sumOverSinThetaSquared +
392  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
393  sumOverSinThetaSquared +
394  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverSinTheta *
395  sumOverSinTheta -
396  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverSinTheta *
397  sumOverSinTheta -
398  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
399  sumOverSinTheta +
400  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSinTheta *
401  sumOverSinTheta -
402  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
403  sumOverSinTheta +
404  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
405  sumOverSinTheta +
406  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
407  sumOverCosThetaSinTheta -
408  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
409  sumOverCosThetaSinTheta) /
410  alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
411 
412  // y deviations of the upper z endfaces (in mm)
413  theResult.GetParameter(aHalfbarrel, 1, 2).first =
414  (sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
415  sumOverSinThetaSquared -
416  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
417  sumOverSinThetaSquared -
418  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
419  sumOverSinThetaSquared +
420  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
421  sumOverSinThetaSquared +
422  sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverSinTheta *
423  sumOverSinTheta -
424  sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverSinTheta *
425  sumOverSinTheta -
426  sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
427  sumOverSinTheta +
428  sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSinTheta *
429  sumOverSinTheta -
430  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
431  sumOverSinTheta +
432  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
433  sumOverSinTheta +
434  sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
435  sumOverCosThetaSinTheta -
436  sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
437  sumOverCosThetaSinTheta) /
438  alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
439  }
440 
441  // another loop is needed here to calculate some terms for the beam parameters
442  double vsumA = 0., vsumB = 0., vsumC = 0., vsumD = 0., vsumE = 0., vsumF = 0.;
443  for (unsigned int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel) {
444  vsumA += theResult.GetParameter(aHalfbarrel, 1, 2).first * termA.at(aHalfbarrel) +
445  theResult.GetParameter(aHalfbarrel, 0, 2).first * termB.at(aHalfbarrel);
446  vsumB += theResult.GetParameter(aHalfbarrel, 1, 1).first * termA.at(aHalfbarrel) +
447  theResult.GetParameter(aHalfbarrel, 0, 1).first * termB.at(aHalfbarrel);
448  vsumC += beamRadii.at(aHalfbarrel) * (theResult.GetParameter(aHalfbarrel, 1, 0).first * termA.at(aHalfbarrel) +
449  theResult.GetParameter(aHalfbarrel, 0, 0).first * termB.at(aHalfbarrel));
450  vsumD += theResult.GetParameter(aHalfbarrel, 1, 2).first * termC.at(aHalfbarrel) +
451  theResult.GetParameter(aHalfbarrel, 0, 2).first * termD.at(aHalfbarrel);
452  vsumE += theResult.GetParameter(aHalfbarrel, 1, 1).first * termC.at(aHalfbarrel) +
453  theResult.GetParameter(aHalfbarrel, 0, 1).first * termD.at(aHalfbarrel);
454  vsumF += beamRadii.at(aHalfbarrel) * (theResult.GetParameter(aHalfbarrel, 1, 0).first * termC.at(aHalfbarrel) +
455  theResult.GetParameter(aHalfbarrel, 0, 0).first * termD.at(aHalfbarrel));
456  }
457 
458  // calculate the beam parameters
459  for (unsigned int beam = 0; beam < 8; ++beam) {
460  // parameter A, defined at lower z
461  theResult.GetBeamParameter(beam, 0).first =
462  (cos(beamPhiPositions.at(beam)) * vsumA - sin(beamPhiPositions.at(beam)) * vsumB - vsumC +
463  sumOverPhiZTPrime.at(beam) * sumOverZTPrimeZTPrimePrime -
464  sumOverPhiZTPrimePrime.at(beam) * sumOverZTPrimeSquared) /
465  beamDenominator;
466 
468  std::cout << "BBBBBBBB: " << cos(beamPhiPositions.at(beam)) * vsumA << " "
469  << -1. * sin(beamPhiPositions.at(beam)) * vsumB << " " << -1. * vsumC << " "
470  << sumOverPhiZTPrime.at(beam) * sumOverZTPrimeZTPrimePrime -
471  sumOverPhiZTPrimePrime.at(beam) * sumOverZTPrimeSquared
472  << " " << beamDenominator << std::endl;
474 
475  // parameter B, defined at upper z
476  theResult.GetBeamParameter(beam, 1).first =
477  (cos(beamPhiPositions.at(beam)) * vsumD - sin(beamPhiPositions.at(beam)) * vsumE - vsumF +
478  sumOverPhiZTPrimePrime.at(beam) * sumOverZTPrimeZTPrimePrime -
479  sumOverPhiZTPrime.at(beam) * sumOverZTPrimePrimeSquared) /
480  beamDenominator;
481  }
482 
483  return theResult;
484 }
485 
490 double LASAlignmentTubeAlgorithm::GetTIBTOBAlignmentParameterCorrection(
491  int det,
492  int beam,
493  int pos,
494  LASGlobalData<LASCoordinateSet>& nominalCoordinates,
495  LASBarrelAlignmentParameterSet& alignmentParameters) {
496  // INITIALIZATION;
497  // ALL THIS IS DUPLICATED FOR THE MOMENT, SHOULD FINALLY BE CALCULATED ONLY ONCE
498  // AND HARD CODED NUMBERS SHOULD CENTRALLY BE IMPORTED FROM src/LASConstants.h
499 
500  // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
501  // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
502  std::vector<std::vector<std::vector<double> > > endFaceZPositions(
503  4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
504  endFaceZPositions.at(0).at(0).at(0) = 1322.5; // TEC+, *, disk1 ///
505  endFaceZPositions.at(0).at(0).at(1) = 2667.5; // TEC+, *, disk9 /// SIDE INFORMATION
506  endFaceZPositions.at(1).at(0).at(0) = -2667.5; // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
507  endFaceZPositions.at(1).at(0).at(1) = -1322.5; // TEC-, *, disk1 ///
508  endFaceZPositions.at(2).at(1).at(0) = -700.; // TIB, -, outer
509  endFaceZPositions.at(2).at(1).at(1) = -300.; // TIB, -, inner
510  endFaceZPositions.at(2).at(0).at(0) = 300.; // TIB, +, inner
511  endFaceZPositions.at(2).at(0).at(1) = 700.; // TIB, +, outer
512  endFaceZPositions.at(3).at(1).at(0) = -1090.; // TOB, -, outer
513  endFaceZPositions.at(3).at(1).at(1) = -300.; // TOB, -, inner
514  endFaceZPositions.at(3).at(0).at(0) = 300.; // TOB, +, inner
515  endFaceZPositions.at(3).at(0).at(1) = 1090.; // TOB, +, outer
516 
517  // the z positions of the virtual planes at which the beam parameters are measured
518  std::vector<double> disk9EndFaceZPositions(2, 0.);
519  disk9EndFaceZPositions.at(0) = -2667.5; // TEC- disk9
520  disk9EndFaceZPositions.at(1) = 2667.5; // TEC+ disk9
521 
522  // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
523  const int theSide = pos < 3 ? 0 : 1;
524 
525  // define the halfbarrel number from det/side
526  const int halfbarrel = det == 2 ? det + theSide : det + 1 + theSide; // TIB:TOB
527 
528  // this is the path the beam has to travel radially after being reflected
529  // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
530  const double radialOffset = det == 2 ? 50. : 36.;
531 
532  // phi positions of the AT beams in rad
533  const double phiPositions[8] = {0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};
534  std::vector<double> beamPhiPositions(8, 0.);
535  for (unsigned int aBeam = 0; aBeam < 8; ++aBeam)
536  beamPhiPositions.at(aBeam) = phiPositions[aBeam];
537 
538  // the radii of the alignment tube beams for each halfbarrel.
539  // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
540  // in TIB/TOB modules these radii differ from the beam radius..
541  // ..due to the radial offsets (after the semitransparent mirrors)
542  const double radii[6] = {564., 564., 514., 514., 600., 600.};
543  std::vector<double> beamRadii(6, 0.);
544  for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel)
545  beamRadii.at(aHalfbarrel) = radii[aHalfbarrel];
546 
547  // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
548  double detReducedZ[2] = {0., 0.};
549  // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
550  double beamReducedZ[2] = {0., 0.};
551 
552  // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
553  detReducedZ[0] = nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() -
554  endFaceZPositions.at(det).at(theSide).at(0); // = zPrime
555  detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
556  detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
557  nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ(); // = zPrimePrime
558  detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
559 
560  // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
561  beamReducedZ[0] = (nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset) -
562  disk9EndFaceZPositions.at(0); // = ZTPrime
563  beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
564  beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
565  (nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset); // ZTPrimePrime
566  beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
567 
568  // the correction to phi from the endcap algorithm;
569  // it is defined such that the correction is to be subtracted ///////////////////////////////// ???
570  double phiCorrection = 0.;
571 
572  // contribution from phi rotation of first end face
573  phiCorrection += detReducedZ[1] * alignmentParameters.GetParameter(halfbarrel, 0, 0).first;
574 
575  // contribution from phi rotation of second end face
576  phiCorrection += detReducedZ[0] * alignmentParameters.GetParameter(halfbarrel, 1, 0).first;
577 
578  // contribution from translation along x of first endface
579  phiCorrection += detReducedZ[1] * sin(beamPhiPositions.at(beam)) *
580  alignmentParameters.GetParameter(halfbarrel, 0, 1).first / beamRadii.at(halfbarrel);
581 
582  // contribution from translation along x of second endface
583  phiCorrection += detReducedZ[0] * sin(beamPhiPositions.at(beam)) *
584  alignmentParameters.GetParameter(halfbarrel, 1, 1).first / beamRadii.at(halfbarrel);
585 
586  // contribution from translation along y of first endface
587  phiCorrection -= detReducedZ[1] * cos(beamPhiPositions.at(beam)) *
588  alignmentParameters.GetParameter(halfbarrel, 0, 2).first / beamRadii.at(halfbarrel);
589 
590  // contribution from translation along y of second endface
591  phiCorrection -= detReducedZ[0] * cos(beamPhiPositions.at(beam)) *
592  alignmentParameters.GetParameter(halfbarrel, 1, 2).first / beamRadii.at(halfbarrel);
593 
594  // contribution from beam parameters;
595  // originally, the contribution in meter is proportional to the radius of the beams: beamRadii.at( 0 )
596  // the additional factor: beamRadii.at( halfbarrel ) converts from meter to radian on the module
597  phiCorrection += beamReducedZ[1] * alignmentParameters.GetBeamParameter(beam, 0).first * beamRadii.at(0) /
598  beamRadii.at(halfbarrel);
599  phiCorrection += beamReducedZ[0] * alignmentParameters.GetBeamParameter(beam, 1).first * beamRadii.at(0) /
600  beamRadii.at(halfbarrel);
601 
602  return phiCorrection;
603 }
604 
609 double LASAlignmentTubeAlgorithm::GetTEC2TECAlignmentParameterCorrection(
610  int det,
611  int beam,
612  int disk,
613  LASGlobalData<LASCoordinateSet>& nominalCoordinates,
614  LASBarrelAlignmentParameterSet& alignmentParameters) {
615  // INITIALIZATION;
616  // ALL THIS IS DUPLICATED FOR THE MOMENT, SHOULD FINALLY BE CALCULATED ONLY ONCE
617  // AND HARD CODED NUMBERS SHOULD CENTRALLY BE IMPORTED FROM src/LASConstants.h
618 
619  // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
620  // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
621  std::vector<std::vector<std::vector<double> > > endFaceZPositions(
622  4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
623  endFaceZPositions.at(0).at(0).at(0) = 1322.5; // TEC+, *, disk1 ///
624  endFaceZPositions.at(0).at(0).at(1) = 2667.5; // TEC+, *, disk9 /// SIDE INFORMATION
625  endFaceZPositions.at(1).at(0).at(0) = -2667.5; // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
626  endFaceZPositions.at(1).at(0).at(1) = -1322.5; // TEC-, *, disk1 ///
627  endFaceZPositions.at(2).at(1).at(0) = -700.; // TIB, -, outer
628  endFaceZPositions.at(2).at(1).at(1) = -300.; // TIB, -, inner
629  endFaceZPositions.at(2).at(0).at(0) = 300.; // TIB, +, inner
630  endFaceZPositions.at(2).at(0).at(1) = 700.; // TIB, +, outer
631  endFaceZPositions.at(3).at(1).at(0) = -1090.; // TOB, -, outer
632  endFaceZPositions.at(3).at(1).at(1) = -300.; // TOB, -, inner
633  endFaceZPositions.at(3).at(0).at(0) = 300.; // TOB, +, inner
634  endFaceZPositions.at(3).at(0).at(1) = 1090.; // TOB, +, outer
635 
636  // the z positions of the virtual planes at which the beam parameters are measured
637  std::vector<double> disk9EndFaceZPositions(2, 0.);
638  disk9EndFaceZPositions.at(0) = -2667.5; // TEC- disk9
639  disk9EndFaceZPositions.at(1) = 2667.5; // TEC+ disk9
640 
641  // for the tec, the halfbarrel numbers are equal to the det numbers...
642  const int halfbarrel = det;
643 
644  // ...so there's no side distinction for the TEC
645  const int theSide = 0;
646 
647  // also, there's no radial offset for the TEC
648  const double radialOffset = 0.;
649 
650  // phi positions of the AT beams in rad
651  const double phiPositions[8] = {0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};
652  std::vector<double> beamPhiPositions(8, 0.);
653  for (unsigned int aBeam = 0; aBeam < 8; ++aBeam)
654  beamPhiPositions.at(aBeam) = phiPositions[aBeam];
655 
656  // the radii of the alignment tube beams for each halfbarrel.
657  // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
658  // in TIB/TOB modules these radii differ from the beam radius..
659  // ..due to the radial offsets (after the semitransparent mirrors)
660  const double radii[6] = {564., 564., 514., 514., 600., 600.};
661  std::vector<double> beamRadii(6, 0.);
662  for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel)
663  beamRadii.at(aHalfbarrel) = radii[aHalfbarrel];
664 
665  // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
666  double detReducedZ[2] = {0., 0.};
667  // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
668  double beamReducedZ[2] = {0., 0.};
669 
670  // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
671  detReducedZ[0] = nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() -
672  endFaceZPositions.at(det).at(theSide).at(0); // = zPrime
673  detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
674  detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
675  nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ(); // = zPrimePrime
676  detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
677 
678  // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
679  beamReducedZ[0] = (nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset) -
680  disk9EndFaceZPositions.at(0); // = ZTPrime
681  beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
682  beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
683  (nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset); // ZTPrimePrime
684  beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
685 
686  // the correction to phi from the endcap algorithm;
687  // it is defined such that the correction is to be subtracted ///////////////////////////////// ???
688  double phiCorrection = 0.;
689 
690  // contribution from phi rotation of first end face
691  phiCorrection += detReducedZ[1] * alignmentParameters.GetParameter(halfbarrel, 0, 0).first;
692 
693  // contribution from phi rotation of second end face
694  phiCorrection += detReducedZ[0] * alignmentParameters.GetParameter(halfbarrel, 1, 0).first;
695 
696  // contribution from translation along x of first endface
697  phiCorrection += detReducedZ[1] * sin(beamPhiPositions.at(beam)) *
698  alignmentParameters.GetParameter(halfbarrel, 0, 1).first / beamRadii.at(halfbarrel);
699 
700  // contribution from translation along x of second endface
701  phiCorrection += detReducedZ[0] * sin(beamPhiPositions.at(beam)) *
702  alignmentParameters.GetParameter(halfbarrel, 1, 1).first / beamRadii.at(halfbarrel);
703 
704  // contribution from translation along y of first endface
705  phiCorrection -= detReducedZ[1] * cos(beamPhiPositions.at(beam)) *
706  alignmentParameters.GetParameter(halfbarrel, 0, 2).first / beamRadii.at(halfbarrel);
707 
708  // contribution from translation along y of second endface
709  phiCorrection -= detReducedZ[0] * cos(beamPhiPositions.at(beam)) *
710  alignmentParameters.GetParameter(halfbarrel, 1, 2).first / beamRadii.at(halfbarrel);
711 
712  // contribution from beam parameters;
713  // originally, the contribution in meter is proportional to the radius of the beams: beamRadii.at( 0 )
714  // the additional factor: beamRadii.at( halfbarrel ) converts from meter to radian on the module
715  phiCorrection += beamReducedZ[1] * alignmentParameters.GetBeamParameter(beam, 0).first * beamRadii.at(0) /
716  beamRadii.at(halfbarrel);
717  phiCorrection += beamReducedZ[0] * alignmentParameters.GetBeamParameter(beam, 1).first * beamRadii.at(0) /
718  beamRadii.at(halfbarrel);
719 
720  return phiCorrection;
721 }
722 
732 void LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile(const char* filename,
733  LASGlobalData<LASCoordinateSet>& measuredCoordinates,
734  LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
735  std::ifstream file(filename);
736  if (file.bad()) {
737  std::cerr << " [LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile] ** ERROR: cannot open file \"" << filename
738  << "\"." << std::endl;
739  return;
740  }
741 
742  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
743  "@@@@@@@@@@@"
744  << std::endl;
745  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
746  "@@@@@@@@@@@"
747  << std::endl;
748  std::cerr << " [LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile] ** WARNING: you are reading a fake measurement "
749  "from a file!"
750  << std::endl;
751  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
752  "@@@@@@@@@@@"
753  << std::endl;
754  std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
755  "@@@@@@@@@@@"
756  << std::endl;
757 
758  // the measured coordinates will finally be overwritten;
759  // first, set them to the nominal values
760  measuredCoordinates = nominalCoordinates;
761 
762  // and put large errors on all values;
763  {
764  LASGlobalLoop moduleLoop;
765  int det, ring, beam, disk, pos;
766 
767  det = 0;
768  ring = 0;
769  beam = 0;
770  disk = 0;
771  do {
772  measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhiError(1000.);
773  } while (moduleLoop.TECLoop(det, ring, beam, disk));
774 
775  det = 2;
776  beam = 0;
777  pos = 0;
778  do {
779  measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhiError(1000.);
780  } while (moduleLoop.TIBTOBLoop(det, beam, pos));
781 
782  det = 0;
783  beam = 0;
784  disk = 0;
785  do {
786  measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhiError(1000.);
787  } while (moduleLoop.TEC2TECLoop(det, beam, disk));
788  }
789 
790  // buffers for read-in
791  int det, beam, z, ring;
792  double phi, phiError;
793 
794  while (!file.eof()) {
795  file >> det;
796  if (file.eof())
797  break; // do not read the last line twice, do not fill trash if file empty
798 
799  file >> beam;
800  file >> z;
801  file >> ring;
802  file >> phi;
803  file >> phiError;
804 
805  if (det > 1) { // TIB/TOB
806  measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhi(phi);
807  measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhiError(phiError);
808  } else { // TEC or TEC(at)
809  if (ring > -1) { // TEC
810  measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhi(phi);
811  measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhiError(phiError);
812  } else { // TEC(at)
813  measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhi(phi);
814  measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhiError(phiError);
815  }
816  }
817  }
818 
819  file.close();
820 }
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
ALPAKA_ACCELERATOR_NAMESPACE::brokenline::radii
results radii
Definition: BrokenLine.h:180
LASCoordinateSet::SetPhiError
void SetPhiError(double aPhiError)
Definition: LASCoordinateSet.h:32
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
LASCoordinateSet::GetZ
double GetZ(void) const
Definition: LASCoordinateSet.h:22
LASGlobalLoop::TEC2TECLoop
bool TEC2TECLoop(int &, int &, int &) const
Definition: LASGlobalLoop.cc:74
EcalCondDBWriter_cfi.beam
beam
Definition: EcalCondDBWriter_cfi.py:45
run3scouting_cff.phiError
phiError
Definition: run3scouting_cff.py:149
LASGlobalLoop::TECLoop
bool TECLoop(int &, int &, int &, int &) const
Definition: LASGlobalLoop.cc:13
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
LASBarrelAlignmentParameterSet::GetBeamParameter
std::pair< double, double > & GetBeamParameter(int aBeam, int aParameter)
Definition: LASBarrelAlignmentParameterSet.cc:75
LASGlobalData::GetTIBTOBEntry
T & GetTIBTOBEntry(int subdetector, int beam, int tibTobPosition)
Definition: LASGlobalData.h:107
LASAlignmentTubeAlgorithm::GetTEC2TECAlignmentParameterCorrection
double GetTEC2TECAlignmentParameterCorrection(int, int, int, LASGlobalData< LASCoordinateSet > &, LASBarrelAlignmentParameterSet &)
Definition: LASAlignmentTubeAlgorithm.cc:609
LASGlobalLoop::TIBTOBLoop
bool TIBTOBLoop(int &, int &, int &) const
Definition: LASGlobalLoop.cc:46
LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile
void ReadMisalignmentFromFile(const char *, LASGlobalData< LASCoordinateSet > &, LASGlobalData< LASCoordinateSet > &)
Definition: LASAlignmentTubeAlgorithm.cc:732
LASCoordinateSet::GetPhi
double GetPhi(void) const
Definition: LASCoordinateSet.h:18
LASGlobalData::GetTEC2TECEntry
T & GetTEC2TECEntry(int subdetector, int beam, int tecDisk)
Definition: LASGlobalData.h:129
LASAlignmentTubeAlgorithm::CalculateParameters
LASBarrelAlignmentParameterSet CalculateParameters(LASGlobalData< LASCoordinateSet > &, LASGlobalData< LASCoordinateSet > &)
Definition: LASAlignmentTubeAlgorithm.cc:12
LASGlobalData::GetTECEntry
T & GetTECEntry(int subdetector, int tecRing, int beam, int tecDisk)
Definition: LASGlobalData.h:84
DMR_cfg.cerr
cerr
Definition: DMR_cfg.py:78
gather_cfg.cout
cout
Definition: gather_cfg.py:144
LASAlignmentTubeAlgorithm::GetTIBTOBAlignmentParameterCorrection
double GetTIBTOBAlignmentParameterCorrection(int, int, int, LASGlobalData< LASCoordinateSet > &, LASBarrelAlignmentParameterSet &)
Definition: LASAlignmentTubeAlgorithm.cc:490
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
Generated for CMSSW Reference Manual by doxygen 1.8.14