update

Author: benjaminhuth

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -240,6 +240,9 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
       return EigenStepperError::StepInvalid;
     }
   }
+  
+  std::cout << "EigenStepper initial h: " << initialH << "\n";
+  std::cout << "EigenStepper h: " << h << "\n";
 
   // Update the track parameters according to the equations of motion
   state.stepping.pars.template segment<3>(eFreePos0) +=

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -26,6 +26,7 @@
 #include "Acts/Propagator/detail/LoopStepperUtils.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/GaussianMixtureReduction.hpp"
+#include "Acts/Propagator/detail/MultiStepperUtils.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Result.hpp"
 
@@ -478,7 +479,7 @@ class MultiEigenStepperLoop
     }
   }
 
-  /// Updates the components in the multistepper
+ /// Updates the components in the multistepper
   ///
   /// @param [in,out] state  The stepping state (thread-local cache)
   /// @param [in] surface The surface we are on

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -8,6 +8,7 @@
 
 #include "Acts/Propagator/MultiEigenStepperLoop.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "Acts/Propagator/detail/MultiStepperUtils.hpp"
 
 namespace Acts {
 
@@ -17,76 +18,15 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
     const FreeToBoundCorrection& freeToBoundCorrection) const
     -> Result<BoundState> {
   assert(!state.components.empty());
-
-  std::vector<std::tuple<double, BoundVector, Covariance>> cmps;
-  cmps.reserve(numberComponents(state));
-  double accumulatedPathLength = 0.0;
-
-  for (auto i = 0ul; i < numberComponents(state); ++i) {
-    auto& cmpState = state.components[i].state;
-
-    // Force the component to be on the surface
-    // This needs to be done because of the `averageOnSurface`-option of the
-    // `MultiStepperSurfaceReached`-Aborter, which can be configured to end the
-    // propagation when the mean of all components reached the destination
-    // surface. Thus, it is not garantueed that all states are actually
-    // onSurface.
-    cmpState.pars.template segment<3>(eFreePos0) =
-        surface
-            .intersect(state.geoContext,
-                       cmpState.pars.template segment<3>(eFreePos0),
-                       cmpState.pars.template segment<3>(eFreeDir0), false)
-            .closest()
-            .position();
-
-    auto bs = SingleStepper::boundState(cmpState, surface, transportCov,
-                                        freeToBoundCorrection);
-
-    if (bs.ok()) {
-      const auto& btp = std::get<BoundTrackParameters>(*bs);
-      cmps.emplace_back(
-          state.components[i].weight, btp.parameters(),
-          btp.covariance().value_or(Acts::BoundSquareMatrix::Zero()));
-      accumulatedPathLength +=
-          std::get<double>(*bs) * state.components[i].weight;
-    }
-  }
-
-  if (cmps.empty()) {
-    return MultiStepperError::AllComponentsConversionToBoundFailed;
-  }
-
-  return BoundState{MultiComponentBoundTrackParameters(
-                        surface.getSharedPtr(), cmps, state.particleHypothesis),
-                    Jacobian::Zero(), accumulatedPathLength};
+  return detail::multiComponentBoundState(*this, state, surface, transportCov, freeToBoundCorrection);
 }
 
 template <typename E, typename R, typename A>
 auto MultiEigenStepperLoop<E, R, A>::curvilinearState(State& state,
                                                       bool transportCov) const
     -> CurvilinearState {
   assert(!state.components.empty());
-
-  std::vector<
-      std::tuple<double, Vector4, Vector3, ActsScalar, BoundSquareMatrix>>
-      cmps;
-  cmps.reserve(numberComponents(state));
-  double accumulatedPathLength = 0.0;
-
-  for (auto i = 0ul; i < numberComponents(state); ++i) {
-    const auto [cp, jac, pl] = SingleStepper::curvilinearState(
-        state.components[i].state, transportCov);
-
-    cmps.emplace_back(state.components[i].weight,
-                      cp.fourPosition(state.geoContext), cp.direction(),
-                      (cp.charge() / cp.absoluteMomentum()),
-                      cp.covariance().value_or(BoundSquareMatrix::Zero()));
-    accumulatedPathLength += state.components[i].weight * pl;
-  }
-
-  return CurvilinearState{
-      MultiComponentCurvilinearTrackParameters(cmps, state.particleHypothesis),
-      Jacobian::Zero(), accumulatedPathLength};
+  return detail::multiComponentCurvilinearState(*this, state, transportCov);
 }
 
 template <typename E, typename R, typename A>

Core/include/Acts/Propagator/MultiEigenStepperSIMD.hpp

