diff --git a/Core/include/Acts/Propagator/Propagator.hpp b/Core/include/Acts/Propagator/Propagator.hpp index 6a38e5d2a8f..d3664985fbc 100644 --- a/Core/include/Acts/Propagator/Propagator.hpp +++ b/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 - namespace Acts { /// Common simplified base interface for propagators. diff --git a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp index e4618296b50..4790c3483be 100644 --- a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp +++ b/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 estimateTrackParamsFromSeedAtSurface( + const GeometryContext& gctx, const MagneticFieldContext& mctx, + const Surface& surface, const Vector3& sp0, const Vector3& sp1, + const Vector3& sp2, double timeSp0, + const std::shared_ptr& bField, + 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,52 @@ 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 -Result 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::failure(lpResult.error()); +Result estimateTrackParamsFromSeedAtSurface( + const GeometryContext& gctx, const MagneticFieldContext& mctx, + const Surface& surface, spacepoint_range_t spRange, + const std::shared_ptr& bField, + 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 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(); } - 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::success(params); + Result paramsResult = estimateTrackParamsFromSeedAtSurface( + gctx, mctx, surface, spPositions[0], spPositions[1], spPositions[2], + spTimes[0].value_or(0), bField, propagator, logger); + return paramsResult; } /// Configuration for the estimation of the covariance matrix of the track @@ -189,3 +218,10 @@ BoundMatrix estimateTrackParamCovariance( bool hasTime); } // namespace Acts + +namespace std { +// register with STL +template <> +struct is_error_code_enum + : std::true_type {}; +} // namespace std diff --git a/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp b/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp index b071d5dd7e6..754b5acafdf 100644 --- a/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp +++ b/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 +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(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(e), c}; +} + Acts::FreeVector Acts::estimateTrackParamsFromSeed(const Vector3& sp0, const Vector3& sp1, const Vector3& sp2, @@ -129,3 +168,81 @@ Acts::BoundMatrix Acts::estimateTrackParamCovariance( return result; } + +Acts::Result Acts::estimateTrackParamsFromSeedAtSurface( + const GeometryContext& gctx, const MagneticFieldContext& mctx, + const Surface& surface, const Vector3& sp0, const Vector3& sp1, + const Vector3& sp2, double timeSp0, + const std::shared_ptr& bField, + const BasePropagator* propagator, const Acts::Logger& logger) { + MagneticFieldProvider::Cache bFieldCache = bField->makeCache(mctx); + + Result bFieldResult = bField->getField(sp0, bFieldCache); + if (!bFieldResult.ok()) { + ACTS_INFO("The magnetic field failed."); + return Result::failure(bFieldResult.error()); + } + const Vector3& bFieldVector = bFieldResult.value(); + + FreeVector freeParams = + estimateTrackParamsFromSeed(sp0, sp1, sp2, bFieldVector); + freeParams[eFreeTime] = timeSp0; + + 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."); + return Result::failure( + EstimateTrackParamsFromSeedError::SurfaceIntersectionFailed); + } + + Direction propagatorDirection = + Direction::fromScalarZeroAsPositive(surfaceIntersection.pathLength()); + + PropagatorPlainOptions propagatorOptions(gctx, mctx); + propagatorOptions.direction = propagatorDirection; + + Result 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::failure( + EstimateTrackParamsFromSeedError::PropagationFailed); + } + + return Result::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(); + + return Result::success(params); +} diff --git a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp index f2711ce1d1a..b8dbd2f99c1 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp +++ b/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()); diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/PrototracksToParameters.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/PrototracksToParameters.cpp index 53fe1cdbdad..4da4c6e829b 100644 --- a/Examples/Algorithms/TrackFindingExaTrkX/src/PrototracksToParameters.cpp +++ b/Examples/Algorithms/TrackFindingExaTrkX/src/PrototracksToParameters.cpp @@ -55,7 +55,6 @@ PrototracksToParameters::~PrototracksToParameters() {} ProcessCode PrototracksToParameters::execute( const AlgorithmContext &ctx) const { - auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext); const auto &sps = m_inputSpacePoints(ctx); auto prototracks = m_inputProtoTracks(ctx); @@ -154,21 +153,9 @@ ProcessCode PrototracksToParameters::execute( .geometryId(); const auto &surface = *m_cfg.geometry->findSurface(geoId); - 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 is too small " << field.norm()); - continue; - } - - auto parsResult = Acts::estimateTrackParamsFromSeed( - ctx.geoContext, seed.sp(), surface, field); + auto parsResult = Acts::estimateTrackParamsFromSeedAtSurface( + ctx.geoContext, ctx.magFieldContext, surface, seed.sp(), + m_cfg.magneticField); if (!parsResult.ok()) { ACTS_WARNING("Skip track because of bad params"); } diff --git a/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp b/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp index 45dbbd9d0bb..6b9811d3bba 100644 --- a/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp +++ b/Tests/UnitTests/Core/Seeding/EstimateTrackParamsFromSeedTest.cpp @@ -97,7 +97,7 @@ BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) { }); const Vector3 bField(0, 0, 2._T); auto field = std::make_shared(bField); - ConstantFieldStepper stepper(std::move(field)); + ConstantFieldStepper stepper(field); ConstantFieldPropagator propagator(std::move(stepper), std::move(navigator)); @@ -174,8 +174,8 @@ BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) { BOOST_CHECK(!estFreeParams.hasNaN()); // Test the bound track parameters estimator - auto estFullParamsResult = estimateTrackParamsFromSeed( - geoCtx, spacePointPtrs, *bottomSurface, bField); + auto estFullParamsResult = estimateTrackParamsFromSeedAtSurface( + geoCtx, magCtx, *bottomSurface, spacePointPtrs, field); BOOST_CHECK(estFullParamsResult.ok()); const auto& estFullParams = estFullParamsResult.value(); BOOST_TEST_INFO(