Acts/Propagator/MultiEigenStepperLoop.ipp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Propagator/MultiEigenStepperLoop.hpp"
0010 #include "Acts/Propagator/MultiStepperError.hpp"
0011 #include "Acts/Utilities/Logger.hpp"
0012
0013 namespace Acts {
0014
0015 template <typename E, typename R>
0016 auto MultiEigenStepperLoop<E, R>::boundState(
0017     State& state, const Surface& surface, bool transportCov,
0018     const FreeToBoundCorrection& freeToBoundCorrection) const
0019     -> Result<BoundState> {
0020   assert(!state.components.empty());
0021
0022   std::vector<std::tuple<double, BoundVector, Covariance>> cmps;
0023   cmps.reserve(numberComponents(state));
0024   double accumulatedPathLength = 0.0;
0025
0026   for (auto i = 0ul; i < numberComponents(state); ++i) {
0027     auto& cmpState = state.components[i].state;
0028
0029     // Force the component to be on the surface
0030     // This needs to be done because of the `averageOnSurface`-option of the
0031     // `MultiStepperSurfaceReached`-Aborter, which can be configured to end the
0032     // propagation when the mean of all components reached the destination
0033     // surface. Thus, it is not garantueed that all states are actually
0034     // onSurface.
0035     cmpState.pars.template segment<3>(eFreePos0) =
0036         surface
0037             .intersect(state.options.geoContext,
0038                        cmpState.pars.template segment<3>(eFreePos0),
0039                        cmpState.pars.template segment<3>(eFreeDir0),
0040                        BoundaryTolerance::Infinite())
0041             .closest()
0042             .position();
0043
0044     auto bs = SingleStepper::boundState(cmpState, surface, transportCov,
0045                                         freeToBoundCorrection);
0046
0047     if (bs.ok()) {
0048       const auto& btp = std::get<BoundTrackParameters>(*bs);
0049       cmps.emplace_back(
0050           state.components[i].weight, btp.parameters(),
0051           btp.covariance().value_or(Acts::BoundSquareMatrix::Zero()));
0052       accumulatedPathLength +=
0053           std::get<double>(*bs) * state.components[i].weight;
0054     }
0055   }
0056
0057   if (cmps.empty()) {
0058     return MultiStepperError::AllComponentsConversionToBoundFailed;
0059   }
0060
0061   return BoundState{MultiComponentBoundTrackParameters(
0062                         surface.getSharedPtr(), cmps, state.particleHypothesis),
0063                     Jacobian::Zero(), accumulatedPathLength};
0064 }
0065
0066 template <typename E, typename R>
0067 auto MultiEigenStepperLoop<E, R>::curvilinearState(
0068     State& state, bool transportCov) const -> CurvilinearState {
0069   assert(!state.components.empty());
0070
0071   std::vector<std::tuple<double, Vector4, Vector3, double, BoundSquareMatrix>>
0072       cmps;
0073   cmps.reserve(numberComponents(state));
0074   double accumulatedPathLength = 0.0;
0075
0076   for (auto i = 0ul; i < numberComponents(state); ++i) {
0077     const auto [cp, jac, pl] = SingleStepper::curvilinearState(
0078         state.components[i].state, transportCov);
0079
0080     cmps.emplace_back(state.components[i].weight,
0081                       cp.fourPosition(state.options.geoContext), cp.direction(),
0082                       cp.qOverP(),
0083                       cp.covariance().value_or(BoundSquareMatrix::Zero()));
0084     accumulatedPathLength += state.components[i].weight * pl;
0085   }
0086
0087   return CurvilinearState{
0088       MultiComponentCurvilinearTrackParameters(cmps, state.particleHypothesis),
0089       Jacobian::Zero(), accumulatedPathLength};
0090 }
0091
0092 template <typename E, typename R>
0093 template <typename propagator_state_t, typename navigator_t>
0094 Result<double> MultiEigenStepperLoop<E, R>::step(
0095     propagator_state_t& state, const navigator_t& navigator) const {
0096   using Status = Acts::IntersectionStatus;
0097
0098   State& stepping = state.stepping;
0099   auto& components = stepping.components;
0100   const Logger& logger = *m_logger;
0101
0102   // Update step count
0103   stepping.steps++;
0104
0105   // Check if we abort because of m_stepLimitAfterFirstComponentOnSurface
0106   if (stepping.stepCounterAfterFirstComponentOnSurface) {
0107     (*stepping.stepCounterAfterFirstComponentOnSurface)++;
0108
0109     // If the limit is reached, remove all components which are not on a
0110     // surface, reweight the components, perform no step and return 0
0111     if (*stepping.stepCounterAfterFirstComponentOnSurface >=
0112         m_stepLimitAfterFirstComponentOnSurface) {
0113       for (auto& cmp : components) {
0114         if (cmp.status != Status::onSurface) {
0115           cmp.status = Status::unreachable;
0116         }
0117       }
0118
0119       ACTS_VERBOSE("Stepper performed "
0120                    << m_stepLimitAfterFirstComponentOnSurface
0121                    << " steps after the first component hit a surface.");
0122       ACTS_VERBOSE(
0123           "-> remove all components not on a surface, perform no step");
0124
0125       removeMissedComponents(stepping);
0126       reweightComponents(stepping);
0127
0128       ACTS_VERBOSE(components.size()
0129                    << " components left after removing missed components");
0130
0131       stepping.stepCounterAfterFirstComponentOnSurface.reset();
0132
0133       return 0.0;
0134     }
0135   }
0136
0137   // Flag indicating if we need to reweight in the end
0138   bool reweightNecessary = false;
0139
0140   // If at least one component is on a surface, we can remove all missed
0141   // components before the step. If not, we must keep them for the case that all
0142   // components miss and we need to retarget
0143   const auto cmpsOnSurface = std::count_if(
0144       components.cbegin(), components.cend(),
0145       [&](auto& cmp) { return cmp.status == IntersectionStatus::onSurface; });
0146
0147   if (cmpsOnSurface > 0) {
0148     removeMissedComponents(stepping);
0149     reweightNecessary = true;
0150   }
0151
0152   // Loop over all components and collect results in vector, write some
0153   // summary information to a stringstream
0154   SmallVector<std::optional<Result<double>>> results;
0155   double accumulatedPathLength = 0.0;
0156   std::size_t errorSteps = 0;
0157
0158   // Type of the proxy single propagation2 state
0159   using ThisSinglePropState =
0160       detail::SinglePropState<SingleState, decltype(state.navigation),
0161                               decltype(state.options),
0162                               decltype(state.geoContext)>;
0163
0164   // Lambda that performs the step for a component and returns false if the step
0165   // went ok and true if there was an error
0166   auto errorInStep = [&](auto& component) {
0167     if (component.status == Status::onSurface) {
0168       // We need to add these, so the propagation does not fail if we have only
0169       // components on surfaces and failing states
0170       results.emplace_back(std::nullopt);
0171       return false;
0172     }
0173
0174     ThisSinglePropState single_state(component.state, state.navigation,
0175                                      state.options, state.geoContext);
0176
0177     results.emplace_back(SingleStepper::step(single_state, navigator));
0178
0179     if (results.back()->ok()) {
0180       accumulatedPathLength += component.weight * results.back()->value();
0181       return false;
0182     } else {
0183       ++errorSteps;
0184       reweightNecessary = true;
0185       return true;
0186     }
0187   };
0188
0189   // Loop over components and remove errorous components
0190   components.erase(
0191       std::remove_if(components.begin(), components.end(), errorInStep),
0192       components.end());
0193
0194   // Reweight if necessary
0195   if (reweightNecessary) {
0196     reweightComponents(stepping);
0197   }
0198
0199   // Print the result vector to a string so we can log it
0200   auto summary = [](auto& result_vec) {
0201     std::stringstream ss;
0202     for (auto& optRes : result_vec) {
0203       if (!optRes) {
0204         ss << "on surface | ";
0205       } else if (optRes->ok()) {
0206         ss << optRes->value() << " | ";
0207       } else {
0208         ss << optRes->error() << " | ";
0209       }
0210     }
0211     auto str = ss.str();
0212     str.resize(str.size() - 3);
0213     return str;
0214   };
0215
0216   // Print the summary
0217   if (errorSteps == 0) {
0218     ACTS_VERBOSE("Performed steps: " << summary(results));
0219   } else {
0220     ACTS_WARNING("Performed steps with errors: " << summary(results));
0221   }
0222
0223   // Return error if there is no ok result
0224   if (components.empty()) {
0225     return MultiStepperError::AllComponentsSteppingError;
0226   }
0227
0228   // Invalidate the component status after each step
0229   for (auto& cmp : components) {
0230     cmp.status = Status::unreachable;
0231   }
0232
0233   // Return the weighted accumulated path length of all successful steps
0234   stepping.pathAccumulated += accumulatedPathLength;
0235   return accumulatedPathLength;
0236 }
0237
0238 }  // namespace Acts