@@ -14,7 +14,7 @@
 #include "Acts/Propagator/detail/SimdHelpers.hpp"
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Units.hpp"
-#include "Acts/EventData/MultiComponentBoundTrackParameters.hpp"
+#include "Acts/EventData/MultiComponentTrackParameters.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
 #include "Acts/Propagator/DefaultExtension.hpp"
@@ -27,6 +27,7 @@
 #include "Acts/Propagator/detail/SimdHelpers.hpp"
 #include "Acts/Propagator/detail/SimdStepperUtils.hpp"
 #include "Acts/Propagator/detail/SteppingHelper.hpp"
+#include "Acts/Propagator/detail/MultiStepperUtils.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Result.hpp"
 
@@ -124,16 +125,17 @@ class MultiEigenStepperSIMD
   using SingleState = typename SingleStepper::State;
 
   /// @brief Use the definitions from the Single-stepper
-  using typename SingleStepper::BoundState;
   using typename SingleStepper::Covariance;
-  using typename SingleStepper::CurvilinearState;
   using typename SingleStepper::Jacobian;
 
   /// @brief The reducer type
   using Reducer = component_reducer_t;
 
   /// @brief How many components can this stepper manage?
   static constexpr int maxComponents = NComponents;
+  
+  using BoundState = detail::MultiStepperBoundState;
+  using CurvilinearState = detail::MultiStepperCurvilinearState;
 
   /// @brief SIMD typedefs
   using SimdScalar = SimdType<NComponents>;
@@ -369,11 +371,13 @@ class MultiEigenStepperSIMD
   void update(State& state, const Surface& /*surface*/, iterator_t begin,
               iterator_t end, const copy_t& copy) const {
     assert(std::distance(begin, end) <= NComponents);
+    ACTS_VERBOSE("distance " << std::distance(begin, end));
 
     std::size_t i = 0;
     auto it = begin;
     for (; it != end; ++i, ++it) {
       auto proxy = ComponentProxy{state, i};
+      ACTS_VERBOSE("copy to component " << i);
       copy(*it, proxy);
     }
 
@@ -499,9 +503,9 @@ class MultiEigenStepperSIMD
   /// @param surface [in] The surface provided
   /// @param bcheck [in] The boundary check for this status update
   Intersection3D::Status updateSurfaceStatus(
-      State& state, const Surface& /*surface*/, Direction /*navDir*/,
-      const BoundaryCheck& /*bcheck*/,
-      const Logger& /*logger*/ = getDummyLogger()) const;
+      State& state, const Surface& surface, std::uint8_t index, Direction navDir,
+      const BoundaryCheck& bcheck,
+      const Logger& logger = getDummyLogger()) const;
 
   /// Update step size
   ///
@@ -608,13 +612,10 @@ class MultiEigenStepperSIMD
   ///   - the stepwise jacobian towards it (from last bound)
   ///   - and the path length (from start - for ordering)
   Result<BoundState> boundState(
-      State& /*state*/, const Surface& /*surface*/,
-      bool /*transportCov*/ = true,
-      const FreeToBoundCorrection& /*freeToBoundCorrection*/ =
-          FreeToBoundCorrection(false)) const {
-    throw std::runtime_error(
-        "'boundState' not yet implemented for MultiEigenStepper");
-  }
+      State& state, const Surface& surface,
+      bool transportCov = true,
+      const FreeToBoundCorrection& freeToBoundCorrection =
+          FreeToBoundCorrection(false)) const;
 
   /// Create and return a curvilinear state at the current position
   ///
