diff --git a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp index 6d6a62516a7..5d6b4dbb97c 100644 --- a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp +++ b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp @@ -14,15 +14,95 @@ #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/MathHelpers.hpp" +#include "Acts/Utilities/Zip.hpp" #include #include #include #include #include +#include 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. +/// +/// 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: +/// 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 sp0 is the bottom space point +/// @param sp1 is the middle space point +/// @param sp2 is the top space point +/// @param bField is the magnetic field vector +/// +/// @return the free parameters +FreeVector estimateTrackParamsFromSeed(const Vector3& sp0, const Vector3& sp1, + const Vector3& sp2, + const Vector3& bField); + +/// 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. +/// +/// 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 spRange is the range of space points +/// @param bField is the magnetic field vector +/// +/// @return the free parameters +template +FreeVector estimateTrackParamsFromSeed(spacepoint_range_t spRange, + const Vector3& bField) { + // Check the number of provided space points + if (spRange.size() != 3) { + throw std::invalid_argument( + "There should be exactly three space points provided."); + } + + // The global positions of the bottom, middle and space points + std::array spPositions = {Vector3::Zero(), Vector3::Zero(), + Vector3::Zero()}; + std::array, 3> spTimes = {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 (auto [sp, spPosition, spTime] : + Acts::zip(spRange, spPositions, spTimes)) { + if (sp == nullptr) { + throw std::invalid_argument("Empty space point found."); + } + spPosition = Vector3(sp->x(), sp->y(), sp->z()); + spTime = sp->t(); + } + + FreeVector params = estimateTrackParamsFromSeed( + spPositions[0], spPositions[1], spPositions[2], bField); + params[eFreeTime] = spTimes[0].value_or(0); + return params; +} + /// Estimate the full track parameters from three space points /// /// This method is based on the conformal map transformation. It estimates the diff --git a/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp b/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp index eebcd1c5760..055691bbb2f 100644 --- a/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp +++ b/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp @@ -12,6 +12,84 @@ #include +Acts::FreeVector Acts::estimateTrackParamsFromSeed(const Vector3& sp0, + const Vector3& sp1, + const Vector3& sp2, + const Vector3& bField) { + // 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 = sp1 - sp0; + 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(sp0); + // 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() * sp1; + Vector3 local2 = transform.inverse() * sp2; + + // 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); + + 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(); + + // Initialize the free parameters vector + FreeVector params = FreeVector::Zero(); + + // The bottom space point position + params.segment<3>(eFreePos0) = sp0; + + // The estimated direction + params.segment<3>(eFreeDir0) = direction; + + // 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 / (bField.norm() * R); + // The estimated q/p in [GeV/c]^-1 + params[eFreeQOverP] = qOverPt / fastHypot(1., invTanTheta); + + return params; +} + Acts::BoundMatrix Acts::estimateTrackParamCovariance( const EstimateTrackParamCovarianceConfig& config, const BoundVector& params, bool hasTime) { diff --git a/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp b/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp index 6161a8bdd33..e39a6a1a900 100644 --- a/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp +++ b/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp @@ -9,10 +9,8 @@ #include #include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/Direction.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/GenericCurvilinearTrackParameters.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/EventData/detail/TestSourceLink.hpp" #include "Acts/Geometry/GeometryContext.hpp" @@ -29,7 +27,6 @@ #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp" #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp" -#include "Acts/Utilities/CalibrationContext.hpp" #include "Acts/Utilities/Logger.hpp" #include @@ -39,7 +36,6 @@ #include #include #include -#include #include #include #include @@ -99,8 +95,8 @@ BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) { true, // material false // passive }); - auto field = - std::make_shared(Acts::Vector3(0.0, 0.0, 2._T)); + const Vector3 bField(0, 0, 2._T); + auto field = std::make_shared(bField); ConstantFieldStepper stepper(std::move(field)); ConstantFieldPropagator propagator(std::move(stepper), std::move(navigator)); @@ -172,10 +168,15 @@ BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) { spacePointPtrs.begin(), [](const auto& sp) { return &sp.second; }); - // Test the full track parameters estimator + // Test the free track parameters estimator + FreeVector estFreeParams = + estimateTrackParamsFromSeed(spacePointPtrs, bField); + BOOST_CHECK(!estFreeParams.hasNaN()); + + // Test the bound track parameters estimator auto fullParamsOpt = estimateTrackParamsFromSeed( geoCtx, spacePointPtrs.begin(), spacePointPtrs.end(), - *bottomSurface, Vector3(0, 0, 2._T), 0.1_T, *logger); + *bottomSurface, bField, 0.1_T, *logger); BOOST_REQUIRE(fullParamsOpt.has_value()); const auto& estFullParams = fullParamsOpt.value(); BOOST_TEST_INFO(