TrackFitting/src/RefittingAlgorithm.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/TrackFitting/RefittingAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/EventData/BoundTrackParameters.hpp"
0015 #include "Acts/EventData/MultiTrajectory.hpp"
0016 #include "Acts/EventData/SourceLink.hpp"
0017 #include "Acts/EventData/TrackProxy.hpp"
0018 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0019 #include "Acts/EventData/VectorTrackContainer.hpp"
0020 #include "Acts/Surfaces/PerigeeSurface.hpp"
0021 #include "Acts/Surfaces/Surface.hpp"
0022 #include "Acts/Utilities/Result.hpp"
0023 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0024 #include "ActsExamples/TrackFitting/RefittingCalibrator.hpp"
0025 #include "ActsExamples/TrackFitting/TrackFitterFunction.hpp"
0026
0027 #include <algorithm>
0028 #include <array>
0029 #include <cstdint>
0030 #include <memory>
0031 #include <optional>
0032 #include <ostream>
0033 #include <stdexcept>
0034 #include <system_error>
0035 #include <utility>
0036 #include <vector>
0037
0038 namespace ActsExamples {
0039
0040 RefittingAlgorithm::RefittingAlgorithm(
0041     Config config, std::unique_ptr<const Acts::Logger> logger)
0042     : IAlgorithm("RefittingAlgorithm", std::move(logger)),
0043       m_cfg(std::move(config)) {
0044   if (m_cfg.inputTracks.empty()) {
0045     throw std::invalid_argument("Missing input tracks collection");
0046   }
0047   if (m_cfg.outputTracks.empty()) {
0048     throw std::invalid_argument("Missing output tracks collection");
0049   }
0050
0051   m_inputTracks.initialize(m_cfg.inputTracks);
0052   m_outputTracks.initialize(m_cfg.outputTracks);
0053 }
0054
0055 ProcessCode RefittingAlgorithm::execute(const AlgorithmContext& ctx) const {
0056   const auto& inputTracks = m_inputTracks(ctx);
0057
0058   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0059   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0060   TrackContainer tracks(trackContainer, trackStateContainer);
0061
0062   std::shared_ptr<const Acts::Surface> perigeeSurface =
0063       Acts::Surface::makeShared<Acts::PerigeeSurface>(
0064           Acts::Vector3{0., 0., 0.});
0065
0066   // The following code is only necessary if a beamspot constraint is in use but
0067   // unguarded for lifetime and simplicity reasons. The Core KF does not support
0068   // a beamspot constraint by itself but it can be achieved by injecting a fake
0069   // measurement on the perigee surface and processing it first in the fitter.
0070   // This can lead to problems if the beamspot constraint moves the initial
0071   // parameters too far and we start to miss following surfaces and
0072   // measurements. But in a refitting scenario with increased bounds and direct
0073   // navigation it *should* work OK.
0074
0075   auto beamSpotVectorTrackStateContainer =
0076       std::make_shared<Acts::VectorMultiTrajectory>();
0077   auto beamSpotTrackState = beamSpotVectorTrackStateContainer->makeTrackState();
0078
0079   const Acts::Vector2 beamSpotMeasValue{0., 0.};
0080
0081   beamSpotTrackState.setReferenceSurface(perigeeSurface);
0082
0083   if (m_cfg.beamSpotConstraint.has_value()) {
0084     ACTS_DEBUG("Using provided beam spot constraint matrix");
0085     beamSpotTrackState.allocateCalibrated(beamSpotMeasValue,
0086                                           m_cfg.beamSpotConstraint.value());
0087   } else {
0088     ACTS_DEBUG("No beam spot constraint provided, using zero matrix");
0089     beamSpotTrackState.allocateCalibrated(beamSpotMeasValue,
0090                                           Acts::SquareMatrix2::Zero());
0091   }
0092
0093   beamSpotTrackState.setProjectorSubspaceIndices(
0094       std::array<std::uint8_t, 2>{Acts::eBoundLoc0, Acts::eBoundLoc1});
0095
0096   Acts::SourceLink testSL{42};
0097   beamSpotTrackState.setUncalibratedSourceLink(std::move(testSL));
0098
0099   auto beamSpotConstVectorTrackStateContainer =
0100       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0101           std::move(*beamSpotVectorTrackStateContainer));
0102
0103   auto beamSpotConstTrackState =
0104       beamSpotConstVectorTrackStateContainer->getTrackState(
0105           beamSpotTrackState.index());
0106
0107   // Perform the fit for each input track
0108   std::vector<Acts::SourceLink> trackSourceLinks;
0109   std::vector<const Acts::Surface*> surfSequence;
0110   RefittingCalibrator calibrator;
0111
0112   auto itrack = 0ul;
0113   for (const auto& track : inputTracks) {
0114     // Check if you are not in picking mode
0115     ++itrack;
0116     if (m_cfg.pickTrack > -1 &&
0117         static_cast<std::size_t>(m_cfg.pickTrack) != itrack - 1) {
0118       continue;
0119     }
0120
0121     if (!track.hasReferenceSurface()) {
0122       ACTS_VERBOSE("Skip track " << itrack << ": missing ref surface");
0123       continue;
0124     }
0125
0126     // Reusing the track parameters on the beam spot perigee is only valid if
0127     // the track is already expressed with respect to it.
0128     const auto& refSurface = track.referenceSurface();
0129     const bool useBeamSpotConstraint = m_cfg.beamSpotConstraint.has_value();
0130     if (useBeamSpotConstraint && *perigeeSurface != refSurface) {
0131       ACTS_ERROR("Track " << itrack
0132                           << " is not parametrized on the beam spot perigee "
0133                              "surface, cannot apply the beam spot constraint");
0134       return ProcessCode::ABORT;
0135     }
0136
0137     // The fitter matches measurements by surface pointer, so with a beamspot
0138     // constraint the fit has to start on the very surface that carries it.
0139     const std::shared_ptr<const Acts::Surface> initialSurface =
0140         useBeamSpotConstraint ? perigeeSurface : refSurface.getSharedPtr();
0141
0142     TrackFitterFunction::GeneralFitterOptions options{
0143         ctx.recoGeoContext,
0144         ctx.magFieldContext,
0145         ctx.calibContext,
0146         perigeeSurface.get(),
0147         Acts::PropagatorPlainOptions(ctx.recoGeoContext, ctx.magFieldContext),
0148         true};
0149
0150     Acts::BoundTrackParameters initialParams(initialSurface, track.parameters(),
0151                                              track.covariance(),
0152                                              track.particleHypothesis());
0153
0154     if (initialParams.covariance()) {
0155       for (auto i = 0ul; i < m_cfg.initialVarInflation.size(); ++i) {
0156         (*initialParams.covariance())(i, i) *= m_cfg.initialVarInflation.at(i);
0157       }
0158     }
0159
0160     trackSourceLinks.clear();
0161     surfSequence.clear();
0162
0163     for (auto state : track.trackStatesReversed()) {
0164       surfSequence.push_back(&state.referenceSurface());
0165
0166       if (!state.hasCalibrated()) {
0167         continue;
0168       }
0169
0170       auto sl = RefittingCalibrator::RefittingSourceLink{state};
0171       trackSourceLinks.push_back(Acts::SourceLink{sl});
0172     }
0173
0174     if (surfSequence.empty()) {
0175       ACTS_DEBUG("Empty track " << itrack << " found.");
0176       continue;
0177     }
0178
0179     if (useBeamSpotConstraint) {
0180       RefittingCalibrator::RefittingSourceLink beamSpotSL{
0181           beamSpotConstTrackState};
0182       trackSourceLinks.emplace_back(Acts::SourceLink{beamSpotSL});
0183     }
0184
0185     std::ranges::reverse(surfSequence);
0186
0187     ACTS_VERBOSE("Initial parameters: "
0188                  << initialParams.fourPosition(ctx.recoGeoContext).transpose()
0189                  << " -> " << initialParams.direction().transpose());
0190
0191     ACTS_DEBUG("Invoke direct fitter for track " << itrack);
0192     auto result = (*m_cfg.fit)(trackSourceLinks, initialParams, options,
0193                                calibrator, surfSequence, tracks);
0194
0195     if (result.ok()) {
0196       // Get the fit output object
0197       const auto& refittedTrack = result.value();
0198       if (refittedTrack.hasReferenceSurface()) {
0199         ACTS_VERBOSE("Refitted parameters for track " << itrack);
0200         ACTS_VERBOSE("  " << track.parameters().transpose());
0201         ACTS_VERBOSE("Measurements: " << refittedTrack.nMeasurements());
0202         ACTS_VERBOSE("Outliers: " << refittedTrack.nOutliers());
0203       } else {
0204         ACTS_DEBUG("No refitted parameters for track " << itrack);
0205       }
0206     } else {
0207       ACTS_DEBUG("Fit failed for event "
0208                  << ctx.eventNumber << " track " << itrack << " with error: "
0209                  << result.error() << ", " << result.error().message());
0210     }
0211     ++itrack;
0212   }
0213
0214   ACTS_DEBUG("Fitted tracks: " << trackContainer->size());
0215
0216   if (logger().doPrint(Acts::Logging::DEBUG)) {
0217     std::stringstream ss;
0218     trackStateContainer->statistics().toStream(ss);
0219     ACTS_DEBUG(ss.str());
0220   }
0221
0222   ConstTrackContainer constTracks{
0223       std::make_shared<Acts::ConstVectorTrackContainer>(
0224           std::move(*trackContainer)),
0225       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0226           std::move(*trackStateContainer))};
0227
0228   m_outputTracks(ctx, std::move(constTracks));
0229   return ProcessCode::SUCCESS;
0230 }
0231
0232 }  // namespace ActsExamples