CommonAlignmentParametrization/src/ParametersToParametersDerivatives.cc

0001 /** \file ParametersToParametersDerivatives.cc
0002  *
0003  *  $Date: 2010/12/09 19:53:42 $
0004  *  $Revision: 1.1 $
0005  */
0006
0007 #include "Alignment/CommonAlignmentParametrization/interface/ParametersToParametersDerivatives.h"
0008
0009 #include "Alignment/CommonAlignment/interface/Alignable.h"
0010 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
0011 #include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
0012 #include "Alignment/CommonAlignmentParametrization/interface/BowedSurfaceAlignmentDerivatives.h"
0013 #include "Alignment/CommonAlignmentParametrization/interface/FrameToFrameDerivative.h"
0014 #include "Alignment/CommonAlignmentParametrization/interface/TwoBowedSurfacesAlignmentParameters.h"
0015 #include "CondFormats/Alignment/interface/Definitions.h"
0016
0017 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 // already in header:
0020 // #include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
0021
0022 //_________________________________________________________________________________________________
0023 ParametersToParametersDerivatives ::ParametersToParametersDerivatives(const Alignable &component,
0024                                                                       const Alignable &mother)
0025     : isOK_(component.alignmentParameters() && mother.alignmentParameters()) {
0026   if (isOK_) {
0027     isOK_ =
0028         this->init(component, component.alignmentParameters()->type(), mother, mother.alignmentParameters()->type());
0029   }
0030 }
0031
0032 //_________________________________________________________________________________________________
0033 bool ParametersToParametersDerivatives::init(const Alignable &component,
0034                                              int typeComponent,
0035                                              const Alignable &mother,
0036                                              int typeMother) {
0037   using namespace AlignmentParametersFactory;  // for kRigidBody etc.
0038   if ((typeMother == kRigidBody || typeMother == kRigidBody4D) &&
0039       (typeComponent == kRigidBody || typeComponent == kRigidBody4D)) {
0040     return this->initRigidRigid(component, mother);
0041   } else if ((typeMother == kRigidBody || typeMother == kRigidBody4D) && typeComponent == kBowedSurface) {
0042     return this->initBowedRigid(component, mother);
0043   } else if ((typeMother == kRigidBody || typeMother == kRigidBody4D) && typeComponent == kTwoBowedSurfaces) {
0044     return this->init2BowedRigid(component, mother);
0045   } else {
0046     // missing: mother with bows and component without, i.e. having 'common' bow
0047     // parameters
0048     edm::LogError("Alignment") << "@SUB=ParametersToParametersDerivatives::init"
0049                                << "Mother " << parametersTypeName(parametersType(typeMother)) << ", component "
0050                                << parametersTypeName(parametersType(typeComponent)) << ": not supported.";
0051     return false;
0052   }
0053 }
0054
0055 //_________________________________________________________________________________________________
0056 bool ParametersToParametersDerivatives::initRigidRigid(const Alignable &component, const Alignable &mother) {
0057   // See G. Flucke's presentation from  20 Feb 2007
0058   // https://indico.cern.ch/contributionDisplay.py?contribId=15&sessionId=1&confId=10930
0059   // and C. Kleinwort's one from 14 Feb 2013
0060   // https://indico.cern.ch/contributionDisplay.py?contribId=14&sessionId=1&confId=224472
0061
0062   FrameToFrameDerivative f2f;
0063   // frame2frame returns dcomponent/dmother for both being rigid body, so we
0064   // have to invert
0065   AlgebraicMatrix66 m(asSMatrix<6, 6>(f2f.frameToFrameDerivative(&component, &mother)));
0066
0067   if (m.Invert()) {  // now matrix is d(rigid_mother)/d(rigid_component)
0068     // copy to TMatrix
0069     derivatives_.ResizeTo(6, 6);
0070     derivatives_.SetMatrixArray(m.begin());
0071     return true;
0072   } else {
0073     return false;
0074   }
0075 }
0076
0077 //_________________________________________________________________________________________________
0078 bool ParametersToParametersDerivatives::initBowedRigid(const Alignable &component, const Alignable &mother) {
0079   // component is bowed surface, mother rigid body
0080   FrameToFrameDerivative f2f;
0081   // frame2frame returns dcomponent/dmother for both being rigid body, so we
0082   // have to invert
0083   AlgebraicMatrix66 rigM2rigC(asSMatrix<6, 6>(f2f.frameToFrameDerivative(&component, &mother)));
0084   if (rigM2rigC.Invert()) {  // now matrix is d(rigid_mother)/d(rigid_component)
0085     const double halfWidth = 0.5 * component.surface().width();
0086     const double halfLength = 0.5 * component.surface().length();
0087     const AlgebraicMatrix69 m(this->dRigid_dBowed(rigM2rigC, halfWidth, halfLength));
0088
0089     // copy to TMatrix
0090     derivatives_.ResizeTo(6, 9);
0091     derivatives_.SetMatrixArray(m.begin());
0092
0093     return true;
0094   } else {
0095     return false;
0096   }
0097 }
0098
0099 //_________________________________________________________________________________________________
0100 bool ParametersToParametersDerivatives::init2BowedRigid(const Alignable &component, const Alignable &mother) {
0101   // component is two bowed surfaces, mother rigid body
0102   const TwoBowedSurfacesAlignmentParameters *aliPar =
0103       dynamic_cast<TwoBowedSurfacesAlignmentParameters *>(component.alignmentParameters());
0104
0105   if (!aliPar) {
0106     edm::LogError("Alignment") << "@SUB=ParametersToParametersDerivatives::init2BowedRigid"
0107                                << "dynamic_cast to TwoBowedSurfacesAlignmentParameters failed.";
0108     return false;
0109   }
0110
0111   // We treat the two surfaces as independent objects, i.e.
0112   // 1) get the global position of each surface, depending on the ySplit value,
0113   const double ySplit = aliPar->ySplit();
0114   const double halfWidth = 0.5 * component.surface().width();
0115   const double halfLength = 0.5 * component.surface().length();
0116   const double halfLength1 = 0.5 * (halfLength + ySplit);
0117   const double halfLength2 = 0.5 * (halfLength - ySplit);
0118   const double yM1 = 0.5 * (ySplit - halfLength);  // y_mean of surface 1
0119   const double yM2 = yM1 + halfLength;             // y_mean of surface 2
0120   // The sensor positions and orientations could be adjusted using
0121   // TwoBowedSurfacesDeformation attached to the component,
0122   // but that should be 2nd order effect.
0123   const align::GlobalPoint posSurf1(component.surface().toGlobal(align::LocalPoint(0., yM1, 0.)));
0124   const align::GlobalPoint posSurf2(component.surface().toGlobal(align::LocalPoint(0., yM2, 0.)));
0125
0126   // 2) get derivatives for both,
0127   // getDerivative(..) returns dcomponent/dmother for both being rigid body
0128   FrameToFrameDerivative f2fMaker;
0129   AlgebraicMatrix66 f2fSurf1(
0130       f2fMaker.getDerivative(component.globalRotation(), mother.globalRotation(), posSurf1, mother.globalPosition()));
0131   AlgebraicMatrix66 f2fSurf2(
0132       f2fMaker.getDerivative(component.globalRotation(), mother.globalRotation(), posSurf2, mother.globalPosition()));
0133   // We have to invert matrices to get d(rigid_mother)/d(rigid_component):
0134   if (!f2fSurf1.Invert() || !f2fSurf2.Invert())
0135     return false;  // bail out if bad inversion
0136   // Now get d(rigid_mother)/d(bowed_component):
0137   const AlgebraicMatrix69 derivs1(this->dRigid_dBowed(f2fSurf1, halfWidth, halfLength1));
0138   const AlgebraicMatrix69 derivs2(this->dRigid_dBowed(f2fSurf2, halfWidth, halfLength2));
0139
0140   // 3) fill the common matrix by merging the two.
0141   typedef ROOT::Math::SMatrix<double, 6, 18, ROOT::Math::MatRepStd<double, 6, 18>> AlgebraicMatrix6_18;
0142   AlgebraicMatrix6_18 derivs;
0143   derivs.Place_at(derivs1, 0, 0);  // left half
0144   derivs.Place_at(derivs2, 0, 9);  // right half
0145
0146   // copy to TMatrix
0147   derivatives_.ResizeTo(6, 18);
0148   derivatives_.SetMatrixArray(derivs.begin());
0149
0150   return true;
0151 }
0152
0153 //_________________________________________________________________________________________________
0154 ParametersToParametersDerivatives::AlgebraicMatrix69 ParametersToParametersDerivatives::dRigid_dBowed(
0155     const AlgebraicMatrix66 &dRigidM2dRigidC, double halfWidth, double halfLength) {
0156   typedef BowedSurfaceAlignmentDerivatives BowedDerivs;
0157   const double gammaScale = BowedDerivs::gammaScale(2. * halfWidth, 2. * halfLength);
0158
0159   // 'dRigidM2dRigidC' is dmother/dcomponent for both being rigid body
0160   // final matrix will be dmother/dcomponent for mother as rigid body, component
0161   // with bows 1st index (row) is parameter of the mother (0..5), 2nd index
0162   // (column) that of component (0..8):
0163   AlgebraicMatrix69 derivs;
0164   if (0. == gammaScale || 0. == halfWidth || 0. == halfLength) {
0165     isOK_ = false;  // bad input - we would have to devide by that in the following!
0166     edm::LogError("Alignment") << "@SUB=ParametersToParametersDerivatives::dRigid_dBowed"
0167                                << "Some zero length as input.";
0168     return derivs;
0169   }
0170
0171   for (unsigned int iRow = 0; iRow < AlgebraicMatrix69::kRows; ++iRow) {
0172     // loop on 6 rigid body parameters of mother
0173     // First copy the common rigid body part, e.g.:
0174     // (0,0): du_moth/du_comp, (0,1): dv_moth/du_comp, (1,2): dw_moth/dv_comp
0175     for (unsigned int iCol = 0; iCol < 3; ++iCol) {  // 3 movements of component
0176       derivs(iRow, iCol) = dRigidM2dRigidC(iRow, iCol);
0177     }
0178
0179     // Now we have to take care of order, signs and scales for rotation-like
0180     // parameters, see CMS AN-2011/531: slopeX = w10 = -halfWidth * beta
0181     // => dpar_m/dslopeX_comp = dpar_m/d(-hw * beta_comp) =
0182     // -(dpar_m/dbeta_comp)/hw
0183     derivs(iRow, 3) = -dRigidM2dRigidC(iRow, 4) / halfWidth;
0184     // slopeY = w10 = +halfLength * alpha
0185     // => dpar_m/dslopeY_comp = dpar_m/d(+hl * alpha_comp) =
0186     // (dpar_m/dalpha_comp)/hl
0187     derivs(iRow, 4) = dRigidM2dRigidC(iRow, 3) / halfLength;
0188     // rotZ = gammaScale * gamma
0189     // => dpar_m/drotZ_comp = dpar_m/d(gamma_comp * gscale) =
0190     // (dpar_m/dgamma)/gscale
0191     derivs(iRow, 5) = dRigidM2dRigidC(iRow, 5) / gammaScale;
0192
0193     // Finally, sensor internals like their curvatures have no influence on
0194     // mother:
0195     for (unsigned int iCol = 6; iCol < AlgebraicMatrix69::kCols; ++iCol) {
0196       derivs(iRow, iCol) = 0.;  // 3 sagittae of component
0197     }
0198   }
0199
0200   return derivs;
0201 }
0202
0203 //_________________________________________________________________________________________________
0204 double ParametersToParametersDerivatives::operator()(unsigned int indParMother, unsigned int indParComp) const {
0205   // Do range checks?
0206   return derivatives_(indParMother, indParComp);
0207 }