Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-14 07:42:08

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