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