Skip to content

Commit

Permalink
perf: use analytical solution for straight tracks during impact point…
Browse files Browse the repository at this point in the history 
… estimation (#2506)

We use the Newton method to find the 3D PCA for helical tracks. Straight tracks are currently approximated by setting the helix radius to a very large value.

In this PR, the approximation is replaced by an analytical solution (which exists for straight tracks). This should lead to a better performance and (slightly) more precise results. To check the analytical solution, a unit test is added.

A derivation of the analytical solution can be found in the updated reference: [Track_Linearization_and_3D_PCA.pdf](https://github.com/acts-project/acts/files/12794276/Track_Linearization_and_3D_PCA.pdf)
  • Loading branch information
@felix-russo
felix-russo authored Oct 7, 2023
1 parent 020a076 commit d8c447c
Show file tree
Hide file tree
Showing 3 changed files with 111 additions and 37 deletions.
8 changes: 3 additions & 5 deletions Core/include/Acts/Vertexing/ImpactPointEstimator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ struct ImpactParametersAndSigma {
///
/// @brief Estimator for impact point calculations
/// A description of the underlying mathematics can be found here:
/// https://github.com/acts-project/acts/pull/2414
/// https://github.com/acts-project/acts/pull/2506
/// TODO: Upload reference at a better place
template <typename input_track_t, typename propagator_t,
typename propagator_options_t = PropagatorOptions<>>
Expand Down Expand Up @@ -81,10 +81,6 @@ class ImpactPointEstimator {
int maxIterations = 20;
/// Desired precision of deltaPhi in Newton method
double precision = 1.e-10;
/// Minimum q/p value
double minQoP = 1e-15;
/// Maximum curvature value
double maxRho = 1e+15;
};

/// @brief Constructor
Expand Down Expand Up @@ -147,6 +143,8 @@ class ImpactPointEstimator {
/// 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.
/// @note For straight tracks we use an analytical solution; for helical
/// tracks we use the Newton method.
///
/// @tparam nDim Number of dimensions used to compute compatibility
/// @note If nDim = 3 we only consider spatial dimensions; if nDim = 4, we
Expand Down
80 changes: 60 additions & 20 deletions Core/include/Acts/Vertexing/ImpactPointEstimator.ipp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,13 +241,65 @@ Acts::ImpactPointEstimator<input_track_t, propagator_t, propagator_options_t>::
// Reference point R
Vector3 refPoint = trkParams.referenceSurface().center(gctx);

// Perigee parameters (parameters of 2D PCA)
// Extract charge-related particle parameters
double absoluteCharge = trkParams.particleHypothesis().absoluteCharge();
double qOvP = trkParams.parameters()[BoundIndices::eBoundQOverP];

// Z-component of the B field at the reference position.
// Note that we assume a constant B field here!
auto fieldRes = m_cfg.bField->getField(refPoint, state.fieldCache);
if (!fieldRes.ok()) {
return fieldRes.error();
}
double bZ = (*fieldRes)[eZ];

// The particle moves on a straight trajectory if its charge is 0 or if there
// is no B field. In that case, the 3D PCA can be calculated analytically, see
// Sec 3.2 of the reference.
if (absoluteCharge == 0. or bZ == 0.) {
// Momentum direction (constant for straight tracks)
Vector3 momDirStraightTrack = trkParams.direction();

// Current position on the track
Vector3 positionOnTrack = trkParams.position(gctx);

// Distance between positionOnTrack and the 3D PCA
ActsScalar distanceToPca =
(vtxPos.template head<3>() - positionOnTrack).dot(momDirStraightTrack);

// 3D PCA
ActsVector<nDim> pcaStraightTrack;
pcaStraightTrack.template head<3>() =
positionOnTrack + distanceToPca * momDirStraightTrack;
if constexpr (nDim == 4) {
// Track time at positionOnTrack
double timeOnTrack = trkParams.parameters()[BoundIndices::eBoundTime];

ActsScalar m0 = trkParams.particleHypothesis().mass();
ActsScalar p = std::abs(absoluteCharge / qOvP);

// Speed in units of c
ActsScalar beta = p / std::hypot(p, m0);

pcaStraightTrack[3] = timeOnTrack + distanceToPca / beta;
}

// Vector pointing from the vertex position to the 3D PCA
ActsVector<nDim> deltaRStraightTrack = pcaStraightTrack - vtxPos;

return std::make_pair(deltaRStraightTrack, momDirStraightTrack);
}

// Charged particles in a constant B field follow a helical trajectory. In
// that case, we calculate the 3D PCA using the Newton method, see Sec 4.2 in
// the reference.

// Spatial Perigee parameters (i.e., spatial parameters of 2D PCA)
double d0 = trkParams.parameters()[BoundIndices::eBoundLoc0];
double z0 = trkParams.parameters()[BoundIndices::eBoundLoc1];
// Momentum angles at 2D PCA
double phiP = trkParams.parameters()[BoundIndices::eBoundPhi];
double theta = trkParams.parameters()[BoundIndices::eBoundTheta];
double qOvP = trkParams.parameters()[BoundIndices::eBoundQOverP];
double tP = trkParams.parameters()[BoundIndices::eBoundTime];
// Functions of the polar angle theta for later use
double sinTheta = std::sin(theta);
double cotTheta = 1. / std::tan(theta);
Expand All @@ -256,22 +308,8 @@ Acts::ImpactPointEstimator<input_track_t, propagator_t, propagator_options_t>::
// Note that phi corresponds to phiV in the reference.
double phi = phiP;

// B-field z-component at the reference position.
// Note that we assume a constant B field here!
auto fieldRes = m_cfg.bField->getField(refPoint, state.fieldCache);
if (!fieldRes.ok()) {
return fieldRes.error();
}
double bZ = (*fieldRes)[eZ];

// Signed radius of the helix on which the particle moves
double rho = 0;
// Curvature is infinite w/o b field
if (bZ == 0. || std::abs(qOvP) < m_cfg.minQoP) {
rho = m_cfg.maxRho;
} else {
rho = sinTheta * (1. / qOvP) / bZ;
}
double rho = sinTheta * (1. / qOvP) / bZ;

// Position of the helix center.
// We can set the z-position to a convenient value since it is not fixed by
Expand Down Expand Up @@ -308,9 +346,11 @@ Acts::ImpactPointEstimator<input_track_t, propagator_t, propagator_options_t>::
helixCenter + rho * Vector3(-sinPhi, cosPhi, -cotTheta * phi);

if constexpr (nDim == 4) {
// Time at the 2D PCA P
double tP = trkParams.parameters()[BoundIndices::eBoundTime];

ActsScalar m0 = trkParams.particleHypothesis().mass();
ActsScalar p =
std::abs(trkParams.particleHypothesis().absoluteCharge() / qOvP);
ActsScalar p = std::abs(absoluteCharge / qOvP);

// Speed in units of c
ActsScalar beta = p / std::hypot(p, m0);
Expand Down
60 changes: 48 additions & 12 deletions Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include "Acts/MagneticField/NullBField.hpp"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Propagator/StraightLineStepper.hpp"
#include "Acts/Propagator/detail/VoidPropagatorComponents.hpp"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Surfaces/Surface.hpp"
Expand Down Expand Up @@ -54,10 +55,13 @@ using namespace Acts::UnitLiterals;
using Acts::VectorHelpers::makeVector4;

using MagneticField = Acts::ConstantBField;
using StraightPropagator = Acts::Propagator<StraightLineStepper>;
using Stepper = Acts::EigenStepper<>;
using Propagator = Acts::Propagator<Stepper>;
using Estimator =
Acts::ImpactPointEstimator<Acts::BoundTrackParameters, Propagator>;
using StraightLineEstimator =
Acts::ImpactPointEstimator<Acts::BoundTrackParameters, StraightPropagator>;

const Acts::GeometryContext geoContext;
const Acts::MagneticFieldContext magFieldContext;
Expand Down Expand Up @@ -95,7 +99,7 @@ Estimator makeEstimator(double bZ) {
Estimator::Config cfg(field,
std::make_shared<Propagator>(
std::move(stepper), detail::VoidNavigator(),
getDefaultLogger("Prop", Logging::Level::VERBOSE)));
getDefaultLogger("Prop", Logging::Level::WARNING)));
return Estimator(cfg);
}

Expand Down Expand Up @@ -132,7 +136,7 @@ BOOST_AUTO_TEST_SUITE(VertexingImpactPointEstimator)

// Check `calculateDistance`, `estimate3DImpactParameters`, and
// `getVertexCompatibility`.
BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3d, tracks, d0,
BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3D, tracks, d0,
l0, t0, phi, theta, p, q) {
auto particleHypothesis = ParticleHypothesis::pion();

Expand Down Expand Up @@ -193,12 +197,24 @@ BOOST_DATA_TEST_CASE(SingleTrackDistanceParametersCompatibility3d, tracks, d0,
BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p,
q, vx0, vy0, vz0, vt0) {
using Propagator = Acts::Propagator<Stepper>;
using StraightPropagator = Acts::Propagator<StraightLineStepper>;

// Set up quantities for constant B field
auto field = std::make_shared<MagneticField>(Vector3(0, 0, 2_T));
Stepper stepper(field);
auto propagator = std::make_shared<Propagator>(std::move(stepper));
Estimator::Config cfg(field, propagator);
Estimator ipEstimator(cfg);
Estimator::State state(magFieldCache());
Estimator::State ipState(magFieldCache());

// Set up quantities for B = 0
auto zeroField = std::make_shared<MagneticField>(Vector3(0, 0, 0));
StraightLineStepper straightLineStepper;
auto straightLinePropagator =
std::make_shared<StraightPropagator>(straightLineStepper);
StraightLineEstimator::Config zeroFieldCfg(zeroField, straightLinePropagator);
StraightLineEstimator zeroFieldIPEstimator(zeroFieldCfg);
StraightLineEstimator::State zeroFieldIPState(magFieldCache());

// Vertex position and vertex object
Vector4 vtxPos(vx0, vy0, vz0, vt0);
Expand Down Expand Up @@ -240,21 +256,29 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p,
// Perigee surface at vertex position
auto refPerigeeSurface = Surface::makeShared<PerigeeSurface>(refPoint);

// Propagate to the 2D PCA of the reference point
// Set up the propagator options (they are the same with and without B field)
PropagatorOptions pOptions(geoContext, magFieldContext);
auto intersection = refPerigeeSurface
->intersect(geoContext, params.position(geoContext),
params.direction(), false)
.closest();
pOptions.direction =
Direction::fromScalarZeroAsPositive(intersection.pathLength());

// Propagate to the 2D PCA of the reference point in a constant B field
auto result = propagator->propagate(params, *refPerigeeSurface, pOptions);
BOOST_CHECK(result.ok());

const auto& refParams = *result->endParameters;

// Propagate to the 2D PCA of the reference point when B = 0
auto zeroFieldResult =
straightLinePropagator->propagate(params, *refPerigeeSurface, pOptions);
BOOST_CHECK(zeroFieldResult.ok());
const auto& zeroFieldRefParams = *zeroFieldResult->endParameters;

BOOST_TEST_CONTEXT(
"Check time at 2D PCA (i.e., function getImpactParameters)") {
"Check time at 2D PCA (i.e., function getImpactParameters) for helical "
"tracks") {
// Calculate impact parameters
auto ipParams = ipEstimator
.getImpactParameters(refParams, vtx, geoContext,
Expand All @@ -270,14 +294,14 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p,
1e-3);
}

BOOST_TEST_CONTEXT(
"Check time at 3D PCA (i.e., function getDistanceAndMomentum)") {
auto checkGetDistanceAndMomentum = [&vtxPos, &corrMomDir, corrTimeDiff](
const auto& ipe, const auto& rParams,
auto& state) {
// Find 4D distance and momentum of the track at the vertex starting from
// the perigee representation at the reference position
auto distAndMom =
ipEstimator
.getDistanceAndMomentum<4>(geoContext, refParams, vtxPos, state)
.value();
auto distAndMom = ipe.template getDistanceAndMomentum<4>(
geoContext, rParams, vtxPos, state)
.value();

ActsVector<4> distVec = distAndMom.first;
Vector3 momDir = distAndMom.second;
Expand All @@ -293,6 +317,18 @@ BOOST_DATA_TEST_CASE(TimeAtPca, tracksWithoutIPs* vertices, t0, phi, theta, p,
// Momentum direction should correspond to the momentum direction at the
// vertex
CHECK_CLOSE_OR_SMALL(momDir, corrMomDir, 1e-5, 1e-4);
};

BOOST_TEST_CONTEXT(
"Check time at 3D PCA (i.e., function getDistanceAndMomentum) for "
"straight tracks") {
checkGetDistanceAndMomentum(zeroFieldIPEstimator, zeroFieldRefParams,
zeroFieldIPState);
}
BOOST_TEST_CONTEXT(
"Check time at 3D PCA (i.e., function getDistanceAndMomentum) for "
"helical tracks") {
checkGetDistanceAndMomentum(ipEstimator, refParams, ipState);
}
}

Expand Down

0 comments on commit d8c447c

Please sign in to comment.