00001 #include <cmath>
00002
00003 #include <Math/Functions.h>
00004 #include <Math/SVector.h>
00005 #include <Math/SMatrix.h>
00006
00007 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
00008 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00009 #include "DataFormats/VertexReco/interface/Vertex.h"
00010
00011 #include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
00012
00013 using namespace reco;
00014
00015 Measurement1D
00016 SecondaryVertex::computeDist2d(const Vertex &pv, const Vertex &sv,
00017 const GlobalVector &direction, bool withPVError)
00018 {
00019 typedef ROOT::Math::SVector<double, 2> SVector2;
00020 typedef ROOT::Math::SMatrix<double, 2, 2,
00021 ROOT::Math::MatRepSym<double, 2> > SMatrixSym2;
00022
00023 SMatrixSym2 cov = sv.covariance().Sub<SMatrixSym2>(0, 0);
00024 if (withPVError)
00025 cov += pv.covariance().Sub<SMatrixSym2>(0, 0);
00026
00027 SVector2 vector(sv.position().X() - pv.position().X(),
00028 sv.position().Y() - pv.position().Y());
00029
00030 double dist = ROOT::Math::Mag(vector);
00031 double error = std::sqrt(ROOT::Math::Similarity(cov, vector)) / dist;
00032
00033 if ((vector[0] * direction.x() +
00034 vector[1] * direction.y()) < 0.0)
00035 dist = -dist;
00036
00037 return Measurement1D(dist, error);
00038 }
00039
00040 Measurement1D
00041 SecondaryVertex::computeDist3d(const Vertex &pv, const Vertex &sv,
00042 const GlobalVector &direction, bool withPVError)
00043 {
00044 typedef ROOT::Math::SVector<double, 3> SVector3;
00045 typedef ROOT::Math::SMatrix<double, 3, 3,
00046 ROOT::Math::MatRepSym<double, 3> > SMatrixSym3;
00047
00048 SMatrixSym3 cov = sv.covariance();
00049 if (withPVError)
00050 cov += pv.covariance();
00051
00052 SVector3 vector(sv.position().X() - pv.position().X(),
00053 sv.position().Y() - pv.position().Y(),
00054 sv.position().Z() - pv.position().Z());
00055
00056 double dist = ROOT::Math::Mag(vector);
00057 double error = std::sqrt(ROOT::Math::Similarity(cov, vector)) / dist;
00058
00059 if ((vector[0] * direction.x() +
00060 vector[1] * direction.y() +
00061 vector[2] * direction.z()) < 0.0)
00062 dist = -dist;
00063
00064 return Measurement1D(dist, error);
00065 }
1.5.4