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     return std::make_unique<Seeder>(seederConfig);
0110   }
0111
0112   if (m_cfg.seedFinder == SeedFinder::GaussianSeeder) {
0113     using Seeder = Acts::TrackDensityVertexFinder;
0114     Acts::GaussianTrackDensity::Config trkDensityCfg;
0115     trkDensityCfg.extractParameters
0116         .connect<&Acts::InputTrack::extractParameters>();
0117     return std::make_unique<Seeder>(
0118         Seeder::Config{Acts::GaussianTrackDensity(trkDensityCfg)});
0119   }
0120
0121   if (m_cfg.seedFinder == SeedFinder::AdaptiveGridSeeder) {
0122     // Set up track density used during vertex seeding
0123     Acts::AdaptiveGridTrackDensity::Config trkDensityCfg;
0124     // Bin extent in z-direction
0125     trkDensityCfg.spatialBinExtent = m_cfg.spatialBinExtent;
0126     // Bin extent in t-direction
0127     trkDensityCfg.temporalBinExtent = m_cfg.temporalBinExtent;
0128     trkDensityCfg.useTime = m_cfg.useTime;
0129     Acts::AdaptiveGridTrackDensity trkDensity(trkDensityCfg);
0130
0131     // Set up vertex seeder and finder
0132     using Seeder = Acts::AdaptiveGridDensityVertexFinder;
0133     Seeder::Config seederConfig(trkDensity);
0134     seederConfig.extractParameters
0135         .connect<&Acts::InputTrack::extractParameters>();
0136     return std::make_unique<Seeder>(seederConfig);
0137   }
0138
0139   throw std::invalid_argument("Unknown seed finder");
0140 }
0141
0142 Acts::AdaptiveMultiVertexFinder
0143 AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder(
0144     std::shared_ptr<const Acts::IVertexFinder> seedFinder) const {
0145   // Set up deterministic annealing with user-defined temperatures
0146   Acts::AnnealingUtility annealingUtility(m_cfg.annealingConfig);
0147
0148   // Set up the vertex fitter with user-defined annealing
0149   Fitter::Config fitterCfg(m_ipEstimator);
0150   fitterCfg.annealingTool = annealingUtility;
0151   fitterCfg.minWeight = m_cfg.minWeight;
0152   fitterCfg.doSmoothing = m_cfg.doSmoothing;
0153   fitterCfg.useTime = m_cfg.useTime;
0154   fitterCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0155   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&m_linearizer);
0156   Fitter fitter(std::move(fitterCfg),
0157                 logger().cloneWithSuffix("AdaptiveMultiVertexFitter"));
0158
0159   Acts::AdaptiveMultiVertexFinder::Config finderConfig(
0160       std::move(fitter), std::move(seedFinder), m_ipEstimator, m_cfg.bField);
0161   // Set the initial variance of the 4D vertex position. Since time is on a
0162   // numerical scale, we have to provide a greater value in the corresponding
0163   // dimension.
0164   finderConfig.initialVariances = m_cfg.initialVariances;
0165   finderConfig.tracksMaxZinterval = m_cfg.tracksMaxZinterval;
0166   finderConfig.maxIterations = m_cfg.maxIterations;
0167   finderConfig.useTime = m_cfg.useTime;
0168   // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
0169   finderConfig.tracksMaxSignificance = 5;
0170   // This should be used consistently with and without time
0171   finderConfig.doFullSplitting = m_cfg.doFullSplitting;
0172   // 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
0173   finderConfig.maxMergeVertexSignificance = 3;
0174   if (m_cfg.useTime) {
0175     // When using time, we have an extra contribution to the chi2 by the time
0176     // coordinate. We thus need to increase tracksMaxSignificance (i.e., the
0177     // maximum chi2 that a track can have to be associated with a vertex).
0178     // Using the same p-value for 3 dof instead of 2.
0179     // 6.7 corresponds to a p-value of ~0.92 using `chi2(x=6.7,ndf=3)`
0180     finderConfig.tracksMaxSignificance = 6.7;
0181     // Using the same p-value for 2 dof instead of 1.
0182     // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
0183     finderConfig.maxMergeVertexSignificance = 5;
0184   }
0185
0186   finderConfig.extractParameters
0187       .template connect<&Acts::InputTrack::extractParameters>();
0188
0189   if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
0190     finderConfig.doNotBreakWhileSeeding = true;
0191   }
0192
0193   finderConfig.tracksMaxSignificance =
0194       m_cfg.tracksMaxSignificance.value_or(finderConfig.tracksMaxSignificance);
0195   finderConfig.maxMergeVertexSignificance =
0196       m_cfg.maxMergeVertexSignificance.value_or(
0197           finderConfig.maxMergeVertexSignificance);
0198
0199   // Instantiate the finder
0200   return Acts::AdaptiveMultiVertexFinder(std::move(finderConfig),
0201                                          logger().clone());
0202 }
0203
0204 ProcessCode AdaptiveMultiVertexFinderAlgorithm::execute(
0205     const AlgorithmContext& ctx) const {
0206   const auto& inputTrackParameters = m_inputTrackParameters(ctx);
0207
0208   auto inputTracks = makeInputTracks(inputTrackParameters);
0209
0210   if (inputTrackParameters.size() != inputTracks.size()) {
0211     ACTS_ERROR("Input track containers do not align: "
0212                << inputTrackParameters.size() << " != " << inputTracks.size());
0213   }
0214
0215   for (const auto& trk : inputTrackParameters) {
0216     if (trk.covariance() && trk.covariance()->determinant() <= 0) {
0217       // actually we should consider this as an error but I do not want the CI
0218       // to fail
0219       ACTS_WARNING("input track " << trk << " has det(cov) = "
0220                                   << trk.covariance()->determinant());
0221     }
0222   }
0223
0224   // The vertex finder state
0225   auto state = m_vertexFinder.makeState(ctx.magFieldContext);
0226
0227   // In case of the truth seeder, we need to wire the truth vertices into the
0228   // vertex finder
0229   if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
0230     const auto& truthParticles = m_inputTruthParticles(ctx);
0231     const auto& truthVertices = m_inputTruthVertices(ctx);
0232
0233     auto& vertexSeederState =
0234         state.as<Acts::AdaptiveMultiVertexFinder::State>()
0235             .seedFinderState.as<TruthVertexSeeder::State>();
0236
0237     std::map<SimVertexBarcode, std::size_t> vertexParticleCount;
0238
0239     for (const auto& truthVertex : truthVertices) {
0240       // Skip secondary vertices
0241       if (truthVertex.vertexId().vertexSecondary() != 0) {
0242         continue;
0243       }
0244       vertexSeederState.truthVertices.push_back(truthVertex);
0245
0246       // Count the number of particles associated with each vertex
0247       std::size_t particleCount = 0;
0248       for (const auto& particle : truthParticles) {
0249         if (static_cast<SimVertexBarcode>(particle.particleId().vertexId()) ==
0250             truthVertex.vertexId()) {
0251           ++particleCount;
0252         }
0253       }
0254       vertexParticleCount[truthVertex.vertexId()] = particleCount;
0255     }
0256
0257     // sort by number of particles
0258     std::ranges::sort(vertexSeederState.truthVertices, {},
0259                       [&vertexParticleCount](const auto& v) {
0260                         return vertexParticleCount[v.vertexId()];
0261                       });
0262
0263     ACTS_INFO("Got " << truthVertices.size() << " truth vertices and selected "
0264                      << vertexSeederState.truthVertices.size() << " in event");
0265   }
0266
0267   // Default vertexing options, this is where e.g. a constraint could be set
0268   Options finderOpts(ctx.geoContext, ctx.magFieldContext);
0269
0270   VertexCollection vertices;
0271
0272   if (inputTrackParameters.empty()) {
0273     ACTS_DEBUG("Empty track parameter collection found, skipping vertexing");
0274   } else {
0275     ACTS_DEBUG("Have " << inputTrackParameters.size()
0276                        << " input track parameters, running vertexing");
0277     // find vertices
0278     auto result = m_vertexFinder.find(inputTracks, finderOpts, state);
0279
0280     if (result.ok()) {
0281       vertices = std::move(result.value());
0282     } else {
0283       ACTS_ERROR("Error in vertex finder: " << result.error().message());
0284     }
0285   }
0286
0287   // show some debug output
0288   ACTS_DEBUG("Found " << vertices.size() << " vertices in event");
0289   for (const auto& vtx : vertices) {
0290     ACTS_DEBUG("Found vertex at " << vtx.fullPosition().transpose() << " with "
0291                                   << vtx.tracks().size() << " tracks.");
0292   }
0293
0294   // store proto vertices extracted from the found vertices
0295   m_outputProtoVertices(ctx, makeProtoVertices(inputTracks, vertices));
0296
0297   // store found vertices
0298   m_outputVertices(ctx, std::move(vertices));
0299
0300   return ProcessCode::SUCCESS;
0301 }
0302
0303 }  // namespace ActsExamples