diff --git a/DataFormats/PatCandidates/interface/PackedCandidate.h b/DataFormats/PatCandidates/interface/PackedCandidate.h index 47e79e956dfe7..a8fc0ab95fd35 100644 --- a/DataFormats/PatCandidates/interface/PackedCandidate.h +++ b/DataFormats/PatCandidates/interface/PackedCandidate.h @@ -787,6 +787,10 @@ namespace pat { unpackTrk(); return *track_; } + /// Return reference to a pseudo track made with candidate kinematics, + /// parameterized error for eta,phi,pt and full IP covariance + /// and the coviriance matrix is forced to be positive definite according to BPH recommandations + virtual const reco::Track pseudoPosDefTrack() const; /// return a pointer to the track if present. otherwise, return a null pointer const reco::Track *bestTrack() const override { diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index e58f05d83edd2..0271b99782a36 100644 --- a/DataFormats/PatCandidates/src/PackedCandidate.cc +++ b/DataFormats/PatCandidates/src/PackedCandidate.cc @@ -5,6 +5,9 @@ #include "DataFormats/SiStripDetId/interface/StripSubdetector.h" #include "DataFormats/Math/interface/liblogintpack.h" + +#include "TMatrixDSym.h" +#include "TVectorD.h" using namespace logintpack; CovarianceParameterization pat::PackedCandidate::covarianceParameterization_; @@ -105,6 +108,7 @@ void pat::PackedCandidate::unpackCovariance() const { unpackCovarianceElement(*m, packedCovariance_.dxydz, 3, 4); unpackCovarianceElement(*m, packedCovariance_.dlambdadz, 1, 4); unpackCovarianceElement(*m, packedCovariance_.dphidxy, 2, 3); + reco::TrackBase::CovarianceMatrix *expected = nullptr; if (m_.compare_exchange_strong(expected, m.get())) { m.release(); @@ -161,6 +165,53 @@ float pat::PackedCandidate::dz(const Point &p) const { ((vertex_.load()->X() - p.X()) * std::cos(phi) + (vertex_.load()->Y() - p.Y()) * std::sin(phi)) * pzpt; } +const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const { + // perform the regular unpacking of the track + if (!track_) + unpackTrk(); + + //calculate the determinant and verify positivity + double det = 0; + bool notPosDef = (!(*m_).Sub(0, 0).Det(det) || det < 0) || + (!(*m_).Sub(0, 0).Det(det) || det < 0) || + (!(*m_).Sub(0, 0).Det(det) || det < 0) || (!(*m_).Det(det) || det < 0); + + if (notPosDef) { + reco::TrackBase::CovarianceMatrix m(*m_); + //if not positive-definite, alter values to allow for pos-def + TMatrixDSym eigenCov(5); + for (int i = 0; i < 5; i++) { + for (int j = 0; j < 5; j++) { + if (std::isnan((m)(i, j)) || std::isinf((m)(i, j))) + eigenCov(i, j) = 1e-6; + else + eigenCov(i, j) = (m)(i, j); + } + } + TVectorD eigenValues(5); + eigenCov.EigenVectors(eigenValues); + double minEigenValue = eigenValues.Min(); + double delta = 1e-6; + if (minEigenValue < 0) { + for (int i = 0; i < 5; i++) + m(i, i) += delta - minEigenValue; + } + + // make a track object with pos def covariance matrix + return reco::Track(normalizedChi2_ * (*track_).ndof(), + (*track_).ndof(), + *vertex_, + (*track_).momentum(), + (*track_).charge(), + m, + reco::TrackBase::undefAlgorithm, + reco::TrackBase::loose); + } else { + // just return a copy of the unpacked track + return reco::Track(*track_); + } +} + void pat::PackedCandidate::unpackTrk() const { maybeUnpackBoth(); math::RhoEtaPhiVector p3(ptTrk(), etaAtVtx(), phiAtVtx());