diff --git a/DataFormats/PatCandidates/interface/PackedCandidate.h b/DataFormats/PatCandidates/interface/PackedCandidate.h index 5118a70f6ff5c..2e448dea9faf8 100644 --- a/DataFormats/PatCandidates/interface/PackedCandidate.h +++ b/DataFormats/PatCandidates/interface/PackedCandidate.h @@ -709,6 +709,11 @@ class PackedCandidate : public reco::Candidate { 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 { if (packedHits_ != 0 || packedLayers_ != 0) { diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index 2f966c2f3da14..ad759af18d153 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_; @@ -178,6 +181,53 @@ float pat::PackedCandidate::dz(const Point &p) const { 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());