From 59d8fc1d8865a1dcbe6e99a249339a31030a5294 Mon Sep 17 00:00:00 2001
From: Lucy Lewitt <lucy.lewitt@cern.ch>
Date: Tue, 28 May 2024 15:13:36 +0200
Subject: [PATCH] Add a dedicated chi2 method.
Add a dedicated chi2 method.
---
.../EventPrimitivesCovarianceHelpers.h | 7 ++++
.../share/test_CovarianceHelpers.ref | 30 +++++++++++++++
.../test/test_CovarianceHelpers.cxx | 37 +++++++++++++++++++
.../src/GsfMeasurementUpdator.cxx | 18 ++++-----
.../TrkMeasurementUpdator/KalmanUpdator.h | 3 +-
.../src/KalmanUpdator.cxx | 6 +--
.../src/KalmanUpdatorSMatrix.cxx | 5 ++-
.../src/KalmanWeightUpdator.cxx | 9 ++---
.../src/FastVertexFitter.cxx | 4 +-
9 files changed, 95 insertions(+), 24 deletions(-)
diff --git a/Event/EventPrimitives/EventPrimitives/EventPrimitivesCovarianceHelpers.h b/Event/EventPrimitives/EventPrimitives/EventPrimitivesCovarianceHelpers.h
index eb1f3cd6ff00..d02ddab1cbba 100644
--- a/Event/EventPrimitives/EventPrimitives/EventPrimitivesCovarianceHelpers.h
+++ b/Event/EventPrimitives/EventPrimitives/EventPrimitivesCovarianceHelpers.h
@@ -214,6 +214,13 @@ inline bool isPositiveDefiniteSlow(const Amg::MatrixX& mat) {
}
return true;
}
+// Chi squared statistic of a linear regression as defined in Frühwirth,
+// section 3.2.1, equation 3.19. Precision is defined as the inverse of
+// of the covariance.
+template <typename T, typename U>
+inline double chi2(const T& precision, const U& residual, const int sign = 1) {
+ return sign * residual.transpose() * precision * residual;
+}
} // namespace Amg
diff --git a/Event/EventPrimitives/share/test_CovarianceHelpers.ref b/Event/EventPrimitives/share/test_CovarianceHelpers.ref
index d6535642fb42..6b485163da14 100644
--- a/Event/EventPrimitives/share/test_CovarianceHelpers.ref
+++ b/Event/EventPrimitives/share/test_CovarianceHelpers.ref
@@ -12,6 +12,12 @@ Cholesky isPositiveSemiDefinite 0
Cholesky isPositiveDefinite 0
All diagonal Elements are >=0 1
All diagonal Elements are >0 1
+ 50977.5
+ -3154.7
+ 191.7
+-3.12793
+ 3.7e-05
+ChiSquared is 1.33992e+14
Testing dynamic Matrix A
50977.5 -3154.7 191.7 -3.12793 3.7e-05
@@ -25,6 +31,12 @@ Cholesky isPositiveSemiDefinite 0
Cholesky isPositiveDefinite 0
All diagonal Elements are >=0 1
All diagonal Elements are >0 1
+ 50977.5
+ -3154.7
+ 191.7
+-3.12793
+ 3.7e-05
+ChiSquared is 1.33992e+14
Testing Matrix B
36 0 -0.0222 -0.0001 0
@@ -38,6 +50,12 @@ Cholesky isPositiveSemiDefinite 0
Cholesky isPositiveDefinite 0
All diagonal Elements are >=0 0
All diagonal Elements are >0 0
+ 36
+ 0
+-0.0222
+-0.0001
+ 0
+ChiSquared is 43584
Testing Zero Matrix
0 0 0 0 0
@@ -51,6 +69,12 @@ Cholesky isPositiveSemiDefinite 1
Cholesky isPositiveDefinite 0
All diagonal Elements are >=0 1
All diagonal Elements are >0 0
+0
+0
+0
+0
+0
+ChiSquared is 0
Testing Dynamic Zero Matrix
0 0 0 0 0
@@ -64,3 +88,9 @@ Cholesky isPositiveSemiDefinite 1
Cholesky isPositiveDefinite 0
All diagonal Elements are >=0 1
All diagonal Elements are >0 0
+0
+0
+0
+0
+0
+ChiSquared is 0
diff --git a/Event/EventPrimitives/test/test_CovarianceHelpers.cxx b/Event/EventPrimitives/test/test_CovarianceHelpers.cxx
index 5edb72353111..0dfd99fd3f65 100644
--- a/Event/EventPrimitives/test/test_CovarianceHelpers.cxx
+++ b/Event/EventPrimitives/test/test_CovarianceHelpers.cxx
@@ -36,6 +36,13 @@ int main() {
<< Amg::hasPositiveOrZeroDiagElems(A) << '\n';
std::cout << "All diagonal Elements are >0 " << Amg::hasPositiveDiagElems(A)
<< '\n';
+
+ AmgVector(5) V;
+ V << 50977.455023, -3154.699348, 191.699597, -3.127933, 0.000037;
+
+ std::cout << V << '\n';
+ std::cout << "ChiSquared is " << Amg::chi2(A, V)
+ << '\n';
}
std::cout << '\n' << "Testing dynamic Matrix A" << '\n';
@@ -60,6 +67,14 @@ int main() {
<< Amg::hasPositiveOrZeroDiagElems(A) << '\n';
std::cout << "All diagonal Elements are >0 " << Amg::hasPositiveDiagElems(A)
<< '\n';
+
+ Amg::VectorX V;
+ V.resize(5);
+ V << 50977.455023, -3154.699348, 191.699597, -3.127933, 0.000037;
+
+ std::cout << V << '\n';
+ std::cout << "ChiSquared is " << Amg::chi2(A, V)
+ << '\n';
}
std::cout << '\n' << "Testing Matrix B" << '\n';
@@ -82,6 +97,13 @@ int main() {
<< Amg::hasPositiveOrZeroDiagElems(B) << '\n';
std::cout << "All diagonal Elements are >0 " << Amg::hasPositiveDiagElems(B)
<< '\n';
+
+ AmgVector(5) V;
+ V << 36.0000, 0.0000, -0.0222, -0.0001, 0.0000;
+
+ std::cout << V << '\n';
+ std::cout << "ChiSquared is " << Amg::chi2(B, V)
+ << '\n';
}
std::cout << '\n' << "Testing Zero Matrix" << '\n';
@@ -101,6 +123,13 @@ int main() {
<< Amg::hasPositiveOrZeroDiagElems(zero) << '\n';
std::cout << "All diagonal Elements are >0 "
<< Amg::hasPositiveDiagElems(zero) << '\n';
+
+ AmgVector(5) zeroV;
+ zeroV.setZero();
+
+ std::cout << zeroV << '\n';
+ std::cout << "ChiSquared is " << Amg::chi2(zero, zeroV)
+ << '\n';
}
std::cout << '\n' << "Testing Dynamic Zero Matrix" << '\n';
@@ -121,5 +150,13 @@ int main() {
<< Amg::hasPositiveOrZeroDiagElems(zero) << '\n';
std::cout << "All diagonal Elements are >0 "
<< Amg::hasPositiveDiagElems(zero) << '\n';
+
+ Amg::VectorX zeroV;
+ zeroV.resize(5);
+ zeroV.setZero();
+
+ std::cout << zeroV << '\n';
+ std::cout << "ChiSquared is " << Amg::chi2(zero, zeroV)
+ << '\n';
}
}
diff --git a/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/GsfMeasurementUpdator.cxx b/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/GsfMeasurementUpdator.cxx
index 488bc9de8d2e..a02d9f7697f0 100644
--- a/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/GsfMeasurementUpdator.cxx
+++ b/Tracking/TrkFitter/TrkGaussianSumFilterUtils/src/GsfMeasurementUpdator.cxx
@@ -8,6 +8,8 @@
* @brief Implementation code for Gsf Measurement Update
*/
+#include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
+//
#include "TrkGaussianSumFilterUtils/GsfMeasurementUpdator.h"
#include "TrkGaussianSumFilterUtils/GsfConstants.h"
#include "TrkGaussianSumFilterUtils/MultiComponentStateAssembler.h"
@@ -252,8 +254,7 @@ calculateFilterStep_5D(Trk::TrackParameters& TP,
// --- compute filtered covariance matrix
const AmgSymMatrix(5) newCov =
trkCov.similarity(M) + sign * measCov.similarity(K);
- double chiSquared =
- (sign > 0) ? r.transpose() * R * r : r.transpose() * (-R) * r;
+ const double chiSquared = Amg::chi2(R, r, sign);
// create the FQSonSurface
fQ = Trk::FitQualityOnSurface(chiSquared, 5);
TP.updateParameters(newPar, newCov);
@@ -299,8 +300,7 @@ calculateFilterStep_T(Trk::TrackParameters& TP,
// C = M * trkCov * M.T() +/- K * covRio * K.T()
const AmgSymMatrix(5) newCov =
M * trkCov * M.transpose() + sign * K * measCov * K.transpose();
- const double chiSquared =
- (sign > 0) ? r.transpose() * R * r : r.transpose() * (-R) * r;
+ const double chiSquared = Amg::chi2(R, r, sign);
// create the FQSonSurface
fQ = Trk::FitQualityOnSurface(chiSquared, DIM);
// In place update of parameters
@@ -429,10 +429,7 @@ makeChi2_T(Trk::FitQualityOnSurface& updatedFitQoS,
AmgSymMatrix(DIM) R = sign * projection_T<DIM>(trkCov, paramKey);
R += covPar;
// calcualte the chi2 value
- double chiSquared = 0.0;
- if (R.determinant() != 0.0) {
- chiSquared = r.transpose() * R.inverse() * r;
- }
+ const double chiSquared = Amg::chi2(R.inverse(), r);
updatedFitQoS = Trk::FitQualityOnSurface(chiSquared, DIM);
return true;
}
@@ -503,8 +500,7 @@ calculateWeight_T(const Trk::TrackParameters* componentTrackParameters,
return { 0, 0 };
}
// Compute Chi2
- return std::pair<double, double>(
- det, (1. / (double)DIM) * ((r.transpose() * R.inverse() * r)(0, 0)));
+ return {det, (1. / (double)DIM) * Amg::chi2(R.inverse(), r)};
}
std::pair<double, double>
@@ -545,7 +541,7 @@ calculateWeight_2D_3(const Trk::TrackParameters* componentTrackParameters,
return { 0, 0 };
}
// Compute Chi2
- return { det, 0.5 * ((r.transpose() * R.inverse() * r)(0, 0)) };
+ return {det, 0.5 * Amg::chi2(R.inverse(), r)};
}
Trk::MultiComponentState
diff --git a/Tracking/TrkTools/TrkMeasurementUpdator/TrkMeasurementUpdator/KalmanUpdator.h b/Tracking/TrkTools/TrkMeasurementUpdator/TrkMeasurementUpdator/KalmanUpdator.h
index a2e7ca32b488..213db1a4ae11 100755
--- a/Tracking/TrkTools/TrkMeasurementUpdator/TrkMeasurementUpdator/KalmanUpdator.h
+++ b/Tracking/TrkTools/TrkMeasurementUpdator/TrkMeasurementUpdator/KalmanUpdator.h
@@ -15,6 +15,7 @@
#define TRK_KALMANUPDATOR_H
#include "AthenaBaseComps/AthAlgTool.h"
+#include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
#include "EventPrimitives/EventPrimitives.h"
#include "GaudiKernel/MsgStream.h"
#include "GeoPrimitives/GeoPrimitives.h"
@@ -331,7 +332,7 @@ KalmanUpdator::makeChi2Object(const Amg::VectorX& residual,
return FitQualityOnSurface(0.0, (int)covRio.cols());
}
// get chi2 = r.T() * R^-1 * r
- double chiSquared = residual.transpose() * R.inverse() * residual;
+ const double chiSquared = Amg::chi2(R.inverse(), residual);
ATH_MSG_VERBOSE("-U- fitQuality of " << (sign > 0 ? "predicted" : "updated")
<< " state, chi2 :" << chiSquared
<< " / ndof= " << covRio.cols());
diff --git a/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdator.cxx b/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdator.cxx
index e6ba4d44c3b0..caa1ae0799f5 100755
--- a/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdator.cxx
+++ b/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdator.cxx
@@ -250,8 +250,7 @@ std::unique_ptr<Trk::TrackParameters> Trk::KalmanUpdator::combineStates (const T
}
// compute fit quality
- double chiSquared = r.transpose()*R_inv*r;
- fitQoS = new FitQualityOnSurface(chiSquared, 5);
+ fitQoS = new FitQualityOnSurface(Amg::chi2(R_inv, r), 5);
// return cloned version of Track Parameters (MeasuredPerigee, MeasuredAtA...)
auto comb =
@@ -421,8 +420,7 @@ Trk::KalmanUpdator::predictedStateFitQuality (const Trk::TrackParameters& one,
Amg::VectorX r = two.parameters() - one.parameters();
AmgSymMatrix(5) R = (covTrkOne + covTrkTwo).inverse();
// chi2 calculation
- double chiSquared = r.transpose()*R*r;
- return {chiSquared, 5};
+ return {Amg::chi2(R, r), 5};
}
std::vector<double> Trk::KalmanUpdator::initialErrors() const {
diff --git a/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdatorSMatrix.cxx b/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdatorSMatrix.cxx
index dc1c3a691a82..f61861724784 100755
--- a/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdatorSMatrix.cxx
+++ b/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanUpdatorSMatrix.cxx
@@ -14,7 +14,10 @@
// Wolfgang Liebig <http://consult.cern.ch/xwho/people/54608>
///////////////////////////////////////////////////////////////////
+#include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
+
#include "TrkMeasurementUpdator/KalmanUpdatorSMatrix.h"
+
#include "TrkEventPrimitives/ParamDefs.h"
#include "TrkEventPrimitives/LocalParameters.h"
#include "TrkEventPrimitives/FitQualityOnSurface.h"
@@ -1293,7 +1296,7 @@ Trk::FitQualityOnSurface Trk::KalmanUpdatorSMatrix::makeChi2Object(const Amg::Ve
chiSquared = 0.f;
} else {
// get chi2 = r.T() * R^-1 * r
- chiSquared = residual.transpose() * R.inverse() * residual;
+ chiSquared = Amg::chi2(R.inverse(), residual);
ATH_MSG_VERBOSE( "-U- fitQuality of "<< (sign>0?"predicted":"updated")
<<" state, chi2 :" << chiSquared << " / ndof= " << covRio.cols() );
}
diff --git a/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanWeightUpdator.cxx b/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanWeightUpdator.cxx
index d44ff414c32d..3ff9873ba4bd 100755
--- a/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanWeightUpdator.cxx
+++ b/Tracking/TrkTools/TrkMeasurementUpdator/src/KalmanWeightUpdator.cxx
@@ -18,6 +18,7 @@
#include "TrkEventPrimitives/FitQualityOnSurface.h"
#include "TrkParameters/TrackParameters.h"
+#include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
// constructor
Trk::KalmanWeightUpdator::KalmanWeightUpdator(const std::string& t,const std::string& n,const IInterface* p) :
@@ -431,8 +432,7 @@ Trk::KalmanWeightUpdator::predictedStateFitQuality ( const Trk::TrackParamete
AmgSymMatrix(5) R = (covTrkOne + covTrkTwo).inverse();
- double chiSquared = r.transpose()*R*r;
- return {chiSquared, 5};
+ return {Amg::chi2(R, r), 5};
}
std::vector<double> Trk::KalmanWeightUpdator::initialErrors() const {
@@ -675,8 +675,7 @@ Trk::FitQualityOnSurface Trk::KalmanWeightUpdator::makeChi2Object( Amg::VectorX
// attention: similarity(vector) and similarity(matrix) are defined differently in CLHEP:
// similarity(const Amg::VectorX &v): Returns v.T()*s*v
// similarity(const HepSymMatrix &m1): Returns m1*s*m1.T()
- double chiSquared = residual1.transpose()*weight1*residual1;
- chiSquared += residual2.transpose()*weight2*residual2;
+ const double chiSquared = Amg::chi2(weight1, residual1) + Amg::chi2(weight2, residual2);
// number of degree of freedom added
int numberDoF = weight1.cols();
@@ -694,7 +693,7 @@ Trk::FitQualityOnSurface Trk::KalmanWeightUpdator::makeChi2Object( Amg::VectorX
{ // sign: -1 = updated, +1 = predicted parameters.
Amg::MatrixX R = covRio + sign* covTrk.similarity(H);
// get chi2 = r.T() * R^-1 * r
- double chiSquared = residual.transpose()*R.inverse()*residual;
+ const double chiSquared = Amg::chi2(R.inverse(), residual);
ATH_MSG_VERBOSE( "-U- fitQuality of "<< (sign>0?"predicted":"updated")
<<" state, chi2 :" << chiSquared << " / ndof= " << covRio.cols() );
diff --git a/Tracking/TrkVertexFitter/TrkVertexBilloirTools/src/FastVertexFitter.cxx b/Tracking/TrkVertexFitter/TrkVertexBilloirTools/src/FastVertexFitter.cxx
index 9a1d93b16fc8..b7d67c52b939 100755
--- a/Tracking/TrkVertexFitter/TrkVertexBilloirTools/src/FastVertexFitter.cxx
+++ b/Tracking/TrkVertexFitter/TrkVertexBilloirTools/src/FastVertexFitter.cxx
@@ -15,6 +15,7 @@
#include "TrkTrack/LinkToTrack.h"
#include "TrkParticleBase/LinkToTrackParticleBase.h"
#include "TrkParticleBase/TrackParticleBase.h"
+#include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
#include "EventPrimitives/EventPrimitives.h"
#include "TrkLinks/LinkToXAODTrackParticle.h"
//xAOD includes
@@ -260,8 +261,7 @@ namespace Trk
deltaTrk[0] = delta_V[0] - ( firstStartingPoint.position() [0] - linPoint [0] );
deltaTrk[1] = delta_V[1] - ( firstStartingPoint.position() [1] - linPoint [1] );
deltaTrk[2] = delta_V[2] - ( firstStartingPoint.position() [2] - linPoint [2] );
- double chi2FromConstraint ( ( deltaTrk.transpose() * firstStartingPoint.covariancePosition().inverse().eval() * deltaTrk ) [0] );
- chi2New += chi2FromConstraint;
+ chi2New += Amg::chi2(firstStartingPoint.covariancePosition().inverse().eval(), deltaTrk);
}
/* assign new linearization point (= new vertex position in global frame) */
--
GitLab