File indexing completed on 2025-11-05 08:54:01
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/MeasurementHelpers.hpp"
0014 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0015 #include "Acts/EventData/TrackContainerFrontendConcept.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/EventData/TrackProxyConcept.hpp"
0018 #include "Acts/EventData/TrackStateProxyConcept.hpp"
0019 #include "Acts/EventData/TrackStateType.hpp"
0020 #include "Acts/Geometry/GeometryContext.hpp"
0021 #include "Acts/Propagator/StandardAborters.hpp"
0022 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0023 #include "Acts/Surfaces/Surface.hpp"
0024 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0025 #include "Acts/Utilities/Logger.hpp"
0026 #include "Acts/Utilities/Result.hpp"
0027
0028 #include <utility>
0029
0030 namespace Acts {
0031
0032
0033 enum class TrackExtrapolationStrategy {
0034
0035 first,
0036
0037 last,
0038
0039
0040 firstOrLast,
0041 };
0042
0043
0044 enum class TrackExtrapolationError {
0045 CompatibleTrackStateNotFound = 1,
0046 ReferenceSurfaceUnreachable = 2,
0047 };
0048
0049
0050
0051
0052 std::error_code make_error_code(TrackExtrapolationError e);
0053
0054
0055
0056
0057 template <TrackProxyConcept track_proxy_t>
0058 Result<typename track_proxy_t::ConstTrackStateProxy> findFirstMeasurementState(
0059 const track_proxy_t &track) {
0060 using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0061
0062
0063
0064 std::optional<TrackStateProxy> firstMeasurementOpt;
0065
0066 for (const auto &trackState : track.trackStatesReversed()) {
0067 bool isMeasurement =
0068 trackState.typeFlags().test(TrackStateFlag::MeasurementFlag);
0069 bool isOutlier = trackState.typeFlags().test(TrackStateFlag::OutlierFlag);
0070
0071 if (isMeasurement && !isOutlier) {
0072 firstMeasurementOpt = trackState;
0073 }
0074 }
0075
0076 if (firstMeasurementOpt) {
0077 return Result<TrackStateProxy>::success(*firstMeasurementOpt);
0078 }
0079
0080 return Result<TrackStateProxy>::failure(
0081 TrackExtrapolationError::CompatibleTrackStateNotFound);
0082 }
0083
0084
0085
0086
0087 template <TrackProxyConcept track_proxy_t>
0088 Result<typename track_proxy_t::ConstTrackStateProxy> findLastMeasurementState(
0089 const track_proxy_t &track) {
0090 using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0091
0092 for (const auto &trackState : track.trackStatesReversed()) {
0093 bool isMeasurement =
0094 trackState.typeFlags().test(TrackStateFlag::MeasurementFlag);
0095 bool isOutlier = trackState.typeFlags().test(TrackStateFlag::OutlierFlag);
0096
0097 if (isMeasurement && !isOutlier) {
0098 return TrackStateProxy{trackState};
0099 }
0100 }
0101
0102 return Result<TrackStateProxy>::failure(
0103 TrackExtrapolationError::CompatibleTrackStateNotFound);
0104 }
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117 template <TrackProxyConcept track_proxy_t,
0118 typename smoother_t = GainMatrixSmoother>
0119 Result<void> smoothTrack(
0120 const GeometryContext &geoContext, track_proxy_t &track,
0121 const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO),
0122 smoother_t smoother = GainMatrixSmoother()) {
0123 auto &trackContainer = track.container();
0124 auto &trackStateContainer = trackContainer.trackStateContainer();
0125
0126 auto last = findLastMeasurementState(track);
0127 if (!last.ok()) {
0128 ACTS_ERROR("no last track state found");
0129 return last.error();
0130 }
0131
0132 auto smoothingResult =
0133 smoother(geoContext, trackStateContainer, last->index(), logger);
0134
0135 if (!smoothingResult.ok()) {
0136 ACTS_ERROR("Smoothing track " << track.index() << " failed with error "
0137 << smoothingResult.error());
0138 return smoothingResult.error();
0139 }
0140
0141 return Result<void>::success();
0142 }
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153 template <TrackContainerFrontend track_container_t>
0154 Result<void> smoothTracks(
0155 const GeometryContext &geoContext, const track_container_t &trackContainer,
0156 const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO)) {
0157 Result<void> result = Result<void>::success();
0158
0159 for (const auto &track : trackContainer) {
0160 auto smoothingResult = smoothTrack(geoContext, track, logger);
0161
0162
0163 if (!smoothingResult.ok() && result.ok()) {
0164 result = smoothingResult.error();
0165 }
0166 }
0167
0168 return result;
0169 }
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183 template <TrackProxyConcept track_proxy_t>
0184 Result<std::pair<typename track_proxy_t::ConstTrackStateProxy, double>>
0185 findTrackStateForExtrapolation(
0186 const GeometryContext &geoContext, const track_proxy_t &track,
0187 const Surface &referenceSurface, TrackExtrapolationStrategy strategy,
0188 const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0189 Logging::INFO)) {
0190 using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0191
0192 auto intersect = [&](const TrackStateProxy &state) -> Intersection3D {
0193 assert(state.hasSmoothed() || state.hasFiltered());
0194
0195 FreeVector freeVector;
0196 if (state.hasSmoothed()) {
0197 freeVector = MultiTrajectoryHelpers::freeSmoothed(geoContext, state);
0198 } else {
0199 freeVector = MultiTrajectoryHelpers::freeFiltered(geoContext, state);
0200 }
0201
0202 return referenceSurface
0203 .intersect(geoContext, freeVector.template segment<3>(eFreePos0),
0204 freeVector.template segment<3>(eFreeDir0),
0205 BoundaryTolerance::None(), s_onSurfaceTolerance)
0206 .closest();
0207 };
0208
0209 switch (strategy) {
0210 case TrackExtrapolationStrategy::first: {
0211 ACTS_VERBOSE("looking for first track state");
0212
0213 auto first = findFirstMeasurementState(track);
0214 if (!first.ok()) {
0215 ACTS_ERROR("no first track state found");
0216 return first.error();
0217 }
0218
0219 Intersection3D intersection = intersect(*first);
0220 if (!intersection.isValid()) {
0221 ACTS_ERROR("no intersection found");
0222 return Result<std::pair<TrackStateProxy, double>>::failure(
0223 TrackExtrapolationError::ReferenceSurfaceUnreachable);
0224 }
0225
0226 ACTS_VERBOSE("found intersection at " << intersection.pathLength());
0227 return std::pair(*first, intersection.pathLength());
0228 }
0229
0230 case TrackExtrapolationStrategy::last: {
0231 ACTS_VERBOSE("looking for last track state");
0232
0233 auto last = findLastMeasurementState(track);
0234 if (!last.ok()) {
0235 ACTS_ERROR("no last track state found");
0236 return last.error();
0237 }
0238
0239 Intersection3D intersection = intersect(*last);
0240 if (!intersection.isValid()) {
0241 ACTS_ERROR("no intersection found");
0242 return Result<std::pair<TrackStateProxy, double>>::failure(
0243 TrackExtrapolationError::ReferenceSurfaceUnreachable);
0244 }
0245
0246 ACTS_VERBOSE("found intersection at " << intersection.pathLength());
0247 return std::pair(*last, intersection.pathLength());
0248 }
0249
0250 case TrackExtrapolationStrategy::firstOrLast: {
0251 ACTS_VERBOSE("looking for first or last track state");
0252
0253 auto first = findFirstMeasurementState(track);
0254 if (!first.ok()) {
0255 ACTS_ERROR("no first track state found");
0256 return first.error();
0257 }
0258
0259 auto last = findLastMeasurementState(track);
0260 if (!last.ok()) {
0261 ACTS_ERROR("no last track state found");
0262 return last.error();
0263 }
0264
0265 Intersection3D intersectionFirst = intersect(*first);
0266 Intersection3D intersectionLast = intersect(*last);
0267
0268 double absDistanceFirst = std::abs(intersectionFirst.pathLength());
0269 double absDistanceLast = std::abs(intersectionLast.pathLength());
0270
0271 if (intersectionFirst.isValid() && absDistanceFirst <= absDistanceLast) {
0272 ACTS_VERBOSE("using first track state with intersection at "
0273 << intersectionFirst.pathLength());
0274 return std::pair(*first, intersectionFirst.pathLength());
0275 }
0276
0277 if (intersectionLast.isValid() && absDistanceLast <= absDistanceFirst) {
0278 ACTS_VERBOSE("using last track state with intersection at "
0279 << intersectionLast.pathLength());
0280 return std::pair(*last, intersectionLast.pathLength());
0281 }
0282
0283 ACTS_ERROR("no intersection found");
0284 return Result<std::pair<TrackStateProxy, double>>::failure(
0285 TrackExtrapolationError::ReferenceSurfaceUnreachable);
0286 }
0287 }
0288
0289
0290 return Result<std::pair<TrackStateProxy, double>>::failure(
0291 TrackExtrapolationError::CompatibleTrackStateNotFound);
0292 }
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308 template <TrackProxyConcept track_proxy_t, typename propagator_t,
0309 typename propagator_options_t>
0310 Result<void> extrapolateTrackToReferenceSurface(
0311 track_proxy_t &track, const Surface &referenceSurface,
0312 const propagator_t &propagator, propagator_options_t options,
0313 TrackExtrapolationStrategy strategy,
0314 const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0315 Logging::INFO)) {
0316 auto findResult = findTrackStateForExtrapolation(
0317 options.geoContext, track, referenceSurface, strategy, logger);
0318
0319 if (!findResult.ok()) {
0320 ACTS_ERROR("failed to find track state for extrapolation");
0321 return findResult.error();
0322 }
0323
0324 auto &[trackState, distance] = *findResult;
0325
0326 options.direction = Direction::fromScalarZeroAsPositive(distance);
0327
0328 BoundTrackParameters parameters = track.createParametersFromState(trackState);
0329 ACTS_VERBOSE("extrapolating track to reference surface at distance "
0330 << distance << " with direction " << options.direction
0331 << " with starting parameters " << parameters);
0332
0333 auto propagateResult =
0334 propagator.template propagate<BoundTrackParameters, propagator_options_t,
0335 ForcedSurfaceReached>(
0336 parameters, referenceSurface, options);
0337
0338 if (!propagateResult.ok()) {
0339 ACTS_ERROR("failed to extrapolate track: " << propagateResult.error());
0340 return propagateResult.error();
0341 }
0342
0343 track.setReferenceSurface(referenceSurface.getSharedPtr());
0344 track.parameters() = propagateResult->endParameters.value().parameters();
0345 track.covariance() =
0346 propagateResult->endParameters.value().covariance().value();
0347
0348 return Result<void>::success();
0349 }
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365 template <TrackContainerFrontend track_container_t, typename propagator_t,
0366 typename propagator_options_t>
0367 Result<void> extrapolateTracksToReferenceSurface(
0368 const track_container_t &trackContainer, const Surface &referenceSurface,
0369 const propagator_t &propagator, propagator_options_t options,
0370 TrackExtrapolationStrategy strategy,
0371 const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0372 Logging::INFO)) {
0373 Result<void> result = Result<void>::success();
0374
0375 for (const auto &track : trackContainer) {
0376 auto extrapolateResult = extrapolateTrackToReferenceSurface(
0377 track, referenceSurface, propagator, options, strategy, logger);
0378
0379
0380 if (!extrapolateResult.ok() && result.ok()) {
0381 result = extrapolateResult.error();
0382 }
0383 }
0384
0385 return result;
0386 }
0387
0388
0389
0390
0391
0392 template <TrackProxyConcept track_proxy_t>
0393 void calculateTrackQuantities(track_proxy_t track)
0394 requires(!track_proxy_t::ReadOnly)
0395 {
0396 track.chi2() = 0;
0397 track.nDoF() = 0;
0398
0399 track.nHoles() = 0;
0400 track.nMeasurements() = 0;
0401 track.nSharedHits() = 0;
0402 track.nOutliers() = 0;
0403
0404 for (const auto &trackState : track.trackStatesReversed()) {
0405 ConstTrackStateType typeFlags = trackState.typeFlags();
0406
0407 if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) {
0408 track.nHoles()++;
0409 } else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
0410 track.nOutliers()++;
0411 } else if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
0412 if (typeFlags.test(Acts::TrackStateFlag::SharedHitFlag)) {
0413 track.nSharedHits()++;
0414 }
0415 track.nMeasurements()++;
0416 track.chi2() += trackState.chi2();
0417 track.nDoF() += trackState.calibratedSize();
0418 }
0419 }
0420 }
0421
0422
0423
0424
0425
0426
0427
0428
0429 template <TrackProxyConcept track_proxy_t>
0430 void trimTrackFront(track_proxy_t track, bool trimHoles, bool trimOutliers,
0431 bool trimMaterial, bool trimOtherNoneMeasurement)
0432 requires(!track_proxy_t::ReadOnly)
0433 {
0434 using TrackStateProxy = typename track_proxy_t::TrackStateProxy;
0435
0436
0437
0438 std::optional<TrackStateProxy> front;
0439
0440 for (TrackStateProxy trackState : track.trackStatesReversed()) {
0441 TrackStateType typeFlags = trackState.typeFlags();
0442 bool isHole = typeFlags.test(TrackStateFlag::HoleFlag);
0443 bool isOutlier = typeFlags.test(TrackStateFlag::OutlierFlag);
0444 bool isMaterial = typeFlags.test(TrackStateFlag::MaterialFlag) &&
0445 !typeFlags.test(TrackStateFlag::MeasurementFlag);
0446 bool isOtherNoneMeasurement =
0447 !typeFlags.test(TrackStateFlag::MeasurementFlag) && !isHole &&
0448 !isOutlier && !isMaterial;
0449 if (trimHoles && isHole) {
0450 continue;
0451 }
0452 if (trimOutliers && isOutlier) {
0453 continue;
0454 }
0455 if (trimMaterial && isMaterial) {
0456 continue;
0457 }
0458 if (trimOtherNoneMeasurement && isOtherNoneMeasurement) {
0459 continue;
0460 }
0461
0462 front = trackState;
0463 }
0464
0465 if (front.has_value()) {
0466 front.value().previous() = TrackStateProxy::kInvalid;
0467 }
0468 }
0469
0470
0471
0472
0473
0474
0475
0476
0477 template <TrackProxyConcept track_proxy_t>
0478 void trimTrackBack(track_proxy_t track, bool trimHoles, bool trimOutliers,
0479 bool trimMaterial, bool trimOtherNoneMeasurement)
0480 requires(!track_proxy_t::ReadOnly)
0481 {
0482 using TrackStateProxy = typename track_proxy_t::TrackStateProxy;
0483
0484 std::optional<TrackStateProxy> back;
0485
0486 for (TrackStateProxy trackState : track.trackStatesReversed()) {
0487 back = trackState;
0488
0489 TrackStateType typeFlags = trackState.typeFlags();
0490 bool isHole = typeFlags.test(TrackStateFlag::HoleFlag);
0491 bool isOutlier = typeFlags.test(TrackStateFlag::OutlierFlag);
0492 bool isMaterial = typeFlags.test(TrackStateFlag::MaterialFlag) &&
0493 !typeFlags.test(TrackStateFlag::MeasurementFlag);
0494 bool isOtherNoneMeasurement =
0495 !typeFlags.test(TrackStateFlag::MeasurementFlag) && !isHole &&
0496 !isOutlier && !isMaterial;
0497 if (trimHoles && isHole) {
0498 continue;
0499 }
0500 if (trimOutliers && isOutlier) {
0501 continue;
0502 }
0503 if (trimMaterial && isMaterial) {
0504 continue;
0505 }
0506 if (trimOtherNoneMeasurement && isOtherNoneMeasurement) {
0507 continue;
0508 }
0509
0510 break;
0511 }
0512
0513 if (back.has_value()) {
0514 track.tipIndex() = back.value().index();
0515 }
0516 }
0517
0518
0519
0520
0521
0522
0523
0524
0525 template <TrackProxyConcept track_proxy_t>
0526 void trimTrack(track_proxy_t track, bool trimHoles, bool trimOutliers,
0527 bool trimMaterial, bool trimOtherNoneMeasurement)
0528 requires(!track_proxy_t::ReadOnly)
0529 {
0530 trimTrackFront(track, trimHoles, trimOutliers, trimMaterial,
0531 trimOtherNoneMeasurement);
0532 trimTrackBack(track, trimHoles, trimOutliers, trimMaterial,
0533 trimOtherNoneMeasurement);
0534 }
0535
0536
0537
0538
0539
0540
0541 template <std::size_t nMeasurementDim,
0542 TrackStateProxyConcept track_state_proxy_t>
0543 std::pair<ActsVector<nMeasurementDim>, ActsSquareMatrix<nMeasurementDim>>
0544 calculatePredictedResidual(track_state_proxy_t trackState) {
0545 using MeasurementVector = ActsVector<nMeasurementDim>;
0546 using MeasurementMatrix = ActsSquareMatrix<nMeasurementDim>;
0547
0548 if (!trackState.hasPredicted()) {
0549 throw std::invalid_argument("track state has no predicted parameters");
0550 }
0551 if (!trackState.hasCalibrated()) {
0552 throw std::invalid_argument("track state has no calibrated parameters");
0553 }
0554
0555 auto subspaceHelper =
0556 trackState.template projectorSubspaceHelper<nMeasurementDim>();
0557
0558 auto measurement = trackState.template calibrated<nMeasurementDim>();
0559 auto measurementCovariance =
0560 trackState.template calibratedCovariance<nMeasurementDim>();
0561 MeasurementVector predicted =
0562 subspaceHelper.projectVector(trackState.predicted());
0563 MeasurementMatrix predictedCovariance =
0564 subspaceHelper.projectMatrix(trackState.predictedCovariance());
0565
0566 MeasurementVector residual = measurement - predicted;
0567 MeasurementMatrix residualCovariance =
0568 measurementCovariance + predictedCovariance;
0569
0570 return {residual, residualCovariance};
0571 }
0572
0573
0574
0575
0576
0577
0578 template <std::size_t nMeasurementDim,
0579 TrackStateProxyConcept track_state_proxy_t>
0580 std::pair<ActsVector<nMeasurementDim>, ActsSquareMatrix<nMeasurementDim>>
0581 calculateFilteredResidual(track_state_proxy_t trackState) {
0582 using MeasurementVector = ActsVector<nMeasurementDim>;
0583 using MeasurementMatrix = ActsSquareMatrix<nMeasurementDim>;
0584
0585 if (!trackState.hasFiltered()) {
0586 throw std::invalid_argument("track state has no filtered parameters");
0587 }
0588 if (!trackState.hasCalibrated()) {
0589 throw std::invalid_argument("track state has no calibrated parameters");
0590 }
0591
0592 auto subspaceHelper =
0593 trackState.template projectorSubspaceHelper<nMeasurementDim>();
0594
0595 auto measurement = trackState.template calibrated<nMeasurementDim>();
0596 auto measurementCovariance =
0597 trackState.template calibratedCovariance<nMeasurementDim>();
0598 MeasurementVector filtered =
0599 subspaceHelper.projectVector(trackState.filtered());
0600 MeasurementMatrix filteredCovariance =
0601 subspaceHelper.projectMatrix(trackState.filteredCovariance());
0602
0603 MeasurementVector residual = measurement - filtered;
0604 MeasurementMatrix residualCovariance =
0605 measurementCovariance - filteredCovariance;
0606
0607 return {residual, residualCovariance};
0608 }
0609
0610
0611
0612
0613
0614
0615 template <std::size_t nMeasurementDim,
0616 TrackStateProxyConcept track_state_proxy_t>
0617 std::pair<ActsVector<nMeasurementDim>, ActsSquareMatrix<nMeasurementDim>>
0618 calculateSmoothedResidual(track_state_proxy_t trackState) {
0619 using MeasurementVector = ActsVector<nMeasurementDim>;
0620 using MeasurementMatrix = ActsSquareMatrix<nMeasurementDim>;
0621
0622 if (!trackState.hasSmoothed()) {
0623 throw std::invalid_argument("track state has no smoothed parameters");
0624 }
0625 if (!trackState.hasCalibrated()) {
0626 throw std::invalid_argument("track state has no calibrated parameters");
0627 }
0628
0629 auto subspaceHelper =
0630 trackState.template projectorSubspaceHelper<nMeasurementDim>();
0631
0632 auto measurement = trackState.template calibrated<nMeasurementDim>();
0633 auto measurementCovariance =
0634 trackState.template calibratedCovariance<nMeasurementDim>();
0635 MeasurementVector smoothed =
0636 subspaceHelper.projectVector(trackState.smoothed());
0637 MeasurementMatrix smoothedCovariance =
0638 subspaceHelper.projectMatrix(trackState.smoothedCovariance());
0639
0640 MeasurementVector residual = measurement - smoothed;
0641 MeasurementMatrix residualCovariance =
0642 measurementCovariance - smoothedCovariance;
0643
0644 return {residual, residualCovariance};
0645 }
0646
0647
0648
0649
0650
0651 template <TrackStateProxyConcept track_state_proxy_t>
0652 double calculatePredictedChi2(track_state_proxy_t trackState) {
0653 if (!trackState.hasPredicted()) {
0654 throw std::invalid_argument("track state has no predicted parameters");
0655 }
0656 if (!trackState.hasCalibrated()) {
0657 throw std::invalid_argument("track state has no calibrated parameters");
0658 }
0659
0660 return visit_measurement(
0661 trackState.calibratedSize(),
0662 [&]<std::size_t measdim>(
0663 std::integral_constant<std::size_t, measdim>) -> double {
0664 auto [residual, residualCovariance] =
0665 calculatePredictedResidual<measdim>(trackState);
0666
0667 return (residual.transpose() * residualCovariance.inverse() * residual)
0668 .eval()(0, 0);
0669 });
0670 }
0671
0672
0673
0674
0675
0676 template <TrackStateProxyConcept track_state_proxy_t>
0677 double calculateFilteredChi2(track_state_proxy_t trackState) {
0678 if (!trackState.hasFiltered()) {
0679 throw std::invalid_argument("track state has no filtered parameters");
0680 }
0681 if (!trackState.hasCalibrated()) {
0682 throw std::invalid_argument("track state has no calibrated parameters");
0683 }
0684
0685 return visit_measurement(
0686 trackState.calibratedSize(),
0687 [&]<std::size_t measdim>(
0688 std::integral_constant<std::size_t, measdim>) -> double {
0689 auto [residual, residualCovariance] =
0690 calculateFilteredResidual<measdim>(trackState);
0691
0692 return (residual.transpose() * residualCovariance.inverse() * residual)
0693 .eval()(0, 0);
0694 });
0695 }
0696
0697
0698
0699
0700
0701 template <TrackStateProxyConcept track_state_proxy_t>
0702 double calculateSmoothedChi2(track_state_proxy_t trackState) {
0703 if (!trackState.hasSmoothed()) {
0704 throw std::invalid_argument("track state has no smoothed parameters");
0705 }
0706 if (!trackState.hasCalibrated()) {
0707 throw std::invalid_argument("track state has no calibrated parameters");
0708 }
0709
0710 return visit_measurement(
0711 trackState.calibratedSize(),
0712 [&]<std::size_t measdim>(
0713 std::integral_constant<std::size_t, measdim>) -> double {
0714 auto [residual, residualCovariance] =
0715 calculateSmoothedResidual<measdim>(trackState);
0716
0717 return (residual.transpose() * residualCovariance.inverse() * residual)
0718 .eval()(0, 0);
0719 });
0720 }
0721
0722
0723
0724
0725
0726
0727
0728 template <TrackStateProxyConcept track_state_proxy_t>
0729 std::pair<BoundVector, BoundMatrix> calculateUnbiasedParametersCovariance(
0730 track_state_proxy_t trackState) {
0731 if (!trackState.hasSmoothed()) {
0732 throw std::invalid_argument("track state has no smoothed parameters");
0733 }
0734 if (!trackState.hasCalibrated()) {
0735 throw std::invalid_argument("track state has no calibrated parameters");
0736 }
0737
0738 return visit_measurement(
0739 trackState.calibratedSize(),
0740 [&]<std::size_t measdim>(std::integral_constant<std::size_t, measdim>) {
0741 FixedBoundSubspaceHelper<measdim> subspaceHelper =
0742 trackState.template projectorSubspaceHelper<measdim>();
0743
0744
0745 auto H = subspaceHelper.projector();
0746 auto s = trackState.smoothed();
0747 auto C = trackState.smoothedCovariance();
0748 auto m = trackState.template calibrated<measdim>();
0749 auto V = trackState.template calibratedCovariance<measdim>();
0750 auto K =
0751 (C * H.transpose() * (H * C * H.transpose() - V).inverse()).eval();
0752 BoundVector unbiasedParamsVec = s + K * (m - H * s);
0753 BoundMatrix unbiasedParamsCov = C - K * H * C;
0754 return std::make_pair(unbiasedParamsVec, unbiasedParamsCov);
0755 });
0756 }
0757
0758 }
0759
0760 namespace std {
0761
0762 template <>
0763 struct is_error_code_enum<Acts::TrackExtrapolationError> : std::true_type {};
0764 }