Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020-2021 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 http://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/EigenStepper.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/Framework/AlgorithmContext.hpp"
0029 #include "ActsExamples/Framework/ProcessCode.hpp"
0030
0031 #include <memory>
0032 #include <optional>
0033 #include <ostream>
0034 #include <stdexcept>
0035 #include <system_error>
0036
0037 #include "VertexingHelpers.hpp"
0038
0039 ActsExamples::AdaptiveMultiVertexFinderAlgorithm::
0040     AdaptiveMultiVertexFinderAlgorithm(const Config& config,
0041                                        Acts::Logging::Level level)
0042     : ActsExamples::IAlgorithm("AdaptiveMultiVertexFinder", level),
0043       m_cfg(config),
0044       m_propagator{[&]() {
0045         // Set up EigenStepper
0046         Acts::EigenStepper<> stepper(m_cfg.bField);
0047
0048         // Set up the propagator
0049         return std::make_shared<Propagator>(stepper);
0050       }()},
0051       m_ipEstimator{[&]() {
0052         // Set up ImpactPointEstimator
0053         Acts::ImpactPointEstimator::Config ipEstimatorCfg(m_cfg.bField,
0054                                                           m_propagator);
0055         return Acts::ImpactPointEstimator(
0056             ipEstimatorCfg, logger().cloneWithSuffix("ImpactPointEstimator"));
0057       }()},
0058       m_linearizer{[&] {
0059         // Set up the helical track linearizer
0060         Linearizer::Config ltConfig;
0061         ltConfig.bField = m_cfg.bField;
0062         ltConfig.propagator = m_propagator;
0063         return Linearizer(ltConfig,
0064                           logger().cloneWithSuffix("HelicalTrackLinearizer"));
0065       }()},
0066       m_vertexFinder{makeVertexFinder()} {
0067   if (m_cfg.inputTrackParameters.empty()) {
0068     throw std::invalid_argument("Missing input track parameter collection");
0069   }
0070   if (m_cfg.outputProtoVertices.empty()) {
0071     throw std::invalid_argument("Missing output proto vertices collection");
0072   }
0073   if (m_cfg.outputVertices.empty()) {
0074     throw std::invalid_argument("Missing output vertices collection");
0075   }
0076
0077   m_inputTrackParameters.initialize(m_cfg.inputTrackParameters);
0078   m_outputProtoVertices.initialize(m_cfg.outputProtoVertices);
0079   m_outputVertices.initialize(m_cfg.outputVertices);
0080 }
0081
0082 auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const
0083     -> Acts::AdaptiveMultiVertexFinder {
0084   std::shared_ptr<const Acts::IVertexFinder> seedFinder;
0085   if (m_cfg.seedFinder == SeedFinder::GaussianSeeder) {
0086     using Seeder = Acts::TrackDensityVertexFinder;
0087     Acts::GaussianTrackDensity::Config trkDensityCfg;
0088     trkDensityCfg.extractParameters
0089         .connect<&Acts::InputTrack::extractParameters>();
0090     seedFinder = std::make_shared<Seeder>(Seeder::Config{trkDensityCfg});
0091   } else if (m_cfg.seedFinder == SeedFinder::AdaptiveGridSeeder) {
0092     // Set up track density used during vertex seeding
0093     Acts::AdaptiveGridTrackDensity::Config trkDensityCfg;
0094     // Bin extent in z-direction
0095     trkDensityCfg.spatialBinExtent = 15. * Acts::UnitConstants::um;
0096     // Bin extent in t-direction
0097     trkDensityCfg.temporalBinExtent = 19. * Acts::UnitConstants::mm;
0098     trkDensityCfg.useTime = m_cfg.useTime;
0099     Acts::AdaptiveGridTrackDensity trkDensity(trkDensityCfg);
0100
0101     // Set up vertex seeder and finder
0102     using Seeder = Acts::AdaptiveGridDensityVertexFinder;
0103     Seeder::Config seederConfig(trkDensity);
0104     seederConfig.extractParameters
0105         .connect<&Acts::InputTrack::extractParameters>();
0106     seedFinder = std::make_shared<Seeder>(seederConfig);
0107   } else {
0108     throw std::invalid_argument("Unknown seed finder");
0109   }
0110
0111   // Set up deterministic annealing with user-defined temperatures
0112   Acts::AnnealingUtility::Config annealingConfig;
0113   annealingConfig.setOfTemperatures = {1.};
0114   Acts::AnnealingUtility annealingUtility(annealingConfig);
0115
0116   // Set up the vertex fitter with user-defined annealing
0117   Fitter::Config fitterCfg(m_ipEstimator);
0118   fitterCfg.annealingTool = annealingUtility;
0119   fitterCfg.minWeight = 0.001;
0120   fitterCfg.doSmoothing = true;
0121   fitterCfg.useTime = m_cfg.useTime;
0122   fitterCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
0123   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&m_linearizer);
0124   Fitter fitter(std::move(fitterCfg),
0125                 logger().cloneWithSuffix("AdaptiveMultiVertexFitter"));
0126
0127   Acts::AdaptiveMultiVertexFinder::Config finderConfig(
0128       std::move(fitter), seedFinder, m_ipEstimator, m_cfg.bField);
0129   // Set the initial variance of the 4D vertex position. Since time is on a
0130   // numerical scale, we have to provide a greater value in the corresponding
0131   // dimension.
0132   finderConfig.initialVariances << 1e+2, 1e+2, 1e+2, 1e+8;
0133   finderConfig.tracksMaxZinterval = 1. * Acts::UnitConstants::mm;
0134   finderConfig.maxIterations = 200;
0135   finderConfig.useTime = m_cfg.useTime;
0136   // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
0137   finderConfig.tracksMaxSignificance = 5;
0138   // This should be used consistently with and without time
0139   finderConfig.doFullSplitting = false;
0140   // 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
0141   finderConfig.maxMergeVertexSignificance = 3;
0142   if (m_cfg.useTime) {
0143     // When using time, we have an extra contribution to the chi2 by the time
0144     // coordinate. We thus need to increase tracksMaxSignificance (i.e., the
0145     // maximum chi2 that a track can have to be associated with a vertex).
0146     // Using the same p-value for 3 dof instead of 2.
0147     // 6.7 corresponds to a p-value of ~0.92 using `chi2(x=6.7,ndf=3)`
0148     finderConfig.tracksMaxSignificance = 6.7;
0149     // Using the same p-value for 2 dof instead of 1.
0150     // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
0151     finderConfig.maxMergeVertexSignificance = 5;
0152   }
0153   finderConfig.extractParameters
0154       .template connect<&Acts::InputTrack::extractParameters>();
0155
0156   // Instantiate the finder
0157   return Acts::AdaptiveMultiVertexFinder(std::move(finderConfig),
0158                                          logger().clone());
0159 }
0160
0161 ActsExamples::ProcessCode
0162 ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
0163     const ActsExamples::AlgorithmContext& ctx) const {
0164   // retrieve input tracks and convert into the expected format
0165
0166   const auto& inputTrackParameters = m_inputTrackParameters(ctx);
0167   // TODO change this from pointers to tracks parameters to actual tracks
0168   auto inputTracks = makeInputTracks(inputTrackParameters);
0169
0170   if (inputTrackParameters.size() != inputTracks.size()) {
0171     ACTS_ERROR("Input track containers do not align: "
0172                << inputTrackParameters.size() << " != " << inputTracks.size());
0173   }
0174
0175   for (const auto& trk : inputTrackParameters) {
0176     if (trk.covariance() && trk.covariance()->determinant() <= 0) {
0177       // actually we should consider this as an error but I do not want the CI
0178       // to fail
0179       ACTS_WARNING("input track " << trk << " has det(cov) = "
0180                                   << trk.covariance()->determinant());
0181     }
0182   }
0183
0184   //////////////////////////////////////////////
0185   /* Full tutorial example code for reference */
0186   //////////////////////////////////////////////
0187
0188   // The vertex finder state
0189   auto state = m_vertexFinder.makeState(ctx.magFieldContext);
0190
0191   // Default vertexing options, this is where e.g. a constraint could be set
0192   Options finderOpts(ctx.geoContext, ctx.magFieldContext);
0193
0194   VertexCollection vertices;
0195
0196   if (inputTrackParameters.empty()) {
0197     ACTS_DEBUG("Empty track parameter collection found, skipping vertexing");
0198   } else {
0199     ACTS_DEBUG("Have " << inputTrackParameters.size()
0200                        << " input track parameters, running vertexing");
0201     // find vertices
0202     auto result = m_vertexFinder.find(inputTracks, finderOpts, state);
0203
0204     if (result.ok()) {
0205       vertices = std::move(result.value());
0206     } else {
0207       ACTS_ERROR("Error in vertex finder: " << result.error().message());
0208     }
0209   }
0210
0211   // show some debug output
0212   ACTS_INFO("Found " << vertices.size() << " vertices in event");
0213   for (const auto& vtx : vertices) {
0214     ACTS_DEBUG("Found vertex at " << vtx.fullPosition().transpose() << " with "
0215                                   << vtx.tracks().size() << " tracks.");
0216   }
0217
0218   // store proto vertices extracted from the found vertices
0219   m_outputProtoVertices(ctx, makeProtoVertices(inputTracks, vertices));
0220
0221   // store found vertices
0222   m_outputVertices(ctx, std::move(vertices));
0223
0224   return ActsExamples::ProcessCode::SUCCESS;
0225 }