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