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