Skip to content

Commit

Permalink
feat: Add time to impact point estimation (#2414)
Browse files Browse the repository at this point in the history 
In this PR I:
- add the option of calculating the time at the 2D and the 3D point of closest approach (PCA) of a track to a vertex.
- implement unit tests for the above feature

If further refactor `ImpactPointEstimator` by
- providing a reference ([Track_Linearization_and_3D_PCA.pdf](https://github.com/acts-project/acts/files/12513641/Track_Linearization_and_3D_PCA.pdf)) 
- changing variable names to match the reference
- removing unused pieces of code
- adding explanatory comments
- adding TODO's where I think we could improve the code further
  • Loading branch information
@felix-russo
felix-russo authored Sep 29, 2023
1 parent 753dbec commit 8b146dc
Show file tree
Hide file tree
Showing 8 changed files with 539 additions and 290 deletions.
8 changes: 4 additions & 4 deletions Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ auto Acts::AdaptiveMultiVertexFinder<vfitter_t, sfinder_t>::getIPSignificance(
newVtx.setFullCovariance(SquareMatrix4::Zero());
}

auto estRes = m_cfg.ipEstimator.estimateImpactParameters(
auto estRes = m_cfg.ipEstimator.getImpactParameters(
m_extractParameters(*track), newVtx, vertexingOptions.geoContext,
vertexingOptions.magFieldContext);
if (!estRes.ok()) {
Expand All @@ -213,9 +213,9 @@ auto Acts::AdaptiveMultiVertexFinder<vfitter_t, sfinder_t>::getIPSignificance(
ImpactParametersAndSigma ipas = *estRes;

double significance = 0.;
if (ipas.sigmad0 > 0 && ipas.sigmaz0 > 0) {
significance = std::sqrt(std::pow(ipas.IPd0 / ipas.sigmad0, 2) +
std::pow(ipas.IPz0 / ipas.sigmaz0, 2));
if (ipas.sigmaD0 > 0 && ipas.sigmaZ0 > 0) {
significance = std::sqrt(std::pow(ipas.d0 / ipas.sigmaD0, 2) +
std::pow(ipas.z0 / ipas.sigmaZ0, 2));
}

return significance;
Expand Down
2 changes: 1 addition & 1 deletion Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,7 +239,7 @@ Acts::AdaptiveMultiVertexFitter<input_track_t, linearizer_t>::
currentVtxInfo.ip3dParams.emplace(trk, res.value());
}
// Set compatibility with current vertex
auto compRes = m_cfg.ipEst.get3dVertexCompatibility(
auto compRes = m_cfg.ipEst.template getVertexCompatibility<3>(
vertexingOptions.geoContext, &(currentVtxInfo.ip3dParams.at(trk)),
VectorHelpers::position(currentVtxInfo.oldPosition));
if (!compRes.ok()) {
Expand Down
154 changes: 85 additions & 69 deletions Core/include/Acts/Vertexing/ImpactPointEstimator.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
// This file is part of the Acts project.
//
// Copyright (C) 2019 CERN for the benefit of the Acts project
// Copyright (C) 2019-2023 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
Expand All @@ -21,20 +21,27 @@
namespace Acts {

struct ImpactParametersAndSigma {
double IPd0 = 0.;
double IPz0 = 0.;
double IPz0SinTheta = 0.;
double sigmad0 = 0.;
double sigmaz0 = 0.;
double sigmaz0SinTheta = 0.;
double PVsigmad0 = 0.;
double PVsigmaz0 = 0.;
double PVsigmaz0SinTheta = 0.;
// Impact parameters ...
double d0 = 0.;
double z0 = 0.;
// ... and their standard deviations wrt a vertex, e.g.:
// sigmaD0 = sqrt(Var(X) + Var(Y) + Var(d0)),
// where X and Y are the x- and y-coordinate of the vertex
double sigmaD0 = 0.;
double sigmaZ0 = 0.;
// Absolute difference in time between the vertex and the track at the 2D PCA
// ...
std::optional<double> deltaT = std::nullopt;
// ... and standard deviation wrt a vertex
std::optional<double> sigmaDeltaT = std::nullopt;
};

/// @class ImpactPointEstimator
///
/// @brief Estimator for impact point calculations
/// A description of the underlying mathematics can be found here:
/// https://github.com/acts-project/acts/pull/2414
/// TODO: Upload reference at a better place
template <typename input_track_t, typename propagator_t,
typename propagator_options_t = PropagatorOptions<>>
class ImpactPointEstimator {
Expand Down Expand Up @@ -72,7 +79,7 @@ class ImpactPointEstimator {
std::shared_ptr<const propagator_t> propagator;
/// Max. number of iterations in Newton method
int maxIterations = 20;
/// Desired precision in deltaPhi in Newton method
/// Desired precision of deltaPhi in Newton method
double precision = 1.e-10;
/// Minimum q/p value
double minQoP = 1e-15;
Expand All @@ -85,63 +92,88 @@ class ImpactPointEstimator {
/// @param cfg Configuration object
ImpactPointEstimator(const Config& cfg) : m_cfg(cfg) {}

/// @brief Calculates 3D distance between a track and a 3D point
/// @brief Calculates 3D distance between a track and a vertex
///
/// @param gctx The geometry context
/// @param trkParams Track parameters
/// @param vtxPos Position to calculate distance to
/// @param vtxPos 3D position to calculate the distance to
/// @param state The state object
///
/// @return Distance
Result<double> calculate3dDistance(const GeometryContext& gctx,
const BoundTrackParameters& trkParams,
const Vector3& vtxPos, State& state) const;

/// @brief Creates track parameters bound to plane
/// at point of closest approach in 3d to given
/// reference position. The parameters and errors
/// are defined on the plane intersecting the track
/// at point of closest approach, with track orthogonal
/// to the plane and center of the plane defined as the
/// given reference point (vertex).
Result<double> calculateDistance(const GeometryContext& gctx,
const BoundTrackParameters& trkParams,
const Vector3& vtxPos, State& state) const;

/// @brief Estimates the track parameters at the 3D PCA (i.e., a point of
/// minimal 3D distance) to a vertex. The track parameters are defined wrt a
/// reference plane that has its origin at the vertex position and whose
/// z-axis points in the direction of the track momentum. The plane's x-axis
/// points approximately from the vertex to the 3D PCA (it is only approximate
/// because we force it to be orthogonal to the z-axis). The y-axis is
/// calculated as a cross product between x- and z-axis.
///
/// @param gctx The geometry context
/// @param mctx The magnetic field context
/// @param trkParams Track parameters
/// @param vtxPos Reference position (vertex)
/// @param state The state object
///
/// @return New track params
/// @return Track parameters at the 3D PCA
Result<BoundTrackParameters> estimate3DImpactParameters(
const GeometryContext& gctx, const Acts::MagneticFieldContext& mctx,
const BoundTrackParameters& trkParams, const Vector3& vtxPos,
State& state) const;

/// @brief Estimates the compatibility of a
/// track to a vertex position based on the 3d
/// distance between the track and the vertex
/// @brief Estimates the compatibility of a track to a vertex based on their
/// 3D (if nDim = 3) or 4D (if nDim = 4) distance and the track covariance.
/// @note Confusingly, a *smaller* compatibility means that a track is *more*
/// compatible.
///
/// @tparam nDim Number of dimensions used to compute compatibility
/// @note If nDim = 3 we only consider spatial dimensions; if nDim = 4, we
/// also consider time. Other values are not allowed.
/// @param gctx The Geometry context
/// @param trkParams Track parameters at point of closest
/// approach in 3d as retrieved by estimate3DImpactParameters
/// approach in 3D as retrieved by estimate3DImpactParameters
/// @param vertexPos The vertex position
///
/// @return The compatibility value
Result<double> get3dVertexCompatibility(const GeometryContext& gctx,
const BoundTrackParameters* trkParams,
const Vector3& vertexPos) const;

/// @brief Estimates the impact parameters and their errors of a given
/// track w.r.t. a vertex by propagating the trajectory state
/// towards the vertex position.
template <unsigned int nDim>
Result<double> getVertexCompatibility(
const GeometryContext& gctx, const BoundTrackParameters* trkParams,
const ActsVector<nDim>& vertexPos) const;

/// @brief Calculate the distance between a track and a vertex by finding the
/// corresponding 3D PCA. Returns also the momentum direction at the 3D PCA.
/// The template parameter nDim determines whether we calculate the 3D
/// distance (nDim = 3) or the 4D distance (nDim = 4) to the 3D PCA.
///
/// @tparam nDim Number of dimensions used to compute compatibility
/// @note If nDim = 3 we only consider spatial dimensions; if nDim = 4, we
/// also consider time. Other values are not allowed.
/// @param gctx Geometry context
/// @param trkParams Track parameters
/// @param vtxPos Vertex position
/// @param state The state object
template <unsigned int nDim>
Result<std::pair<Acts::ActsVector<nDim>, Acts::Vector3>>
getDistanceAndMomentum(const GeometryContext& gctx,
const BoundTrackParameters& trkParams,
const ActsVector<nDim>& vtxPos, State& state) const;

/// @brief Calculates the impact parameters of a track w.r.t. a vertex. The
/// corresponding errors are approximated by summing the variances of the
/// track and the vertex.
///
/// @param track Track to estimate IP from
/// @param vtx Vertex the track belongs to
/// @param track Track whose impact parameters are calculated
/// @param vtx Vertex corresponding to the track
/// @param gctx The geometry context
/// @param mctx The magnetic field context
Result<ImpactParametersAndSigma> estimateImpactParameters(
/// @param calculateTimeIP If true, the difference in time is computed
Result<ImpactParametersAndSigma> getImpactParameters(
const BoundTrackParameters& track, const Vertex<input_track_t>& vtx,
const GeometryContext& gctx, const MagneticFieldContext& mctx) const;
const GeometryContext& gctx, const MagneticFieldContext& mctx,
bool calculateTimeIP = false) const;

/// @brief Estimates the sign of the 2D and Z lifetime of a given track
/// w.r.t. a vertex and a direction (e.g. a jet direction)
Expand All @@ -154,8 +186,8 @@ class ImpactPointEstimator {
/// @param gctx The geometry context
/// @param mctx The magnetic field context
///
/// @return A pair holding the sign for the 2D an Z lifetimes
Result<std::pair<double, double>> getLifetimesSignOfTrack(
/// @return A pair holding the sign for the 2D and Z lifetimes
Result<std::pair<double, double>> getLifetimeSignOfTrack(
const BoundTrackParameters& track, const Vertex<input_track_t>& vtx,
const Acts::Vector3& direction, const GeometryContext& gctx,
const MagneticFieldContext& mctx) const;
Expand All @@ -182,34 +214,18 @@ class ImpactPointEstimator {
/// @brief Performs a Newton approximation to retrieve a point
/// of closest approach in 3D to a reference position
///
/// @param trkPos Initial position
/// @param vtxPos Reference position
/// @param phi Phi along the helix which will be changed by
/// the Newton method
/// @param theta Track theta
/// @param r Helix radius
/// @param helixCenter Position of the helix center
/// @param vtxPos Vertex position
/// @param phi Azimuthal momentum angle
/// @note Modifying phi corresponds to moving along the track. This function
/// optimizes phi until we reach a 3D PCA.
/// @param theta Polar momentum angle (constant along the track)
/// @param rho Signed helix radius
///
/// @return New phi value
Result<double> performNewtonApproximation(const Vector3& trkPos,
const Vector3& vtxPos, double phi,
double theta, double r) const;

/// @brief Helper function to calculate relative
/// distance between track and vtxPos and the
/// direction of the momentum
///
/// @param gctx The geometry context
/// @param trkParams Track parameters
/// @param vtxPos The vertex position
/// @param deltaR Relative position between
/// track and vtxPos, to be determined by method
/// @param momDir Momentum direction, to be
/// determined by method
/// @param state The state object
Result<void> getDistanceAndMomentum(const GeometryContext& gctx,
const BoundTrackParameters& trkParams,
const Vector3& vtxPos, Vector3& deltaR,
Vector3& momDir, State& state) const;
/// @return Phi value at 3D PCA
Result<double> performNewtonOptimization(const Vector3& helixCenter,
const Vector3& vtxPos, double phi,
double theta, double rho) const;
};

} // namespace Acts
Expand Down
Loading

0 comments on commit 8b146dc

Please sign in to comment.