Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-14 08:02:31

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