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