Skip to content

Commit

Permalink
refactor: Add nonlinear correction in KF (#1233)
Browse files Browse the repository at this point in the history
This PR refactors a few methods of the stepper related with transforming from free to bound parameters. It adds the option for correction of non-linear effects during the transformation using a few sigma points (https://doi.org/10.1117/12.280797).
  • Loading branch information
XiaocongAi authored May 9, 2022
1 parent 0b8796e commit 538fd0b
Show file tree
Hide file tree
Showing 27 changed files with 647 additions and 96 deletions.
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
// This file is part of the Acts project.
//
// Copyright (C) 2021-2022 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Definitions/Common.hpp"
#include "Acts/EventData/detail/TransformationFreeToBound.hpp"
#include "Acts/Utilities/Logger.hpp"

namespace Acts {

/// @brief Free to bound transformation Correction configuration class
///
struct FreeToBoundCorrection {
/// Apply correction or not
bool apply = false;

/// UKF tuning parameters
ActsScalar alpha = 0.1;
ActsScalar beta = 2;

/// The cutoff of incident angles cosine for correction
ActsScalar cosIncidentAngleMinCutoff = 1e-5;
ActsScalar cosIncidentAngleMaxCutoff = 0.99500417;

/// Default constructor
FreeToBoundCorrection() = default;

/// Construct from boolean and UKF parameters (alpha, beta)
///
/// @param apply_ Wheter to apply correction
/// @param alpha_ The UKF tuning parameter alpha
/// @param beta_ The UKF tuning parameter beta
FreeToBoundCorrection(bool apply_, ActsScalar alpha_, ActsScalar beta_);

/// Construct from boolean only
///
/// @param apply_ Wheter to apply correction
explicit FreeToBoundCorrection(bool apply_);

/// Return boolean for applying correction or not
operator bool() const;
};

namespace detail {

/// @brief Corrected free to bound transform class based on covariance matrix sqrt root in UKF: https://doi.org/10.1117/12.280797
///
struct CorrectedFreeToBoundTransformer {
/// Construct from boolean, UKF tuning parameters (alpha, beta) and incident
/// angle cutoff for correction
///
/// @param alpha The UKF tuning parameter alpha
/// @param beta The UKF tuning parameter beta
/// @param cosIncidentAngleMinCutoff The cosine of max incident angle
/// @param cosIncidentAngleMaxCutoff The cosine of min incident angle
CorrectedFreeToBoundTransformer(ActsScalar alpha, ActsScalar beta,
ActsScalar cosIncidentAngleMinCutoff,
ActsScalar cosIncidentAngleMaxCutoff);

/// Construct from a FreeToBoundCorrection
///
/// @param freeToBoundCorrection The freeToBoundCorrection object
CorrectedFreeToBoundTransformer(
const FreeToBoundCorrection& freeToBoundCorrection);

/// Default constructors
CorrectedFreeToBoundTransformer() = default;
CorrectedFreeToBoundTransformer(const CorrectedFreeToBoundTransformer&) =
default;
CorrectedFreeToBoundTransformer(CorrectedFreeToBoundTransformer&&) = default;
CorrectedFreeToBoundTransformer& operator=(
const CorrectedFreeToBoundTransformer&) = default;
CorrectedFreeToBoundTransformer& operator=(
CorrectedFreeToBoundTransformer&&) = default;

/// Get the non-linearity corrected bound parameters and its covariance
///
/// @param freeParams The free parameters vector
/// @param freeCovariance The free parameters covariance
/// @param Surface The surface of the bound parameters being represented
/// @param geoContext The geometry context
/// @param navDir The navigation direction
/// @param logger The logger
std::optional<std::tuple<BoundVector, BoundSymMatrix>> operator()(
const FreeVector& freeParams, const FreeSymMatrix& freeCovariance,
const Surface& surface, const GeometryContext& geoContext,
NavigationDirection navDir = NavigationDirection::Forward,
LoggerWrapper logger = getDummyLogger()) const;

private:
/// The parameters to tune the weight in UKF (0 < alpha <=1)
ActsScalar m_alpha = 0.1;
ActsScalar m_beta = 2;

/// The maximum incident angle (i.e. minimum cos incident angle) cutoff for
/// correction
ActsScalar m_cosIncidentAngleMinCutoff = 1e-5;

/// The minimum incident angle (i.e. maximum cos incident angle) cutoff for
/// correction, note cos(0.1) = 0.99500417
ActsScalar m_cosIncidentAngleMaxCutoff = 0.99500417;
};

} // namespace detail
} // namespace Acts
15 changes: 11 additions & 4 deletions Core/include/Acts/Propagator/AtlasStepper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
#include "Acts/EventData/detail/TransformationBoundToFree.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
Expand Down Expand Up @@ -580,13 +581,16 @@ class AtlasStepper {
/// @param [in] state State that will be presented as @c BoundState
/// @param [in] surface The surface to which we bind the state
/// @param [in] transportCov Flag steering covariance transport
/// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
///
/// @return A bound state:
/// - the parameters at the surface
/// - the stepwise jacobian towards it
/// - and the path length (from start - for ordering)
Result<BoundState> boundState(State& state, const Surface& surface,
bool transportCov = true) const {
Result<BoundState> boundState(
State& state, const Surface& surface, bool transportCov = true,
const FreeToBoundCorrection& freeToBoundCorrection =
FreeToBoundCorrection(false)) const {
// the convert method invalidates the state (in case it's reused)
state.state_ready = false;
// extract state information
Expand All @@ -604,7 +608,7 @@ class AtlasStepper {
// The transport of the covariance
std::optional<Covariance> covOpt = std::nullopt;
if (state.covTransport && transportCov) {
transportCovarianceToBound(state, surface);
transportCovarianceToBound(state, surface, freeToBoundCorrection);
}
if (state.cov != Covariance::Zero()) {
covOpt = state.cov;
Expand Down Expand Up @@ -871,7 +875,10 @@ class AtlasStepper {
///
/// @param [in,out] state State of the stepper
/// @param [in] surface is the surface to which the covariance is forwarded to
void transportCovarianceToBound(State& state, const Surface& surface) const {
void transportCovarianceToBound(
State& state, const Surface& surface,
const FreeToBoundCorrection& /*freeToBoundCorrection*/ =
FreeToBoundCorrection(false)) const {
Acts::Vector3 gp(state.pVector[0], state.pVector[1], state.pVector[2]);
Acts::Vector3 mom(state.pVector[4], state.pVector[5], state.pVector[6]);
mom /= std::abs(state.pVector[7]);
Expand Down
13 changes: 10 additions & 3 deletions Core/include/Acts/Propagator/EigenStepper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -313,13 +313,16 @@ class EigenStepper {
/// @param [in] state State that will be presented as @c BoundState
/// @param [in] surface The surface to which we bind the state
/// @param [in] transportCov Flag steering covariance transport
/// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
///
/// @return A bound state:
/// - the parameters at the surface
/// - the stepwise jacobian towards it (from last bound)
/// - and the path length (from start - for ordering)
Result<BoundState> boundState(State& state, const Surface& surface,
bool transportCov = true) const;
Result<BoundState> boundState(
State& state, const Surface& surface, bool transportCov = true,
const FreeToBoundCorrection& freeToBoundCorrection =
FreeToBoundCorrection(false)) const;

/// Create and return a curvilinear state at the current position
///
Expand Down Expand Up @@ -372,8 +375,12 @@ class EigenStepper {
///
/// @param [in,out] state State of the stepper
/// @param [in] surface is the surface to which the covariance is forwarded to
/// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
/// @note no check is done if the position is actually on the surface
void transportCovarianceToBound(State& state, const Surface& surface) const;
void transportCovarianceToBound(
State& state, const Surface& surface,
const FreeToBoundCorrection& freeToBoundCorrection =
FreeToBoundCorrection(false)) const;

/// Perform a Runge-Kutta track parameter propagation step
///
Expand Down
14 changes: 9 additions & 5 deletions Core/include/Acts/Propagator/EigenStepper.ipp
Original file line number Diff line number Diff line change
Expand Up @@ -49,13 +49,15 @@ void Acts::EigenStepper<E, A>::resetState(State& state,
}

template <typename E, typename A>
auto Acts::EigenStepper<E, A>::boundState(State& state, const Surface& surface,
bool transportCov) const
auto Acts::EigenStepper<E, A>::boundState(
State& state, const Surface& surface, bool transportCov,
const FreeToBoundCorrection& freeToBoundCorrection) const
-> Result<BoundState> {
return detail::boundState(
state.geoContext, state.cov, state.jacobian, state.jacTransport,
state.derivative, state.jacToGlobal, state.pars,
state.covTransport && transportCov, state.pathAccumulated, surface);
state.covTransport && transportCov, state.pathAccumulated, surface,
freeToBoundCorrection);
}

template <typename E, typename A>
Expand Down Expand Up @@ -100,10 +102,12 @@ void Acts::EigenStepper<E, A>::transportCovarianceToCurvilinear(

template <typename E, typename A>
void Acts::EigenStepper<E, A>::transportCovarianceToBound(
State& state, const Surface& surface) const {
State& state, const Surface& surface,
const FreeToBoundCorrection& freeToBoundCorrection) const {
detail::transportCovarianceToBound(
state.geoContext.get(), state.cov, state.jacobian, state.jacTransport,
state.derivative, state.jacToGlobal, state.pars, surface);
state.derivative, state.jacToGlobal, state.pars, surface,
freeToBoundCorrection);
}

template <typename E, typename A>
Expand Down
29 changes: 20 additions & 9 deletions Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include "Acts/Definitions/Units.hpp"
#include "Acts/EventData/MultiComponentBoundTrackParameters.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
#include "Acts/MagneticField/MagneticFieldProvider.hpp"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/EigenStepperError.hpp"
Expand Down Expand Up @@ -430,11 +431,13 @@ class MultiEigenStepperLoop
cmp.state, state.navigation, state.options, state.geoContext);
}

Result<BoundState> boundState(const Surface& surface, bool transportCov) {
return detail::boundState(all_state.geoContext, cov(), jacobian(),
jacTransport(), derivative(), jacToGlobal(),
pars(), all_state.covTransport && transportCov,
cmp.state.pathAccumulated, surface);
Result<BoundState> boundState(
const Surface& surface, bool transportCov,
const FreeToBoundCorrection& freeToBoundCorrection) {
return detail::boundState(
all_state.geoContext, cov(), jacobian(), jacTransport(), derivative(),
jacToGlobal(), pars(), all_state.covTransport && transportCov,
cmp.state.pathAccumulated, surface, freeToBoundCorrection);
}

void update(const FreeVector& freeParams, const BoundVector& boundParams,
Expand Down Expand Up @@ -761,13 +764,16 @@ class MultiEigenStepperLoop
/// @param [in] state State that will be presented as @c BoundState
/// @param [in] surface The surface to which we bind the state
/// @param [in] transportCov Flag steering covariance transport
/// @param [in] freeToBoundCorrection Flag steering non-linear correction during global to local correction
///
/// @return A bound state:
/// - the parameters at the surface
/// - the stepwise jacobian towards it (from last bound)
/// - and the path length (from start - for ordering)
Result<BoundState> boundState(State& state, const Surface& surface,
bool transportCov = true) const;
Result<BoundState> boundState(
State& state, const Surface& surface, bool transportCov = true,
const FreeToBoundCorrection& freeToBoundCorrection =
FreeToBoundCorrection(false)) const;

/// Create and return a curvilinear state at the current position
///
Expand Down Expand Up @@ -805,11 +811,16 @@ class MultiEigenStepperLoop
///
/// @param [in,out] state State of the stepper
/// @param [in] surface is the surface to which the covariance is forwarded
/// @param [in] freeToBoundCorrection Flag steering non-linear correction during global to local correction
/// to
/// @note no check is done if the position is actually on the surface
void transportCovarianceToBound(State& state, const Surface& surface) const {
void transportCovarianceToBound(
State& state, const Surface& surface,
const FreeToBoundCorrection& freeToBoundCorrection =
FreeToBoundCorrection(false)) const {
for (auto& component : state.components) {
SingleStepper::transportCovarianceToBound(component.state, surface);
SingleStepper::transportCovarianceToBound(component.state, surface,
freeToBoundCorrection);
}
}

Expand Down
12 changes: 6 additions & 6 deletions Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp
Original file line number Diff line number Diff line change
Expand Up @@ -9,21 +9,21 @@
namespace Acts {

template <typename E, typename R, typename A>
auto MultiEigenStepperLoop<E, R, A>::boundState(State& state,
const Surface& surface,
bool transportCov) const
auto MultiEigenStepperLoop<E, R, A>::boundState(
State& state, const Surface& surface, bool transportCov,
const FreeToBoundCorrection& freeToBoundCorrection) const
-> Result<BoundState> {
if (numberComponents(state) == 1) {
return SingleStepper::boundState(state.components.front().state, surface,
transportCov);
} else { // Do the combination
transportCov, freeToBoundCorrection);
} else { // Do the combinatio
SmallVector<std::pair<double, BoundTrackParameters>> states;
double accumulatedPathLength = 0.0;
int failedBoundTransforms = 0;

for (auto i = 0ul; i < numberComponents(state); ++i) {
auto bs = SingleStepper::boundState(state.components[i].state, surface,
transportCov);
transportCov, freeToBoundCorrection);

if (bs.ok()) {
states.push_back(
Expand Down
5 changes: 3 additions & 2 deletions Core/include/Acts/Propagator/StepperConcept.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
#include "Acts/Propagator/ConstrainedStep.hpp"
#include "Acts/Surfaces/BoundaryCheck.hpp"
#include "Acts/Surfaces/Surface.hpp"
Expand Down Expand Up @@ -114,12 +115,12 @@ constexpr bool MultiStepperStateConcept= require<
static_assert(time_exists, "time method not found");
constexpr static bool overstep_exists = has_method<const S, double, overstep_t, const state&>;
static_assert(overstep_exists, "overstepLimit method not found");
constexpr static bool bound_state_method_exists= has_method<const S, Result<typename S::BoundState>, bound_state_method_t, state&, const Surface&, bool>;
constexpr static bool bound_state_method_exists= has_method<const S, Result<typename S::BoundState>, bound_state_method_t, state&, const Surface&, bool, const FreeToBoundCorrection&>;
static_assert(bound_state_method_exists, "boundState method not found");
constexpr static bool curvilinear_state_method_exists = has_method<const S, typename S::CurvilinearState, curvilinear_state_method_t, state&, bool>;
static_assert(curvilinear_state_method_exists, "curvilinearState method not found");
constexpr static bool covariance_transport_exists = require<has_method<const S, void, covariance_transport_curvilinear_t, state&>,
has_method<const S, void, covariance_transport_bound_t, state&, const Surface&>>;
has_method<const S, void, covariance_transport_bound_t, state&, const Surface&, const FreeToBoundCorrection&>>;
static_assert(covariance_transport_exists, "covarianceTransport method not found");
constexpr static bool update_surface_exists = has_method<const S, Intersection3D::Status, update_surface_status_t, state&, const Surface&, const BoundaryCheck&, LoggerWrapper>;
static_assert(update_surface_exists, "updateSurfaceStatus method not found");
Expand Down
14 changes: 11 additions & 3 deletions Core/include/Acts/Propagator/StraightLineStepper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/MagneticField/NullBField.hpp"
Expand Down Expand Up @@ -283,13 +284,16 @@ class StraightLineStepper {
/// @param [in] state State that will be presented as @c BoundState
/// @param [in] surface The surface to which we bind the state
/// @param [in] transportCov Flag steering covariance transport
/// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
///
/// @return A bound state:
/// - the parameters at the surface
/// - the stepwise jacobian towards it (from last bound)
/// - and the path length (from start - for ordering)
Result<BoundState> boundState(State& state, const Surface& surface,
bool transportCov = true) const;
Result<BoundState> boundState(
State& state, const Surface& surface, bool transportCov = true,
const FreeToBoundCorrection& freeToBoundCorrection =
FreeToBoundCorrection(false)) const;

/// Create and return a curvilinear state at the current position
///
Expand Down Expand Up @@ -343,8 +347,12 @@ class StraightLineStepper {
/// @param [in] surface is the surface to which the covariance is
/// forwarded to
/// @note no check is done if the position is actually on the surface
/// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
///
void transportCovarianceToBound(State& state, const Surface& surface) const;
void transportCovarianceToBound(
State& state, const Surface& surface,
const FreeToBoundCorrection& freeToBoundCorrection =
FreeToBoundCorrection(false)) const;

/// Perform a straight line propagation step
///
Expand Down
Loading

0 comments on commit 538fd0b

Please sign in to comment.