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   m_inputFreeParameters.maybeInitialize(m_cfg.inputFreeParameters);
0134
0135   m_outputTrackParameters.initialize(m_cfg.outputTrackParameters);
0136   m_outputSeeds.maybeInitialize(m_cfg.outputSeeds);
0137   m_outputTracks.maybeInitialize(m_cfg.outputProtoTracks);
0138 }
0139
0140 ProcessCode TrackParamsEstimationAlgorithm::execute(
0141     const AlgorithmContext& ctx) const {
0142   auto const& seeds = m_inputSeeds(ctx);
0143   ACTS_VERBOSE("Read " << seeds.size() << " seeds");
0144
0145   TrackParametersContainer trackParameters;
0146   trackParameters.reserve(seeds.size());
0147
0148   SeedContainer outputSeeds;
0149   if (m_outputSeeds.isInitialized()) {
0150     outputSeeds.assignSpacePointContainer(seeds.spacePointContainer());
0151     outputSeeds.reserve(seeds.size());
0152   }
0153
0154   const ProtoTrackContainer* inputTracks = nullptr;
0155   ProtoTrackContainer outputTracks;
0156   if (m_inputTracks.isInitialized() && m_outputTracks.isInitialized()) {
0157     const auto& inputTracksRef = m_inputTracks(ctx);
0158     if (seeds.size() != inputTracksRef.size()) {
0159       ACTS_FATAL("Inconsistent number of seeds and proto tracks");
0160       return ProcessCode::ABORT;
0161     }
0162     inputTracks = &inputTracksRef;
0163     outputTracks.reserve(seeds.size());
0164   }
0165
0166   const std::vector<Acts::ParticleHypothesis>* inputParticleHypotheses =
0167       nullptr;
0168   if (m_inputParticleHypotheses.isInitialized()) {
0169     const auto& inputParticleHypothesesRef = m_inputParticleHypotheses(ctx);
0170     if (seeds.size() != inputParticleHypothesesRef.size()) {
0171       ACTS_FATAL("Inconsistent number of seeds and particle hypotheses");
0172       return ProcessCode::ABORT;
0173     }
0174     inputParticleHypotheses = &inputParticleHypothesesRef;
0175   }
0176
0177   const std::vector<Acts::FreeVector>* inputFreeParameters = nullptr;
0178   if (m_inputFreeParameters.isInitialized()) {
0179     const auto& inputFreeParametersRef = m_inputFreeParameters(ctx);
0180     if (seeds.size() != inputFreeParametersRef.size()) {
0181       ACTS_FATAL("Inconsistent number of seeds and free parameters");
0182       return ProcessCode::ABORT;
0183     }
0184     inputFreeParameters = &inputFreeParametersRef;
0185   }
0186
0187   auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext);
0188
0189   IndexSourceLink::SurfaceAccessor surfaceAccessor{*m_cfg.trackingGeometry};
0190
0191   const SpacePointContainer& spacePoints = seeds.spacePointContainer();
0192
0193   const SeedPropagator propagator(Acts::SympyStepper(m_cfg.magneticField),
0194                                   Acts::VoidNavigator(),
0195                                   logger().cloneWithSuffix("Propagator"));
0196   const SeedPropagatorOptions propagatorOptions(ctx.recoGeoContext,
0197                                                 ctx.magFieldContext);
0198
0199   // reused across seeds by the helix fit
0200   std::array<SpacePointIndex, 3> triplet{};
0201   std::vector<Acts::Vector3> positions;
0202   std::vector<double> weights;
0203
0204   struct {
0205     std::size_t selection = 0;
0206     std::size_t fit = 0;
0207     std::size_t degenerate = 0;
0208     std::size_t transport = 0;
0209
0210     std::size_t total() const {
0211       return selection + fit + degenerate + transport;
0212     }
0213   } skipped;
0214
0215   // Loop over all found seeds to estimate track parameters
0216   for (std::size_t iseed = 0; iseed < seeds.size(); ++iseed) {
0217     const auto& seed = seeds[iseed];
0218     if (seed.spacePoints().size() < 3) {
0219       ACTS_WARNING("Seed " << iseed << " has less than 3 space points, skip");
0220       continue;
0221     }
0222
0223     std::span<const SpacePointIndex> selected = seed.spacePointIndices();
0224
0225     // given free parameters replace the fit, unless they have no momentum
0226     const bool hasFreeParams =
0227         inputFreeParameters != nullptr &&
0228         inputFreeParameters->at(iseed)[Acts::eFreeQOverP] != 0;
0229
0230     // the triplet selections pick three space points
0231     if (m_cfg.spacePointSelection != SeedSpacePointSelection::All) {
0232       const std::optional<std::array<SpacePointIndex, 3>> selectedTriplet =
0233           selectSeedSpacePoints(spacePoints, seed.spacePointIndices(),
0234                                 m_cfg.spacePointSelection,
0235                                 m_cfg.minTransverseDistance);
0236       // a seed with given free parameters keeps all of its space points
0237       if (!selectedTriplet.has_value() && !hasFreeParams) {
0238         ACTS_DEBUG("Seed " << iseed << " failed space point selection, skip");
0239         ++skipped.selection;
0240         continue;
0241       }
0242       if (selectedTriplet.has_value()) {
0243         triplet = *selectedTriplet;
0244         selected = triplet;
0245       }
0246     }
0247
0248     // Get the bottom space point and its reference surface
0249     const ConstSpacePointProxy bottomSp = spacePoints.at(selected.front());
0250     if (bottomSp.sourceLinks().empty()) {
0251       ACTS_WARNING("Missing source link in the space point");
0252       continue;
0253     }
0254
0255     const Acts::Vector3 bottomSpVec{bottomSp.x(), bottomSp.y(), bottomSp.z()};
0256
0257     const Acts::SourceLink& bottomSourceLink = bottomSp.sourceLinks()[0];
0258     const Acts::Surface* bottomSurface = surfaceAccessor(bottomSourceLink);
0259     if (bottomSurface == nullptr) {
0260       ACTS_WARNING(
0261           "Surface from source link is not found in the tracking geometry");
0262       continue;
0263     }
0264
0265     // Get the magnetic field at the bottom space point
0266     const auto fieldRes = m_cfg.magneticField->getField(bottomSpVec, bCache);
0267     if (!fieldRes.ok()) {
0268       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0269       return ProcessCode::ABORT;
0270     }
0271     const Acts::Vector3& field = *fieldRes;
0272
0273     if (field.norm() < m_cfg.bFieldMin) {
0274       ACTS_WARNING("Magnetic field at seed " << iseed << " is too small "
0275                                              << field.norm());
0276       continue;
0277     }
0278
0279     positions.clear();
0280     for (const SpacePointIndex index : selected) {
0281       const ConstSpacePointProxy sp = spacePoints.at(index);
0282       positions.emplace_back(sp.x(), sp.y(), sp.z());
0283     }
0284
0285     if (m_cfg.spacePointWeight) {
0286       weights.clear();
0287       for (const Acts::Vector3& position : positions) {
0288         weights.emplace_back(m_cfg.spacePointWeight(position));
0289       }
0290     }
0291
0292     const double t0 = std::isnan(bottomSp.time()) ? 0.0 : bottomSp.time();
0293
0294     Acts::Result<Acts::FreeVector> freeParams =
0295         hasFreeParams ? Acts::Result<Acts::FreeVector>::success(
0296                             inputFreeParameters->at(iseed))
0297                       : estimateFreeParams(positions, field, t0, weights,
0298                                            m_cfg.geometricRefineIterations);
0299     if (hasFreeParams) {
0300       (*freeParams)[Acts::eFreeTime] = t0;
0301     }
0302     if (!freeParams.ok()) {
0303       ACTS_DEBUG("Seed " << iseed << " could not be fitted: "
0304                          << freeParams.error().message());
0305       ++skipped.fit;
0306       continue;
0307     }
0308     // degenerate space points make the estimate not a number, a straight line
0309     // fit leaves q/p at zero, which the covariance transport cannot take
0310     if (!freeParams->allFinite() || (*freeParams)[Acts::eFreeQOverP] == 0) {
0311       ACTS_DEBUG("Seed " << iseed << " has a degenerate estimate, skip");
0312       ++skipped.degenerate;
0313       continue;
0314     }
0315
0316     const Acts::ParticleHypothesis hypothesis =
0317         inputParticleHypotheses != nullptr ? inputParticleHypotheses->at(iseed)
0318                                            : m_cfg.particleHypothesis;
0319
0320     Acts::Result<Acts::BoundVector> boundParams = transportToSurface(
0321         propagator, propagatorOptions, *freeParams, *bottomSurface, hypothesis);
0322     // given free parameters which do not reach the surface give way to the fit
0323     if (!boundParams.ok() && hasFreeParams) {
0324       freeParams = estimateFreeParams(positions, field, t0, weights,
0325                                       m_cfg.geometricRefineIterations);
0326       if (freeParams.ok() && freeParams->allFinite() &&
0327           (*freeParams)[Acts::eFreeQOverP] != 0) {
0328         boundParams =
0329             transportToSurface(propagator, propagatorOptions, *freeParams,
0330                                *bottomSurface, hypothesis);
0331       }
0332     }
0333     if (!boundParams.ok()) {
0334       ACTS_DEBUG("Seed " << iseed
0335                          << " could not be transported to the surface of its "
0336                             "first space point: "
0337                          << boundParams.error().message());
0338       ++skipped.transport;
0339       continue;
0340     }
0341
0342     Acts::EstimateTrackParamCovarianceConfig config{
0343         .initialSigmas =
0344             Eigen::Map<const Acts::BoundVector>{m_cfg.initialSigmas.data()},
0345         .initialSigmaQoverPt = m_cfg.initialSigmaQoverPt,
0346         .initialSigmaPtRel = m_cfg.initialSigmaPtRel,
0347         .initialVarInflation = Eigen::Map<const Acts::BoundVector>{
0348             m_cfg.initialVarInflation.data()}};
0349
0350     const Acts::BoundMatrix cov = Acts::estimateTrackParamCovariance(
0351         config, *boundParams, !std::isnan(bottomSp.time()));
0352
0353     const TrackParameters& trackParams = trackParameters.emplace_back(
0354         bottomSurface->getSharedPtr(), *boundParams, cov, hypothesis);
0355     ACTS_VERBOSE("Estimated track parameters: " << trackParams);
0356     if (m_outputSeeds.isInitialized()) {
0357       auto newSp = outputSeeds.createSeed();
0358       // TODO copy shorthand
0359       newSp.assignSpacePointIndices(seed.spacePointIndices());
0360       newSp.quality() = seed.quality();
0361       newSp.vertexZ() = seed.vertexZ();
0362     }
0363     if (m_outputTracks.isInitialized() && inputTracks != nullptr) {
0364       outputTracks.push_back(inputTracks->at(iseed));
0365     }
0366   }
0367
0368   ACTS_DEBUG("Estimated " << trackParameters.size() << " track parameters from "
0369                           << seeds.size() << " seeds");
0370   if (skipped.total() > 0) {
0371     ACTS_DEBUG(
0372         "Skipped " << skipped.selection
0373                    << " seeds without a space point selection, " << skipped.fit
0374                    << " without a fit, " << skipped.degenerate
0375                    << " with a degenerate estimate and " << skipped.transport
0376                    << " without a transport to the surface");
0377   }
0378
0379   m_outputTrackParameters(ctx, std::move(trackParameters));
0380   if (m_outputSeeds.isInitialized()) {
0381     m_outputSeeds(ctx, std::move(outputSeeds));
0382   }
0383
0384   if (m_outputTracks.isInitialized()) {
0385     m_outputTracks(ctx, std::move(outputTracks));
0386   }
0387
0388   return ProcessCode::SUCCESS;
0389 }
0390
0391 }  // namespace ActsExamples