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