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: qOverP and absoluteMomentum for Stepper interface #2253

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.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion Core/include/Acts/Navigation/NextNavigator.hpp
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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