CMS 3D CMS Logo
List of all members | Public Member Functions
LASEndcapAlgorithm Class Reference

#include <LASEndcapAlgorithm.h>

Public Member Functions

LASEndcapAlignmentParameterSet CalculateParameters (LASGlobalData< LASCoordinateSet > &, LASGlobalData< LASCoordinateSet > &)
 
double GetAlignmentParameterCorrection (int, int, int, int, LASGlobalData< LASCoordinateSet > &, LASEndcapAlignmentParameterSet &)
 
 LASEndcapAlgorithm ()
 

Detailed Description

calculate parameters for both endcaps from measurement

TODO:

Definition at line 25 of file LASEndcapAlgorithm.h.

Constructor & Destructor Documentation

◆ LASEndcapAlgorithm()

LASEndcapAlgorithm::LASEndcapAlgorithm ( )

Definition at line 7 of file LASEndcapAlgorithm.cc.

7 {}

Member Function Documentation

◆ CalculateParameters()

LASEndcapAlignmentParameterSet LASEndcapAlgorithm::CalculateParameters ( LASGlobalData< LASCoordinateSet > &  measuredCoordinates,
LASGlobalData< LASCoordinateSet > &  nominalCoordinates 
)

implementation of the analytical solution for the endcap; described in bruno's thesis (Appendix B): http://darwin.bth.rwth-aachen.de/opus3/volltexte/2002/348 but extended with the beams' phi positions

Definition at line 15 of file LASEndcapAlgorithm.cc.

