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.

Sign up for GitHub

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: Disentangle stepper state creation and initialization #4069

Merged
merged 3 commits into from
Feb 5, 2025
Merged

refactor: Disentangle stepper state creation and initialization #4069

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5a11b34
refactor!: Disentangle stepper state creation and initialization
andiwand Feb 4, 2025
b9f1e88
try workaround compiler error
andiwand Feb 4, 2025
bd7a592
try again
andiwand Feb 4, 2025
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
77 changes: 36 additions & 41 deletions Core/include/Acts/Propagator/AtlasStepper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@
#include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp"

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Common.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/EventData/TransformationHelpers.hpp"
Expand Down Expand Up @@ -137,33 +139,53 @@ class AtlasStepper {

explicit AtlasStepper(const Config& config) : m_bField(config.bField) {}

State makeState(const Options& options,
const BoundTrackParameters& par) const {
State makeState(const Options& options) const {
State state{options, m_bField->makeCache(options.magFieldContext)};
return state;
}

void initialize(State& state, const BoundTrackParameters& par) const {
initialize(state, par.parameters(), par.covariance(),
par.particleHypothesis(), par.referenceSurface());
}

void initialize(State& state, const BoundVector& boundParams,
const std::optional<BoundMatrix>& cov,
ParticleHypothesis particleHypothesis,
const Surface& surface) const {
state.particleHypothesis = particleHypothesis;

state.particleHypothesis = par.particleHypothesis();
state.pathAccumulated = 0;
state.nSteps = 0;
state.nStepTrials = 0;
state.stepSize = ConstrainedStep(state.options.maxStepSize);
state.previousStepSize = 0;
state.statistics = StepperStatistics();

// The rest of this constructor is copy&paste of AtlasStepper::update() -
// this is a nasty but working solution for the stepper state without
// functions

const auto pos = par.position(options.geoContext);
const auto Vp = par.parameters();
const auto& Vp = boundParams;

double Sf = std::sin(Vp[eBoundPhi]);
double Cf = std::cos(Vp[eBoundPhi]);
double Se = std::sin(Vp[eBoundTheta]);
double Ce = std::cos(Vp[eBoundTheta]);

const Vector3 dir = {Cf * Se, Sf * Se, Ce};
const auto pos = surface.localToGlobal(
state.options.geoContext, boundParams.segment<2>(eBoundLoc0), dir);

double* pVector = state.pVector.data();

pVector[0] = pos[ePos0];
pVector[1] = pos[ePos1];
pVector[2] = pos[ePos2];
pVector[3] = par.time();
pVector[4] = Cf * Se;
pVector[5] = Sf * Se;
pVector[6] = Ce;
pVector[3] = boundParams[eBoundTime];
pVector[4] = dir[ePos0];
pVector[5] = dir[ePos1];
pVector[6] = dir[ePos2];
pVector[7] = Vp[eBoundQOverP];

// @todo: remove magic numbers - is that the charge ?
Expand All @@ -173,13 +195,13 @@ class AtlasStepper {
}

// prepare the jacobian if we have a covariance
if (par.covariance()) {
state.covTransport = cov.has_value();
if (state.covTransport) {
// copy the covariance matrix
state.covariance = new BoundSquareMatrix(*par.covariance());
state.covTransport = true;
state.covariance = new BoundSquareMatrix(*cov);
state.useJacobian = true;
const auto transform = par.referenceSurface().referenceFrame(
options.geoContext, pos, par.direction());
const auto transform =
surface.referenceFrame(state.options.geoContext, pos, dir);

pVector[8] = transform(0, eBoundLoc0);
pVector[16] = transform(0, eBoundLoc1);
Expand Down Expand Up @@ -243,7 +265,6 @@ class AtlasStepper {
pVector[59] = 0.;

// special treatment for surface types
const auto& surface = par.referenceSurface();
// the disc needs polar coordinate adaptations
if (surface.type() == Surface::Disc) {
double lCf = std::cos(Vp[1]);
Expand Down Expand Up @@ -305,34 +326,8 @@ class AtlasStepper {
}
}

state.stepSize = ConstrainedStep(options.maxStepSize);

// now declare the state as ready
state.state_ready = true;

return state;
}

/// @brief Resets the state
///
/// @param [in, out] state State of the stepper
/// @param [in] boundParams Parameters in bound parametrisation
/// @param [in] cov Covariance matrix
/// @param [in] surface Reset state will be on this surface
/// @param [in] stepSize Step size
void resetState(
State& state, const BoundVector& boundParams,
const BoundSquareMatrix& cov, const Surface& surface,
const double stepSize = std::numeric_limits<double>::max()) const {
// Update the stepping state
update(state,
transformBoundToFreeParameters(surface, state.options.geoContext,
boundParams),
boundParams, cov, surface);
state.stepSize = ConstrainedStep(stepSize);
state.pathAccumulated = 0.;

setIdentityJacobian(state);
}

/// Get the field for the stepping
Expand Down
21 changes: 7 additions & 14 deletions Core/include/Acts/Propagator/EigenStepper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,6 @@
#include "Acts/Utilities/Intersection.hpp"
#include "Acts/Utilities/Result.hpp"

#include <limits>
#include <type_traits>

namespace Acts {
Expand Down Expand Up @@ -150,20 +149,14 @@ class EigenStepper {
/// @param [in] config The configuration of the stepper
explicit EigenStepper(const Config& config) : m_bField(config.bField) {}

State makeState(const Options& options,
const BoundTrackParameters& par) const;
State makeState(const Options& options) const;

/// @brief Resets the state
///
/// @param [in, out] state State of the stepper
/// @param [in] boundParams Parameters in bound parametrisation
/// @param [in] cov Covariance matrix
/// @param [in] surface The reference surface of the bound parameters
/// @param [in] stepSize Step size
void resetState(
State& state, const BoundVector& boundParams,
const BoundSquareMatrix& cov, const Surface& surface,
const double stepSize = std::numeric_limits<double>::max()) const;
void initialize(State& state, const BoundTrackParameters& par) const;

void initialize(State& state, const BoundVector& boundParams,
const std::optional<BoundMatrix>& cov,
ParticleHypothesis particleHypothesis,
const Surface& surface) const;

/// Get the field for the stepping, it checks first if the access is still
/// within the Cell, and updates the cell if necessary.
Expand Down
72 changes: 32 additions & 40 deletions Core/include/Acts/Propagator/EigenStepper.ipp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,57 +20,49 @@ Acts::EigenStepper<E>::EigenStepper(
: m_bField(std::move(bField)) {}

template <typename E>
auto Acts::EigenStepper<E>::makeState(
const Options& options, const BoundTrackParameters& par) const -> State {
auto Acts::EigenStepper<E>::makeState(const Options& options) const -> State {
State state{options, m_bField->makeCache(options.magFieldContext)};

state.particleHypothesis = par.particleHypothesis();

Vector3 position = par.position(options.geoContext);
Vector3 direction = par.direction();
state.pars.template segment<3>(eFreePos0) = position;
state.pars.template segment<3>(eFreeDir0) = direction;
state.pars[eFreeTime] = par.time();
state.pars[eFreeQOverP] = par.parameters()[eBoundQOverP];

// Init the jacobian matrix if needed
if (par.covariance()) {
// Get the reference surface for navigation
const auto& surface = par.referenceSurface();
// set the covariance transport flag to true and copy
state.covTransport = true;
state.cov = BoundSquareMatrix(*par.covariance());
state.jacToGlobal =
surface.boundToFreeJacobian(options.geoContext, position, direction);
}

state.stepSize = ConstrainedStep(options.maxStepSize);

return state;
}

template <typename E>
void Acts::EigenStepper<E>::resetState(State& state,
void Acts::EigenStepper<E>::initialize(State& state,
const BoundTrackParameters& par) const {
initialize(state, par.parameters(), par.covariance(),
par.particleHypothesis(), par.referenceSurface());
}

template <typename E>
void Acts::EigenStepper<E>::initialize(State& state,
const BoundVector& boundParams,
const BoundSquareMatrix& cov,
const Surface& surface,
const double stepSize) const {
const std::optional<BoundMatrix>& cov,
ParticleHypothesis particleHypothesis,
const Surface& surface) const {
FreeVector freeParams = transformBoundToFreeParameters(
surface, state.options.geoContext, boundParams);

// Update the stepping state
state.particleHypothesis = particleHypothesis;

state.pathAccumulated = 0;
state.nSteps = 0;
state.nStepTrials = 0;
state.stepSize = ConstrainedStep(state.options.maxStepSize);
state.previousStepSize = 0;
state.statistics = StepperStatistics();

state.pars = freeParams;
state.cov = cov;
state.stepSize = ConstrainedStep(stepSize);
state.pathAccumulated = 0.;

// Reinitialize the stepping jacobian
state.jacToGlobal = surface.boundToFreeJacobian(
state.options.geoContext, freeParams.template segment<3>(eFreePos0),
freeParams.template segment<3>(eFreeDir0));
state.jacobian = BoundMatrix::Identity();
state.jacTransport = FreeMatrix::Identity();
state.derivative = FreeVector::Zero();
// Init the jacobian matrix if needed
state.covTransport = cov.has_value();
if (state.covTransport) {
state.cov = *cov;
state.jacToGlobal = surface.boundToFreeJacobian(
state.options.geoContext, freeParams.segment<3>(eFreePos0),
freeParams.segment<3>(eFreeDir0));
state.jacobian = BoundMatrix::Identity();
state.jacTransport = FreeMatrix::Identity();
state.derivative = FreeVector::Zero();
}
}

template <typename E>
Expand Down
44 changes: 17 additions & 27 deletions Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,6 +198,9 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
SingleState state;
double weight;
IntersectionStatus status;

Component(SingleState state_, double weight_, IntersectionStatus status_)
: state(std::move(state_)), weight(weight_), status(status_) {}
};

Options options;
Expand Down Expand Up @@ -240,49 +243,34 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
getDefaultLogger("GSF", Logging::INFO))
: EigenStepper<extension_t>(config), m_logger(std::move(logger)) {}

/// Construct and initialize a state
State makeState(const Options& options,
State makeState(const Options& options) const {
State state(options);
return state;
}

void initialize(State& state,
const MultiComponentBoundTrackParameters& par) const {
if (par.components().empty()) {
throw std::invalid_argument(
"Cannot construct MultiEigenStepperLoop::State with empty "
"multi-component parameters");
}

State state(options);

state.particleHypothesis = par.particleHypothesis();

const auto surface = par.referenceSurface().getSharedPtr();

for (auto i = 0ul; i < par.components().size(); ++i) {
const auto& [weight, singlePars] = par[i];
state.components.push_back({SingleStepper::makeState(options, singlePars),
weight, IntersectionStatus::onSurface});
auto& cmp =
state.components.emplace_back(SingleStepper::makeState(state.options),
weight, IntersectionStatus::onSurface);
SingleStepper::initialize(cmp.state, singlePars);
}

if (std::get<2>(par.components().front())) {
state.covTransport = true;
}

return state;
}

/// @brief Resets the state
///
/// @param [in, out] state State of the stepper
/// @param [in] boundParams Parameters in bound parametrisation
/// @param [in] cov Covariance matrix
/// @param [in] surface The reference surface of the bound parameters
/// @param [in] stepSize Step size
void resetState(
State& state, const BoundVector& boundParams,
const BoundSquareMatrix& cov, const Surface& surface,
const double stepSize = std::numeric_limits<double>::max()) const {
for (auto& component : state.components) {
SingleStepper::resetState(component.state, boundParams, cov, surface,
stepSize);
}
}

/// A proxy struct which allows access to a single component of the
Expand Down Expand Up @@ -420,8 +408,10 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
Result<ComponentProxy> addComponent(State& state,
const BoundTrackParameters& pars,
double weight) const {
state.components.push_back({SingleStepper::makeState(state.options, pars),
weight, IntersectionStatus::onSurface});
auto& cmp =
state.components.emplace_back(SingleStepper::makeState(state.options),
weight, IntersectionStatus::onSurface);
SingleStepper::initialize(cmp.state, pars);

return ComponentProxy{state.components.back(), state};
}
Expand Down
5 changes: 3 additions & 2 deletions Core/include/Acts/Propagator/Propagator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -386,8 +386,9 @@ class Propagator final
private:
const Logger& logger() const { return *m_logger; }

template <typename propagator_state_t, typename path_aborter_t>
void initialize(propagator_state_t& state) const;
template <typename propagator_state_t, typename parameters_t,
typename path_aborter_t>
void initialize(propagator_state_t& state, const parameters_t& start) const;

template <typename propagator_state_t, typename result_t>
void moveStateToResult(propagator_state_t& state, result_t& result) const;
Expand Down
Loading
Loading