Skip to content

Commit

Permalink
fix: enforce calculation of impact parameters wrt vertex position (#2524
Browse files Browse the repository at this point in the history 
)

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.
  • Loading branch information
@felix-russo
felix-russo authored Oct 26, 2023
1 parent 8ce9fd3 commit 7a2383c
Show file tree
Hide file tree
Showing 4 changed files with 60 additions and 68 deletions.
2 changes: 1 addition & 1 deletion Core/include/Acts/Vertexing/AMVFInfo.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ struct VertexInfo {
// Vector of all tracks that are currently assigned to vertex
std::vector<const input_track_t*> trackLinks;

std::map<const input_track_t*, const BoundTrackParameters> impactParams3D;
std::map<const input_track_t*, BoundTrackParameters> impactParams3D;
};

} // namespace Acts
9 changes: 4 additions & 5 deletions Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,15 +241,14 @@ class AdaptiveMultiVertexFitter {
Vertex<InputTrack_t>* vtx,
const std::vector<Vertex<InputTrack_t>*>& 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<void> prepareVertexForFit(
Result<void> updateImpactParams3D(
State& state, Vertex<InputTrack_t>* vtx,
const VertexingOptions<InputTrack_t>& vertexingOptions) const;

Expand Down
65 changes: 33 additions & 32 deletions Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,9 +68,6 @@ Acts::AdaptiveMultiVertexFitter<input_track_t, linearizer_t>::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
Expand Down Expand Up @@ -124,12 +121,6 @@ Acts::AdaptiveMultiVertexFitter<input_track_t, linearizer_t>::addVtxToFit(

std::vector<Vertex<input_track_t>*> 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<Vertex<input_track_t>*> lastIterAddedVertices = {&newVertex};
// List of vertices added in current iteration
Expand Down Expand Up @@ -189,25 +180,47 @@ bool Acts::AdaptiveMultiVertexFitter<input_track_t, linearizer_t>::
}

template <typename input_track_t, typename linearizer_t>
Acts::Result<void> Acts::
AdaptiveMultiVertexFitter<input_track_t, linearizer_t>::prepareVertexForFit(
Acts::Result<void>
Acts::AdaptiveMultiVertexFitter<input_track_t, linearizer_t>::
updateImpactParams3D(
State& state, Vertex<input_track_t>* vtx,
const VertexingOptions<input_track_t>& 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 {};
}
Expand All @@ -218,26 +231,14 @@ Acts::AdaptiveMultiVertexFitter<input_track_t, linearizer_t>::
setAllVertexCompatibilities(
State& state, Vertex<input_track_t>* vtx,
const VertexingOptions<input_track_t>& vertexingOptions) const {
VertexInfo<input_track_t>& 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<input_track_t>& 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<double> compatibilityResult(0.);
if (m_cfg.useTime) {
compatibilityResult = m_cfg.ipEst.template getVertexCompatibility<4>(
Expand Down
52 changes: 22 additions & 30 deletions Core/include/Acts/Vertexing/ImpactPointEstimator.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,43 +116,22 @@ Acts::ImpactPointEstimator<input_track_t, propagator_t, propagator_options_t>::
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<nDim - 1> subCovMat;
if constexpr (nDim == 3) {
subCovMat = trkParams->spatialImpactParameterCovariance().value();
} else {
subCovMat = trkParams->impactParameterCovariance().value();
}
ActsSquareMatrix<nDim - 1> 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<nDim - 1> localVertexCoords;
localVertexCoords.template head<2>() =
Vector2(originToVertex.dot(xAxis), originToVertex.dot(yAxis));
localVertexCoords.template head<2>() = Vector2(0., 0.);

ActsVector<nDim - 1> localTrackCoords;
localTrackCoords.template head<2>() =
Expand All @@ -167,6 +146,19 @@ Acts::ImpactPointEstimator<input_track_t, propagator_t, propagator_options_t>::
// residual
ActsVector<nDim - 1> 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<nDim - 1> subCovMat;
if constexpr (nDim == 3) {
subCovMat = trkParams->spatialImpactParameterCovariance().value();
} else {
subCovMat = trkParams->impactParameterCovariance().value();
}
ActsSquareMatrix<nDim - 1> weight = subCovMat.inverse();

// return chi2
return residual.dot(weight * residual);
}
Expand Down

0 comments on commit 7a2383c

Please sign in to comment.