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