From f09e7d910664d617403ca98ba5d0a7fc73e3e8ea Mon Sep 17 00:00:00 2001 From: russofel Date: Wed, 29 Nov 2023 16:36:58 +0100 Subject: [PATCH 01/12] introduce dummy template param --- .../Vertexing/AdaptiveMultiVertexFitter.ipp | 6 +++--- .../Vertexing/KalmanVertexTrackUpdater.ipp | 4 ++-- .../Acts/Vertexing/KalmanVertexUpdater.hpp | 20 +++++++++++-------- .../Acts/Vertexing/KalmanVertexUpdater.ipp | 20 +++++++++---------- .../KalmanVertexTrackUpdaterTests.cpp | 2 +- .../Vertexing/KalmanVertexUpdaterTests.cpp | 4 ++-- 6 files changed, 30 insertions(+), 26 deletions(-) diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp index 2c2a3969180..cf71f0c211d 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp @@ -308,8 +308,8 @@ Acts::Result Acts:: trkAtVtx.isLinearized = true; } // Update the vertex with the new track - KalmanVertexUpdater::updateVertexWithTrack(*vtx, - trkAtVtx); + KalmanVertexUpdater::updateVertexWithTrack(*vtx, + trkAtVtx); } else { ACTS_VERBOSE("Track weight too low. Skip track."); } @@ -369,7 +369,7 @@ void Acts::AdaptiveMultiVertexFitter< for (const auto trk : state.vtxInfoMap[vtx].trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); if (trkAtVtx.trackWeight > m_cfg.minWeight) { - KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); + KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); } } } diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index 8ff1cb323c6..4c80ea69800 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -11,7 +11,7 @@ #include "Acts/Utilities/detail/periodic.hpp" #include "Acts/Vertexing/VertexingError.hpp" -template +template void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, const Vertex& vtx) { // Check if linearized state exists @@ -44,7 +44,7 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, const ActsVector<5> constTerm = linTrack.constantTerm.head<5>(); // Cache object filled with zeros - KalmanVertexUpdater::Cache cache; + KalmanVertexUpdater::Cache cache; // Calculate the update of the vertex position when the track is removed. This // might be unintuitive, but it is needed to compute a symmetric chi2. diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp index 29d4fe7920d..fd135e4d1e1 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp @@ -35,7 +35,9 @@ namespace KalmanVertexUpdater { /// Section 5.3.5 /// Cache object, the comments indicate the names of the variables in Ref. (1) +template struct Cache { + using VertexPosition = ActsVector; // \tilde{x_k} Vector3 newVertexPos = Vector3::Zero(); // C_k @@ -55,7 +57,7 @@ struct Cache { /// /// @param vtx Vertex to be updated /// @param trk Track to be used for updating the vertex -template +template void updateVertexWithTrack(Vertex& vtx, TrackAtVertex& trk); @@ -69,10 +71,11 @@ void updateVertexWithTrack(Vertex& vtx, /// @note Tracks are removed during the smoothing procedure to compute /// the chi2 of the track wrt the updated vertex position /// @param[out] cache A cache to store the results of this function -template +template void calculateUpdate(const Acts::Vertex& vtx, const Acts::LinearizedTrack& linTrack, - const double trackWeight, const int sign, Cache& cache); + const double trackWeight, const int sign, + Cache& cache); namespace detail { @@ -83,9 +86,9 @@ namespace detail { /// @param cache Cache containing updated vertex position /// /// @return Chi2 -template +template double vertexPositionChi2Update(const Vertex& oldVtx, - const Cache& cache); + const Cache& cache); /// @brief Calculates chi2 of refitted track parameters /// w.r.t. updated vertex @@ -95,8 +98,9 @@ double vertexPositionChi2Update(const Vertex& oldVtx, /// this function /// /// @return Chi2 -template -double trackParametersChi2(const LinearizedTrack& linTrack, const Cache& cache); +template +double trackParametersChi2(const LinearizedTrack& linTrack, + const Cache& cache); /// @brief Adds or removes (depending on `sign`) tracks from vertex /// and updates the vertex @@ -106,7 +110,7 @@ double trackParametersChi2(const LinearizedTrack& linTrack, const Cache& cache); /// @param sign +1 (add track) or -1 (remove track) /// @note Tracks are removed during the smoothing procedure to compute /// the chi2 of the track wrt the updated vertex position -template +template void update(Vertex& vtx, TrackAtVertex& trk, int sign); } // Namespace detail diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp index 8138eb5e7a8..2957fb46bc6 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp @@ -10,19 +10,19 @@ #include -template +template void Acts::KalmanVertexUpdater::updateVertexWithTrack( Vertex& vtx, TrackAtVertex& trk) { - detail::update(vtx, trk, 1); + detail::update(vtx, trk, 1); } -template +template void Acts::KalmanVertexUpdater::detail::update( Vertex& vtx, TrackAtVertex& trk, int sign) { double trackWeight = trk.trackWeight; // Set up cache where entire content is set to 0 - Cache cache; + Cache cache; // Calculate update and save result in cache calculateUpdate(vtx, trk.linearizedState, trackWeight, sign, cache); @@ -65,11 +65,11 @@ void Acts::KalmanVertexUpdater::detail::update( } } -template +template void Acts::KalmanVertexUpdater::calculateUpdate( const Acts::Vertex& vtx, const Acts::LinearizedTrack& linTrack, const double trackWeight, - const int sign, Cache& cache) { + const int sign, Cache& cache) { // Retrieve variables from the track linearization. The comments indicate the // corresponding symbol used in Ref. (1). // A_k @@ -118,18 +118,18 @@ void Acts::KalmanVertexUpdater::calculateUpdate( (trkParams - constTerm)); } -template +template double Acts::KalmanVertexUpdater::detail::vertexPositionChi2Update( - const Vertex& oldVtx, const Cache& cache) { + const Vertex& oldVtx, const Cache& cache) { Vector3 posDiff = cache.newVertexPos - oldVtx.position(); // Calculate and return corresponding chi2 return posDiff.transpose() * (cache.oldVertexWeight * posDiff); } -template +template double Acts::KalmanVertexUpdater::detail::trackParametersChi2( - const LinearizedTrack& linTrack, const Cache& cache) { + const LinearizedTrack& linTrack, const Cache& cache) { // A_k const ActsMatrix<5, 3> posJac = linTrack.positionJacobian.block<5, 3>(0, 0); // B_k diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp index 871899777ac..83d5fcb8935 100644 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp @@ -178,7 +178,7 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) { Vertex vtx(vtxPos); // Update trkAtVertex with assumption of originating from vtx - KalmanVertexTrackUpdater::update(trkAtVtx, vtx); + KalmanVertexTrackUpdater::update(trkAtVtx, vtx); // The old distance double oldDistance = diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp index 0360756c319..6d2ddea556d 100644 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp @@ -172,8 +172,8 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) { vtx.setFullCovariance(SquareMatrix4::Identity() * 0.01); // Update trkAtVertex with assumption of originating from vtx - KalmanVertexUpdater::updateVertexWithTrack(vtx, - trkAtVtx); + KalmanVertexUpdater::updateVertexWithTrack( + vtx, trkAtVtx); if (debug) { std::cout << "Old vertex position: " << vtxPos << std::endl; From b176406019b9e4adf9a64b2c4c7eb76205d95c73 Mon Sep 17 00:00:00 2001 From: russofel Date: Wed, 29 Nov 2023 17:26:32 +0100 Subject: [PATCH 02/12] implement templatng for KalmanVertexUpdater --- .../Vertexing/KalmanVertexTrackUpdater.ipp | 11 +-- .../Acts/Vertexing/KalmanVertexUpdater.hpp | 11 +-- .../Acts/Vertexing/KalmanVertexUpdater.ipp | 79 +++++++++++++------ 3 files changed, 66 insertions(+), 35 deletions(-) diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index 4c80ea69800..3707b0b14de 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -36,8 +36,8 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, const ActsMatrix<5, 3> momJac = linTrack.momentumJacobian.block<5, 3>(0, 0); // p_k const ActsVector<5> trkParams = linTrack.parametersAtPCA.head<5>(); - // TODO we could use `linTrack.weightAtPCA` but only if we would use time - // G_k + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time G_k const ActsSquareMatrix<5> trkParamWeight = linTrack.covarianceAtPCA.block<5, 5>(0, 0).inverse(); // c_k @@ -71,15 +71,16 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, -vtxCov * posJac.transpose() * trkParamWeight * momJac * cache.wMat; // Difference in positions. cache.newVertexPos corresponds to \tilde{x_k^{n*}} in Ref. (1). - Vector3 posDiff = vtxPos - cache.newVertexPos; + Vector3 posDiff = vtxPos - cache.newVertexPos.template head<3>(); // r_k^n ActsVector<5> paramDiff = trkParams - (constTerm + posJac * vtxPos + momJac * newTrkMomentum); // New chi2 to be set later - double chi2 = posDiff.dot(cache.newVertexWeight * posDiff) + - paramDiff.dot(trkParamWeight * paramDiff); + double chi2 = + posDiff.dot(cache.newVertexWeight.template block<3, 3>(0, 0) * posDiff) + + paramDiff.dot(trkParamWeight * paramDiff); Acts::BoundMatrix fullPerTrackCov = detail::createFullTrackCovariance( cache.wMat, crossCovVP, vtxCov, newTrkParams); diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp index fd135e4d1e1..d084c6e9888 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp @@ -37,15 +37,16 @@ namespace KalmanVertexUpdater { /// Cache object, the comments indicate the names of the variables in Ref. (1) template struct Cache { - using VertexPosition = ActsVector; + using VertexVector = ActsVector; + using VertexMatrix = ActsSquareMatrix; // \tilde{x_k} - Vector3 newVertexPos = Vector3::Zero(); + VertexVector newVertexPos = VertexVector::Zero(); // C_k - SquareMatrix3 newVertexCov = SquareMatrix3::Zero(); + VertexMatrix newVertexCov = VertexMatrix::Zero(); // C_k^-1 - SquareMatrix3 newVertexWeight = SquareMatrix3::Zero(); + VertexMatrix newVertexWeight = VertexMatrix::Zero(); // C_{k-1}^-1 - SquareMatrix3 oldVertexWeight = SquareMatrix3::Zero(); + VertexMatrix oldVertexWeight = VertexMatrix::Zero(); // W_k SquareMatrix3 wMat = SquareMatrix3::Zero(); }; diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp index 2957fb46bc6..a62e668de8f 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp @@ -19,6 +19,12 @@ void Acts::KalmanVertexUpdater::updateVertexWithTrack( template void Acts::KalmanVertexUpdater::detail::update( Vertex& vtx, TrackAtVertex& trk, int sign) { + if constexpr (nDimVertex != 3 && nDimVertex != 4) { + throw std::invalid_argument( + "The vertex dimension must either be 3 (when fitting the spatial " + "coordinates) or 4 (when fitting the spatial coordinates + time)."); + } + double trackWeight = trk.trackWeight; // Set up cache where entire content is set to 0 @@ -47,8 +53,13 @@ void Acts::KalmanVertexUpdater::detail::update( ndf += sign * trackWeight * 2.; // Updating the vertex - vtx.setPosition(cache.newVertexPos); - vtx.setCovariance(cache.newVertexCov); + if constexpr (nDimVertex == 3) { + vtx.setPosition(cache.newVertexPos); + vtx.setCovariance(cache.newVertexCov); + } else if constexpr (nDimVertex == 4) { + vtx.setFullPosition(cache.newVertexPos); + vtx.setFullCovariance(cache.newVertexCov); + } vtx.setFitQuality(chi2, ndf); if (sign == 1) { @@ -70,31 +81,38 @@ void Acts::KalmanVertexUpdater::calculateUpdate( const Acts::Vertex& vtx, const Acts::LinearizedTrack& linTrack, const double trackWeight, const int sign, Cache& cache) { + constexpr unsigned int nParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // Retrieve variables from the track linearization. The comments indicate the // corresponding symbol used in Ref. (1). // A_k - const ActsMatrix<5, 3> posJac = linTrack.positionJacobian.block<5, 3>(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix<5, 3> momJac = linTrack.momentumJacobian.block<5, 3>(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ActsVector<5> trkParams = linTrack.parametersAtPCA.head<5>(); + const ParameterVector trkParams = linTrack.parametersAtPCA.head(); // c_k - const ActsVector<5> constTerm = linTrack.constantTerm.head<5>(); - // TODO we could use `linTrack.weightAtPCA` but only if we would use time - // G_k - // Note that, when removing a track, G_k -> - G_k, see Ref. (1). + const ParameterVector constTerm = linTrack.constantTerm.head(); + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time G_k Note that, when removing a track, G_k -> - G_k, see Ref. (1). // Further note that, as we use the weighted formalism, the track weight // matrix (i.e., the inverse track covariance matrix) should be multiplied // with the track weight from the AMVF formalism. Here, we choose to // consider these two multiplicative factors directly in the updates of // newVertexWeight and newVertexPos. - const ActsSquareMatrix<5> trkParamWeight = - linTrack.covarianceAtPCA.block<5, 5>(0, 0).inverse(); + const ParameterMatrix trkParamWeight = + linTrack.covarianceAtPCA.block(0, 0).inverse(); // Retrieve current position of the vertex and its current weight matrix - const Vector3& oldVtxPos = vtx.position(); + const ActsVector oldVtxPos = + vtx.fullPosition().template head(); // C_{k-1}^-1 - cache.oldVertexWeight = (vtx.covariance()).inverse(); + cache.oldVertexWeight = + (vtx.fullCovariance().template block(0, 0)) + .inverse(); // W_k cache.wMat = (momJac.transpose() * (trkParamWeight * momJac)).inverse(); @@ -121,7 +139,8 @@ void Acts::KalmanVertexUpdater::calculateUpdate( template double Acts::KalmanVertexUpdater::detail::vertexPositionChi2Update( const Vertex& oldVtx, const Cache& cache) { - Vector3 posDiff = cache.newVertexPos - oldVtx.position(); + ActsVector posDiff = + cache.newVertexPos - oldVtx.fullPosition().template head(); // Calculate and return corresponding chi2 return posDiff.transpose() * (cache.oldVertexWeight * posDiff); @@ -130,21 +149,31 @@ double Acts::KalmanVertexUpdater::detail::vertexPositionChi2Update( template double Acts::KalmanVertexUpdater::detail::trackParametersChi2( const LinearizedTrack& linTrack, const Cache& cache) { + constexpr unsigned int nParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // A_k - const ActsMatrix<5, 3> posJac = linTrack.positionJacobian.block<5, 3>(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix<5, 3> momJac = linTrack.momentumJacobian.block<5, 3>(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ActsVector<5> trkParams = linTrack.parametersAtPCA.head<5>(); + const ParameterVector trkParams = linTrack.parametersAtPCA.head(); // c_k - const ActsVector<5> constTerm = linTrack.constantTerm.head<5>(); - // TODO we could use `linTrack.weightAtPCA` but only if we would use time - // G_k - const ActsSquareMatrix<5> trkParamWeight = - linTrack.covarianceAtPCA.block<5, 5>(0, 0).inverse(); + const ParameterVector constTerm = linTrack.constantTerm.head(); + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time G_k Note that, when removing a track, G_k -> - G_k, see Ref. (1). + // Further note that, as we use the weighted formalism, the track weight + // matrix (i.e., the inverse track covariance matrix) should be multiplied + // with the track weight from the AMVF formalism. Here, we choose to + // consider these two multiplicative factors directly in the updates of + // newVertexWeight and newVertexPos. + const ParameterMatrix trkParamWeight = + linTrack.covarianceAtPCA.block(0, 0).inverse(); // A_k * \tilde{x_k} - const ActsVector<5> posJacVtxPos = posJac * cache.newVertexPos; + const ParameterVector posJacVtxPos = posJac * cache.newVertexPos; // \tilde{q_k} Vector3 newTrkMom = cache.wMat * momJac.transpose() * trkParamWeight * @@ -161,11 +190,11 @@ double Acts::KalmanVertexUpdater::detail::trackParametersChi2( */ // \tilde{p_k} - ActsVector<5> linearizedTrackParameters = + ParameterVector linearizedTrackParameters = constTerm + posJacVtxPos + momJac * newTrkMom; // r_k - ActsVector<5> paramDiff = trkParams - linearizedTrackParameters; + ParameterVector paramDiff = trkParams - linearizedTrackParameters; // Return chi2 return paramDiff.transpose() * (trkParamWeight * paramDiff); From f7db3ac37941ede204cf47a89d77edadf0e8f0d3 Mon Sep 17 00:00:00 2001 From: russofel Date: Wed, 29 Nov 2023 19:38:45 +0100 Subject: [PATCH 03/12] implement templating for KalmanVertexTrackUpdater --- .../Vertexing/KalmanVertexTrackUpdater.hpp | 11 +-- .../Vertexing/KalmanVertexTrackUpdater.ipp | 88 ++++++++++++------- .../Acts/Vertexing/KalmanVertexUpdater.ipp | 12 ++- 3 files changed, 69 insertions(+), 42 deletions(-) diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp index 72b5f9b668d..28707603ce3 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp @@ -24,7 +24,7 @@ namespace KalmanVertexTrackUpdater { /// /// @param track Track to update /// @param vtx Vertex `track` belongs to -template +template void update(TrackAtVertex& track, const Vertex& vtx); @@ -37,10 +37,11 @@ namespace detail { /// @param newTrkCov New track covariance matrixs /// @param vtxCov Vertex covariance matrix /// @param newTrkParams New track parameter -inline BoundMatrix createFullTrackCovariance(const SquareMatrix3& wMat, - const SquareMatrix3& crossCovVP, - const SquareMatrix3& vtxCov, - const BoundVector& newTrkParams); +template +inline BoundMatrix createFullTrackCovariance( + const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, + const ActsSquareMatrix& vtxCov, + const BoundVector& newTrkParams); } // Namespace detail diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index 3707b0b14de..57826a46e8a 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -14,6 +14,11 @@ template void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, const Vertex& vtx) { + using VertexVector = ActsVector; + using VertexMatrix = ActsSquareMatrix; + constexpr unsigned int nParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // Check if linearized state exists if (!track.isLinearized) { throw std::invalid_argument("TrackAtVertex object must be linearized."); @@ -21,27 +26,30 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, // Extract vertex position and covariance // \tilde{x_n} - const Vector3& vtxPos = vtx.position(); + const VertexVector vtxPos = vtx.fullPosition().template head(); // C_n - const SquareMatrix3& vtxCov = vtx.covariance(); + const VertexMatrix vtxCov = + vtx.fullCovariance().template block(0, 0); // Get the linearized track const LinearizedTrack& linTrack = track.linearizedState; - // Retrieve variables from the track linearization. The comments indicate the - // corresponding symbol used in the Ref. (1). + // corresponding symbol used in Ref. (1). // A_k - const ActsMatrix<5, 3> posJac = linTrack.positionJacobian.block<5, 3>(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix<5, 3> momJac = linTrack.momentumJacobian.block<5, 3>(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ActsVector<5> trkParams = linTrack.parametersAtPCA.head<5>(); - // TODO we could use `linTrack.weightAtPCA` but only if we would always fit - // time G_k - const ActsSquareMatrix<5> trkParamWeight = - linTrack.covarianceAtPCA.block<5, 5>(0, 0).inverse(); + const ParameterVector trkParams = linTrack.parametersAtPCA.head(); // c_k - const ActsVector<5> constTerm = linTrack.constantTerm.head<5>(); + const ParameterVector constTerm = linTrack.constantTerm.head(); + // TODO we could use `linTrack.weightAtPCA` but only if we would always fit + // time. + // G_k + const ParameterMatrix trkParamWeight = + linTrack.covarianceAtPCA.block(0, 0).inverse(); // Cache object filled with zeros KalmanVertexUpdater::Cache cache; @@ -67,14 +75,14 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, newTrkParams(BoundIndices::eBoundQOverP) = newTrkMomentum(2); // qOverP // E_k^n - const SquareMatrix3 crossCovVP = + const ActsMatrix crossCovVP = -vtxCov * posJac.transpose() * trkParamWeight * momJac * cache.wMat; // Difference in positions. cache.newVertexPos corresponds to \tilde{x_k^{n*}} in Ref. (1). - Vector3 posDiff = vtxPos - cache.newVertexPos.template head<3>(); + VertexVector posDiff = vtxPos - cache.newVertexPos.template head<3>(); // r_k^n - ActsVector<5> paramDiff = + ParameterVector paramDiff = trkParams - (constTerm + posJac * vtxPos + momJac * newTrkMomentum); // New chi2 to be set later @@ -82,8 +90,9 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, posDiff.dot(cache.newVertexWeight.template block<3, 3>(0, 0) * posDiff) + paramDiff.dot(trkParamWeight * paramDiff); - Acts::BoundMatrix fullPerTrackCov = detail::createFullTrackCovariance( - cache.wMat, crossCovVP, vtxCov, newTrkParams); + Acts::BoundMatrix fullPerTrackCov = + detail::createFullTrackCovariance(cache.wMat, crossCovVP, + vtxCov, newTrkParams); // Create new refitted parameters std::shared_ptr perigeeSurface = @@ -101,25 +110,41 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, return; } +template inline Acts::BoundMatrix Acts::KalmanVertexTrackUpdater::detail::createFullTrackCovariance( - const SquareMatrix3& wMat, const SquareMatrix3& crossCovVP, - const SquareMatrix3& vtxCov, const BoundVector& newTrkParams) { + const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, + const ActsSquareMatrix& vtxCov, + const BoundVector& newTrkParams) { // D_k^n ActsSquareMatrix<3> momCov = wMat + crossCovVP.transpose() * vtxCov.inverse() * crossCovVP; - // Full (x,y,z,phi, theta, q/p) covariance matrix - // To be made 7d again after switching to (x,y,z,phi, theta, q/p, t) - ActsSquareMatrix<6> fullTrkCov(ActsSquareMatrix<6>::Zero()); - - fullTrkCov.block<3, 3>(0, 0) = vtxCov; - fullTrkCov.block<3, 3>(0, 3) = crossCovVP; - fullTrkCov.block<3, 3>(3, 0) = crossCovVP.transpose(); - fullTrkCov.block<3, 3>(3, 3) = momCov; + // Full x, y, z, phi, theta, q/p, and, optionally, t covariance matrix. Note + // that we call this set of parameters here free even though that is not quite + // correct (free parameters corresponds to x, y, z, t, px, py, pz) + constexpr unsigned int nFreeParams = nDimVertex + 3; + ActsSquareMatrix fullTrkCov( + ActsSquareMatrix::Zero()); + + fullTrkCov.template block<3, 3>(0, 0) = vtxCov.template block<3, 3>(0, 0); + fullTrkCov.template block<3, 3>(0, 3) = crossCovVP.template block<3, 3>(0, 0); + fullTrkCov.template block<3, 3>(3, 0) = + crossCovVP.template block<3, 3>(0, 0).transpose(); + fullTrkCov.template block<3, 3>(3, 3) = momCov; + + // Fill time cross-covariances + if constexpr (nFreeParams == 7) { + fullTrkCov.block<3, 1>(0, 6) = vtxCov.block<3, 1>(0, 3); + fullTrkCov.block<3, 1>(3, 6) = crossCovVP.block<3, 1>(0, 3); + fullTrkCov.block<1, 3>(6, 0) = vtxCov.block<1, 3>(3, 0); + fullTrkCov.block<1, 3>(6, 3) = crossCovVP.block<1, 3>(3, 0); + } // Combined track jacobian - ActsMatrix<5, 6> trkJac(ActsMatrix<5, 6>::Zero()); + constexpr unsigned int nBoundParams = nDimVertex + 2; + ActsMatrix trkJac( + ActsMatrix::Zero()); // First row trkJac(0, 0) = -std::sin(newTrkParams[2]); @@ -131,11 +156,14 @@ Acts::KalmanVertexTrackUpdater::detail::createFullTrackCovariance( trkJac(1, 0) = -trkJac(0, 1) / tanTheta; trkJac(1, 1) = trkJac(0, 0) / tanTheta; - trkJac.block<4, 4>(1, 2) = ActsMatrix<4, 4>::Identity(); + // Dimension of the part of the Jacobian that is an identity matrix + constexpr unsigned int nDimIdentity = nFreeParams - 2; + trkJac.template block(1, 2) = + ActsMatrix::Identity(); // Full perigee track covariance BoundMatrix fullPerTrackCov(BoundMatrix::Identity()); - fullPerTrackCov.block<5, 5>(0, 0) = + fullPerTrackCov.block(0, 0) = (trkJac * (fullTrkCov * trkJac.transpose())); return fullPerTrackCov; diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp index a62e668de8f..b478c0fc5af 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp @@ -97,7 +97,9 @@ void Acts::KalmanVertexUpdater::calculateUpdate( // c_k const ParameterVector constTerm = linTrack.constantTerm.head(); // TODO we could use `linTrack.weightAtPCA` but only if we would always fit - // time G_k Note that, when removing a track, G_k -> - G_k, see Ref. (1). + // time. + // G_k + // Note that, when removing a track, G_k -> - G_k, see Ref. (1). // Further note that, as we use the weighted formalism, the track weight // matrix (i.e., the inverse track covariance matrix) should be multiplied // with the track weight from the AMVF formalism. Here, we choose to @@ -163,12 +165,8 @@ double Acts::KalmanVertexUpdater::detail::trackParametersChi2( // c_k const ParameterVector constTerm = linTrack.constantTerm.head(); // TODO we could use `linTrack.weightAtPCA` but only if we would always fit - // time G_k Note that, when removing a track, G_k -> - G_k, see Ref. (1). - // Further note that, as we use the weighted formalism, the track weight - // matrix (i.e., the inverse track covariance matrix) should be multiplied - // with the track weight from the AMVF formalism. Here, we choose to - // consider these two multiplicative factors directly in the updates of - // newVertexWeight and newVertexPos. + // time. + // G_k const ParameterMatrix trkParamWeight = linTrack.covarianceAtPCA.block(0, 0).inverse(); From 35444da85e73922d85b754e5d10f1936ab5a799c Mon Sep 17 00:00:00 2001 From: russofel Date: Wed, 29 Nov 2023 20:18:09 +0100 Subject: [PATCH 04/12] bugs --- .../Vertexing/KalmanVertexTrackUpdater.ipp | 19 ++++++++++++------- .../Acts/Vertexing/KalmanVertexUpdater.ipp | 6 +++--- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index 57826a46e8a..ed9b950bf71 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -79,7 +79,8 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, -vtxCov * posJac.transpose() * trkParamWeight * momJac * cache.wMat; // Difference in positions. cache.newVertexPos corresponds to \tilde{x_k^{n*}} in Ref. (1). - VertexVector posDiff = vtxPos - cache.newVertexPos.template head<3>(); + VertexVector posDiff = + vtxPos - cache.newVertexPos.template head(); // r_k^n ParameterVector paramDiff = @@ -87,7 +88,9 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, // New chi2 to be set later double chi2 = - posDiff.dot(cache.newVertexWeight.template block<3, 3>(0, 0) * posDiff) + + posDiff.dot( + cache.newVertexWeight.template block(0, 0) * + posDiff) + paramDiff.dot(trkParamWeight * paramDiff); Acts::BoundMatrix fullPerTrackCov = @@ -96,7 +99,7 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, // Create new refitted parameters std::shared_ptr perigeeSurface = - Surface::makeShared(vtxPos); + Surface::makeShared(vtxPos.template head<3>()); BoundTrackParameters refittedPerigee = BoundTrackParameters( perigeeSurface, newTrkParams, std::move(fullPerTrackCov), @@ -135,10 +138,12 @@ Acts::KalmanVertexTrackUpdater::detail::createFullTrackCovariance( // Fill time cross-covariances if constexpr (nFreeParams == 7) { - fullTrkCov.block<3, 1>(0, 6) = vtxCov.block<3, 1>(0, 3); - fullTrkCov.block<3, 1>(3, 6) = crossCovVP.block<3, 1>(0, 3); - fullTrkCov.block<1, 3>(6, 0) = vtxCov.block<1, 3>(3, 0); - fullTrkCov.block<1, 3>(6, 3) = crossCovVP.block<1, 3>(3, 0); + fullTrkCov.template block<3, 1>(0, 6) = vtxCov.template block<3, 1>(0, 3); + fullTrkCov.template block<3, 1>(3, 6) = + crossCovVP.template block<3, 1>(0, 3); + fullTrkCov.template block<1, 3>(6, 0) = vtxCov.template block<1, 3>(3, 0); + fullTrkCov.template block<1, 3>(6, 3) = + crossCovVP.template block<1, 3>(3, 0); } // Combined track jacobian diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp index b478c0fc5af..1cddefab3b9 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp @@ -120,9 +120,9 @@ void Acts::KalmanVertexUpdater::calculateUpdate( cache.wMat = (momJac.transpose() * (trkParamWeight * momJac)).inverse(); // G_k^B = G_k - G_k*B_k*W_k*B_k^(T)*G_k - ActsSquareMatrix<5> gBMat = - trkParamWeight - trkParamWeight * momJac * cache.wMat * - momJac.transpose() * trkParamWeight; + ParameterMatrix gBMat = trkParamWeight - trkParamWeight * momJac * + cache.wMat * momJac.transpose() * + trkParamWeight; // C_k^-1 cache.newVertexWeight = cache.oldVertexWeight + sign * trackWeight * From 90f281d45e36fb6673afc12aa8b1bb3700902af0 Mon Sep 17 00:00:00 2001 From: russofel Date: Wed, 29 Nov 2023 20:18:27 +0100 Subject: [PATCH 05/12] fit time if useTime==true --- .../Vertexing/AdaptiveMultiVertexFitter.ipp | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp index cf71f0c211d..b464800d617 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp @@ -307,9 +307,16 @@ Acts::Result Acts:: trkAtVtx.linearizedState = *result; trkAtVtx.isLinearized = true; } - // Update the vertex with the new track - KalmanVertexUpdater::updateVertexWithTrack(*vtx, - trkAtVtx); + // Update the vertex with the new track. The second template argument + // corresponds to the number of fitted vertex dimensions (i.e., 3 if we + // only fit spatial coordinates and 4 if we also fit time) + if (m_cfg.useTime) { + KalmanVertexUpdater::updateVertexWithTrack( + *vtx, trkAtVtx); + } else { + KalmanVertexUpdater::updateVertexWithTrack( + *vtx, trkAtVtx); + } } else { ACTS_VERBOSE("Track weight too low. Skip track."); } @@ -369,7 +376,11 @@ void Acts::AdaptiveMultiVertexFitter< for (const auto trk : state.vtxInfoMap[vtx].trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); if (trkAtVtx.trackWeight > m_cfg.minWeight) { - KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); + if (m_cfg.useTime) { + KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); + } else { + KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); + } } } } From 37c6ffdb603ae26cf11b4b4ce6bf15bedc407f1d Mon Sep 17 00:00:00 2001 From: russofel Date: Thu, 30 Nov 2023 13:44:57 +0100 Subject: [PATCH 06/12] remove bug and small improvements --- .../Vertexing/KalmanVertexTrackUpdater.hpp | 21 ++++-- .../Vertexing/KalmanVertexTrackUpdater.ipp | 64 ++++++++++--------- 2 files changed, 47 insertions(+), 38 deletions(-) diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp index 28707603ce3..3868183ac22 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp @@ -18,7 +18,14 @@ namespace Acts { namespace KalmanVertexTrackUpdater { /// KalmanVertexTrackUpdater -/// +/// Based on +/// Ref. (1): +/// R. Frühwirth et al. +/// Vertex reconstruction and track bundling at the lep collider using robust +/// algorithms +/// Computer Physics Comm.: 96 (1996) 189 +/// Chapter 2.1 + /// @brief Refits a single track with the knowledge of /// the vertex it has originated from /// @@ -30,15 +37,15 @@ void update(TrackAtVertex& track, namespace detail { -/// @brief reates a new covariance matrix for the -/// refitted track parameters +/// @brief Calculates a covariance matrix for the refitted track parameters /// -/// @param sMat Track ovariance in momentum space -/// @param newTrkCov New track covariance matrixs +/// @param wMat W_k matrix from Ref. (1) +/// @param crossCovVP Cross-covariance matrix between vertex position and track +/// momentum /// @param vtxCov Vertex covariance matrix -/// @param newTrkParams New track parameter +/// @param newTrkParams Refitted track parameters template -inline BoundMatrix createFullTrackCovariance( +inline BoundMatrix calculateTrackCovariance( const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, const ActsSquareMatrix& vtxCov, const BoundVector& newTrkParams); diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index ed9b950bf71..305fc373f49 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -93,17 +93,16 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, posDiff) + paramDiff.dot(trkParamWeight * paramDiff); - Acts::BoundMatrix fullPerTrackCov = - detail::createFullTrackCovariance(cache.wMat, crossCovVP, - vtxCov, newTrkParams); + Acts::BoundMatrix newTrackCov = detail::calculateTrackCovariance( + cache.wMat, crossCovVP, vtxCov, newTrkParams); // Create new refitted parameters std::shared_ptr perigeeSurface = Surface::makeShared(vtxPos.template head<3>()); - BoundTrackParameters refittedPerigee = BoundTrackParameters( - perigeeSurface, newTrkParams, std::move(fullPerTrackCov), - track.fittedParams.particleHypothesis()); + BoundTrackParameters refittedPerigee = + BoundTrackParameters(perigeeSurface, newTrkParams, std::move(newTrackCov), + track.fittedParams.particleHypothesis()); // Set new properties track.fittedParams = refittedPerigee; @@ -115,7 +114,7 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, template inline Acts::BoundMatrix -Acts::KalmanVertexTrackUpdater::detail::createFullTrackCovariance( +Acts::KalmanVertexTrackUpdater::detail::calculateTrackCovariance( const SquareMatrix3& wMat, const ActsMatrix& crossCovVP, const ActsSquareMatrix& vtxCov, const BoundVector& newTrkParams) { @@ -124,52 +123,55 @@ Acts::KalmanVertexTrackUpdater::detail::createFullTrackCovariance( wMat + crossCovVP.transpose() * vtxCov.inverse() * crossCovVP; // Full x, y, z, phi, theta, q/p, and, optionally, t covariance matrix. Note - // that we call this set of parameters here free even though that is not quite - // correct (free parameters corresponds to x, y, z, t, px, py, pz) + // that we call this set of parameters "free" in the following even though + // that is not quite correct (free parameters correspond to + // x, y, z, t, px, py, pz) constexpr unsigned int nFreeParams = nDimVertex + 3; - ActsSquareMatrix fullTrkCov( + ActsSquareMatrix freeTrkCov( ActsSquareMatrix::Zero()); - fullTrkCov.template block<3, 3>(0, 0) = vtxCov.template block<3, 3>(0, 0); - fullTrkCov.template block<3, 3>(0, 3) = crossCovVP.template block<3, 3>(0, 0); - fullTrkCov.template block<3, 3>(3, 0) = + freeTrkCov.template block<3, 3>(0, 0) = vtxCov.template block<3, 3>(0, 0); + freeTrkCov.template block<3, 3>(0, 3) = crossCovVP.template block<3, 3>(0, 0); + freeTrkCov.template block<3, 3>(3, 0) = crossCovVP.template block<3, 3>(0, 0).transpose(); - fullTrkCov.template block<3, 3>(3, 3) = momCov; + freeTrkCov.template block<3, 3>(3, 3) = momCov; - // Fill time cross-covariances + // Fill time (cross-)covariances if constexpr (nFreeParams == 7) { - fullTrkCov.template block<3, 1>(0, 6) = vtxCov.template block<3, 1>(0, 3); - fullTrkCov.template block<3, 1>(3, 6) = - crossCovVP.template block<3, 1>(0, 3); - fullTrkCov.template block<1, 3>(6, 0) = vtxCov.template block<1, 3>(3, 0); - fullTrkCov.template block<1, 3>(6, 3) = + freeTrkCov.template block<3, 1>(0, 6) = vtxCov.template block<3, 1>(0, 3); + freeTrkCov.template block<3, 1>(3, 6) = + crossCovVP.template block<1, 3>(3, 0).transpose(); + freeTrkCov.template block<1, 3>(6, 0) = vtxCov.template block<1, 3>(3, 0); + freeTrkCov.template block<1, 3>(6, 3) = crossCovVP.template block<1, 3>(3, 0); + freeTrkCov(6, 6) = vtxCov(3, 3); } - // Combined track jacobian + // Jacobian relating "free" and bound track parameters constexpr unsigned int nBoundParams = nDimVertex + 2; - ActsMatrix trkJac( + ActsMatrix freeToBoundJac( ActsMatrix::Zero()); + // TODO: Jacobian is not quite correct // First row - trkJac(0, 0) = -std::sin(newTrkParams[2]); - trkJac(0, 1) = std::cos(newTrkParams[2]); + freeToBoundJac(0, 0) = -std::sin(newTrkParams[2]); + freeToBoundJac(0, 1) = std::cos(newTrkParams[2]); double tanTheta = std::tan(newTrkParams[3]); // Second row - trkJac(1, 0) = -trkJac(0, 1) / tanTheta; - trkJac(1, 1) = trkJac(0, 0) / tanTheta; + freeToBoundJac(1, 0) = -freeToBoundJac(0, 1) / tanTheta; + freeToBoundJac(1, 1) = freeToBoundJac(0, 0) / tanTheta; // Dimension of the part of the Jacobian that is an identity matrix constexpr unsigned int nDimIdentity = nFreeParams - 2; - trkJac.template block(1, 2) = + freeToBoundJac.template block(1, 2) = ActsMatrix::Identity(); // Full perigee track covariance - BoundMatrix fullPerTrackCov(BoundMatrix::Identity()); - fullPerTrackCov.block(0, 0) = - (trkJac * (fullTrkCov * trkJac.transpose())); + BoundMatrix boundTrackCov(BoundMatrix::Identity()); + boundTrackCov.block(0, 0) = + (freeToBoundJac * (freeTrkCov * freeToBoundJac.transpose())); - return fullPerTrackCov; + return boundTrackCov; } From c8343f7734adcf3a8a382c200ea5aee3a76a501e Mon Sep 17 00:00:00 2001 From: russofel Date: Thu, 30 Nov 2023 13:45:11 +0100 Subject: [PATCH 07/12] add unit test --- .../AdaptiveMultiVertexFitterTests.cpp | 140 ++++++++++++++++++ 1 file changed, 140 insertions(+) diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp index a6a01858ba2..9dc2c43b111 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp @@ -82,12 +82,18 @@ std::uniform_real_distribution phiDist(-M_PI, M_PI); std::uniform_real_distribution thetaDist(1.0, M_PI - 1.0); // Track charge helper distribution std::uniform_real_distribution qDist(-1, 1); +// Distribution of track time (relative to vertex time). Values are unrealistic +// and only used for testing purposes. +std::uniform_real_distribution relTDist(-4_ps, 4_ps); // Track IP resolution distribution std::uniform_real_distribution resIPDist(0., 100._um); // Track angular distribution std::uniform_real_distribution resAngDist(0., 0.1); // Track q/p resolution distribution std::uniform_real_distribution resQoPDist(-0.1, 0.1); +// Track time resolution distribution. Values are unrealistic and only used for +// testing purposes. +std::uniform_real_distribution resTDist(0_ps, 8_ps); /// @brief Unit test for AdaptiveMultiVertexFitter /// @@ -314,6 +320,140 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test) { << vtxList.at(2).fullPosition()); } +/// @brief Unit test for fitting a 4D vertex position +/// +BOOST_AUTO_TEST_CASE(time_fitting) { + // Set up RNG + int mySeed = 31415; + std::mt19937 gen(mySeed); + + // Set up constant B-Field + auto bField = std::make_shared(Vector3{0.0, 0.0, 1_T}); + + // Set up EigenStepper + EigenStepper<> stepper(bField); + + // Set up propagator with void navigator + auto propagator = std::make_shared(stepper); + + VertexingOptions vertexingOptions(geoContext, + magFieldContext); + + // IP 3D Estimator + using IPEstimator = ImpactPointEstimator; + + IPEstimator::Config ip3dEstCfg(bField, propagator); + IPEstimator ip3dEst(ip3dEstCfg); + + AdaptiveMultiVertexFitter::Config fitterCfg( + ip3dEst); + + // Linearizer for BoundTrackParameters type test + Linearizer::Config ltConfig(bField, propagator); + Linearizer linearizer(ltConfig); + + // Test smoothing + fitterCfg.doSmoothing = true; + // Do time fit + fitterCfg.useTime = true; + + AdaptiveMultiVertexFitter fitter( + std::move(fitterCfg)); + + // Vertex position + double trueVtxTime = 40.0_ps; + Vector3 trueVtxPos(-0.15_mm, -0.1_mm, -1.5_mm); + + // Seed position of the vertex + Vector4 vtxSeedPos(0.0_mm, 0.0_mm, -1.4_mm, 0.0_ps); + + Vertex vtx(vtxSeedPos); + // Set initial covariance matrix to a large value + SquareMatrix4 initialCovariance(SquareMatrix4::Identity() * 1e+8); + vtx.setFullCovariance(initialCovariance); + + // Vector of all tracks + std::vector trks; + + unsigned int nTracks = 4; + for (unsigned int _ = 0; _ < nTracks; _++) { + // Construct positive or negative charge randomly + double q = qDist(gen) < 0 ? -1. : 1.; + + // Track resolution + double resD0 = resIPDist(gen); + double resZ0 = resIPDist(gen); + double resPh = resAngDist(gen); + double resTh = resAngDist(gen); + double resQp = resQoPDist(gen); + double resT = resTDist(gen); + + // Random diagonal covariance matrix + Covariance covMat; + + covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0., + 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0., + 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., resT * resT; + + // Random track parameters + BoundTrackParameters::ParametersVector paramVec; + paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen), + q / pTDist(gen), trueVtxTime + relTDist(gen); + + std::shared_ptr perigeeSurface = + Surface::makeShared(trueVtxPos); + + trks.emplace_back(perigeeSurface, paramVec, std::move(covMat), + ParticleHypothesis::pion()); + } + + std::vector trksPtr; + for (const auto& trk : trks) { + trksPtr.push_back(&trk); + } + + // Prepare fitter state + AdaptiveMultiVertexFitter::State state( + *bField, magFieldContext); + + for (const auto& trk : trks) { + ACTS_DEBUG("Track parameters:\n" << trk); + // Index of current vertex + state.vtxInfoMap[&vtx].trackLinks.push_back(&trk); + state.tracksAtVerticesMap.insert( + std::make_pair(std::make_pair(&trk, &vtx), + TrackAtVertex(1., trk, &trk))); + } + + state.addVertexToMultiMap(vtx); + + auto res = fitter.addVtxToFit(state, vtx, linearizer, vertexingOptions); + + BOOST_CHECK(res.ok()); + + ACTS_DEBUG("Truth vertex position: " << trueVtxPos.transpose()); + ACTS_DEBUG("Fitted vertex position: " << vtx.position().transpose()); + + ACTS_DEBUG("Truth vertex time: " << trueVtxTime); + ACTS_DEBUG("Fitted vertex time: " << vtx.time()); + + // Check that true vertex position and time approximately recovered + CHECK_CLOSE_ABS(trueVtxPos, vtx.position(), 60_um); + CHECK_CLOSE_ABS(trueVtxTime, vtx.time(), 2_ps); + + const SquareMatrix4& vtxCov = vtx.fullCovariance(); + + ACTS_DEBUG("Vertex 4D covariance after the fit:\n" << vtxCov); + + // Check that variances of the vertex position/time are positive + for (std::size_t i = 0; i <= 3; i++) { + BOOST_CHECK(vtxCov(i, i) > 0.); + } + + // Check that the covariance matrix is approximately symmetric + CHECK_CLOSE_ABS(vtxCov - vtxCov.transpose(), SquareMatrix4::Zero(), 1e-5); +} + /// @brief Unit test for AdaptiveMultiVertexFitter /// based on Athena unit test, i.e. same setting and /// test values are used here From 676f9bccc305522ba40240489bf6aade4acc5e8f Mon Sep 17 00:00:00 2001 From: russofel Date: Thu, 30 Nov 2023 14:00:01 +0100 Subject: [PATCH 08/12] document template params --- .../Vertexing/KalmanVertexTrackUpdater.hpp | 7 ++++++ .../Acts/Vertexing/KalmanVertexUpdater.hpp | 23 +++++++++++++++++++ 2 files changed, 30 insertions(+) diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp index 3868183ac22..e90de63943d 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp @@ -29,6 +29,10 @@ namespace KalmanVertexTrackUpdater { /// @brief Refits a single track with the knowledge of /// the vertex it has originated from /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param track Track to update /// @param vtx Vertex `track` belongs to template @@ -39,6 +43,9 @@ namespace detail { /// @brief Calculates a covariance matrix for the refitted track parameters /// +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param wMat W_k matrix from Ref. (1) /// @param crossCovVP Cross-covariance matrix between vertex position and track /// momentum diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp index d084c6e9888..c7a2b587145 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.hpp @@ -35,6 +35,8 @@ namespace KalmanVertexUpdater { /// Section 5.3.5 /// Cache object, the comments indicate the names of the variables in Ref. (1) +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). template struct Cache { using VertexVector = ActsVector; @@ -56,6 +58,10 @@ struct Cache { /// However, it does not add the track to the TrackAtVertex list. This to be /// done manually after calling the method. /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param vtx Vertex to be updated /// @param trk Track to be used for updating the vertex template @@ -65,6 +71,11 @@ void updateVertexWithTrack(Vertex& vtx, /// @brief Calculates updated vertex position and covariance as well as the /// updated track momentum when adding/removing linTrack. Saves the result in /// cache. +/// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param vtx Vertex /// @param linTrack Linearized track to be added or removed /// @param trackWeight Track weight @@ -83,6 +94,10 @@ namespace detail { /// @brief Calculates the update of the vertex position chi2 after /// adding/removing the track /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param oldVtx Vertex before the track was added/removed /// @param cache Cache containing updated vertex position /// @@ -94,6 +109,10 @@ double vertexPositionChi2Update(const Vertex& oldVtx, /// @brief Calculates chi2 of refitted track parameters /// w.r.t. updated vertex /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param linTrack Linearized version of track /// @param cache Cache containing some quantities needed in /// this function @@ -106,6 +125,10 @@ double trackParametersChi2(const LinearizedTrack& linTrack, /// @brief Adds or removes (depending on `sign`) tracks from vertex /// and updates the vertex /// +/// @tparam input_track_t Track object type +/// @tparam nDimVertex number of dimensions of the vertex. Can be 3 (if we only +/// fit its spatial coordinates) or 4 (if we also fit time). +/// /// @param vtx Vertex to be updated /// @param trk Track to be added to/removed from vtx /// @param sign +1 (add track) or -1 (remove track) From 20a33b5cc1e12d4ae7af6e2cc2429d9ba7afc99d Mon Sep 17 00:00:00 2001 From: russofel Date: Thu, 30 Nov 2023 14:03:07 +0100 Subject: [PATCH 09/12] improve comments --- Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp index b464800d617..89bb09ba2a2 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp @@ -309,7 +309,7 @@ Acts::Result Acts:: } // Update the vertex with the new track. The second template argument // corresponds to the number of fitted vertex dimensions (i.e., 3 if we - // only fit spatial coordinates and 4 if we also fit time) + // only fit spatial coordinates and 4 if we also fit time). if (m_cfg.useTime) { KalmanVertexUpdater::updateVertexWithTrack( *vtx, trkAtVtx); @@ -376,6 +376,10 @@ void Acts::AdaptiveMultiVertexFitter< for (const auto trk : state.vtxInfoMap[vtx].trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); if (trkAtVtx.trackWeight > m_cfg.minWeight) { + // Update the new track under the assumption that it originates at the + // vertex. The second template argument corresponds to the number of + // fitted vertex dimensions (i.e., 3 if we only fit spatial coordinates + // and 4 if we also fit time). if (m_cfg.useTime) { KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); } else { From f589947f3d38da25f7af78896605a15dcf7acdc5 Mon Sep 17 00:00:00 2001 From: russofel Date: Thu, 30 Nov 2023 14:05:58 +0100 Subject: [PATCH 10/12] nParams -> nBoundParams --- .../Vertexing/KalmanVertexTrackUpdater.ipp | 22 +++++----- .../Acts/Vertexing/KalmanVertexUpdater.ipp | 44 ++++++++++--------- 2 files changed, 36 insertions(+), 30 deletions(-) diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index 305fc373f49..a21309e9519 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -16,9 +16,9 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, const Vertex& vtx) { using VertexVector = ActsVector; using VertexMatrix = ActsSquareMatrix; - constexpr unsigned int nParams = nDimVertex + 2; - using ParameterVector = ActsVector; - using ParameterMatrix = ActsSquareMatrix; + constexpr unsigned int nBoundParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // Check if linearized state exists if (!track.isLinearized) { throw std::invalid_argument("TrackAtVertex object must be linearized."); @@ -36,20 +36,22 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, // Retrieve variables from the track linearization. The comments indicate the // corresponding symbol used in Ref. (1). // A_k - const ActsMatrix posJac = - linTrack.positionJacobian.block(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix momJac = - linTrack.momentumJacobian.block(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ParameterVector trkParams = linTrack.parametersAtPCA.head(); + const ParameterVector trkParams = + linTrack.parametersAtPCA.head(); // c_k - const ParameterVector constTerm = linTrack.constantTerm.head(); + const ParameterVector constTerm = linTrack.constantTerm.head(); // TODO we could use `linTrack.weightAtPCA` but only if we would always fit // time. // G_k const ParameterMatrix trkParamWeight = - linTrack.covarianceAtPCA.block(0, 0).inverse(); + linTrack.covarianceAtPCA.block(0, 0) + .inverse(); // Cache object filled with zeros KalmanVertexUpdater::Cache cache; diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp index 1cddefab3b9..b1886315fd3 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp @@ -81,21 +81,22 @@ void Acts::KalmanVertexUpdater::calculateUpdate( const Acts::Vertex& vtx, const Acts::LinearizedTrack& linTrack, const double trackWeight, const int sign, Cache& cache) { - constexpr unsigned int nParams = nDimVertex + 2; - using ParameterVector = ActsVector; - using ParameterMatrix = ActsSquareMatrix; + constexpr unsigned int nBoundParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // Retrieve variables from the track linearization. The comments indicate the // corresponding symbol used in Ref. (1). // A_k - const ActsMatrix posJac = - linTrack.positionJacobian.block(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix momJac = - linTrack.momentumJacobian.block(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ParameterVector trkParams = linTrack.parametersAtPCA.head(); + const ParameterVector trkParams = + linTrack.parametersAtPCA.head(); // c_k - const ParameterVector constTerm = linTrack.constantTerm.head(); + const ParameterVector constTerm = linTrack.constantTerm.head(); // TODO we could use `linTrack.weightAtPCA` but only if we would always fit // time. // G_k @@ -106,7 +107,8 @@ void Acts::KalmanVertexUpdater::calculateUpdate( // consider these two multiplicative factors directly in the updates of // newVertexWeight and newVertexPos. const ParameterMatrix trkParamWeight = - linTrack.covarianceAtPCA.block(0, 0).inverse(); + linTrack.covarianceAtPCA.block(0, 0) + .inverse(); // Retrieve current position of the vertex and its current weight matrix const ActsVector oldVtxPos = @@ -151,24 +153,26 @@ double Acts::KalmanVertexUpdater::detail::vertexPositionChi2Update( template double Acts::KalmanVertexUpdater::detail::trackParametersChi2( const LinearizedTrack& linTrack, const Cache& cache) { - constexpr unsigned int nParams = nDimVertex + 2; - using ParameterVector = ActsVector; - using ParameterMatrix = ActsSquareMatrix; + constexpr unsigned int nBoundParams = nDimVertex + 2; + using ParameterVector = ActsVector; + using ParameterMatrix = ActsSquareMatrix; // A_k - const ActsMatrix posJac = - linTrack.positionJacobian.block(0, 0); + const ActsMatrix posJac = + linTrack.positionJacobian.block(0, 0); // B_k - const ActsMatrix momJac = - linTrack.momentumJacobian.block(0, 0); + const ActsMatrix momJac = + linTrack.momentumJacobian.block(0, 0); // p_k - const ParameterVector trkParams = linTrack.parametersAtPCA.head(); + const ParameterVector trkParams = + linTrack.parametersAtPCA.head(); // c_k - const ParameterVector constTerm = linTrack.constantTerm.head(); + const ParameterVector constTerm = linTrack.constantTerm.head(); // TODO we could use `linTrack.weightAtPCA` but only if we would always fit // time. // G_k const ParameterMatrix trkParamWeight = - linTrack.covarianceAtPCA.block(0, 0).inverse(); + linTrack.covarianceAtPCA.block(0, 0) + .inverse(); // A_k * \tilde{x_k} const ParameterVector posJacVtxPos = posJac * cache.newVertexPos; From feb695e05d79ad7aa29506bbdcc2154e31fce8f8 Mon Sep 17 00:00:00 2001 From: russofel Date: Thu, 30 Nov 2023 14:09:21 +0100 Subject: [PATCH 11/12] throw invalid_argument if nDimVertex != 3 or 4 --- Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index a21309e9519..5a183704a75 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -14,6 +14,12 @@ template void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, const Vertex& vtx) { + if constexpr (nDimVertex != 3 && nDimVertex != 4) { + throw std::invalid_argument( + "The vertex dimension must either be 3 (when fitting the spatial " + "coordinates) or 4 (when fitting the spatial coordinates + time)."); + } + using VertexVector = ActsVector; using VertexMatrix = ActsSquareMatrix; constexpr unsigned int nBoundParams = nDimVertex + 2; From 92031bebfca3a1adcceb8914e1c42bbf47f67eab Mon Sep 17 00:00:00 2001 From: Felix Russo <72298366+felix-russo@users.noreply.github.com> Date: Thu, 30 Nov 2023 18:04:19 +0100 Subject: [PATCH 12/12] Apply suggestions from code review Co-authored-by: Alexander J. Pfleger <70842573+AJPfleger@users.noreply.github.com> --- .../Vertexing/AdaptiveMultiVertexFitterTests.cpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp index 9dc2c43b111..7b23de2f213 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp @@ -391,9 +391,15 @@ BOOST_AUTO_TEST_CASE(time_fitting) { // Random diagonal covariance matrix Covariance covMat; - covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0., - 0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0., - 0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., resT * resT; + // clang-format off + covMat << + resD0 * resD0, 0., 0., 0., 0., 0., + 0., resZ0 * resZ0, 0., 0., 0., 0., + 0., 0., resPh * resPh, 0., 0., 0., + 0., 0., 0., resTh * resTh, 0., 0., + 0., 0., 0., 0., resQp * resQp, 0., + 0., 0., 0., 0., 0., resT * resT; + // clang-format on // Random track parameters BoundTrackParameters::ParametersVector paramVec; @@ -447,7 +453,7 @@ BOOST_AUTO_TEST_CASE(time_fitting) { // Check that variances of the vertex position/time are positive for (std::size_t i = 0; i <= 3; i++) { - BOOST_CHECK(vtxCov(i, i) > 0.); + BOOST_CHECK_GT(vtxCov(i, i), 0.); } // Check that the covariance matrix is approximately symmetric