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