File indexing completed on 2026-09-08 08:38:33
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Common.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/BoundTrackParameters.hpp"
0014 #include "Acts/EventData/MultiTrajectory.hpp"
0015 #include "Acts/EventData/TrackContainerFrontendConcept.hpp"
0016 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0017 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0018 #include "Acts/Geometry/GeometryContext.hpp"
0019 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0020 #include "Acts/Propagator/ActorList.hpp"
0021 #include "Acts/Propagator/DirectNavigator.hpp"
0022 #include "Acts/Propagator/PropagatorOptions.hpp"
0023 #include "Acts/Propagator/StandardAborters.hpp"
0024 #include "Acts/Propagator/detail/LoopProtection.hpp"
0025 #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
0026 #include "Acts/Utilities/CalibrationContext.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Utilities/Result.hpp"
0029
0030 #include <memory>
0031 #include <vector>
0032
0033 namespace Acts::Experimental {
0034
0035
0036 template <typename traj_t>
0037 struct ReferenceTrajectoryBuilderOptions {
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047 ReferenceTrajectoryBuilderOptions(
0048 const GeometryContext& gctx, const MagneticFieldContext& mctx,
0049 const PropagatorPlainOptions& pOptions, const Surface* tSurface = nullptr,
0050 bool mScattering = true, bool eLoss = true,
0051 const FreeToBoundCorrection& freeToBoundCorrection_ =
0052 FreeToBoundCorrection(false))
0053 : geoContext(gctx),
0054 magFieldContext(mctx),
0055 propagatorPlainOptions(pOptions),
0056 referenceSurface(tSurface),
0057 multipleScattering(mScattering),
0058 energyLoss(eLoss),
0059 freeToBoundCorrection(freeToBoundCorrection_) {}
0060
0061
0062 std::reference_wrapper<const GeometryContext> geoContext;
0063
0064 std::reference_wrapper<const MagneticFieldContext> magFieldContext;
0065
0066
0067 PropagatorPlainOptions propagatorPlainOptions;
0068
0069
0070 const Surface* referenceSurface = nullptr;
0071
0072
0073 bool multipleScattering = true;
0074
0075
0076 bool energyLoss = true;
0077
0078
0079
0080 FreeToBoundCorrection freeToBoundCorrection;
0081 };
0082
0083
0084 template <typename traj_t>
0085 struct ReferenceTrajectoryBuilderResult {
0086
0087 traj_t* trajectory{nullptr};
0088
0089
0090 std::size_t lastTrackStateIndex = kTrackIndexInvalid;
0091
0092
0093
0094 std::optional<BoundTrackParameters> referenceParameters;
0095
0096
0097 bool finished = false;
0098
0099
0100 PathLimitReached pathLimitReached;
0101 };
0102
0103
0104 template <typename propagator_t, typename traj_t>
0105 class ReferenceTrajectoryBuilder {
0106
0107 using NavigatorType = typename propagator_t::Navigator;
0108
0109
0110 static constexpr bool isDirectNavigator =
0111 std::is_same_v<NavigatorType, DirectNavigator>;
0112
0113
0114 using ResultType = ReferenceTrajectoryBuilderResult<traj_t>;
0115
0116 public:
0117
0118 using Options = ReferenceTrajectoryBuilderOptions<traj_t>;
0119
0120
0121 using TrackStateProxy = typename traj_t::TrackStateProxy;
0122
0123
0124 using ConstTrackStateProxy = typename traj_t::ConstTrackStateProxy;
0125
0126
0127
0128
0129 explicit ReferenceTrajectoryBuilder(
0130 propagator_t pPropagator,
0131 std::unique_ptr<const Logger> _logger =
0132 getDefaultLogger("ReferenceTrajectoryBuilder", Logging::INFO))
0133 : m_propagator(std::move(pPropagator)),
0134 m_logger{std::move(_logger)},
0135 m_actorLogger{m_logger->cloneWithSuffix("Actor")} {}
0136
0137 private:
0138
0139 propagator_t m_propagator;
0140
0141
0142 std::unique_ptr<const Logger> m_logger;
0143
0144 std::unique_ptr<const Logger> m_actorLogger;
0145
0146
0147 const Logger& logger() const { return *m_logger; }
0148
0149 class Actor {
0150 public:
0151
0152 using result_type = ResultType;
0153
0154
0155 SurfaceReached targetReached{std::numeric_limits<double>::lowest()};
0156
0157
0158 bool multipleScattering = true;
0159
0160
0161 bool energyLoss = true;
0162
0163
0164
0165 FreeToBoundCorrection freeToBoundCorrection;
0166
0167
0168 EndOfWorldReached endOfWorldReached;
0169
0170
0171 VolumeConstraintAborter volumeConstraintAborter;
0172
0173
0174 const Logger* actorLogger{nullptr};
0175
0176
0177 const Logger& logger() const { return *actorLogger; }
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189 template <typename propagator_state_t, typename stepper_t,
0190 typename navigator_t>
0191 Result<void> act(propagator_state_t& state, const stepper_t& stepper,
0192 const navigator_t& navigator, result_type& result,
0193 const Logger& ) const {
0194 assert(result.trajectory && "No MultiTrajectory set");
0195
0196 if (result.finished) {
0197 return Result<void>::success();
0198 }
0199
0200 ACTS_VERBOSE("ReferenceTrajectory step at pos: "
0201 << stepper.position(state.stepping).transpose()
0202 << " dir: " << stepper.direction(state.stepping).transpose()
0203 << " momentum: "
0204 << stepper.absoluteMomentum(state.stepping));
0205
0206 if (result.pathLimitReached.internalLimit ==
0207 std::numeric_limits<double>::max()) {
0208 detail::setupLoopProtection(state, stepper, result.pathLimitReached,
0209 true, logger());
0210 }
0211
0212 if (const Surface* surface = navigator.currentSurface(state.navigation);
0213 surface != nullptr) {
0214 ACTS_VERBOSE("Handle Surface " << surface->geometryId() << " "
0215 << state.options.direction);
0216 auto res = handleSurface(*surface, state, stepper, navigator, result);
0217 if (!res.ok()) {
0218 ACTS_DEBUG("Error in " << state.options.direction
0219 << " filter: " << res.error());
0220 return res.error();
0221 }
0222 }
0223
0224 const bool isEndOfWorldReached =
0225 endOfWorldReached.checkAbort(state, stepper, navigator, logger());
0226 const bool isVolumeConstraintReached = volumeConstraintAborter.checkAbort(
0227 state, stepper, navigator, logger());
0228 const bool isPathLimitReached = result.pathLimitReached.checkAbort(
0229 state, stepper, navigator, logger());
0230 const bool isTargetReached =
0231 targetReached.checkAbort(state, stepper, navigator, logger());
0232 if (isEndOfWorldReached || isVolumeConstraintReached ||
0233 isPathLimitReached || isTargetReached) {
0234 ACTS_VERBOSE(
0235 "Finalizing reference trajectory: "
0236 << (isEndOfWorldReached ? "end of world reached; " : "")
0237 << (isVolumeConstraintReached ? "volume constraint reached; " : "")
0238 << (isPathLimitReached ? "path limit reached; " : "")
0239 << (isTargetReached ? "target surface reached; " : ""));
0240
0241 if (isTargetReached) {
0242 ACTS_VERBOSE("Setting parameters at target surface");
0243
0244 auto res = stepper.boundState(state.stepping, *targetReached.surface);
0245 if (!res.ok()) {
0246 ACTS_DEBUG("Error while acquiring bound state for target surface: "
0247 << res.error() << " " << res.error().message());
0248 return res.error();
0249 } else {
0250 const auto& [boundParams, jacobian, pathLength] = *res;
0251 result.referenceParameters = boundParams;
0252 }
0253 }
0254
0255 result.finished = true;
0256 }
0257
0258 return Result<void>::success();
0259 }
0260
0261 template <typename propagator_state_t, typename stepper_t,
0262 typename navigator_t>
0263 bool checkAbort(propagator_state_t& , const stepper_t& ,
0264 const navigator_t& , const result_type& result,
0265 const Logger& ) const {
0266 return result.finished;
0267 }
0268
0269 template <typename propagator_state_t, typename stepper_t,
0270 typename navigator_t>
0271 Result<void> handleSurface(const Surface& surface,
0272 propagator_state_t& state,
0273 const stepper_t& stepper,
0274 const navigator_t& navigator,
0275 result_type& result) const {
0276 stepper.transportCovarianceToBound(state.stepping, surface,
0277 freeToBoundCorrection);
0278
0279 const Result<detail::PointwiseMaterialEffects> materialInteractionPreRes =
0280 detail::performMaterialInteraction(
0281 state, stepper, surface,
0282 detail::determineMaterialUpdateMode(
0283 state, navigator, MaterialUpdateMode::PreUpdate),
0284 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0285 logger());
0286 if (!materialInteractionPreRes.ok()) {
0287 ACTS_DEBUG("Error during pre-update material interaction: "
0288 << materialInteractionPreRes.error());
0289 return materialInteractionPreRes.error();
0290 }
0291
0292 TrackStatePropMask mask =
0293 TrackStatePropMask::Predicted | TrackStatePropMask::Jacobian;
0294 TrackStateProxy trackStateProxy =
0295 result.trajectory->makeTrackState(mask, result.lastTrackStateIndex);
0296
0297 ConstTrackStateProxy trackStateProxyConst{trackStateProxy};
0298
0299 trackStateProxy.setReferenceSurface(surface.getSharedPtr());
0300 auto res = stepper.boundState(state.stepping, surface, false,
0301 freeToBoundCorrection);
0302 if (!res.ok()) {
0303 ACTS_DEBUG("Propagate to surface " << surface.geometryId()
0304 << " failed: " << res.error());
0305 return res.error();
0306 }
0307 const auto& [boundParams, jacobian, pathLength] = *res;
0308
0309 trackStateProxy.predicted() = boundParams.parameters();
0310 trackStateProxy.predictedCovariance() = state.stepping.cov;
0311 trackStateProxy.jacobian() = jacobian;
0312 trackStateProxy.pathLength() = pathLength;
0313
0314 auto typeFlags = trackStateProxy.typeFlags();
0315 typeFlags.setHasParameters();
0316 if (surface.hasMaterial()) {
0317 typeFlags.setHasMaterial();
0318 }
0319
0320 result.lastTrackStateIndex = trackStateProxy.index();
0321
0322 const Result<detail::PointwiseMaterialEffects>
0323 materialInteractionPostRes = detail::performMaterialInteraction(
0324 state, stepper, surface,
0325 detail::determineMaterialUpdateMode(
0326 state, navigator, MaterialUpdateMode::PostUpdate),
0327 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0328 logger());
0329 if (!materialInteractionPostRes.ok()) {
0330 ACTS_DEBUG("Error during post-update material interaction: "
0331 << materialInteractionPostRes.error());
0332 return materialInteractionPostRes.error();
0333 }
0334
0335 return Result<void>::success();
0336 }
0337 };
0338
0339 public:
0340
0341
0342
0343
0344
0345
0346
0347 template <TrackContainerFrontend track_container_t>
0348 Result<typename track_container_t::TrackProxy> build(
0349 const BoundTrackParameters& sParameters, const Options& actorOptions,
0350 track_container_t& trackContainer) const {
0351 auto propagatorOptions = makePropagatorOptions(
0352 actorOptions, nullptr, actorOptions.referenceSurface);
0353 return buildImpl(sParameters, propagatorOptions, trackContainer);
0354 }
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365 template <TrackContainerFrontend track_container_t>
0366 Result<typename track_container_t::TrackProxy> build(
0367 const BoundTrackParameters& sParameters, const Options& actorOptions,
0368 const std::vector<const Surface*>& sSequence,
0369 track_container_t& trackContainer) const
0370 requires isDirectNavigator
0371 {
0372 auto propagatorOptions = makePropagatorOptions(
0373 actorOptions, &sSequence, actorOptions.referenceSurface);
0374 return buildImpl(sParameters, propagatorOptions, trackContainer);
0375 }
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385 template <typename track_proxy_t, typename source_link_range_t>
0386 void attachSourceLinks(
0387 track_proxy_t trackProxy, const source_link_range_t& sourceLinkRange,
0388 const SourceLinkSurfaceAccessor& surfaceAccessor) const {
0389 const std::size_t nMeasurements = std::ranges::distance(sourceLinkRange);
0390
0391 ACTS_VERBOSE("Preparing " << nMeasurements << " input measurements");
0392 std::unordered_map<const Surface*, SourceLink> inputMeasurements;
0393 for (const SourceLink& sourceLink : sourceLinkRange) {
0394 const Surface* surface = surfaceAccessor(sourceLink);
0395 inputMeasurements.try_emplace(surface, sourceLink);
0396 }
0397
0398 for (auto trackState : trackProxy.trackStates()) {
0399 if (!trackState.hasReferenceSurface()) {
0400 continue;
0401 }
0402 const Surface& surface = trackState.referenceSurface();
0403
0404 if (!surface.isSensitive()) {
0405 continue;
0406 }
0407
0408 auto typeFlagsMap = trackState.typeFlags();
0409
0410 if (const auto it = inputMeasurements.find(&surface);
0411 it == inputMeasurements.end()) {
0412 typeFlagsMap.setIsHole();
0413 } else {
0414 SourceLink sourceLink = it->second;
0415
0416 trackState.setUncalibratedSourceLink(std::move(sourceLink));
0417 typeFlagsMap.setHasMeasurement();
0418 }
0419 }
0420 }
0421
0422
0423 using Calibrator =
0424 Delegate<void(const GeometryContext&, const CalibrationContext&,
0425 const SourceLink&, TrackStateProxy)>;
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435 template <typename track_proxy_t>
0436 void calibrateMeasurements(const GeometryContext& geoContext,
0437 const CalibrationContext& calibrationContext,
0438 track_proxy_t trackProxy,
0439 const Calibrator& calibrator) const {
0440 for (auto trackState : trackProxy.trackStates()) {
0441 if (!trackState.typeFlags().hasMeasurement()) {
0442 continue;
0443 }
0444
0445 trackState.addComponents(TrackStatePropMask::Calibrated);
0446
0447 const SourceLink& sourceLink = trackState.getUncalibratedSourceLink();
0448 calibrator(geoContext, calibrationContext, sourceLink, trackState);
0449 }
0450 }
0451
0452
0453 using Updater = Delegate<Result<void>(const GeometryContext&, TrackStateProxy,
0454 const Logger&)>;
0455
0456
0457
0458
0459
0460
0461
0462
0463 template <typename track_proxy_t>
0464 Result<void> filter(const GeometryContext& geoContext,
0465 track_proxy_t trackProxy, const Updater& updater) const {
0466 std::optional<TrackStateProxy> lastTrackState;
0467 BoundVector accumulatedBoundDeltas = BoundVector::Zero();
0468 BoundMatrix predictedCovariance = BoundMatrix::Zero();
0469
0470 for (auto trackState : trackProxy.trackStates()) {
0471 if (lastTrackState.has_value()) {
0472
0473
0474
0475 accumulatedBoundDeltas = trackState.jacobian() * accumulatedBoundDeltas;
0476 trackState.predicted() += accumulatedBoundDeltas;
0477
0478 predictedCovariance = trackState.jacobian() * predictedCovariance *
0479 trackState.jacobian().transpose();
0480
0481 {
0482 const FreeVector freeParams = transformBoundToFreeParameters(
0483 trackState.referenceSurface(), geoContext,
0484 trackState.predicted());
0485
0486 const Result<MaterialSlab> materialSlabRes =
0487 detail::evaluateMaterialSlab(
0488 geoContext, trackState.referenceSurface(),
0489 Direction::Forward(), freeParams.segment<3>(eFreePos0),
0490 freeParams.segment<3>(eFreeDir0),
0491 MaterialUpdateMode::PreUpdate);
0492 if (!materialSlabRes.ok()) {
0493 ACTS_DEBUG(
0494 "Error evaluating material slab: " << materialSlabRes.error());
0495 return materialSlabRes.error();
0496 }
0497 const MaterialSlab materialSlab = materialSlabRes.value();
0498
0499 const detail::PointwiseMaterialEffects materialEffects =
0500 detail::computeMaterialEffects(
0501 materialSlab, trackProxy.particleHypothesis(),
0502 freeParams.segment<3>(eFreeDir0), freeParams[eFreeQOverP],
0503 true, true, true);
0504
0505 predictedCovariance(eBoundPhi, eBoundPhi) +=
0506 materialEffects.variancePhi;
0507 predictedCovariance(eBoundTheta, eBoundTheta) +=
0508 materialEffects.varianceTheta;
0509 predictedCovariance(eBoundQOverP, eBoundQOverP) +=
0510 materialEffects.varianceQoverP;
0511 }
0512
0513 trackState.predictedCovariance() = predictedCovariance;
0514 }
0515
0516 if (!trackState.typeFlags().hasMeasurement()) {
0517 trackState.shareFrom(trackState, TrackStatePropMask::Predicted,
0518 TrackStatePropMask::Filtered);
0519 } else {
0520 trackState.addComponents(TrackStatePropMask::Filtered);
0521
0522 const Result<void> updateResult =
0523 updater(geoContext, trackState, logger());
0524 if (!updateResult.ok()) {
0525 ACTS_DEBUG("Error in filter: " << updateResult.error());
0526 return updateResult.error();
0527 }
0528 }
0529
0530 lastTrackState = trackState;
0531 accumulatedBoundDeltas += trackState.filtered() - trackState.predicted();
0532 predictedCovariance = trackState.filteredCovariance();
0533
0534 {
0535 const FreeVector freeParams = transformBoundToFreeParameters(
0536 trackState.referenceSurface(), geoContext, trackState.filtered());
0537
0538 const Result<MaterialSlab> materialSlabRes =
0539 detail::evaluateMaterialSlab(
0540 geoContext, trackState.referenceSurface(), Direction::Forward(),
0541 freeParams.segment<3>(eFreePos0),
0542 freeParams.segment<3>(eFreeDir0),
0543 MaterialUpdateMode::PostUpdate);
0544 if (!materialSlabRes.ok()) {
0545 ACTS_DEBUG(
0546 "Error evaluating material slab: " << materialSlabRes.error());
0547 return materialSlabRes.error();
0548 }
0549 const MaterialSlab materialSlab = materialSlabRes.value();
0550
0551 const detail::PointwiseMaterialEffects materialEffects =
0552 detail::computeMaterialEffects(
0553 materialSlab, trackProxy.particleHypothesis(),
0554 freeParams.segment<3>(eFreeDir0), freeParams[eFreeQOverP], true,
0555 true, true);
0556
0557 predictedCovariance(eBoundPhi, eBoundPhi) +=
0558 materialEffects.variancePhi;
0559 predictedCovariance(eBoundTheta, eBoundTheta) +=
0560 materialEffects.varianceTheta;
0561 predictedCovariance(eBoundQOverP, eBoundQOverP) +=
0562 materialEffects.varianceQoverP;
0563 }
0564 }
0565
0566 return Result<void>::success();
0567 }
0568
0569 private:
0570 auto makePropagatorOptions(const Options& actorOptions,
0571 const std::vector<const Surface*>* sSequence,
0572 const Surface* targetSurface) const {
0573 using Actors = ActorList<Actor>;
0574 using PropagatorOptions = typename propagator_t::template Options<Actors>;
0575
0576 PropagatorOptions propagatorOptions(actorOptions.geoContext,
0577 actorOptions.magFieldContext);
0578 propagatorOptions.setPlainOptions(actorOptions.propagatorPlainOptions);
0579
0580 if constexpr (!isDirectNavigator) {
0581 if (sSequence != nullptr) {
0582 for (const Surface* surface : *sSequence) {
0583 propagatorOptions.navigation.appendExternalSurface(*surface);
0584 }
0585 }
0586 } else {
0587 assert(sSequence != nullptr &&
0588 "DirectNavigator requires a surface sequence for "
0589 "ReferenceTrajectory");
0590 propagatorOptions.navigation.externalSurfaces = *sSequence;
0591 }
0592
0593 auto& actor = propagatorOptions.actorList.template get<Actor>();
0594 actor.targetReached.surface = targetSurface;
0595 actor.multipleScattering = actorOptions.multipleScattering;
0596 actor.energyLoss = actorOptions.energyLoss;
0597 actor.freeToBoundCorrection = actorOptions.freeToBoundCorrection;
0598 actor.actorLogger = m_actorLogger.get();
0599
0600 return propagatorOptions;
0601 }
0602
0603 template <typename propagator_options_t,
0604 TrackContainerFrontend track_container_t>
0605 auto buildImpl(const BoundTrackParameters& sParameters,
0606 const propagator_options_t& propagatorOptions,
0607 track_container_t& trackContainer) const
0608 -> Result<typename track_container_t::TrackProxy> {
0609 auto propagatorState = m_propagator.makeState(propagatorOptions);
0610
0611 auto propagatorInitResult =
0612 m_propagator.initialize(propagatorState, sParameters);
0613 if (!propagatorInitResult.ok()) {
0614 ACTS_DEBUG("Propagation initialization failed: "
0615 << propagatorInitResult.error());
0616 return propagatorInitResult.error();
0617 }
0618
0619 auto& actorResult = propagatorState.template get<ResultType>();
0620 actorResult.trajectory = &trackContainer.trackStateContainer();
0621
0622 auto result = m_propagator.propagate(propagatorState);
0623
0624 if (!result.ok()) {
0625 ACTS_DEBUG("Propagation failed: " << result.error());
0626 return result.error();
0627 }
0628
0629 auto track = trackContainer.makeTrack();
0630 track.tipIndex() = actorResult.lastTrackStateIndex;
0631 if (actorResult.referenceParameters.has_value()) {
0632 const auto& params = *actorResult.referenceParameters;
0633 track.parameters() = params.parameters();
0634 track.covariance() = params.covariance().value();
0635 track.setReferenceSurface(params.referenceSurface().getSharedPtr());
0636 }
0637
0638 track.linkForward();
0639
0640 return track;
0641 }
0642 };
0643
0644 }