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