TrackFinding/src/TrackParamsEstimationAlgorithm.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/TrackFinding/TrackParamsEstimationAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/Seed.hpp"
0014 #include "Acts/Geometry/GeometryIdentifier.hpp"
0015 #include "Acts/Geometry/TrackingGeometry.hpp"
0016 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0017 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0020 #include "ActsExamples/EventData/SimSpacePoint.hpp"
0021 #include "ActsExamples/EventData/Track.hpp"
0022 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0023
0024 #include <cstddef>
0025 #include <optional>
0026 #include <ostream>
0027 #include <stdexcept>
0028 #include <utility>
0029 #include <vector>
0030
0031 namespace ActsExamples {
0032
0033 TrackParamsEstimationAlgorithm::TrackParamsEstimationAlgorithm(
0034     TrackParamsEstimationAlgorithm::Config cfg, Acts::Logging::Level lvl)
0035     : IAlgorithm("TrackParamsEstimationAlgorithm", lvl), m_cfg(std::move(cfg)) {
0036   if (m_cfg.inputSeeds.empty()) {
0037     throw std::invalid_argument("Missing seeds input collection");
0038   }
0039   if (m_cfg.outputTrackParameters.empty()) {
0040     throw std::invalid_argument("Missing track parameters output collection");
0041   }
0042   if (!m_cfg.trackingGeometry) {
0043     throw std::invalid_argument("Missing tracking geometry");
0044   }
0045   if (!m_cfg.magneticField) {
0046     throw std::invalid_argument("Missing magnetic field");
0047   }
0048
0049   m_inputSeeds.initialize(m_cfg.inputSeeds);
0050   m_inputTracks.maybeInitialize(m_cfg.inputProtoTracks);
0051
0052   m_outputTrackParameters.initialize(m_cfg.outputTrackParameters);
0053   m_outputSeeds.maybeInitialize(m_cfg.outputSeeds);
0054   m_outputTracks.maybeInitialize(m_cfg.outputProtoTracks);
0055 }
0056
0057 ProcessCode TrackParamsEstimationAlgorithm::execute(
0058     const AlgorithmContext& ctx) const {
0059   auto const& seeds = m_inputSeeds(ctx);
0060   ACTS_VERBOSE("Read " << seeds.size() << " seeds");
0061
0062   TrackParametersContainer trackParameters;
0063   trackParameters.reserve(seeds.size());
0064
0065   SimSeedContainer outputSeeds;
0066   if (m_outputSeeds.isInitialized()) {
0067     outputSeeds.reserve(seeds.size());
0068   }
0069
0070   const ProtoTrackContainer* inputTracks = nullptr;
0071   ProtoTrackContainer outputTracks;
0072   if (m_inputTracks.isInitialized() && m_outputTracks.isInitialized()) {
0073     const auto& inputTracksRef = m_inputTracks(ctx);
0074     if (seeds.size() != inputTracksRef.size()) {
0075       ACTS_FATAL("Inconsistent number of seeds and prototracks");
0076       return ProcessCode::ABORT;
0077     }
0078     inputTracks = &inputTracksRef;
0079     outputTracks.reserve(seeds.size());
0080   }
0081
0082   auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext);
0083
0084   IndexSourceLink::SurfaceAccessor surfaceAccessor{*m_cfg.trackingGeometry};
0085
0086   // Loop over all found seeds to estimate track parameters
0087   for (std::size_t iseed = 0; iseed < seeds.size(); ++iseed) {
0088     const auto& seed = seeds[iseed];
0089     // Get the bottom space point and its reference surface
0090     const auto& bottomSP = seed.sp().front();
0091     if (bottomSP->sourceLinks().empty()) {
0092       ACTS_WARNING("Missing source link in the space point");
0093       continue;
0094     }
0095     const auto& sourceLink = bottomSP->sourceLinks()[0];
0096     const Acts::Surface* surface = surfaceAccessor(sourceLink);
0097
0098     if (surface == nullptr) {
0099       ACTS_WARNING(
0100           "Surface from source link is not found in the tracking geometry");
0101       continue;
0102     }
0103
0104     // Get the magnetic field at the bottom space point
0105     auto fieldRes = m_cfg.magneticField->getField(
0106         {bottomSP->x(), bottomSP->y(), bottomSP->z()}, bCache);
0107     if (!fieldRes.ok()) {
0108       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0109       return ProcessCode::ABORT;
0110     }
0111     Acts::Vector3 field = *fieldRes;
0112
0113     if (field.norm() < m_cfg.bFieldMin) {
0114       ACTS_WARNING("Magnetic field at seed " << iseed << " is too small "
0115                                              << field.norm());
0116       continue;
0117     }
0118
0119     // Estimate the track parameters from seed
0120     const auto paramsResult = Acts::estimateTrackParamsFromSeed(
0121         ctx.geoContext, seed.sp(), *surface, field);
0122     if (!paramsResult.ok()) {
0123       ACTS_WARNING("Skip track because param estimation failed "
0124                    << paramsResult.error());
0125       continue;
0126     }
0127     const auto& params = *paramsResult;
0128
0129     Acts::EstimateTrackParamCovarianceConfig config{
0130         .initialSigmas =
0131             Eigen::Map<const Acts::BoundVector>{m_cfg.initialSigmas.data()},
0132         .initialSigmaPtRel = m_cfg.initialSigmaPtRel,
0133         .initialVarInflation = Eigen::Map<const Acts::BoundVector>{
0134             m_cfg.initialVarInflation.data()}};
0135
0136     Acts::BoundSquareMatrix cov = Acts::estimateTrackParamCovariance(
0137         config, params, bottomSP->t().has_value());
0138
0139     trackParameters.emplace_back(surface->getSharedPtr(), params, cov,
0140                                  m_cfg.particleHypothesis);
0141     if (m_outputSeeds.isInitialized()) {
0142       outputSeeds.push_back(seed);
0143     }
0144     if (m_outputTracks.isInitialized() && inputTracks != nullptr) {
0145       outputTracks.push_back(inputTracks->at(iseed));
0146     }
0147   }
0148
0149   ACTS_VERBOSE("Estimated " << trackParameters.size() << " track parameters");
0150
0151   m_outputTrackParameters(ctx, std::move(trackParameters));
0152   if (m_outputSeeds.isInitialized()) {
0153     m_outputSeeds(ctx, std::move(outputSeeds));
0154   }
0155
0156   if (m_outputTracks.isInitialized()) {
0157     m_outputTracks(ctx, std::move(outputTracks));
0158   }
0159
0160   return ProcessCode::SUCCESS;
0161 }
0162
0163 }  // namespace ActsExamples