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