@@ -630,57 +631,7 @@ class MultiEigenStepperSIMD
   ///   - and the path length (from start - for ordering)
   /// TODO reformulate with reducer functions
   CurvilinearState curvilinearState(State& /*state*/,
-                                    bool /*transportCov*/ = true) const {
-    throw std::runtime_error("not implemented in Multi Stepper");
-    // // std optional because CurvilinearState is not default constructable
-    // std::array<std::optional<CurvilinearState>, NComponents> states;
-    //
-    // // Compute all states
-    // for (auto i = 0ul; i < NComponents; ++i) {
-    //   FreeVector pars, derivative;
-    //   FreeMatrix jacTransport;
-    //
-    //   for (auto j = 0ul; j < eFreeSize; ++j) {
-    //     pars[j] = state.pars[j][i];
-    //   }
-    //
-    //   for (auto j = 0ul; j < eFreeSize; ++j)
-    //     for (auto k = 0ul; k < eFreeSize; ++k) {
-    //       jacTransport(j, k) = state.jacTransport(j, k)[i];
-    //     }
-    //
-    //   states[i] = detail::curvilinearState(
-    //       state.covs[i], state.jacobians[i], jacTransport, derivative,
-    //       state.jacToGlobals[i], pars, state.covTransport && transportCov,
-    //       state.pathAccumulated);
-    // }
-    //
-    // // Sum everything up
-    // Vector4 pos = Vector4::Zero();
-    // Vector3 dir = Vector3::Zero();
-    // double p = 0., pathlen = 0.;
-    // Jacobian jac = Jacobian::Zero();
-    //
-    // for (const auto& curvstate : states) {
-    //   const auto& curvpars =
-    //   std::get<CurvilinearTrackParameters>(*curvstate); pos +=
-    //   curvpars.fourPosition(state.geoContext); dir +=
-    //   curvpars.unitDirection(); p += curvpars.absoluteabsoluteMomentum();
-    //   pathlen += std::get<double>(*curvstate);
-    //   jac += std::get<Jacobian>(*curvstate);
-    // }
-    //
-    // // Average over all
-    // double q = std::get<double>(*states.front());
-    // pos /= NComponents;
-    // dir.normalize();
-    // p /= NComponents;
-    // pathlen /= NComponents;
-    // jac /= NComponents;
-    //
-    // return CurvilinearState{CurvilinearTrackParameters(pos, dir, p, q), jac,
-    //                         pathlen};
-  }
+                                    bool /*transportCov*/ = true) const;
 
   /// Method for on-demand transport of the covariance
   /// to a new curvilinear frame at current  position,

Core/include/Acts/Propagator/MultiEigenStepperSIMD.ipp

