Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/tools/output_test_stream.hpp>
0011 #include <boost/test/unit_test.hpp>
0012
0013 #include "Acts/Definitions/Algebra.hpp"
0014 #include "Acts/Definitions/Direction.hpp"
0015 #include "Acts/Definitions/TrackParametrization.hpp"
0016 #include "Acts/Definitions/Units.hpp"
0017 #include "Acts/EventData/Charge.hpp"
0018 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0019 #include "Acts/EventData/TrackParameters.hpp"
0020 #include "Acts/Geometry/GeometryContext.hpp"
0021 #include "Acts/Geometry/GeometryIdentifier.hpp"
0022 #include "Acts/MagneticField/ConstantBField.hpp"
0023 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0024 #include "Acts/Propagator/EigenStepper.hpp"
0025 #include "Acts/Propagator/Propagator.hpp"
0026 #include "Acts/Surfaces/PerigeeSurface.hpp"
0027 #include "Acts/Surfaces/Surface.hpp"
0028 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0029 #include "Acts/Utilities/AnnealingUtility.hpp"
0030 #include "Acts/Utilities/Logger.hpp"
0031 #include "Acts/Utilities/Result.hpp"
0032 #include "Acts/Vertexing/AMVFInfo.hpp"
0033 #include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp"
0034 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
0035 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0036 #include "Acts/Vertexing/TrackAtVertex.hpp"
0037 #include "Acts/Vertexing/Vertex.hpp"
0038 #include "Acts/Vertexing/VertexingOptions.hpp"
0039
0040 #include <algorithm>
0041 #include <array>
0042 #include <cmath>
0043 #include <iostream>
0044 #include <map>
0045 #include <memory>
0046 #include <random>
0047 #include <tuple>
0048 #include <type_traits>
0049 #include <utility>
0050 #include <vector>
0051
0052 namespace Acts::Test {
0053
0054 using namespace Acts::UnitLiterals;
0055 using Acts::VectorHelpers::makeVector4;
0056
0057 // Set up logger
0058 ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("AMVFitterTests", Acts::Logging::INFO))
0059
0060 using Covariance = BoundSquareMatrix;
0061 using Propagator = Acts::Propagator<EigenStepper<>>;
0062 using Linearizer = HelicalTrackLinearizer;
0063
0064 // Create a test context
0065 GeometryContext geoContext = GeometryContext();
0066 MagneticFieldContext magFieldContext = MagneticFieldContext();
0067
0068 // Vertex x/y position distribution
0069 std::uniform_real_distribution<double> vXYDist(-0.1_mm, 0.1_mm);
0070 // Vertex z position distribution
0071 std::uniform_real_distribution<double> vZDist(-20_mm, 20_mm);
0072 // Track d0 distribution
0073 std::uniform_real_distribution<double> d0Dist(-0.01_mm, 0.01_mm);
0074 // Track z0 distribution
0075 std::uniform_real_distribution<double> z0Dist(-0.2_mm, 0.2_mm);
0076 // Track pT distribution
0077 std::uniform_real_distribution<double> pTDist(1._GeV, 30._GeV);
0078 // Track phi distribution
0079 std::uniform_real_distribution<double> phiDist(-M_PI, M_PI);
0080 // Track theta distribution
0081 std::uniform_real_distribution<double> thetaDist(1.0, M_PI - 1.0);
0082 // Track charge helper distribution
0083 std::uniform_real_distribution<double> qDist(-1, 1);
0084 // Distribution of track time (relative to vertex time). Values are unrealistic
0085 // and only used for testing purposes.
0086 std::uniform_real_distribution<double> relTDist(-4_ps, 4_ps);
0087 // Track IP resolution distribution
0088 std::uniform_real_distribution<double> resIPDist(0., 100._um);
0089 // Track angular distribution
0090 std::uniform_real_distribution<double> resAngDist(0., 0.1);
0091 // Track q/p resolution distribution
0092 std::uniform_real_distribution<double> resQoPDist(-0.1, 0.1);
0093 // Track time resolution distribution. Values are unrealistic and only used for
0094 // testing purposes.
0095 std::uniform_real_distribution<double> resTDist(0_ps, 8_ps);
0096
0097 /// @brief Unit test for AdaptiveMultiVertexFitter
0098 ///
0099 BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test) {
0100   // Set up RNG
0101   int mySeed = 31415;
0102   std::mt19937 gen(mySeed);
0103
0104   // Set up constant B-Field
0105   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0106
0107   // Set up EigenStepper
0108   EigenStepper<> stepper(bField);
0109
0110   // Set up propagator with void navigator
0111   auto propagator = std::make_shared<Propagator>(stepper);
0112
0113   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0114
0115   // IP 3D Estimator
0116   ImpactPointEstimator::Config ip3dEstCfg(bField, propagator);
0117   ImpactPointEstimator ip3dEst(ip3dEstCfg);
0118
0119   // Linearizer for BoundTrackParameters type test
0120   Linearizer::Config ltConfig;
0121   ltConfig.bField = bField;
0122   ltConfig.propagator = propagator;
0123   Linearizer linearizer(ltConfig);
0124
0125   AdaptiveMultiVertexFitter::Config fitterCfg(ip3dEst);
0126   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0127
0128   // Test smoothing
0129   fitterCfg.doSmoothing = true;
0130   fitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0131
0132   AdaptiveMultiVertexFitter fitter(std::move(fitterCfg));
0133
0134   // Create positions of three vertices, two of which (1 and 2) are
0135   // close to one another and will share a common track later
0136   Vector3 vtxPos1(-0.15_mm, -0.1_mm, -1.5_mm);
0137   Vector3 vtxPos2(-0.1_mm, -0.15_mm, -3._mm);
0138   Vector3 vtxPos3(0.2_mm, 0.2_mm, 10._mm);
0139
0140   std::vector<Vector3> vtxPosVec{vtxPos1, vtxPos2, vtxPos3};
0141
0142   // Resolutions, use the same for all tracks
0143   double resD0 = resIPDist(gen);
0144   double resZ0 = resIPDist(gen);
0145   double resPh = resAngDist(gen);
0146   double resTh = resAngDist(gen);
0147   double resQp = resQoPDist(gen);
0148
0149   std::vector<Vertex> vtxList;
0150   for (auto& vtxPos : vtxPosVec) {
0151     Vertex vtx(vtxPos);
0152     // Set some vertex covariance
0153     SquareMatrix4 posCovariance(SquareMatrix4::Identity());
0154     vtx.setFullCovariance(posCovariance);
0155     // Add to vertex list
0156     vtxList.push_back(vtx);
0157   }
0158
0159   std::vector<Vertex*> vtxPtrList;
0160   ACTS_DEBUG("All vertices in test case:");
0161   int cv = 0;
0162   for (auto& vtx : vtxList) {
0163     cv++;
0164     ACTS_DEBUG("\t" << cv << ". vertex ptr: " << &vtx);
0165     vtxPtrList.push_back(&vtx);
0166   }
0167
0168   std::vector<BoundTrackParameters> allTracks;
0169
0170   unsigned int nTracksPerVtx = 4;
0171   // Construct nTracksPerVtx * 3 (3 vertices) random track emerging
0172   // from vicinity of vertex positions
0173   for (unsigned int iTrack = 0; iTrack < nTracksPerVtx * vtxPosVec.size();
0174        iTrack++) {
0175     // Construct positive or negative charge randomly
0176     double q = qDist(gen) < 0 ? -1. : 1.;
0177
0178     // Fill vector of track objects with simple covariance matrix
0179     Covariance covMat;
0180
0181     covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0.,
0182         0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0.,
0183         0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;
0184
0185     // Index of current vertex
0186     int vtxIdx = (int)(iTrack / nTracksPerVtx);
0187
0188     // Construct random track parameters
0189     BoundTrackParameters::ParametersVector paramVec;
0190     paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen),
0191         q / pTDist(gen), 0.;
0192
0193     std::shared_ptr<PerigeeSurface> perigeeSurface =
0194         Surface::makeShared<PerigeeSurface>(vtxPosVec[vtxIdx]);
0195
0196     allTracks.emplace_back(perigeeSurface, paramVec, std::move(covMat),
0197                            ParticleHypothesis::pion());
0198   }
0199
0200   int ct = 0;
0201   ACTS_DEBUG("All tracks in test case:");
0202   for (auto& trk : allTracks) {
0203     ct++;
0204     ACTS_DEBUG("\t" << ct << ". track ptr: " << &trk);
0205   }
0206
0207   AdaptiveMultiVertexFitter::State state(*bField, magFieldContext);
0208
0209   for (unsigned int iTrack = 0; iTrack < nTracksPerVtx * vtxPosVec.size();
0210        iTrack++) {
0211     // Index of current vertex
0212     int vtxIdx = (int)(iTrack / nTracksPerVtx);
0213
0214     InputTrack inputTrack{&allTracks[iTrack]};
0215
0216     state.vtxInfoMap[&(vtxList[vtxIdx])].trackLinks.push_back(inputTrack);
0217     state.tracksAtVerticesMap.insert(
0218         std::make_pair(std::make_pair(inputTrack, &(vtxList[vtxIdx])),
0219                        TrackAtVertex(1., allTracks[iTrack], inputTrack)));
0220
0221     // Use first track also for second vertex to let vtx1 and vtx2
0222     // share this track
0223     if (iTrack == 0) {
0224       state.vtxInfoMap[&(vtxList.at(1))].trackLinks.push_back(inputTrack);
0225       state.tracksAtVerticesMap.insert(
0226           std::make_pair(std::make_pair(inputTrack, &(vtxList.at(1))),
0227                          TrackAtVertex(1., allTracks[iTrack], inputTrack)));
0228     }
0229   }
0230
0231   for (auto& vtx : vtxPtrList) {
0232     state.addVertexToMultiMap(*vtx);
0233     ACTS_DEBUG("Vertex, with ptr: " << vtx);
0234     for (auto& trk : state.vtxInfoMap[vtx].trackLinks) {
0235       ACTS_DEBUG("\t track ptr: " << trk);
0236     }
0237   }
0238
0239   ACTS_DEBUG("Checking all vertices linked to a single track:");
0240   for (auto& trk : allTracks) {
0241     ACTS_DEBUG("Track with ptr: " << &trk);
0242     auto range = state.trackToVerticesMultiMap.equal_range(InputTrack{&trk});
0243     for (auto vtxIter = range.first; vtxIter != range.second; ++vtxIter) {
0244       ACTS_DEBUG("\t used by vertex: " << vtxIter->second);
0245     }
0246   }
0247
0248   // Copy vertex seeds from state.vertexCollection to new
0249   // list in order to be able to compare later
0250   std::vector<Vertex> seedListCopy = vtxList;
0251
0252   auto res1 = fitter.addVtxToFit(state, vtxList.at(0), vertexingOptions);
0253   ACTS_DEBUG("Tracks linked to each vertex AFTER fit:");
0254   int c = 0;
0255   for (auto& vtx : vtxPtrList) {
0256     c++;
0257     ACTS_DEBUG(c << ". vertex, with ptr: " << vtx);
0258     for (const auto& trk : state.vtxInfoMap[vtx].trackLinks) {
0259       ACTS_DEBUG("\t track ptr: " << trk);
0260     }
0261   }
0262
0263   ACTS_DEBUG("Checking all vertices linked to a single track AFTER fit:");
0264   for (auto& trk : allTracks) {
0265     ACTS_DEBUG("Track with ptr: " << &trk);
0266     auto range = state.trackToVerticesMultiMap.equal_range(InputTrack{&trk});
0267     for (auto vtxIter = range.first; vtxIter != range.second; ++vtxIter) {
0268       ACTS_DEBUG("\t used by vertex: " << vtxIter->second);
0269     }
0270   }
0271
0272   BOOST_CHECK(res1.ok());
0273
0274   ACTS_DEBUG("Vertex positions after fit of vertex 1 and 2:");
0275   for (std::size_t vtxIter = 0; vtxIter < 3; vtxIter++) {
0276     ACTS_DEBUG("Vtx " << vtxIter + 1 << ", seed position:\n "
0277                       << seedListCopy.at(vtxIter).fullPosition()
0278                       << "\nFitted position:\n "
0279                       << vtxList.at(vtxIter).fullPosition());
0280   }
0281
0282   // After fit of first vertex, only first and second vertex seed
0283   // should have been modified while third vertex should remain untouched
0284   BOOST_CHECK_NE(vtxList.at(0).fullPosition(),
0285                  seedListCopy.at(0).fullPosition());
0286   BOOST_CHECK_NE(vtxList.at(1).fullPosition(),
0287                  seedListCopy.at(1).fullPosition());
0288   BOOST_CHECK_EQUAL(vtxList.at(2).fullPosition(),
0289                     seedListCopy.at(2).fullPosition());
0290
0291   CHECK_CLOSE_ABS(vtxList.at(0).fullPosition(),
0292                   seedListCopy.at(0).fullPosition(), 1_mm);
0293   CHECK_CLOSE_ABS(vtxList.at(1).fullPosition(),
0294                   seedListCopy.at(1).fullPosition(), 1_mm);
0295
0296   auto res2 = fitter.addVtxToFit(state, vtxList.at(2), vertexingOptions);
0297   BOOST_CHECK(res2.ok());
0298
0299   // Now also the third vertex should have been modified and fitted
0300   BOOST_CHECK_NE(vtxList.at(2).fullPosition(),
0301                  seedListCopy.at(2).fullPosition());
0302   CHECK_CLOSE_ABS(vtxList.at(2).fullPosition(),
0303                   seedListCopy.at(2).fullPosition(), 1_mm);
0304
0305   ACTS_DEBUG("Vertex positions after fit of vertex 3:");
0306   ACTS_DEBUG("Vtx 1, seed position:\n " << seedListCopy.at(0).fullPosition()
0307                                         << "\nFitted position:\n "
0308                                         << vtxList.at(0).fullPosition());
0309   ACTS_DEBUG("Vtx 2, seed position:\n " << seedListCopy.at(1).fullPosition()
0310                                         << "\nFitted position:\n "
0311                                         << vtxList.at(1).fullPosition());
0312   ACTS_DEBUG("Vtx 3, seed position:\n " << seedListCopy.at(2).fullPosition()
0313                                         << "\nFitted position:\n "
0314                                         << vtxList.at(2).fullPosition());
0315 }
0316
0317 /// @brief Unit test for fitting a 4D vertex position
0318 ///
0319 BOOST_AUTO_TEST_CASE(time_fitting) {
0320   // Set up RNG
0321   int mySeed = 31415;
0322   std::mt19937 gen(mySeed);
0323
0324   // Set up constant B-Field
0325   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 1_T});
0326
0327   // Set up EigenStepper
0328   EigenStepper<> stepper(bField);
0329
0330   // Set up propagator with void navigator
0331   auto propagator = std::make_shared<Propagator>(stepper);
0332
0333   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0334
0335   ImpactPointEstimator::Config ip3dEstCfg(bField, propagator);
0336   ImpactPointEstimator ip3dEst(ip3dEstCfg);
0337
0338   AdaptiveMultiVertexFitter::Config fitterCfg(ip3dEst);
0339
0340   // Linearizer for BoundTrackParameters type test
0341   Linearizer::Config ltConfig;
0342   ltConfig.bField = bField;
0343   ltConfig.propagator = propagator;
0344   Linearizer linearizer(ltConfig);
0345
0346   // Test smoothing
0347   fitterCfg.doSmoothing = true;
0348   // Do time fit
0349   fitterCfg.useTime = true;
0350   fitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0351   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0352
0353   AdaptiveMultiVertexFitter fitter(std::move(fitterCfg));
0354
0355   // Vertex position
0356   double trueVtxTime = 40.0_ps;
0357   Vector3 trueVtxPos(-0.15_mm, -0.1_mm, -1.5_mm);
0358
0359   // Seed position of the vertex
0360   Vector4 vtxSeedPos(0.0_mm, 0.0_mm, -1.4_mm, 0.0_ps);
0361
0362   Vertex vtx(vtxSeedPos);
0363   // Set initial covariance matrix to a large value
0364   SquareMatrix4 initialCovariance(SquareMatrix4::Identity() * 1e+8);
0365   vtx.setFullCovariance(initialCovariance);
0366
0367   // Vector of all tracks
0368   std::vector<BoundTrackParameters> trks;
0369
0370   unsigned int nTracks = 4;
0371   for (unsigned int _ = 0; _ < nTracks; _++) {
0372     // Construct positive or negative charge randomly
0373     double q = qDist(gen) < 0 ? -1. : 1.;
0374
0375     // Track resolution
0376     double resD0 = resIPDist(gen);
0377     double resZ0 = resIPDist(gen);
0378     double resPh = resAngDist(gen);
0379     double resTh = resAngDist(gen);
0380     double resQp = resQoPDist(gen);
0381     double resT = resTDist(gen);
0382
0383     // Random diagonal covariance matrix
0384     Covariance covMat;
0385
0386     // clang-format off
0387     covMat <<
0388       resD0 * resD0, 0., 0., 0., 0., 0.,
0389       0., resZ0 * resZ0, 0., 0., 0., 0.,
0390       0., 0., resPh * resPh, 0., 0., 0.,
0391       0., 0., 0., resTh * resTh, 0., 0.,
0392       0., 0., 0., 0., resQp * resQp, 0.,
0393       0., 0., 0., 0., 0., resT * resT;
0394     // clang-format on
0395
0396     // Random track parameters
0397     BoundTrackParameters::ParametersVector paramVec;
0398     paramVec << d0Dist(gen), z0Dist(gen), phiDist(gen), thetaDist(gen),
0399         q / pTDist(gen), trueVtxTime + relTDist(gen);
0400
0401     std::shared_ptr<PerigeeSurface> perigeeSurface =
0402         Surface::makeShared<PerigeeSurface>(trueVtxPos);
0403
0404     trks.emplace_back(perigeeSurface, paramVec, std::move(covMat),
0405                       ParticleHypothesis::pion());
0406   }
0407
0408   std::vector<const BoundTrackParameters*> trksPtr;
0409   for (const auto& trk : trks) {
0410     trksPtr.push_back(&trk);
0411   }
0412
0413   // Prepare fitter state
0414   AdaptiveMultiVertexFitter::State state(*bField, magFieldContext);
0415
0416   for (const auto& trk : trks) {
0417     ACTS_DEBUG("Track parameters:\n" << trk);
0418     // Index of current vertex
0419     state.vtxInfoMap[&vtx].trackLinks.push_back(InputTrack{&trk});
0420     state.tracksAtVerticesMap.insert(
0421         std::make_pair(std::make_pair(InputTrack{&trk}, &vtx),
0422                        TrackAtVertex(1., trk, InputTrack{&trk})));
0423   }
0424
0425   state.addVertexToMultiMap(vtx);
0426
0427   auto res = fitter.addVtxToFit(state, vtx, vertexingOptions);
0428
0429   BOOST_CHECK(res.ok());
0430
0431   ACTS_DEBUG("Truth vertex position:  " << trueVtxPos.transpose());
0432   ACTS_DEBUG("Fitted vertex position: " << vtx.position().transpose());
0433
0434   ACTS_DEBUG("Truth vertex time:  " << trueVtxTime);
0435   ACTS_DEBUG("Fitted vertex time: " << vtx.time());
0436
0437   // Check that true vertex position and time approximately recovered
0438   CHECK_CLOSE_ABS(trueVtxPos, vtx.position(), 60_um);
0439   CHECK_CLOSE_ABS(trueVtxTime, vtx.time(), 2_ps);
0440
0441   const SquareMatrix4& vtxCov = vtx.fullCovariance();
0442
0443   ACTS_DEBUG("Vertex 4D covariance after the fit:\n" << vtxCov);
0444
0445   // Check that variances of the vertex position/time are positive
0446   for (std::size_t i = 0; i <= 3; i++) {
0447     BOOST_CHECK_GT(vtxCov(i, i), 0.);
0448   }
0449
0450   // Check that the covariance matrix is approximately symmetric
0451   CHECK_CLOSE_ABS(vtxCov - vtxCov.transpose(), SquareMatrix4::Zero(), 1e-5);
0452 }
0453
0454 /// @brief Unit test for AdaptiveMultiVertexFitter
0455 /// based on Athena unit test, i.e. same setting and
0456 /// test values are used here
0457 BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) {
0458   // Set up constant B-Field
0459   auto bField = std::make_shared<ConstantBField>(Vector3{0.0, 0.0, 2_T});
0460
0461   // Set up EigenStepper
0462   // EigenStepper<> stepper(bField);
0463   EigenStepper<> stepper(bField);
0464
0465   // Set up propagator with void navigator
0466   auto propagator = std::make_shared<Propagator>(stepper);
0467
0468   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0469
0470   ImpactPointEstimator::Config ip3dEstCfg(bField, propagator);
0471   ImpactPointEstimator ip3dEst(ip3dEstCfg);
0472
0473   std::vector<double> temperatures(1, 3.);
0474   AnnealingUtility::Config annealingConfig;
0475   annealingConfig.setOfTemperatures = temperatures;
0476   AnnealingUtility annealingUtility(annealingConfig);
0477
0478   AdaptiveMultiVertexFitter::Config fitterCfg(ip3dEst);
0479
0480   fitterCfg.annealingTool = annealingUtility;
0481   fitterCfg.extractParameters.connect<&InputTrack::extractParameters>();
0482
0483   // Linearizer for BoundTrackParameters type test
0484   Linearizer::Config ltConfig;
0485   ltConfig.bField = bField;
0486   ltConfig.propagator = propagator;
0487   Linearizer linearizer(ltConfig);
0488
0489   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
0490
0491   // Test smoothing
0492   // fitterCfg.doSmoothing = true;
0493
0494   AdaptiveMultiVertexFitter fitter(std::move(fitterCfg));
0495
0496   // Create first vector of tracks
0497   Vector3 pos1a(0.5_mm, -0.5_mm, 2.4_mm);
0498   Vector3 mom1a(1000_MeV, 0_MeV, -500_MeV);
0499   Vector3 pos1b(0.5_mm, -0.5_mm, 3.5_mm);
0500   Vector3 mom1b(0_MeV, 1000_MeV, 500_MeV);
0501   Vector3 pos1c(-0.2_mm, 0.1_mm, 3.4_mm);
0502   Vector3 mom1c(-50_MeV, 180_MeV, 300_MeV);
0503
0504   Vector3 pos1d(-0.1_mm, 0.3_mm, 3.0_mm);
0505   Vector3 mom1d(-80_MeV, 480_MeV, -100_MeV);
0506   Vector3 pos1e(-0.01_mm, 0.01_mm, 2.9_mm);
0507   Vector3 mom1e(-600_MeV, 10_MeV, 210_MeV);
0508
0509   Vector3 pos1f(-0.07_mm, 0.03_mm, 2.5_mm);
0510   Vector3 mom1f(240_MeV, 110_MeV, 150_MeV);
0511
0512   // Start creating some track parameters
0513   Covariance covMat1;
0514   covMat1 << 1_mm * 1_mm, 0, 0., 0, 0., 0, 0, 1_mm * 1_mm, 0, 0., 0, 0, 0., 0,
0515       0.1, 0, 0, 0, 0, 0., 0, 0.1, 0, 0, 0., 0, 0, 0, 1. / (10_GeV * 10_GeV), 0,
0516       0, 0, 0, 0, 0, 1_ns;
0517
0518   std::vector<BoundTrackParameters> params1 = {
0519       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos1a),
0520                                    geoContext, makeVector4(pos1a, 0),
0521                                    mom1a.normalized(), 1_e / mom1a.norm(),
0522                                    covMat1, ParticleHypothesis::pion())
0523           .value(),
0524       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos1b),
0525                                    geoContext, makeVector4(pos1b, 0),
0526                                    mom1b.normalized(), -1_e / mom1b.norm(),
0527                                    covMat1, ParticleHypothesis::pion())
0528           .value(),
0529       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos1c),
0530                                    geoContext, makeVector4(pos1c, 0),
0531                                    mom1c.normalized(), 1_e / mom1c.norm(),
0532                                    covMat1, ParticleHypothesis::pion())
0533           .value(),
0534       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos1d),
0535                                    geoContext, makeVector4(pos1d, 0),
0536                                    mom1d.normalized(), -1_e / mom1d.norm(),
0537                                    covMat1, ParticleHypothesis::pion())
0538           .value(),
0539       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos1e),
0540                                    geoContext, makeVector4(pos1e, 0),
0541                                    mom1e.normalized(), 1_e / mom1e.norm(),
0542                                    covMat1, ParticleHypothesis::pion())
0543           .value(),
0544       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos1f),
0545                                    geoContext, makeVector4(pos1f, 0),
0546                                    mom1f.normalized(), -1_e / mom1f.norm(),
0547                                    covMat1, ParticleHypothesis::pion())
0548           .value(),
0549   };
0550
0551   // Create second vector of tracks
0552   Vector3 pos2a(0.2_mm, 0_mm, -4.9_mm);
0553   Vector3 mom2a(5000_MeV, 30_MeV, 200_MeV);
0554   Vector3 pos2b(-0.5_mm, 0.1_mm, -5.1_mm);
0555   Vector3 mom2b(800_MeV, 1200_MeV, 200_MeV);
0556   Vector3 pos2c(0.05_mm, -0.5_mm, -4.7_mm);
0557   Vector3 mom2c(400_MeV, -300_MeV, -200_MeV);
0558
0559   // Define covariance as used in athena unit test
0560   Covariance covMat2 = covMat1;
0561
0562   std::vector<BoundTrackParameters> params2 = {
0563       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos2a),
0564                                    geoContext, makeVector4(pos2a, 0),
0565                                    mom2a.normalized(), 1_e / mom2a.norm(),
0566                                    covMat2, ParticleHypothesis::pion())
0567           .value(),
0568       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos2b),
0569                                    geoContext, makeVector4(pos2b, 0),
0570                                    mom2b.normalized(), -1_e / mom2b.norm(),
0571                                    covMat2, ParticleHypothesis::pion())
0572           .value(),
0573       BoundTrackParameters::create(Surface::makeShared<PerigeeSurface>(pos2c),
0574                                    geoContext, makeVector4(pos2c, 0),
0575                                    mom2c.normalized(), -1_e / mom2c.norm(),
0576                                    covMat2, ParticleHypothesis::pion())
0577           .value(),
0578   };
0579
0580   AdaptiveMultiVertexFitter::State state(*bField, magFieldContext);
0581
0582   // The constraint vertex position covariance
0583   SquareMatrix4 covConstr(SquareMatrix4::Identity());
0584   covConstr = covConstr * 1e+8;
0585   covConstr(3, 3) = 0.;
0586
0587   // Prepare first vertex
0588   Vector3 vtxPos1(0.15_mm, 0.15_mm, 2.9_mm);
0589   Vertex vtx1(vtxPos1);
0590
0591   // Add to vertex list
0592   state.vertexCollection.push_back(&vtx1);
0593
0594   // The constraint vtx for vtx1
0595   Vertex vtx1Constr(vtxPos1);
0596   vtx1Constr.setFullCovariance(covConstr);
0597   vtx1Constr.setFitQuality(0, -3);
0598
0599   // Prepare vtx info for fitter
0600   VertexInfo vtxInfo1;
0601   vtxInfo1.seedPosition = vtxInfo1.linPoint;
0602   vtxInfo1.linPoint.setZero();
0603   vtxInfo1.linPoint.head<3>() = vtxPos1;
0604   vtxInfo1.constraint = std::move(vtx1Constr);
0605   vtxInfo1.oldPosition = vtxInfo1.linPoint;
0606
0607   for (const auto& trk : params1) {
0608     vtxInfo1.trackLinks.push_back(InputTrack{&trk});
0609     state.tracksAtVerticesMap.insert(
0610         std::make_pair(std::make_pair(InputTrack{&trk}, &vtx1),
0611                        TrackAtVertex(1.5, trk, InputTrack{&trk})));
0612   }
0613
0614   // Prepare second vertex
0615   Vector3 vtxPos2(0.3_mm, -0.2_mm, -4.8_mm);
0616   Vertex vtx2(vtxPos2);
0617
0618   // Add to vertex list
0619   state.vertexCollection.push_back(&vtx2);
0620
0621   // The constraint vtx for vtx2
0622   Vertex vtx2Constr(vtxPos2);
0623   vtx2Constr.setFullCovariance(covConstr);
0624   vtx2Constr.setFitQuality(0, -3);
0625
0626   // Prepare vtx info for fitter
0627   VertexInfo vtxInfo2;
0628   vtxInfo2.linPoint.setZero();
0629   vtxInfo2.linPoint.head<3>() = vtxPos2;
0630   vtxInfo2.constraint = std::move(vtx2Constr);
0631   vtxInfo2.oldPosition = vtxInfo2.linPoint;
0632   vtxInfo2.seedPosition = vtxInfo2.linPoint;
0633
0634   for (const auto& trk : params2) {
0635     vtxInfo2.trackLinks.push_back(InputTrack{&trk});
0636     state.tracksAtVerticesMap.insert(
0637         std::make_pair(std::make_pair(InputTrack{&trk}, &vtx2),
0638                        TrackAtVertex(1.5, trk, InputTrack{&trk})));
0639   }
0640
0641   state.vtxInfoMap[&vtx1] = std::move(vtxInfo1);
0642   state.vtxInfoMap[&vtx2] = std::move(vtxInfo2);
0643
0644   state.addVertexToMultiMap(vtx1);
0645   state.addVertexToMultiMap(vtx2);
0646
0647   // Fit vertices
0648   fitter.fit(state, vertexingOptions);
0649
0650   auto vtx1Fitted = state.vertexCollection.at(0);
0651   auto vtx1PosFitted = vtx1Fitted->position();
0652   auto vtx1CovFitted = vtx1Fitted->covariance();
0653   auto trks1 = state.vtxInfoMap.at(vtx1Fitted).trackLinks;
0654   auto vtx1FQ = vtx1Fitted->fitQuality();
0655
0656   auto vtx2Fitted = state.vertexCollection.at(1);
0657   auto vtx2PosFitted = vtx2Fitted->position();
0658   auto vtx2CovFitted = vtx2Fitted->covariance();
0659   auto trks2 = state.vtxInfoMap.at(vtx2Fitted).trackLinks;
0660   auto vtx2FQ = vtx2Fitted->fitQuality();
0661
0662   // Vertex 1
0663   ACTS_DEBUG("Vertex 1, position: " << vtx1PosFitted);
0664   ACTS_DEBUG("Vertex 1, covariance: " << vtx1CovFitted);
0665   for (const auto& trk : trks1) {
0666     auto& trkAtVtx =
0667         state.tracksAtVerticesMap.at(std::make_pair(trk, vtx1Fitted));
0668     ACTS_DEBUG("\tTrack weight:" << trkAtVtx.trackWeight);
0669   }
0670   ACTS_DEBUG("Vertex 1, chi2: " << vtx1FQ.first);
0671   ACTS_DEBUG("Vertex 1, ndf: " << vtx1FQ.second);
0672
0673   // Vertex 2
0674   ACTS_DEBUG("Vertex 2, position: " << vtx2PosFitted);
0675   ACTS_DEBUG("Vertex 2, covariance: " << vtx2CovFitted);
0676   for (const auto& trk : trks2) {
0677     auto& trkAtVtx =
0678         state.tracksAtVerticesMap.at(std::make_pair(trk, vtx2Fitted));
0679     ACTS_DEBUG("\tTrack weight:" << trkAtVtx.trackWeight);
0680   }
0681   ACTS_DEBUG("Vertex 2, chi2: " << vtx2FQ.first);
0682   ACTS_DEBUG("Vertex 2, ndf: " << vtx2FQ.second);
0683
0684   // Expected values from Athena implementation
0685   // Vertex 1
0686   const Vector3 expVtx1Pos(0.077_mm, -0.189_mm, 2.924_mm);
0687
0688   // Helper matrix to create const expVtx1Cov below
0689   SquareMatrix3 expVtx1Cov;
0690   expVtx1Cov << 0.329, 0.016, -0.035, 0.016, 0.250, 0.085, -0.035, 0.085, 0.242;
0691
0692   ActsVector<6> expVtx1TrkWeights;
0693   expVtx1TrkWeights << 0.8128, 0.7994, 0.8164, 0.8165, 0.8165, 0.8119;
0694   const double expVtx1chi2 = 0.9812;
0695   const double expVtx1ndf = 6.7474;
0696
0697   // Vertex 2
0698   const Vector3 expVtx2Pos(-0.443_mm, -0.044_mm, -4.829_mm);
0699   // Helper matrix to create const expVtx2Cov below
0700   SquareMatrix3 expVtx2Cov;
0701   expVtx2Cov << 1.088, 0.028, -0.066, 0.028, 0.643, 0.073, -0.066, 0.073, 0.435;
0702
0703   const Vector3 expVtx2TrkWeights(0.8172, 0.8150, 0.8137);
0704   const double expVtx2chi2 = 0.2114;
0705   const double expVtx2ndf = 1.8920;
0706
0707   // Compare the results
0708   // Vertex 1
0709   CHECK_CLOSE_ABS(vtx1PosFitted, expVtx1Pos, 0.001_mm);
0710   CHECK_CLOSE_ABS(vtx1CovFitted, expVtx1Cov, 0.001_mm);
0711   int trkCount = 0;
0712   for (const auto& trk : trks1) {
0713     auto& trkAtVtx =
0714         state.tracksAtVerticesMap.at(std::make_pair(trk, vtx1Fitted));
0715     CHECK_CLOSE_ABS(trkAtVtx.trackWeight, expVtx1TrkWeights[trkCount], 0.001);
0716     trkCount++;
0717   }
0718   CHECK_CLOSE_ABS(vtx1FQ.first, expVtx1chi2, 0.001);
0719   CHECK_CLOSE_ABS(vtx1FQ.second, expVtx1ndf, 0.001);
0720
0721   // Vertex 2
0722   CHECK_CLOSE_ABS(vtx2PosFitted, expVtx2Pos, 0.001_mm);
0723   CHECK_CLOSE_ABS(vtx2CovFitted, expVtx2Cov, 0.001_mm);
0724   trkCount = 0;
0725   for (const auto& trk : trks2) {
0726     auto& trkAtVtx =
0727         state.tracksAtVerticesMap.at(std::make_pair(trk, vtx2Fitted));
0728     CHECK_CLOSE_ABS(trkAtVtx.trackWeight, expVtx2TrkWeights[trkCount], 0.001);
0729     trkCount++;
0730   }
0731   CHECK_CLOSE_ABS(vtx2FQ.first, expVtx2chi2, 0.001);
0732   CHECK_CLOSE_ABS(vtx2FQ.second, expVtx2ndf, 0.001);
0733 }
0734
0735 }  // namespace Acts::Test