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/Geometry/GeometryIdentifier.hpp"
0015 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Utilities/Logger.hpp"
0018 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0019 #include "ActsExamples/EventData/SpacePoint.hpp"
0020 #include "ActsExamples/EventData/Track.hpp"
0021 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0022
0023 #include <array>
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     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0035     : IAlgorithm("TrackParamsEstimationAlgorithm", std::move(logger)),
0036       m_cfg(cfg) {
0037   if (m_cfg.inputSeeds.empty()) {
0038     throw std::invalid_argument("Missing seeds input collection");
0039   }
0040   if (m_cfg.outputTrackParameters.empty()) {
0041     throw std::invalid_argument("Missing track parameters output collection");
0042   }
0043   if (!m_cfg.trackingGeometry) {
0044     throw std::invalid_argument("Missing tracking geometry");
0045   }
0046   if (!m_cfg.magneticField) {
0047     throw std::invalid_argument("Missing magnetic field");
0048   }
0049
0050   m_inputSeeds.initialize(m_cfg.inputSeeds);
0051   m_inputTracks.maybeInitialize(m_cfg.inputProtoTracks);
0052   m_inputParticleHypotheses.maybeInitialize(m_cfg.inputParticleHypotheses);
0053
0054   m_outputTrackParameters.initialize(m_cfg.outputTrackParameters);
0055   m_outputSeeds.maybeInitialize(m_cfg.outputSeeds);
0056   m_outputTracks.maybeInitialize(m_cfg.outputProtoTracks);
0057 }
0058
0059 ProcessCode TrackParamsEstimationAlgorithm::execute(
0060     const AlgorithmContext& ctx) const {
0061   auto const& seeds = m_inputSeeds(ctx);
0062   ACTS_VERBOSE("Read " << seeds.size() << " seeds");
0063
0064   TrackParametersContainer trackParameters;
0065   trackParameters.reserve(seeds.size());
0066
0067   SeedContainer outputSeeds;
0068   if (m_outputSeeds.isInitialized()) {
0069     outputSeeds.assignSpacePointContainer(seeds.spacePointContainer());
0070     outputSeeds.reserve(seeds.size());
0071   }
0072
0073   const ProtoTrackContainer* inputTracks = nullptr;
0074   ProtoTrackContainer outputTracks;
0075   if (m_inputTracks.isInitialized() && m_outputTracks.isInitialized()) {
0076     const auto& inputTracksRef = m_inputTracks(ctx);
0077     if (seeds.size() != inputTracksRef.size()) {
0078       ACTS_FATAL("Inconsistent number of seeds and proto tracks");
0079       return ProcessCode::ABORT;
0080     }
0081     inputTracks = &inputTracksRef;
0082     outputTracks.reserve(seeds.size());
0083   }
0084
0085   const std::vector<Acts::ParticleHypothesis>* inputParticleHypotheses =
0086       nullptr;
0087   if (m_inputParticleHypotheses.isInitialized()) {
0088     const auto& inputParticleHypothesesRef = m_inputParticleHypotheses(ctx);
0089     if (seeds.size() != inputParticleHypothesesRef.size()) {
0090       ACTS_FATAL("Inconsistent number of seeds and particle hypotheses");
0091       return ProcessCode::ABORT;
0092     }
0093     inputParticleHypotheses = &inputParticleHypothesesRef;
0094   }
0095
0096   auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext);
0097
0098   IndexSourceLink::SurfaceAccessor surfaceAccessor{*m_cfg.trackingGeometry};
0099
0100   const SpacePointContainer& spacePoints = seeds.spacePointContainer();
0101
0102   // Loop over all found seeds to estimate track parameters
0103   for (std::size_t iseed = 0; iseed < seeds.size(); ++iseed) {
0104     const auto& seed = seeds[iseed];
0105     if (seed.spacePoints().size() < 3) {
0106       ACTS_WARNING("Seed " << iseed << " has less than 3 space points, skip");
0107       continue;
0108     }
0109
0110     const std::optional<std::array<SpacePointIndex, 3>> selected =
0111         selectSeedSpacePoints(spacePoints, seed.spacePointIndices(),
0112                               m_cfg.spacePointSelection);
0113     if (!selected.has_value()) {
0114       ACTS_DEBUG("Seed " << iseed << " has no space point selection, skip");
0115       continue;
0116     }
0117
0118     // Get the bottom space point and its reference surface
0119     const ConstSpacePointProxy bottomSp = spacePoints.at((*selected)[0]);
0120     const ConstSpacePointProxy middleSp = spacePoints.at((*selected)[1]);
0121     const ConstSpacePointProxy topSp = spacePoints.at((*selected)[2]);
0122     if (bottomSp.sourceLinks().empty()) {
0123       ACTS_WARNING("Missing source link in the space point");
0124       continue;
0125     }
0126
0127     const Acts::Vector3 bottomSpVec{bottomSp.x(), bottomSp.y(), bottomSp.z()};
0128     const Acts::Vector3 middleSpVec{middleSp.x(), middleSp.y(), middleSp.z()};
0129     const Acts::Vector3 topSpVec{topSp.x(), topSp.y(), topSp.z()};
0130
0131     const Acts::SourceLink& bottomSourceLink = bottomSp.sourceLinks()[0];
0132     const Acts::Surface* bottomSurface = surfaceAccessor(bottomSourceLink);
0133     if (bottomSurface == nullptr) {
0134       ACTS_WARNING(
0135           "Surface from source link is not found in the tracking geometry");
0136       continue;
0137     }
0138
0139     // Get the magnetic field at the bottom space point
0140     const auto fieldRes = m_cfg.magneticField->getField(bottomSpVec, bCache);
0141     if (!fieldRes.ok()) {
0142       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0143       return ProcessCode::ABORT;
0144     }
0145     const Acts::Vector3& field = *fieldRes;
0146
0147     if (field.norm() < m_cfg.bFieldMin) {
0148       ACTS_WARNING("Magnetic field at seed " << iseed << " is too small "
0149                                              << field.norm());
0150       continue;
0151     }
0152
0153     // Estimate the track parameters from seed
0154     Acts::Result<Acts::BoundVector> boundParams =
0155         Acts::estimateTrackParamsFromSeed(
0156             ctx.recoGeoContext, *bottomSurface, bottomSpVec,
0157             std::isnan(bottomSp.time()) ? 0.0 : bottomSp.time(), middleSpVec,
0158             topSpVec, field);
0159     if (!boundParams.ok()) {
0160       ACTS_WARNING("Failed to estimate track parameters from seed: "
0161                    << boundParams.error().message());
0162       continue;
0163     }
0164
0165     Acts::EstimateTrackParamCovarianceConfig config{
0166         .initialSigmas =
0167             Eigen::Map<const Acts::BoundVector>{m_cfg.initialSigmas.data()},
0168         .initialSigmaQoverPt = m_cfg.initialSigmaQoverPt,
0169         .initialSigmaPtRel = m_cfg.initialSigmaPtRel,
0170         .initialVarInflation = Eigen::Map<const Acts::BoundVector>{
0171             m_cfg.initialVarInflation.data()}};
0172
0173     const Acts::BoundMatrix cov = Acts::estimateTrackParamCovariance(
0174         config, *boundParams, !std::isnan(bottomSp.time()));
0175
0176     const Acts::ParticleHypothesis hypothesis =
0177         inputParticleHypotheses != nullptr ? inputParticleHypotheses->at(iseed)
0178                                            : m_cfg.particleHypothesis;
0179
0180     const TrackParameters& trackParams = trackParameters.emplace_back(
0181         bottomSurface->getSharedPtr(), *boundParams, cov, hypothesis);
0182     ACTS_VERBOSE("Estimated track parameters: " << trackParams);
0183     if (m_outputSeeds.isInitialized()) {
0184       auto newSp = outputSeeds.createSeed();
0185       // TODO copy shorthand
0186       newSp.assignSpacePointIndices(seed.spacePointIndices());
0187       newSp.quality() = seed.quality();
0188       newSp.vertexZ() = seed.vertexZ();
0189     }
0190     if (m_outputTracks.isInitialized() && inputTracks != nullptr) {
0191       outputTracks.push_back(inputTracks->at(iseed));
0192     }
0193   }
0194
0195   ACTS_DEBUG("Estimated " << trackParameters.size() << " track parameters");
0196
0197   m_outputTrackParameters(ctx, std::move(trackParameters));
0198   if (m_outputSeeds.isInitialized()) {
0199     m_outputSeeds(ctx, std::move(outputSeeds));
0200   }
0201
0202   if (m_outputTracks.isInitialized()) {
0203     m_outputTracks(ctx, std::move(outputTracks));
0204   }
0205
0206   return ProcessCode::SUCCESS;
0207 }
0208
0209 }  // namespace ActsExamples