Option to make track covariance Pos-def in packedcandidate

cms-sw · Jan 18, 2023 · bae1ece · bae1ece
1 parent 7b7506a
commit bae1ece
Show file tree

Hide file tree

Showing 2 changed files with 55 additions and 0 deletions.
diff --git 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
@@ -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<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) {
+    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());