feat: Free parameter estimation from seed

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 <array>
 #include <cmath>
 #include <iostream>
 #include <iterator>
 #include <optional>
+#include <stdexcept>
 
 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 <std::ranges::range spacepoint_range_t>
+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<Vector3, 3> spPositions = {Vector3::Zero(), Vector3::Zero(),
+                                        Vector3::Zero()};
+  std::array<std::optional<double>, 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

Core/src/Seeding/EstimateTrackParamsFromSeed.cpp

@@ -12,6 +12,84 @@
 
 #include <numbers>
 
+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) {

Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp

@@ -9,10 +9,8 @@
 #include <boost/test/unit_test.hpp>
 
 #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 <algorithm>
@@ -39,7 +36,6 @@
 #include <map>
 #include <memory>
 #include <optional>
-#include <ostream>
 #include <random>
 #include <utility>
 #include <vector>
@@ -99,8 +95,8 @@ BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) {
       true,  // material
       false  // passive
   });
-  auto field =
-      std::make_shared<Acts::ConstantBField>(Acts::Vector3(0.0, 0.0, 2._T));
+  const Vector3 bField(0, 0, 2._T);
+  auto field = std::make_shared<Acts::ConstantBField>(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(