File indexing completed on 2025-08-28 08:11:33
0001
0002
0003
0004
0005
0006
0007
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
0032
0033
0034
0035
0036
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
0104 state.steps++;
0105
0106
0107 if (state.stepCounterAfterFirstComponentOnSurface) {
0108 (*state.stepCounterAfterFirstComponentOnSurface)++;
0109
0110
0111
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
0139 bool reweightNecessary = false;
0140
0141
0142
0143
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
0154
0155 SmallVector<std::optional<Result<double>>> results;
0156 double accumulatedPathLength = 0.0;
0157 std::size_t errorSteps = 0;
0158
0159
0160
0161 auto errorInStep = [this, &results, propDir, material, &accumulatedPathLength,
0162 &errorSteps, &reweightNecessary](auto& component) {
0163 if (component.status == Status::onSurface) {
0164
0165
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
0184 components.erase(
0185 std::remove_if(components.begin(), components.end(), errorInStep),
0186 components.end());
0187
0188
0189 if (reweightNecessary) {
0190 reweightComponents(state);
0191 }
0192
0193
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
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
0218 if (components.empty()) {
0219 return MultiStepperError::AllComponentsSteppingError;
0220 }
0221
0222
0223 for (auto& cmp : components) {
0224 cmp.status = Status::unreachable;
0225 }
0226
0227
0228 state.pathAccumulated += accumulatedPathLength;
0229 return accumulatedPathLength;
0230 }
0231
0232 }