From ff30771aa67aac23fd62d30dc5e7eaf957f03835 Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Fri, 8 Nov 2024 16:43:19 +0100 Subject: [PATCH] refactor(gx2f): container to represent the mathematical gx2f-system (#3806) --- .../TrackFitting/GlobalChiSquareFitter.hpp | 215 ++++++++++-------- .../TrackFitting/GlobalChiSquareFitter.cpp | 13 +- 2 files changed, 131 insertions(+), 97 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index 968843fa479..c5fb8483aa7 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -283,6 +283,84 @@ struct ScatteringProperties { bool m_materialIsValid; }; +/// @brief A container to manage all properties of a gx2f system +/// +/// This struct manages the mathematical infrastructure for the gx2f. It +/// initializes and maintains the extended aMatrix and extended bVector. +struct Gx2fSystem { + public: + /// @brief Constructor to initialize matrices and vectors to zero based on specified dimensions. + /// + /// @param nDims Number of dimensions for the extended matrix and vector. + explicit Gx2fSystem(std::size_t nDims) + : m_nDims{nDims}, + m_aMatrix{Eigen::MatrixXd::Zero(nDims, nDims)}, + m_bVector{Eigen::VectorXd::Zero(nDims)} {} + + // Accessor for nDims (const reference). + std::size_t nDims() const { return m_nDims; } + + // Accessor for chi2 + double chi2() const { return m_chi2; } + + // Modifier for chi2 + double& chi2() { return m_chi2; } + + // Accessor for the matrix. + const Eigen::MatrixXd& aMatrix() const { return m_aMatrix; } + + // Accessor for a modifiable reference to the matrix. + Eigen::MatrixXd& aMatrix() { return m_aMatrix; } + + // Accessor for the vector. + const Eigen::VectorXd& bVector() const { return m_bVector; } + + // Accessor for a modifiable reference to the vector. + Eigen::VectorXd& bVector() { return m_bVector; } + + // Accessor for NDF + std::size_t ndf() const { return m_ndf; } + + // Modifier for NDF + std::size_t& ndf() { return m_ndf; } + + // It automatically deduces if we want to fit e.g. q/p and adjusts itself + // later. We have only 3 cases, because we always have l0, l1, phi, theta: + // - 4: no magnetic field -> q/p is empty + // - 5: no time measurement -> time is not fittable + // - 6: full fit + std::size_t findRequiredNdf() { + std::size_t ndfSystem = 0; + if (m_aMatrix(4, 4) == 0) { + ndfSystem = 4; + } else if (m_aMatrix(5, 5) == 0) { + ndfSystem = 5; + } else { + ndfSystem = 6; + } + + return ndfSystem; + } + + bool isWellDefined() { return m_ndf > findRequiredNdf(); } + + private: + /// Number of dimensions of the (extended) system + std::size_t m_nDims; + + /// Sum of chi-squared values. + double m_chi2 = 0.; + + /// Extended matrix for accumulation. + Eigen::MatrixXd m_aMatrix; + + /// Extended vector for accumulation. + Eigen::VectorXd m_bVector; + + /// Number of degrees of freedom of the system + std::size_t m_ndf = 0u; +}; + /// @brief Process measurements and fill the aMatrix and bVector /// /// The function processes each measurement for the GX2F Actor fitting process. @@ -292,15 +370,12 @@ struct ScatteringProperties { /// @tparam kMeasDim Number of dimensions of the measurement /// @tparam track_state_t The type of the track state /// -/// @param aMatrixExtended The aMatrix sums over the second derivatives -/// @param bVectorExtended The bVector sums over the first derivatives -/// @param chi2sum The total chi2 of the system +/// @param extendedSystem All parameters of the current equation system /// @param jacobianFromStart The Jacobian matrix from start to the current state /// @param trackState The track state to analyse /// @param logger A logger instance template -void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, - Eigen::VectorXd& bVectorExtended, double& chi2sum, +void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem, const std::vector& jacobianFromStart, const track_state_t& trackState, const Logger& logger) { @@ -317,13 +392,12 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, } // Create an extended Jacobian. This one contains only eBoundSize rows, - // because the rest is irrelevant + // because the rest is irrelevant. We fill it in the next steps // TODO make dimsExtendedParams template with unrolling - const std::size_t dimsExtendedParams = aMatrixExtended.rows(); // We create an empty Jacobian and fill it in the next steps Eigen::MatrixXd extendedJacobian = - Eigen::MatrixXd::Zero(eBoundSize, dimsExtendedParams); + Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims()); // This part of the Jacobian comes from the material-less propagation extendedJacobian.topLeftCorner() = @@ -360,13 +434,14 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, const ActsVector residual = measurement - projPredicted; // Finally contribute to chi2sum, aMatrix, and bVector - chi2sum += (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); + extendedSystem.chi2() += + (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); - aMatrixExtended += + extendedSystem.aMatrix() += (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) .eval(); - bVectorExtended += + extendedSystem.bVector() += (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) .eval() .transpose(); @@ -410,17 +485,14 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended, /// /// @tparam track_state_t The type of the track state /// -/// @param aMatrixExtended The aMatrix sums over the second derivatives -/// @param bVectorExtended The bVector sums over the first derivatives -/// @param chi2sum The total chi2 of the system +/// @param extendedSystem All parameters of the current equation system /// @param nMaterialsHandled How many materials we already handled. Used for the offset. /// @param scatteringMap The scattering map, containing all scattering angles and covariances /// @param trackState The track state to analyse /// @param logger A logger instance template void addMaterialToGx2fSums( - Eigen::MatrixXd& aMatrixExtended, Eigen::VectorXd& bVectorExtended, - double& chi2sum, const std::size_t nMaterialsHandled, + Gx2fSystem& extendedSystem, const std::size_t nMaterialsHandled, const std::unordered_map& scatteringMap, const track_state_t& trackState, const Logger& logger) { @@ -444,18 +516,18 @@ void addMaterialToGx2fSums( const ActsScalar invCov = scatteringMapId->second.invCovarianceMaterial(); // Phi contribution - aMatrixExtended(deltaPosition, deltaPosition) += + extendedSystem.aMatrix()(deltaPosition, deltaPosition) += invCov * sinThetaLoc * sinThetaLoc; - bVectorExtended(deltaPosition, 0) -= + extendedSystem.bVector()(deltaPosition, 0) -= invCov * scatteringAngles[eBoundPhi] * sinThetaLoc; - chi2sum += invCov * scatteringAngles[eBoundPhi] * sinThetaLoc * - scatteringAngles[eBoundPhi] * sinThetaLoc; + extendedSystem.chi2() += invCov * scatteringAngles[eBoundPhi] * sinThetaLoc * + scatteringAngles[eBoundPhi] * sinThetaLoc; // Theta Contribution - aMatrixExtended(deltaPosition + 1, deltaPosition + 1) += invCov; - bVectorExtended(deltaPosition + 1, 0) -= + extendedSystem.aMatrix()(deltaPosition + 1, deltaPosition + 1) += invCov; + extendedSystem.bVector()(deltaPosition + 1, 0) -= invCov * scatteringAngles[eBoundTheta]; - chi2sum += + extendedSystem.chi2() += invCov * scatteringAngles[eBoundTheta] * scatteringAngles[eBoundTheta]; ACTS_VERBOSE( @@ -495,18 +567,14 @@ void addMaterialToGx2fSums( /// @tparam track_proxy_t The type of the track proxy /// /// @param track A mutable track proxy to operate on -/// @param aMatrixExtended The aMatrix, summing over second derivatives for the fitting system -/// @param bVectorExtended The bVector, summing over first derivatives for the fitting system -/// @param chi2sum Accumulated chi2 value of the system -/// @param countNdf The number of degrees of freedom counted so far +/// @param extendedSystem All parameters of the current equation system /// @param multipleScattering Flag to consider multiple scattering in the calculation /// @param scatteringMap Map of geometry identifiers to scattering properties, containing all scattering angles and covariances /// @param geoIdVector A vector to store geometry identifiers for tracking processed elements /// @param logger A logger instance template void fillGx2fSystem( - const track_proxy_t track, Eigen::MatrixXd& aMatrixExtended, - Eigen::VectorXd& bVectorExtended, double& chi2sum, std::size_t& countNdf, + const track_proxy_t track, Gx2fSystem& extendedSystem, const bool multipleScattering, const std::unordered_map& scatteringMap, @@ -559,11 +627,11 @@ void fillGx2fSystem( "Found measurement with less than 1 or more than 6 dimension(s)."); } - countNdf += measDim; + extendedSystem.ndf() += measDim; visit_measurement(measDim, [&](auto N) { - addMeasurementToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum, - jacobianFromStart, trackState, logger); + addMeasurementToGx2fSums(extendedSystem, jacobianFromStart, + trackState, logger); }); } @@ -574,9 +642,8 @@ void fillGx2fSystem( jacobianFromStart.emplace_back(BoundMatrix::Identity()); // Add the material contribution to the system - addMaterialToGx2fSums(aMatrixExtended, bVectorExtended, chi2sum, - geoIdVector.size(), scatteringMap, trackState, - logger); + addMaterialToGx2fSums(extendedSystem, geoIdVector.size(), scatteringMap, + trackState, logger); geoIdVector.emplace_back(geoId); } @@ -592,13 +659,11 @@ void fillGx2fSystem( /// no qop/time fit) /// /// @param fullCovariancePredicted The covariance matrix to update -/// @param aMatrixExtended The matrix containing the coefficients of the linear system. -/// @param ndfSystem The number of degrees of freedom, determining the size of meaning full block +/// @param extendedSystem All parameters of the current equation system /// /// @return deltaParams The calculated delta parameters. void updateGx2fCovarianceParams(BoundMatrix& fullCovariancePredicted, - Eigen::MatrixXd& aMatrixExtended, - const std::size_t ndfSystem); + Gx2fSystem& extendedSystem); /// Global Chi Square fitter (GX2F) implementation. /// @@ -1128,8 +1193,6 @@ class Gx2Fitter { BoundVector deltaParams = BoundVector::Zero(); double chi2sum = 0; double oldChi2sum = std::numeric_limits::max(); - BoundMatrix aMatrix = BoundMatrix::Zero(); - BoundVector bVector = BoundVector::Zero(); // We need to create a temporary track container. We create several times a // new track and delete it after updating the parameters. However, if we @@ -1145,10 +1208,6 @@ class Gx2Fitter { // and used for the final track std::size_t tipIndex = Acts::MultiTrajectoryTraits::kInvalid; - // Here we will store, the ndf of the system. It automatically deduces if we - // want to fit e.g. q/p and adjusts itself later. - std::size_t ndfSystem = std::numeric_limits::max(); - // The scatteringMap stores for each visited surface their scattering // properties std::unordered_map scatteringMap; @@ -1242,9 +1301,6 @@ class Gx2Fitter { track.tipIndex() = tipIndex; track.linkForward(); - // This goes up for each measurement (for each dimension) - std::size_t countNdf = 0; - // Count the material surfaces, to set up the system. In the multiple // scattering case, we need to extend our system. std::size_t nMaterialSurfaces = 0; @@ -1278,12 +1334,9 @@ class Gx2Fitter { // dimensions for the scattering angles. const std::size_t dimsExtendedParams = eBoundSize + 2 * nMaterialSurfaces; - // Set to zero before filling - chi2sum = 0; - Eigen::MatrixXd aMatrixExtended = - Eigen::MatrixXd::Zero(dimsExtendedParams, dimsExtendedParams); - Eigen::VectorXd bVectorExtended = - Eigen::VectorXd::Zero(dimsExtendedParams); + // System that we fill with the information gathered by the actor and + // evaluate later + Gx2fSystem extendedSystem{dimsExtendedParams}; // This vector stores the IDs for each visited material in order. We use // it later for updating the scattering angles. We cannot use @@ -1291,22 +1344,10 @@ class Gx2Fitter { // all stored material in each propagation. std::vector geoIdVector; - fillGx2fSystem(track, aMatrixExtended, bVectorExtended, chi2sum, countNdf, - multipleScattering, scatteringMap, geoIdVector, - *m_addToSumLogger); - - // Get required number of degrees of freedom ndfSystem. - // We have only 3 cases, because we always have l0, l1, phi, theta - // 4: no magnetic field -> q/p is empty - // 5: no time measurement -> time not fittable - // 6: full fit - if (aMatrixExtended(4, 4) == 0) { - ndfSystem = 4; - } else if (aMatrixExtended(5, 5) == 0) { - ndfSystem = 5; - } else { - ndfSystem = 6; - } + fillGx2fSystem(track, extendedSystem, multipleScattering, scatteringMap, + geoIdVector, *m_addToSumLogger); + + chi2sum = extendedSystem.chi2(); // This check takes into account the evaluated dimensions of the // measurements. To fit, we need at least NDF+1 measurements. However, we @@ -1317,50 +1358,45 @@ class Gx2Fitter { // We skip the check during the first iteration, since we cannot guarantee // to hit all/enough measurement surfaces with the initial parameter // guess. - if ((nUpdate > 0) && (ndfSystem + 1 > countNdf)) { + if ((nUpdate > 0) && !extendedSystem.isWellDefined()) { ACTS_INFO("Not enough measurements. Require " - << ndfSystem + 1 << ", but only " << countNdf - << " could be used."); + << extendedSystem.findRequiredNdf() + 1 << ", but only " + << extendedSystem.ndf() << " could be used."); return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements; } - // get back the Bound vector components - aMatrix = aMatrixExtended.topLeftCorner().eval(); - bVector = bVectorExtended.topLeftCorner().eval(); - // calculate delta params [a] * delta = b Eigen::VectorXd deltaParamsExtended = - aMatrixExtended.colPivHouseholderQr().solve(bVectorExtended); + extendedSystem.aMatrix().colPivHouseholderQr().solve( + extendedSystem.bVector()); deltaParams = deltaParamsExtended.topLeftCorner().eval(); ACTS_VERBOSE("aMatrix:\n" - << aMatrix << "\n" + << extendedSystem.aMatrix() << "\n" << "bVector:\n" - << bVector << "\n" + << extendedSystem.bVector() << "\n" << "deltaParams:\n" << deltaParams << "\n" << "deltaParamsExtended:\n" << deltaParamsExtended << "\n" << "oldChi2sum = " << oldChi2sum << "\n" - << "chi2sum = " << chi2sum); + << "chi2sum = " << extendedSystem.chi2()); if ((gx2fOptions.relChi2changeCutOff != 0) && (nUpdate > 0) && - (std::abs(chi2sum / oldChi2sum - 1) < + (std::abs(extendedSystem.chi2() / oldChi2sum - 1) < gx2fOptions.relChi2changeCutOff)) { ACTS_INFO("Abort with relChi2changeCutOff after " << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax << " iterations."); - updateGx2fCovarianceParams(fullCovariancePredicted, aMatrixExtended, - ndfSystem); + updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem); break; } - if (chi2sum > oldChi2sum + 1e-5) { + if (extendedSystem.chi2() > oldChi2sum + 1e-5) { ACTS_DEBUG("chi2 not converging monotonically"); - updateGx2fCovarianceParams(fullCovariancePredicted, aMatrixExtended, - ndfSystem); + updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem); break; } @@ -1378,8 +1414,7 @@ class Gx2Fitter { return Experimental::GlobalChiSquareFitterError::DidNotConverge; } - updateGx2fCovarianceParams(fullCovariancePredicted, aMatrixExtended, - ndfSystem); + updateGx2fCovarianceParams(fullCovariancePredicted, extendedSystem); break; } @@ -1397,7 +1432,7 @@ class Gx2Fitter { } } - oldChi2sum = chi2sum; + oldChi2sum = extendedSystem.chi2(); } ACTS_DEBUG("Finished to iterate"); ACTS_VERBOSE("Final parameters: " << params.parameters().transpose()); diff --git a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp index 0784322308b..de61ee090d4 100644 --- a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp +++ b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp @@ -11,18 +11,17 @@ #include "Acts/Definitions/TrackParametrization.hpp" void Acts::Experimental::updateGx2fCovarianceParams( - BoundMatrix& fullCovariancePredicted, Eigen::MatrixXd& aMatrixExtended, - const std::size_t ndfSystem) { + BoundMatrix& fullCovariancePredicted, Gx2fSystem& extendedSystem) { // make invertible - for (int i = 0; i < aMatrixExtended.rows(); ++i) { - if (aMatrixExtended(i, i) == 0.) { - aMatrixExtended(i, i) = 1.; + for (std::size_t i = 0; i < extendedSystem.nDims(); ++i) { + if (extendedSystem.aMatrix()(i, i) == 0.) { + extendedSystem.aMatrix()(i, i) = 1.; } } - visit_measurement(ndfSystem, [&](auto N) { + visit_measurement(extendedSystem.findRequiredNdf(), [&](auto N) { fullCovariancePredicted.topLeftCorner() = - aMatrixExtended.inverse().topLeftCorner(); + extendedSystem.aMatrix().inverse().topLeftCorner(); }); return;