From ee60abe1d3b49d1fb9d1c16022e8ee653c582ea1 Mon Sep 17 00:00:00 2001 From: jean-roch Date: Tue, 4 Oct 2022 13:40:50 +0200 Subject: [PATCH 1/8] introduce a forcePosDef option in unpacking the pseudo track of PackedCandidate --- .../PatCandidates/interface/PackedCandidate.h | 16 ++++--- .../PatCandidates/src/PackedCandidate.cc | 46 +++++++++++++++++-- 2 files changed, 53 insertions(+), 9 deletions(-) diff --git a/DataFormats/PatCandidates/interface/PackedCandidate.h b/DataFormats/PatCandidates/interface/PackedCandidate.h index 47e79e956dfe7..625b954f739c6 100644 --- a/DataFormats/PatCandidates/interface/PackedCandidate.h +++ b/DataFormats/PatCandidates/interface/PackedCandidate.h @@ -782,11 +782,15 @@ namespace pat { /// Return reference to a pseudo track made with candidate kinematics, /// parameterized error for eta,phi,pt and full IP covariance - virtual const reco::Track &pseudoTrack() const { + virtual const reco::Track &pseudoTrack(bool forcePosDef=false) const { if (!track_) - unpackTrk(); + unpackTrk(forcePosDef); return *track_; } + virtual void resetPseudoTrack() const { + delete m_.exchange(nullptr); + delete track_.exchange(nullptr); + } /// return a pointer to the track if present. otherwise, return a null pointer const reco::Track *bestTrack() const override { @@ -1019,7 +1023,7 @@ namespace pat { void packVtx(bool unpackAfterwards = true); void unpackVtx() const; void packCovariance(const reco::TrackBase::CovarianceMatrix &m, bool unpackAfterwards = true); - void unpackCovariance() const; + void unpackCovariance(bool forcePosDef = true) const; void maybeUnpackBoth() const { if (!p4c_) unpack(); @@ -1030,9 +1034,9 @@ namespace pat { if (!track_) unpackTrk(); } - void maybeUnpackCovariance() const { + void maybeUnpackCovariance(bool forcePosDef=false) const { if (!m_) - unpackCovariance(); + unpackCovariance(forcePosDef); } void packBoth() { pack(false); @@ -1044,7 +1048,7 @@ namespace pat { unpackVtx(); } // do it this way, so that we don't loose precision on the angles before // computing dxy,dz - void unpackTrk() const; + void unpackTrk(bool forcePosDef=false) const; uint8_t packedPuppiweight_; int8_t packedPuppiweightNoLepDiff_; // storing the DIFFERENCE of (all - "no diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index e58f05d83edd2..a4fa3365c855e 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_; @@ -89,7 +92,7 @@ void pat::PackedCandidate::packCovariance(const reco::TrackBase::CovarianceMatri unpackCovariance(); } -void pat::PackedCandidate::unpackCovariance() const { +void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const { const CovarianceParameterization &p = covarianceParameterization(); if (p.isValid()) { auto m = std::make_unique(); @@ -105,6 +108,43 @@ void pat::PackedCandidate::unpackCovariance() const { unpackCovarianceElement(*m, packedCovariance_.dxydz, 3, 4); unpackCovarianceElement(*m, packedCovariance_.dlambdadz, 1, 4); unpackCovarianceElement(*m, packedCovariance_.dphidxy, 2, 3); + + if (forcePosDef){ + // std::cout<<"unpacking enforcing positive-definite cov matrix"<Det(det_); + // std::cout<<"during unpacking, the determinant is: "<Det(det_); + // std::cout<<" the determinant of the corrected covariance matrix is: "<X() - p.X()) * std::cos(phi) + (vertex_.load()->Y() - p.Y()) * std::sin(phi)) * pzpt; } -void pat::PackedCandidate::unpackTrk() const { +void pat::PackedCandidate::unpackTrk(bool forcePosDef) const { maybeUnpackBoth(); math::RhoEtaPhiVector p3(ptTrk(), etaAtVtx(), phiAtVtx()); - maybeUnpackCovariance(); + maybeUnpackCovariance(forcePosDef); int numberOfStripLayers = stripLayersWithMeasurement(), numberOfPixelLayers = pixelLayersWithMeasurement(); int numberOfPixelHits = this->numberOfPixelHits(); int numberOfHits = this->numberOfHits(); From 0fb9dea085240afd926868397ce9dc1ba254dc4b Mon Sep 17 00:00:00 2001 From: jean-roch Date: Tue, 4 Oct 2022 14:46:35 +0200 Subject: [PATCH 2/8] code formatting --- .../PatCandidates/interface/PackedCandidate.h | 6 +- .../PatCandidates/src/PackedCandidate.cc | 56 +++++++++---------- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/DataFormats/PatCandidates/interface/PackedCandidate.h b/DataFormats/PatCandidates/interface/PackedCandidate.h index 625b954f739c6..c2a0c709f7e66 100644 --- a/DataFormats/PatCandidates/interface/PackedCandidate.h +++ b/DataFormats/PatCandidates/interface/PackedCandidate.h @@ -782,7 +782,7 @@ namespace pat { /// Return reference to a pseudo track made with candidate kinematics, /// parameterized error for eta,phi,pt and full IP covariance - virtual const reco::Track &pseudoTrack(bool forcePosDef=false) const { + virtual const reco::Track &pseudoTrack(bool forcePosDef = false) const { if (!track_) unpackTrk(forcePosDef); return *track_; @@ -1034,7 +1034,7 @@ namespace pat { if (!track_) unpackTrk(); } - void maybeUnpackCovariance(bool forcePosDef=false) const { + void maybeUnpackCovariance(bool forcePosDef = false) const { if (!m_) unpackCovariance(forcePosDef); } @@ -1048,7 +1048,7 @@ namespace pat { unpackVtx(); } // do it this way, so that we don't loose precision on the angles before // computing dxy,dz - void unpackTrk(bool forcePosDef=false) const; + void unpackTrk(bool forcePosDef = false) const; uint8_t packedPuppiweight_; int8_t packedPuppiweightNoLepDiff_; // storing the DIFFERENCE of (all - "no diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index a4fa3365c855e..c40dedebced11 100644 --- a/DataFormats/PatCandidates/src/PackedCandidate.cc +++ b/DataFormats/PatCandidates/src/PackedCandidate.cc @@ -109,39 +109,39 @@ void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const { unpackCovarianceElement(*m, packedCovariance_.dlambdadz, 1, 4); unpackCovarianceElement(*m, packedCovariance_.dphidxy, 2, 3); - if (forcePosDef){ + if (forcePosDef) { // std::cout<<"unpacking enforcing positive-definite cov matrix"<Det(det_); // std::cout<<"during unpacking, the determinant is: "<Det(det_); - // std::cout<<" the determinant of the corrected covariance matrix is: "<Det(det_); + // std::cout<<" the determinant of the corrected covariance matrix is: "< Date: Tue, 4 Oct 2022 17:03:18 +0200 Subject: [PATCH 3/8] the default is meant to be false --- DataFormats/PatCandidates/interface/PackedCandidate.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DataFormats/PatCandidates/interface/PackedCandidate.h b/DataFormats/PatCandidates/interface/PackedCandidate.h index c2a0c709f7e66..91774d854d4ef 100644 --- a/DataFormats/PatCandidates/interface/PackedCandidate.h +++ b/DataFormats/PatCandidates/interface/PackedCandidate.h @@ -1023,7 +1023,7 @@ namespace pat { void packVtx(bool unpackAfterwards = true); void unpackVtx() const; void packCovariance(const reco::TrackBase::CovarianceMatrix &m, bool unpackAfterwards = true); - void unpackCovariance(bool forcePosDef = true) const; + void unpackCovariance(bool forcePosDef = false) const; void maybeUnpackBoth() const { if (!p4c_) unpack(); From 5a743ccb8d283c64fed76dab1946d1c1462fdffd Mon Sep 17 00:00:00 2001 From: jean-roch Date: Tue, 4 Oct 2022 17:06:08 +0200 Subject: [PATCH 4/8] code naming. Sylvester test for posdef --- .../PatCandidates/src/PackedCandidate.cc | 42 +++++++++---------- 1 file changed, 20 insertions(+), 22 deletions(-) diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index c40dedebced11..c4e7ac3832919 100644 --- a/DataFormats/PatCandidates/src/PackedCandidate.cc +++ b/DataFormats/PatCandidates/src/PackedCandidate.cc @@ -112,36 +112,34 @@ void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const { if (forcePosDef) { // std::cout<<"unpacking enforcing positive-definite cov matrix"<Det(det_); - // std::cout<<"during unpacking, the determinant is: "<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); + // std::cout<<"during unpacking, the determinant is: "<Det(det_); - // std::cout<<" the determinant of the corrected covariance matrix is: "<Det(det); + // std::cout<<" the determinant of the corrected covariance matrix is: "< Date: Tue, 4 Oct 2022 17:21:58 +0200 Subject: [PATCH 5/8] code formatting --- DataFormats/PatCandidates/src/PackedCandidate.cc | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index c4e7ac3832919..d8d79cecd6f67 100644 --- a/DataFormats/PatCandidates/src/PackedCandidate.cc +++ b/DataFormats/PatCandidates/src/PackedCandidate.cc @@ -113,10 +113,9 @@ void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const { // std::cout<<"unpacking enforcing positive-definite cov matrix"<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); + 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); // std::cout<<"during unpacking, the determinant is: "<Det(det); + // computed_ = m->Det(det); // std::cout<<" the determinant of the corrected covariance matrix is: "< Date: Fri, 7 Oct 2022 11:23:24 +0200 Subject: [PATCH 6/8] moving to separate copy of pseudo track --- .../PatCandidates/interface/PackedCandidate.h | 20 ++-- .../PatCandidates/src/PackedCandidate.cc | 93 ++++++++++++------- 2 files changed, 67 insertions(+), 46 deletions(-) diff --git a/DataFormats/PatCandidates/interface/PackedCandidate.h b/DataFormats/PatCandidates/interface/PackedCandidate.h index 91774d854d4ef..a8fc0ab95fd35 100644 --- a/DataFormats/PatCandidates/interface/PackedCandidate.h +++ b/DataFormats/PatCandidates/interface/PackedCandidate.h @@ -782,15 +782,15 @@ namespace pat { /// Return reference to a pseudo track made with candidate kinematics, /// parameterized error for eta,phi,pt and full IP covariance - virtual const reco::Track &pseudoTrack(bool forcePosDef = false) const { + virtual const reco::Track &pseudoTrack() const { if (!track_) - unpackTrk(forcePosDef); + unpackTrk(); return *track_; } - virtual void resetPseudoTrack() const { - delete m_.exchange(nullptr); - delete track_.exchange(nullptr); - } + /// 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 { @@ -1023,7 +1023,7 @@ namespace pat { void packVtx(bool unpackAfterwards = true); void unpackVtx() const; void packCovariance(const reco::TrackBase::CovarianceMatrix &m, bool unpackAfterwards = true); - void unpackCovariance(bool forcePosDef = false) const; + void unpackCovariance() const; void maybeUnpackBoth() const { if (!p4c_) unpack(); @@ -1034,9 +1034,9 @@ namespace pat { if (!track_) unpackTrk(); } - void maybeUnpackCovariance(bool forcePosDef = false) const { + void maybeUnpackCovariance() const { if (!m_) - unpackCovariance(forcePosDef); + unpackCovariance(); } void packBoth() { pack(false); @@ -1048,7 +1048,7 @@ namespace pat { unpackVtx(); } // do it this way, so that we don't loose precision on the angles before // computing dxy,dz - void unpackTrk(bool forcePosDef = false) const; + void unpackTrk() const; uint8_t packedPuppiweight_; int8_t packedPuppiweightNoLepDiff_; // storing the DIFFERENCE of (all - "no diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index d8d79cecd6f67..46f32f37d4b48 100644 --- a/DataFormats/PatCandidates/src/PackedCandidate.cc +++ b/DataFormats/PatCandidates/src/PackedCandidate.cc @@ -92,7 +92,7 @@ void pat::PackedCandidate::packCovariance(const reco::TrackBase::CovarianceMatri unpackCovariance(); } -void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const { +void pat::PackedCandidate::unpackCovariance() const { const CovarianceParameterization &p = covarianceParameterization(); if (p.isValid()) { auto m = std::make_unique(); @@ -109,39 +109,6 @@ void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const { unpackCovarianceElement(*m, packedCovariance_.dlambdadz, 1, 4); unpackCovarianceElement(*m, packedCovariance_.dphidxy, 2, 3); - if (forcePosDef) { - // std::cout<<"unpacking enforcing positive-definite cov matrix"<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); - // std::cout<<"during unpacking, the determinant is: "<Det(det); - // std::cout<<" the determinant of the corrected covariance matrix is: "<X() - p.X()) * std::cos(phi) + (vertex_.load()->Y() - p.Y()) * std::sin(phi)) * pzpt; } -void pat::PackedCandidate::unpackTrk(bool forcePosDef) const { +const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const { + //if (posdeftrack_) return *posdeftrack_; + // perform the regular unpacking of the track + if (!track_) + unpackTrk(); + + // std::cout<<"unpacking enforcing positive-definite cov matrix"<(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) { + std::cout << "during unpacking, the determinant is: " << det << std::endl; + 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; + } + + 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); + std::cout << " the determinant of the corrected covariance matrix is: " << det << std::endl; + // 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()); - maybeUnpackCovariance(forcePosDef); + maybeUnpackCovariance(); int numberOfStripLayers = stripLayersWithMeasurement(), numberOfPixelLayers = pixelLayersWithMeasurement(); int numberOfPixelHits = this->numberOfPixelHits(); int numberOfHits = this->numberOfHits(); From 3f8ecfd8dbd3e9796959e61da24e1df0c0f55746 Mon Sep 17 00:00:00 2001 From: jean-roch Date: Fri, 7 Oct 2022 14:05:18 +0200 Subject: [PATCH 7/8] remove unused var and cout --- DataFormats/PatCandidates/src/PackedCandidate.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index 46f32f37d4b48..0d7ee67a616e0 100644 --- a/DataFormats/PatCandidates/src/PackedCandidate.cc +++ b/DataFormats/PatCandidates/src/PackedCandidate.cc @@ -179,7 +179,7 @@ const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const { (!(*m_).Sub(0, 0).Det(det) || det < 0) || (!(*m_).Det(det) || det < 0); if (notPosDef) { - std::cout << "during unpacking, the determinant is: " << det << std::endl; + // std::cout << "during unpacking, the determinant is: " << det << std::endl; reco::TrackBase::CovarianceMatrix m(*m_); //if not positive-definite, alter values to allow for pos-def TMatrixDSym eigenCov(5); @@ -200,10 +200,10 @@ const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const { m(i, i) += delta - minEigenValue; } - 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); - std::cout << " the determinant of the corrected covariance matrix is: " << det << std::endl; + // 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); + // std::cout << " the determinant of the corrected covariance matrix is: " << det << std::endl; // make a track object with pos def covariance matrix return reco::Track(normalizedChi2_ * (*track_).ndof(), (*track_).ndof(), From 5c28160dc3d95197e943bd2ec9f7fb89ff5bf4e0 Mon Sep 17 00:00:00 2001 From: jean-roch Date: Fri, 14 Oct 2022 10:37:36 +0200 Subject: [PATCH 8/8] removed commented-out statements --- DataFormats/PatCandidates/src/PackedCandidate.cc | 7 ------- 1 file changed, 7 deletions(-) diff --git a/DataFormats/PatCandidates/src/PackedCandidate.cc b/DataFormats/PatCandidates/src/PackedCandidate.cc index 0d7ee67a616e0..0271b99782a36 100644 --- a/DataFormats/PatCandidates/src/PackedCandidate.cc +++ b/DataFormats/PatCandidates/src/PackedCandidate.cc @@ -166,12 +166,10 @@ float pat::PackedCandidate::dz(const Point &p) const { } const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const { - //if (posdeftrack_) return *posdeftrack_; // perform the regular unpacking of the track if (!track_) unpackTrk(); - // std::cout<<"unpacking enforcing positive-definite cov matrix"<(0, 0).Det(det) || det < 0) || @@ -179,7 +177,6 @@ const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const { (!(*m_).Sub(0, 0).Det(det) || det < 0) || (!(*m_).Det(det) || det < 0); if (notPosDef) { - // std::cout << "during unpacking, the determinant is: " << det << std::endl; reco::TrackBase::CovarianceMatrix m(*m_); //if not positive-definite, alter values to allow for pos-def TMatrixDSym eigenCov(5); @@ -200,10 +197,6 @@ const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const { m(i, i) += delta - minEigenValue; } - // 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); - // std::cout << " the determinant of the corrected covariance matrix is: " << det << std::endl; // make a track object with pos def covariance matrix return reco::Track(normalizedChi2_ * (*track_).ndof(), (*track_).ndof(),