VertexTools/src/VertexDistanceXY.cc

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