VertexTools/src/VertexDistance3D.cc

0001 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0002 #include <cfloat>
0003
0004 using namespace reco;
0005
0006 Measurement1D VertexDistance3D::signedDistance(const Vertex& vtx1,
0007                                                const Vertex& vtx2,
0008                                                const GlobalVector& momentum) const {
0009   Measurement1D unsignedDistance = distance(vtx1, vtx2);
0010   Basic3DVector<float> diff = Basic3DVector<float>(vtx2.position()) - Basic3DVector<float>(vtx1.position());
0011   //   Basic3DVector<float> (vtx2 - vtx1);
0012   if ((momentum.x() * diff.x() + momentum.y() * diff.y() * momentum.z() * diff.z()) < 0)
0013     return Measurement1D(-1.0 * unsignedDistance.value(), unsignedDistance.error());
0014   return unsignedDistance;
0015 }
0016
0017 Measurement1D VertexDistance3D::distance(const GlobalPoint& vtx1Position,
0018                                          const GlobalError& vtx1PositionError,
0019                                          const GlobalPoint& vtx2Position,
0020                                          const GlobalError& vtx2PositionError) const {
0021   AlgebraicSymMatrix33 error = vtx1PositionError.matrix() + vtx2PositionError.matrix();
0022   GlobalVector diff = vtx1Position - vtx2Position;
0023   AlgebraicVector3 vDiff;
0024   vDiff[0] = diff.x();
0025   vDiff[1] = diff.y();
0026   vDiff[2] = diff.z();
0027
0028   double dist = diff.mag();
0029
0030   double err2 = ROOT::Math::Similarity(error, vDiff);
0031   double err = 0.;
0032   if (dist != 0)
0033     err = sqrt(err2) / dist;
0034
0035   return Measurement1D(dist, err);
0036 }
0037
0038 float VertexDistance3D::compatibility(const GlobalPoint& vtx1Position,
0039                                       const GlobalError& vtx1PositionError,
0040                                       const GlobalPoint& vtx2Position,
0041                                       const GlobalError& vtx2PositionError) const {
0042   // error matrix of residuals
0043   AlgebraicSymMatrix33 err1 = vtx1PositionError.matrix();
0044   AlgebraicSymMatrix33 err2 = vtx2PositionError.matrix();
0045   AlgebraicSymMatrix33 error = err1 + err2;
0046   if (error == theNullMatrix)
0047     return FLT_MAX;
0048
0049   // position residuals
0050   GlobalVector diff = vtx2Position - vtx1Position;
0051   AlgebraicVector3 vDiff;
0052   vDiff[0] = diff.x();
0053   vDiff[1] = diff.y();
0054   vDiff[2] = diff.z();
0055
0056   // Invert error matrix of residuals
0057   bool ifail = !error.InvertChol();
0058   if (ifail) {
0059     throw cms::Exception("VertexDistance3D::matrix inversion problem");
0060   }
0061
0062   return ROOT::Math::Similarity(error, vDiff);
0063 }