LaserAlignment/src/LASEndcapAlgorithm.cc

0001
0002 #include "Alignment/LaserAlignment/interface/LASEndcapAlgorithm.h"
0003
0004 ///
0005 ///
0006 ///
0007 LASEndcapAlgorithm::LASEndcapAlgorithm() {}
0008
0009 ///
0010 /// implementation of the analytical solution for the endcap;
0011 /// described in bruno's thesis (Appendix B):
0012 /// http://darwin.bth.rwth-aachen.de/opus3/volltexte/2002/348
0013 /// but extended with the beams' phi positions
0014 ///
0015 LASEndcapAlignmentParameterSet LASEndcapAlgorithm::CalculateParameters(
0016     LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
0017   std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Starting." << std::endl;
0018
0019   // loop object
0020   LASGlobalLoop globalLoop;
0021   int det, ring, beam, disk;
0022
0023   // vector containing the z positions of the disks in mm;
0024   // outer vector: TEC+/-, inner vector: 9 disks
0025   const double zPositions[9] = {1322.5, 1462.5, 1602.5, 1742.5, 1882.5, 2057.5, 2247.5, 2452.5, 2667.5};
0026   std::vector<std::vector<double> > diskZPositions(2, std::vector<double>(9, 0.));
0027   for (det = 0; det < 2; ++det) {
0028     for (disk = 0; disk < 9; ++disk) {
0029       // sign depends on side of course
0030       diskZPositions.at(det).at(disk) = (det == 0 ? zPositions[disk] : -1. * zPositions[disk]);
0031     }
0032   }
0033
0034   // vector containing the phi positions of the beam in rad;
0035   // outer vector: TEC+/-, inner vector: 8 beams
0036   const double phiPositions[8] = {0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486};
0037   std::vector<std::vector<double> > beamPhiPositions(2, std::vector<double>(8, 0.));
0038   for (det = 0; det < 2; ++det) {
0039     for (beam = 0; beam < 8; ++beam) {
0040       beamPhiPositions.at(det).at(beam) = phiPositions[beam];
0041     }
0042   }
0043
0044   // vector containing the radius of ring4,ring6 = (0,1)
0045   std::vector<double> radius(2, 0.);
0046   radius.at(0) = 564.;
0047   radius.at(1) = 840.;
0048
0049   // constants
0050   const double endcapLength = 1345.;  // mm
0051
0052   // now come some sums which are later used in the formulas for the parameters.
0053   // the rings are implicitly summed over, however, this brings a little complication:
0054   // to make the calculation of the parameters independent of the ring (=radius),
0055   // we define some of the sums twice, once for phi and once for y=r*phi
0056
0057   // sum over r*phi or phi for each endcap and for each disk (both rings)
0058   // outer vector: TEC+/-, inner vector: 9 disks
0059   std::vector<std::vector<double> > sumOverY(2, std::vector<double>(9, 0.));
0060   std::vector<std::vector<double> > sumOverPhi(2, std::vector<double>(9, 0.));
0061
0062   // sum over phi for each endcap and for each beam (both rings)
0063   // outer vector: TEC+/-, inner vector: 8 beams
0064   std::vector<std::vector<double> > kSumOverPhi(2, std::vector<double>(8, 0.));
0065
0066   // double sum over r*phi or phi, for each endcap (both rings)
0067   // outer vector: TEC+/-
0068   std::vector<double> doubleSumOverY(2, 0.);
0069   std::vector<double> doubleSumOverPhi(2, 0.);
0070
0071   // sum over r*phi*z or phi*z, for each endcap and for each beam (both rings)
0072   // outer vector: TEC+/-, inner vector: 8 beams
0073   std::vector<std::vector<double> > kSumOverPhiZ(2, std::vector<double>(8, 0.));
0074
0075   // sum over r*phi*z or phi*z, for each endcap (both rings)
0076   // outer vector: TEC+/-
0077   std::vector<double> doubleSumOverYZ(2, 0.);
0078   std::vector<double> doubleSumOverPhiZ(2, 0.);
0079
0080   // sum over sin(phi_nominal)*R*phi for each endcap and for each disk (both rings)
0081   std::vector<std::vector<double> > sumOverSinThetaY(2, std::vector<double>(9, 0.));
0082
0083   // sum over cos(phi_nominal)*R*phi for each endcap and for each disk (both rings)
0084   std::vector<std::vector<double> > sumOverCosThetaY(2, std::vector<double>(9, 0.));
0085
0086   // double sum over sin or cos(phi_nominal)*phi, for each endcap
0087   std::vector<double> doubleSumOverSinThetaY(2, 0.);
0088   std::vector<double> doubleSumOverCosThetaY(2, 0.);
0089
0090   // double sum over sin or cos(phi_nominal)*phi*z, for each endcap
0091   std::vector<double> doubleSumOverSinThetaYZ(2, 0.);
0092   std::vector<double> doubleSumOverCosThetaYZ(2, 0.);
0093
0094   // sum over z values / sum over z^2 values
0095   std::vector<double> sumOverZ(2, 0.);
0096   std::vector<double> sumOverZSquared(2, 0.);
0097
0098   // now calculate the sums
0099   det = 0;
0100   ring = 0;
0101   beam = 0;
0102   disk = 0;
0103   do {
0104     if (ring == 1)
0105       continue;  //################################################################################################### BOUND TO RING6
0106     // current radius, depends on the ring
0107     const double radius = nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetR();
0108
0109     // residual in r*phi (in the formulas this corresponds to y_ik)
0110     const double residual = (measuredCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi() -
0111                              nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
0112                             radius;
0113
0114     sumOverY.at(det).at(disk) += residual;
0115     sumOverPhi.at(det).at(disk) += residual / radius;
0116     kSumOverPhi.at(det).at(beam) += residual / radius;
0117
0118     doubleSumOverY.at(det) += residual;
0119     doubleSumOverPhi.at(det) += residual / radius;
0120
0121     kSumOverPhiZ.at(det).at(beam) += diskZPositions.at(det).at(disk) * residual / radius;
0122
0123     doubleSumOverYZ.at(det) += diskZPositions.at(det).at(disk) * residual;
0124     doubleSumOverPhiZ.at(det) += diskZPositions.at(det).at(disk) * residual / radius;
0125
0126     sumOverSinThetaY.at(det).at(disk) += sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
0127     sumOverCosThetaY.at(det).at(disk) += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
0128
0129     doubleSumOverSinThetaY.at(det) += sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
0130     doubleSumOverCosThetaY.at(det) += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) * residual;
0131
0132     doubleSumOverSinThetaYZ.at(det) += sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
0133                                        diskZPositions.at(det).at(disk) * residual;
0134     doubleSumOverCosThetaYZ.at(det) += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
0135                                        diskZPositions.at(det).at(disk) * residual;
0136
0137   } while (globalLoop.TECLoop(det, ring, beam, disk));
0138
0139   // disk-wise sums
0140   for (disk = 0; disk < 9; ++disk) {
0141     sumOverZ.at(0) += diskZPositions.at(0).at(disk);
0142     sumOverZ.at(1) += diskZPositions.at(1).at(disk);
0143     sumOverZSquared.at(0) += pow(diskZPositions.at(0).at(disk), 2);
0144     sumOverZSquared.at(1) += pow(diskZPositions.at(1).at(disk), 2);
0145   }
0146
0147   // now we can calculate the parameters for both TECs simultaneously,
0148   // so they're all vectors( 2 ) for TEC+/- (global parameters), or dim 2*9 (disk parameters),
0149   // or dim 2*8 (beam parameters)
0150
0151   // define them..
0152
0153   // deltaPhi_0
0154   std::vector<double> deltaPhi0(2, 0.);
0155
0156   // deltaPhi_t
0157   std::vector<double> deltaPhiT(2, 0.);
0158
0159   // deltaPhi_k (k=0..8)
0160   std::vector<std::vector<double> > deltaPhiK(2, std::vector<double>(9, 0.));
0161
0162   // deltaX_0
0163   std::vector<double> deltaX0(2, 0.);
0164
0165   // deltaX_t
0166   std::vector<double> deltaXT(2, 0.);
0167
0168   // deltaX_k (k=0..8)
0169   std::vector<std::vector<double> > deltaXK(2, std::vector<double>(9, 0.));
0170
0171   // deltaY_0
0172   std::vector<double> deltaY0(2, 0.);
0173
0174   // deltaY_t
0175   std::vector<double> deltaYT(2, 0.);
0176
0177   // deltaY_k (k=0..8)
0178   std::vector<std::vector<double> > deltaYK(2, std::vector<double>(9, 0.));
0179
0180   // beam parameters: deltaTheta_A, deltaTheta_B (i=0..7)
0181   std::vector<std::vector<double> > deltaTA(2, std::vector<double>(8, 0.));
0182   std::vector<std::vector<double> > deltaTB(2, std::vector<double>(8, 0.));
0183
0184   // ..and fill them;
0185   // the additional divisors "/ 2." come from the fact that we average over both rings
0186   for (det = 0; det < 2; ++det) {  // TEC+/- loop
0187
0188     // deltaPhi_0
0189     // here we use the phi-sums to eliminate the radius
0190     deltaPhi0.at(det) =
0191         (sumOverZSquared.at(det) * doubleSumOverPhi.at(det) - sumOverZ.at(det) * doubleSumOverPhiZ.at(det)) /
0192         (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)));  // / 2.; // @@@@@@@
0193
0194     // deltaPhi_t
0195     // again use the phi-sums
0196     deltaPhiT.at(det) = endcapLength * (9. * doubleSumOverPhiZ.at(det) - sumOverZ.at(det) * doubleSumOverPhi.at(det)) /
0197                         (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)));  // / 2.; // @@@@@@@
0198
0199     // deltaPhi_k (k=0..8)
0200     // again use the phi-sums
0201     for (disk = 0; disk < 9; ++disk) {
0202       deltaPhiK.at(det).at(disk) = (-1. * diskZPositions.at(det).at(disk) * deltaPhiT.at(det) / endcapLength) -
0203                                    (deltaPhi0.at(det)) - sumOverPhi.at(det).at(disk) / 8.;  // / 2.; // @@@@@@@
0204     }
0205
0206     // deltaX_0
0207     deltaX0.at(det) = 2. *
0208                       (sumOverZ.at(det) * doubleSumOverSinThetaYZ.at(det) -
0209                        sumOverZSquared.at(det) * doubleSumOverSinThetaY.at(det)) /
0210                       (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)));  // / 2.; // @@@@@@@
0211
0212     // deltaX_t
0213     deltaXT.at(det) = 2. * endcapLength *
0214                       (sumOverZ.at(det) * doubleSumOverSinThetaY.at(det) - 9. * doubleSumOverSinThetaYZ.at(det)) /
0215                       (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)));  // / 2.; // @@@@@@@
0216
0217     // deltaX_k (k=0..8)
0218     for (disk = 0; disk < 9; ++disk) {
0219       deltaXK.at(det).at(disk) = (-1. * diskZPositions.at(det).at(disk) * deltaXT.at(det) / endcapLength) -
0220                                  (deltaX0.at(det)) + 2. * sumOverSinThetaY.at(det).at(disk) / 8.;  // / 2.; // @@@@@@@
0221     }
0222
0223     // deltaY_0
0224     deltaY0.at(det) = 2. *
0225                       (sumOverZSquared.at(det) * doubleSumOverCosThetaY.at(det) -
0226                        sumOverZ.at(det) * doubleSumOverCosThetaYZ.at(det)) /
0227                       (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)));  // / 2.; // @@@@@@@
0228
0229     // deltaY_t
0230     deltaYT.at(det) = 2. * endcapLength *
0231                       (9. * doubleSumOverCosThetaYZ.at(det) - sumOverZ.at(det) * doubleSumOverCosThetaY.at(det)) /
0232                       (8. * (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)));  // / 2.; // @@@@@@@
0233
0234     // deltaY_k (k=0..8)
0235     for (disk = 0; disk < 9; ++disk) {
0236       deltaYK.at(det).at(disk) = (-1. * diskZPositions.at(det).at(disk) * deltaYT.at(det) / endcapLength) -
0237                                  (deltaY0.at(det)) - 2. * sumOverCosThetaY.at(det).at(disk) / 8.;  // / 2.; // @@@@@@@
0238     }
0239
0240     // the beam parameters deltaTA & deltaTB
0241     for (beam = 0; beam < 8; ++beam) {
0242       deltaTA.at(det).at(beam) =
0243           deltaPhi0.at(det) -
0244           (kSumOverPhi.at(det).at(beam) * sumOverZSquared.at(det) - kSumOverPhiZ.at(det).at(beam) * sumOverZ.at(det)) /
0245               (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)) +
0246           (cos(beamPhiPositions.at(det).at(beam)) * deltaY0.at(det) -
0247            sin(beamPhiPositions.at(det).at(beam)) * deltaX0.at(det)) /
0248               radius.at(0);  // for ring 4..
0249       // + ( cos( beamPhiPositions.at( det ).at( beam ) ) * deltaY0.at( det ) - sin( beamPhiPositions.at( det ).at( beam ) ) * deltaX0.at( det ) ) / radius.at( 1 ); // ...and ring 6
0250
0251       deltaTB.at(det).at(beam) =
0252           -1. * deltaPhiT.at(det) - deltaPhi0.at(det) -
0253           (kSumOverPhi.at(det).at(beam) * sumOverZ.at(det) - 9. * kSumOverPhiZ.at(det).at(beam)) * endcapLength /
0254               (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)) -
0255           (kSumOverPhiZ.at(det).at(beam) * sumOverZ.at(det) - kSumOverPhi.at(det).at(beam) * sumOverZSquared.at(det)) /
0256               (pow(sumOverZ.at(det), 2) - 9. * sumOverZSquared.at(det)) +
0257           ((deltaXT.at(det) + deltaX0.at(det)) * sin(beamPhiPositions.at(det).at(beam)) -
0258            (deltaYT.at(det) + deltaY0.at(det)) * cos(beamPhiPositions.at(det).at(beam))) /
0259               radius.at(0);  // for ring4..
0260       // + ( ( deltaXT.at( det ) + deltaX0.at( det ) ) * sin( beamPhiPositions.at( det ).at( beam ) ) - ( deltaYT.at( det ) + deltaY0.at( det ) ) * cos( beamPhiPositions.at( det ).at( beam ) ) )
0261       // / radius.at( 1 ); // ..and ring6
0262     }
0263   }
0264
0265   // fill the result
0266   LASEndcapAlignmentParameterSet theResult;
0267
0268   // for the moment we fill only the values, not the errors
0269
0270   // disk parameters
0271   for (det = 0; det < 2; ++det) {
0272     for (disk = 0; disk < 9; ++disk) {
0273       // the rotation parameters: deltaPhi_k
0274       theResult.GetDiskParameter(det, disk, 0).first = deltaPhiK.at(det).at(disk);
0275       // the x offsets: deltaX_k
0276       theResult.GetDiskParameter(det, disk, 1).first = deltaXK.at(det).at(disk);
0277       // the y offsets: deltaY_k
0278       theResult.GetDiskParameter(det, disk, 2).first = deltaYK.at(det).at(disk);
0279     }
0280   }
0281
0282   // global parameters
0283   for (int det = 0; det < 2; ++det) {
0284     theResult.GetGlobalParameter(det, 0).first = deltaPhi0.at(det);
0285     theResult.GetGlobalParameter(det, 1).first = deltaPhiT.at(det);
0286     theResult.GetGlobalParameter(det, 2).first = deltaX0.at(det);
0287     theResult.GetGlobalParameter(det, 3).first = deltaXT.at(det);
0288     theResult.GetGlobalParameter(det, 4).first = deltaY0.at(det);
0289     theResult.GetGlobalParameter(det, 5).first = deltaYT.at(det);
0290   }
0291
0292   // beam parameters
0293   for (int det = 0; det < 2; ++det) {
0294     for (int beam = 0; beam < 8; ++beam) {
0295       theResult.GetBeamParameter(det, 1 /*R6*/, beam, 0).first =
0296           deltaTA.at(det).at(beam);  ////////////////////////////////////////////
0297       theResult.GetBeamParameter(det, 1 /*R6*/, beam, 1).first =
0298           deltaTB.at(det).at(beam);  ////////////////////////////////////////////
0299     }
0300   }
0301
0302   std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Done." << std::endl;
0303
0304   return (theResult);
0305 }
0306
0307 ///
0308 /// for a given set of endcap alignment parameters "endcapParameters",
0309 /// this function returns the global phi offset from nominalPosition
0310 /// for a module specified by (det,ring,beam,disk)
0311 ///
0312 double LASEndcapAlgorithm::GetAlignmentParameterCorrection(int det,
0313                                                            int ring,
0314                                                            int beam,
0315                                                            int disk,
0316                                                            LASGlobalData<LASCoordinateSet>& nominalCoordinates,
0317                                                            LASEndcapAlignmentParameterSet& endcapParameters) {
0318   // ring dependent radius, to be softcoded...
0319   const double radius = ring == 0 ? 564. : 840.;
0320   const double endcapLength = 1345.;  // mm
0321
0322   // the correction to phi from the endcap algorithm;
0323   // it is defined such that the correction is to be subtracted
0324   double phiCorrection = 0.;
0325
0326   // plain disk phi
0327   phiCorrection += endcapParameters.GetDiskParameter(det, disk, 0).first;
0328
0329   // phi component from x deviation
0330   phiCorrection -= sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
0331                    endcapParameters.GetDiskParameter(det, disk, 1).first;
0332
0333   // phi component from y deviation
0334   phiCorrection += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
0335                    endcapParameters.GetDiskParameter(det, disk, 2).first;
0336
0337   // phi correction from global phi
0338   phiCorrection += endcapParameters.GetGlobalParameter(det, 0).first;
0339
0340   // correction from global x deviation
0341   phiCorrection -= sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
0342                    endcapParameters.GetGlobalParameter(det, 2).first;
0343
0344   // correction from global y deviation
0345   phiCorrection += cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) / radius *
0346                    endcapParameters.GetGlobalParameter(det, 4).first;
0347
0348   // correction from global torsion
0349   phiCorrection += nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength *
0350                    endcapParameters.GetGlobalParameter(det, 1).first;
0351
0352   // correction from global x shear
0353   phiCorrection -= nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength / radius *
0354                    sin(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
0355                    endcapParameters.GetGlobalParameter(det, 3).first;
0356
0357   // correction from global y shear
0358   phiCorrection += nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength / radius *
0359                    cos(nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetPhi()) *
0360                    endcapParameters.GetGlobalParameter(det, 5).first;
0361
0362   // correction from beam parameters
0363   phiCorrection += (nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength - 1.) *
0364                    endcapParameters.GetBeamParameter(det, 1, beam, 0).first;
0365   phiCorrection += nominalCoordinates.GetTECEntry(det, ring, beam, disk).GetZ() / endcapLength *
0366                    endcapParameters.GetBeamParameter(det, 1, beam, 1).first;
0367
0368   return phiCorrection;
0369 }