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