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