TrackFinding/src/TrackParamsEstimationAlgorithm.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021-2024 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/ParticleHypothesis.hpp"
0014 #include "Acts/EventData/SourceLink.hpp"
0015 #include "Acts/Geometry/GeometryIdentifier.hpp"
0016 #include "Acts/Geometry/TrackingGeometry.hpp"
0017 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0018 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0019 #include "Acts/Seeding/Seed.hpp"
0020 #include "Acts/Surfaces/Surface.hpp"
0021 #include "Acts/Utilities/Result.hpp"
0022 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0023 #include "ActsExamples/EventData/SimSpacePoint.hpp"
0024 #include "ActsExamples/EventData/Track.hpp"
0025 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0026
0027 #include <cmath>
0028 #include <cstddef>
0029 #include <optional>
0030 #include <ostream>
0031 #include <stdexcept>
0032 #include <system_error>
0033 #include <utility>
0034 #include <vector>
0035
0036 namespace ActsExamples {
0037
0038 namespace {
0039
0040 Acts::BoundSquareMatrix makeInitialCovariance(
0041     const TrackParamsEstimationAlgorithm::Config& config,
0042     const Acts::BoundVector& params, const SimSpacePoint& sp) {
0043   Acts::BoundSquareMatrix result = Acts::BoundSquareMatrix::Zero();
0044
0045   for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
0046     double sigma = config.initialSigmas[i];
0047
0048     // Add momentum dependent uncertainties
0049     sigma +=
0050         config.initialSimgaQoverPCoefficients[i] * params[Acts::eBoundQOverP];
0051
0052     double var = sigma * sigma;
0053
0054     // Inflate the time uncertainty if no time measurement is available
0055     if (i == Acts::eBoundTime && !sp.t().has_value()) {
0056       var *= config.noTimeVarInflation;
0057     }
0058
0059     // Inflate the initial covariance
0060     var *= config.initialVarInflation[i];
0061
0062     result(i, i) = var;
0063   }
0064
0065   return result;
0066 }
0067
0068 }  // namespace
0069
0070 ActsExamples::TrackParamsEstimationAlgorithm::TrackParamsEstimationAlgorithm(
0071     ActsExamples::TrackParamsEstimationAlgorithm::Config cfg,
0072     Acts::Logging::Level lvl)
0073     : ActsExamples::IAlgorithm("TrackParamsEstimationAlgorithm", lvl),
0074       m_cfg(std::move(cfg)) {
0075   if (m_cfg.inputSeeds.empty()) {
0076     throw std::invalid_argument("Missing seeds input collection");
0077   }
0078   if (m_cfg.outputTrackParameters.empty()) {
0079     throw std::invalid_argument("Missing track parameters output collection");
0080   }
0081   if (!m_cfg.trackingGeometry) {
0082     throw std::invalid_argument("Missing tracking geometry");
0083   }
0084   if (!m_cfg.magneticField) {
0085     throw std::invalid_argument("Missing magnetic field");
0086   }
0087
0088   m_inputSeeds.initialize(m_cfg.inputSeeds);
0089   m_inputTracks.maybeInitialize(m_cfg.inputProtoTracks);
0090
0091   m_outputTrackParameters.initialize(m_cfg.outputTrackParameters);
0092   m_outputSeeds.maybeInitialize(m_cfg.outputSeeds);
0093   m_outputTracks.maybeInitialize(m_cfg.outputProtoTracks);
0094 }
0095
0096 ActsExamples::ProcessCode ActsExamples::TrackParamsEstimationAlgorithm::execute(
0097     const ActsExamples::AlgorithmContext& ctx) const {
0098   auto const& seeds = m_inputSeeds(ctx);
0099   ACTS_VERBOSE("Read " << seeds.size() << " seeds");
0100
0101   TrackParametersContainer trackParameters;
0102   trackParameters.reserve(seeds.size());
0103
0104   SimSeedContainer outputSeeds;
0105   if (m_outputSeeds.isInitialized()) {
0106     outputSeeds.reserve(seeds.size());
0107   }
0108
0109   const ProtoTrackContainer* inputTracks = nullptr;
0110   ProtoTrackContainer outputTracks;
0111   if (m_inputTracks.isInitialized() && m_outputTracks.isInitialized()) {
0112     const auto& inputTracksRef = m_inputTracks(ctx);
0113     if (seeds.size() != inputTracksRef.size()) {
0114       ACTS_FATAL("Inconsistent number of seeds and prototracks");
0115       return ProcessCode::ABORT;
0116     }
0117     inputTracks = &inputTracksRef;
0118     outputTracks.reserve(seeds.size());
0119   }
0120
0121   auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext);
0122
0123   IndexSourceLink::SurfaceAccessor surfaceAccessor{*m_cfg.trackingGeometry};
0124
0125   // Loop over all found seeds to estimate track parameters
0126   for (std::size_t iseed = 0; iseed < seeds.size(); ++iseed) {
0127     const auto& seed = seeds[iseed];
0128     // Get the bottom space point and its reference surface
0129     const auto bottomSP = seed.sp().front();
0130     if (bottomSP->sourceLinks().empty()) {
0131       ACTS_WARNING("Missing source link in the space point")
0132       continue;
0133     }
0134     const auto& sourceLink = bottomSP->sourceLinks()[0];
0135     const Acts::Surface* surface = surfaceAccessor(sourceLink);
0136
0137     if (surface == nullptr) {
0138       ACTS_WARNING(
0139           "Surface from source link is not found in the tracking geometry");
0140       continue;
0141     }
0142
0143     // Get the magnetic field at the bottom space point
0144     auto fieldRes = m_cfg.magneticField->getField(
0145         {bottomSP->x(), bottomSP->y(), bottomSP->z()}, bCache);
0146     if (!fieldRes.ok()) {
0147       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0148       return ProcessCode::ABORT;
0149     }
0150     Acts::Vector3 field = *fieldRes;
0151
0152     // Estimate the track parameters from seed
0153     auto optParams = Acts::estimateTrackParamsFromSeed(
0154         ctx.geoContext, seed.sp().begin(), seed.sp().end(), *surface, field,
0155         m_cfg.bFieldMin, logger());
0156     if (!optParams.has_value()) {
0157       ACTS_WARNING("Estimation of track parameters for seed " << iseed
0158                                                               << " failed.");
0159       continue;
0160     }
0161
0162     const auto& params = optParams.value();
0163
0164     Acts::BoundSquareMatrix cov =
0165         makeInitialCovariance(m_cfg, params, *bottomSP);
0166
0167     trackParameters.emplace_back(surface->getSharedPtr(), params, cov,
0168                                  m_cfg.particleHypothesis);
0169     if (m_outputSeeds.isInitialized()) {
0170       outputSeeds.push_back(seed);
0171     }
0172     if (m_outputTracks.isInitialized() && inputTracks != nullptr) {
0173       outputTracks.push_back(inputTracks->at(iseed));
0174     }
0175   }
0176
0177   ACTS_VERBOSE("Estimated " << trackParameters.size() << " track parameters");
0178
0179   m_outputTrackParameters(ctx, std::move(trackParameters));
0180   if (m_outputSeeds.isInitialized()) {
0181     m_outputSeeds(ctx, std::move(outputSeeds));
0182   }
0183
0184   if (m_outputTracks.isInitialized()) {
0185     m_outputTracks(ctx, std::move(outputTracks));
0186   }
0187
0188   return ProcessCode::SUCCESS;
0189 }
0190 }  // namespace ActsExamples