Merge pull request #2592 from frmeier/patch_hierachyconstraints62XSLHC

Patch hierarchyconstraints to CMSSW_6_2_X_SHLC

Alignment/CommonAlignmentParametrization/interface/FrameToFrameDerivative.h

@@ -6,7 +6,10 @@
 
 /// \class FrameToFrameDerivative
 ///
-/// class for calculating derivatives between different frames
+/// Class for calculating the jacobian d_object/d_composedObject
+/// for the rigid body parametrisation of both, i.e. the derivatives
+/// expressing the influence of u, v, w, alpha, beta, gamma of the
+/// composedObject on u, v, w, alpha, beta, gamma of its component 'object'.
 ///
 ///  $Date: 2010/12/14 01:02:34 $
 ///  $Revision: 1.7 $
@@ -18,12 +21,25 @@ class FrameToFrameDerivative
 {
   public:
 
-  /// Return the derivative DeltaFrame(object)/DeltaFrame(composedobject) 
+  /// Return the derivative DeltaFrame(object)/DeltaFrame(composedObject),
+  /// i.e. a 6x6 matrix:
+  ///
+  /// / du/du_c du/dv_c du/dw_c du/da_c du/db_c du/dg_c |
+  /// | dv/du_c dv/dv_c dv/dw_c dv/da_c dv/db_c dv/dg_c |
+  /// | dw/du_c dw/dv_c dw/dw_c dw/da_c dw/db_c dw/dg_c |
+  /// | da/du_c da/dv_c da/dw_c da/da_c da/db_c da/dg_c |
+  /// | db/du_c db/dv_c db/dw_c db/da_c db/db_c db/dg_c |
+  /// \ dg/du_c dg/dv_c dg/dw_c dg/da_c dg/db_c dg/dg_c /
+  ///
+  /// where u, v, w, a, b, g are shifts and rotations of the object
+  /// and u_c, v_c, w_c, a_c, b_c, g_c those of the composed object.
+
   AlgebraicMatrix frameToFrameDerivative(const Alignable* object,
 					 const Alignable* composedObject) const;
 
   /// Calculates derivatives DeltaFrame(object)/DeltaFrame(composedobject) 
-  /// using their positions and orientations.
+  /// using their positions and orientations, see method frameToFrameDerivative(..)
+  /// for definition.
   /// As a new method it gets a new interface avoiding CLHEP that should anyway be
   /// replaced by SMatrix at some point...
   AlgebraicMatrix66 getDerivative(const align::RotationType &objectRot,

Alignment/CommonAlignmentParametrization/interface/KarimakiAlignmentDerivatives.h

@@ -17,7 +17,15 @@ class KarimakiAlignmentDerivatives
 {
 public:
 
-  /// Returns 6x2 jacobian matrix
+  /// Returns 6x2 jacobian matrix of derivatives of residuals in x and y
+  /// with respect to rigid body aligment parameters:
+  ///
+  /// / dr_x/du  dr_y/du |
+  /// | dr_x/dv  dr_y/dv |
+  /// | dr_x/dw  dr_y/dw |
+  /// | dr_x/da  dr_y/da |
+  /// | dr_x/db  dr_y/db |
+  /// \ dr_x/dg  dr_y/dg /
   AlgebraicMatrix operator()(const TrajectoryStateOnSurface &tsos) const;
 
 };

Alignment/CommonAlignmentParametrization/interface/ParametersToParametersDerivatives.h

@@ -3,11 +3,33 @@
 
 /// \class ParametersToParametersDerivatives
 ///
-/// Class for calculating derivatives for hierarchies between different kind
-/// of alignment parameters (note that not all combinations might be supported!),
-/// needed e.g. to formulate constraints to remove the additional degrees of
-/// freedom introduced if larger structure and their components are aligned
-/// simultaneously.
+/// Class for getting the jacobian d_mother/d_component for various kinds
+/// of alignment parametrisations, i.e. the derivatives expressing the influence
+/// of the parameters of the 'component' on the parameters of its 'mother'.
+/// This is needed e.g. to formulate constraints to remove the additional
+/// degrees of freedom introduced if larger structures and their components
+/// are aligned simultaneously.
+/// The jacobian matrix is
+///
+/// / dp1_l/dp1_i dp1_l/dp2_i  ...  dp1_l/dpn_i |
+/// | dp2_l/dp1_i dp2_l/dp2_i  ...  dp2_l/dpn_i |
+/// |      .           .                 .      |
+/// |      .           .                 .      |
+/// |      .           .                 .      |
+/// \ dpm_l/dpm_i dpm_l/dpm_i  ...  dpm_l/dpn_i /
+///
+/// where 
+/// p1_l, p2_l, ..., pn_l are the n parameters of the composite 'mother' object
+/// and
+/// p1_i, p2_i, ..., pm_i are the m parameters of its component.
+///
+/// Note that not all combinations of parameters are supported:
+/// Please check method isOK() before accessing the derivatives via 
+/// operator(unsigned int indParMother, unsigned int indParComp).
+///
+/// Currently these parameters are supported:
+/// - mother: rigid body parameters,
+/// - component: rigid body, bowed surface or two bowed surfaces parameters.
 ///
 ///  $Date: 2010/12/14 01:08:25 $
 ///  $Revision: 1.2 $
@@ -26,8 +48,9 @@ class ParametersToParametersDerivatives
   /// Indicate whether able to provide the derivatives.
   bool isOK() const { return isOK_;}
 
-  /// Return the derivative DeltaParam(object)/DeltaParam(composedobject), indices start with 0.
-  /// But check isOK() first!
+  /// Return the derivative DeltaParam(mother)/DeltaParam(component).
+  /// Indices start with 0 - but check isOK() first!
+  /// See class description about matrix.
   double operator() (unsigned int indParMother, unsigned int indParComp) const; 
 
   // Not this - would make the internals public:
@@ -45,8 +68,10 @@ class ParametersToParametersDerivatives
   bool init2BowedRigid(const Alignable &component, const Alignable &mother);
 
   typedef ROOT::Math::SMatrix<double,6,9,ROOT::Math::MatRepStd<double,6,9> > AlgebraicMatrix69;
-  AlgebraicMatrix69 dBowed_dRigid(const AlgebraicMatrix66 &f2f,
-				  double halfWidth, double halfLength) const;
+  /// from d(rigid_mother)/d(rigid_component) to d(rigid_mother)/d(bowed_component)
+  /// for bad input (length or width zero), set object to invalid: isOK_ = false
+  AlgebraicMatrix69 dRigid_dBowed(const AlgebraicMatrix66 &dRigidM2dRigidC,
+				  double halfWidth, double halfLength);
 
   /// data members
   bool            isOK_; /// can we provide the desired?

Alignment/CommonAlignmentParametrization/src/BeamSpotAlignmentParameters.cc

@@ -10,7 +10,6 @@
 #include "Alignment/CommonAlignment/interface/Utilities.h"
 #include "Alignment/CommonAlignment/interface/Alignable.h"
 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
-#include "Alignment/CommonAlignmentParametrization/interface/FrameToFrameDerivative.h"
 #include "Alignment/CommonAlignmentParametrization/interface/BeamSpotAlignmentDerivatives.h"
 #include "Alignment/CommonAlignmentParametrization/interface/AlignmentParametersFactory.h"
 #include "CondFormats/Alignment/interface/Definitions.h"
@@ -97,10 +96,11 @@ BeamSpotAlignmentParameters::derivatives( const TrajectoryStateOnSurface &tsos,
 
   if (ali == alidet) { // same alignable => same frame
     return BeamSpotAlignmentDerivatives()(tsos);
-  } else { // different alignable => transform into correct frame
-    const AlgebraicMatrix deriv = BeamSpotAlignmentDerivatives()(tsos);
-    FrameToFrameDerivative ftfd;
-    return ftfd.frameToFrameDerivative(alidet, ali) * deriv;
+  } else {
+    throw cms::Exception("MisMatch")
+      << "BeamSpotAlignmentParameters::derivatives: The hit alignable must match the "
+      << "aligned one, i.e. these parameters make only sense for AlignableBeamSpot.\n";
+    return AlgebraicMatrix(N_PARAM, 2); // please compiler
   }
 }
 

Alignment/CommonAlignmentParametrization/src/FrameToFrameDerivative.cc

@@ -96,8 +96,7 @@ FrameToFrameDerivative::getDerivative(const align::RotationType &objectRot,
   derivative[5][4] = derivBB[2][1];
   derivative[5][5] = derivBB[2][2];
 
-  return(derivative.T());
-
+  return derivative;
 }
 
 

Alignment/CommonAlignmentParametrization/src/ParametersToParametersDerivatives.cc

@@ -59,33 +59,46 @@ bool ParametersToParametersDerivatives::init(const Alignable &component, int typ
 bool ParametersToParametersDerivatives::initRigidRigid(const Alignable &component,
 						       const Alignable &mother)
 {
-  // simply frame to frame!
-  FrameToFrameDerivative f2fDerivMaker;
-  AlgebraicMatrix66 m(asSMatrix<6,6>(f2fDerivMaker.frameToFrameDerivative(&component, &mother)));
-
-  // copy to TMatrix
-  derivatives_.ResizeTo(6,6);
-  derivatives_.SetMatrixArray(m.begin());
-
-  return true;
+  // See G. Flucke's presentation from  20 Feb 2007
+  // https://indico.cern.ch/contributionDisplay.py?contribId=15&sessionId=1&confId=10930
+  // and C. Kleinwort's one from 14 Feb 2013
+  // https://indico.cern.ch/contributionDisplay.py?contribId=14&sessionId=1&confId=224472
+
+  FrameToFrameDerivative f2f;
+  // frame2frame returns dcomponent/dmother for both being rigid body, so we have to invert
+  AlgebraicMatrix66 m(asSMatrix<6,6>(f2f.frameToFrameDerivative(&component, &mother)));
+
+  if (m.Invert()) { // now matrix is d(rigid_mother)/d(rigid_component)
+    // copy to TMatrix
+    derivatives_.ResizeTo(6,6);
+    derivatives_.SetMatrixArray(m.begin());
+    return true;
+  } else {
+    return false;
+  }
 }
 
 //_________________________________________________________________________________________________
 bool ParametersToParametersDerivatives::initBowedRigid(const Alignable &component,
 						       const Alignable &mother)
 {
   // component is bowed surface, mother rigid body
-  FrameToFrameDerivative f2fMaker;
-  const AlgebraicMatrix66 f2f(asSMatrix<6,6>(f2fMaker.frameToFrameDerivative(&component,&mother)));
-  const double halfWidth  = 0.5 * component.surface().width();
-  const double halfLength = 0.5 * component.surface().length();
-  const AlgebraicMatrix69 m(this->dBowed_dRigid(f2f, halfWidth, halfLength));
-
-  // copy to TMatrix
-  derivatives_.ResizeTo(6,9);
-  derivatives_.SetMatrixArray(m.begin());
-
-  return true;
+  FrameToFrameDerivative f2f;
+  // frame2frame returns dcomponent/dmother for both being rigid body, so we have to invert
+  AlgebraicMatrix66 rigM2rigC(asSMatrix<6,6>(f2f.frameToFrameDerivative(&component,&mother)));
+  if (rigM2rigC.Invert()) { // now matrix is d(rigid_mother)/d(rigid_component)
+    const double halfWidth  = 0.5 * component.surface().width();
+    const double halfLength = 0.5 * component.surface().length();
+    const AlgebraicMatrix69 m(this->dRigid_dBowed(rigM2rigC, halfWidth, halfLength));
+
+    // copy to TMatrix
+    derivatives_.ResizeTo(6,9);
+    derivatives_.SetMatrixArray(m.begin());
+
+    return true;
+  } else {
+    return false;
+  }
 }
 
 //_________________________________________________________________________________________________
@@ -118,15 +131,19 @@ bool ParametersToParametersDerivatives::init2BowedRigid(const Alignable &compone
   const align::GlobalPoint posSurf2(component.surface().toGlobal(align::LocalPoint(0.,yM2,0.)));
 
   // 2) get derivatives for both,
+  // getDerivative(..) returns dcomponent/dmother for both being rigid body
   FrameToFrameDerivative f2fMaker;
-  const AlgebraicMatrix66 f2fSurf1(f2fMaker.getDerivative(component.globalRotation(),
-                                                          mother.globalRotation(),
-                                                          posSurf1, mother.globalPosition()));
-  const AlgebraicMatrix66 f2fSurf2(f2fMaker.getDerivative(component.globalRotation(),
-                                                          mother.globalRotation(),
-                                                          posSurf2, mother.globalPosition()));
-  const AlgebraicMatrix69 derivs1(this->dBowed_dRigid(f2fSurf1, halfWidth, halfLength1));
-  const AlgebraicMatrix69 derivs2(this->dBowed_dRigid(f2fSurf2, halfWidth, halfLength2));
+  AlgebraicMatrix66 f2fSurf1(f2fMaker.getDerivative(component.globalRotation(),
+                                                    mother.globalRotation(),
+                                                    posSurf1, mother.globalPosition()));
+  AlgebraicMatrix66 f2fSurf2(f2fMaker.getDerivative(component.globalRotation(),
+                                                    mother.globalRotation(),
+                                                    posSurf2, mother.globalPosition()));
+  // We have to invert matrices to get d(rigid_mother)/d(rigid_component):
+  if (!f2fSurf1.Invert() || !f2fSurf2.Invert()) return false; // bail out if bad inversion
+  // Now get d(rigid_mother)/d(bowed_component):
+  const AlgebraicMatrix69 derivs1(this->dRigid_dBowed(f2fSurf1, halfWidth, halfLength1));
+  const AlgebraicMatrix69 derivs2(this->dRigid_dBowed(f2fSurf2, halfWidth, halfLength2));
 
   // 3) fill the common matrix by merging the two.
   typedef ROOT::Math::SMatrix<double,6,18,ROOT::Math::MatRepStd<double,6,18> > AlgebraicMatrix6_18;
@@ -143,36 +160,47 @@ bool ParametersToParametersDerivatives::init2BowedRigid(const Alignable &compone
 
 //_________________________________________________________________________________________________
 ParametersToParametersDerivatives::AlgebraicMatrix69
-ParametersToParametersDerivatives::dBowed_dRigid(const AlgebraicMatrix66 &f2f,
-						 double halfWidth, double halfLength) const
+ParametersToParametersDerivatives::dRigid_dBowed(const AlgebraicMatrix66 &dRigidM2dRigidC,
+						 double halfWidth, double halfLength)
 {
   typedef BowedSurfaceAlignmentDerivatives BowedDerivs;
   const double gammaScale = BowedDerivs::gammaScale(2.*halfWidth, 2.*halfLength);
 
-  // 1st index (column) is parameter of the mother (<6),
-  // 2nd index (row) that of component (<9):
+  // 'dRigidM2dRigidC' is dmother/dcomponent for both being rigid body
+  // final matrix will be dmother/dcomponent for mother as rigid body, component with bows
+  // 1st index (row) is parameter of the mother (0..5),
+  // 2nd index (column) that of component (0..8):
   AlgebraicMatrix69 derivs;
+  if (0. == gammaScale || 0. == halfWidth || 0. == halfLength) {
+    isOK_ = false; // bad input - we would have to devide by that in the following!
+    edm::LogError("Alignment") << "@SUB=ParametersToParametersDerivatives::dRigid_dBowed"
+			       << "Some zero length as input.";
+    return derivs;
+  }
 
-  for (unsigned int iRow = 0; iRow < 6; ++iRow) { // 6 rigid body parameters of mother
+  for (unsigned int iRow = 0; iRow < AlgebraicMatrix69::kRows; ++iRow) {
+    // loop on 6 rigid body parameters of mother
     // First copy the common rigid body part, e.g.:
-    // - (0,0): du_comp/du_moth
-    // - (0,1): dv_comp/du_moth
-    // - (1,2): dw_comp/dv_moth
+    // (0,0): du_moth/du_comp, (0,1): dv_moth/du_comp, (1,2): dw_moth/dv_comp
     for (unsigned int iCol = 0; iCol < 3; ++iCol) { // 3 movements of component
-      derivs(iRow, iCol) = f2f(iRow, iCol);  
+      derivs(iRow, iCol) = dRigidM2dRigidC(iRow, iCol);  
     }
 
-    // Now we have to take care of order and scales for rotation-like parameters:
-    // slopeX -> halfWidth * beta
-    derivs(iRow, 3) = halfWidth  * f2f(iRow, 4); // = dslopeX_c/dpar_m = hw * db_c/dpar_m
-    // slopeY -> halfLength * alpha
-    derivs(iRow, 4) = halfLength * f2f(iRow, 3); // = dslopeY_c/dpar_m = hl * da_c/dpar_m
-    // rotZ -> gammaScale * gamma
-    derivs(iRow, 5) = gammaScale * f2f(iRow, 5); // = drotZ_c/dpar_m = gscale * dg_c/dpar_m
-
-    // Finally, movements and rotations have no influence on surface internals:
-    for (unsigned int iCol = 6; iCol < 9; ++iCol) { // 3 sagittae of component
-      derivs(iRow, iCol) = 0.;  
+    // Now we have to take care of order, signs and scales for rotation-like parameters,
+    // see CMS AN-2011/531:
+    // slopeX = w10 = -halfWidth * beta
+    // => dpar_m/dslopeX_comp = dpar_m/d(-hw * beta_comp) = -(dpar_m/dbeta_comp)/hw
+    derivs(iRow, 3) = -dRigidM2dRigidC(iRow, 4)/halfWidth; 
+    // slopeY = w10 = +halfLength * alpha
+    // => dpar_m/dslopeY_comp = dpar_m/d(+hl * alpha_comp) = (dpar_m/dalpha_comp)/hl
+    derivs(iRow, 4) = dRigidM2dRigidC(iRow, 3)/halfLength;
+    // rotZ = gammaScale * gamma
+    // => dpar_m/drotZ_comp = dpar_m/d(gamma_comp * gscale) = (dpar_m/dgamma)/gscale
+    derivs(iRow, 5) = dRigidM2dRigidC(iRow, 5)/gammaScale;
+
+    // Finally, sensor internals like their curvatures have no influence on mother:
+    for (unsigned int iCol = 6; iCol < AlgebraicMatrix69::kCols; ++iCol) {
+      derivs(iRow, iCol) = 0.; // 3 sagittae of component
     }
   }
 

Alignment/CommonAlignmentParametrization/src/RigidBodyAlignmentParameters.cc

@@ -93,7 +93,7 @@ RigidBodyAlignmentParameters::derivatives( const TrajectoryStateOnSurface &tsos,
   } else { // different alignable => transform into correct frame
     const AlgebraicMatrix deriv = KarimakiAlignmentDerivatives()(tsos);
     FrameToFrameDerivative ftfd;
-    return ftfd.frameToFrameDerivative(alidet, ali) * deriv;
+    return ftfd.frameToFrameDerivative(alidet, ali).T() * deriv;
   }
 }
 

Alignment/CommonAlignmentParametrization/src/RigidBodyAlignmentParameters4D.cc

@@ -29,7 +29,7 @@ RigidBodyAlignmentParameters4D::derivatives( const TrajectoryStateOnSurface &tso
   } else { // different alignable => transform into correct frame
     const AlgebraicMatrix deriv = SegmentAlignmentDerivatives4D()(tsos);
     FrameToFrameDerivative ftfd;
-    return ftfd.frameToFrameDerivative(alidet, ali) * deriv;
+    return ftfd.frameToFrameDerivative(alidet, ali).T() * deriv;
   }
 }