File indexing completed on 2025-10-29 07:55:45
0001 
0002 
0003 
0004 
0005 
0006 
0007 
0008 
0009 #include <boost/test/unit_test.hpp>
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "Acts/Definitions/Direction.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/EventData/Charge.hpp"
0017 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0018 #include "Acts/EventData/TrackParameters.hpp"
0019 #include "Acts/Geometry/GeometryContext.hpp"
0020 #include "Acts/Geometry/GeometryIdentifier.hpp"
0021 #include "Acts/MagneticField/ConstantBField.hpp"
0022 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0023 #include "Acts/Propagator/EigenStepper.hpp"
0024 #include "Acts/Propagator/Propagator.hpp"
0025 #include "Acts/Surfaces/PerigeeSurface.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "Acts/Utilities/Result.hpp"
0028 #include "Acts/Vertexing/FullBilloirVertexFitter.hpp"
0029 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
0030 #include "Acts/Vertexing/IVertexFinder.hpp"
0031 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0032 #include "Acts/Vertexing/IterativeVertexFinder.hpp"
0033 #include "Acts/Vertexing/TrackAtVertex.hpp"
0034 #include "Acts/Vertexing/Vertex.hpp"
0035 #include "Acts/Vertexing/VertexingOptions.hpp"
0036 #include "Acts/Vertexing/ZScanVertexFinder.hpp"
0037 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0038 
0039 #include <algorithm>
0040 #include <array>
0041 #include <cmath>
0042 #include <functional>
0043 #include <iostream>
0044 #include <iterator>
0045 #include <memory>
0046 #include <numbers>
0047 #include <optional>
0048 #include <random>
0049 #include <string>
0050 #include <system_error>
0051 #include <tuple>
0052 #include <type_traits>
0053 #include <utility>
0054 #include <vector>
0055 
0056 #include "VertexingDataHelper.hpp"
0057 
0058 using namespace Acts;
0059 using namespace Acts::UnitLiterals;
0060 
0061 namespace ActsTests {
0062 
0063 using Covariance = BoundSquareMatrix;
0064 using Propagator = Acts::Propagator<EigenStepper<>>;
0065 using Linearizer = HelicalTrackLinearizer;
0066 
0067 
0068 GeometryContext geoContext = GeometryContext();
0069 MagneticFieldContext magFieldContext = MagneticFieldContext();
0070 
0071 const std::string toolString = "IVF";
0072 
0073 
0074 std::uniform_real_distribution<double> vXYDist(-0.1_mm, 0.1_mm);
0075 
0076 std::uniform_real_distribution<double> vZDist(-20_mm, 20_mm);
0077 
0078 std::uniform_real_distribution<double> d0Dist(-0.01_mm, 0.01_mm);
0079 
0080 std::uniform_real_distribution<double> z0Dist(-0.2_mm, 0.2_mm);
0081 
0082 std::uniform_real_distribution<double> pTDist(0.4_GeV, 10_GeV);
0083 
0084 std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0085                                                std::numbers::pi);
0086 
0087 std::uniform_real_distribution<double> thetaDist(1., std::numbers::pi - 1.);
0088 
0089 std::uniform_real_distribution<double> qDist(-1, 1);
0090 
0091 std::uniform_real_distribution<double> resIPDist(0., 100_um);
0092 
0093 std::uniform_real_distribution<double> resAngDist(0., 0.1);
0094 
0095 std::uniform_real_distribution<double> resQoPDist(-0.01, 0.01);
0096 
0097 std::uniform_int_distribution<std::uint32_t> nVertexDist(1, 6);
0098 
0099 std::uniform_int_distribution<std::uint32_t> nTracksDist(5, 15);
0100 
0101 
0102 struct InputTrackStub {
0103   explicit InputTrackStub(const BoundTrackParameters& params)
0104       : m_parameters(params) {}
0105 
0106   const BoundTrackParameters& parameters() const { return m_parameters; }
0107 
0108   
0109 
0110  private:
0111   BoundTrackParameters m_parameters;
0112 };
0113 
0114 BOOST_AUTO_TEST_SUITE(VertexingSuite)
0115 
0116 
0117 
0118 BOOST_AUTO_TEST_CASE(iterative_finder_test) {
0119   bool debug = false;
0120 
0121   
0122   int mySeed = 31415;
0123   std::mt19937 gen(mySeed);
0124 
0125   
0126   unsigned int nEvents = 5;  
0127 
0128   
0129   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0130 
0131   for (unsigned int iEvent = 0; iEvent < nEvents; ++iEvent) {
0132     
0133     EigenStepper<> stepper(bField);
0134 
0135     
0136     auto propagator = std::make_shared<Propagator>(stepper);
0137 
0138     
0139     Linearizer::Config ltConfig;
0140     ltConfig.bField = bField;
0141     ltConfig.propagator = propagator;
0142     Linearizer linearizer(ltConfig);
0143 
0144     using BilloirFitter = FullBilloirVertexFitter;
0145 
0146     
0147     BilloirFitter::Config vertexFitterCfg;
0148     vertexFitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0149     vertexFitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(
0150         &linearizer);
0151 
0152     BilloirFitter bFitter(vertexFitterCfg);
0153 
0154     ImpactPointEstimator::Config ipEstimatorCfg(bField, propagator);
0155     ImpactPointEstimator ipEstimator(ipEstimatorCfg);
0156 
0157     using ZScanSeedFinder = ZScanVertexFinder;
0158 
0159     ZScanSeedFinder::Config seedFinderCfg(ipEstimator);
0160     seedFinderCfg.extractParameters.connect<&InputTrack::extractParameters>();
0161 
0162     auto sFinder = std::make_shared<ZScanSeedFinder>(seedFinderCfg);
0163 
0164     
0165 
0166     IterativeVertexFinder::Config cfg(std::move(bFitter), std::move(sFinder),
0167                                       ipEstimator);
0168     cfg.field = bField;
0169     cfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0170 
0171     cfg.reassignTracksAfterFirstFit = true;
0172     cfg.extractParameters.connect<&InputTrack::extractParameters>();
0173 
0174     IterativeVertexFinder finder(std::move(cfg));
0175     IVertexFinder::State state{
0176         IterativeVertexFinder::State(*bField, magFieldContext)};
0177 
0178     
0179     std::vector<std::unique_ptr<const BoundTrackParameters>> tracks;
0180 
0181     std::vector<InputTrack> inputTracks;
0182 
0183     
0184     std::vector<Vertex> trueVertices;
0185 
0186     
0187     std::uint32_t nVertices = nVertexDist(gen);
0188     for (std::uint32_t iVertex = 0; iVertex < nVertices; ++iVertex) {
0189       
0190       std::uint32_t nTracks = nTracksDist(gen);
0191 
0192       if (debug) {
0193         std::cout << "Event " << iEvent << ", Vertex " << iVertex << "/"
0194                   << nVertices << " with " << nTracks << " tracks."
0195                   << std::endl;
0196       }
0197       
0198       std::shared_ptr<PerigeeSurface> perigeeSurface =
0199           Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0200 
0201       
0202       double x = vXYDist(gen);
0203       double y = vXYDist(gen);
0204       double z = vZDist(gen);
0205 
0206       
0207       Vertex trueV(Vector3(x, y, z));
0208       std::vector<TrackAtVertex> tracksAtTrueVtx;
0209 
0210       
0211       double d0_v = std::hypot(x, y);
0212       double z0_v = z;
0213 
0214       
0215       
0216       for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
0217         
0218         double q = qDist(gen) < 0 ? -1. : 1.;
0219 
0220         
0221         BoundVector paramVec;
0222         double z0track = z0_v + z0Dist(gen);
0223         paramVec << d0_v + d0Dist(gen), z0track, phiDist(gen), thetaDist(gen),
0224             q / pTDist(gen), 0.;
0225 
0226         
0227         double res_d0 = resIPDist(gen);
0228         double res_z0 = resIPDist(gen);
0229         double res_ph = resAngDist(gen);
0230         double res_th = resAngDist(gen);
0231         double res_qp = resQoPDist(gen);
0232 
0233         
0234         Covariance covMat;
0235         covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0.,
0236             0., 0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0.,
0237             res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0.,
0238             0., 0., 0., 0., 1.;
0239         auto params =
0240             BoundTrackParameters(perigeeSurface, paramVec, std::move(covMat),
0241                                  ParticleHypothesis::pion());
0242 
0243         tracks.push_back(std::make_unique<BoundTrackParameters>(params));
0244 
0245         TrackAtVertex trAtVt(0., params, InputTrack{tracks.back().get()});
0246         tracksAtTrueVtx.push_back(trAtVt);
0247       }
0248 
0249       trueV.setTracksAtVertex(tracksAtTrueVtx);
0250       trueVertices.push_back(trueV);
0251 
0252     }  
0253 
0254     
0255     std::shuffle(std::begin(tracks), std::end(tracks), gen);
0256 
0257     for (const auto& trk : tracks) {
0258       inputTracks.emplace_back(trk.get());
0259     }
0260 
0261     VertexingOptions vertexingOptions(geoContext, magFieldContext);
0262 
0263     
0264     auto res = finder.find(inputTracks, vertexingOptions, state);
0265 
0266     if (!res.ok()) {
0267       BOOST_FAIL(res.error().message());
0268     }
0269 
0270     
0271     auto vertexCollection = *res;
0272 
0273     
0274     CHECK_CLOSE_ABS(vertexCollection.size(), nVertices, 2);
0275 
0276     if (debug) {
0277       std::cout << "########## RESULT: ########## Event " << iEvent
0278                 << std::endl;
0279       std::cout << "Number of true vertices: " << nVertices << std::endl;
0280       std::cout << "Number of reco vertices: " << vertexCollection.size()
0281                 << std::endl;
0282 
0283       int count = 1;
0284       std::cout << "----- True vertices -----" << std::endl;
0285       for (const auto& vertex : trueVertices) {
0286         Vector3 pos = vertex.position();
0287         std::cout << count << ". True Vertex:\t Position:"
0288                   << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
0289                   << std::endl;
0290         std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
0291                   << std::endl;
0292         count++;
0293       }
0294       std::cout << "----- Reco vertices -----" << std::endl;
0295       count = 1;
0296       for (const auto& vertex : vertexCollection) {
0297         Vector3 pos = vertex.position();
0298         std::cout << count << ". Reco Vertex:\t Position:"
0299                   << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
0300                   << std::endl;
0301         std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
0302                   << std::endl;
0303         count++;
0304       }
0305     }
0306 
0307     
0308     bool allVerticesFound = true;
0309     for (const auto& trueVertex : trueVertices) {
0310       Vector4 truePos = trueVertex.fullPosition();
0311       bool currentVertexFound = false;
0312       for (const auto& recoVertex : vertexCollection) {
0313         Vector4 recoPos = recoVertex.fullPosition();
0314         
0315         double zDistance = std::abs(truePos[eZ] - recoPos[eZ]);
0316         if (zDistance < 2_mm) {
0317           currentVertexFound = true;
0318         }
0319       }
0320       if (!currentVertexFound) {
0321         allVerticesFound = false;
0322       }
0323     }
0324 
0325     
0326     BOOST_CHECK(allVerticesFound);
0327   }
0328 }
0329 
0330 
0331 
0332 
0333 
0334 BOOST_AUTO_TEST_CASE(iterative_finder_test_user_track_type) {
0335   bool debug = false;
0336 
0337   
0338   int mySeed = 31415;
0339   std::mt19937 gen(mySeed);
0340 
0341   
0342   unsigned int nEvents = 5;  
0343 
0344   
0345   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0346 
0347   for (unsigned int iEvent = 0; iEvent < nEvents; ++iEvent) {
0348     
0349     EigenStepper<> stepper(bField);
0350 
0351     
0352     auto propagator = std::make_shared<Propagator>(stepper);
0353 
0354     
0355     Linearizer::Config ltConfigUT;
0356     ltConfigUT.bField = bField;
0357     ltConfigUT.propagator = propagator;
0358     Linearizer linearizer(ltConfigUT);
0359 
0360     
0361     using BilloirFitter = FullBilloirVertexFitter;
0362 
0363     
0364     
0365     auto extractParameters = [](const InputTrack& params) {
0366       return params.as<InputTrackStub>()->parameters();
0367     };
0368 
0369     
0370     BilloirFitter::Config vertexFitterCfg;
0371     vertexFitterCfg.extractParameters.connect(extractParameters);
0372     vertexFitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(
0373         &linearizer);
0374 
0375     BilloirFitter bFitter(vertexFitterCfg);
0376 
0377     ImpactPointEstimator::Config ipEstimatorCfg(bField, propagator);
0378     ImpactPointEstimator ipEstimator(ipEstimatorCfg);
0379 
0380     using ZScanSeedFinder = ZScanVertexFinder;
0381     ZScanSeedFinder::Config seedFinderCfg(ipEstimator);
0382     seedFinderCfg.extractParameters.connect(extractParameters);
0383 
0384     auto sFinder = std::make_shared<ZScanSeedFinder>(seedFinderCfg);
0385 
0386     
0387     IterativeVertexFinder::Config cfg(std::move(bFitter), std::move(sFinder),
0388                                       ipEstimator);
0389     cfg.reassignTracksAfterFirstFit = true;
0390     cfg.extractParameters.connect(extractParameters);
0391     cfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0392     cfg.field = bField;
0393 
0394     IterativeVertexFinder finder(std::move(cfg));
0395     IVertexFinder::State state{
0396         IterativeVertexFinder::State(*bField, magFieldContext)};
0397 
0398     
0399     std::vector<std::unique_ptr<const InputTrackStub>> tracks;
0400 
0401     std::vector<InputTrack> inputTracks;
0402 
0403     
0404     std::vector<Vertex> trueVertices;
0405 
0406     
0407     std::uint32_t nVertices = nVertexDist(gen);
0408     for (std::uint32_t iVertex = 0; iVertex < nVertices; ++iVertex) {
0409       
0410       std::uint32_t nTracks = nTracksDist(gen);
0411 
0412       if (debug) {
0413         std::cout << "Event " << iEvent << ", Vertex " << iVertex << "/"
0414                   << nVertices << " with " << nTracks << " tracks."
0415                   << std::endl;
0416       }
0417       
0418       std::shared_ptr<PerigeeSurface> perigeeSurface =
0419           Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0420 
0421       
0422       double x = vXYDist(gen);
0423       double y = vXYDist(gen);
0424       double z = vZDist(gen);
0425 
0426       
0427       Vertex trueV(Vector3(x, y, z));
0428       std::vector<TrackAtVertex> tracksAtTrueVtx;
0429 
0430       
0431       double d0_v = std::hypot(x, y);
0432       double z0_v = z;
0433 
0434       
0435       
0436       for (std::uint32_t iTrack = 0; iTrack < nTracks; iTrack++) {
0437         
0438         double q = qDist(gen) < 0 ? -1. : 1.;
0439 
0440         
0441         BoundVector paramVec;
0442         double z0track = z0_v + z0Dist(gen);
0443         paramVec << d0_v + d0Dist(gen), z0track, phiDist(gen), thetaDist(gen),
0444             q / pTDist(gen), 0.;
0445 
0446         
0447         double res_d0 = resIPDist(gen);
0448         double res_z0 = resIPDist(gen);
0449         double res_ph = resAngDist(gen);
0450         double res_th = resAngDist(gen);
0451         double res_qp = resQoPDist(gen);
0452 
0453         
0454         Covariance covMat;
0455 
0456         covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0.,
0457             0., 0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0.,
0458             res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0.,
0459             0., 0., 0., 0., 1.;
0460         auto params =
0461             BoundTrackParameters(perigeeSurface, paramVec, std::move(covMat),
0462                                  ParticleHypothesis::pion());
0463 
0464         tracks.push_back(std::make_unique<InputTrackStub>(params));
0465 
0466         TrackAtVertex trAtVt(0., params, InputTrack{tracks.back().get()});
0467         tracksAtTrueVtx.push_back(trAtVt);
0468       }
0469 
0470       trueV.setTracksAtVertex(tracksAtTrueVtx);
0471       trueVertices.push_back(trueV);
0472 
0473     }  
0474 
0475     
0476     std::shuffle(std::begin(tracks), std::end(tracks), gen);
0477 
0478     for (const auto& trk : tracks) {
0479       inputTracks.emplace_back(trk.get());
0480     }
0481 
0482     VertexingOptions vertexingOptionsUT(geoContext, magFieldContext);
0483 
0484     
0485     auto res = finder.find(inputTracks, vertexingOptionsUT, state);
0486 
0487     if (!res.ok()) {
0488       BOOST_FAIL(res.error().message());
0489     }
0490 
0491     
0492     auto vertexCollectionUT = *res;
0493 
0494     
0495     CHECK_CLOSE_ABS(vertexCollectionUT.size(), nVertices, 2);
0496 
0497     if (debug) {
0498       std::cout << "########## RESULT: ########## Event " << iEvent
0499                 << std::endl;
0500       std::cout << "Number of true vertices: " << nVertices << std::endl;
0501       std::cout << "Number of reco vertices: " << vertexCollectionUT.size()
0502                 << std::endl;
0503 
0504       int count = 1;
0505       std::cout << "----- True vertices -----" << std::endl;
0506       for (const auto& vertex : trueVertices) {
0507         Vector3 pos = vertex.position();
0508         std::cout << count << ". True Vertex:\t Position:"
0509                   << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
0510                   << std::endl;
0511         std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
0512                   << std::endl;
0513         count++;
0514       }
0515       std::cout << "----- Reco vertices -----" << std::endl;
0516       count = 1;
0517       for (const auto& vertex : vertexCollectionUT) {
0518         Vector3 pos = vertex.position();
0519         std::cout << count << ". Reco Vertex:\t Position:"
0520                   << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
0521                   << std::endl;
0522         std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
0523                   << std::endl;
0524         count++;
0525       }
0526     }
0527 
0528     
0529     bool allVerticesFound = true;
0530     for (const auto& trueVertex : trueVertices) {
0531       Vector4 truePos = trueVertex.fullPosition();
0532       bool currentVertexFound = false;
0533       for (const auto& recoVertex : vertexCollectionUT) {
0534         Vector4 recoPos = recoVertex.fullPosition();
0535         
0536         double zDistance = std::abs(truePos[eZ] - recoPos[eZ]);
0537         if (zDistance < 2_mm) {
0538           currentVertexFound = true;
0539         }
0540       }
0541       if (!currentVertexFound) {
0542         allVerticesFound = false;
0543       }
0544     }
0545 
0546     
0547     BOOST_CHECK(allVerticesFound);
0548   }
0549 }
0550 
0551 
0552 
0553 
0554 BOOST_AUTO_TEST_CASE(iterative_finder_test_athena_reference) {
0555   
0556   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 2_T});
0557 
0558   
0559   EigenStepper<> stepper(bField);
0560 
0561   
0562   auto propagator = std::make_shared<Propagator>(stepper);
0563 
0564   
0565   Linearizer::Config ltConfig;
0566   ltConfig.bField = bField;
0567   ltConfig.propagator = propagator;
0568   Linearizer linearizer(ltConfig);
0569 
0570   using BilloirFitter = FullBilloirVertexFitter;
0571 
0572   
0573   BilloirFitter::Config vertexFitterCfg;
0574   vertexFitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0575   vertexFitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(
0576       &linearizer);
0577 
0578   BilloirFitter bFitter(vertexFitterCfg);
0579 
0580   ImpactPointEstimator::Config ipEstimatorCfg(bField, propagator);
0581   ImpactPointEstimator ipEstimator(ipEstimatorCfg);
0582 
0583   using ZScanSeedFinder = ZScanVertexFinder;
0584 
0585   ZScanSeedFinder::Config seedFinderCfg(ipEstimator);
0586   seedFinderCfg.extractParameters.connect<&InputTrack::extractParameters>();
0587 
0588   auto sFinder = std::make_shared<ZScanSeedFinder>(seedFinderCfg);
0589 
0590   IterativeVertexFinder::Config cfg(std::move(bFitter), std::move(sFinder),
0591                                     ipEstimator);
0592   cfg.maxVertices = 200;
0593   cfg.maximumChi2cutForSeeding = 49;
0594   cfg.significanceCutSeeding = 12;
0595   cfg.extractParameters.connect<&InputTrack::extractParameters>();
0596   cfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0597   cfg.field = bField;
0598 
0599   IterativeVertexFinder finder(std::move(cfg));
0600   IVertexFinder::State state{
0601       IterativeVertexFinder::State(*bField, magFieldContext)};
0602 
0603   auto csvData = readTracksAndVertexCSV(toolString);
0604   auto tracks = std::get<TracksData>(csvData);
0605 
0606   std::vector<InputTrack> inputTracks;
0607   for (const auto& trk : tracks) {
0608     inputTracks.emplace_back(&trk);
0609   }
0610 
0611   Vertex beamSpot = std::get<BeamSpotData>(csvData);
0612   
0613   SquareMatrix4 fullCovariance = SquareMatrix4::Zero();
0614   fullCovariance.topLeftCorner<3, 3>() = beamSpot.covariance();
0615   fullCovariance(eTime, eTime) = 100_ns;
0616   beamSpot.setFullCovariance(fullCovariance);
0617   VertexingOptions vertexingOptions(geoContext, magFieldContext, beamSpot);
0618 
0619   
0620   auto findResult = finder.find(inputTracks, vertexingOptions, state);
0621 
0622   BOOST_CHECK(findResult.ok());
0623 
0624   if (!findResult.ok()) {
0625     std::cout << findResult.error().message() << std::endl;
0626   }
0627 
0628   
0629   
0630 
0631   
0632   
0633   
0634   
0635 
0636   
0637   
0638   
0639   
0640   
0641   
0642   
0643   
0644   
0645   
0646   
0647   
0648 }
0649 
0650 BOOST_AUTO_TEST_SUITE_END()
0651 
0652 }