CommonAlignmentParametrization/src/TwoBowedSurfacesAlignmentParameters.cc

0001 /** \file TwoBowedSurfacesAlignmentParameters.cc
0002  *
0003  *  Version    : $Revision: 1.2 $
0004  *  last update: $Date: 2010/11/23 13:50:50 $
0005  *  by         : $Author: flucke $
0006  */
0007
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0012
0013 #include "Alignment/CommonAlignment/interface/Alignable.h"
0014 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
0015 #include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
0016 #include "Alignment/CommonAlignmentParametrization/interface/KarimakiAlignmentDerivatives.h"
0017 #include "CondFormats/Alignment/interface/Definitions.h"
0018 //#include "DataFormats/Alignment/interface/SurfaceDeformation.h"
0019 #include "Geometry/CommonTopologies/interface/TwoBowedSurfacesDeformation.h"
0020
0021 // This class's header
0022 #include "Alignment/CommonAlignmentParametrization/interface/TwoBowedSurfacesAlignmentParameters.h"
0023
0024 #include <cmath>
0025 #include <iostream>
0026 //_________________________________________________________________________________________________
0027 TwoBowedSurfacesAlignmentParameters::TwoBowedSurfacesAlignmentParameters(Alignable *ali)
0028     : AlignmentParameters(ali, AlgebraicVector(N_PARAM), AlgebraicSymMatrix(N_PARAM, 0)),
0029       ySplit_(this->ySplitFromAlignable(ali)) {}
0030
0031 //_________________________________________________________________________________________________
0032 TwoBowedSurfacesAlignmentParameters ::TwoBowedSurfacesAlignmentParameters(Alignable *alignable,
0033                                                                           const AlgebraicVector &parameters,
0034                                                                           const AlgebraicSymMatrix &covMatrix)
0035     : AlignmentParameters(alignable, parameters, covMatrix), ySplit_(this->ySplitFromAlignable(alignable)) {
0036   if (parameters.num_row() != N_PARAM) {
0037     throw cms::Exception("BadParameters") << "in TwoBowedSurfacesAlignmentParameters(): " << parameters.num_row()
0038                                           << " instead of " << N_PARAM << " parameters.";
0039   }
0040 }
0041
0042 //_________________________________________________________________________________________________
0043 TwoBowedSurfacesAlignmentParameters ::TwoBowedSurfacesAlignmentParameters(Alignable *alignable,
0044                                                                           const AlgebraicVector &parameters,
0045                                                                           const AlgebraicSymMatrix &covMatrix,
0046                                                                           const std::vector<bool> &selection)
0047     : AlignmentParameters(alignable, parameters, covMatrix, selection), ySplit_(this->ySplitFromAlignable(alignable)) {
0048   if (parameters.num_row() != N_PARAM) {
0049     throw cms::Exception("BadParameters") << "in TwoBowedSurfacesAlignmentParameters(): " << parameters.num_row()
0050                                           << " instead of " << N_PARAM << " parameters.";
0051   }
0052 }
0053
0054 //_________________________________________________________________________________________________
0055 TwoBowedSurfacesAlignmentParameters *TwoBowedSurfacesAlignmentParameters::clone(
0056     const AlgebraicVector &parameters, const AlgebraicSymMatrix &covMatrix) const {
0057   TwoBowedSurfacesAlignmentParameters *rbap =
0058       new TwoBowedSurfacesAlignmentParameters(this->alignable(), parameters, covMatrix, selector());
0059
0060   if (this->userVariables())
0061     rbap->setUserVariables(this->userVariables()->clone());
0062   rbap->setValid(this->isValid());
0063
0064   return rbap;
0065 }
0066
0067 //_________________________________________________________________________________________________
0068 TwoBowedSurfacesAlignmentParameters *TwoBowedSurfacesAlignmentParameters::cloneFromSelected(
0069     const AlgebraicVector &parameters, const AlgebraicSymMatrix &covMatrix) const {
0070   return this->clone(this->expandVector(parameters, this->selector()),
0071                      this->expandSymMatrix(covMatrix, this->selector()));
0072 }
0073
0074 //_________________________________________________________________________________________________
0075 AlgebraicMatrix TwoBowedSurfacesAlignmentParameters::derivatives(const TrajectoryStateOnSurface &tsos,
0076                                                                  const AlignableDetOrUnitPtr &alidet) const {
0077   const Alignable *ali = this->alignable();  // Alignable of these parameters
0078   AlgebraicMatrix result(N_PARAM, 2);        // initialised with zeros
0079
0080   if (ali == alidet) {
0081     const AlignableSurface &surf = ali->surface();
0082
0083     // matrix of dimension BowedDerivs::N_PARAM x 2
0084     const AlgebraicMatrix derivs(BowedDerivs()(tsos, surf.width(), surf.length(), true, ySplit_));  // split at ySplit_!
0085
0086     // Parameters belong to surface part with y < ySplit_ or y >= ySplit_?
0087     const double localY = tsos.localParameters().mixedFormatVector()[4];
0088     const unsigned int indexOffset = (localY < ySplit_ ? 0 : dx2);
0089     // Copy derivatives to relevant part of result
0090     for (unsigned int i = BowedDerivs::dx; i < BowedDerivs::N_PARAM; ++i) {
0091       result[indexOffset + i][0] = derivs[i][0];
0092       result[indexOffset + i][1] = derivs[i][1];
0093     }
0094   } else {
0095     // The following is even more difficult for
0096     // TwoBowedSurfacesAlignmentParameters than for
0097     // BowedSurfaceAlignmentParameters where this text comes from:
0098     //
0099     // We could give this a meaning by applying frame-to-frame derivatives
0100     // to the rigid body part of the parameters (be careful that alpha ~=
0101     // dslopeY and beta ~= -dslopeX, but with changed scale!) and keep the
0102     // surface structure parameters untouched in local meaning. In this way we
0103     // could do higher level alignment and determine 'average' surface
0104     // structures for the components.
0105     throw cms::Exception("MisMatch") << "TwoBowedSurfacesAlignmentParameters::derivatives: The hit "
0106                                         "alignable must match the "
0107                                      << "aligned one (i.e. bowed surface parameters cannot be used for "
0108                                         "composed alignables)\n";
0109   }
0110
0111   return result;
0112 }
0113
0114 //_________________________________________________________________________________________________
0115 void TwoBowedSurfacesAlignmentParameters::apply() {
0116   Alignable *alignable = this->alignable();
0117   if (!alignable) {
0118     throw cms::Exception("BadParameters") << "TwoBowedSurfacesAlignmentParameters::apply: parameters without "
0119                                              "alignable";
0120   }
0121
0122   // Some repeatedly needed variables
0123   const AlignableSurface &surface = alignable->surface();
0124   const double halfLength = surface.length() * 0.5;         // full module
0125   const double halfLength1 = (halfLength + ySplit_) * 0.5;  // low-y surface
0126   const double halfLength2 = (halfLength - ySplit_) * 0.5;  // high-y surface
0127
0128   // first copy the parameters into separate parts for the two surfaces
0129   const AlgebraicVector &params = theData->parameters();
0130   std::vector<double> rigidBowPar1(BowedDerivs::N_PARAM);  // 1st surface (y <  ySplit_)
0131   std::vector<double> rigidBowPar2(BowedDerivs::N_PARAM);  // 2nd surface (y >= ySplit_)
0132   for (unsigned int i = 0; i < BowedDerivs::N_PARAM; ++i) {
0133     rigidBowPar1[i] = params[i];
0134     rigidBowPar2[i] = params[i + BowedDerivs::N_PARAM];
0135   }
0136   // Now adjust slopes to angles, note that dslopeX <-> -beta & dslopeY <->
0137   // alpha, see BowedSurfaceAlignmentParameters::rotation(): FIXME: use atan?
0138   rigidBowPar1[3] = params[dslopeY1] / halfLength1;               // alpha1
0139   rigidBowPar2[3] = params[dslopeY2] / halfLength2;               // alpha2
0140   rigidBowPar1[4] = -params[dslopeX1] / (surface.width() * 0.5);  // beta1
0141   rigidBowPar2[4] = -params[dslopeX2] / (surface.width() * 0.5);  // beta2
0142   // gamma is simply scaled
0143   const double gammaScale1 = BowedDerivs::gammaScale(surface.width(), 2.0 * halfLength1);
0144   rigidBowPar1[5] = params[drotZ1] / gammaScale1;
0145   //   const double gammaScale2 = std::sqrt(halfLength2 * halfLength2
0146   //                       + surface.width() * surface.width()/4.);
0147   const double gammaScale2 = BowedDerivs::gammaScale(surface.width(), 2.0 * halfLength2);
0148   rigidBowPar2[5] = params[drotZ2] / gammaScale2;
0149
0150   // Get rigid body rotations of full module as mean of the two surfaces:
0151   align::EulerAngles angles(3);  // to become 'common' rotation in local frame
0152   for (unsigned int i = 0; i < 3; ++i) {
0153     angles[i] = (rigidBowPar1[i + 3] + rigidBowPar2[i + 3]) * 0.5;
0154   }
0155   // Module rotations are around other axes than the one we determined,
0156   // so we have to correct that the surfaces are shifted by the rotation around
0157   // the module axis - in linear approximation just an additional shift:
0158   const double yMean1 = -halfLength + halfLength1;  // y of alpha1 rotation axis in module frame
0159   const double yMean2 = halfLength - halfLength2;   // y of alpha2 rotation axis in module frame
0160   rigidBowPar1[dz1] -= angles[0] * yMean1;          // correct w1 for alpha
0161   rigidBowPar2[dz1] -= angles[0] * yMean2;          // correct w2 for alpha
0162   // Nothing for beta1/2 since anyway both around the y-axis of the module.
0163   rigidBowPar1[dx1] += angles[2] * yMean1;  // correct x1 for gamma
0164   rigidBowPar2[dx1] += angles[2] * yMean2;  // correct x1 for gamma
0165
0166   // Get rigid body shifts of full module as mean of the two surfaces:
0167   const align::LocalVector shift((rigidBowPar1[dx1] + rigidBowPar2[dx1]) * 0.5,   // dx1!
0168                                  (rigidBowPar1[dy1] + rigidBowPar2[dy1]) * 0.5,   // dy1!
0169                                  (rigidBowPar1[dz1] + rigidBowPar2[dz1]) * 0.5);  // dz1!
0170   // Apply module shift and rotation:
0171   alignable->move(surface.toGlobal(shift));
0172   // original code:
0173   //  alignable->rotateInLocalFrame( align::toMatrix(angles) );
0174   // correct for rounding errors:
0175   align::RotationType rot(surface.toGlobal(align::toMatrix(angles)));
0176   align::rectify(rot);
0177   alignable->rotateInGlobalFrame(rot);
0178
0179   // only update the surface deformations if they were selected for alignment
0180   if (selector()[dsagittaX1] || selector()[dsagittaXY1] || selector()[dsagittaY1] || selector()[dsagittaX2] ||
0181       selector()[dsagittaXY2] || selector()[dsagittaY2]) {
0182     // Fill surface structures with mean bows and half differences for all
0183     // parameters:
0184     std::vector<align::Scalar> deformations;
0185     deformations.reserve(13);
0186     // first part: average bows
0187     deformations.push_back((params[dsagittaX1] + params[dsagittaX2]) * 0.5);
0188     deformations.push_back((params[dsagittaXY1] + params[dsagittaXY2]) * 0.5);
0189     deformations.push_back((params[dsagittaY1] + params[dsagittaY2]) * 0.5);
0190     // second part: half difference of all corrections
0191     for (unsigned int i = 0; i < BowedDerivs::N_PARAM; ++i) {
0192       // sign means that we have to apply e.g.
0193       // - sagittaX for sensor 1: deformations[0] + deformations[9]
0194       // - sagittaX for sensor 2: deformations[0] - deformations[9]
0195       // - additional dx for sensor 1:  deformations[3]
0196       // - additional dx for sensor 2: -deformations[3]
0197       deformations.push_back((rigidBowPar1[i] - rigidBowPar2[i]) * 0.5);
0198     }
0199     // finally: keep track of where we have split the module
0200     deformations.push_back(ySplit_);  // index is 12
0201
0202     const TwoBowedSurfacesDeformation deform{deformations};
0203
0204     // FIXME: true to propagate down?
0205     //        Needed for hierarchy with common deformation parameter,
0206     //        but that is not possible now anyway.
0207     alignable->addSurfaceDeformation(&deform, false);
0208   }
0209 }
0210
0211 //_________________________________________________________________________________________________
0212 int TwoBowedSurfacesAlignmentParameters::type() const { return AlignmentParametersFactory::kTwoBowedSurfaces; }
0213
0214 //_________________________________________________________________________________________________
0215 void TwoBowedSurfacesAlignmentParameters::print() const {
0216   std::cout << "Contents of TwoBowedSurfacesAlignmentParameters:"
0217             << "\nParameters: " << theData->parameters() << "\nCovariance: " << theData->covariance() << std::endl;
0218 }
0219
0220 //_________________________________________________________________________________________________
0221 double TwoBowedSurfacesAlignmentParameters::ySplitFromAlignable(const Alignable *ali) const {
0222   if (!ali)
0223     return 0.;
0224
0225   const align::PositionType pos(ali->globalPosition());
0226   const double r = pos.perp();
0227
0228   // The returned numbers for TEC are calculated as stated below from
0229   // what is found in CMS-NOTE 2003/20.
0230   // Note that at that time it was planned to use ST sensors for the outer TEC
0231   // while in the end there are only a few of them in the tracker - the others
0232   // are HPK. No idea whether there are subtle changes in geometry. The numbers
0233   // have been cross checked with y-residuals in data, see
0234   // https://hypernews.cern.ch/HyperNews/CMS/get/recoTracking/1018/1/1/2/1/1/1/2/1.html.
0235
0236   if (r < 58.) {  // Pixel, TIB, TID or TEC ring 1-4
0237     edm::LogError("Alignment") << "@SUB=TwoBowedSurfacesAlignmentParameters::ySplitFromAlignable"
0238                                << "There are no split modules for radii < 58, but r = " << r;
0239     return 0.;
0240   } else if (fabs(pos.z()) < 118.) {  // TOB
0241     return 0.;
0242   } else if (r > 90.) {  // TEC ring 7
0243     // W7a Height active= 106.900mm (Tab. 2) (but 106.926 mm p. 40!?)
0244     // W7a R min active = 888.400mm (Paragraph 5.5.7)
0245     // W7a R max active = W7a R min active + W7a Height active =
0246     //                  = 888.400mm + 106.900mm = 995.300mm
0247     // W7b Height active=  94.900mm (Tab. 2)  (but 94.876 mm p. 41!?)
0248     // W7b R min active = 998.252mm (Paragraph 5.5.8)
0249     // W7b R max active = 998.252mm + 94.900mm = 1093.152mm
0250     // mean radius module = 0.5*(1093.152mm + 888.400mm) = 990.776mm
0251     // mean radius gap = 0.5*(998.252mm + 995.300mm) = 996.776mm
0252     // ySplit = (mean radius gap - mean radius module) // local y and r have
0253     // same directions!
0254     //        = 996.776mm - 990.776mm = 6mm
0255     return 0.6;
0256   } else if (r > 75.) {  // TEC ring 6
0257     // W6a Height active= 96.100mm (Tab. 2) (but 96.136 mm p. 38!?)
0258     // W6a R min active = 727.000mm (Paragraph 5.5.5)
0259     // W6a R max active = W6a R min active + W6a Height active =
0260     //                  = 727.000mm + 96.100mm = 823.100mm
0261     // W6b Height active= 84.900mm (Tab. 2) (but 84.936 mm p. 39!?)
0262     // W6b R min active = 826.060mm (Paragraph 5.5.6)
0263     // W6b R max active = 826.060mm + 84.900mm = 910.960mm
0264     // mean radius module = 0.5*(910.960mm + 727.000mm) = 818.980mm
0265     // mean radius gap = 0.5*(826.060mm + 823.100mm) = 824.580mm
0266     //          -1: local y and r have opposite directions!
0267     // ySplit = -1*(mean radius gap - mean radius module)
0268     //        = -1*(824.580mm - 818.980mm) = -5.6mm
0269     return -0.56;
0270   } else {  // TEC ring 5 - smaller radii alreay excluded before
0271     // W5a Height active= 81.200mm (Tab. 2) (but 81.169 mm p. 36!?)
0272     // W5a R min active = 603.200mm (Paragraph 5.5.3)
0273     // W5a R max active = W5a R min active + W5a Height active =
0274     //                  = 603.200mm + 81.200mm = 684.400mm
0275     // W5b Height active= 63.200mm (Tab. 2) (63.198 mm on p. 37)
0276     // W5b R min active = 687.293mm (Abschnitt 5.5.4 der note)
0277     // W5b R max active = 687.293mm + 63.200mm = 750.493mm
0278     // mean radius module = 0.5*(750.493mm + 603.200mm) = 676.8465mm
0279     // mean radius gap = 0.5*(687.293mm + 684.400mm) = 685.8465mm
0280     //          -1: local y and r have opposite directions!
0281     // ySplit = -1*(mean radius gap - mean radius module)
0282     //        = -1*(685.8465mm - 676.8465mm) = -9mm
0283     return -0.9;
0284   }
0285   //  return 0.;
0286 }