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/Direction.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/EventData/BoundTrackParameters.hpp"
0015 #include "Acts/EventData/ParticleHypothesis.hpp"
0016 #include "Acts/Geometry/GeometryIdentifier.hpp"
0017 #include "Acts/Propagator/Propagator.hpp"
0018 #include "Acts/Propagator/PropagatorError.hpp"
0019 #include "Acts/Propagator/SympyStepper.hpp"
0020 #include "Acts/Propagator/VoidNavigator.hpp"
0021 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0022 #include "Acts/Surfaces/Surface.hpp"
0023 #include "Acts/Utilities/Intersection.hpp"
0024 #include "Acts/Utilities/Logger.hpp"
0025 #include "Acts/Utilities/VectorHelpers.hpp"
0026 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0027 #include "ActsExamples/EventData/SpacePoint.hpp"
0028 #include "ActsExamples/EventData/Track.hpp"
0029 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0030
0031 #include <array>
0032 #include <cmath>
0033 #include <cstddef>
0034 #include <optional>
0035 #include <ostream>
0036 #include <span>
0037 #include <stdexcept>
0038 #include <utility>
0039 #include <vector>
0040
0041 namespace ActsExamples {
0042
0043 namespace {
0044
0045 // the step is at most a space point's distance from its own module, so there is
0046 // nothing in between to navigate
0047 using SeedPropagator =
0048     Acts::Propagator<Acts::SympyStepper, Acts::VoidNavigator>;
0049 using SeedPropagatorOptions = SeedPropagator::Options<>;
0050
0051 Acts::Result<Acts::FreeVector> estimateFreeParams(
0052     std::span<const Acts::Vector3> positions, const Acts::Vector3& bField,
0053     double t0, std::span<const double> weights, std::size_t refineIterations) {
0054   // three positions determine a helix exactly, more are fitted
0055   if (positions.size() == 3) {
0056     return Acts::Result<Acts::FreeVector>::success(
0057         Acts::estimateTrackParamsFromSeed(positions[0], t0, positions[1],
0058                                           positions[2], bField));
0059   }
0060
0061   return Acts::estimateTrackParamsFromSpacePoints(positions, bField, t0,
0062                                                   refineIterations, weights);
0063 }
0064
0065 // a space point does not have to lie on the surface of its own measurement
0066 Acts::Result<Acts::BoundVector> transportToSurface(
0067     const SeedPropagator& propagator, const SeedPropagatorOptions& options,
0068     const Acts::FreeVector& freeParams, const Acts::Surface& surface,
0069     const Acts::ParticleHypothesis& hypothesis) {
0070   const Acts::Vector3 direction = freeParams.segment<3>(Acts::eFreeDir0);
0071
0072   // only for the direction to propagate in, so the bounds do not matter
0073   const Acts::Intersection3D intersection =
0074       surface
0075           .intersect(options.geoContext, freeParams.segment<3>(Acts::eFreePos0),
0076                      direction)
0077           .closest();
0078   if (!intersection.isValid()) {
0079     return Acts::Result<Acts::BoundVector>::failure(
0080         Acts::PropagatorError::Failure);
0081   }
0082
0083   SeedPropagatorOptions surfaceOptions = options;
0084   surfaceOptions.direction =
0085       Acts::Direction::fromScalarZeroAsPositive(intersection.pathLength());
0086
0087   const Acts::BoundTrackParameters start =
0088       Acts::BoundTrackParameters::createCurvilinear(
0089           freeParams.segment<4>(Acts::eFreePos0), direction,
0090           freeParams[Acts::eFreeQOverP], std::nullopt, hypothesis);
0091
0092   auto result = propagator.propagate(start, surface, surfaceOptions);
0093   if (!result.ok()) {
0094     return Acts::Result<Acts::BoundVector>::failure(result.error());
0095   }
0096   return Acts::Result<Acts::BoundVector>::success(
0097       result->endParameters.value().parameters());
0098 }
0099
0100 }  // namespace
0101
0102 TrackParamsEstimationAlgorithm::SpacePointWeight
0103 TrackParamsEstimationAlgorithm::inverseRadiusPowerWeight(double exponent) {
0104   return [exponent](const Acts::Vector3& position) {
0105     const double r = Acts::VectorHelpers::perp(position);
0106     if (r <= 0) {
0107       return 1.;
0108     }
0109     return std::pow(r, -exponent);
0110   };
0111 }
0112
0113 TrackParamsEstimationAlgorithm::TrackParamsEstimationAlgorithm(
0114     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0115     : IAlgorithm("TrackParamsEstimationAlgorithm", std::move(logger)),
0116       m_cfg(cfg) {
0117   if (m_cfg.inputSeeds.empty()) {
0118     throw std::invalid_argument("Missing seeds input collection");
0119   }
0120   if (m_cfg.outputTrackParameters.empty()) {
0121     throw std::invalid_argument("Missing track parameters output collection");
0122   }
0123   if (!m_cfg.trackingGeometry) {
0124     throw std::invalid_argument("Missing tracking geometry");
0125   }
0126   if (!m_cfg.magneticField) {
0127     throw std::invalid_argument("Missing magnetic field");
0128   }
0129
0130   m_inputSeeds.initialize(m_cfg.inputSeeds);
0131   m_inputTracks.maybeInitialize(m_cfg.inputProtoTracks);
0132   m_inputParticleHypotheses.maybeInitialize(m_cfg.inputParticleHypotheses);
0133
0134   m_outputTrackParameters.initialize(m_cfg.outputTrackParameters);
0135   m_outputSeeds.maybeInitialize(m_cfg.outputSeeds);
0136   m_outputTracks.maybeInitialize(m_cfg.outputProtoTracks);
0137 }
0138
0139 ProcessCode TrackParamsEstimationAlgorithm::execute(
0140     const AlgorithmContext& ctx) const {
0141   auto const& seeds = m_inputSeeds(ctx);
0142   ACTS_VERBOSE("Read " << seeds.size() << " seeds");
0143
0144   TrackParametersContainer trackParameters;
0145   trackParameters.reserve(seeds.size());
0146
0147   SeedContainer outputSeeds;
0148   if (m_outputSeeds.isInitialized()) {
0149     outputSeeds.assignSpacePointContainer(seeds.spacePointContainer());
0150     outputSeeds.reserve(seeds.size());
0151   }
0152
0153   const ProtoTrackContainer* inputTracks = nullptr;
0154   ProtoTrackContainer outputTracks;
0155   if (m_inputTracks.isInitialized() && m_outputTracks.isInitialized()) {
0156     const auto& inputTracksRef = m_inputTracks(ctx);
0157     if (seeds.size() != inputTracksRef.size()) {
0158       ACTS_FATAL("Inconsistent number of seeds and proto tracks");
0159       return ProcessCode::ABORT;
0160     }
0161     inputTracks = &inputTracksRef;
0162     outputTracks.reserve(seeds.size());
0163   }
0164
0165   const std::vector<Acts::ParticleHypothesis>* inputParticleHypotheses =
0166       nullptr;
0167   if (m_inputParticleHypotheses.isInitialized()) {
0168     const auto& inputParticleHypothesesRef = m_inputParticleHypotheses(ctx);
0169     if (seeds.size() != inputParticleHypothesesRef.size()) {
0170       ACTS_FATAL("Inconsistent number of seeds and particle hypotheses");
0171       return ProcessCode::ABORT;
0172     }
0173     inputParticleHypotheses = &inputParticleHypothesesRef;
0174   }
0175
0176   auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext);
0177
0178   IndexSourceLink::SurfaceAccessor surfaceAccessor{*m_cfg.trackingGeometry};
0179
0180   const SpacePointContainer& spacePoints = seeds.spacePointContainer();
0181
0182   const SeedPropagator propagator(Acts::SympyStepper(m_cfg.magneticField),
0183                                   Acts::VoidNavigator(),
0184                                   logger().cloneWithSuffix("Propagator"));
0185   const SeedPropagatorOptions propagatorOptions(ctx.recoGeoContext,
0186                                                 ctx.magFieldContext);
0187
0188   // reused across seeds by the helix fit
0189   std::array<SpacePointIndex, 3> triplet{};
0190   std::vector<Acts::Vector3> positions;
0191   std::vector<double> weights;
0192
0193   struct {
0194     std::size_t selection = 0;
0195     std::size_t fit = 0;
0196     std::size_t degenerate = 0;
0197     std::size_t transport = 0;
0198
0199     std::size_t total() const {
0200       return selection + fit + degenerate + transport;
0201     }
0202   } skipped;
0203
0204   // Loop over all found seeds to estimate track parameters
0205   for (std::size_t iseed = 0; iseed < seeds.size(); ++iseed) {
0206     const auto& seed = seeds[iseed];
0207     if (seed.spacePoints().size() < 3) {
0208       ACTS_WARNING("Seed " << iseed << " has less than 3 space points, skip");
0209       continue;
0210     }
0211
0212     std::span<const SpacePointIndex> selected = seed.spacePointIndices();
0213
0214     // the triplet selections pick three space points
0215     if (m_cfg.spacePointSelection != SeedSpacePointSelection::All) {
0216       const std::optional<std::array<SpacePointIndex, 3>> selectedTriplet =
0217           selectSeedSpacePoints(spacePoints, seed.spacePointIndices(),
0218                                 m_cfg.spacePointSelection,
0219                                 m_cfg.minTransverseDistance);
0220       if (!selectedTriplet.has_value()) {
0221         ACTS_DEBUG("Seed " << iseed << " failed space point selection, skip");
0222         ++skipped.selection;
0223         continue;
0224       }
0225       triplet = *selectedTriplet;
0226       selected = triplet;
0227     }
0228
0229     // Get the bottom space point and its reference surface
0230     const ConstSpacePointProxy bottomSp = spacePoints.at(selected.front());
0231     if (bottomSp.sourceLinks().empty()) {
0232       ACTS_WARNING("Missing source link in the space point");
0233       continue;
0234     }
0235
0236     const Acts::Vector3 bottomSpVec{bottomSp.x(), bottomSp.y(), bottomSp.z()};
0237
0238     const Acts::SourceLink& bottomSourceLink = bottomSp.sourceLinks()[0];
0239     const Acts::Surface* bottomSurface = surfaceAccessor(bottomSourceLink);
0240     if (bottomSurface == nullptr) {
0241       ACTS_WARNING(
0242           "Surface from source link is not found in the tracking geometry");
0243       continue;
0244     }
0245
0246     // Get the magnetic field at the bottom space point
0247     const auto fieldRes = m_cfg.magneticField->getField(bottomSpVec, bCache);
0248     if (!fieldRes.ok()) {
0249       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0250       return ProcessCode::ABORT;
0251     }
0252     const Acts::Vector3& field = *fieldRes;
0253
0254     if (field.norm() < m_cfg.bFieldMin) {
0255       ACTS_WARNING("Magnetic field at seed " << iseed << " is too small "
0256                                              << field.norm());
0257       continue;
0258     }
0259
0260     positions.clear();
0261     for (const SpacePointIndex index : selected) {
0262       const ConstSpacePointProxy sp = spacePoints.at(index);
0263       positions.emplace_back(sp.x(), sp.y(), sp.z());
0264     }
0265
0266     if (m_cfg.spacePointWeight) {
0267       weights.clear();
0268       for (const Acts::Vector3& position : positions) {
0269         weights.emplace_back(m_cfg.spacePointWeight(position));
0270       }
0271     }
0272
0273     const double t0 = std::isnan(bottomSp.time()) ? 0.0 : bottomSp.time();
0274
0275     const Acts::Result<Acts::FreeVector> freeParams = estimateFreeParams(
0276         positions, field, t0, weights, m_cfg.geometricRefineIterations);
0277     if (!freeParams.ok()) {
0278       ACTS_DEBUG("Seed " << iseed << " could not be fitted: "
0279                          << freeParams.error().message());
0280       ++skipped.fit;
0281       continue;
0282     }
0283     // degenerate space points make the estimate not a number, a straight line
0284     // fit leaves q/p at zero, which the covariance transport cannot take
0285     if (!freeParams->allFinite() || (*freeParams)[Acts::eFreeQOverP] == 0) {
0286       ACTS_DEBUG("Seed " << iseed << " has a degenerate estimate, skip");
0287       ++skipped.degenerate;
0288       continue;
0289     }
0290
0291     const Acts::ParticleHypothesis hypothesis =
0292         inputParticleHypotheses != nullptr ? inputParticleHypotheses->at(iseed)
0293                                            : m_cfg.particleHypothesis;
0294
0295     const Acts::Result<Acts::BoundVector> boundParams = transportToSurface(
0296         propagator, propagatorOptions, *freeParams, *bottomSurface, hypothesis);
0297     if (!boundParams.ok()) {
0298       ACTS_DEBUG("Seed " << iseed
0299                          << " could not be transported to the surface of its "
0300                             "first space point: "
0301                          << boundParams.error().message());
0302       ++skipped.transport;
0303       continue;
0304     }
0305
0306     Acts::EstimateTrackParamCovarianceConfig config{
0307         .initialSigmas =
0308             Eigen::Map<const Acts::BoundVector>{m_cfg.initialSigmas.data()},
0309         .initialSigmaQoverPt = m_cfg.initialSigmaQoverPt,
0310         .initialSigmaPtRel = m_cfg.initialSigmaPtRel,
0311         .initialVarInflation = Eigen::Map<const Acts::BoundVector>{
0312             m_cfg.initialVarInflation.data()}};
0313
0314     const Acts::BoundMatrix cov = Acts::estimateTrackParamCovariance(
0315         config, *boundParams, !std::isnan(bottomSp.time()));
0316
0317     const TrackParameters& trackParams = trackParameters.emplace_back(
0318         bottomSurface->getSharedPtr(), *boundParams, cov, hypothesis);
0319     ACTS_VERBOSE("Estimated track parameters: " << trackParams);
0320     if (m_outputSeeds.isInitialized()) {
0321       auto newSp = outputSeeds.createSeed();
0322       // TODO copy shorthand
0323       newSp.assignSpacePointIndices(seed.spacePointIndices());
0324       newSp.quality() = seed.quality();
0325       newSp.vertexZ() = seed.vertexZ();
0326     }
0327     if (m_outputTracks.isInitialized() && inputTracks != nullptr) {
0328       outputTracks.push_back(inputTracks->at(iseed));
0329     }
0330   }
0331
0332   ACTS_DEBUG("Estimated " << trackParameters.size() << " track parameters from "
0333                           << seeds.size() << " seeds");
0334   if (skipped.total() > 0) {
0335     ACTS_DEBUG(
0336         "Skipped " << skipped.selection
0337                    << " seeds without a space point selection, " << skipped.fit
0338                    << " without a fit, " << skipped.degenerate
0339                    << " with a degenerate estimate and " << skipped.transport
0340                    << " without a transport to the surface");
0341   }
0342
0343   m_outputTrackParameters(ctx, std::move(trackParameters));
0344   if (m_outputSeeds.isInitialized()) {
0345     m_outputSeeds(ctx, std::move(outputSeeds));
0346   }
0347
0348   if (m_outputTracks.isInitialized()) {
0349     m_outputTracks(ctx, std::move(outputTracks));
0350   }
0351
0352   return ProcessCode::SUCCESS;
0353 }
0354
0355 }  // namespace ActsExamples