Core/Vertexing/GridDensityVertexFinderTests.cpp

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