feat!: Estimate parameters on surface also if bottom SP is not on surface

Core/include/Acts/Propagator/Propagator.hpp

@@ -13,13 +13,8 @@
 #include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp"
 // clang-format on
 
-#include "Acts/Definitions/Algebra.hpp"
-#include "Acts/Definitions/PdgParticle.hpp"
-#include "Acts/Definitions/Units.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/EventData/TrackParametersConcept.hpp"
-#include "Acts/Geometry/GeometryContext.hpp"
-#include "Acts/MagneticField/MagneticFieldContext.hpp"
 #include "Acts/Propagator/ActorList.hpp"
 #include "Acts/Propagator/PropagatorOptions.hpp"
 #include "Acts/Propagator/PropagatorResult.hpp"
@@ -32,8 +27,6 @@
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
 
-#include <optional>
-
 namespace Acts {
 
 /// Common simplified base interface for propagators.

Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp

@@ -11,7 +11,6 @@
 #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/Zip.hpp"
 
@@ -21,6 +20,19 @@
 
 namespace Acts {
 
+class GeometryContext;
+class MagneticFieldContext;
+class MagneticFieldProvider;
+class BasePropagator;
+
+enum class EstimateTrackParamsFromSeedError {
+  // ensure all values are non-zero
+  SurfaceIntersectionFailed = 1,
+  PropagationFailed,
+};
+
+std::error_code make_error_code(EstimateTrackParamsFromSeedError e);
+
 /// Estimate the full track parameters from three space points
 ///
 /// This method is based on the conformal map transformation. It estimates the
@@ -50,15 +62,36 @@ 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.
+/// 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
+/// See the free parameter version for more details.
+///
+/// @tparam spacepoint_iterator_t  The type of space point iterator
+///
+/// @param gctx the geometry context
+/// @param mctx the magnetic field context
+/// @param surface the surface where the bound parameters should be expressed
+/// @param sp0 the bottom space point
+/// @param sp1 the middle space point
+/// @param sp2 the top space point
+/// @param bField the magnetic field provider
+/// @param propagator is the propagator used in case the bottom spacepoint is not on the surface
+/// @param logger the logger to be used
+///
+/// @return bound parameters result
+Result<BoundVector> estimateTrackParamsFromSeedAtSurface(
+    const GeometryContext& gctx, const MagneticFieldContext& mctx,
+    const Surface& surface, const Vector3& sp0, const Vector3& sp1,
+    const Vector3& sp2,
+    const std::shared_ptr<const MagneticFieldProvider>& bField,
+    const std::shared_ptr<const BasePropagator>& propagator = nullptr,
+    const Acts::Logger& logger = getDummyLogger());
+
+/// Estimate the full track parameters from three space points
+///
+/// See the free parameter version with the explicit spacepoints for more
+/// details.
 ///
 /// @tparam spacepoint_iterator_t  The type of space point iterator
 ///
@@ -99,56 +132,55 @@ FreeVector estimateTrackParamsFromSeed(spacepoint_range_t spRange,
 
 /// Estimate the full track parameters from three space points
 ///
-/// 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.
-///
-/// 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
+/// See the free parameter version for more details.
 ///
 /// @tparam spacepoint_iterator_t  The type of space point iterator
 ///
-/// @param gctx is the geometry context
-/// @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 gctx the geometry context
+/// @param mctx the magnetic field context
+/// @param surface the surface where the bound parameters should be expressed
+/// @param spRange the range of space points (bottom, middle, top)
+/// @param bField the magnetic field provider
+/// @param propagator is the propagator used in case the bottom spacepoint is not on the surface
+/// @param logger the logger to be used
 ///
-/// @return bound parameters
+/// @return bound parameters result
 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);
-
-  const auto* sp0 = *spRange.begin();
-  Vector3 origin = Vector3(sp0->x(), sp0->y(), sp0->z());
-  Vector3 direction = freeParams.segment<3>(eFreeDir0);
-
-  BoundVector params = BoundVector::Zero();
-  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, origin, direction);
-  if (!lpResult.ok()) {
-    return Result<BoundVector>::failure(lpResult.error());
+Result<BoundVector> estimateTrackParamsFromSeedAtSurface(
+    const GeometryContext& gctx, const MagneticFieldContext& mctx,
+    const Surface& surface, spacepoint_range_t spRange,
+    const std::shared_ptr<const MagneticFieldProvider>& bField,
+    const std::shared_ptr<const BasePropagator>& propagator = nullptr,
+    const Acts::Logger& logger = getDummyLogger()) {
+  // 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();
   }
-  Vector2 bottomLocalPos = lpResult.value();
-  // The estimated loc0 and loc1
-  params[eBoundLoc0] = bottomLocalPos.x();
-  params[eBoundLoc1] = bottomLocalPos.y();
-  params[eBoundTime] = sp0->t().value_or(0);
 
-  return Result<BoundVector>::success(params);
+  Result<BoundVector> paramsResult = estimateTrackParamsFromSeedAtSurface(
+      gctx, mctx, surface, spPositions[0], spPositions[1], spPositions[2],
+      bField, propagator, logger);
+  if (paramsResult.ok()) {
+    paramsResult.value()[eBoundTime] = spTimes[0].value_or(0);
+  }
+  return paramsResult;
 }
 
 /// Configuration for the estimation of the covariance matrix of the track
@@ -189,3 +221,10 @@ BoundMatrix estimateTrackParamCovariance(
     bool hasTime);
 
 }  // namespace Acts
