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