Skip to content

Commit

Permalink
refactor: qOverP and absoluteMomentum for Stepper interface (#2253)
Browse files Browse the repository at this point in the history
With this PR I try to align the naming in the Stepper interface to `qOverP` and `absoluteMomentum` which is also my goal for the track parameters.

I pulled these changes out of #2181.
  • Loading branch information
andiwand committed Jun 30, 2023
1 parent 7b0e2ed commit f22a278
Show file tree
Hide file tree
Showing 26 changed files with 106 additions and 85 deletions.
2 changes: 1 addition & 1 deletion Core/include/Acts/Navigation/NextNavigator.hpp
Expand Up @@ -411,7 +411,7 @@ class NextNavigator {
NavigationState& nState) const {
nState.position = stepper.position(state.stepping);
nState.direction = stepper.direction(state.stepping);
nState.absMomentum = stepper.momentum(state.stepping);
nState.absMomentum = stepper.absoluteMomentum(state.stepping);
auto fieldResult = stepper.getField(state.stepping, nState.position);
if (!fieldResult.ok()) {
ACTS_ERROR(volInfo(state) << posInfo(state, stepper)
Expand Down
14 changes: 8 additions & 6 deletions Core/include/Acts/Propagator/AtlasStepper.hpp
Expand Up @@ -353,14 +353,16 @@ class AtlasStepper {
return Vector3(state.pVector[4], state.pVector[5], state.pVector[6]);
}

double qop(const State& state) const { return state.pVector[7]; }
double qOverP(const State& state) const { return state.pVector[7]; }

double momentum(const State& state) const {
return 1. / std::abs(qop(state));
double absoluteMomentum(const State& state) const {
return 1. / std::abs(qOverP(state));
}

/// Charge access
double charge(const State& state) const { return qop(state) > 0. ? 1. : -1.; }
double charge(const State& state) const {
return qOverP(state) > 0. ? 1. : -1.;
}

/// Overstep limit
double overstepLimit(const State& /*state*/) const { return m_overstepLimit; }
Expand Down Expand Up @@ -1248,7 +1250,7 @@ class AtlasStepper {

// Evaluate the time propagation
double dtds =
std::hypot(1, state.options.mass / momentum(state.stepping));
std::hypot(1, state.options.mass / absoluteMomentum(state.stepping));
state.stepping.pVector[3] += h * dtds;
state.stepping.pVector[59] = dtds;
state.stepping.field = f;
Expand All @@ -1257,7 +1259,7 @@ class AtlasStepper {
if (Jac) {
double dtdl = h * state.options.mass * state.options.mass *
charge(state.stepping) /
(momentum(state.stepping) * dtds);
(absoluteMomentum(state.stepping) * dtds);
state.stepping.pVector[43] += dtdl;

// Jacobian calculation
Expand Down
12 changes: 6 additions & 6 deletions Core/include/Acts/Propagator/EigenStepper.hpp
Expand Up @@ -216,21 +216,21 @@ class EigenStepper {
/// QoP direction accessor
///
/// @param state [in] The stepping state (thread-local cache)
double qop(const State& state) const { return state.pars[eFreeQOverP]; }
double qOverP(const State& state) const { return state.pars[eFreeQOverP]; }

/// Absolute momentum accessor
///
/// @param state [in] The stepping state (thread-local cache)
double momentum(const State& state) const {
double absoluteMomentum(const State& state) const {
auto q = charge(state);
return std::abs((q == 0 ? 1 : q) / qop(state));
return std::abs((q == 0 ? 1 : q) / qOverP(state));
}

/// Charge access
///
/// @param state [in] The stepping state (thread-local cache)
double charge(const State& state) const {
return std::copysign(state.absCharge, qop(state));
return std::copysign(state.absCharge, qOverP(state));
}

/// Time access
Expand Down Expand Up @@ -363,10 +363,10 @@ class EigenStepper {
/// @param [in,out] state State object that will be updated
/// @param [in] uposition the updated position
/// @param [in] udirection the updated direction
/// @param [in] qop the updated qop value
/// @param [in] qOverP the updated qOverP value
/// @param [in] time the updated time value
void update(State& state, const Vector3& uposition, const Vector3& udirection,
double qop, double time) const;
double qOverP, double time) const;

/// Method for on-demand transport of the covariance
/// to a new curvilinear frame at current position,
Expand Down
4 changes: 2 additions & 2 deletions Core/include/Acts/Propagator/EigenStepper.ipp
Expand Up @@ -84,12 +84,12 @@ void Acts::EigenStepper<E, A>::update(State& state,

template <typename E, typename A>
void Acts::EigenStepper<E, A>::update(State& state, const Vector3& uposition,
const Vector3& udirection, double qop,
const Vector3& udirection, double qOverP,
double time) const {
state.pars.template segment<3>(eFreePos0) = uposition;
state.pars.template segment<3>(eFreeDir0) = udirection;
state.pars[eFreeTime] = time;
state.pars[eFreeQOverP] = qop;
state.pars[eFreeQOverP] = qOverP;
}

template <typename E, typename A>
Expand Down
11 changes: 9 additions & 2 deletions Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
Expand Up @@ -66,7 +66,7 @@ struct WeightedComponentReducerLoop {
// TODO: Maybe we can cache this value and only update it when the parameters
// change
template <typename stepper_state_t>
static ActsScalar qop(const stepper_state_t& s) {
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 {
Expand Down Expand Up @@ -627,10 +627,17 @@ class MultiEigenStepperLoop
return Reducer::direction(state);
}

/// QoP access
///
/// @param state [in] The stepping state (thread-local cache)
double qOverP(const State& state) const { return Reducer::qOverP(state); }

/// Absolute momentum accessor
///
/// @param state [in] The stepping state (thread-local cache)
double momentum(const State& state) const { return Reducer::momentum(state); }
double absoluteMomentum(const State& state) const {
return Reducer::momentum(state);
}

/// Charge access
///
Expand Down
2 changes: 1 addition & 1 deletion Core/include/Acts/Propagator/StandardAborters.hpp
Expand Up @@ -224,7 +224,7 @@ struct ParticleStopped {
bool operator()(propagator_state_t& state, const stepper_t& stepper,
const navigator_t& navigator,
const Logger& /*logger*/) const {
if (stepper.momentum(state.stepping) > 0) {
if (stepper.absoluteMomentum(state.stepping) > 0) {
return false;
}
navigator.targetReached(state.navigation, true);
Expand Down
12 changes: 8 additions & 4 deletions Core/include/Acts/Propagator/StepperConcept.hpp
Expand Up @@ -40,7 +40,8 @@ METHOD_TRAIT(reset_state_t, resetState);
METHOD_TRAIT(get_field_t, getField);
METHOD_TRAIT(position_t, position);
METHOD_TRAIT(direction_t, direction);
METHOD_TRAIT(momentum_t, momentum);
METHOD_TRAIT(qop_t, qOverP);
METHOD_TRAIT(absolute_momentum_t, absoluteMomentum);
METHOD_TRAIT(charge_t, charge);
METHOD_TRAIT(time_t, time);
METHOD_TRAIT(overstep_t, overstepLimit);
Expand Down Expand Up @@ -107,8 +108,10 @@ constexpr bool MultiStepperStateConcept= require<
static_assert(position_exists, "position method not found");
constexpr static bool direction_exists = has_method<const S, Vector3, direction_t, const state&>;
static_assert(direction_exists, "direction method not found");
constexpr static bool momentum_exists = has_method<const S, double, momentum_t, const state&>;
static_assert(momentum_exists, "momentum method not found");
constexpr static bool qop_exists = has_method<const S, double, qop_t, const state&>;
static_assert(qop_exists, "qOverP method not found");
constexpr static bool absolute_momentum_exists = has_method<const S, double, absolute_momentum_t, const state&>;
static_assert(absolute_momentum_exists, "absoluteMomentum method not found");
constexpr static bool charge_exists = has_method<const S, double, charge_t, const state&>;
static_assert(charge_exists, "charge method not found");
constexpr static bool time_exists = has_method<const S, double, time_t, const state&>;
Expand Down Expand Up @@ -140,7 +143,8 @@ constexpr bool MultiStepperStateConcept= require<
curvilinear_state_exists,
position_exists,
direction_exists,
momentum_exists,
qop_exists,
absolute_momentum_exists,
charge_exists,
time_exists,
bound_state_method_exists,
Expand Down
10 changes: 5 additions & 5 deletions Core/include/Acts/Propagator/StraightLineStepper.hpp
Expand Up @@ -196,21 +196,21 @@ class StraightLineStepper {
/// QoP direction accessor
///
/// @param state [in] The stepping state (thread-local cache)
double qop(const State& state) const { return state.pars[eFreeQOverP]; }
double qOverP(const State& state) const { return state.pars[eFreeQOverP]; }

/// Absolute momentum accessor
///
/// @param state [in] The stepping state (thread-local cache)
double momentum(const State& state) const {
double absoluteMomentum(const State& state) const {
auto q = charge(state);
return std::abs((q == 0 ? 1 : q) / qop(state));
return std::abs((q == 0 ? 1 : q) / qOverP(state));
}

/// Charge access
///
/// @param state [in] The stepping state (thread-local cache)
double charge(const State& state) const {
return std::copysign(state.absCharge, qop(state));
return std::copysign(state.absCharge, qOverP(state));
}

/// Time access
Expand Down Expand Up @@ -387,7 +387,7 @@ class StraightLineStepper {
const navigator_t& /*navigator*/) const {
// use the adjusted step size
const auto h = state.stepping.stepSize.value();
const double p = momentum(state.stepping);
const double p = absoluteMomentum(state.stepping);
// time propagates along distance as 1/b = sqrt(1 + m²/p²)
const auto dtds = std::hypot(1., state.options.mass / p);
// Update the track parameters according to the equations of motion
Expand Down
14 changes: 7 additions & 7 deletions Core/include/Acts/Propagator/detail/GenericDefaultExtension.hpp
Expand Up @@ -62,7 +62,7 @@ struct GenericDefaultExtension {
const navigator_t& /*navigator*/, ThisVector3& knew,
const Vector3& bField, std::array<Scalar, 4>& kQoP, const int i = 0,
const double h = 0., const ThisVector3& kprev = ThisVector3::Zero()) {
auto qop = stepper.qop(state.stepping);
auto qop = stepper.qOverP(state.stepping);
// First step does not rely on previous data
if (i == 0) {
knew = qop * stepper.direction(state.stepping).cross(bField);
Expand Down Expand Up @@ -137,7 +137,7 @@ struct GenericDefaultExtension {
/// = sqrt(m^2/p^2 + c^{-2}) with the mass m and the momentum p.
using std::hypot;
auto derivative =
hypot(1, state.options.mass / stepper.momentum(state.stepping));
hypot(1, state.options.mass / stepper.absoluteMomentum(state.stepping));
state.stepping.pars[eFreeTime] += h * derivative;
if (state.stepping.covTransport) {
state.stepping.derivative(3) = derivative;
Expand Down Expand Up @@ -182,7 +182,7 @@ struct GenericDefaultExtension {

auto& sd = state.stepping.stepData;
auto dir = stepper.direction(state.stepping);
auto qop = stepper.qop(state.stepping);
auto qop = stepper.qOverP(state.stepping);

D = FreeMatrix::Identity();

Expand Down Expand Up @@ -241,10 +241,10 @@ struct GenericDefaultExtension {

dGdL = h / 6. * (dk1dL + 2. * (dk2dL + dk3dL) + dk4dL);

D(3, 7) =
h * state.options.mass * state.options.mass *
stepper.qop(state.stepping) /
std::hypot(1., state.options.mass / stepper.momentum(state.stepping));
D(3, 7) = h * state.options.mass * state.options.mass *
stepper.qOverP(state.stepping) /
std::hypot(1., state.options.mass /
stepper.absoluteMomentum(state.stepping));
return true;
}
};
Expand Down
Expand Up @@ -79,7 +79,8 @@ struct GenericDenseEnvironmentExtension {
const navigator_t& navigator) const {
// Check for valid particle properties
if (stepper.charge(state.stepping) == 0. || state.options.mass == 0. ||
stepper.momentum(state.stepping) < state.options.momentumCutOff) {
stepper.absoluteMomentum(state.stepping) <
state.options.momentumCutOff) {
return 0;
}

Expand Down Expand Up @@ -123,7 +124,7 @@ struct GenericDenseEnvironmentExtension {
auto volumeMaterial = navigator.currentVolumeMaterial(state.navigation);
ThisVector3 position = stepper.position(state.stepping);
material = (volumeMaterial->material(position.template cast<double>()));
initialMomentum = stepper.momentum(state.stepping);
initialMomentum = stepper.absoluteMomentum(state.stepping);
currentMomentum = initialMomentum;
qop[0] = q / initialMomentum;
initializeEnergyLoss(state);
Expand Down Expand Up @@ -173,7 +174,7 @@ struct GenericDenseEnvironmentExtension {
const navigator_t& /*navigator*/, const double h) const {
// Evaluate the new momentum
auto newMomentum =
stepper.momentum(state.stepping) +
stepper.absoluteMomentum(state.stepping) +
(h / 6.) * (dPds[0] + 2. * (dPds[1] + dPds[2]) + dPds[3]);

// Break propagation if momentum becomes below cut-off
Expand Down
2 changes: 1 addition & 1 deletion Core/include/Acts/Propagator/detail/LoopProtection.hpp
Expand Up @@ -39,7 +39,7 @@ void setupLoopProtection(propagator_state_t& state, const stepper_t& stepper,
}

// Transverse component at start is taken for the loop protection
const double p = stepper.momentum(state.stepping);
const double p = stepper.absoluteMomentum(state.stepping);
// Calculate the full helix path
const double helixPath = state.stepping.navDir * 2 * M_PI * p / B;
// And set it as the loop limit if it overwrites the internal limit
Expand Down
Expand Up @@ -79,9 +79,9 @@ struct PointwiseMaterialInteraction {
pos(stepper.position(state.stepping)),
time(stepper.time(state.stepping)),
dir(stepper.direction(state.stepping)),
qOverP(stepper.qop(state.stepping)),
qOverP(stepper.qOverP(state.stepping)),
absQ(std::abs(stepper.charge(state.stepping))),
momentum(stepper.momentum(state.stepping)),
momentum(stepper.absoluteMomentum(state.stepping)),
mass(state.options.mass),
pdg(state.options.absPdgCode),
performCovarianceTransport(state.stepping.covTransport),
Expand Down
4 changes: 2 additions & 2 deletions Core/include/Acts/Propagator/detail/SteppingLogger.hpp
Expand Up @@ -71,8 +71,8 @@ struct SteppingLogger {
Step step;
step.stepSize = state.stepping.stepSize;
step.position = stepper.position(state.stepping);
step.momentum =
stepper.momentum(state.stepping) * stepper.direction(state.stepping);
step.momentum = stepper.absoluteMomentum(state.stepping) *
stepper.direction(state.stepping);

if (navigator.currentSurface(state.navigation) != nullptr) {
// hang on to surface
Expand Down
Expand Up @@ -63,9 +63,9 @@ struct VolumeMaterialInteraction {
pos(stepper.position(state.stepping)),
time(stepper.time(state.stepping)),
dir(stepper.direction(state.stepping)),
qOverP(stepper.qop(state.stepping)),
qOverP(stepper.qOverP(state.stepping)),
absQ(std::abs(stepper.charge(state.stepping))),
momentum(stepper.momentum(state.stepping)),
momentum(stepper.absoluteMomentum(state.stepping)),
mass(state.options.mass),
pdg(state.options.absPdgCode),
performCovarianceTransport(state.stepping.covTransport),
Expand Down
3 changes: 2 additions & 1 deletion Core/include/Acts/TrackFitting/KalmanFitter.hpp
Expand Up @@ -344,7 +344,8 @@ class KalmanFitter {
ACTS_VERBOSE("KalmanFitter step at pos: "
<< stepper.position(state.stepping).transpose()
<< " dir: " << stepper.direction(state.stepping).transpose()
<< " momentum: " << stepper.momentum(state.stepping));
<< " momentum: "
<< stepper.absoluteMomentum(state.stepping));

// Add the measurement surface as external surface to navigator.
// We will try to hit those surface by ignoring boundary checks.
Expand Down
6 changes: 3 additions & 3 deletions Core/include/Acts/TrackFitting/detail/GsfUtils.hpp
Expand Up @@ -123,7 +123,7 @@ class ScopedGsfInfoPrinterAndChecker {
: m_state(state),
m_stepper(stepper),
m_navigator(navigator),
m_p_initial(stepper.momentum(state.stepping)),
m_p_initial(stepper.absoluteMomentum(state.stepping)),
m_logger{logger} {
// Some initial printing
checks(true);
Expand All @@ -132,15 +132,15 @@ class ScopedGsfInfoPrinterAndChecker {
<< stepper.position(state.stepping).transpose()
<< " with direction "
<< stepper.direction(state.stepping).transpose()
<< " and momentum " << stepper.momentum(state.stepping)
<< " and momentum " << stepper.absoluteMomentum(state.stepping)
<< " and charge " << stepper.charge(state.stepping));
ACTS_VERBOSE("Propagation is in " << state.stepping.navDir << " mode");
print_component_stats();
}

~ScopedGsfInfoPrinterAndChecker() {
if (m_navigator.currentSurface(m_state.navigation)) {
const auto p_final = m_stepper.momentum(m_state.stepping);
const auto p_final = m_stepper.absoluteMomentum(m_state.stepping);
ACTS_VERBOSE("Component status at end of step:");
print_component_stats();
ACTS_VERBOSE("Delta Momentum = " << std::setprecision(5)
Expand Down
Expand Up @@ -219,7 +219,7 @@ struct SimulationActor {
return Particle(previous)
.setPosition4(stepper.position(state), stepper.time(state))
.setDirection(stepper.direction(state))
.setAbsoluteMomentum(stepper.momentum(state))
.setAbsoluteMomentum(stepper.absoluteMomentum(state))
.setProperTime(properTime);
}

Expand Down
Expand Up @@ -94,8 +94,8 @@ struct AutodiffExtensionWrapper {
// A fake stepper
struct FakeStepper {
auto charge(const FakeStepperState& s) const { return s.q; }
auto qop(const FakeStepperState& s) const { return s.pars(eFreeQOverP); }
auto momentum(const FakeStepperState& s) const {
auto qOverP(const FakeStepperState& s) const { return s.pars(eFreeQOverP); }
auto absoluteMomentum(const FakeStepperState& s) const {
return s.q / s.pars(eFreeQOverP);
}
auto direction(const FakeStepperState& s) const {
Expand Down Expand Up @@ -128,7 +128,7 @@ struct AutodiffExtensionWrapper {
initial_params.segment<3>(eFreeDir0) = stepper.direction(state.stepping);
initial_params(eFreeQOverP) =
(fstate.stepping.q != 0. ? fstate.stepping.q : 1.) /
stepper.momentum(state.stepping);
stepper.absoluteMomentum(state.stepping);

const auto& sd = state.stepping.stepData;

Expand Down

0 comments on commit f22a278

Please sign in to comment.