Acts/Propagator/MultiEigenStepperLoop.ipp

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