Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.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/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0014 #include "Acts/Propagator/SympyStepper.hpp"
0015 #include "Acts/Utilities/AnnealingUtility.hpp"
0016 #include "Acts/Utilities/Logger.hpp"
0017 #include "Acts/Utilities/Result.hpp"
0018 #include "Acts/Vertexing/AdaptiveGridTrackDensity.hpp"
0019 #include "Acts/Vertexing/AdaptiveMultiVertexFinder.hpp"
0020 #include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp"
0021 #include "Acts/Vertexing/GaussianTrackDensity.hpp"
0022 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0023 #include "Acts/Vertexing/TrackAtVertex.hpp"
0024 #include "Acts/Vertexing/TrackDensityVertexFinder.hpp"
0025 #include "Acts/Vertexing/Vertex.hpp"
0026 #include "Acts/Vertexing/VertexingOptions.hpp"
0027 #include "ActsExamples/EventData/ProtoVertex.hpp"
0028 #include "ActsExamples/EventData/SimParticle.hpp"
0029 #include "ActsExamples/EventData/SimVertex.hpp"
0030 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0031 #include "ActsExamples/Framework/ProcessCode.hpp"
0032 #include "ActsExamples/TruthTracking/TruthVertexFinder.hpp"
0033
0034 #include <algorithm>
0035 #include <memory>
0036 #include <optional>
0037 #include <ostream>
0038 #include <stdexcept>
0039 #include <system_error>
0040 #include <unordered_map>
0041 #include <utility>
0042
0043 #include "TruthVertexSeeder.hpp"
0044 #include "VertexingHelpers.hpp"
0045
0046 namespace ActsExamples {
0047
0048 AdaptiveMultiVertexFinderAlgorithm::AdaptiveMultiVertexFinderAlgorithm(
0049     const Config& config, Acts::Logging::Level level)
0050     : IAlgorithm("AdaptiveMultiVertexFinder", level),
0051       m_cfg(config),
0052       m_propagator{[&]() {
0053         // Set up SympyStepper
0054         Acts::SympyStepper stepper(m_cfg.bField);
0055
0056         // Set up the propagator
0057         return std::make_shared<Propagator>(stepper);
0058       }()},
0059       m_ipEstimator{[&]() {
0060         // Set up ImpactPointEstimator
0061         Acts::ImpactPointEstimator::Config ipEstimatorCfg(m_cfg.bField,
0062                                                           m_propagator);
0063         return Acts::ImpactPointEstimator(
0064             ipEstimatorCfg, logger().cloneWithSuffix("ImpactPointEstimator"));
0065       }()},
0066       m_linearizer{[&] {
0067         // Set up the helical track linearizer
0068         Linearizer::Config ltConfig;
0069         ltConfig.bField = m_cfg.bField;
0070         ltConfig.propagator = m_propagator;
0071         return Linearizer(ltConfig,
0072                           logger().cloneWithSuffix("HelicalTrackLinearizer"));
0073       }()},
0074       m_vertexSeeder{makeVertexSeeder()},
0075       m_vertexFinder{makeVertexFinder(m_vertexSeeder)} {
0076   if (m_cfg.inputTrackParameters.empty()) {
0077     throw std::invalid_argument("Missing input track parameter collection");
0078   }
0079   if (m_cfg.outputProtoVertices.empty()) {
0080     throw std::invalid_argument("Missing output proto vertices collection");
0081   }
0082   if (m_cfg.outputVertices.empty()) {
0083     throw std::invalid_argument("Missing output vertices collection");
0084   }
0085   if (m_cfg.seedFinder == SeedFinder::TruthSeeder &&
0086       m_cfg.inputTruthVertices.empty()) {
0087     throw std::invalid_argument("Missing input truth vertex collection");
0088   }
0089
0090   // Sanitize the configuration
0091   if (m_cfg.seedFinder != SeedFinder::TruthSeeder &&
0092       (!m_cfg.inputTruthParticles.empty() ||
0093        !m_cfg.inputTruthVertices.empty())) {
0094     ACTS_INFO("Ignoring truth input as seed finder is not TruthSeeder");
0095     m_cfg.inputTruthVertices.clear();
0096     m_cfg.inputTruthVertices.clear();
0097   }
0098
0099   m_inputTrackParameters.initialize(m_cfg.inputTrackParameters);
0100   m_inputTruthParticles.maybeInitialize(m_cfg.inputTruthParticles);
0101   m_inputTruthVertices.maybeInitialize(m_cfg.inputTruthVertices);
0102   m_outputProtoVertices.initialize(m_cfg.outputProtoVertices);
0103   m_outputVertices.initialize(m_cfg.outputVertices);
0104 }
0105
0106 std::unique_ptr<Acts::IVertexFinder>
0107 AdaptiveMultiVertexFinderAlgorithm::makeVertexSeeder() const {
0108   if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
0109     using Seeder = ActsExamples::TruthVertexSeeder;
0110     Seeder::Config seederConfig;
0111     seederConfig.useXY = false;
0112     seederConfig.useTime = m_cfg.useTime;
0113     return std::make_unique<Seeder>(seederConfig);
0114   }
0115
0116   if (m_cfg.seedFinder == SeedFinder::GaussianSeeder) {
0117     using Seeder = Acts::TrackDensityVertexFinder;
0118     Acts::GaussianTrackDensity::Config trkDensityCfg;
0119     trkDensityCfg.extractParameters
0120         .connect<&Acts::InputTrack::extractParameters>();
0121     return std::make_unique<Seeder>(Seeder::Config{trkDensityCfg});
0122   }
0123
0124   if (m_cfg.seedFinder == SeedFinder::AdaptiveGridSeeder) {
0125     // Set up track density used during vertex seeding
0126     Acts::AdaptiveGridTrackDensity::Config trkDensityCfg;
0127     // Bin extent in z-direction
0128     trkDensityCfg.spatialBinExtent = m_cfg.spatialBinExtent;
0129     // Bin extent in t-direction
0130     trkDensityCfg.temporalBinExtent = m_cfg.temporalBinExtent;
0131     trkDensityCfg.useTime = m_cfg.useTime;
0132     Acts::AdaptiveGridTrackDensity trkDensity(trkDensityCfg);
0133
0134     // Set up vertex seeder and finder
0135     using Seeder = Acts::AdaptiveGridDensityVertexFinder;
0136     Seeder::Config seederConfig(trkDensity);
0137     seederConfig.extractParameters
0138         .connect<&Acts::InputTrack::extractParameters>();
0139     return std::make_unique<Seeder>(seederConfig);
0140   }
0141
0142   throw std::invalid_argument("Unknown seed finder");
0143 }
0144
0145 Acts::AdaptiveMultiVertexFinder
0146 AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder(
0147     std::shared_ptr<const Acts::IVertexFinder> seedFinder) const {
0148   // Set up deterministic annealing with user-defined temperatures
0149   Acts::AnnealingUtility annealingUtility(m_cfg.annealingConfig);
0150
0151   // Set up the vertex fitter with user-defined annealing
0152   Fitter::Config fitterCfg(m_ipEstimator);
0153   fitterCfg.annealingTool = annealingUtility;
0154   fitterCfg.minWeight = m_cfg.minWeight;
0155   fitterCfg.doSmoothing = m_cfg.doSmoothing;
0156   fitterCfg.useTime = m_cfg.useTime;
0157   fitterCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0158   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&m_linearizer);
0159   Fitter fitter(std::move(fitterCfg),
0160                 logger().cloneWithSuffix("AdaptiveMultiVertexFitter"));
0161
0162   Acts::AdaptiveMultiVertexFinder::Config finderConfig(
0163       std::move(fitter), std::move(seedFinder), m_ipEstimator, m_cfg.bField);
0164   // Set the initial variance of the 4D vertex position. Since time is on a
0165   // numerical scale, we have to provide a greater value in the corresponding
0166   // dimension.
0167   finderConfig.initialVariances = m_cfg.initialVariances;
0168   finderConfig.tracksMaxZinterval = m_cfg.tracksMaxZinterval;
0169   finderConfig.maxIterations = m_cfg.maxIterations;
0170   finderConfig.useTime = m_cfg.useTime;
0171   // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
0172   finderConfig.tracksMaxSignificance = 5;
0173   // This should be used consistently with and without time
0174   finderConfig.doFullSplitting = m_cfg.doFullSplitting;
0175   // 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
0176   finderConfig.maxMergeVertexSignificance = 3;
0177   if (m_cfg.useTime) {
0178     // When using time, we have an extra contribution to the chi2 by the time
0179     // coordinate. We thus need to increase tracksMaxSignificance (i.e., the
0180     // maximum chi2 that a track can have to be associated with a vertex).
0181     // Using the same p-value for 3 dof instead of 2.
0182     // 6.7 corresponds to a p-value of ~0.92 using `chi2(x=6.7,ndf=3)`
0183     finderConfig.tracksMaxSignificance = 6.7;
0184     // Using the same p-value for 2 dof instead of 1.
0185     // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
0186     finderConfig.maxMergeVertexSignificance = 5;
0187   }
0188
0189   finderConfig.extractParameters
0190       .template connect<&Acts::InputTrack::extractParameters>();
0191
0192   if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
0193     finderConfig.doNotBreakWhileSeeding = true;
0194   }
0195
0196   finderConfig.tracksMaxSignificance =
0197       m_cfg.tracksMaxSignificance.value_or(finderConfig.tracksMaxSignificance);
0198   finderConfig.maxMergeVertexSignificance =
0199       m_cfg.maxMergeVertexSignificance.value_or(
0200           finderConfig.maxMergeVertexSignificance);
0201
0202   // Instantiate the finder
0203   return Acts::AdaptiveMultiVertexFinder(std::move(finderConfig),
0204                                          logger().clone());
0205 }
0206
0207 ProcessCode AdaptiveMultiVertexFinderAlgorithm::execute(
0208     const AlgorithmContext& ctx) const {
0209   const auto& inputTrackParameters = m_inputTrackParameters(ctx);
0210
0211   auto inputTracks = makeInputTracks(inputTrackParameters);
0212
0213   if (inputTrackParameters.size() != inputTracks.size()) {
0214     ACTS_ERROR("Input track containers do not align: "
0215                << inputTrackParameters.size() << " != " << inputTracks.size());
0216   }
0217
0218   for (const auto& trk : inputTrackParameters) {
0219     if (trk.covariance() && trk.covariance()->determinant() <= 0) {
0220       // actually we should consider this as an error but I do not want the CI
0221       // to fail
0222       ACTS_WARNING("input track " << trk << " has det(cov) = "
0223                                   << trk.covariance()->determinant());
0224     }
0225   }
0226
0227   // The vertex finder state
0228   auto state = m_vertexFinder.makeState(ctx.magFieldContext);
0229
0230   // In case of the truth seeder, we need to wire the truth vertices into the
0231   // vertex finder
0232   if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
0233     const auto& truthParticles = m_inputTruthParticles(ctx);
0234     const auto& truthVertices = m_inputTruthVertices(ctx);
0235
0236     auto& vertexSeederState =
0237         state.as<Acts::AdaptiveMultiVertexFinder::State>()
0238             .seedFinderState.as<TruthVertexSeeder::State>();
0239
0240     std::map<SimVertexBarcode, std::size_t> vertexParticleCount;
0241
0242     for (const auto& truthVertex : truthVertices) {
0243       // Skip secondary vertices
0244       if (truthVertex.vertexId().vertexSecondary() != 0) {
0245         continue;
0246       }
0247       vertexSeederState.truthVertices.push_back(truthVertex);
0248
0249       // Count the number of particles associated with each vertex
0250       std::size_t particleCount = 0;
0251       for (const auto& particle : truthParticles) {
0252         if (static_cast<SimVertexBarcode>(particle.particleId().vertexId()) ==
0253             truthVertex.vertexId()) {
0254           ++particleCount;
0255         }
0256       }
0257       vertexParticleCount[truthVertex.vertexId()] = particleCount;
0258     }
0259
0260     // sort by number of particles
0261     std::ranges::sort(vertexSeederState.truthVertices, {},
0262                       [&vertexParticleCount](const auto& v) {
0263                         return vertexParticleCount[v.vertexId()];
0264                       });
0265
0266     ACTS_INFO("Got " << truthVertices.size() << " truth vertices and selected "
0267                      << vertexSeederState.truthVertices.size() << " in event");
0268   }
0269
0270   // Default vertexing options, this is where e.g. a constraint could be set
0271   Options finderOpts(ctx.geoContext, ctx.magFieldContext);
0272
0273   VertexCollection vertices;
0274
0275   if (inputTrackParameters.empty()) {
0276     ACTS_DEBUG("Empty track parameter collection found, skipping vertexing");
0277   } else {
0278     ACTS_DEBUG("Have " << inputTrackParameters.size()
0279                        << " input track parameters, running vertexing");
0280     // find vertices
0281     auto result = m_vertexFinder.find(inputTracks, finderOpts, state);
0282
0283     if (result.ok()) {
0284       vertices = std::move(result.value());
0285     } else {
0286       ACTS_ERROR("Error in vertex finder: " << result.error().message());
0287     }
0288   }
0289
0290   // show some debug output
0291   ACTS_INFO("Found " << vertices.size() << " vertices in event");
0292   for (const auto& vtx : vertices) {
0293     ACTS_DEBUG("Found vertex at " << vtx.fullPosition().transpose() << " with "
0294                                   << vtx.tracks().size() << " tracks.");
0295   }
0296
0297   // store proto vertices extracted from the found vertices
0298   m_outputProtoVertices(ctx, makeProtoVertices(inputTracks, vertices));
0299
0300   // store found vertices
0301   m_outputVertices(ctx, std::move(vertices));
0302
0303   return ProcessCode::SUCCESS;
0304 }
0305
0306 }  // namespace ActsExamples