CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
LASEndcapAlgorithm.cc
Go to the documentation of this file.
1 
3 
4 
5 
10 }
11 
12 
13 
14 
15 
23  LASGlobalData<LASCoordinateSet>& nominalCoordinates ) {
24 
25  std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Starting." << std::endl;
26 
27  // loop object
28  LASGlobalLoop globalLoop;
29  int det, ring, beam, disk;
30 
31  // vector containing the z positions of the disks in mm;
32  // outer vector: TEC+/-, inner vector: 9 disks
33  const double zPositions[9] = { 1322.5, 1462.5, 1602.5, 1742.5, 1882.5, 2057.5, 2247.5, 2452.5, 2667.5 };
34  std::vector<std::vector<double> > diskZPositions( 2, std::vector<double>( 9, 0. ) );
35  for( det = 0; det < 2; ++det ) {
36  for( disk = 0; disk < 9; ++disk ) {
37  // sign depends on side of course
38  diskZPositions.at( det ).at( disk ) = ( det==0 ? zPositions[disk] : -1. * zPositions[disk] );
39  }
40  }
41 
42  // vector containing the phi positions of the beam in rad;
43  // outer vector: TEC+/-, inner vector: 8 beams
44  const double phiPositions[8] = { 0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486 };
45  std::vector<std::vector<double> > beamPhiPositions( 2, std::vector<double>( 8, 0. ) );
46  for( det = 0; det < 2; ++ det ) {
47  for( beam = 0; beam < 8; ++beam ) {
48  beamPhiPositions.at( det ).at( beam ) = phiPositions[beam];
49  }
50  }
51 
52  // vector containing the radius of ring4,ring6 = (0,1)
53  std::vector<double> radius( 2, 0. );
54  radius.at( 0 ) = 564.; radius.at( 1 ) = 840.;
55 
56  // constants
57  const double endcapLength = 1345.; // mm
58 
59  // now come some sums which are later used in the formulas for the parameters.
60  // the rings are implicitly summed over, however, this brings a little complication:
61  // to make the calculation of the parameters independent of the ring (=radius),
62  // we define some of the sums twice, once for phi and once for y=r*phi
63 
64  // sum over r*phi or phi for each endcap and for each disk (both rings)
65  // outer vector: TEC+/-, inner vector: 9 disks
66  std::vector<std::vector<double> > sumOverY( 2, std::vector<double>( 9, 0. ) );
67  std::vector<std::vector<double> > sumOverPhi( 2, std::vector<double>( 9, 0. ) );
68 
69  // sum over phi for each endcap and for each beam (both rings)
70  // outer vector: TEC+/-, inner vector: 8 beams
71  std::vector<std::vector<double> > kSumOverPhi( 2, std::vector<double>( 8, 0. ) );
72 
73  // double sum over r*phi or phi, for each endcap (both rings)
74  // outer vector: TEC+/-
75  std::vector<double> doubleSumOverY( 2, 0. );
76  std::vector<double> doubleSumOverPhi( 2, 0. );
77 
78  // sum over r*phi*z or phi*z, for each endcap and for each beam (both rings)
79  // outer vector: TEC+/-, inner vector: 8 beams
80  std::vector<std::vector<double> > kSumOverPhiZ( 2, std::vector<double>( 8, 0. ) );
81 
82  // sum over r*phi*z or phi*z, for each endcap (both rings)
83  // outer vector: TEC+/-
84  std::vector<double> doubleSumOverYZ( 2, 0. );
85  std::vector<double> doubleSumOverPhiZ( 2, 0. );
86 
87  // sum over sin(phi_nominal)*R*phi for each endcap and for each disk (both rings)
88  std::vector<std::vector<double> > sumOverSinThetaY( 2, std::vector<double>( 9, 0. ) );
89 
90  // sum over cos(phi_nominal)*R*phi for each endcap and for each disk (both rings)
91  std::vector<std::vector<double> > sumOverCosThetaY( 2, std::vector<double>( 9, 0. ) );
92 
93  // double sum over sin or cos(phi_nominal)*phi, for each endcap
94  std::vector<double> doubleSumOverSinThetaY( 2, 0. );
95  std::vector<double> doubleSumOverCosThetaY( 2, 0. );
96 
97  // double sum over sin or cos(phi_nominal)*phi*z, for each endcap
98  std::vector<double> doubleSumOverSinThetaYZ( 2, 0. );
99  std::vector<double> doubleSumOverCosThetaYZ( 2, 0. );
100 
101  // sum over z values / sum over z^2 values
102  std::vector<double> sumOverZ( 2, 0. );
103  std::vector<double> sumOverZSquared( 2, 0. );
104 
105 
106  // now calculate the sums
107  det = 0; ring = 0; beam = 0; disk = 0;
108  do {
109  if( ring == 1 ) continue; //################################################################################################### BOUND TO RING6
110  // current radius, depends on the ring
111  const double radius = nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetR();
112 
113  // residual in r*phi (in the formulas this corresponds to y_ik)
114  const double residual = ( measuredCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() - nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * radius;
115 
116  sumOverY.at( det ).at( disk ) += residual;
117  sumOverPhi.at( det ).at( disk ) += residual / radius;
118  kSumOverPhi.at( det ).at( beam ) += residual / radius;
119 
120  doubleSumOverY.at( det ) += residual;
121  doubleSumOverPhi.at( det ) += residual / radius;
122 
123  kSumOverPhiZ.at( det ).at( beam ) += diskZPositions.at( det ).at( disk ) * residual / radius;
124 
125  doubleSumOverYZ.at( det ) += diskZPositions.at( det ).at( disk ) * residual;
126  doubleSumOverPhiZ.at( det ) += diskZPositions.at( det ).at( disk ) * residual / radius;
127 
128  sumOverSinThetaY.at( det ).at( disk ) += sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
129  sumOverCosThetaY.at( det ).at( disk ) += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
130 
131  doubleSumOverSinThetaY.at( det ) += sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
132  doubleSumOverCosThetaY.at( det ) += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
133 
134  doubleSumOverSinThetaYZ.at( det ) += sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * diskZPositions.at( det ).at( disk ) * residual;
135  doubleSumOverCosThetaYZ.at( det ) += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * diskZPositions.at( det ).at( disk ) * residual;
136 
137  } while( globalLoop.TECLoop( det, ring, beam, disk ) );
138 
139 
140  // disk-wise sums
141  for( disk = 0; disk < 9; ++disk ) {
142  sumOverZ.at( 0 ) += diskZPositions.at( 0 ).at( disk ); sumOverZ.at( 1 ) += diskZPositions.at( 1 ).at( disk );
143  sumOverZSquared.at( 0 ) += pow( diskZPositions.at( 0 ).at( disk ), 2 ); sumOverZSquared.at( 1 ) += pow( diskZPositions.at( 1 ).at( disk ), 2 );
144  }
145 
146 
147  // now we can calculate the parameters for both TECs simultaneously,
148  // so they're all vectors( 2 ) for TEC+/- (global parameters), or dim 2*9 (disk parameters),
149  // or dim 2*8 (beam parameters)
150 
151  // define them..
152 
153  // deltaPhi_0
154  std::vector<double> deltaPhi0( 2, 0. );
155 
156  // deltaPhi_t
157  std::vector<double> deltaPhiT( 2, 0. );
158 
159  // deltaPhi_k (k=0..8)
160  std::vector<std::vector<double> > deltaPhiK( 2, std::vector<double>( 9, 0. ) );
161 
162  // deltaX_0
163  std::vector<double> deltaX0( 2, 0. );
164 
165  // deltaX_t
166  std::vector<double> deltaXT( 2, 0. );
167 
168  // deltaX_k (k=0..8)
169  std::vector<std::vector<double> > deltaXK( 2, std::vector<double>( 9, 0. ) );
170 
171  // deltaY_0
172  std::vector<double> deltaY0( 2, 0. );
173 
174  // deltaY_t
175  std::vector<double> deltaYT( 2, 0. );
176 
177  // deltaY_k (k=0..8)
178  std::vector<std::vector<double> > deltaYK( 2, std::vector<double>( 9, 0. ) );
179 
180  // beam parameters: deltaTheta_A, deltaTheta_B (i=0..7)
181  std::vector<std::vector<double> > deltaTA( 2, std::vector<double>( 8, 0. ) );
182  std::vector<std::vector<double> > deltaTB( 2, std::vector<double>( 8, 0. ) );
183 
184 
185 
186  // ..and fill them;
187  // the additional divisors "/ 2." come from the fact that we average over both rings
188  for( det = 0; det < 2; ++det ) { // TEC+/- loop
189 
190  // deltaPhi_0
191  // here we use the phi-sums to eliminate the radius
192  deltaPhi0.at( det ) = ( sumOverZSquared.at( det ) * doubleSumOverPhi.at( det ) - sumOverZ.at( det ) * doubleSumOverPhiZ.at( det ) )
193  / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) ); // / 2.; // @@@@@@@
194 
195  // deltaPhi_t
196  // again use the phi-sums
197  deltaPhiT.at( det ) = endcapLength * ( 9. * doubleSumOverPhiZ.at( det ) - sumOverZ.at( det ) * doubleSumOverPhi.at( det ) )
198  / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) ); // / 2.; // @@@@@@@
199 
200  // deltaPhi_k (k=0..8)
201  // again use the phi-sums
202  for( disk = 0; disk < 9; ++disk ) {
203  deltaPhiK.at( det ).at( disk ) = ( -1. * diskZPositions.at( det ).at( disk ) * deltaPhiT.at( det ) / endcapLength )
204  - ( deltaPhi0.at( det ) ) - sumOverPhi.at( det ).at( disk ) / 8.;// / 2.; // @@@@@@@
205  }
206 
207  // deltaX_0
208  deltaX0.at( det ) = 2. * ( sumOverZ.at( det ) * doubleSumOverSinThetaYZ.at( det ) - sumOverZSquared.at( det ) * doubleSumOverSinThetaY.at( det ) )
209  / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
210 
211  // deltaX_t
212  deltaXT.at( det ) = 2. * endcapLength * ( sumOverZ.at( det ) * doubleSumOverSinThetaY.at( det ) - 9. * doubleSumOverSinThetaYZ.at( det ) )
213  / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
214 
215  // deltaX_k (k=0..8)
216  for( disk = 0; disk < 9; ++disk ) {
217  deltaXK.at( det ).at( disk ) = ( -1. * diskZPositions.at( det ).at( disk ) * deltaXT.at( det ) / endcapLength )
218  - ( deltaX0.at( det ) ) + 2. * sumOverSinThetaY.at( det ).at( disk ) / 8.;// / 2.; // @@@@@@@
219  }
220 
221  // deltaY_0
222  deltaY0.at( det ) = 2. * ( sumOverZSquared.at( det ) * doubleSumOverCosThetaY.at( det ) - sumOverZ.at( det ) * doubleSumOverCosThetaYZ.at( det ) )
223  / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
224 
225  // deltaY_t
226  deltaYT.at( det ) = 2. * endcapLength * ( 9. * doubleSumOverCosThetaYZ.at( det ) - sumOverZ.at( det ) * doubleSumOverCosThetaY.at( det ) )
227  / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
228 
229  // deltaY_k (k=0..8)
230  for( disk = 0; disk < 9; ++disk ) {
231  deltaYK.at( det ).at( disk ) = ( -1. * diskZPositions.at( det ).at( disk ) * deltaYT.at( det ) / endcapLength )
232  - ( deltaY0.at( det ) ) - 2. * sumOverCosThetaY.at( det ).at( disk ) / 8.;// / 2.; // @@@@@@@
233  }
234 
235 
236  // the beam parameters deltaTA & deltaTB
237  for( beam = 0; beam < 8; ++beam ) {
238 
239  deltaTA.at( det ).at( beam ) = deltaPhi0.at( det )
240  - ( kSumOverPhi.at( det ).at( beam ) * sumOverZSquared.at( det ) - kSumOverPhiZ.at( det ).at( beam ) * sumOverZ.at( det ) )
241  / ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) )
242  + ( cos( beamPhiPositions.at( det ).at( beam ) ) * deltaY0.at( det ) - sin( beamPhiPositions.at( det ).at( beam ) ) * deltaX0.at( det ) ) / radius.at( 0 ); // for ring 4..
243  // + ( cos( beamPhiPositions.at( det ).at( beam ) ) * deltaY0.at( det ) - sin( beamPhiPositions.at( det ).at( beam ) ) * deltaX0.at( det ) ) / radius.at( 1 ); // ...and ring 6
244 
245 
246  deltaTB.at( det ).at( beam ) = -1. * deltaPhiT.at( det ) - deltaPhi0.at( det )
247  - ( kSumOverPhi.at( det ).at( beam ) * sumOverZ.at( det ) - 9. * kSumOverPhiZ.at( det ).at( beam ) )
248  * endcapLength / ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) )
249  - ( kSumOverPhiZ.at( det ).at( beam ) * sumOverZ.at( det ) - kSumOverPhi.at( det ).at( beam ) * sumOverZSquared.at( det ) )
250  / ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) )
251  + ( ( deltaXT.at( det ) + deltaX0.at( det ) ) * sin( beamPhiPositions.at( det ).at( beam ) ) - ( deltaYT.at( det ) + deltaY0.at( det ) ) * cos( beamPhiPositions.at( det ).at( beam ) ) )
252  / radius.at( 0 ); // for ring4..
253  // + ( ( deltaXT.at( det ) + deltaX0.at( det ) ) * sin( beamPhiPositions.at( det ).at( beam ) ) - ( deltaYT.at( det ) + deltaY0.at( det ) ) * cos( beamPhiPositions.at( det ).at( beam ) ) )
254  // / radius.at( 1 ); // ..and ring6
255 
256  }
257 
258 
259  }
260 
261 
262  // fill the result
264 
265  // for the moment we fill only the values, not the errors
266 
267  // disk parameters
268  for( det = 0; det < 2; ++det ) {
269  for( disk = 0; disk < 9; ++disk ) {
270  // the rotation parameters: deltaPhi_k
271  theResult.GetDiskParameter( det, disk, 0 ).first = deltaPhiK.at( det ).at( disk );
272  // the x offsets: deltaX_k
273  theResult.GetDiskParameter( det, disk, 1 ).first = deltaXK.at( det ).at( disk );
274  // the y offsets: deltaY_k
275  theResult.GetDiskParameter( det, disk, 2 ).first = deltaYK.at( det ).at( disk );
276  }
277  }
278 
279  // global parameters
280  for( int det = 0; det < 2; ++det ) {
281  theResult.GetGlobalParameter( det, 0 ).first = deltaPhi0.at( det );
282  theResult.GetGlobalParameter( det, 1 ).first = deltaPhiT.at( det );
283  theResult.GetGlobalParameter( det, 2 ).first = deltaX0.at( det );
284  theResult.GetGlobalParameter( det, 3 ).first = deltaXT.at( det );
285  theResult.GetGlobalParameter( det, 4 ).first = deltaY0.at( det );
286  theResult.GetGlobalParameter( det, 5 ).first = deltaYT.at( det );
287  }
288 
289  // beam parameters
290  for( int det = 0; det < 2; ++det ) {
291  for( int beam = 0; beam < 8; ++beam ) {
292  theResult.GetBeamParameter( det, 1 /*R6*/, beam, 0 ).first = deltaTA.at( det ).at( beam );
293  theResult.GetBeamParameter( det, 1 /*R6*/, beam, 1 ).first = deltaTB.at( det ).at( beam );
294  }
295  }
296 
297  std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Done." << std::endl;
298 
299  return( theResult );
300 
301 }
302 
303 
304 
305 
306 
312 double LASEndcapAlgorithm::GetAlignmentParameterCorrection( int det, int ring, int beam, int disk, LASGlobalData<LASCoordinateSet>& nominalCoordinates, LASEndcapAlignmentParameterSet& endcapParameters ) {
313 
314  // ring dependent radius, to be softcoded...
315  const double radius = ring==0 ? 564. : 840.;
316  const double endcapLength = 1345.; // mm
317 
318  // the correction to phi from the endcap algorithm;
319  // it is defined such that the correction is to be subtracted
320  double phiCorrection = 0.;
321 
322  // plain disk phi
323  phiCorrection += endcapParameters.GetDiskParameter( det, disk, 0 ).first;
324 
325  // phi component from x deviation
326  phiCorrection -= sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetDiskParameter( det, disk, 1 ).first;
327 
328  // phi component from y deviation
329  phiCorrection += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetDiskParameter( det, disk, 2 ).first;
330 
331  // phi correction from global phi
332  phiCorrection += endcapParameters.GetGlobalParameter( det, 0 ).first;
333 
334  // correction from global x deviation
335  phiCorrection -= sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetGlobalParameter( det, 2 ).first;
336 
337  // correction from global y deviation
338  phiCorrection += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetGlobalParameter( det, 4 ).first;
339 
340  // correction from global torsion
341  phiCorrection += nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength * endcapParameters.GetGlobalParameter( det, 1 ).first;
342 
343  // correction from global x shear
344  phiCorrection -= nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength / radius *
345  sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * endcapParameters.GetGlobalParameter( det, 3 ).first;
346 
347  // correction from global y shear
348  phiCorrection += nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength / radius *
349  cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * endcapParameters.GetGlobalParameter( det, 5 ).first;
350 
351  // correction from beam parameters
352  phiCorrection += ( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength - 1. ) * endcapParameters.GetBeamParameter( det, 1, beam, 0 ).first;
353  phiCorrection += nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength * endcapParameters.GetBeamParameter( det, 1, beam, 1 ).first;
354 
355  return phiCorrection;
356 
357 }
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::pair< double, double > & GetGlobalParameter(int aSubdetector, int aParameter)
double GetPhi(void) const
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::pair< double, double > & GetDiskParameter(int aSubdetector, int aDisk, int aParameter)
LASEndcapAlignmentParameterSet CalculateParameters(LASGlobalData< LASCoordinateSet > &, LASGlobalData< LASCoordinateSet > &)
std::pair< double, double > & GetBeamParameter(int aSubdetector, int aRing, int aBeam, int aParameter)
T & GetTECEntry(int subdetector, int tecRing, int beam, int tecDisk)
Definition: LASGlobalData.h:91
double GetAlignmentParameterCorrection(int, int, int, int, LASGlobalData< LASCoordinateSet > &, LASEndcapAlignmentParameterSet &)
double GetZ(void) const
bool TECLoop(int &, int &, int &, int &) const
tuple cout
Definition: gather_cfg.py:121
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
double GetR(void) const