File indexing completed on 2026-07-16 08:07:34
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/EventData/MultiComponentTrackParameters.hpp"
0015 #include "Acts/EventData/TrackParameters.hpp"
0016 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0017 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0018 #include "Acts/Material/IVolumeMaterial.hpp"
0019 #include "Acts/Propagator/ConstrainedStep.hpp"
0020 #include "Acts/Propagator/StepperConcept.hpp"
0021 #include "Acts/Propagator/StepperOptions.hpp"
0022 #include "Acts/Propagator/StepperStatistics.hpp"
0023 #include "Acts/Propagator/detail/LoopStepperUtils.hpp"
0024 #include "Acts/Propagator/detail/SteppingHelper.hpp"
0025 #include "Acts/Surfaces/Surface.hpp"
0026 #include "Acts/Utilities/Intersection.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Utilities/Result.hpp"
0029
0030 #include <algorithm>
0031 #include <cmath>
0032 #include <cstddef>
0033 #include <limits>
0034 #include <sstream>
0035 #include <vector>
0036
0037 #include <boost/container/small_vector.hpp>
0038
0039 namespace Acts {
0040
0041 namespace detail {
0042
0043 struct MaxMomentumComponent {
0044 template <typename component_range_t>
0045 auto operator()(const component_range_t& cmps) const {
0046 return std::ranges::max_element(cmps, [&](const auto& a, const auto& b) {
0047 return std::abs(a.state.pars[eFreeQOverP]) >
0048 std::abs(b.state.pars[eFreeQOverP]);
0049 });
0050 }
0051 };
0052
0053 struct MaxWeightComponent {
0054 template <typename component_range_t>
0055 auto operator()(const component_range_t& cmps) {
0056 return std::ranges::max_element(cmps, [&](const auto& a, const auto& b) {
0057 return a.weight < b.weight;
0058 });
0059 }
0060 };
0061
0062 template <typename component_chooser_t>
0063 struct SingleComponentReducer {
0064 template <typename stepper_state_t>
0065 static Vector3 position(const stepper_state_t& s) {
0066 return component_chooser_t{}(s.components)
0067 ->state.pars.template segment<3>(eFreePos0);
0068 }
0069
0070 template <typename stepper_state_t>
0071 static Vector3 direction(const stepper_state_t& s) {
0072 return component_chooser_t{}(s.components)
0073 ->state.pars.template segment<3>(eFreeDir0);
0074 }
0075
0076 template <typename stepper_state_t>
0077 static double qOverP(const stepper_state_t& s) {
0078 const auto cmp = component_chooser_t{}(s.components);
0079 return cmp->state.pars[eFreeQOverP];
0080 }
0081
0082 template <typename stepper_state_t>
0083 static double absoluteMomentum(const stepper_state_t& s) {
0084 const auto cmp = component_chooser_t{}(s.components);
0085 return s.particleHypothesis.extractMomentum(cmp->state.pars[eFreeQOverP]);
0086 }
0087
0088 template <typename stepper_state_t>
0089 static Vector3 momentum(const stepper_state_t& s) {
0090 const auto cmp = component_chooser_t{}(s.components);
0091 return s.particleHypothesis.extractMomentum(cmp->state.pars[eFreeQOverP]) *
0092 cmp->state.pars.template segment<3>(eFreeDir0);
0093 }
0094
0095 template <typename stepper_state_t>
0096 static double charge(const stepper_state_t& s) {
0097 const auto cmp = component_chooser_t{}(s.components);
0098 return s.particleHypothesis.extractCharge(cmp->state.pars[eFreeQOverP]);
0099 }
0100
0101 template <typename stepper_state_t>
0102 static double time(const stepper_state_t& s) {
0103 return component_chooser_t{}(s.components)->state.pars[eFreeTime];
0104 }
0105
0106 template <typename stepper_state_t>
0107 static FreeVector pars(const stepper_state_t& s) {
0108 return component_chooser_t{}(s.components)->state.pars;
0109 }
0110
0111 template <typename stepper_state_t>
0112 static FreeVector cov(const stepper_state_t& s) {
0113 return component_chooser_t{}(s.components)->state.cov;
0114 }
0115 };
0116
0117 }
0118
0119 using MaxMomentumReducerLoop =
0120 detail::SingleComponentReducer<detail::MaxMomentumComponent>;
0121 using MaxWeightReducerLoop =
0122 detail::SingleComponentReducer<detail::MaxWeightComponent>;
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134 template <Concepts::SingleStepper single_stepper_t,
0135 typename component_reducer_t = MaxWeightReducerLoop>
0136 class MultiStepperLoop : public single_stepper_t {
0137
0138
0139 std::size_t m_stepLimitAfterFirstComponentOnSurface = 50;
0140
0141
0142 std::unique_ptr<const Acts::Logger> m_logger;
0143
0144
0145
0146 template <typename T>
0147 using SmallVector = boost::container::small_vector<T, 16>;
0148
0149 public:
0150
0151 using SingleStepper = single_stepper_t;
0152
0153
0154 using SingleOptions = typename SingleStepper::Options;
0155
0156
0157 using SingleState = typename SingleStepper::State;
0158
0159
0160 using SingleConfig = typename SingleStepper::Config;
0161
0162
0163 using typename SingleStepper::Covariance;
0164 using typename SingleStepper::Jacobian;
0165
0166
0167 using BoundState =
0168 std::tuple<MultiComponentBoundTrackParameters, Jacobian, double>;
0169
0170
0171 using Reducer = component_reducer_t;
0172
0173
0174 static constexpr int maxComponents = std::numeric_limits<int>::max();
0175
0176 struct Config : public SingleStepper::Config {
0177
0178
0179 std::size_t stepLimitAfterFirstComponentOnSurface = 50;
0180 };
0181
0182 struct Options : public SingleOptions {
0183 using SingleOptions::SingleOptions;
0184 };
0185
0186 struct State {
0187
0188 struct Component {
0189
0190 SingleState state;
0191
0192 double weight;
0193
0194 IntersectionStatus status;
0195
0196
0197
0198
0199
0200 Component(SingleState state_, double weight_, IntersectionStatus status_)
0201 : state(std::move(state_)), weight(weight_), status(status_) {}
0202 };
0203
0204 Options options;
0205
0206
0207 ParticleHypothesis particleHypothesis = ParticleHypothesis::pion();
0208
0209
0210 SmallVector<Component> components;
0211
0212 bool covTransport = false;
0213 double pathAccumulated = 0.;
0214 std::size_t steps = 0;
0215
0216
0217
0218 std::optional<std::size_t> stepCounterAfterFirstComponentOnSurface;
0219
0220
0221 StepperStatistics statistics;
0222
0223
0224
0225
0226
0227
0228 explicit State(const Options& optionsIn) : options(optionsIn) {}
0229 };
0230
0231
0232
0233
0234
0235 explicit MultiStepperLoop(std::shared_ptr<const MagneticFieldProvider> bField,
0236 std::unique_ptr<const Logger> logger =
0237 getDefaultLogger("GSF", Logging::INFO))
0238 : SingleStepper(std::move(bField)), m_logger(std::move(logger)) {}
0239
0240
0241
0242
0243 explicit MultiStepperLoop(const Config& config,
0244 std::unique_ptr<const Logger> logger =
0245 getDefaultLogger("MultiStepperLoop",
0246 Logging::INFO))
0247 : SingleStepper(config),
0248 m_stepLimitAfterFirstComponentOnSurface(
0249 config.stepLimitAfterFirstComponentOnSurface),
0250 m_logger(std::move(logger)) {}
0251
0252
0253
0254
0255 State makeState(const Options& options) const {
0256 State state(options);
0257 return state;
0258 }
0259
0260
0261
0262
0263 void initialize(State& state,
0264 const MultiComponentBoundTrackParameters& par) const {
0265 if (par.components().empty()) {
0266 throw std::invalid_argument(
0267 "Cannot construct MultiEigenStepperLoop::State with empty "
0268 "multi-component parameters");
0269 }
0270
0271 state.particleHypothesis = par.particleHypothesis();
0272
0273 const auto surface = par.referenceSurface().getSharedPtr();
0274
0275 for (auto i = 0ul; i < par.components().size(); ++i) {
0276 const auto& [weight, singlePars] = par[i];
0277 auto& cmp =
0278 state.components.emplace_back(SingleStepper::makeState(state.options),
0279 weight, IntersectionStatus::onSurface);
0280 SingleStepper::initialize(cmp.state, singlePars);
0281 }
0282
0283 if (std::get<2>(par.components().front())) {
0284 state.covTransport = true;
0285 }
0286 }
0287
0288
0289
0290
0291
0292 using ConstComponentProxy =
0293 detail::LoopComponentProxyBase<const typename State::Component,
0294 MultiStepperLoop>;
0295
0296
0297
0298
0299
0300 using ComponentProxy =
0301 detail::LoopComponentProxy<typename State::Component, MultiStepperLoop>;
0302
0303
0304
0305
0306
0307
0308
0309 auto componentIterable(State& state) const {
0310 struct Iterator {
0311 using difference_type [[maybe_unused]] = std::ptrdiff_t;
0312 using value_type [[maybe_unused]] = ComponentProxy;
0313 using reference [[maybe_unused]] = ComponentProxy;
0314 using pointer [[maybe_unused]] = void;
0315 using iterator_category [[maybe_unused]] = std::forward_iterator_tag;
0316
0317 typename decltype(state.components)::iterator it;
0318 const State& s;
0319
0320
0321 auto& operator++() { ++it; return *this; }
0322 auto operator==(const Iterator& other) const { return it == other.it; }
0323 auto operator*() const { return ComponentProxy(*it, s); }
0324
0325 };
0326
0327 struct Iterable {
0328 State& s;
0329
0330
0331 auto begin() { return Iterator{s.components.begin(), s}; }
0332 auto end() { return Iterator{s.components.end(), s}; }
0333
0334 };
0335
0336 return Iterable{state};
0337 }
0338
0339
0340
0341
0342
0343
0344
0345 auto constComponentIterable(const State& state) const {
0346 struct ConstIterator {
0347 using difference_type [[maybe_unused]] = std::ptrdiff_t;
0348 using value_type [[maybe_unused]] = ConstComponentProxy;
0349 using reference [[maybe_unused]] = ConstComponentProxy;
0350 using pointer [[maybe_unused]] = void;
0351 using iterator_category [[maybe_unused]] = std::forward_iterator_tag;
0352
0353 typename decltype(state.components)::const_iterator it;
0354 const State& s;
0355
0356
0357 auto& operator++() { ++it; return *this; }
0358 auto operator==(const ConstIterator& other) const { return it == other.it; }
0359 auto operator*() const { return ConstComponentProxy{*it}; }
0360
0361 };
0362
0363 struct Iterable {
0364 const State& s;
0365
0366
0367 auto begin() const { return ConstIterator{s.components.cbegin(), s}; }
0368 auto end() const { return ConstIterator{s.components.cend(), s}; }
0369
0370 };
0371
0372 return Iterable{state};
0373 }
0374
0375
0376
0377
0378
0379 std::size_t numberComponents(const State& state) const {
0380 return state.components.size();
0381 }
0382
0383
0384
0385
0386 void removeMissedComponents(State& state) const {
0387 auto [beg, end] =
0388 std::ranges::remove_if(state.components, [](const auto& cmp) {
0389 return cmp.status == IntersectionStatus::unreachable;
0390 });
0391
0392 state.components.erase(beg, end);
0393 }
0394
0395
0396
0397
0398 void reweightComponents(State& state) const {
0399 double sumOfWeights = 0.0;
0400 for (const auto& cmp : state.components) {
0401 sumOfWeights += cmp.weight;
0402 }
0403 for (auto& cmp : state.components) {
0404 cmp.weight /= sumOfWeights;
0405 }
0406 }
0407
0408
0409
0410
0411 void clearComponents(State& state) const { state.components.clear(); }
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426 Result<ComponentProxy> addComponent(State& state,
0427 const BoundTrackParameters& pars,
0428 double weight) const {
0429 auto& cmp =
0430 state.components.emplace_back(SingleStepper::makeState(state.options),
0431 weight, IntersectionStatus::onSurface);
0432 SingleStepper::initialize(cmp.state, pars);
0433
0434 return ComponentProxy{state.components.back(), state};
0435 }
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446 Result<Vector3> getField(State& state, const Vector3& pos) const {
0447 return SingleStepper::getField(state.components.front().state, pos);
0448 }
0449
0450
0451
0452
0453
0454 Vector3 position(const State& state) const {
0455 return Reducer::position(state);
0456 }
0457
0458
0459
0460
0461
0462 Vector3 direction(const State& state) const {
0463 return Reducer::direction(state);
0464 }
0465
0466
0467
0468
0469
0470 double qOverP(const State& state) const { return Reducer::qOverP(state); }
0471
0472
0473
0474
0475
0476 double absoluteMomentum(const State& state) const {
0477 return Reducer::absoluteMomentum(state);
0478 }
0479
0480
0481
0482
0483
0484 Vector3 momentum(const State& state) const {
0485 return Reducer::momentum(state);
0486 }
0487
0488
0489
0490
0491
0492 double charge(const State& state) const { return Reducer::charge(state); }
0493
0494
0495
0496
0497
0498 ParticleHypothesis particleHypothesis(const State& state) const {
0499 return state.particleHypothesis;
0500 }
0501
0502
0503
0504
0505
0506 double time(const State& state) const { return Reducer::time(state); }
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522 IntersectionStatus updateSurfaceStatus(
0523 State& state, const Surface& surface, std::uint8_t index,
0524 Direction navDir, const BoundaryTolerance& boundaryTolerance,
0525 double surfaceTolerance, ConstrainedStep::Type stype,
0526 const Logger& logger = getDummyLogger()) const {
0527 using Status = IntersectionStatus;
0528
0529 std::array<int, 3> counts = {0, 0, 0};
0530
0531 for (auto& component : state.components) {
0532 component.status = detail::updateSingleSurfaceStatus<SingleStepper>(
0533 *this, component.state, surface, index, navDir, boundaryTolerance,
0534 surfaceTolerance, stype, logger);
0535 ++counts[static_cast<std::size_t>(component.status)];
0536 }
0537
0538
0539
0540
0541 if (counts[static_cast<std::size_t>(Status::onSurface)] > 0) {
0542 removeMissedComponents(state);
0543 reweightComponents(state);
0544 }
0545
0546 ACTS_VERBOSE("Component status wrt "
0547 << surface.geometryId() << " at {"
0548 << surface.center(state.options.geoContext).transpose()
0549 << "}:\t" << [&]() {
0550 std::stringstream ss;
0551 for (auto& component : state.components) {
0552 ss << component.status << "\t";
0553 }
0554 return ss.str();
0555 }());
0556
0557
0558
0559 if (!state.stepCounterAfterFirstComponentOnSurface &&
0560 counts[static_cast<std::size_t>(Status::onSurface)] > 0 &&
0561 counts[static_cast<std::size_t>(Status::reachable)] > 0) {
0562 state.stepCounterAfterFirstComponentOnSurface = 0;
0563 ACTS_VERBOSE("started stepCounterAfterFirstComponentOnSurface");
0564 }
0565
0566
0567
0568
0569 if (state.stepCounterAfterFirstComponentOnSurface &&
0570 counts[static_cast<std::size_t>(Status::onSurface)] == 0) {
0571 state.stepCounterAfterFirstComponentOnSurface.reset();
0572 ACTS_VERBOSE("switch off stepCounterAfterFirstComponentOnSurface");
0573 }
0574
0575
0576
0577
0578 if (counts[static_cast<std::size_t>(Status::reachable)] > 0) {
0579 return Status::reachable;
0580 } else if (counts[static_cast<std::size_t>(Status::onSurface)] > 0) {
0581 state.stepCounterAfterFirstComponentOnSurface.reset();
0582 return Status::onSurface;
0583 } else {
0584 return Status::unreachable;
0585 }
0586 }
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599 template <typename object_intersection_t>
0600 void updateStepSize(State& state, const object_intersection_t& oIntersection,
0601 Direction direction, ConstrainedStep::Type stype) const {
0602 const Surface& surface = *oIntersection.object();
0603
0604 for (auto& component : state.components) {
0605 auto intersection = surface.intersect(
0606 component.state.options.geoContext,
0607 SingleStepper::position(component.state),
0608 direction * SingleStepper::direction(component.state),
0609 BoundaryTolerance::None())[oIntersection.index()];
0610
0611 SingleStepper::updateStepSize(component.state, intersection, direction,
0612 stype);
0613 }
0614 }
0615
0616
0617
0618
0619
0620
0621 void updateStepSize(State& state, double stepSize,
0622 ConstrainedStep::Type stype) const {
0623 for (auto& component : state.components) {
0624 SingleStepper::updateStepSize(component.state, stepSize, stype);
0625 }
0626 }
0627
0628
0629
0630
0631
0632
0633
0634
0635
0636 double getStepSize(const State& state, ConstrainedStep::Type stype) const {
0637 return std::ranges::min_element(
0638 state.components,
0639 [=](const auto& a, const auto& b) {
0640 return std::abs(a.state.stepSize.value(stype)) <
0641 std::abs(b.state.stepSize.value(stype));
0642 })
0643 ->state.stepSize.value(stype);
0644 }
0645
0646
0647
0648
0649
0650 void releaseStepSize(State& state, ConstrainedStep::Type stype) const {
0651 for (auto& component : state.components) {
0652 SingleStepper::releaseStepSize(component.state, stype);
0653 }
0654 }
0655
0656
0657
0658
0659
0660 std::string outputStepSize(const State& state) const {
0661 std::stringstream ss;
0662 for (const auto& component : state.components) {
0663 ss << component.state.stepSize.toString() << " || ";
0664 }
0665
0666 return ss.str();
0667 }
0668
0669
0670
0671
0672
0673
0674
0675
0676
0677
0678
0679
0680
0681
0682
0683
0684
0685
0686
0687
0688
0689 Result<BoundState> boundState(
0690 State& state, const Surface& surface, bool transportCov = true,
0691 const FreeToBoundCorrection& freeToBoundCorrection =
0692 FreeToBoundCorrection(false)) const;
0693
0694
0695
0696
0697
0698
0699
0700
0701 bool prepareCurvilinearState(State& state) const {
0702 static_cast<void>(state);
0703 return true;
0704 }
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720 BoundState curvilinearState(State& state, bool transportCov = true) const;
0721
0722
0723
0724
0725
0726
0727 void transportCovarianceToCurvilinear(State& state) const {
0728 for (auto& component : state.components) {
0729 SingleStepper::transportCovarianceToCurvilinear(component.state);
0730 }
0731 }
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743
0744 void transportCovarianceToBound(
0745 State& state, const Surface& surface,
0746 const FreeToBoundCorrection& freeToBoundCorrection =
0747 FreeToBoundCorrection(false)) const {
0748 for (auto& component : state.components) {
0749 SingleStepper::transportCovarianceToBound(component.state, surface,
0750 freeToBoundCorrection);
0751 }
0752 }
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764 Result<double> step(State& state, Direction propDir,
0765 const IVolumeMaterial* material) const;
0766 };
0767
0768 }
0769
0770 #include "MultiStepperLoop.ipp"