Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-08-28 08:11:33

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