16  {
17  std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Starting." << std::endl;
18 
19  // loop object
20  LASGlobalLoop globalLoop;
21  int det, ring, beam, disk;
22 
23  // vector containing the z positions of the disks in mm;
24  // outer vector: TEC+/-, inner vector: 9 disks
25  const double zPositions[9] = {1322.5, 1462.5, 1602.5, 1742.5, 1882.5, 2057.5, 2247.5, 2452.5, 2667.5};
26  std::vector<std::vector<double> > diskZPositions(2, std::vector<double>(9, 0.));
27  for (det = 0; det < 2; ++det) {
28  for (disk = 0; disk < 9; ++disk) {
29  // sign depends on side of course
30  diskZPositions.at(det).at(disk) = (det == 0 ? zPositions[disk] : -1. * zPositions[disk]);
31  }
32  }
33 
34  // vector containing the phi positions of the beam in rad;
35  // outer vector: TEC+/-, inner vector: 8 beams
36  const double phiPositions[8] = {0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486};
37  std::vector<std::vector<double> > beamPhiPositions(2, std::vector<double>(8, 0.));
38  for (det = 0; det < 2; ++det) {
39  for (beam = 0; beam < 8; ++beam) {
40  beamPhiPositions.at(det).at(beam) = phiPositions[beam];
41  }
42  }
43 
44  // vector containing the radius of ring4,ring6 = (0,1)
45  std::vector<double> radius(2, 0.);
46  radius.at(0) = 564.;
47  radius.at(1) = 840.;
48 
49  // constants
50  const double endcapLength = 1345.; // mm
51 
52  // now come some sums which are later used in the formulas for the parameters.
53  // the rings are implicitly summed over, however, this brings a little complication:
54  // to make the calculation of the parameters independent of the ring (=radius),
55  // we define some of the sums twice, once for phi and once for y=r*phi
56 
57  // sum over r*phi or phi for each endcap and for each disk (both rings)
58  // outer vector: TEC+/-, inner vector: 9 disks
59  std::vector<std::vector<double> > sumOverY(2, std::vector<double>(9, 0.));
60  std::vector<std::vector<double> > sumOverPhi(2, std::vector<double>(9, 0.));
61 
62  // sum over phi for each endcap and for each beam (both rings)
63  // outer vector: TEC+/-, inner vector: 8 beams
64  std::vector<std::vector<double> > kSumOverPhi(2, std::vector<double>(8, 0.));
65 
66  // double sum over r*phi or phi, for each endcap (both rings)
67  // outer vector: TEC+/-
68  std::vector<double> doubleSumOverY(2, 0.);
69  std::vector<double> doubleSumOverPhi(2, 0.);
70 
71  // sum over r*phi*z or phi*z, for each endcap and for each beam (both rings)
72  // outer vector: TEC+/-, inner vector: 8 beams
73  std::vector<std::vector<double> > kSumOverPhiZ(2, std::vector<double>(8, 0.));
74 
75  // sum over r*phi*z or phi*z, for each endcap (both rings)
76  // outer vector: TEC+/-
77  std::vector<double> doubleSumOverYZ(2, 0.);
78  std::vector<double> doubleSumOverPhiZ(2, 0.);
79 
80  // sum over sin(phi_nominal)*R*phi for each endcap and for each disk (both rings)
81  std::vector<std::vector<double> > sumOverSinThetaY(2, std::vector<double>(9, 0.));
82 
83  // sum over cos(phi_nominal)*R*phi for each endcap and for each disk (both rings)
84  std::vector<std::vector<double> > sumOverCosThetaY(2, std::vector<double>(9, 0.));
85 
86  // double sum over sin or cos(phi_nominal)*phi, for each endcap
87  std::vector<double> doubleSumOverSinThetaY(2, 0.);
88  std::vector<double> doubleSumOverCosThetaY(2, 0.);
89 
90  // double sum over sin or cos(phi_nominal)*phi*z, for each endcap
91  std::vector<double> doubleSumOverSinThetaYZ(2, 0.);
92  std::vector<double> doubleSumOverCosThetaYZ(2, 0.);
93 
94  // sum over z values / sum over z^2 values
95  std::vector<double> sumOverZ(2, 0.);
96  std::vector<double> sumOverZSquared(2, 0.);
97 
98  // now calculate the sums
99  det = 0;
100  ring = 0;
101  beam = 0;
102  disk = 0;
103  do {
104  if (ring == 1)
105  continue; //################################################################################################### BOUND TO RING6
106  // current radius, depends on the ring
107  const double radius = nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetR();
108 
109  // residual in r*phi (in the formulas this corresponds to y_ik)
110  const double residual = (measuredCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi() -
111  nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
112  radius;
113 
114  sumOverY.at(det).at(disk) += residual;
115  sumOverPhi.at(det).at(disk) += residual / radius;
116  kSumOverPhi.at(det).at(beam) += residual / radius;
117 
118  doubleSumOverY.at(det) += residual;
119  doubleSumOverPhi.at(det) += residual / radius;
120 
121  kSumOverPhiZ.at(det).at(beam) += diskZPositions.at(det).at(disk) * residual / radius;
122 
123  doubleSumOverYZ.at(det) += diskZPositions.at(det).at(disk) * residual;
124  doubleSumOverPhiZ.at(det) += diskZPositions.at(det).at(disk) * residual / radius;
125 
126  sumOverSinThetaY.at(det).at(disk) += sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
127  sumOverCosThetaY.at(det).at(disk) += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
128 
129  doubleSumOverSinThetaY.at(det) += sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
130  doubleSumOverCosThetaY.at(det) += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
131 
132  doubleSumOverSinThetaYZ.at(det) += sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
133  diskZPositions.at(det).at(disk) * residual;
134  doubleSumOverCosThetaYZ.at(det) += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
135  diskZPositions.at(det).at(disk) * residual;
136 
137  } while (globalLoop.TECLoop(det, ring, beam, disk));
138 
139  // disk-wise sums
140  for (disk = 0; disk < 9; ++disk) {
141  sumOverZ.at(0) += diskZPositions.at(0).at(disk);
142  sumOverZ.at(1) += diskZPositions.at(1).at(disk);
143  sumOverZSquared.at(0) += pow(diskZPositions.at(0).at(disk), 2);
144  sumOverZSquared.at(1) += pow(diskZPositions.at(1).at(disk), 2);
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  // ..and fill them;
185  // the additional divisors "/ 2." come from the fact that we average over both rings
186  for (det = 0; det < 2; ++det) { // TEC+/- loop
187 
188  // deltaPhi_0
189  // here we use the phi-sums to eliminate the radius
190  deltaPhi0.at(det) =
191  (sumOverZSquared.at(det) * doubleSumOverPhi.at(det) - sumOverZ.at(det) * doubleSumOverPhiZ.at(det)) /
192  (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det))); // / 2.; // @@@@@@@
193 
194  // deltaPhi_t
195  // again use the phi-sums
196  deltaPhiT.at(det) = endcapLength * (9. * doubleSumOverPhiZ.at(det) - sumOverZ.at(det) * doubleSumOverPhi.at(det)) /
197  (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det))); // / 2.; // @@@@@@@
198 
199  // deltaPhi_k (k=0..8)
200  // again use the phi-sums
201  for (disk = 0; disk < 9; ++disk) {
202  deltaPhiK.at(det).at(disk) = (-1. * diskZPositions.at(det).at(disk) * deltaPhiT.at(det) / endcapLength) -
203  (deltaPhi0.at(det)) - sumOverPhi.at(det).at(disk) / 8.; // / 2.; // @@@@@@@
204  }
205 
206  // deltaX_0
207  deltaX0.at(det) = 2. *
208  (sumOverZ.at(det) * doubleSumOverSinThetaYZ.at(det) -
209  sumOverZSquared.at(det) * doubleSumOverSinThetaY.at(det)) /
210  (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det))); // / 2.; // @@@@@@@
211 
212  // deltaX_t
213  deltaXT.at(det) = 2. * endcapLength *
214  (sumOverZ.at(det) * doubleSumOverSinThetaY.at(det) - 9. * doubleSumOverSinThetaYZ.at(det)) /
215  (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det))); // / 2.; // @@@@@@@
216 
217  // deltaX_k (k=0..8)
218  for (disk = 0; disk < 9; ++disk) {
219  deltaXK.at(det).at(disk) = (-1. * diskZPositions.at(det).at(disk) * deltaXT.at(det) / endcapLength) -
220  (deltaX0.at(det)) + 2. * sumOverSinThetaY.at(det).at(disk) / 8.; // / 2.; // @@@@@@@
221  }
222 
223  // deltaY_0
224  deltaY0.at(det) = 2. *
225  (sumOverZSquared.at(det) * doubleSumOverCosThetaY.at(det) -
226  sumOverZ.at(det) * doubleSumOverCosThetaYZ.at(det)) /
227  (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det))); // / 2.; // @@@@@@@
228 
229  // deltaY_t
230  deltaYT.at(det) = 2. * endcapLength *
231  (9. * doubleSumOverCosThetaYZ.at(det) - sumOverZ.at(det) * doubleSumOverCosThetaY.at(det)) /
232  (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det))); // / 2.; // @@@@@@@
233 
234  // deltaY_k (k=0..8)
235  for (disk = 0; disk < 9; ++disk) {
236  deltaYK.at(det).at(disk) = (-1. * diskZPositions.at(det).at(disk) * deltaYT.at(det) / endcapLength) -
237  (deltaY0.at(det)) - 2. * sumOverCosThetaY.at(det).at(disk) / 8.; // / 2.; // @@@@@@@
238  }
239 
240  // the beam parameters deltaTA & deltaTB
241  for (beam = 0; beam < 8; ++beam) {
242  deltaTA.at(det).at(beam) =
243  deltaPhi0.at(det) -
244  (kSumOverPhi.at(det).at(beam) * sumOverZSquared.at(det) - kSumOverPhiZ.at(det).at(beam) * sumOverZ.at(det)) /
245  (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)) +
246  (cos(beamPhiPositions.at(det).at(beam)) * deltaY0.at(det) -
247  sin(beamPhiPositions.at(det).at(beam)) * deltaX0.at(det)) /
248  radius.at(0); // for ring 4..
249  // + ( cos( beamPhiPositions.at( det ).at( beam ) ) * deltaY0.at( det ) - sin( beamPhiPositions.at( det ).at( beam ) ) * deltaX0.at( det ) ) / radius.at( 1 ); // ...and ring 6
250 
251  deltaTB.at(det).at(beam) =
252  -1. * deltaPhiT.at(det) - deltaPhi0.at(det) -
253  (kSumOverPhi.at(det).at(beam) * sumOverZ.at(det) - 9. * kSumOverPhiZ.at(det).at(beam)) * endcapLength /
254  (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)) -
255  (kSumOverPhiZ.at(det).at(beam) * sumOverZ.at(det) - kSumOverPhi.at(det).at(beam) * sumOverZSquared.at(det)) /
256  (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)) +
257  ((deltaXT.at(det) + deltaX0.at(det)) * sin(beamPhiPositions.at(det).at(beam)) -
258  (deltaYT.at(det) + deltaY0.at(det)) * cos(beamPhiPositions.at(det).at(beam))) /
259  radius.at(0); // for ring4..
260  // + ( ( deltaXT.at( det ) + deltaX0.at( det ) ) * sin( beamPhiPositions.at( det ).at( beam ) ) - ( deltaYT.at( det ) + deltaY0.at( det ) ) * cos( beamPhiPositions.at( det ).at( beam ) ) )
261  // / radius.at( 1 ); // ..and ring6
262  }
263  }
264 
265  // fill the result
266  LASEndcapAlignmentParameterSet theResult;
267 
268  // for the moment we fill only the values, not the errors
269 
270  // disk parameters
271  for (det = 0; det < 2; ++det) {
272  for (disk = 0; disk < 9; ++disk) {
273  // the rotation parameters: deltaPhi_k
274  theResult.GetDiskParameter(det, disk, 0).first = deltaPhiK.at(det).at(disk);
275  // the x offsets: deltaX_k
276  theResult.GetDiskParameter(det, disk, 1).first = deltaXK.at(det).at(disk);
277  // the y offsets: deltaY_k
278  theResult.GetDiskParameter(det, disk, 2).first = deltaYK.at(det).at(disk);
279  }
280  }
281 
282  // global parameters
283  for (int det = 0; det < 2; ++det) {
284  theResult.GetGlobalParameter(det, 0).first = deltaPhi0.at(det);
285  theResult.GetGlobalParameter(det, 1).first = deltaPhiT.at(det);
286  theResult.GetGlobalParameter(det, 2).first = deltaX0.at(det);
287  theResult.GetGlobalParameter(det, 3).first = deltaXT.at(det);
288  theResult.GetGlobalParameter(det, 4).first = deltaY0.at(det);
289  theResult.GetGlobalParameter(det, 5).first = deltaYT.at(det);
290  }
291 
292  // beam parameters
293  for (int det = 0; det < 2; ++det) {
294  for (int beam = 0; beam < 8; ++beam) {
295  theResult.GetBeamParameter(det, 1 /*R6*/, beam, 0).first =
296  deltaTA.at(det).at(beam);
297  theResult.GetBeamParameter(det, 1 /*R6*/, beam, 1).first =
298  deltaTB.at(det).at(beam);
299  }
300  }
301 
302  std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Done." << std::endl;
303 
304  return (theResult);
305 }

