Core/Vertexing/GridDensityVertexFinderTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
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/TrackParameters.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0018 #include "Acts/Surfaces/PerigeeSurface.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Utilities/Intersection.hpp"
0021 #include "Acts/Utilities/Result.hpp"
0022 #include "Acts/Utilities/UnitVectors.hpp"
0023 #include "Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp"
0024 #include "Acts/Vertexing/AdaptiveGridTrackDensity.hpp"
0025 #include "Acts/Vertexing/GaussianGridTrackDensity.hpp"
0026 #include "Acts/Vertexing/GridDensityVertexFinder.hpp"
0027 #include "Acts/Vertexing/IVertexFinder.hpp"
0028 #include "Acts/Vertexing/Vertex.hpp"
0029 #include "Acts/Vertexing/VertexingOptions.hpp"
0030 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0031
0032 #include <iostream>
0033 #include <memory>
0034 #include <numbers>
0035 #include <random>
0036 #include <system_error>
0037 #include <vector>
0038
0039 using namespace Acts;
0040 using namespace Acts::UnitLiterals;
0041 using Acts::VectorHelpers::makeVector4;
0042
0043 namespace ActsTests {
0044
0045 using Covariance = BoundSquareMatrix;
0046
0047 // Create a test context
0048 GeometryContext geoContext = GeometryContext();
0049 MagneticFieldContext magFieldContext = MagneticFieldContext();
0050
0051 const double zVertexPos1 = 12.;
0052 const double zVertexPos2 = -3.;
0053 // x position
0054 std::normal_distribution<double> xdist(1_mm, 0.1_mm);
0055 // y position
0056 std::normal_distribution<double> ydist(-0.7_mm, 0.1_mm);
0057 // z1 position
0058 std::normal_distribution<double> z1dist(zVertexPos1 * 1_mm, 1_mm);
0059 // z2 position
0060 std::normal_distribution<double> z2dist(zVertexPos2 * 1_mm, 0.5_mm);
0061 // Track pT distribution
0062 std::uniform_real_distribution<double> pTDist(0.1_GeV, 100_GeV);
0063 // Track phi distribution
0064 std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0065                                                std::numbers::pi);
0066 // Track eta distribution
0067 std::uniform_real_distribution<double> etaDist(-4., 4.);
0068
0069 BOOST_AUTO_TEST_SUITE(VertexingSuite)
0070 ///
0071 /// @brief Unit test for GridDensityVertexFinder without caching
0072 /// of track density values
0073 ///
0074 BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_test) {
0075   bool debugMode = false;
0076
0077   // Note that the AdaptiveGridTrackDensity and the GaussianGridTrackDensity
0078   // only furnish exactly the same results for uneven mainGridSize, where the
0079   // binnings of the two track densities align. For even mainGridSize the
0080   // binning of GaussianGridTrackDensity is:
0081   // ..., [-binSize, 0), [0, binSize), ...
0082   // and the binning of AdaptiveGridTrackDensity is:
0083   // ..., [-0.5*binSize, 0.5*binSize), [0.5*binSize, 1.5*binSize), ...
0084   // This is because the AdaptiveGridTrackDensity always has 0 as a bin center.
0085   // As a consequence of these different binnings, results would be shifted for
0086   // binSize/2 if mainGridSize is even.
0087   const int mainGridSize = 3001;
0088   const int trkGridSize = 35;
0089
0090   Covariance covMat = Covariance::Identity();
0091
0092   // Perigee surface for track parameters
0093   Vector3 pos0{0, 0, 0};
0094   std::shared_ptr<PerigeeSurface> perigeeSurface =
0095       Surface::makeShared<PerigeeSurface>(pos0);
0096
0097   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0098
0099   using Finder1 = GridDensityVertexFinder;
0100   GaussianGridTrackDensity::Config gDensityConfig(100, mainGridSize,
0101                                                   trkGridSize);
0102   GaussianGridTrackDensity gDensity(gDensityConfig);
0103   Finder1::Config cfg1(gDensity);
0104   cfg1.cacheGridStateForTrackRemoval = false;
0105   cfg1.extractParameters.connect<&InputTrack::extractParameters>();
0106   Finder1 finder1(cfg1);
0107   IVertexFinder::State state1 = finder1.makeState(magFieldContext);
0108
0109   // Use custom grid density here with same bin size as Finder1
0110   AdaptiveGridTrackDensity::Config adaptiveDensityConfig;
0111   adaptiveDensityConfig.spatialTrkGridSizeRange = {trkGridSize, trkGridSize};
0112   adaptiveDensityConfig.spatialBinExtent = 2. / 30.01 * 1_mm;
0113   AdaptiveGridTrackDensity adaptiveDensity(adaptiveDensityConfig);
0114
0115   using Finder2 = AdaptiveGridDensityVertexFinder;
0116   Finder2::Config cfg2(adaptiveDensity);
0117   cfg2.cacheGridStateForTrackRemoval = false;
0118   cfg2.extractParameters.connect<&InputTrack::extractParameters>();
0119   Finder2 finder2(cfg2);
0120   IVertexFinder::State state2 = finder2.makeState(magFieldContext);
0121
0122   int mySeed = 31415;
0123   std::mt19937 gen(mySeed);
0124   unsigned int nTracks = 200;
0125
0126   std::vector<BoundTrackParameters> trackVec;
0127   trackVec.reserve(nTracks);
0128
0129   // Create nTracks tracks for test case
0130   for (unsigned int i = 0; i < nTracks; i++) {
0131     // The position of the particle
0132     Vector3 pos(xdist(gen), ydist(gen), 0);
0133
0134     // Create momentum and charge of track
0135     double pt = pTDist(gen);
0136     double phi = phiDist(gen);
0137     double eta = etaDist(gen);
0138     double charge = std::copysign(1., etaDist(gen));
0139
0140     // project the position on the surface
0141     Vector3 direction = makeDirectionFromPhiEta(phi, eta);
0142     Intersection3D intersection =
0143         perigeeSurface->intersect(geoContext, pos, direction).closest();
0144     pos = intersection.position();
0145
0146     // Produce most of the tracks at near z1 position,
0147     // some near z2. Highest track density then expected at z1
0148     pos[eZ] = ((i % 4) == 0) ? z2dist(gen) : z1dist(gen);
0149
0150     trackVec.push_back(BoundTrackParameters::create(
0151                            geoContext, perigeeSurface, makeVector4(pos, 0),
0152                            direction, charge / pt, covMat,
0153                            ParticleHypothesis::pion())
0154                            .value());
0155   }
0156
0157   std::vector<InputTrack> inputTracks;
0158   for (const auto& trk : trackVec) {
0159     inputTracks.emplace_back(&trk);
0160   }
0161
0162   auto res1 = finder1.find(inputTracks, vertexingOptions, state1);
0163   if (!res1.ok()) {
0164     std::cout << res1.error().message() << std::endl;
0165   }
0166
0167   auto res2 = finder2.find(inputTracks, vertexingOptions, state2);
0168   if (!res2.ok()) {
0169     std::cout << res2.error().message() << std::endl;
0170   }
0171
0172   double zResult1 = 0;
0173   if (res1.ok()) {
0174     BOOST_CHECK(!(*res1).empty());
0175     Vector3 result1 = (*res1).back().position();
0176     if (debugMode) {
0177       std::cout << "Vertex position result 1: " << result1.transpose()
0178                 << std::endl;
0179     }
0180     CHECK_CLOSE_ABS(result1[eZ], zVertexPos1, 1_mm);
0181     zResult1 = result1[eZ];
0182   }
0183
0184   double zResult2 = 0;
0185   if (res2.ok()) {
0186     BOOST_CHECK(!(*res2).empty());
0187     Vector3 result2 = (*res2).back().position();
0188     if (debugMode) {
0189       std::cout << "Vertex position result 2: " << result2.transpose()
0190                 << std::endl;
0191     }
0192     CHECK_CLOSE_ABS(result2[eZ], zVertexPos1, 1_mm);
0193     zResult2 = result2[eZ];
0194   }
0195
0196   // Both finders should give same results
0197   CHECK_CLOSE_REL(zResult1, zResult2, 1e-5);
0198 }
0199
0200 BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_track_caching_test) {
0201   bool debugMode = false;
0202
0203   const int mainGridSize = 3001;
0204   const int trkGridSize = 35;
0205
0206   Covariance covMat = Covariance::Identity();
0207
0208   // Perigee surface for track parameters
0209   Vector3 pos0{0, 0, 0};
0210   std::shared_ptr<PerigeeSurface> perigeeSurface =
0211       Surface::makeShared<PerigeeSurface>(pos0);
0212
0213   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0214
0215   using Finder1 = GridDensityVertexFinder;
0216   using GridDensity = GaussianGridTrackDensity;
0217
0218   // Use custom grid density here
0219   GridDensity::Config densityConfig(100_mm, mainGridSize, trkGridSize);
0220   densityConfig.useHighestSumZPosition = true;
0221   GridDensity density(densityConfig);
0222
0223   Finder1::Config cfg(density);
0224   cfg.cacheGridStateForTrackRemoval = true;
0225   cfg.extractParameters.connect<&InputTrack::extractParameters>();
0226   Finder1 finder1(cfg);
0227
0228   // Use custom grid density here with same bin size as Finder1
0229   AdaptiveGridTrackDensity::Config adaptiveDensityConfig;
0230   adaptiveDensityConfig.spatialTrkGridSizeRange = {trkGridSize, trkGridSize};
0231   adaptiveDensityConfig.spatialBinExtent = 2. / 30.01 * 1_mm;
0232   adaptiveDensityConfig.useHighestSumZPosition = true;
0233   AdaptiveGridTrackDensity adaptiveDensity(adaptiveDensityConfig);
0234
0235   using Finder2 = AdaptiveGridDensityVertexFinder;
0236   Finder2::Config cfg2(adaptiveDensity);
0237   cfg2.cacheGridStateForTrackRemoval = true;
0238   cfg2.extractParameters.connect<&InputTrack::extractParameters>();
0239   Finder2 finder2(cfg2);
0240
0241   int mySeed = 31415;
0242   std::mt19937 gen(mySeed);
0243   unsigned int nTracks = 200;
0244
0245   std::vector<BoundTrackParameters> trackVec;
0246   trackVec.reserve(nTracks);
0247
0248   // Create nTracks tracks for test case
0249   for (unsigned int i = 0; i < nTracks; i++) {
0250     // The position of the particle
0251     Vector3 pos(xdist(gen), ydist(gen), 0);
0252
0253     // Create momentum and charge of track
0254     double pt = pTDist(gen);
0255     double phi = phiDist(gen);
0256     double eta = etaDist(gen);
0257     double charge = std::copysign(1., etaDist(gen));
0258
0259     // project the position on the surface
0260     Vector3 direction = makeDirectionFromPhiEta(phi, eta);
0261     Intersection3D intersection =
0262         perigeeSurface->intersect(geoContext, pos, direction).closest();
0263     pos = intersection.position();
0264
0265     // Produce most of the tracks at near z1 position,
0266     // some near z2. Highest track density then expected at z1
0267     pos[eZ] = ((i % 4) == 0) ? z2dist(gen) : z1dist(gen);
0268
0269     trackVec.push_back(BoundTrackParameters::create(
0270                            geoContext, perigeeSurface, makeVector4(pos, 0),
0271                            direction, charge / pt, covMat,
0272                            ParticleHypothesis::pion())
0273                            .value());
0274   }
0275
0276   std::vector<InputTrack> inputTracks;
0277   for (const auto& trk : trackVec) {
0278     inputTracks.emplace_back(&trk);
0279   }
0280
0281   IVertexFinder::State state1 = finder1.makeState(magFieldContext);
0282   IVertexFinder::State state2 = finder2.makeState(magFieldContext);
0283
0284   double zResult1 = 0;
0285   double zResult2 = 0;
0286
0287   auto res1 = finder1.find(inputTracks, vertexingOptions, state1);
0288   if (!res1.ok()) {
0289     std::cout << res1.error().message() << std::endl;
0290   }
0291   if (res1.ok()) {
0292     BOOST_CHECK(!(*res1).empty());
0293     Vector3 result = (*res1).back().position();
0294     if (debugMode) {
0295       std::cout << "Vertex position after first fill 1: " << result.transpose()
0296                 << std::endl;
0297     }
0298     CHECK_CLOSE_ABS(result[eZ], zVertexPos1, 1_mm);
0299     zResult1 = result[eZ];
0300   }
0301
0302   auto res2 = finder2.find(inputTracks, vertexingOptions, state2);
0303   if (!res2.ok()) {
0304     std::cout << res2.error().message() << std::endl;
0305   }
0306   if (res2.ok()) {
0307     BOOST_CHECK(!(*res2).empty());
0308     Vector3 result = (*res2).back().position();
0309     if (debugMode) {
0310       std::cout << "Vertex position after first fill 2: " << result.transpose()
0311                 << std::endl;
0312     }
0313     CHECK_CLOSE_ABS(result[eZ], zVertexPos1, 1_mm);
0314     zResult2 = result[eZ];
0315   }
0316
0317   CHECK_CLOSE_REL(zResult1, zResult2, 1e-5);
0318
0319   int trkCount = 0;
0320   std::vector<InputTrack> removedTracks;
0321   for (const auto& trk : trackVec) {
0322     if ((trkCount % 4) != 0) {
0323       removedTracks.emplace_back(&trk);
0324     }
0325     trkCount++;
0326   }
0327
0328   state1.as<Finder1::State>().tracksToRemove = removedTracks;
0329   state2.as<Finder2::State>().tracksToRemove = removedTracks;
0330
0331   auto res3 = finder1.find(inputTracks, vertexingOptions, state1);
0332   if (!res3.ok()) {
0333     std::cout << res3.error().message() << std::endl;
0334   }
0335   if (res3.ok()) {
0336     BOOST_CHECK(!(*res3).empty());
0337     Vector3 result = (*res3).back().position();
0338     if (debugMode) {
0339       std::cout
0340           << "Vertex position after removing tracks near first density peak 1: "
0341           << result << std::endl;
0342     }
0343     CHECK_CLOSE_ABS(result[eZ], zVertexPos2, 1_mm);
0344     zResult1 = result[eZ];
0345   }
0346
0347   auto res4 = finder2.find(inputTracks, vertexingOptions, state2);
0348   if (!res4.ok()) {
0349     std::cout << res4.error().message() << std::endl;
0350   }
0351   if (res4.ok()) {
0352     BOOST_CHECK(!(*res4).empty());
0353     Vector3 result = (*res4).back().position();
0354     if (debugMode) {
0355       std::cout
0356           << "Vertex position after removing tracks near first density peak 2: "
0357           << result << std::endl;
0358     }
0359     CHECK_CLOSE_ABS(result[eZ], zVertexPos2, 1_mm);
0360     zResult2 = result[eZ];
0361   }
0362
0363   CHECK_CLOSE_REL(zResult1, zResult2, 1e-5);
0364 }
0365
0366 ///
0367 /// @brief Unit test for GridDensityVertexFinder with seed with estimation
0368 ///
0369 BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_seed_width_test) {
0370   bool debugMode = false;
0371
0372   const int mainGridSize = 3001;
0373   const int trkGridSize = 35;
0374
0375   Covariance covMat = Covariance::Identity();
0376
0377   // Perigee surface for track parameters
0378   Vector3 pos0{0, 0, 0};
0379   std::shared_ptr<PerigeeSurface> perigeeSurface =
0380       Surface::makeShared<PerigeeSurface>(pos0);
0381
0382   VertexingOptions vertexingOptions(geoContext, magFieldContext);
0383   Vertex constraintVtx;
0384   constraintVtx.setCovariance(SquareMatrix3::Identity());
0385   vertexingOptions.constraint = constraintVtx;
0386
0387   using Finder1 = GridDensityVertexFinder;
0388   GaussianGridTrackDensity::Config gDensityConfig(100, mainGridSize,
0389                                                   trkGridSize);
0390   GaussianGridTrackDensity gDensity(gDensityConfig);
0391   Finder1::Config cfg1(gDensity);
0392   cfg1.cacheGridStateForTrackRemoval = false;
0393   cfg1.estimateSeedWidth = true;
0394   cfg1.extractParameters.connect<&InputTrack::extractParameters>();
0395   Finder1 finder1(cfg1);
0396   IVertexFinder::State state1 = finder1.makeState(magFieldContext);
0397
0398   // Use custom grid density here with same bin size as Finder1
0399   AdaptiveGridTrackDensity::Config adaptiveDensityConfig;
0400   adaptiveDensityConfig.spatialTrkGridSizeRange = {trkGridSize, trkGridSize};
0401   adaptiveDensityConfig.spatialBinExtent = 2. / 30.01 * 1_mm;
0402   AdaptiveGridTrackDensity adaptiveDensity(adaptiveDensityConfig);
0403
0404   using Finder2 = AdaptiveGridDensityVertexFinder;
0405   Finder2::Config cfg2(adaptiveDensity);
0406   cfg2.cacheGridStateForTrackRemoval = false;
0407   cfg2.estimateSeedWidth = true;
0408   cfg2.extractParameters.connect<&InputTrack::extractParameters>();
0409   Finder2 finder2(cfg2);
0410   IVertexFinder::State state2 = finder2.makeState(magFieldContext);
0411
0412   int mySeed = 31415;
0413   std::mt19937 gen(mySeed);
0414   unsigned int nTracks = 20;
0415
0416   std::vector<BoundTrackParameters> trackVec;
0417   trackVec.reserve(nTracks);
0418
0419   // Create nTracks tracks for test case
0420   for (unsigned int i = 0; i < nTracks; i++) {
0421     // The position of the particle
0422     Vector3 pos(xdist(gen), ydist(gen), 0);
0423
0424     // Create momentum and charge of track
0425     double pt = pTDist(gen);
0426     double phi = phiDist(gen);
0427     double eta = etaDist(gen);
0428     double charge = std::copysign(1., etaDist(gen));
0429
0430     // project the position on the surface
0431     Vector3 direction = makeDirectionFromPhiEta(phi, eta);
0432     Intersection3D intersection =
0433         perigeeSurface->intersect(geoContext, pos, direction).closest();
0434     pos = intersection.position();
0435
0436     pos[eZ] = z1dist(gen);
0437
0438     trackVec.push_back(BoundTrackParameters::create(
0439                            geoContext, perigeeSurface, makeVector4(pos, 0),
0440                            direction, charge / pt, covMat,
0441                            ParticleHypothesis::pion())
0442                            .value());
0443   }
0444
0445   std::vector<InputTrack> inputTracks;
0446   for (const auto& trk : trackVec) {
0447     inputTracks.emplace_back(&trk);
0448   }
0449
0450   // Test finder 1
0451   auto res1 = finder1.find(inputTracks, vertexingOptions, state1);
0452   if (!res1.ok()) {
0453     std::cout << res1.error().message() << std::endl;
0454   }
0455
0456   double covZZ1 = 0;
0457   if (res1.ok()) {
0458     BOOST_CHECK(!(*res1).empty());
0459     SquareMatrix3 cov = (*res1).back().covariance();
0460     BOOST_CHECK_NE(constraintVtx.covariance(), cov);
0461     BOOST_CHECK_NE(cov(eZ, eZ), 0.);
0462     covZZ1 = cov(eZ, eZ);
0463     if (debugMode) {
0464       std::cout << "Estimated z-seed width 1: " << cov(eZ, eZ) << std::endl;
0465     }
0466   }
0467
0468   // Test finder 2
0469   auto res2 = finder2.find(inputTracks, vertexingOptions, state2);
0470   if (!res2.ok()) {
0471     std::cout << res2.error().message() << std::endl;
0472   }
0473
0474   double covZZ2 = 0;
0475   if (res2.ok()) {
0476     BOOST_CHECK(!(*res2).empty());
0477     SquareMatrix3 cov = (*res2).back().covariance();
0478     BOOST_CHECK_NE(constraintVtx.covariance(), cov);
0479     BOOST_CHECK_NE(cov(eZ, eZ), 0.);
0480     covZZ2 = cov(eZ, eZ);
0481     if (debugMode) {
0482       std::cout << "Estimated z-seed width 2: " << cov(eZ, eZ) << std::endl;
0483     }
0484   }
0485
0486   // Test for same seed width
0487   CHECK_CLOSE_ABS(covZZ1, covZZ2, 1e-4);
0488 }
0489
0490 BOOST_AUTO_TEST_SUITE_END()
0491
0492 }  // namespace ActsTests