+
+namespace std {
+// register with STL
+template <>
+struct is_error_code_enum<Acts::EstimateTrackParamsFromSeedError>
+    : std::true_type {};
+}  // namespace std

Core/src/Seeding/EstimateTrackParamsFromSeed.cpp

@@ -8,11 +8,50 @@
 
 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
 
+#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Propagator/PropagatorOptions.hpp"
+#include "Acts/Propagator/VoidNavigator.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
 
 #include <numbers>
 
+namespace {
+
+class EstimateTrackParamsFromSeedErrorCategory : public std::error_category {
+ public:
+  // Return a short descriptive name for the category.
+  const char* name() const noexcept final {
+    return "EstimateTrackParamsFromSeedError";
+  }
+
+  // Return what each enum means in text.
+  std::string message(int c) const final {
+    using Acts::EstimateTrackParamsFromSeedError;
+
+    switch (static_cast<EstimateTrackParamsFromSeedError>(c)) {
+      case EstimateTrackParamsFromSeedError::SurfaceIntersectionFailed:
+        return "Surface intersection failed";
+      case EstimateTrackParamsFromSeedError::PropagationFailed:
+        return "Propagation failed";
+      default:
+        return "unknown";
+    }
+  }
+};
+
+}  // namespace
+
+std::error_code Acts::make_error_code(
+    Acts::EstimateTrackParamsFromSeedError e) {
+  static EstimateTrackParamsFromSeedErrorCategory c;
+  return {static_cast<int>(e), c};
+}
+
 Acts::FreeVector Acts::estimateTrackParamsFromSeed(const Vector3& sp0,
                                                    const Vector3& sp1,
                                                    const Vector3& sp2,
@@ -129,3 +168,84 @@ Acts::BoundMatrix Acts::estimateTrackParamCovariance(
 
   return result;
 }
+
+Acts::Result<Acts::BoundVector> Acts::estimateTrackParamsFromSeedAtSurface(
+    const GeometryContext& gctx, const MagneticFieldContext& mctx,
+    const Surface& surface, const Vector3& sp0, const Vector3& sp1,
+    const Vector3& sp2,
+    const std::shared_ptr<const MagneticFieldProvider>& bField,
+    const std::shared_ptr<const BasePropagator>& propagator,
+    const Acts::Logger& logger) {
+  MagneticFieldProvider::Cache bFieldCache = bField->makeCache(mctx);
+
+  Result<Vector3> bFieldResult = bField->getField(sp0, bFieldCache);
+  if (!bFieldResult.ok()) {
+    ACTS_INFO("The magnetic field failed.");
+    return Result<BoundVector>::failure(bFieldResult.error());
+  }
+  const Vector3& bFieldVector = bFieldResult.value();
+
+  FreeVector freeParams =
+      estimateTrackParamsFromSeed(sp0, sp1, sp2, bFieldVector);
+
+  Vector3 origin = sp0;
+  Vector3 direction = freeParams.segment<3>(eFreeDir0);
+  double qOverPt = freeParams[eFreeQOverP];
+
+  auto lpResult = surface.globalToLocal(gctx, origin, direction);
+  if (!lpResult.ok()) {
+    // no cov transport matrix is needed here
+    // particle hypothesis does not matter here
+    CurvilinearTrackParameters estimatedParams(freeParams.segment<4>(eFreePos0),
+                                               direction, qOverPt, std::nullopt,
+                                               ParticleHypothesis::pion());
+
+    auto surfaceIntersection =
+        surface.intersect(gctx, origin, direction).closest();
+
+    if (!surfaceIntersection.isValid()) {
+      ACTS_INFO(
+          "The surface does not intersect with the origin and estimated "
+          "direction.");
+      // TODO
+      return Result<BoundVector>::failure(
+          EstimateTrackParamsFromSeedError::SurfaceIntersectionFailed);
+    }
+
+    Direction propagatorDirection =
+        Direction::fromScalarZeroAsPositive(surfaceIntersection.pathLength());
+
+    PropagatorPlainOptions propagatorOptions(gctx, mctx);
+    propagatorOptions.direction = propagatorDirection;
+
+    Result<BoundTrackParameters> result =
+        (propagator != nullptr)
+            ? propagator->propagateToSurface(estimatedParams, surface,
+                                             propagatorOptions)
+            : Propagator(EigenStepper<>(bField), VoidNavigator(),
+                         logger().cloneWithSuffix("Propagator"))
+                  .propagateToSurface(estimatedParams, surface,
+                                      propagatorOptions);
+    if (!result.ok()) {
+      ACTS_INFO("The propagation failed.");
+      return Result<BoundVector>::failure(
+          EstimateTrackParamsFromSeedError::PropagationFailed);
+    }
+
+    return Result<BoundVector>::success(result.value().parameters());
+  }
+
+  BoundVector params = BoundVector::Zero();
+  params[eBoundPhi] = VectorHelpers::phi(direction);
+  params[eBoundTheta] = VectorHelpers::theta(direction);
+  params[eBoundQOverP] = freeParams[eFreeQOverP];
+
+  Vector2 bottomLocalPos = lpResult.value();
+  // The estimated loc0 and loc1
+  params[eBoundLoc0] = bottomLocalPos.x();
+  params[eBoundLoc1] = bottomLocalPos.y();
+  // We need the bottom space point time later
+  params[eBoundTime] = 0;
+
+  return Result<BoundVector>::success(params);
+}

Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp

@@ -8,7 +8,6 @@
 
 #include "ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp"
 
-#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/Seed.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
@@ -101,24 +100,10 @@ ProcessCode TrackParamsEstimationAlgorithm::execute(
       continue;
     }
 
-    // Get the magnetic field at the bottom space point
-    auto fieldRes = m_cfg.magneticField->getField(
-        {bottomSP->x(), bottomSP->y(), bottomSP->z()}, bCache);
-    if (!fieldRes.ok()) {
-      ACTS_ERROR("Field lookup error: " << fieldRes.error());
-      return ProcessCode::ABORT;
-    }
-    Acts::Vector3 field = *fieldRes;
-
-    if (field.norm() < m_cfg.bFieldMin) {
-      ACTS_WARNING("Magnetic field at seed " << iseed << " is too small "
-                                             << field.norm());
-      continue;
-    }
-
     // Estimate the track parameters from seed
-    const auto paramsResult = Acts::estimateTrackParamsFromSeed(
-        ctx.geoContext, seed.sp(), *surface, field);
+    const auto paramsResult = Acts::estimateTrackParamsFromSeedAtSurface(
+        ctx.geoContext, ctx.magFieldContext, *surface, seed.sp(),
+        m_cfg.magneticField);
     if (!paramsResult.ok()) {
       ACTS_WARNING("Skip track because param estimation failed "
                    << paramsResult.error());