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