00001 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
00002 #include <cfloat>
00003
00004
00005 using namespace reco;
00006
00007 Measurement1D
00008 VertexDistance3D::signedDistance(const Vertex& vtx1, const Vertex & vtx2,
00009 const GlobalVector & momentum) const
00010 {
00011 Measurement1D unsignedDistance = distance(vtx1, vtx2);
00012 Basic3DVector<float> diff = Basic3DVector<float> (vtx2.position()) -
00013 Basic3DVector<float> (vtx1.position());
00014
00015 if ((momentum.x()*diff.x() + momentum.y()*diff.y() * momentum.z()*diff.z()) < 0 )
00016 return Measurement1D(-1.0*unsignedDistance.value(),unsignedDistance.error());
00017 return unsignedDistance;
00018 }
00019
00020
00021 Measurement1D
00022 VertexDistance3D::distance(const GlobalPoint & vtx1Position,
00023 const GlobalError & vtx1PositionError,
00024 const GlobalPoint & vtx2Position,
00025 const GlobalError & vtx2PositionError) const
00026 {
00027 AlgebraicSymMatrix error = vtx1PositionError.matrix()
00028 + vtx2PositionError.matrix();
00029 GlobalVector diff = vtx1Position - vtx2Position;
00030 AlgebraicVector vDiff(3);
00031 vDiff[0] = diff.x();
00032 vDiff[1] = diff.y();
00033 vDiff[2] = diff.z();
00034
00035 double dist=diff.mag();
00036
00037 double err2 = error.similarity(vDiff);
00038 double err = 0.;
00039 if (dist != 0) err = sqrt(err2)/dist;
00040
00041 return Measurement1D(dist,err);
00042 }
00043
00044
00045
00046 float
00047 VertexDistance3D::compatibility(const GlobalPoint & vtx1Position,
00048 const GlobalError & vtx1PositionError,
00049 const GlobalPoint & vtx2Position,
00050 const GlobalError & vtx2PositionError) const
00051 {
00052
00053 AlgebraicSymMatrix err1 = vtx1PositionError.matrix();
00054 AlgebraicSymMatrix err2 = vtx2PositionError.matrix();
00055 AlgebraicSymMatrix error = err1 + err2;
00056 if (error == theNullMatrix) return FLT_MAX;
00057
00058
00059 GlobalVector diff = vtx2Position - vtx1Position;
00060 AlgebraicVector vDiff(3);
00061 vDiff[0] = diff.x();
00062 vDiff[1] = diff.y();
00063 vDiff[2] = diff.z();
00064
00065
00066 int ifail;
00067 error.invert(ifail);
00068 if (ifail != 0) {
00069 throw cms::Exception("VertexDistance3D::matrix inversion problem");
00070 }
00071
00072 return error.similarity(vDiff);
00073 }