From 7a2383c6b3d121e8a7707ddddc059e95dff4a0d2 Mon Sep 17 00:00:00 2001 From: felix-russo <72298366+felix-russo@users.noreply.github.com> Date: Thu, 26 Oct 2023 04:21:54 +0200 Subject: [PATCH] fix: enforce calculation of impact parameters wrt vertex position (#2524) Previously, the function `getVertexCompatibility` could be called when the `BoundTrackParameters` were not defined wrt the vertex position. This PR checks that the origin of the track reference surface and the vertex coincide: ![image](https://github.com/acts-project/acts/assets/72298366/8ddeb4c5-e4b7-449c-85c7-1911a5433e10) This condition could not be enforced before because the track parameters were sometimes defined wrt a previous vertex position. I think this is a bug: This PR makes sure that the track parameters are always defined wrt the newest vertex thus making the check above possible. I am not sure about the CPU performance implications of this change, will need to check this. --- Core/include/Acts/Vertexing/AMVFInfo.hpp | 2 +- .../Vertexing/AdaptiveMultiVertexFitter.hpp | 9 ++- .../Vertexing/AdaptiveMultiVertexFitter.ipp | 65 ++++++++++--------- .../Acts/Vertexing/ImpactPointEstimator.ipp | 52 +++++++-------- 4 files changed, 60 insertions(+), 68 deletions(-) diff --git a/Core/include/Acts/Vertexing/AMVFInfo.hpp b/Core/include/Acts/Vertexing/AMVFInfo.hpp index 534c72cdbab..895bed78786 100644 --- a/Core/include/Acts/Vertexing/AMVFInfo.hpp +++ b/Core/include/Acts/Vertexing/AMVFInfo.hpp @@ -48,7 +48,7 @@ struct VertexInfo { // Vector of all tracks that are currently assigned to vertex std::vector trackLinks; - std::map impactParams3D; + std::map impactParams3D; }; } // namespace Acts diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp index c2e23bd9959..77083ca936b 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp @@ -241,15 +241,14 @@ class AdaptiveMultiVertexFitter { Vertex* vtx, const std::vector*>& verticesVec) const; - /// @brief 1) Calls ImpactPointEstimator::estimate3DImpactParameters - /// for all tracks that are associated with vtx (i.e., all elements - /// of the trackLinks vector in the VertexInfo of vtx). - /// 2) Saves the 3D impact parameters in the VertexInfo of vtx. + /// @brief Checks whether the impact parameters of the associated + /// tracks were calculated wrt the vertex position. Updates them + /// if needed. /// /// @param state Vertex fitter state /// @param vtx Vertex object /// @param vertexingOptions Vertexing options - Result prepareVertexForFit( + Result updateImpactParams3D( State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const; diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp index 8f59ac534dc..919f5812d20 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp @@ -68,9 +68,6 @@ Acts::AdaptiveMultiVertexFitter::fitImpl( if (xyDiff.norm() > m_cfg.maxDistToLinPoint) { // Set flag for relinearization vtxInfo.relinearize = true; - // Recalculate the track impact parameters at the current vertex - // position - prepareVertexForFit(state, vtx, vertexingOptions); } // Check if we use the constraint during the vertex fit @@ -124,12 +121,6 @@ Acts::AdaptiveMultiVertexFitter::addVtxToFit( std::vector*> verticesToFit; - // Save the 3D impact parameters of all tracks associated with newVertex. - auto res = prepareVertexForFit(state, &newVertex, vertexingOptions); - if (!res.ok()) { - return res.error(); - } - // List of vertices added in last iteration std::vector*> lastIterAddedVertices = {&newVertex}; // List of vertices added in current iteration @@ -189,25 +180,47 @@ bool Acts::AdaptiveMultiVertexFitter:: } template -Acts::Result Acts:: - AdaptiveMultiVertexFitter::prepareVertexForFit( +Acts::Result +Acts::AdaptiveMultiVertexFitter:: + updateImpactParams3D( State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const { // Vertex info object auto& vtxInfo = state.vtxInfoMap[vtx]; - // Vertex seed position - const Vector3& seedPos = vtxInfo.seedPosition.template head<3>(); + // Vertex position, i.e., point wrt which the impact parameters are estimated + const Vector3& vtxPosition = vtxInfo.oldPosition.template head<3>(); // Loop over all tracks at the vertex for (const auto& trk : vtxInfo.trackLinks) { + // Track parameters + auto trkParams = m_extractParameters(*trk); + + // Origin of the track reference surface + Vector3 surfaceOrigin = trkParams.referenceSurface() + .transform(vertexingOptions.geoContext) + .translation(); + // Skip the impact point estimation if the impact parameters of trk wrt the + // current vertex position were already calculated. + if (surfaceOrigin == vtxPosition && + vtxInfo.impactParams3D.find(trk) != vtxInfo.impactParams3D.end()) { + continue; + } + auto res = m_cfg.ipEst.estimate3DImpactParameters( vertexingOptions.geoContext, vertexingOptions.magFieldContext, - m_extractParameters(*trk), seedPos, state.ipState); + trkParams, vtxPosition, state.ipState); if (!res.ok()) { return res.error(); } - // Save 3D impact parameters of the track - vtxInfo.impactParams3D.emplace(trk, res.value()); + + // Try to create a new map entry. If "trk" already exists as key, the + // corresponding value will not be overwritten and the boolean "inserted" + // will be set to false ... + auto [it, inserted] = vtxInfo.impactParams3D.emplace(trk, *res); + // ... and we have to overwrite manually. + if (not inserted) { + it->second = *res; + } } return {}; } @@ -218,26 +231,14 @@ Acts::AdaptiveMultiVertexFitter:: setAllVertexCompatibilities( State& state, Vertex* vtx, const VertexingOptions& vertexingOptions) const { - VertexInfo& vtxInfo = state.vtxInfoMap[vtx]; + // Update the 3D impact parameters of all tracks + updateImpactParams3D(state, vtx, vertexingOptions); - // Loop over all tracks that are associated with vtx and estimate their + VertexInfo& vtxInfo = state.vtxInfoMap[vtx]; + // Loop over the tracks that are associated with vtx and estimate their // compatibility for (const auto& trk : vtxInfo.trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); - // Recover from cases where linearization point != 0 but - // more tracks were added later on - if (vtxInfo.impactParams3D.find(trk) == vtxInfo.impactParams3D.end()) { - auto res = m_cfg.ipEst.estimate3DImpactParameters( - vertexingOptions.geoContext, vertexingOptions.magFieldContext, - m_extractParameters(*trk), VectorHelpers::position(vtxInfo.linPoint), - state.ipState); - if (!res.ok()) { - return res.error(); - } - // Set impactParams3D for current trackAtVertex - vtxInfo.impactParams3D.emplace(trk, res.value()); - } - // Set compatibility with current vertex Acts::Result compatibilityResult(0.); if (m_cfg.useTime) { compatibilityResult = m_cfg.ipEst.template getVertexCompatibility<4>( diff --git a/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp b/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp index 1fc3b2830d7..c2c3b6a18c5 100644 --- a/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp +++ b/Core/include/Acts/Vertexing/ImpactPointEstimator.ipp @@ -116,43 +116,22 @@ Acts::ImpactPointEstimator:: return VertexingError::EmptyInput; } - // Retrieve weight matrix of the track's local x-, y-, and time-coordinate - // (the latter only if nDim = 4). For this, the covariance needs to be set. - if (not trkParams->covariance().has_value()) { - return VertexingError::NoCovariance; - } - ActsSquareMatrix subCovMat; - if constexpr (nDim == 3) { - subCovMat = trkParams->spatialImpactParameterCovariance().value(); - } else { - subCovMat = trkParams->impactParameterCovariance().value(); - } - ActsSquareMatrix weight = subCovMat.inverse(); - - // Orientation of the surface (i.e., axes of the corresponding coordinate - // system) - RotationMatrix3 surfaceAxes = - trkParams->referenceSurface().transform(gctx).rotation(); // Origin of the surface coordinate system Vector3 surfaceOrigin = trkParams->referenceSurface().transform(gctx).translation(); - // x- and y-axis of the surface coordinate system - Vector3 xAxis = surfaceAxes.col(0); - Vector3 yAxis = surfaceAxes.col(1); - - // Vector pointing from the surface origin to the vertex position - // TODO: The vertex should always be at the surfaceOrigin since the - // track parameters should be obtained by estimate3DImpactParameters. - // Therefore, originToVertex should always be 0, which is currently not the - // case. - Vector3 originToVertex = vertexPos.template head<3>() - surfaceOrigin; + if (vertexPos.template head<3>() != surfaceOrigin) { + throw std::invalid_argument( + "The origin of the track reference surface must be at the vertex " + "position. Please make sure that trkParams was computed using " + "estimate3DImpactParameters."); + } // x-, y-, and possibly time-coordinate of the vertex and the track in the - // surface coordinate system + // surface coordinate system. Since the vertex is at the surface origin, its + // spatial coordinates are 0. ActsVector localVertexCoords; - localVertexCoords.template head<2>() = - Vector2(originToVertex.dot(xAxis), originToVertex.dot(yAxis)); + localVertexCoords.template head<2>() = Vector2(0., 0.); ActsVector localTrackCoords; localTrackCoords.template head<2>() = @@ -167,6 +146,19 @@ Acts::ImpactPointEstimator:: // residual ActsVector residual = localTrackCoords - localVertexCoords; + // Retrieve weight matrix of the track's local x-, y-, and time-coordinate + // (the latter only if nDim = 4). For this, the covariance needs to be set. + if (not trkParams->covariance().has_value()) { + return VertexingError::NoCovariance; + } + ActsSquareMatrix subCovMat; + if constexpr (nDim == 3) { + subCovMat = trkParams->spatialImpactParameterCovariance().value(); + } else { + subCovMat = trkParams->impactParameterCovariance().value(); + } + ActsSquareMatrix weight = subCovMat.inverse(); + // return chi2 return residual.dot(weight * residual); }