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