LaserAlignment/src/LASAlignmentTubeAlgorithm.cc

0001
0002 #include "Alignment/LaserAlignment/interface/LASAlignmentTubeAlgorithm.h"
0003
0004 ///
0005 ///
0006 ///
0007 LASAlignmentTubeAlgorithm::LASAlignmentTubeAlgorithm() {}
0008
0009 ///
0010 ///
0011 ///
0012 LASBarrelAlignmentParameterSet LASAlignmentTubeAlgorithm::CalculateParameters(
0013     LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
0014   std::cout << " [LASAlignmentTubeAlgorithm::CalculateParameters] -- Starting." << std::endl;
0015
0016   // for debugging only
0017   //######################################################################################
0018   //ReadMisalignmentFromFile( "misalign-var.txt", measuredCoordinates, nominalCoordinates );
0019   //######################################################################################
0020
0021   // loop object
0022   LASGlobalLoop globalLoop;
0023   int det, beam, disk, pos;
0024
0025   // phi positions of the AT beams in rad
0026   const double phiPositions[8] = {0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};
0027   std::vector<double> beamPhiPositions(8, 0.);
0028   for (beam = 0; beam < 8; ++beam)
0029     beamPhiPositions.at(beam) = phiPositions[beam];
0030
0031   // the radii of the alignment tube beams for each halfbarrel.
0032   // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
0033   // in TIB/TOB modules these radii differ from the beam radius..
0034   // ..due to the radial offsets (after the semitransparent mirrors)
0035   const double radii[6] = {564., 564., 514., 514., 600., 600.};
0036   std::vector<double> beamRadii(6, 0.);
0037   for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel)
0038     beamRadii.at(aHalfbarrel) = radii[aHalfbarrel];
0039
0040   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
0041   // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
0042   std::vector<std::vector<std::vector<double> > > endFaceZPositions(
0043       4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
0044   endFaceZPositions.at(0).at(0).at(0) = 1322.5;   // TEC+, *, disk1 ///
0045   endFaceZPositions.at(0).at(0).at(1) = 2667.5;   // TEC+, *, disk9 /// SIDE INFORMATION
0046   endFaceZPositions.at(1).at(0).at(0) = -2667.5;  // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
0047   endFaceZPositions.at(1).at(0).at(1) = -1322.5;  // TEC-, *, disk1 ///
0048   endFaceZPositions.at(2).at(1).at(0) = -700.;    // TIB,  -, outer
0049   endFaceZPositions.at(2).at(1).at(1) = -300.;    // TIB,  -, inner
0050   endFaceZPositions.at(2).at(0).at(0) = 300.;     // TIB,  +, inner
0051   endFaceZPositions.at(2).at(0).at(1) = 700.;     // TIB,  +, outer
0052   endFaceZPositions.at(3).at(1).at(0) = -1090.;   // TOB,  -, outer
0053   endFaceZPositions.at(3).at(1).at(1) = -300.;    // TOB,  -, inner
0054   endFaceZPositions.at(3).at(0).at(0) = 300.;     // TOB,  +, inner
0055   endFaceZPositions.at(3).at(0).at(1) = 1090.;    // TOB,  +, outer
0056
0057   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
0058   double detReducedZ[2] = {0., 0.};
0059   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
0060   double beamReducedZ[2] = {0., 0.};
0061
0062   // the z positions of the virtual planes at which the beam parameters are measured
0063   std::vector<double> disk9EndFaceZPositions(2, 0.);
0064   disk9EndFaceZPositions.at(0) = -2667.5;  // TEC- disk9
0065   disk9EndFaceZPositions.at(1) = 2667.5;   // TEC+ disk9
0066
0067   // define sums over measured values to "simplify" the beam parameter formulas
0068
0069   // all these have 6 entries, one for each halfbarrel (TEC+,TEC-,TIB+,TIB-,TOB+,TOB-)
0070   std::vector<double> sumOverPhiZPrime(6, 0.);
0071   std::vector<double> sumOverPhiZPrimePrime(6, 0.);
0072   std::vector<double> sumOverPhiZPrimeSinTheta(6, 0.);
0073   std::vector<double> sumOverPhiZPrimePrimeSinTheta(6, 0.);
0074   std::vector<double> sumOverPhiZPrimeCosTheta(6, 0.);
0075   std::vector<double> sumOverPhiZPrimePrimeCosTheta(6, 0.);
0076
0077   // these have 8 entries, one for each beam
0078   std::vector<double> sumOverPhiZTPrime(8, 0.);
0079   std::vector<double> sumOverPhiZTPrimePrime(8, 0.);
0080
0081   // define sums over nominal values
0082
0083   // all these have 6 entries, one for each halfbarrel (TEC+,TEC-,TIB+,TIB-,TOB+,TOB-)
0084   std::vector<double> sumOverZPrimeSquared(6, 0.);
0085   std::vector<double> sumOverZPrimePrimeSquared(6, 0.);
0086   std::vector<double> sumOverZPrimeZPrimePrime(6, 0.);
0087   std::vector<double> sumOverZPrimeZTPrime(6, 0.);
0088   std::vector<double> sumOverZPrimeZTPrimePrime(6, 0.);
0089   std::vector<double> sumOverZPrimePrimeZTPrime(6, 0.);
0090   std::vector<double> sumOverZPrimePrimeZTPrimePrime(6, 0.);
0091
0092   // all these are scalars
0093   double sumOverZTPrimeSquared = 0.;
0094   double sumOverZTPrimePrimeSquared = 0.;
0095   double sumOverZTPrimeZTPrimePrime = 0.;
0096
0097   // now calculate them for TIBTOB
0098   det = 2;
0099   beam = 0;
0100   pos = 0;
0101   do {
0102     // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
0103     const int theSide = pos < 3 ? 0 : 1;
0104
0105     // define the halfbarrel number from det/side
0106     const int halfbarrel = det == 2 ? det + theSide : det + 1 + theSide;  // TIB:TOB
0107
0108     // this is the path the beam has to travel radially after being reflected
0109     // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
0110     const double radialOffset = det == 2 ? 50. : 36.;
0111
0112     // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
0113     detReducedZ[0] = measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() -
0114                      endFaceZPositions.at(det).at(theSide).at(0);  // = zPrime
0115     detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0116     detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
0117                      measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ();  // = zPrimePrime
0118     detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0119
0120     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
0121     beamReducedZ[0] = (measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset) -
0122                       disk9EndFaceZPositions.at(0);  // = ZTPrime
0123     beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0124     beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
0125                       (measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset);  // ZTPrimePrime
0126     beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0127
0128     // residual in phi (in the formulas this corresponds to y_ik/R)
0129     const double phiResidual = measuredCoordinates.GetTIBTOBEntry(det, beam, pos).GetPhi() -
0130                                nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetPhi();
0131
0132     // sums over measured values
0133     sumOverPhiZPrime.at(halfbarrel) += phiResidual * detReducedZ[0];
0134     sumOverPhiZPrimePrime.at(halfbarrel) += phiResidual * detReducedZ[1];
0135     sumOverPhiZPrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * sin(beamPhiPositions.at(beam));
0136     sumOverPhiZPrimePrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * sin(beamPhiPositions.at(beam));
0137     sumOverPhiZPrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * cos(beamPhiPositions.at(beam));
0138     sumOverPhiZPrimePrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * cos(beamPhiPositions.at(beam));
0139
0140     sumOverPhiZTPrime.at(beam) += phiResidual * beamReducedZ[0];  // note the index change here..
0141     sumOverPhiZTPrimePrime.at(beam) += phiResidual * beamReducedZ[1];
0142
0143     // sums over nominal values
0144     sumOverZPrimeSquared.at(halfbarrel) += pow(detReducedZ[0], 2) / 8.;  // these are defined beam-wise, so: / 8.
0145     sumOverZPrimePrimeSquared.at(halfbarrel) += pow(detReducedZ[1], 2) / 8.;
0146     sumOverZPrimeZPrimePrime.at(halfbarrel) += detReducedZ[0] * detReducedZ[1] / 8.;
0147     sumOverZPrimeZTPrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[0] / 8.;
0148     sumOverZPrimeZTPrimePrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[1] / 8.;
0149     sumOverZPrimePrimeZTPrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[0] / 8.;
0150     sumOverZPrimePrimeZTPrimePrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[1] / 8.;
0151
0152     sumOverZTPrimeSquared += pow(beamReducedZ[0], 2) / 8.;
0153     sumOverZTPrimePrimeSquared += pow(beamReducedZ[1], 2) / 8.;
0154     sumOverZTPrimeZTPrimePrime += beamReducedZ[0] * beamReducedZ[1] / 8.;
0155
0156   } while (globalLoop.TIBTOBLoop(det, beam, pos));
0157
0158   // now for TEC2TEC
0159   det = 0;
0160   beam = 0;
0161   disk = 0;
0162   do {
0163     // for the tec, the halfbarrel numbers are equal to the det numbers...
0164     const int halfbarrel = det;
0165
0166     // ...so there's no side distinction for the TEC
0167     const int theSide = 0;
0168
0169     // also, there's no radial offset for the TEC
0170     const double radialOffset = 0.;
0171
0172     // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
0173     detReducedZ[0] = measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() -
0174                      endFaceZPositions.at(det).at(theSide).at(0);  // = zPrime
0175     detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0176     detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
0177                      measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ();  // = zPrimePrime
0178     detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0179
0180     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
0181     beamReducedZ[0] = (measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset) -
0182                       disk9EndFaceZPositions.at(0);  // = ZTPrime
0183     beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0184     beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
0185                       (measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset);  // ZTPrimePrime
0186     beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0187
0188     // residual in phi (in the formulas this corresponds to y_ik/R)
0189     const double phiResidual = measuredCoordinates.GetTEC2TECEntry(det, beam, disk).GetPhi() -
0190                                nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetPhi();
0191
0192     // sums over measured values
0193     sumOverPhiZPrime.at(halfbarrel) += phiResidual * detReducedZ[0];
0194     sumOverPhiZPrimePrime.at(halfbarrel) += phiResidual * detReducedZ[1];
0195     sumOverPhiZPrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * sin(beamPhiPositions.at(beam));
0196     sumOverPhiZPrimePrimeSinTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * sin(beamPhiPositions.at(beam));
0197     sumOverPhiZPrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[0] * cos(beamPhiPositions.at(beam));
0198     sumOverPhiZPrimePrimeCosTheta.at(halfbarrel) += phiResidual * detReducedZ[1] * cos(beamPhiPositions.at(beam));
0199
0200     sumOverPhiZTPrime.at(beam) += phiResidual * beamReducedZ[0];  // note the index change here..
0201     sumOverPhiZTPrimePrime.at(beam) += phiResidual * beamReducedZ[1];
0202
0203     // sums over nominal values
0204     sumOverZPrimeSquared.at(halfbarrel) += pow(detReducedZ[0], 2) / 8.;  // these are defined beam-wise, so: / 8.
0205     sumOverZPrimePrimeSquared.at(halfbarrel) += pow(detReducedZ[1], 2) / 8.;
0206     sumOverZPrimeZPrimePrime.at(halfbarrel) += detReducedZ[0] * detReducedZ[1] / 8.;
0207     sumOverZPrimeZTPrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[0] / 8.;
0208     sumOverZPrimeZTPrimePrime.at(halfbarrel) += detReducedZ[0] * beamReducedZ[1] / 8.;
0209     sumOverZPrimePrimeZTPrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[0] / 8.;
0210     sumOverZPrimePrimeZTPrimePrime.at(halfbarrel) += detReducedZ[1] * beamReducedZ[1] / 8.;
0211
0212     sumOverZTPrimeSquared += pow(beamReducedZ[0], 2) / 8.;
0213     sumOverZTPrimePrimeSquared += pow(beamReducedZ[1], 2) / 8.;
0214     sumOverZTPrimeZTPrimePrime += beamReducedZ[0] * beamReducedZ[1] / 8.;
0215
0216   } while (globalLoop.TEC2TECLoop(det, beam, disk));
0217
0218   // more "simplification" terms...
0219   // these here are functions of theta and can be calculated directly
0220   double sumOverSinTheta = 0.;
0221   double sumOverCosTheta = 0.;
0222   double sumOverSinThetaSquared = 0.;
0223   double sumOverCosThetaSquared = 0.;
0224   double sumOverCosThetaSinTheta = 0.;
0225   double mixedTrigonometricTerm = 0.;
0226
0227   for (beam = 0; beam < 8; ++beam) {
0228     sumOverSinTheta += sin(beamPhiPositions.at(beam));
0229     sumOverCosTheta += cos(beamPhiPositions.at(beam));
0230     sumOverSinThetaSquared += pow(sin(beamPhiPositions.at(beam)), 2);
0231     sumOverCosThetaSquared += pow(cos(beamPhiPositions.at(beam)), 2);
0232     sumOverCosThetaSinTheta += cos(beamPhiPositions.at(beam)) * sin(beamPhiPositions.at(beam));
0233   }
0234
0235   mixedTrigonometricTerm = 8. * (sumOverCosThetaSquared * sumOverSinThetaSquared - pow(sumOverCosThetaSinTheta, 2)) -
0236                            pow(sumOverCosTheta, 2) * sumOverSinThetaSquared -
0237                            pow(sumOverSinTheta, 2) * sumOverCosThetaSquared +
0238                            2. * sumOverCosTheta * sumOverSinTheta * sumOverCosThetaSinTheta;
0239
0240   // even more shortcuts before we can calculate the parameters
0241   double beamDenominator =
0242       (pow(sumOverZTPrimeZTPrimePrime, 2) - sumOverZTPrimeSquared * sumOverZTPrimePrimeSquared) * beamRadii.at(0);
0243   std::vector<double> alignmentDenominator(6, 0.);
0244   std::vector<double> termA(6, 0.), termB(6, 0.), termC(6, 0.), termD(6, 0.);
0245   for (unsigned int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel) {
0246     alignmentDenominator.at(aHalfbarrel) =
0247         (pow(sumOverZPrimeZPrimePrime.at(aHalfbarrel), 2) -
0248          sumOverZPrimeSquared.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel)) *
0249         mixedTrigonometricTerm;
0250     termA.at(aHalfbarrel) = sumOverZPrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
0251                             sumOverZPrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeSquared;
0252     termB.at(aHalfbarrel) = sumOverZPrimePrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
0253                             sumOverZPrimePrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeSquared;
0254     termC.at(aHalfbarrel) = sumOverZPrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
0255                             sumOverZPrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimePrimeSquared;
0256     termD.at(aHalfbarrel) = sumOverZPrimePrimeZTPrimePrime.at(aHalfbarrel) * sumOverZTPrimeZTPrimePrime -
0257                             sumOverZPrimePrimeZTPrime.at(aHalfbarrel) * sumOverZTPrimePrimeSquared;
0258   }
0259
0260   // have eight alignment tube beams..
0261   const int numberOfBeams = 8;
0262
0263   // that's all for preparation, now it gets ugly:
0264   // calculate the alignment parameters
0265   LASBarrelAlignmentParameterSet theResult;
0266
0267   // can do this in one go for all halfbarrels
0268   for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel) {
0269     // no errors yet
0270
0271     // rotation angles of the lower z endfaces (in rad)
0272     theResult.GetParameter(aHalfbarrel, 0, 0).first =
0273         (sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
0274              sumOverSinThetaSquared -
0275          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
0276              sumOverSinThetaSquared -
0277          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0278              sumOverSinThetaSquared +
0279          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0280              sumOverSinThetaSquared +
0281          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
0282              sumOverSinTheta -
0283          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) *
0284              sumOverCosThetaSinTheta * sumOverSinTheta -
0285          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
0286              sumOverSinTheta +
0287          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
0288              sumOverSinTheta -
0289          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * pow(sumOverCosThetaSinTheta, 2) +
0290          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) *
0291              pow(sumOverCosThetaSinTheta, 2) +
0292          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0293              sumOverCosThetaSinTheta -
0294          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0295              sumOverCosThetaSinTheta) /
0296         alignmentDenominator.at(aHalfbarrel);
0297
0298     // rotation angles of the upper z endfaces (in rad)
0299     theResult.GetParameter(aHalfbarrel, 1, 0).first =
0300         (sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
0301              sumOverSinThetaSquared -
0302          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
0303              sumOverSinThetaSquared -
0304          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0305              sumOverSinThetaSquared +
0306          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0307              sumOverSinThetaSquared +
0308          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) *
0309              sumOverCosThetaSinTheta * sumOverSinTheta -
0310          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) *
0311              sumOverCosThetaSinTheta * sumOverSinTheta -
0312          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) *
0313              sumOverCosThetaSquared * sumOverSinTheta +
0314          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
0315              sumOverSinTheta -
0316          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
0317              sumOverCosThetaSinTheta +
0318          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSinTheta *
0319              sumOverCosThetaSinTheta +
0320          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0321              sumOverCosThetaSinTheta -
0322          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0323              sumOverCosThetaSinTheta) /
0324         alignmentDenominator.at(aHalfbarrel);
0325
0326     // x deviations of the lower z endfaces (in mm)
0327     theResult.GetParameter(aHalfbarrel, 0, 1).first =
0328         -1. *
0329         (sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
0330              sumOverSinTheta -
0331          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
0332              sumOverSinTheta -
0333          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0334              sumOverSinTheta +
0335          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0336              sumOverSinTheta -
0337          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0338              sumOverCosThetaSinTheta +
0339          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0340              sumOverCosThetaSinTheta +
0341          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0342              sumOverCosThetaSinTheta -
0343          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
0344              sumOverCosThetaSinTheta -
0345          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0346              sumOverCosThetaSquared +
0347          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
0348              sumOverCosThetaSquared +
0349          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0350              sumOverCosTheta -
0351          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0352              sumOverCosTheta) /
0353         alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
0354
0355     // x deviations of the upper z endfaces (in mm)
0356     theResult.GetParameter(aHalfbarrel, 1, 1).first =
0357         -1. *
0358         (sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSquared *
0359              sumOverSinTheta -
0360          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSquared *
0361              sumOverSinTheta -
0362          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0363              sumOverSinTheta +
0364          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0365              sumOverSinTheta -
0366          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0367              sumOverCosThetaSinTheta +
0368          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0369              sumOverCosThetaSinTheta +
0370          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0371              sumOverCosThetaSinTheta -
0372          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
0373              sumOverCosThetaSinTheta -
0374          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0375              sumOverCosThetaSquared +
0376          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
0377              sumOverCosThetaSquared +
0378          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0379              sumOverCosTheta -
0380          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0381              sumOverCosTheta) /
0382         alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
0383
0384     // y deviations of the lower z endfaces (in mm)
0385     theResult.GetParameter(aHalfbarrel, 0, 2).first =
0386         (sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0387              sumOverSinThetaSquared -
0388          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0389              sumOverSinThetaSquared -
0390          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0391              sumOverSinThetaSquared +
0392          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
0393              sumOverSinThetaSquared +
0394          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverSinTheta *
0395              sumOverSinTheta -
0396          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverSinTheta *
0397              sumOverSinTheta -
0398          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
0399              sumOverSinTheta +
0400          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosThetaSinTheta *
0401              sumOverSinTheta -
0402          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0403              sumOverSinTheta +
0404          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0405              sumOverSinTheta +
0406          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0407              sumOverCosThetaSinTheta -
0408          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeSquared.at(aHalfbarrel) * numberOfBeams *
0409              sumOverCosThetaSinTheta) /
0410         alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
0411
0412     // y deviations of the upper z endfaces (in mm)
0413     theResult.GetParameter(aHalfbarrel, 1, 2).first =
0414         (sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0415              sumOverSinThetaSquared -
0416          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0417              sumOverSinThetaSquared -
0418          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0419              sumOverSinThetaSquared +
0420          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
0421              sumOverSinThetaSquared +
0422          sumOverPhiZPrimePrimeCosTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverSinTheta *
0423              sumOverSinTheta -
0424          sumOverPhiZPrimeCosTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverSinTheta *
0425              sumOverSinTheta -
0426          sumOverPhiZPrimePrime.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosThetaSinTheta *
0427              sumOverSinTheta +
0428          sumOverPhiZPrime.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosThetaSinTheta *
0429              sumOverSinTheta -
0430          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * sumOverCosTheta *
0431              sumOverSinTheta +
0432          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * sumOverCosTheta *
0433              sumOverSinTheta +
0434          sumOverPhiZPrimePrimeSinTheta.at(aHalfbarrel) * sumOverZPrimeZPrimePrime.at(aHalfbarrel) * numberOfBeams *
0435              sumOverCosThetaSinTheta -
0436          sumOverPhiZPrimeSinTheta.at(aHalfbarrel) * sumOverZPrimePrimeSquared.at(aHalfbarrel) * numberOfBeams *
0437              sumOverCosThetaSinTheta) /
0438         alignmentDenominator.at(aHalfbarrel) * beamRadii.at(aHalfbarrel);
0439   }
0440
0441   // another loop is needed here to calculate some terms for the beam parameters
0442   double vsumA = 0., vsumB = 0., vsumC = 0., vsumD = 0., vsumE = 0., vsumF = 0.;
0443   for (unsigned int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel) {
0444     vsumA += theResult.GetParameter(aHalfbarrel, 1, 2).first * termA.at(aHalfbarrel) +
0445              theResult.GetParameter(aHalfbarrel, 0, 2).first * termB.at(aHalfbarrel);
0446     vsumB += theResult.GetParameter(aHalfbarrel, 1, 1).first * termA.at(aHalfbarrel) +
0447              theResult.GetParameter(aHalfbarrel, 0, 1).first * termB.at(aHalfbarrel);
0448     vsumC += beamRadii.at(aHalfbarrel) * (theResult.GetParameter(aHalfbarrel, 1, 0).first * termA.at(aHalfbarrel) +
0449                                           theResult.GetParameter(aHalfbarrel, 0, 0).first * termB.at(aHalfbarrel));
0450     vsumD += theResult.GetParameter(aHalfbarrel, 1, 2).first * termC.at(aHalfbarrel) +
0451              theResult.GetParameter(aHalfbarrel, 0, 2).first * termD.at(aHalfbarrel);
0452     vsumE += theResult.GetParameter(aHalfbarrel, 1, 1).first * termC.at(aHalfbarrel) +
0453              theResult.GetParameter(aHalfbarrel, 0, 1).first * termD.at(aHalfbarrel);
0454     vsumF += beamRadii.at(aHalfbarrel) * (theResult.GetParameter(aHalfbarrel, 1, 0).first * termC.at(aHalfbarrel) +
0455                                           theResult.GetParameter(aHalfbarrel, 0, 0).first * termD.at(aHalfbarrel));
0456   }
0457
0458   // calculate the beam parameters
0459   for (unsigned int beam = 0; beam < 8; ++beam) {
0460     // parameter A, defined at lower z
0461     theResult.GetBeamParameter(beam, 0).first =
0462         (cos(beamPhiPositions.at(beam)) * vsumA - sin(beamPhiPositions.at(beam)) * vsumB - vsumC +
0463          sumOverPhiZTPrime.at(beam) * sumOverZTPrimeZTPrimePrime -
0464          sumOverPhiZTPrimePrime.at(beam) * sumOverZTPrimeSquared) /
0465         beamDenominator;
0466
0467     ///////////////////////////////////////////////////////////////////////////////////////////////////
0468     std::cout << "BBBBBBBB: " << cos(beamPhiPositions.at(beam)) * vsumA << "  "
0469               << -1. * sin(beamPhiPositions.at(beam)) * vsumB << "  " << -1. * vsumC << "  "
0470               << sumOverPhiZTPrime.at(beam) * sumOverZTPrimeZTPrimePrime -
0471                      sumOverPhiZTPrimePrime.at(beam) * sumOverZTPrimeSquared
0472               << "  " << beamDenominator << std::endl;
0473     ///////////////////////////////////////////////////////////////////////////////////////////////////
0474
0475     // parameter B, defined at upper z
0476     theResult.GetBeamParameter(beam, 1).first =
0477         (cos(beamPhiPositions.at(beam)) * vsumD - sin(beamPhiPositions.at(beam)) * vsumE - vsumF +
0478          sumOverPhiZTPrimePrime.at(beam) * sumOverZTPrimeZTPrimePrime -
0479          sumOverPhiZTPrime.at(beam) * sumOverZTPrimePrimeSquared) /
0480         beamDenominator;
0481   }
0482
0483   return theResult;
0484 }
0485
0486 ///
0487 /// get global phi correction from alignment parameters
0488 /// for an alignment tube module in TIB/TOB
0489 ///
0490 double LASAlignmentTubeAlgorithm::GetTIBTOBAlignmentParameterCorrection(
0491     int det,
0492     int beam,
0493     int pos,
0494     LASGlobalData<LASCoordinateSet>& nominalCoordinates,
0495     LASBarrelAlignmentParameterSet& alignmentParameters) {
0496   // INITIALIZATION;
0497   // ALL THIS IS DUPLICATED FOR THE MOMENT, SHOULD FINALLY BE CALCULATED ONLY ONCE
0498   // AND HARD CODED NUMBERS SHOULD CENTRALLY BE IMPORTED FROM src/LASConstants.h
0499
0500   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
0501   // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
0502   std::vector<std::vector<std::vector<double> > > endFaceZPositions(
0503       4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
0504   endFaceZPositions.at(0).at(0).at(0) = 1322.5;   // TEC+, *, disk1 ///
0505   endFaceZPositions.at(0).at(0).at(1) = 2667.5;   // TEC+, *, disk9 /// SIDE INFORMATION
0506   endFaceZPositions.at(1).at(0).at(0) = -2667.5;  // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
0507   endFaceZPositions.at(1).at(0).at(1) = -1322.5;  // TEC-, *, disk1 ///
0508   endFaceZPositions.at(2).at(1).at(0) = -700.;    // TIB,  -, outer
0509   endFaceZPositions.at(2).at(1).at(1) = -300.;    // TIB,  -, inner
0510   endFaceZPositions.at(2).at(0).at(0) = 300.;     // TIB,  +, inner
0511   endFaceZPositions.at(2).at(0).at(1) = 700.;     // TIB,  +, outer
0512   endFaceZPositions.at(3).at(1).at(0) = -1090.;   // TOB,  -, outer
0513   endFaceZPositions.at(3).at(1).at(1) = -300.;    // TOB,  -, inner
0514   endFaceZPositions.at(3).at(0).at(0) = 300.;     // TOB,  +, inner
0515   endFaceZPositions.at(3).at(0).at(1) = 1090.;    // TOB,  +, outer
0516
0517   // the z positions of the virtual planes at which the beam parameters are measured
0518   std::vector<double> disk9EndFaceZPositions(2, 0.);
0519   disk9EndFaceZPositions.at(0) = -2667.5;  // TEC- disk9
0520   disk9EndFaceZPositions.at(1) = 2667.5;   // TEC+ disk9
0521
0522   // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
0523   const int theSide = pos < 3 ? 0 : 1;
0524
0525   // define the halfbarrel number from det/side
0526   const int halfbarrel = det == 2 ? det + theSide : det + 1 + theSide;  // TIB:TOB
0527
0528   // this is the path the beam has to travel radially after being reflected
0529   // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
0530   const double radialOffset = det == 2 ? 50. : 36.;
0531
0532   // phi positions of the AT beams in rad
0533   const double phiPositions[8] = {0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};
0534   std::vector<double> beamPhiPositions(8, 0.);
0535   for (unsigned int aBeam = 0; aBeam < 8; ++aBeam)
0536     beamPhiPositions.at(aBeam) = phiPositions[aBeam];
0537
0538   // the radii of the alignment tube beams for each halfbarrel.
0539   // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
0540   // in TIB/TOB modules these radii differ from the beam radius..
0541   // ..due to the radial offsets (after the semitransparent mirrors)
0542   const double radii[6] = {564., 564., 514., 514., 600., 600.};
0543   std::vector<double> beamRadii(6, 0.);
0544   for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel)
0545     beamRadii.at(aHalfbarrel) = radii[aHalfbarrel];
0546
0547   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
0548   double detReducedZ[2] = {0., 0.};
0549   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
0550   double beamReducedZ[2] = {0., 0.};
0551
0552   // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
0553   detReducedZ[0] = nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() -
0554                    endFaceZPositions.at(det).at(theSide).at(0);  // = zPrime
0555   detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0556   detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
0557                    nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ();  // = zPrimePrime
0558   detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0559
0560   // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
0561   beamReducedZ[0] = (nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset) -
0562                     disk9EndFaceZPositions.at(0);  // = ZTPrime
0563   beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0564   beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
0565                     (nominalCoordinates.GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset);  // ZTPrimePrime
0566   beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0567
0568   // the correction to phi from the endcap algorithm;
0569   // it is defined such that the correction is to be subtracted ///////////////////////////////// ???
0570   double phiCorrection = 0.;
0571
0572   // contribution from phi rotation of first end face
0573   phiCorrection += detReducedZ[1] * alignmentParameters.GetParameter(halfbarrel, 0, 0).first;
0574
0575   // contribution from phi rotation of second end face
0576   phiCorrection += detReducedZ[0] * alignmentParameters.GetParameter(halfbarrel, 1, 0).first;
0577
0578   // contribution from translation along x of first endface
0579   phiCorrection += detReducedZ[1] * sin(beamPhiPositions.at(beam)) *
0580                    alignmentParameters.GetParameter(halfbarrel, 0, 1).first / beamRadii.at(halfbarrel);
0581
0582   // contribution from translation along x of second endface
0583   phiCorrection += detReducedZ[0] * sin(beamPhiPositions.at(beam)) *
0584                    alignmentParameters.GetParameter(halfbarrel, 1, 1).first / beamRadii.at(halfbarrel);
0585
0586   // contribution from translation along y of first endface
0587   phiCorrection -= detReducedZ[1] * cos(beamPhiPositions.at(beam)) *
0588                    alignmentParameters.GetParameter(halfbarrel, 0, 2).first / beamRadii.at(halfbarrel);
0589
0590   // contribution from translation along y of second endface
0591   phiCorrection -= detReducedZ[0] * cos(beamPhiPositions.at(beam)) *
0592                    alignmentParameters.GetParameter(halfbarrel, 1, 2).first / beamRadii.at(halfbarrel);
0593
0594   // contribution from beam parameters;
0595   // originally, the contribution in meter is proportional to the radius of the beams: beamRadii.at( 0 )
0596   // the additional factor: beamRadii.at( halfbarrel ) converts from meter to radian on the module
0597   phiCorrection += beamReducedZ[1] * alignmentParameters.GetBeamParameter(beam, 0).first * beamRadii.at(0) /
0598                    beamRadii.at(halfbarrel);
0599   phiCorrection += beamReducedZ[0] * alignmentParameters.GetBeamParameter(beam, 1).first * beamRadii.at(0) /
0600                    beamRadii.at(halfbarrel);
0601
0602   return phiCorrection;
0603 }
0604
0605 ///
0606 /// get global phi correction from alignment parameters
0607 /// for an alignment tube module in TEC(AT)
0608 ///
0609 double LASAlignmentTubeAlgorithm::GetTEC2TECAlignmentParameterCorrection(
0610     int det,
0611     int beam,
0612     int disk,
0613     LASGlobalData<LASCoordinateSet>& nominalCoordinates,
0614     LASBarrelAlignmentParameterSet& alignmentParameters) {
0615   // INITIALIZATION;
0616   // ALL THIS IS DUPLICATED FOR THE MOMENT, SHOULD FINALLY BE CALCULATED ONLY ONCE
0617   // AND HARD CODED NUMBERS SHOULD CENTRALLY BE IMPORTED FROM src/LASConstants.h
0618
0619   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
0620   // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
0621   std::vector<std::vector<std::vector<double> > > endFaceZPositions(
0622       4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
0623   endFaceZPositions.at(0).at(0).at(0) = 1322.5;   // TEC+, *, disk1 ///
0624   endFaceZPositions.at(0).at(0).at(1) = 2667.5;   // TEC+, *, disk9 /// SIDE INFORMATION
0625   endFaceZPositions.at(1).at(0).at(0) = -2667.5;  // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
0626   endFaceZPositions.at(1).at(0).at(1) = -1322.5;  // TEC-, *, disk1 ///
0627   endFaceZPositions.at(2).at(1).at(0) = -700.;    // TIB,  -, outer
0628   endFaceZPositions.at(2).at(1).at(1) = -300.;    // TIB,  -, inner
0629   endFaceZPositions.at(2).at(0).at(0) = 300.;     // TIB,  +, inner
0630   endFaceZPositions.at(2).at(0).at(1) = 700.;     // TIB,  +, outer
0631   endFaceZPositions.at(3).at(1).at(0) = -1090.;   // TOB,  -, outer
0632   endFaceZPositions.at(3).at(1).at(1) = -300.;    // TOB,  -, inner
0633   endFaceZPositions.at(3).at(0).at(0) = 300.;     // TOB,  +, inner
0634   endFaceZPositions.at(3).at(0).at(1) = 1090.;    // TOB,  +, outer
0635
0636   // the z positions of the virtual planes at which the beam parameters are measured
0637   std::vector<double> disk9EndFaceZPositions(2, 0.);
0638   disk9EndFaceZPositions.at(0) = -2667.5;  // TEC- disk9
0639   disk9EndFaceZPositions.at(1) = 2667.5;   // TEC+ disk9
0640
0641   // for the tec, the halfbarrel numbers are equal to the det numbers...
0642   const int halfbarrel = det;
0643
0644   // ...so there's no side distinction for the TEC
0645   const int theSide = 0;
0646
0647   // also, there's no radial offset for the TEC
0648   const double radialOffset = 0.;
0649
0650   // phi positions of the AT beams in rad
0651   const double phiPositions[8] = {0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784};
0652   std::vector<double> beamPhiPositions(8, 0.);
0653   for (unsigned int aBeam = 0; aBeam < 8; ++aBeam)
0654     beamPhiPositions.at(aBeam) = phiPositions[aBeam];
0655
0656   // the radii of the alignment tube beams for each halfbarrel.
0657   // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
0658   // in TIB/TOB modules these radii differ from the beam radius..
0659   // ..due to the radial offsets (after the semitransparent mirrors)
0660   const double radii[6] = {564., 564., 514., 514., 600., 600.};
0661   std::vector<double> beamRadii(6, 0.);
0662   for (int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel)
0663     beamRadii.at(aHalfbarrel) = radii[aHalfbarrel];
0664
0665   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
0666   double detReducedZ[2] = {0., 0.};
0667   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
0668   double beamReducedZ[2] = {0., 0.};
0669
0670   // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
0671   detReducedZ[0] = nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() -
0672                    endFaceZPositions.at(det).at(theSide).at(0);  // = zPrime
0673   detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0674   detReducedZ[1] = endFaceZPositions.at(det).at(theSide).at(1) -
0675                    nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ();  // = zPrimePrime
0676   detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0677
0678   // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
0679   beamReducedZ[0] = (nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset) -
0680                     disk9EndFaceZPositions.at(0);  // = ZTPrime
0681   beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0682   beamReducedZ[1] = disk9EndFaceZPositions.at(1) -
0683                     (nominalCoordinates.GetTEC2TECEntry(det, beam, disk).GetZ() - radialOffset);  // ZTPrimePrime
0684   beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0685
0686   // the correction to phi from the endcap algorithm;
0687   // it is defined such that the correction is to be subtracted ///////////////////////////////// ???
0688   double phiCorrection = 0.;
0689
0690   // contribution from phi rotation of first end face
0691   phiCorrection += detReducedZ[1] * alignmentParameters.GetParameter(halfbarrel, 0, 0).first;
0692
0693   // contribution from phi rotation of second end face
0694   phiCorrection += detReducedZ[0] * alignmentParameters.GetParameter(halfbarrel, 1, 0).first;
0695
0696   // contribution from translation along x of first endface
0697   phiCorrection += detReducedZ[1] * sin(beamPhiPositions.at(beam)) *
0698                    alignmentParameters.GetParameter(halfbarrel, 0, 1).first / beamRadii.at(halfbarrel);
0699
0700   // contribution from translation along x of second endface
0701   phiCorrection += detReducedZ[0] * sin(beamPhiPositions.at(beam)) *
0702                    alignmentParameters.GetParameter(halfbarrel, 1, 1).first / beamRadii.at(halfbarrel);
0703
0704   // contribution from translation along y of first endface
0705   phiCorrection -= detReducedZ[1] * cos(beamPhiPositions.at(beam)) *
0706                    alignmentParameters.GetParameter(halfbarrel, 0, 2).first / beamRadii.at(halfbarrel);
0707
0708   // contribution from translation along y of second endface
0709   phiCorrection -= detReducedZ[0] * cos(beamPhiPositions.at(beam)) *
0710                    alignmentParameters.GetParameter(halfbarrel, 1, 2).first / beamRadii.at(halfbarrel);
0711
0712   // contribution from beam parameters;
0713   // originally, the contribution in meter is proportional to the radius of the beams: beamRadii.at( 0 )
0714   // the additional factor: beamRadii.at( halfbarrel ) converts from meter to radian on the module
0715   phiCorrection += beamReducedZ[1] * alignmentParameters.GetBeamParameter(beam, 0).first * beamRadii.at(0) /
0716                    beamRadii.at(halfbarrel);
0717   phiCorrection += beamReducedZ[0] * alignmentParameters.GetBeamParameter(beam, 1).first * beamRadii.at(0) /
0718                    beamRadii.at(halfbarrel);
0719
0720   return phiCorrection;
0721 }
0722
0723 ///
0724 /// allows to push in a simple simulated misalignment for quick internal testing purposes;
0725 /// overwrites LASGlobalData<LASCoordinateSet>& measuredCoordinates;
0726 /// call at beginning of LASBarrelAlgorithm::CalculateParameters method
0727 ///
0728 /// one line per module,
0729 /// format for TEC:              det ring beam disk phi phiErr
0730 /// format for TEC(at) & TIBTOB: det beam   z  "-1" phi phiErr
0731 ///
0732 void LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile(const char* filename,
0733                                                          LASGlobalData<LASCoordinateSet>& measuredCoordinates,
0734                                                          LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
0735   std::ifstream file(filename);
0736   if (file.bad()) {
0737     std::cerr << " [LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile] ** ERROR: cannot open file \"" << filename
0738               << "\"." << std::endl;
0739     return;
0740   }
0741
0742   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0743                "@@@@@@@@@@@"
0744             << std::endl;
0745   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0746                "@@@@@@@@@@@"
0747             << std::endl;
0748   std::cerr << " [LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile] ** WARNING: you are reading a fake measurement "
0749                "from a file!"
0750             << std::endl;
0751   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0752                "@@@@@@@@@@@"
0753             << std::endl;
0754   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0755                "@@@@@@@@@@@"
0756             << std::endl;
0757
0758   // the measured coordinates will finally be overwritten;
0759   // first, set them to the nominal values
0760   measuredCoordinates = nominalCoordinates;
0761
0762   // and put large errors on all values;
0763   {
0764     LASGlobalLoop moduleLoop;
0765     int det, ring, beam, disk, pos;
0766
0767     det = 0;
0768     ring = 0;
0769     beam = 0;
0770     disk = 0;
0771     do {
0772       measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhiError(1000.);
0773     } while (moduleLoop.TECLoop(det, ring, beam, disk));
0774
0775     det = 2;
0776     beam = 0;
0777     pos = 0;
0778     do {
0779       measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhiError(1000.);
0780     } while (moduleLoop.TIBTOBLoop(det, beam, pos));
0781
0782     det = 0;
0783     beam = 0;
0784     disk = 0;
0785     do {
0786       measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhiError(1000.);
0787     } while (moduleLoop.TEC2TECLoop(det, beam, disk));
0788   }
0789
0790   // buffers for read-in
0791   int det, beam, z, ring;
0792   double phi, phiError;
0793
0794   while (!file.eof()) {
0795     file >> det;
0796     if (file.eof())
0797       break;  // do not read the last line twice, do not fill trash if file empty
0798
0799     file >> beam;
0800     file >> z;
0801     file >> ring;
0802     file >> phi;
0803     file >> phiError;
0804
0805     if (det > 1) {  // TIB/TOB
0806       measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhi(phi);
0807       measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhiError(phiError);
0808     } else {            // TEC or TEC(at)
0809       if (ring > -1) {  // TEC
0810         measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhi(phi);
0811         measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhiError(phiError);
0812       } else {  // TEC(at)
0813         measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhi(phi);
0814         measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhiError(phiError);
0815       }
0816     }
0817   }
0818
0819   file.close();
0820 }