@@ -7,6 +7,7 @@
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include "Acts/Utilities/Zip.hpp"
+#include "Acts/Propagator/detail/MultiStepperUtils.hpp"
 
 namespace Acts {
 
@@ -65,8 +66,8 @@ void MultiEigenStepperSIMD<N, AA, BB, CC, DD>::transportCovarianceToCurvilinear(
 template <int N, typename AA, typename BB, typename CC, typename DD>
 Intersection3D::Status
 MultiEigenStepperSIMD<N, AA, BB, CC, DD>::updateSurfaceStatus(
-    State& state, const Surface& surface, Direction navDir,
-    const BoundaryCheck& bcheck, const Logger& l) const {
+    State& state, const Surface& surface, std::uint8_t index, Direction navDir,
+    const BoundaryCheck& bcheck, const Logger& /*l*/) const {
   ACTS_DEBUG("MultiEigenStepperSIMD::updateSurfaceStatus");
 
   std::array<int, 4> counts = {0, 0, 0, 0};
@@ -80,7 +81,7 @@ MultiEigenStepperSIMD<N, AA, BB, CC, DD>::updateSurfaceStatus(
     const auto prevStatus = cmp.status();
 
     cmp.status() = detail::updateSingleSurfaceStatus<SingleProxyStepper>(
-        cmp.singleStepper(*this), state, surface, navDir, bcheck,
+        cmp.singleStepper(*this), state, surface, index, navDir, bcheck,
         s_onSurfaceTolerance, logger()  /*Acts::getDummyLogger()*/);
 
     ACTS_VERBOSE("  cmp" << cmp.index() << ": " << prevStatus << " -> " << cmp.status());
@@ -103,6 +104,21 @@ MultiEigenStepperSIMD<N, AA, BB, CC, DD>::updateSurfaceStatus(
     return Status::missed;
 }
 
+template <int N, typename AA, typename BB, typename CC, typename DD>
+auto MultiEigenStepperSIMD<N, AA, BB, CC, DD>::boundState(
+      State& state, const Surface& surface,
+      bool transportCov,
+      const FreeToBoundCorrection& freeToBoundCorrection) const -> Result<BoundState> {
+  return detail::multiComponentBoundState(*this, state, surface, transportCov, freeToBoundCorrection);
+}
+
+template <int N, typename AA, typename BB, typename CC, typename DD>
+auto MultiEigenStepperSIMD<N, AA, BB, CC, DD>::curvilinearState(
+      State& state, 
+      bool transportCov) const -> CurvilinearState {
+  return detail::multiComponentCurvilinearState(*this, state, transportCov);
+}
+
 // Seperated stepsize estimate from Single eigen stepper
 /// TODO state should be constant, but then magne
 template <int N, typename AA, typename BB, typename CC, typename DD>
@@ -184,7 +200,7 @@ Result<double> MultiEigenStepperSIMD<N, AA, BB, CC, DD>::estimate_step_size(
                                          0.25)),
                  4.);
 
-    current_estimate.setValue(current_estimate.value() * stepSizeScaling);
+    current_estimate.update(current_estimate.value() * stepSizeScaling, ConstrainedStep::actor);
 
     // If step size becomes too small the particle remains at the initial
     // place
@@ -219,9 +235,15 @@ Result<double> MultiEigenStepperSIMD<N, AA, BB, CC, DD>::step(
   const auto dir = multiDirection(stepping);
 
   ACTS_VERBOSE("Pos before step:");
-  ACTS_VERBOSE("  x = " << multiPosition(stepping)[0]);
-  ACTS_VERBOSE("  y = " << multiPosition(stepping)[1]);
-  ACTS_VERBOSE("  z = " << multiPosition(stepping)[2]);
+  ACTS_VERBOSE("  x = " << pos[0]);
+  ACTS_VERBOSE("  y = " << pos[1]);
+  ACTS_VERBOSE("  z = " << pos[2]);
+  
+  
+  ACTS_VERBOSE("Dir before step:");
+  ACTS_VERBOSE("  dx = " << dir[0]);
+  ACTS_VERBOSE("  dy = " << dir[1]);
+  ACTS_VERBOSE("  dz = " << dir[2]);
 
   // First Runge-Kutta point
   sd.B_first = getMultiField(stepping, pos);
@@ -270,7 +292,7 @@ Result<double> MultiEigenStepperSIMD<N, AA, BB, CC, DD>::step(
   // }();
 
   SimdScalar h = 0;
-  for (int i = 0; i < stepping.stepSizes.size(); ++i) {
+  for (auto i = 0ul; i < stepping.stepSizes.size(); ++i) {
     h[i] = stepping.stepSizes[i].value();
   }
 
@@ -285,13 +307,23 @@ Result<double> MultiEigenStepperSIMD<N, AA, BB, CC, DD>::step(
   const SimdScalar half_h = h * SimdScalar(0.5);
 
   // Second Runge-Kutta point
-  const SimdVector3 pos1 = pos + half_h * dir + h2 * 0.125 * sd.k1;
+  const SimdVector3 pos1 = pos + half_h * dir + h2 * SimdScalar(0.125) * sd.k1;
   sd.B_middle = getMultiField(stepping, pos1);
+  
+  ACTS_VERBOSE("Runge-Kutta k1:");
+  ACTS_VERBOSE("  x = " << sd.k1[0]);
+  ACTS_VERBOSE("  y = " << sd.k1[1]);
+  ACTS_VERBOSE("  z = " << sd.k1[2]);
 
   if (!stepping.extension.k2(state, MultiProxyStepper{}, navigator, sd.k2,
                              sd.B_middle, sd.kQoP, half_h, sd.k1)) {
     return EigenStepperError::StepInvalid;
   }
+  
+  ACTS_VERBOSE("Runge-Kutta k2:");
+  ACTS_VERBOSE("  x = " << sd.k2[0]);
+  ACTS_VERBOSE("  y = " << sd.k2[1]);
+  ACTS_VERBOSE("  z = " << sd.k2[2]);
 
   // check for nan
   for (int i = 0; i < 3; ++i) {
@@ -312,12 +344,22 @@ Result<double> MultiEigenStepperSIMD<N, AA, BB, CC, DD>::step(
   // Last Runge-Kutta point
   const SimdVector3 pos2 = pos + h * dir + h2 * 0.5 * sd.k3;
   sd.B_last = getMultiField(stepping, pos2);
+  
+  ACTS_VERBOSE("Runge-Kutta k3:");
+  ACTS_VERBOSE("  x = " << sd.k3[0]);
+  ACTS_VERBOSE("  y = " << sd.k3[1]);
+  ACTS_VERBOSE("  z = " << sd.k3[2]);
 
   if (!stepping.extension.k4(state, MultiProxyStepper{}, navigator, sd.k4,
                              sd.B_last, sd.kQoP, h, sd.k3)) {
     return EigenStepperError::StepInvalid;
   }
-
+  
+  ACTS_VERBOSE("Runge-Kutta k4:");
+  ACTS_VERBOSE("  x = " << sd.k4[0]);
+  ACTS_VERBOSE("  y = " << sd.k4[1]);
+  ACTS_VERBOSE("  z = " << sd.k4[2]);
+  
   // check for nan
   for (int i = 0; i < 3; ++i) {
     assert(!sd.k4[i].isNaN().any() && "k4 contains nan");