Core/Vertexing/AdaptiveMultiVertexFinderTests.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 <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Geometry/GeometryIdentifier.hpp"
0019 #include "Acts/MagneticField/ConstantBField.hpp"
0020 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0021 #include "Acts/Propagator/EigenStepper.hpp"
0022 #include "Acts/Propagator/Propagator.hpp"
0023 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0024 #include "Acts/Utilities/AnnealingUtility.hpp"
0025 #include "Acts/Utilities/Helpers.hpp"
0026 #include "Acts/Utilities/Result.hpp"
0027 #include "Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp"
0028 #include "Acts/Vertexing/AdaptiveMultiVertexFinder.hpp"
0029 #include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp"
0030 #include "Acts/Vertexing/GaussianTrackDensity.hpp"
0031 #include "Acts/Vertexing/GridDensityVertexFinder.hpp"
0032 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
0033 #include "Acts/Vertexing/IVertexFinder.hpp"
0034 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0035 #include "Acts/Vertexing/TrackAtVertex.hpp"
0036 #include "Acts/Vertexing/TrackDensityVertexFinder.hpp"
0037 #include "Acts/Vertexing/Vertex.hpp"
0038 #include "Acts/Vertexing/VertexingOptions.hpp"
0039
0040 #include <algorithm>
0041 #include <array>
0042 #include <chrono>
0043 #include <cmath>
0044 #include <functional>
0045 #include <iostream>
0046 #include <map>
0047 #include <memory>
0048 #include <numbers>
0049 #include <string>
0050 #include <system_error>
0051 #include <tuple>
0052 #include <utility>
0053 #include <vector>
0054
0055 #include "VertexingDataHelper.hpp"
0056
0057 using namespace Acts::UnitLiterals;
0058
0059 namespace Acts::Test {
0060
0061 using Covariance = BoundSquareMatrix;
0062 using Propagator = Acts::Propagator<EigenStepper<>>;
0063 using Linearizer = HelicalTrackLinearizer;
0064
0065 // Create a test context
0066 GeometryContext geoContext = GeometryContext();
0067 MagneticFieldContext magFieldContext = MagneticFieldContext();
0068
0069 const std::string toolString = "AMVF";
0070
0071 /// @brief AMVF test with Gaussian seed finder
0072 BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_test) {
0073   // Set debug mode
0074   bool debugMode = false;
0075   // Set up constant B-Field
0076   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 2_T));
0077
0078   // Set up EigenStepper
0079   // EigenStepper<> stepper(bField);
0080   EigenStepper<> stepper(bField);
0081
0082   // Set up propagator with void navigator
0083   auto propagator = std::make_shared<Propagator>(stepper);
0084
0085   // IP 3D Estimator
0086   ImpactPointEstimator::Config ipEstimatorCfg(bField, propagator);
0087   ImpactPointEstimator ipEstimator(ipEstimatorCfg);
0088
0089   std::vector<double> temperatures{
0090       8., 4., 2., std::numbers::sqrt2, std::sqrt(3. / 2.), 1.};
0091   AnnealingUtility::Config annealingConfig;
0092   annealingConfig.setOfTemperatures = temperatures;
0093   AnnealingUtility annealingUtility(annealingConfig);
0094
0095   using Fitter = AdaptiveMultiVertexFitter;
0096
0097   Fitter::Config fitterCfg(ipEstimator);
0098
0099   fitterCfg.annealingTool = annealingUtility;
0100
0101   // Linearizer for BoundTrackParameters type test
0102   Linearizer::Config ltConfig;
0103   ltConfig.bField = bField;
0104   ltConfig.propagator = propagator;
0105   Linearizer linearizer(ltConfig);
0106
0107   // Test smoothing
0108   fitterCfg.doSmoothing = true;
0109   fitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0110   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0111
0112   Fitter fitter(fitterCfg);
0113
0114   GaussianTrackDensity::Config densityCfg;
0115   densityCfg.extractParameters.connect<&InputTrack::extractParameters>();
0116   auto seedFinder = std::make_shared<TrackDensityVertexFinder>(
0117       TrackDensityVertexFinder::Config{densityCfg});
0118
0119   AdaptiveMultiVertexFinder::Config finderConfig(std::move(fitter), seedFinder,
0120                                                  ipEstimator, bField);
0121   finderConfig.extractParameters.connect<&InputTrack::extractParameters>();
0122
0123   AdaptiveMultiVertexFinder finder(std::move(finderConfig));
0124   IVertexFinder::State state = finder.makeState(magFieldContext);
0125
0126   auto csvData = readTracksAndVertexCSV(toolString);
0127   std::vector<BoundTrackParameters> tracks = std::get<TracksData>(csvData);
0128
0129   if (debugMode) {
0130     std::cout << "Number of tracks in event: " << tracks.size() << std::endl;
0131     int maxCout = 10;
0132     int count = 0;
0133     for (const auto& trk : tracks) {
0134       std::cout << count << ". track: " << std::endl;
0135       std::cout << "params: " << trk << std::endl;
0136       count++;
0137       if (count == maxCout) {
0138         break;
0139       }
0140     }
0141   }
0142
0143   std::vector<InputTrack> inputTracks;
0144   for (const auto& trk : tracks) {
0145     inputTracks.emplace_back(&trk);
0146   }
0147
0148   // TODO: test without using beam spot constraint
0149   Vertex bsConstr = std::get<BeamSpotData>(csvData);
0150   VertexingOptions vertexingOptions(geoContext, magFieldContext, bsConstr);
0151
0152   auto t1 = std::chrono::system_clock::now();
0153   auto findResult = finder.find(inputTracks, vertexingOptions, state);
0154   auto t2 = std::chrono::system_clock::now();
0155
0156   auto timediff =
0157       std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
0158
0159   if (!findResult.ok()) {
0160     std::cout << findResult.error().message() << std::endl;
0161   }
0162
0163   BOOST_CHECK(findResult.ok());
0164
0165   std::vector<Vertex> allVertices = *findResult;
0166
0167   if (debugMode) {
0168     std::cout << "Time needed: " << timediff << " ms." << std::endl;
0169     std::cout << "Number of vertices reconstructed: " << allVertices.size()
0170               << std::endl;
0171
0172     int count = 0;
0173     for (const auto& vtx : allVertices) {
0174       count++;
0175       std::cout << count << ". Vertex at position: " << vtx.position()[0]
0176                 << ", " << vtx.position()[1] << ", " << vtx.position()[2]
0177                 << std::endl;
0178       std::cout << count << ". Vertex with cov: " << vtx.covariance()
0179                 << std::endl;
0180       std::cout << "\t with n tracks: " << vtx.tracks().size() << std::endl;
0181     }
0182   }
0183
0184   // Test expected outcomes from athena implementation
0185   // Number of reconstructed vertices
0186   auto verticesInfo = std::get<VerticesData>(csvData);
0187   const int expNRecoVertices = verticesInfo.size();
0188
0189   BOOST_CHECK_EQUAL(allVertices.size(), expNRecoVertices);
0190
0191   double relTol = 1e-2;
0192   double small = 1e-3;
0193   for (int i = 0; i < expNRecoVertices; i++) {
0194     auto recoVtx = allVertices[i];
0195     auto expVtx = verticesInfo[i];
0196     CHECK_CLOSE_OR_SMALL(recoVtx.position(), expVtx.position, relTol, small);
0197     CHECK_CLOSE_OR_SMALL(recoVtx.covariance(), expVtx.covariance, relTol,
0198                          small);
0199     BOOST_CHECK_EQUAL(recoVtx.tracks().size(), expVtx.nTracks);
0200     CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].trackWeight, expVtx.trk1Weight,
0201                          relTol, small);
0202     CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].vertexCompatibility,
0203                          expVtx.trk1Comp, relTol, small);
0204   }
0205 }
0206
0207 // Dummy user-defined InputTrackStub type
0208 struct InputTrackStub {
0209   InputTrackStub(const BoundTrackParameters& params, int id)
0210       : m_parameters(params), m_id(id) {}
0211
0212   const BoundTrackParameters& parameters() const { return m_parameters; }
0213   // store e.g. link to original objects here
0214
0215   int id() const { return m_id; }
0216
0217  private:
0218   BoundTrackParameters m_parameters;
0219
0220   // Some test track ID
0221   int m_id;
0222 };
0223
0224 /// @brief AMVF test with user-defined input track type
0225 BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_usertype_test) {
0226   // Set debug mode
0227   bool debugMode = false;
0228   // Set up constant B-Field
0229   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 2_T));
0230
0231   // Set up EigenStepper
0232   // EigenStepper<> stepper(bField);
0233   EigenStepper<> stepper(bField);
0234
0235   // Set up propagator with void navigator
0236   auto propagator = std::make_shared<Propagator>(stepper);
0237
0238   // Create a custom std::function to extract BoundTrackParameters from
0239   // user-defined InputTrackStub
0240   auto extractParameters = [](const InputTrack& track) {
0241     return track.as<InputTrackStub>()->parameters();
0242   };
0243
0244   // IP 3D Estimator
0245   ImpactPointEstimator::Config ipEstimatorCfg(bField, propagator);
0246   ImpactPointEstimator ipEstimator(ipEstimatorCfg);
0247
0248   std::vector<double> temperatures{
0249       8., 4., 2., std::numbers::sqrt2, std::sqrt(3. / 2.), 1.};
0250   AnnealingUtility::Config annealingConfig;
0251   annealingConfig.setOfTemperatures = temperatures;
0252   AnnealingUtility annealingUtility(annealingConfig);
0253
0254   using Fitter = AdaptiveMultiVertexFitter;
0255
0256   Fitter::Config fitterCfg(ipEstimator);
0257
0258   fitterCfg.annealingTool = annealingUtility;
0259
0260   // Linearizer
0261   Linearizer::Config ltConfig;
0262   ltConfig.bField = bField;
0263   ltConfig.propagator = propagator;
0264   Linearizer linearizer(ltConfig);
0265
0266   // Test smoothing
0267   fitterCfg.doSmoothing = true;
0268   fitterCfg.extractParameters.connect(extractParameters);
0269   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0270
0271   Fitter fitter(fitterCfg);
0272
0273   GaussianTrackDensity::Config densityCfg;
0274   densityCfg.extractParameters.connect(extractParameters);
0275   auto seedFinder = std::make_shared<TrackDensityVertexFinder>(
0276       TrackDensityVertexFinder::Config{densityCfg});
0277
0278   AdaptiveMultiVertexFinder::Config finderConfig(
0279       std::move(fitter), std::move(seedFinder), ipEstimator, bField);
0280   finderConfig.extractParameters.connect(extractParameters);
0281
0282   AdaptiveMultiVertexFinder finder(std::move(finderConfig));
0283   IVertexFinder::State state = finder.makeState(magFieldContext);
0284
0285   auto csvData = readTracksAndVertexCSV(toolString);
0286   auto tracks = std::get<TracksData>(csvData);
0287
0288   std::vector<InputTrackStub> userTracks;
0289   int idCount = 0;
0290   for (const auto& trk : tracks) {
0291     userTracks.push_back(InputTrackStub(trk, idCount));
0292     idCount++;
0293   }
0294
0295   if (debugMode) {
0296     std::cout << "Number of tracks in event: " << tracks.size() << std::endl;
0297     int maxCout = 10;
0298     int count = 0;
0299     for (const auto& trk : tracks) {
0300       std::cout << count << ". track: " << std::endl;
0301       std::cout << "params: " << trk << std::endl;
0302       count++;
0303       if (count == maxCout) {
0304         break;
0305       }
0306     }
0307   }
0308
0309   std::vector<InputTrack> userInputTracks;
0310   for (const auto& trk : userTracks) {
0311     userInputTracks.emplace_back(&trk);
0312   }
0313
0314   Vertex constraintVtx;
0315   constraintVtx.setPosition(std::get<BeamSpotData>(csvData).position());
0316   constraintVtx.setCovariance(std::get<BeamSpotData>(csvData).covariance());
0317
0318   VertexingOptions vertexingOptions(geoContext, magFieldContext, constraintVtx);
0319
0320   auto findResult = finder.find(userInputTracks, vertexingOptions, state);
0321
0322   if (!findResult.ok()) {
0323     std::cout << findResult.error().message() << std::endl;
0324   }
0325
0326   BOOST_CHECK(findResult.ok());
0327
0328   std::vector<Vertex> allVertices = *findResult;
0329
0330   if (debugMode) {
0331     std::cout << "Number of vertices reconstructed: " << allVertices.size()
0332               << std::endl;
0333
0334     int count = 0;
0335     for (const auto& vtx : allVertices) {
0336       count++;
0337       std::cout << count << ". Vertex at position: " << vtx.position()[0]
0338                 << ", " << vtx.position()[1] << ", " << vtx.position()[2]
0339                 << std::endl;
0340       std::cout << count << ". Vertex with cov: " << vtx.covariance()
0341                 << std::endl;
0342       std::cout << "\t with n tracks: " << vtx.tracks().size() << std::endl;
0343     }
0344     for (auto& trk : allVertices[0].tracks()) {
0345       std::cout << "Track ID at first vertex: "
0346                 << trk.originalParams.as<InputTrackStub>()->id() << std::endl;
0347     }
0348   }
0349
0350   auto verticesInfo = std::get<VerticesData>(csvData);
0351   const int expNRecoVertices = verticesInfo.size();
0352
0353   BOOST_CHECK_EQUAL(allVertices.size(), expNRecoVertices);
0354
0355   double relTol = 1e-2;
0356   double small = 1e-3;
0357   for (int i = 0; i < expNRecoVertices; i++) {
0358     auto recoVtx = allVertices[i];
0359     auto expVtx = verticesInfo[i];
0360     CHECK_CLOSE_OR_SMALL(recoVtx.position(), expVtx.position, relTol, small);
0361     CHECK_CLOSE_OR_SMALL(recoVtx.covariance(), expVtx.covariance, relTol,
0362                          small);
0363     BOOST_CHECK_EQUAL(recoVtx.tracks().size(), expVtx.nTracks);
0364     CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].trackWeight, expVtx.trk1Weight,
0365                          relTol, small);
0366     CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].vertexCompatibility,
0367                          expVtx.trk1Comp, relTol, small);
0368   }
0369 }
0370
0371 /// @brief AMVF test with grid seed finder
0372 BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_grid_seed_finder_test) {
0373   // Set debug mode
0374   bool debugMode = false;
0375   if (debugMode) {
0376     std::cout << "Starting AMVF test with grid seed finder..." << std::endl;
0377   }
0378   // Set up constant B-Field
0379   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 2_T));
0380
0381   // Set up EigenStepper
0382   // EigenStepper<> stepper(bField);
0383   EigenStepper<> stepper(bField);
0384
0385   // Set up propagator with void navigator
0386   auto propagator = std::make_shared<Propagator>(stepper);
0387
0388   // IP Estimator
0389   ImpactPointEstimator::Config ipEstCfg(bField, propagator);
0390   ImpactPointEstimator ipEst(ipEstCfg);
0391
0392   std::vector<double> temperatures{
0393       8., 4., 2., std::numbers::sqrt2, std::sqrt(3. / 2.), 1.};
0394   AnnealingUtility::Config annealingConfig;
0395   annealingConfig.setOfTemperatures = temperatures;
0396   AnnealingUtility annealingUtility(annealingConfig);
0397
0398   using Fitter = AdaptiveMultiVertexFitter;
0399
0400   Fitter::Config fitterCfg(ipEst);
0401
0402   fitterCfg.annealingTool = annealingUtility;
0403
0404   // Linearizer for BoundTrackParameters type test
0405   Linearizer::Config ltConfig;
0406   ltConfig.bField = bField;
0407   ltConfig.propagator = propagator;
0408   Linearizer linearizer(ltConfig);
0409
0410   // Test smoothing
0411   fitterCfg.doSmoothing = true;
0412   fitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0413   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0414
0415   Fitter fitter(fitterCfg);
0416
0417   using SeedFinder = GridDensityVertexFinder;
0418   SeedFinder::Config seedFinderCfg{{{250, 4000, 55}}};
0419   seedFinderCfg.cacheGridStateForTrackRemoval = true;
0420   seedFinderCfg.extractParameters.connect<&InputTrack::extractParameters>();
0421
0422   auto seedFinder = std::make_shared<SeedFinder>(seedFinderCfg);
0423
0424   AdaptiveMultiVertexFinder::Config finderConfig(
0425       std::move(fitter), std::move(seedFinder), ipEst, bField);
0426   finderConfig.extractParameters.connect<&InputTrack::extractParameters>();
0427
0428   AdaptiveMultiVertexFinder finder(std::move(finderConfig));
0429   IVertexFinder::State state = finder.makeState(magFieldContext);
0430
0431   auto csvData = readTracksAndVertexCSV(toolString);
0432   auto tracks = std::get<TracksData>(csvData);
0433
0434   if (debugMode) {
0435     std::cout << "Number of tracks in event: " << tracks.size() << std::endl;
0436     int maxCout = 10;
0437     int count = 0;
0438     for (const auto& trk : tracks) {
0439       std::cout << count << ". track: " << std::endl;
0440       std::cout << "params: " << trk << std::endl;
0441       count++;
0442       if (count == maxCout) {
0443         break;
0444       }
0445     }
0446   }
0447
0448   std::vector<InputTrack> inputTracks;
0449   for (const auto& trk : tracks) {
0450     inputTracks.emplace_back(&trk);
0451   }
0452
0453   // TODO: test using beam spot constraint
0454   Vertex bsConstr = std::get<BeamSpotData>(csvData);
0455   VertexingOptions vertexingOptions(geoContext, magFieldContext, bsConstr);
0456
0457   auto t1 = std::chrono::system_clock::now();
0458   auto findResult = finder.find(inputTracks, vertexingOptions, state);
0459   auto t2 = std::chrono::system_clock::now();
0460
0461   auto timediff =
0462       std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
0463
0464   if (!findResult.ok()) {
0465     std::cout << findResult.error().message() << std::endl;
0466   }
0467
0468   BOOST_CHECK(findResult.ok());
0469
0470   std::vector<Vertex> allVertices = *findResult;
0471
0472   if (debugMode) {
0473     std::cout << "Time needed: " << timediff << " ms." << std::endl;
0474     std::cout << "Number of vertices reconstructed: " << allVertices.size()
0475               << std::endl;
0476
0477     int count = 0;
0478     for (const auto& vtx : allVertices) {
0479       count++;
0480       std::cout << count << ". Vertex at position: " << vtx.position()[0]
0481                 << ", " << vtx.position()[1] << ", " << vtx.position()[2]
0482                 << std::endl;
0483       std::cout << count << ". Vertex with cov: " << vtx.covariance()
0484                 << std::endl;
0485       std::cout << "\t with n tracks: " << vtx.tracks().size() << std::endl;
0486     }
0487   }
0488   // Test expected outcomes from athena implementation
0489   // Number of reconstructed vertices
0490   auto verticesInfo = std::get<VerticesData>(csvData);
0491   const int expNRecoVertices = verticesInfo.size();
0492
0493   BOOST_CHECK_EQUAL(allVertices.size(), expNRecoVertices);
0494   std::vector<bool> vtxFound(expNRecoVertices, false);
0495
0496   for (const auto& vtx : allVertices) {
0497     double vtxZ = vtx.position()[2];
0498     double diffZ = 1e5;
0499     int foundVtxIdx = -1;
0500     for (int i = 0; i < expNRecoVertices; i++) {
0501       if (!vtxFound[i]) {
0502         if (std::abs(vtxZ - verticesInfo[i].position[2]) < diffZ) {
0503           diffZ = std::abs(vtxZ - verticesInfo[i].position[2]);
0504           foundVtxIdx = i;
0505         }
0506       }
0507     }
0508     if (diffZ < 0.5_mm) {
0509       vtxFound[foundVtxIdx] = true;
0510       CHECK_CLOSE_ABS(vtx.tracks().size(), verticesInfo[foundVtxIdx].nTracks,
0511                       1);
0512     }
0513   }
0514   for (bool found : vtxFound) {
0515     BOOST_CHECK_EQUAL(found, true);
0516   }
0517 }
0518
0519 /// @brief AMVF test with adaptive grid seed finder
0520 BOOST_AUTO_TEST_CASE(
0521     adaptive_multi_vertex_finder_adaptive_grid_seed_finder_test) {
0522   // Set debug mode
0523   bool debugMode = false;
0524   if (debugMode) {
0525     std::cout << "Starting AMVF test with adaptive grid seed finder..."
0526               << std::endl;
0527   }
0528   // Set up constant B-Field
0529   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 2_T));
0530
0531   // Set up EigenStepper
0532   // EigenStepper<> stepper(bField);
0533   EigenStepper<> stepper(bField);
0534
0535   // Set up propagator with void navigator
0536   auto propagator = std::make_shared<Propagator>(stepper);
0537
0538   // IP Estimator
0539   ImpactPointEstimator::Config ipEstCfg(bField, propagator);
0540   ImpactPointEstimator ipEst(ipEstCfg);
0541
0542   std::vector<double> temperatures{
0543       8., 4., 2., std::numbers::sqrt2, std::sqrt(3. / 2.), 1.};
0544   AnnealingUtility::Config annealingConfig;
0545   annealingConfig.setOfTemperatures = temperatures;
0546   AnnealingUtility annealingUtility(annealingConfig);
0547
0548   using Fitter = AdaptiveMultiVertexFitter;
0549
0550   Fitter::Config fitterCfg(ipEst);
0551
0552   fitterCfg.annealingTool = annealingUtility;
0553
0554   // Linearizer for BoundTrackParameters type test
0555   Linearizer::Config ltConfig;
0556   ltConfig.bField = bField;
0557   ltConfig.propagator = propagator;
0558   Linearizer linearizer(ltConfig);
0559
0560   // Test smoothing
0561   fitterCfg.doSmoothing = true;
0562   fitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0563   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0564
0565   Fitter fitter(fitterCfg);
0566
0567   // Grid density used during vertex seed finding
0568   AdaptiveGridTrackDensity::Config gridDensityCfg;
0569   // force track to have exactly spatialTrkGridSize spatial bins for testing
0570   // purposes
0571   gridDensityCfg.spatialTrkGridSizeRange = {55, 55};
0572   gridDensityCfg.spatialBinExtent = 0.05;
0573   AdaptiveGridTrackDensity gridDensity(gridDensityCfg);
0574
0575   using SeedFinder = AdaptiveGridDensityVertexFinder;
0576   SeedFinder::Config seedFinderCfg(gridDensity);
0577   seedFinderCfg.cacheGridStateForTrackRemoval = true;
0578   seedFinderCfg.extractParameters.connect<&InputTrack::extractParameters>();
0579
0580   auto seedFinder = std::make_shared<SeedFinder>(seedFinderCfg);
0581
0582   AdaptiveMultiVertexFinder::Config finderConfig(
0583       std::move(fitter), std::move(seedFinder), ipEst, bField);
0584   finderConfig.extractParameters.connect<&InputTrack::extractParameters>();
0585
0586   AdaptiveMultiVertexFinder finder(std::move(finderConfig));
0587   IVertexFinder::State state = finder.makeState(magFieldContext);
0588
0589   auto csvData = readTracksAndVertexCSV(toolString);
0590   auto tracks = std::get<TracksData>(csvData);
0591
0592   if (debugMode) {
0593     std::cout << "Number of tracks in event: " << tracks.size() << std::endl;
0594     int maxCout = 10;
0595     int count = 0;
0596     for (const auto& trk : tracks) {
0597       std::cout << count << ". track: " << std::endl;
0598       std::cout << "params: " << trk << std::endl;
0599       count++;
0600       if (count == maxCout) {
0601         break;
0602       }
0603     }
0604   }
0605
0606   std::vector<InputTrack> inputTracks;
0607   for (const auto& trk : tracks) {
0608     inputTracks.emplace_back(&trk);
0609   }
0610
0611   Vertex bsConstr = std::get<BeamSpotData>(csvData);
0612   VertexingOptions vertexingOptions(geoContext, magFieldContext, bsConstr);
0613
0614   auto t1 = std::chrono::system_clock::now();
0615   auto findResult = finder.find(inputTracks, vertexingOptions, state);
0616   auto t2 = std::chrono::system_clock::now();
0617
0618   auto timediff =
0619       std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
0620
0621   if (!findResult.ok()) {
0622     std::cout << findResult.error().message() << std::endl;
0623   }
0624
0625   BOOST_CHECK(findResult.ok());
0626
0627   std::vector<Vertex> allVertices = *findResult;
0628
0629   if (debugMode) {
0630     std::cout << "Time needed: " << timediff << " ms." << std::endl;
0631     std::cout << "Number of vertices reconstructed: " << allVertices.size()
0632               << std::endl;
0633
0634     int count = 0;
0635     for (const auto& vtx : allVertices) {
0636       count++;
0637       std::cout << count << ". Vertex at position: " << vtx.position()[0]
0638                 << ", " << vtx.position()[1] << ", " << vtx.position()[2]
0639                 << std::endl;
0640       std::cout << count << ". Vertex with cov: " << vtx.covariance()
0641                 << std::endl;
0642       std::cout << "\t with n tracks: " << vtx.tracks().size() << std::endl;
0643     }
0644   }
0645   // Test expected outcomes from athena implementation
0646   // Number of reconstructed vertices
0647   auto verticesInfo = std::get<VerticesData>(csvData);
0648   const int expNRecoVertices = verticesInfo.size();
0649
0650   BOOST_CHECK_EQUAL(allVertices.size(), expNRecoVertices);
0651   std::vector<bool> vtxFound(expNRecoVertices, false);
0652
0653   for (const auto& vtx : allVertices) {
0654     double vtxZ = vtx.position()[2];
0655     double diffZ = 1e5;
0656     int foundVtxIdx = -1;
0657     for (int i = 0; i < expNRecoVertices; i++) {
0658       if (!vtxFound[i]) {
0659         if (std::abs(vtxZ - verticesInfo[i].position[2]) < diffZ) {
0660           diffZ = std::abs(vtxZ - verticesInfo[i].position[2]);
0661           foundVtxIdx = i;
0662         }
0663       }
0664     }
0665     if (diffZ < 0.5_mm) {
0666       vtxFound[foundVtxIdx] = true;
0667       CHECK_CLOSE_ABS(vtx.tracks().size(), verticesInfo[foundVtxIdx].nTracks,
0668                       2);
0669     }
0670   }
0671   for (bool found : vtxFound) {
0672     BOOST_CHECK_EQUAL(found, true);
0673   }
0674 }
0675
0676 }  // namespace Acts::Test