Skip to content

Commit

Permalink
Merge branch 'main' into refactor/avoid-new-delete
Browse files Browse the repository at this point in the history
  • Loading branch information
paulgessinger authored Nov 22, 2024
2 parents 0f71f79 + c3c4bf5 commit c62c22c
Show file tree
Hide file tree
Showing 7 changed files with 55 additions and 163 deletions.
172 changes: 35 additions & 137 deletions Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,17 +8,14 @@

#pragma once

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Acts/Utilities/MathHelpers.hpp"
#include "Acts/Utilities/Zip.hpp"

#include <array>
#include <cmath>
#include <iterator>
#include <optional>
#include <stdexcept>

Expand All @@ -27,17 +24,16 @@ namespace Acts {
/// Estimate the full track parameters from three space points
///
/// This method is based on the conformal map transformation. It estimates the
/// full free track parameters, i.e. (x, y, z, t, dx, dy, dz, q/p) at the
/// bottom space point. The bottom space is assumed to be the first element
/// in the range defined by the iterators. The magnetic field (which might be
/// along any direction) is also necessary for the momentum estimation.
/// full free track parameters, i.e. (x, y, z, t, dx, dy, dz, q/p) at the bottom
/// space point. The bottom space is assumed to be the first element in the
/// range defined by the iterators. The magnetic field (which might be along any
/// direction) is also necessary for the momentum estimation.
///
/// This is a purely spatial estimation, i.e. the time parameter will be set to
/// 0.
///
/// It resembles the method used in ATLAS for the track parameters
/// estimated from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane
/// here:
/// It resembles the method used in ATLAS for the track parameters estimated
/// from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane here:
/// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
///
/// @tparam spacepoint_iterator_t The type of space point iterator
Expand All @@ -55,14 +51,13 @@ FreeVector estimateTrackParamsFromSeed(const Vector3& sp0, const Vector3& sp1,
/// Estimate the full track parameters from three space points
///
/// This method is based on the conformal map transformation. It estimates the
/// full free track parameters, i.e. (x, y, z, t, dx, dy, dz, q/p) at the
/// bottom space point. The bottom space is assumed to be the first element
/// in the range defined by the iterators. The magnetic field (which might be
/// along any direction) is also necessary for the momentum estimation.
/// full free track parameters, i.e. (x, y, z, t, dx, dy, dz, q/p) at the bottom
/// space point. The bottom space is assumed to be the first element in the
/// range defined by the iterators. The magnetic field (which might be along any
/// direction) is also necessary for the momentum estimation.
///
/// It resembles the method used in ATLAS for the track parameters
/// estimated from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane
/// here:
/// It resembles the method used in ATLAS for the track parameters estimated
/// from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane here:
/// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
///
/// @tparam spacepoint_iterator_t The type of space point iterator
Expand Down Expand Up @@ -107,150 +102,53 @@ FreeVector estimateTrackParamsFromSeed(spacepoint_range_t spRange,
/// This method is based on the conformal map transformation. It estimates the
/// full bound track parameters, i.e. (loc0, loc1, phi, theta, q/p, t) at the
/// bottom space point. The bottom space is assumed to be the first element
/// in the range defined by the iterators. It must lie on the surface
/// provided for the representation of the bound track parameters. The magnetic
/// field (which might be along any direction) is also necessary for the
/// momentum estimation.
/// in the range defined by the iterators. It must lie on the surface provided
/// for the representation of the bound track parameters. The magnetic field
/// (which might be along any direction) is also necessary for the momentum
/// estimation.
///
/// It resembles the method used in ATLAS for the track parameters
/// estimated from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane
/// here:
/// It resembles the method used in ATLAS for the track parameters estimated
/// from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane here:
/// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
///
/// @tparam spacepoint_iterator_t The type of space point iterator
///
/// @param gctx is the geometry context
/// @param spBegin is the begin iterator for the space points
/// @param spEnd is the end iterator for the space points
/// @param spRange is the range of space points
/// @param surface is the surface of the bottom space point. The estimated bound
/// track parameters will be represented also at this surface
/// @param bField is the magnetic field vector
/// @param logger A logger instance
///
/// @return optional bound parameters
template <typename spacepoint_iterator_t>
std::optional<BoundVector> estimateTrackParamsFromSeed(
const GeometryContext& gctx, spacepoint_iterator_t spBegin,
spacepoint_iterator_t spEnd, const Surface& surface, const Vector3& bField,
const Acts::Logger& logger = getDummyLogger()) {
// Check the number of provided space points
std::size_t numSP = std::distance(spBegin, spEnd);
if (numSP != 3) {
ACTS_ERROR("There should be exactly three space points provided.");
return std::nullopt;
}

// The global positions of the bottom, middle and space points
std::array<Vector3, 3> spGlobalPositions = {Vector3::Zero(), Vector3::Zero(),
Vector3::Zero()};
std::array<std::optional<float>, 3> spGlobalTimes = {
std::nullopt, std::nullopt, std::nullopt};
// The first, second and third space point are assumed to be bottom, middle
// and top space point, respectively
for (std::size_t isp = 0; isp < 3; ++isp) {
spacepoint_iterator_t it = std::next(spBegin, isp);
if (*it == nullptr) {
ACTS_ERROR("Empty space point found. This should not happen.");
return std::nullopt;
}
const auto& sp = *it;
spGlobalPositions[isp] = Vector3(sp->x(), sp->y(), sp->z());
spGlobalTimes[isp] = sp->t();
}

// Define a new coordinate frame with its origin at the bottom space point, z
// axis long the magnetic field direction and y axis perpendicular to vector
// from the bottom to middle space point. Hence, the projection of the middle
// space point on the transverse plane will be located at the x axis of the
// new frame.
Vector3 relVec = spGlobalPositions[1] - spGlobalPositions[0];
Vector3 newZAxis = bField.normalized();
Vector3 newYAxis = newZAxis.cross(relVec).normalized();
Vector3 newXAxis = newYAxis.cross(newZAxis);
RotationMatrix3 rotation;
rotation.col(0) = newXAxis;
rotation.col(1) = newYAxis;
rotation.col(2) = newZAxis;
// The center of the new frame is at the bottom space point
Translation3 trans(spGlobalPositions[0]);
// The transform which constructs the new frame
Transform3 transform(trans * rotation);

// The coordinate of the middle and top space point in the new frame
Vector3 local1 = transform.inverse() * spGlobalPositions[1];
Vector3 local2 = transform.inverse() * spGlobalPositions[2];

// In the new frame the bottom sp is at the origin, while the middle
// sp in along the x axis. As such, the x-coordinate of the circle is
// at: x-middle / 2.
// The y coordinate can be found by using the straight line passing
// between the mid point between the middle and top sp and perpendicular to
// the line connecting them
Vector2 circleCenter;
circleCenter(0) = 0.5 * local1(0);
/// @return bound parameters
template <std::ranges::range spacepoint_range_t>
Result<BoundVector> estimateTrackParamsFromSeed(const GeometryContext& gctx,
spacepoint_range_t spRange,
const Surface& surface,
const Vector3& bField) {
FreeVector freeParams = estimateTrackParamsFromSeed(spRange, bField);

ActsScalar deltaX21 = local2(0) - local1(0);
ActsScalar sumX21 = local2(0) + local1(0);
// straight line connecting the two points
// y = a * x + c (we don't care about c right now)
// we simply need the slope
// we compute 1./a since this is what we need for the following computation
ActsScalar ia = deltaX21 / local2(1);
// Perpendicular line is then y = -1/a *x + b
// we can evaluate b given we know a already by imposing
// the line passes through P = (0.5 * (x2 + x1), 0.5 * y2)
ActsScalar b = 0.5 * (local2(1) + ia * sumX21);
circleCenter(1) = -ia * circleCenter(0) + b;
// Radius is a signed distance between circleCenter and first sp, which is at
// (0, 0) in the new frame. Sign depends on the slope a (positive vs negative)
int sign = ia > 0 ? -1 : 1;
const ActsScalar R = circleCenter.norm();
ActsScalar invTanTheta =
local2.z() / (2 * R * std::asin(local2.head<2>().norm() / (2 * R)));
// The momentum direction in the new frame (the center of the circle has the
// coordinate (-1.*A/(2*B), 1./(2*B)))
ActsScalar A = -circleCenter(0) / circleCenter(1);
Vector3 transDirection(1., A, fastHypot(1, A) * invTanTheta);
// Transform it back to the original frame
Vector3 direction = rotation * transDirection.normalized();
const auto* sp0 = *spRange.begin();
Vector3 origin = Vector3(sp0->x(), sp0->y(), sp0->z());
Vector3 direction = freeParams.segment<3>(eFreeDir0);

// Initialize the bound parameters vector
BoundVector params = BoundVector::Zero();

// The estimated phi and theta
params[eBoundPhi] = VectorHelpers::phi(direction);
params[eBoundTheta] = VectorHelpers::theta(direction);
params[eBoundQOverP] = freeParams[eFreeQOverP];

// Transform the bottom space point to local coordinates of the provided
// surface
auto lpResult = surface.globalToLocal(gctx, spGlobalPositions[0], direction);
auto lpResult = surface.globalToLocal(gctx, origin, direction);
if (!lpResult.ok()) {
ACTS_ERROR(
"Global to local transformation did not succeed. Please make sure the "
"bottom space point lies on the provided surface.");
return std::nullopt;
return Result<BoundVector>::failure(lpResult.error());
}
Vector2 bottomLocalPos = lpResult.value();
// The estimated loc0 and loc1
params[eBoundLoc0] = bottomLocalPos.x();
params[eBoundLoc1] = bottomLocalPos.y();
params[eBoundTime] = spGlobalTimes[0].value_or(0.);

ActsScalar bFieldStrength = bField.norm();
// The estimated q/pt in [GeV/c]^-1 (note that the pt is the projection of
// momentum on the transverse plane of the new frame)
ActsScalar qOverPt = sign / (bFieldStrength * R);
// The estimated q/p in [GeV/c]^-1
params[eBoundQOverP] = qOverPt / fastHypot(1., invTanTheta);
params[eBoundTime] = sp0->t().value_or(0);

if (params.hasNaN()) {
ACTS_ERROR(
"The NaN value exists at the estimated track parameters from seed with"
<< "\nbottom sp: " << spGlobalPositions[0] << "\nmiddle sp: "
<< spGlobalPositions[1] << "\ntop sp: " << spGlobalPositions[2]);
return std::nullopt;
}
return params;
return Result<BoundVector>::success(params);
}

/// Configuration for the estimation of the covariance matrix of the track
Expand Down
1 change: 1 addition & 0 deletions Core/src/Seeding/EstimateTrackParamsFromSeed.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"

#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Utilities/MathHelpers.hpp"

#include <numbers>

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,10 @@
#include "ActsExamples/EventData/Track.hpp"
#include "ActsExamples/Framework/AlgorithmContext.hpp"

#include <cmath>
#include <cstddef>
#include <optional>
#include <ostream>
#include <stdexcept>
#include <system_error>
#include <utility>
#include <vector>

Expand Down Expand Up @@ -119,16 +117,14 @@ ProcessCode TrackParamsEstimationAlgorithm::execute(
}

// Estimate the track parameters from seed
auto optParams = Acts::estimateTrackParamsFromSeed(
ctx.geoContext, seed.sp().begin(), seed.sp().end(), *surface, field,
logger());
if (!optParams.has_value()) {
ACTS_WARNING("Estimation of track parameters for seed " << iseed
<< " failed.");
const auto paramsResult = Acts::estimateTrackParamsFromSeed(
ctx.geoContext, seed.sp(), *surface, field);
if (!paramsResult.ok()) {
ACTS_WARNING("Skip track because param estimation failed "
<< paramsResult.error());
continue;
}

const auto& params = optParams.value();
const auto& params = *paramsResult;

Acts::EstimateTrackParamCovarianceConfig config{
.initialSigmas =
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -167,17 +167,17 @@ ProcessCode PrototracksToParameters::execute(
continue;
}

auto pars = Acts::estimateTrackParamsFromSeed(
ctx.geoContext, seed.sp().begin(), seed.sp().end(), surface, field);

if (not pars) {
auto parsResult = Acts::estimateTrackParamsFromSeed(
ctx.geoContext, seed.sp(), surface, field);
if (!parsResult.ok()) {
ACTS_WARNING("Skip track because of bad params");
}
const auto &pars = *parsResult;

seededTracks.push_back(track);
seeds.emplace_back(std::move(seed));
parameters.push_back(Acts::BoundTrackParameters(
surface.getSharedPtr(), *pars, m_covariance, m_cfg.particleHypothesis));
surface.getSharedPtr(), pars, m_covariance, m_cfg.particleHypothesis));
}

if (skippedTracks > 0) {
Expand Down
8 changes: 4 additions & 4 deletions Examples/Python/tests/root_file_hashes.txt
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,11 @@ test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 8e0116c
test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e
test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: 82a6744980553e6274df78eea15f0dec22676b1c04e14afc3828bff9bbf5e1b1
test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 06d6ae1d05cb611b19df3c59531997c9b0108f5ef6027d76c4827bd2d9edb921
test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 17c48c5a61b1a5495d91336cdf06f9c24e50d81349c1f31d7c70ffff5810a376
test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: b5805e54030ab8ac80a8c0a764700c65433dc659783fc8ff3b2c96e512a1d045
test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 53a96c225ed78bd7da1491ae3ff4173f0e6c8bffe27e3eab1ab7adfe55426080
test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 2a98d8ec8fae97e18f4661580b71885f558d7222f94bca5edfbf5cdb595021f7
test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f
test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 86be5a086d2a87dfde9320bb880bd0788d733ea9727cb5ee6dc0282ec4be39f4
test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: ffce6a73f16986cb3f0386d4a8c1e0ff6f0b4130b9bb12d1af0eb905d000e3e9
test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 7ec5f4b24fa69bfc9d812a6f224d74d36df2d0a42aef791ae0116f309d6f1433
test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: 1eaae038ced2cc5c757480ca42eab60cdaff14d812c34a807a841267d6bfa110
test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575
test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 35a65e15a6f479f628a96f56ee78e1ac371d71a686ee0c974944d681499fe6bd
test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 3e257de624674fa9a19dcc72598c78c29a52633821acaa56dc2aa39a1395f1b5
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -174,11 +174,10 @@ BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) {
BOOST_CHECK(!estFreeParams.hasNaN());

// Test the bound track parameters estimator
auto fullParamsOpt = estimateTrackParamsFromSeed(
geoCtx, spacePointPtrs.begin(), spacePointPtrs.end(),
*bottomSurface, bField, *logger);
BOOST_REQUIRE(fullParamsOpt.has_value());
const auto& estFullParams = fullParamsOpt.value();
auto estFullParamsResult = estimateTrackParamsFromSeed(
geoCtx, spacePointPtrs, *bottomSurface, bField);
BOOST_CHECK(estFullParamsResult.ok());
const auto& estFullParams = estFullParamsResult.value();
BOOST_TEST_INFO(
"The estimated full track parameters at the bottom space point: "
"\n"
Expand Down
2 changes: 0 additions & 2 deletions docs/known-warnings.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,6 @@
.*Duplicate C\+\+ declaration.*
.*Duplicate ID: .*
.*Duplicate explicit target name: .*
.*doxygenfile: Found multiple matches for file \"SeedFinder.hpp.*
.*Duplicate explicit target name: .*class_acts_1_1_axis_aligned_bounding_box_.*
.*undefined label: .*class_acts_1_1_convex_polygon_bounds_3_01_polygon_dynamic_01_4.*
.*undefined label: .*class_acts_1_1_grid_surface_material_t.*
# I think these are because we use specialization
Expand Down

0 comments on commit c62c22c

Please sign in to comment.