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