Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor!: Remove mean reduction from MultiEigenStepperLoop #3671

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
100 changes: 1 addition & 99 deletions Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -45,104 +45,6 @@ namespace Acts {

using namespace Acts::UnitLiterals;

/// @brief Reducer struct for the Loop MultiEigenStepper which reduces the
/// multicomponent state to simply by summing the weighted values
///
/// @note Usage is not encouraged, since it can lead to navigation failures
/// as the global position might not be on surface, even if all components
/// are on surface
struct WeightedComponentReducerLoop {
template <typename component_range_t>
static Vector3 toVector3(const component_range_t& comps,
const FreeIndices i) {
return std::accumulate(
comps.begin(), comps.end(), Vector3{Vector3::Zero()},
[i](const auto& sum, const auto& cmp) -> Vector3 {
return sum + cmp.weight * cmp.state.pars.template segment<3>(i);
});
}

template <typename stepper_state_t>
static Vector3 position(const stepper_state_t& s) {
return toVector3(s.components, eFreePos0);
}

template <typename stepper_state_t>
static Vector3 direction(const stepper_state_t& s) {
return toVector3(s.components, eFreeDir0).normalized();
}

// TODO: Maybe we can cache this value and only update it when the parameters
// change
template <typename stepper_state_t>
static ActsScalar qOverP(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * cmp.state.pars[eFreeQOverP];
});
}

template <typename stepper_state_t>
static ActsScalar absoluteMomentum(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[&s](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * s.particleHypothesis.extractMomentum(
cmp.state.pars[eFreeQOverP]);
});
}

template <typename stepper_state_t>
static Vector3 momentum(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), Vector3::Zero().eval(),
[&s](const auto& sum, const auto& cmp) -> Vector3 {
return sum + cmp.weight *
s.particleHypothesis.extractMomentum(
cmp.state.pars[eFreeQOverP]) *
cmp.state.pars.template segment<3>(eFreeDir0);
});
}

template <typename stepper_state_t>
static ActsScalar charge(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[&s](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * s.particleHypothesis.extractCharge(
cmp.state.pars[eFreeQOverP]);
});
}

template <typename stepper_state_t>
static ActsScalar time(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * cmp.state.pars[eFreeTime];
});
}

template <typename stepper_state_t>
static FreeVector pars(const stepper_state_t& s) {
return std::accumulate(s.components.begin(), s.components.end(),
FreeVector{FreeVector::Zero()},
[](const auto& sum, const auto& cmp) -> FreeVector {
return sum + cmp.weight * cmp.state.pars;
});
}

template <typename stepper_state_t>
static FreeVector cov(const stepper_state_t& s) {
return std::accumulate(s.components.begin(), s.components.end(),
FreeMatrix{FreeMatrix::Zero()},
[](const auto& sum, const auto& cmp) -> FreeMatrix {
return sum + cmp.weight * cmp.state.cov;
});
}
};

namespace detail {

struct MaxMomentumComponent {
Expand Down Expand Up @@ -240,7 +142,7 @@ using MaxWeightReducerLoop =
/// @tparam small_vector_size A size-hint how much memory should be allocated
/// by the small vector
template <typename extension_t = EigenStepperDefaultExtension,
typename component_reducer_t = WeightedComponentReducerLoop>
typename component_reducer_t = MaxWeightReducerLoop>
class MultiEigenStepperLoop : public EigenStepper<extension_t> {
/// Limits the number of steps after at least one component reached the
/// surface
Expand Down
25 changes: 0 additions & 25 deletions Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -138,31 +138,6 @@ auto makeDefaultBoundPars(bool cov = true, std::size_t n = 4,
//////////////////////
/// Test the reducers
//////////////////////
BOOST_AUTO_TEST_CASE(test_weighted_reducer) {
// Can use this multistepper since we only care about the state which is
// invariant
using MultiState = typename MultiStepperLoop::State;

constexpr std::size_t N = 4;
const auto multi_pars = makeDefaultBoundPars(false, N);

MultiState state(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize);
SingleStepper singleStepper(defaultBField);

WeightedComponentReducerLoop reducer{};

Acts::Vector3 pos = Acts::Vector3::Zero();
Acts::Vector3 dir = Acts::Vector3::Zero();
for (const auto &[sstate, weight, _] : state.components) {
pos += weight * singleStepper.position(sstate);
dir += weight * singleStepper.direction(sstate);
}
dir.normalize();

BOOST_CHECK_EQUAL(reducer.position(state), pos);
BOOST_CHECK_EQUAL(reducer.direction(state), dir);
}

BOOST_AUTO_TEST_CASE(test_max_weight_reducer) {
// Can use this multistepper since we only care about the state which is
// invariant
Expand Down
2 changes: 1 addition & 1 deletion docs/core/reconstruction/track_fitting.md
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ Simplified overview of the GSF algorithm.
:::

### The Multi-Stepper
To implement the GSF, a special stepper is needed, that can handle a multi-component state internally: The {class}`Acts::MultiEigenStepperLoop`, which is based on the {class}`Acts::EigenStepper` and thus shares a lot of code with it. It interfaces to the navigation as one aggregate state to limit the navigation overhead, but internally processes a multi-component state. How this aggregation is performed can be configured via a template parameter, by default weighted average is used ({struct}`Acts::WeightedComponentReducerLoop`).
To implement the GSF, a special stepper is needed, that can handle a multi-component state internally: The {class}`Acts::MultiEigenStepperLoop`, which is based on the {class}`Acts::EigenStepper` and thus shares a lot of code with it. It interfaces to the navigation as one aggregate state to limit the navigation overhead, but internally processes a multi-component state. How this aggregation is performed can be configured via a template parameter, by default maximum weight is used ({type}`Acts::MaxWeightReducerLoop`).

Even though the multi-stepper interface exposes only one aggregate state and thus is compatible with most standard tools, there is a special aborter is required to stop the navigation when the surface is reached, the {struct}`Acts::MultiStepperSurfaceReached`. It checks if all components have reached the target surface already and updates their state accordingly. Optionally, it also can stop the propagation when the aggregate state reaches the surface.

Expand Down
Loading