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 <cstddef>
0024 #include <optional>
0025 #include <ostream>
0026 #include <stdexcept>
0027 #include <utility>
0028 #include <vector>
0029
0030 namespace ActsExamples {
0031
0032 TrackParamsEstimationAlgorithm::TrackParamsEstimationAlgorithm(
0033     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0034     : IAlgorithm("TrackParamsEstimationAlgorithm", std::move(logger)),
0035       m_cfg(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   m_inputParticleHypotheses.maybeInitialize(m_cfg.inputParticleHypotheses);
0052
0053   m_outputTrackParameters.initialize(m_cfg.outputTrackParameters);
0054   m_outputSeeds.maybeInitialize(m_cfg.outputSeeds);
0055   m_outputTracks.maybeInitialize(m_cfg.outputProtoTracks);
0056 }
0057
0058 ProcessCode TrackParamsEstimationAlgorithm::execute(
0059     const AlgorithmContext& ctx) const {
0060   auto const& seeds = m_inputSeeds(ctx);
0061   ACTS_VERBOSE("Read " << seeds.size() << " seeds");
0062
0063   TrackParametersContainer trackParameters;
0064   trackParameters.reserve(seeds.size());
0065
0066   SeedContainer outputSeeds;
0067   if (m_outputSeeds.isInitialized()) {
0068     outputSeeds.assignSpacePointContainer(seeds.spacePointContainer());
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 proto tracks");
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     if (seed.spacePoints().size() < 3) {
0103       ACTS_WARNING("Seed " << iseed << " has less than 3 space points");
0104       continue;
0105     } else if (seed.spacePoints().size() > 3) {
0106       ACTS_WARNING(
0107           "Seed "
0108           << iseed
0109           << " has more than 3 space points, only the first 3 will be used");
0110     }
0111
0112     // Get the bottom space point and its reference surface
0113     const ConstSpacePointProxy bottomSp = seed.spacePoints()[0];
0114     const ConstSpacePointProxy middleSp = seed.spacePoints()[1];
0115     const ConstSpacePointProxy topSp = seed.spacePoints()[2];
0116     if (bottomSp.sourceLinks().empty()) {
0117       ACTS_WARNING("Missing source link in the space point");
0118       continue;
0119     }
0120
0121     const Acts::Vector3 bottomSpVec{bottomSp.x(), bottomSp.y(), bottomSp.z()};
0122     const Acts::Vector3 middleSpVec{middleSp.x(), middleSp.y(), middleSp.z()};
0123     const Acts::Vector3 topSpVec{topSp.x(), topSp.y(), topSp.z()};
0124
0125     const Acts::SourceLink& bottomSourceLink = bottomSp.sourceLinks()[0];
0126     const Acts::Surface* bottomSurface = surfaceAccessor(bottomSourceLink);
0127     if (bottomSurface == nullptr) {
0128       ACTS_WARNING(
0129           "Surface from source link is not found in the tracking geometry");
0130       continue;
0131     }
0132
0133     // Get the magnetic field at the bottom space point
0134     const auto fieldRes = m_cfg.magneticField->getField(bottomSpVec, bCache);
0135     if (!fieldRes.ok()) {
0136       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0137       return ProcessCode::ABORT;
0138     }
0139     const Acts::Vector3& field = *fieldRes;
0140
0141     if (field.norm() < m_cfg.bFieldMin) {
0142       ACTS_WARNING("Magnetic field at seed " << iseed << " is too small "
0143                                              << field.norm());
0144       continue;
0145     }
0146
0147     // Estimate the track parameters from seed
0148     Acts::Result<Acts::BoundVector> boundParams =
0149         Acts::estimateTrackParamsFromSeed(
0150             ctx.geoContext, *bottomSurface, bottomSpVec,
0151             std::isnan(bottomSp.time()) ? 0.0 : bottomSp.time(), middleSpVec,
0152             topSpVec, field);
0153     if (!boundParams.ok()) {
0154       ACTS_WARNING("Failed to estimate track parameters from seed: "
0155                    << boundParams.error().message());
0156       continue;
0157     }
0158
0159     Acts::EstimateTrackParamCovarianceConfig config{
0160         .initialSigmas =
0161             Eigen::Map<const Acts::BoundVector>{m_cfg.initialSigmas.data()},
0162         .initialSigmaQoverPt = m_cfg.initialSigmaQoverPt,
0163         .initialSigmaPtRel = m_cfg.initialSigmaPtRel,
0164         .initialVarInflation = Eigen::Map<const Acts::BoundVector>{
0165             m_cfg.initialVarInflation.data()}};
0166
0167     const Acts::BoundMatrix cov = Acts::estimateTrackParamCovariance(
0168         config, *boundParams, !std::isnan(bottomSp.time()));
0169
0170     const Acts::ParticleHypothesis hypothesis =
0171         inputParticleHypotheses != nullptr ? inputParticleHypotheses->at(iseed)
0172                                            : m_cfg.particleHypothesis;
0173
0174     const TrackParameters& trackParams = trackParameters.emplace_back(
0175         bottomSurface->getSharedPtr(), *boundParams, cov, hypothesis);
0176     ACTS_VERBOSE("Estimated track parameters: " << trackParams);
0177     if (m_outputSeeds.isInitialized()) {
0178       auto newSp = outputSeeds.createSeed();
0179       // TODO copy shorthand
0180       newSp.assignSpacePointIndices(seed.spacePointIndices());
0181       newSp.quality() = seed.quality();
0182       newSp.vertexZ() = seed.vertexZ();
0183     }
0184     if (m_outputTracks.isInitialized() && inputTracks != nullptr) {
0185       outputTracks.push_back(inputTracks->at(iseed));
0186     }
0187   }
0188
0189   ACTS_DEBUG("Estimated " << trackParameters.size() << " track parameters");
0190
0191   m_outputTrackParameters(ctx, std::move(trackParameters));
0192   if (m_outputSeeds.isInitialized()) {
0193     m_outputSeeds(ctx, std::move(outputSeeds));
0194   }
0195
0196   if (m_outputTracks.isInitialized()) {
0197     m_outputTracks(ctx, std::move(outputTracks));
0198   }
0199
0200   return ProcessCode::SUCCESS;
0201 }
0202
0203 }  // namespace ActsExamples