References EcalCondDBWriter_cfi::beam, funct::cos(), gather_cfg::cout, LASEndcapAlignmentParameterSet::GetBeamParameter(), LASEndcapAlignmentParameterSet::GetDiskParameter(), LASEndcapAlignmentParameterSet::GetGlobalParameter(), LASCoordinateSet::GetPhi(), LASCoordinateSet::GetR(), LASGlobalData< T >::GetTECEntry(), funct::pow(), CosmicsPD_Skims::radius, relativeConstraints::ring, funct::sin(), and LASGlobalLoop::TECLoop().

Referenced by LaserAlignment::endRunProduce().

◆ GetAlignmentParameterCorrection()

double LASEndcapAlgorithm::GetAlignmentParameterCorrection ( int  det,
int  ring,
int  beam,
int  disk,
LASGlobalData< LASCoordinateSet > &  nominalCoordinates,
LASEndcapAlignmentParameterSet &  endcapParameters 
)

for a given set of endcap alignment parameters "endcapParameters", this function returns the global phi offset from nominalPosition for a module specified by (det,ring,beam,disk)

Definition at line 312 of file LASEndcapAlgorithm.cc.

317  {
318  // ring dependent radius, to be softcoded...
319  const double radius = ring == 0 ? 564. : 840.;
320  const double endcapLength = 1345.; // mm
321 
322  // the correction to phi from the endcap algorithm;
323  // it is defined such that the correction is to be subtracted
324  double phiCorrection = 0.;
325 
326  // plain disk phi
327  phiCorrection += endcapParameters.GetDiskParameter(det, disk, 0).first;
328 
329  // phi component from x deviation
330  phiCorrection -= sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
331  endcapParameters.GetDiskParameter(det, disk, 1).first;
332 
333  // phi component from y deviation
334  phiCorrection += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
335  endcapParameters.GetDiskParameter(det, disk, 2).first;
336 
337  // phi correction from global phi
338  phiCorrection += endcapParameters.GetGlobalParameter(det, 0).first;
339 
340  // correction from global x deviation
341  phiCorrection -= sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
342  endcapParameters.GetGlobalParameter(det, 2).first;
343 
344  // correction from global y deviation
345  phiCorrection += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
346  endcapParameters.GetGlobalParameter(det, 4).first;
347 
348  // correction from global torsion
349  phiCorrection += nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength *
350  endcapParameters.GetGlobalParameter(det, 1).first;
351 
352  // correction from global x shear
353  phiCorrection -= nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength / radius *
354  sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
355  endcapParameters.GetGlobalParameter(det, 3).first;
356 
357  // correction from global y shear
358  phiCorrection += nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength / radius *
359  cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
360  endcapParameters.GetGlobalParameter(det, 5).first;
361 
362  // correction from beam parameters
363  phiCorrection += (nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength - 1.) *
364  endcapParameters.GetBeamParameter(det, 1, beam, 0).first;
365  phiCorrection += nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength *
366  endcapParameters.GetBeamParameter(det, 1, beam, 1).first;
367 
368  return phiCorrection;
369 }

