From 910e11fcf1c056c4dd5dc7984a39bdb28f5cf7a1 Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Thu, 4 Aug 2022 16:29:28 +0200 Subject: [PATCH] fix: Multi-Stepper stepping error in multiple components + refactoring (#1339) This PR does refactor the `MultiEigenStepperLoop` and fixes a situation, where a few stepping-erros could break the whole propagation. - **Fix:** If some components had stepping errors and all other components where on a surface, the propagation terminated - **Refactor:** - remove internal logger from `MultiEigenStepperLoop` since it is not really configurable from outside - improve `MultiEigenStepperLoop::step` (only construct summary-string in `VERBOSE` mode or on error, remove unnecessary loops and vectors, move things to lambdas, ...) **NOTE:** This does not fix the reason why stepping errors occur in some situations, this (hopefully) will be done in a subsequent PR) (cherry picked from commit 1b973d52cbd015bc0ef746f7a73b195ac5a8ae5e) --- .../Acts/Propagator/MultiEigenStepperLoop.hpp | 11 +- .../Acts/Propagator/MultiEigenStepperLoop.ipp | 115 ++++++++++-------- .../Acts/Propagator/MultiStepperError.hpp | 3 +- Core/src/Propagator/MultiStepperError.cpp | 2 + .../Core/Propagator/MultiStepperTests.cpp | 3 +- 5 files changed, 77 insertions(+), 57 deletions(-) diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp index b30ac589435..20cf12075e9 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp @@ -177,9 +177,6 @@ template , typename auctioneer_t = detail::VoidAuctioneer> class MultiEigenStepperLoop : public EigenStepper { - /// Allows logging in the member functions - const LoggerWrapper logger; - /// Limits the number of steps after at least one component reached the /// surface std::size_t m_stepLimitAfterFirstComponentOnSurface = 50; @@ -209,9 +206,8 @@ class MultiEigenStepperLoop static constexpr int maxComponents = std::numeric_limits::max(); /// Constructor from a magnetic field and a optionally provided Logger - MultiEigenStepperLoop(std::shared_ptr bField, - LoggerWrapper l = getDummyLogger()) - : EigenStepper(bField), logger(l) {} + MultiEigenStepperLoop(std::shared_ptr bField) + : EigenStepper(bField) {} struct State { /// The struct that stores the individual components @@ -608,9 +604,10 @@ class MultiEigenStepperLoop /// @param state [in,out] The stepping state (thread-local cache) /// @param surface [in] The surface provided /// @param bcheck [in] The boundary check for this status update + /// @param logger [in] A @c LoggerWrapper instance Intersection3D::Status updateSurfaceStatus( State& state, const Surface& surface, const BoundaryCheck& bcheck, - LoggerWrapper /*extLogger*/ = getDummyLogger()) const { + LoggerWrapper logger = getDummyLogger()) const { using Status = Intersection3D::Status; std::array counts = {0, 0, 0, 0}; diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp index 66ece75fc24..6611eac0b4c 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp @@ -35,16 +35,14 @@ auto MultiEigenStepperLoop::boundState( } } - if (failedBoundTransforms > 0) { - ACTS_ERROR("Multi component bound state: " - << failedBoundTransforms << " of " << numberComponents(state) - << " transforms failed"); - } - if (states.size() == 0) { return MultiStepperError::AllComponentsConversionToBoundFailed; } + if (failedBoundTransforms > 0) { + return MultiStepperError::SomeComponentsConversionToBoundFailed; + } + // TODO At ATLAS, the final parameters seem to be computed with the mode of // the mixture. At the moment, we use the mean of the mixture here, but // there should be done a comparison sometimes in the future. This could @@ -100,8 +98,29 @@ template template Result MultiEigenStepperLoop::step( propagator_state_t& state) const { + const auto& logger = state.options.logger; State& stepping = state.stepping; + // It is not possible to remove components from the vector, since the + // GSF actor relies on the fact that the ordering and number of + // components does not change + auto invalidateComponent = [](auto& cmp) { + cmp.status = Intersection3D::Status::missed; + cmp.weight = 0.0; + cmp.state.pars.template segment<3>(eFreeDir0) = Vector3::Zero(); + }; + + // Lambda for reweighting the components + auto reweight = [](auto& cmps) { + ActsScalar sumOfWeights = 0.0; + for (const auto& cmp : cmps) { + sumOfWeights += cmp.weight; + } + for (auto& cmp : cmps) { + cmp.weight /= sumOfWeights; + } + }; + // Update step count stepping.steps++; @@ -113,30 +132,19 @@ Result MultiEigenStepperLoop::step( // surface, reweight the components, perform no step and return 0 if (*stepping.stepCounterAfterFirstComponentOnSurface >= m_stepLimitAfterFirstComponentOnSurface) { - auto& cmps = stepping.components; - - // It is not possible to remove components from the vector, since the - // GSF actor relies on the fact that the ordering and number of - // components does not change - for (auto& cmp : cmps) { + for (auto& cmp : stepping.components) { if (cmp.status != Intersection3D::Status::onSurface) { - cmp.status = Intersection3D::Status::missed; - cmp.weight = 0.0; - cmp.state.pars.template segment<3>(eFreeDir0) = Vector3::Zero(); + invalidateComponent(cmp); } } - // Reweight - const auto sum_of_weights = std::accumulate( - cmps.begin(), cmps.end(), ActsScalar{0}, - [](auto sum, const auto& cmp) { return sum + cmp.weight; }); - for (auto& cmp : cmps) { - cmp.weight /= sum_of_weights; - } + reweight(stepping.components); + ACTS_VERBOSE("Stepper performed " + << m_stepLimitAfterFirstComponentOnSurface + << " after the first component hit a surface."); ACTS_VERBOSE( - "hit m_stepLimitAfterFirstComponentOnSurface, " - "perform no step"); + "-> remove all components not on a surface, perform no step"); stepping.stepCounterAfterFirstComponentOnSurface.reset(); @@ -146,15 +154,17 @@ Result MultiEigenStepperLoop::step( // Loop over all components and collect results in vector, write some // summary information to a stringstream - SmallVector> results; - std::stringstream ss; + SmallVector>> results; double accumulatedPathLength = 0.0; + std::size_t errorSteps = 0; for (auto& component : stepping.components) { // We must also propagate missed components for the case that all // components miss the target and we need to re-target if (component.status == Intersection3D::Status::onSurface) { - ss << "cmp skipped\t"; + // We need to add these, so the propagation does not fail if we have only + // components on surfaces and failing states + results.emplace_back(std::nullopt); continue; } @@ -165,41 +175,50 @@ Result MultiEigenStepperLoop::step( ThisSinglePropState single_state(component.state, state.navigation, state.options, state.geoContext); - results.push_back(SingleStepper::step(single_state)); + results.emplace_back(SingleStepper::step(single_state)); - if (results.back().ok()) { - accumulatedPathLength += component.weight * *results.back(); - ss << *results.back() << "\t"; + if (results.back()->ok()) { + accumulatedPathLength += component.weight * results.back()->value(); } else { - ss << "step error: " << results.back().error() << "\t"; + ++errorSteps; + invalidateComponent(component); } } - // Return no component was updated - if (results.empty()) { - return 0.0; + // Since we have invalidated some components, we need to reweight + if (errorSteps > 0) { + reweight(stepping.components); } - // Collect pointers to results which are ok, since Result is not copyable - SmallVector*> ok_results; - for (auto& res : results) { - if (res.ok()) { - ok_results.push_back(&res); + // Print the result vector to a string so we can log it + auto summary = [](auto& result_vec) { + std::stringstream ss; + for (auto& optRes : result_vec) { + if (not optRes) { + ss << "on surface | "; + } else if (optRes->ok()) { + ss << optRes->value() << " | "; + } else { + ss << optRes->error() << " | "; + } } + auto str = ss.str(); + str.resize(str.size() - 3); + return str; + }; + + // Print the summary + if (errorSteps == 0) { + ACTS_VERBOSE("Performed steps: " << summary(results)); + } else { + ACTS_WARNING("Performed steps with errors: " << summary(results)); } // Return error if there is no ok result - if (ok_results.empty()) { + if (errorSteps == results.size()) { return MultiStepperError::AllComponentsSteppingError; } - // Print the summary - if (ok_results.size() == results.size()) { - ACTS_VERBOSE("Performed steps: " << ss.str()); - } else { - ACTS_WARNING("Performed steps with errors: " << ss.str()); - } - // Return the weighted accumulated path length of all successful steps stepping.pathAccumulated += accumulatedPathLength; return accumulatedPathLength; diff --git a/Core/include/Acts/Propagator/MultiStepperError.hpp b/Core/include/Acts/Propagator/MultiStepperError.hpp index 66a1e1c2d7d..f0bd03b535f 100644 --- a/Core/include/Acts/Propagator/MultiStepperError.hpp +++ b/Core/include/Acts/Propagator/MultiStepperError.hpp @@ -18,7 +18,8 @@ enum class MultiStepperError { StateOfMultipleComponentsRequested = 2, AverageTrackLeftCurrentVolume = 3, AllComponentsSteppingError = 4, - AllComponentsConversionToBoundFailed = 5 + AllComponentsConversionToBoundFailed = 5, + SomeComponentsConversionToBoundFailed = 6 }; std::error_code make_error_code(Acts::MultiStepperError e); diff --git a/Core/src/Propagator/MultiStepperError.cpp b/Core/src/Propagator/MultiStepperError.cpp index 789cb75ae62..a1761d48d31 100644 --- a/Core/src/Propagator/MultiStepperError.cpp +++ b/Core/src/Propagator/MultiStepperError.cpp @@ -31,6 +31,8 @@ class MultiStepperErrorCategory : public std::error_category { return "Stepping error occured in all components"; case MultiStepperError::AllComponentsConversionToBoundFailed: return "The conversion to the bound state failed for all components"; + case MultiStepperError::SomeComponentsConversionToBoundFailed: + return "The conversion to the bound state failed for some components"; default: return "unknown"; } diff --git a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp index 517e83319f6..d3c1089729f 100644 --- a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp @@ -12,7 +12,7 @@ #include "Acts/MagneticField/NullBField.hpp" #include "Acts/Propagator/MultiEigenStepperLoop.hpp" #include "Acts/Propagator/MultiStepperAborters.hpp" -#include +#include "Acts/Propagator/Navigator.hpp" using namespace Acts; using namespace Acts::VectorHelpers; @@ -40,6 +40,7 @@ struct Options { double stepSizeCutOff = 0.0; std::size_t maxRungeKuttaStepTrials = 10; double mass = 1.0; + LoggerWrapper logger = Acts::getDummyLogger(); }; struct Navigation {};