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