Skip to content

Commit

Permalink
refactor: Add method to reset jacobian in stepper (#1171)
Browse files Browse the repository at this point in the history
When trying to use the ATLAS stepper with the KF I noticed that it wouldn't work since we use to set the stepping state jacobian to a BoundMatrix::Identity() which isn't compatible with the ATLAS stepper.
This PR add a method to the stepper to set the jacobian to the identity, and thus solve this issue. 
Similarly I had to transform the update to a reset for the propagation to the Perigee as it relied on manually setting some matrix afterward (which was not compatible with the ATLAS stepper).
  • Loading branch information
Corentin-Allaire committed Feb 21, 2022
1 parent 9294d61 commit f9dbc02
Show file tree
Hide file tree
Showing 4 changed files with 66 additions and 13 deletions.
48 changes: 48 additions & 0 deletions Core/include/Acts/Propagator/AtlasStepper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -1376,6 +1376,54 @@ class AtlasStepper {
return h;
}

/// Method that reset the Jacobian to the Identity for when no bound state are
/// available
///
/// @param [in,out] state State of the stepper
void setIdentityJacobian(State& state) const {
state.jacobian[0] = 1.; // dL0/dL0
state.jacobian[1] = 0.; // dL0/dL1
state.jacobian[2] = 0.; // dL0/dPhi
state.jacobian[3] = 0.; // dL0/dThe
state.jacobian[4] = 0.; // dL0/dCM
state.jacobian[5] = 0.; // dL0/dT

state.jacobian[6] = 0.; // dL1/dL0
state.jacobian[7] = 1.; // dL1/dL1
state.jacobian[8] = 0.; // dL1/dPhi
state.jacobian[9] = 0.; // dL1/dThe
state.jacobian[10] = 0.; // dL1/dCM
state.jacobian[11] = 0.; // dL1/dT

state.jacobian[12] = 0.; // dPhi/dL0
state.jacobian[13] = 0.; // dPhi/dL1
state.jacobian[14] = 1.; // dPhi/dPhi
state.jacobian[15] = 0.; // dPhi/dThe
state.jacobian[16] = 0.; // dPhi/dCM
state.jacobian[17] = 0.; // dPhi/dT

state.jacobian[18] = 0.; // dThe/dL0
state.jacobian[19] = 0.; // dThe/dL1
state.jacobian[20] = 0.; // dThe/dPhi
state.jacobian[21] = 1.; // dThe/dThe
state.jacobian[22] = 0.; // dThe/dCM
state.jacobian[23] = 0.; // dThe/dT

state.jacobian[24] = 0.; // dCM /dL0
state.jacobian[25] = 0.; // dCM /dL1
state.jacobian[26] = 0.; // dCM /dPhi
state.jacobian[27] = 0.; // dCM /dTheta
state.jacobian[28] = 1.; // dCM /dCM
state.jacobian[29] = 0.; // dCM/dT

state.jacobian[30] = 0.; // dT/dL0
state.jacobian[31] = 0.; // dT/dL1
state.jacobian[32] = 0.; // dT/dPhi
state.jacobian[33] = 0.; // dT/dThe
state.jacobian[34] = 0.; // dT/dCM
state.jacobian[35] = 1.; // dT/dT
}

private:
std::shared_ptr<const MagneticFieldProvider> m_bField;

Expand Down
6 changes: 6 additions & 0 deletions Core/include/Acts/Propagator/EigenStepper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,12 @@ class EigenStepper {
template <typename propagator_state_t>
Result<double> step(propagator_state_t& state) const;

/// Method that reset the Jacobian to the Identity for when no bound state are
/// available
///
/// @param [in,out] state State of the stepper
void setIdentityJacobian(State& state) const;

protected:
/// Magnetic field inside of the detector
std::shared_ptr<const MagneticFieldProvider> m_bField;
Expand Down
5 changes: 5 additions & 0 deletions Core/include/Acts/Propagator/EigenStepper.ipp
Original file line number Diff line number Diff line change
Expand Up @@ -265,3 +265,8 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
}
return h;
}

template <typename E, typename A>
void Acts::EigenStepper<E, A>::setIdentityJacobian(State& state) const {
state.jacobian = BoundMatrix::Identity();
}
20 changes: 7 additions & 13 deletions Core/include/Acts/TrackFitting/KalmanFitter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -814,7 +814,7 @@ class KalmanFitter {
}
// Not creating bound state here, so need manually reinitialize
// jacobian
state.stepping.jacobian = BoundMatrix::Identity();
stepper.setIdentityJacobian(state.stepping);
if (surface->surfaceMaterial() != nullptr) {
// Update state and stepper with material effects
materialInteractor(surface, state, stepper);
Expand Down Expand Up @@ -972,16 +972,14 @@ class KalmanFitter {
(std::abs(firstIntersection.intersection.pathLength) <=
std::abs(lastIntersection.intersection.pathLength));
if (closerTofirstCreatedState) {
stepper.update(state.stepping, firstParams,
firstCreatedState.smoothed(),
firstCreatedState.smoothedCovariance(),
firstCreatedState.referenceSurface());
stepper.resetState(state.stepping, firstCreatedState.smoothed(),
firstCreatedState.smoothedCovariance(),
firstCreatedState.referenceSurface());
reverseDirection = (firstIntersection.intersection.pathLength < 0);
} else {
stepper.update(state.stepping, lastParams,
lastCreatedMeasurement.smoothed(),
lastCreatedMeasurement.smoothedCovariance(),
lastCreatedMeasurement.referenceSurface());
stepper.resetState(state.stepping, lastCreatedMeasurement.smoothed(),
lastCreatedMeasurement.smoothedCovariance(),
lastCreatedMeasurement.referenceSurface());
reverseDirection = (lastIntersection.intersection.pathLength < 0);
}
const auto& surface = closerTofirstCreatedState
Expand All @@ -1000,10 +998,6 @@ class KalmanFitter {
state.stepping.navDir =
(state.stepping.navDir == forward) ? backward : forward;
}
// Reinitialize the stepping jacobian
state.stepping.jacobian = BoundMatrix::Identity();
state.stepping.jacTransport = FreeMatrix::Identity();
state.stepping.derivative = FreeVector::Zero();
// Reset the step size
state.stepping.stepSize = ConstrainedStep(
state.stepping.navDir * std::abs(state.options.maxStepSize));
Expand Down

0 comments on commit f9dbc02

Please sign in to comment.