References EcalCondDBWriter_cfi::beam, funct::cos(), LASEndcapAlignmentParameterSet::GetBeamParameter(), LASEndcapAlignmentParameterSet::GetDiskParameter(), LASEndcapAlignmentParameterSet::GetGlobalParameter(), LASCoordinateSet::GetPhi(), LASGlobalData< T >::GetTECEntry(), LASCoordinateSet::GetZ(), CosmicsPD_Skims::radius, relativeConstraints::ring, and funct::sin().

Referenced by LaserAlignment::ApplyEndcapMaskingCorrections().

LASEndcapAlignmentParameterSet::GetGlobalParameter
std::pair< double, double > & GetGlobalParameter(int aSubdetector, int aParameter)
Definition: LASEndcapAlignmentParameterSet.cc:97
LASEndcapAlignmentParameterSet
Definition: LASEndcapAlignmentParameterSet.h:25
gather_cfg.cout
cout
Definition: gather_cfg.py:144
LASCoordinateSet::GetZ
double GetZ(void) const
Definition: LASCoordinateSet.h:22
LASGlobalLoop
Definition: LASGlobalLoop.h:24
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
LASGlobalData::GetTECEntry
T & GetTECEntry(int subdetector, int tecRing, int beam, int tecDisk)
Definition: LASGlobalData.h:84
EcalCondDBWriter_cfi.beam
beam
Definition: EcalCondDBWriter_cfi.py:45
LASGlobalLoop::TECLoop
bool TECLoop(int &, int &, int &, int &) const
Definition: LASGlobalLoop.cc:13
LASEndcapAlignmentParameterSet::GetBeamParameter
std::pair< double, double > & GetBeamParameter(int aSubdetector, int aRing, int aBeam, int aParameter)
Definition: LASEndcapAlignmentParameterSet.cc:122
relativeConstraints.ring
ring
Definition: relativeConstraints.py:68
CosmicsPD_Skims.radius
radius
Definition: CosmicsPD_Skims.py:135
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
LASCoordinateSet::GetPhi
double GetPhi(void) const
Definition: LASCoordinateSet.h:18
LASCoordinateSet::GetR
double GetR(void) const
Definition: LASCoordinateSet.h:20
LASEndcapAlignmentParameterSet::GetDiskParameter
std::pair< double, double > & GetDiskParameter(int aSubdetector, int aDisk, int aParameter)
Definition: LASEndcapAlignmentParameterSet.cc:63
Generated for CMSSW Reference Manual by doxygen 1.8.16