From ee60abe1d3b49d1fb9d1c16022e8ee653c582ea1 Mon Sep 17 00:00:00 2001
From: jean-roch <vlimant@cern.ch>
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<reco::TrackBase::CovarianceMatrix>();
@@ -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"<<std::endl;
+ //calculate the determinant and verify positivity
+ double det_(0);
+ bool computed_ = m->Det(det_);
+ // std::cout<<"during unpacking, the determinant is: "<<det_<<std::endl;
+ if (computed_ && det_<0){
+ //if not positive-definite, alter values to allow for pos-def
+ TMatrixDSym COV(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))){
+ COV(i, j) = 1e-6;
+ }else{
+ COV(i,j) = (*m)(i,j);
+ }
+ }
+ }
+ TVectorD EIG(5);
+ COV.EigenVectors(EIG);
+ double MIN_EIG = 100000;
+ double DELTA = 1e-6;
+ for (int i = 0; i < 5; i++){
+ if (EIG(i) < MIN_EIG)
+ MIN_EIG = EIG(i);
+ }
+ if (MIN_EIG<0) {
+ for (int i = 0; i < 5; i++)
+ (*m)(i,i) += DELTA - MIN_EIG;
+ }
+
+ computed_ = m->Det(det_);
+ // std::cout<<" the determinant of the corrected covariance matrix is: "<<det_<<std::endl;
+ }
+ }
+
reco::TrackBase::CovarianceMatrix *expected = nullptr;
if (m_.compare_exchange_strong(expected, m.get())) {
m.release();
@@ -161,10 +201,10 @@ 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;
}
-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 <vlimant@cern.ch>
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"<<std::endl;
//calculate the determinant and verify positivity
double det_(0);
bool computed_ = m->Det(det_);
// std::cout<<"during unpacking, the determinant is: "<<det_<<std::endl;
- if (computed_ && det_<0){
- //if not positive-definite, alter values to allow for pos-def
- TMatrixDSym COV(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))){
- COV(i, j) = 1e-6;
- }else{
- COV(i,j) = (*m)(i,j);
- }
- }
- }
- TVectorD EIG(5);
- COV.EigenVectors(EIG);
- double MIN_EIG = 100000;
- double DELTA = 1e-6;
- for (int i = 0; i < 5; i++){
- if (EIG(i) < MIN_EIG)
- MIN_EIG = EIG(i);
- }
- if (MIN_EIG<0) {
- for (int i = 0; i < 5; i++)
- (*m)(i,i) += DELTA - MIN_EIG;
- }
-
- computed_ = m->Det(det_);
- // std::cout<<" the determinant of the corrected covariance matrix is: "<<det_<<std::endl;
+ if (computed_ && det_ < 0) {
+ //if not positive-definite, alter values to allow for pos-def
+ TMatrixDSym COV(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))) {
+ COV(i, j) = 1e-6;
+ } else {
+ COV(i, j) = (*m)(i, j);
+ }
+ }
+ }
+ TVectorD EIG(5);
+ COV.EigenVectors(EIG);
+ double MIN_EIG = 100000;
+ double DELTA = 1e-6;
+ for (int i = 0; i < 5; i++) {
+ if (EIG(i) < MIN_EIG)
+ MIN_EIG = EIG(i);
+ }
+ if (MIN_EIG < 0) {
+ for (int i = 0; i < 5; i++)
+ (*m)(i, i) += DELTA - MIN_EIG;
+ }
+
+ computed_ = m->Det(det_);
+ // std::cout<<" the determinant of the corrected covariance matrix is: "<<det_<<std::endl;
}
}
From f4b8f113b413d4bc3b143a11d4663eac4024574e Mon Sep 17 00:00:00 2001
From: jean-roch <vlimant@cern.ch>
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 <vlimant@cern.ch>
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"<<std::endl;
//calculate the determinant and verify positivity
- double det_(0);
- bool computed_ = m->Det(det_);
- // std::cout<<"during unpacking, the determinant is: "<<det_<<std::endl;
- if (computed_ && det_ < 0) {
+ double det = 0;
+ bool notPosDef = (!m->Sub<AlgebraicSymMatrix22>(0, 0).Det(det) || det<0)
+ || (!m->Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det<0)
+ || (!m->Sub<AlgebraicSymMatrix44>(0, 0).Det(det) || det<0)
+ || (!m->Det(det) || det<0);
+ // std::cout<<"during unpacking, the determinant is: "<<det<<std::endl;
+ if (notPosDef) {
//if not positive-definite, alter values to allow for pos-def
- TMatrixDSym COV(5);
+ 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))) {
- COV(i, j) = 1e-6;
- } else {
- COV(i, j) = (*m)(i, 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 EIG(5);
- COV.EigenVectors(EIG);
- double MIN_EIG = 100000;
- double DELTA = 1e-6;
- for (int i = 0; i < 5; i++) {
- if (EIG(i) < MIN_EIG)
- MIN_EIG = EIG(i);
- }
- if (MIN_EIG < 0) {
+ 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 - MIN_EIG;
+ (*m)(i, i) += delta - minEigenValue;
}
- computed_ = m->Det(det_);
- // std::cout<<" the determinant of the corrected covariance matrix is: "<<det_<<std::endl;
+ // computed_ = m->Det(det);
+ // std::cout<<" the determinant of the corrected covariance matrix is: "<<det<<std::endl;
}
}
From 23ebfb9697db139d44f05539b6453748ba30b259 Mon Sep 17 00:00:00 2001
From: jean-roch <vlimant@cern.ch>
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"<<std::endl;
//calculate the determinant and verify positivity
double det = 0;
- bool notPosDef = (!m->Sub<AlgebraicSymMatrix22>(0, 0).Det(det) || det<0)
- || (!m->Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det<0)
- || (!m->Sub<AlgebraicSymMatrix44>(0, 0).Det(det) || det<0)
- || (!m->Det(det) || det<0);
+ bool notPosDef = (!m->Sub<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
+ (!m->Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det < 0) ||
+ (!m->Sub<AlgebraicSymMatrix44>(0, 0).Det(det) || det < 0) || (!m->Det(det) || det < 0);
// std::cout<<"during unpacking, the determinant is: "<<det<<std::endl;
if (notPosDef) {
//if not positive-definite, alter values to allow for pos-def
@@ -125,7 +124,7 @@ void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const {
for (int j = 0; j < 5; j++) {
if (std::isnan((*m)(i, j)) || std::isinf((*m)(i, j)))
eigenCov(i, j) = 1e-6;
- else
+ else
eigenCov(i, j) = (*m)(i, j);
}
}
@@ -138,7 +137,7 @@ void pat::PackedCandidate::unpackCovariance(bool forcePosDef) const {
(*m)(i, i) += delta - minEigenValue;
}
- // computed_ = m->Det(det);
+ // computed_ = m->Det(det);
// std::cout<<" the determinant of the corrected covariance matrix is: "<<det<<std::endl;
}
}
From 1622d28c46f3f9cb0d15ecbcec9931c41d84a22d Mon Sep 17 00:00:00 2001
From: jean-roch <vlimant@cern.ch>
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<reco::TrackBase::CovarianceMatrix>();
@@ -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"<<std::endl;
- //calculate the determinant and verify positivity
- double det = 0;
- bool notPosDef = (!m->Sub<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
- (!m->Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det < 0) ||
- (!m->Sub<AlgebraicSymMatrix44>(0, 0).Det(det) || det < 0) || (!m->Det(det) || det < 0);
- // std::cout<<"during unpacking, the determinant is: "<<det<<std::endl;
- if (notPosDef) {
- //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;
- }
-
- // computed_ = m->Det(det);
- // std::cout<<" the determinant of the corrected covariance matrix is: "<<det<<std::endl;
- }
- }
-
reco::TrackBase::CovarianceMatrix *expected = nullptr;
if (m_.compare_exchange_strong(expected, m.get())) {
m.release();
@@ -198,10 +165,64 @@ 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;
}
-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"<<std::endl;
+ //calculate the determinant and verify positivity
+ double det = 0;
+ bool notPosDef = (!(*m_).Sub<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
+ (!(*m_).Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det < 0) ||
+ (!(*m_).Sub<AlgebraicSymMatrix44>(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<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
+ (!m.Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det < 0) ||
+ (!m.Sub<AlgebraicSymMatrix44>(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 <vlimant@cern.ch>
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<AlgebraicSymMatrix44>(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<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
- (!m.Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det < 0) ||
- (!m.Sub<AlgebraicSymMatrix44>(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<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
+ // (!m.Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det < 0) ||
+ // (!m.Sub<AlgebraicSymMatrix44>(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 <vlimant@cern.ch>
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"<<std::endl;
//calculate the determinant and verify positivity
double det = 0;
bool notPosDef = (!(*m_).Sub<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
@@ -179,7 +177,6 @@ const reco::Track pat::PackedCandidate::pseudoPosDefTrack() const {
(!(*m_).Sub<AlgebraicSymMatrix44>(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<AlgebraicSymMatrix22>(0, 0).Det(det) || det < 0) ||
- // (!m.Sub<AlgebraicSymMatrix33>(0, 0).Det(det) || det < 0) ||
- // (!m.Sub<AlgebraicSymMatrix44>(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(),