Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:13:00

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/TrackParametrization.hpp"
0013 #include "Acts/EventData/ParticleHypothesis.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/Surfaces/PerigeeSurface.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0018 #include "Acts/Utilities/Result.hpp"
0019 #include "Acts/Vertexing/AdaptiveGridTrackDensity.hpp"
0020 
0021 #include <memory>
0022 #include <numbers>
0023 #include <optional>
0024 #include <utility>
0025 
0026 using namespace Acts::UnitLiterals;
0027 
0028 namespace Acts::Test {
0029 
0030 using Covariance = BoundSquareMatrix;
0031 
0032 Covariance makeRandomCovariance(int seed = 31415) {
0033   std::srand(seed);
0034   Covariance randMat((Covariance::Random() + 1.5 * Covariance::Identity()) *
0035                      0.05);
0036 
0037   // symmetric covariance matrix
0038   Covariance covMat = 0.5 * (randMat + randMat.transpose());
0039 
0040   return covMat;
0041 }
0042 
0043 BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) {
0044   // Using a large track grid so we can choose a small bin size
0045   const std::uint32_t spatialTrkGridSize = 4001;
0046   // Arbitrary (but small) bin size
0047   const double binExtent = 3.1e-4;
0048   // Arbitrary impact parameters
0049   const double d0 = 0.4;
0050   const double z0 = -0.2;
0051   Vector2 impactParameters{d0, z0};
0052 
0053   Covariance covMat = makeRandomCovariance();
0054   SquareMatrix2 subCovMat = covMat.block<2, 2>(0, 0);
0055   BoundVector paramVec;
0056   paramVec << d0, z0, 0, 0, 0, 0;
0057 
0058   // Create perigee surface
0059   std::shared_ptr<PerigeeSurface> perigeeSurface =
0060       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0061 
0062   BoundTrackParameters params1(perigeeSurface, paramVec, covMat,
0063                                ParticleHypothesis::pion());
0064 
0065   AdaptiveGridTrackDensity::Config cfg;
0066   // force track to have exactly spatialTrkGridSize spatial bins for testing
0067   // purposes
0068   cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0069   cfg.spatialBinExtent = binExtent;
0070   AdaptiveGridTrackDensity grid(cfg);
0071 
0072   // Empty map
0073   AdaptiveGridTrackDensity::DensityMap mainDensityMap;
0074 
0075   // Add track
0076   auto trackDensityMap = grid.addTrack(params1, mainDensityMap);
0077 
0078   double relTol = 1e-5;
0079   double small = 1e-5;
0080 
0081   auto gaussian2D = [&](const Vector2& args, const Vector2& mus,
0082                         const SquareMatrix2& sigmas) {
0083     Vector2 diffs = args - mus;
0084     double coef = 1 / std::sqrt(sigmas.determinant());
0085     double expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs);
0086     return coef * std::exp(expo);
0087   };
0088 
0089   for (const auto& [bin, _] : mainDensityMap) {
0090     // Extract variables for better readability
0091     std::int32_t zBin = bin.first;
0092     float density = mainDensityMap.at(bin);
0093     // Argument for 2D gaussian
0094     Vector2 dzVec{0., grid.getBinCenter(zBin, binExtent)};
0095     // Compute correct density...
0096     float correctDensity = gaussian2D(dzVec, impactParameters, subCovMat);
0097     // ... and check if our result is equivalent
0098     CHECK_CLOSE_OR_SMALL(correctDensity, density, relTol, small);
0099   }
0100 
0101   // Analytical maximum of the Gaussian (can be obtained by expressing the
0102   // exponent as (az - b)^2 + c and noting correctMaxZ = b/a)
0103   double correctMaxZ =
0104       -0.5 * (subCovMat(0, 1) + subCovMat(1, 0)) / subCovMat(0, 0) * d0 + z0;
0105   // Analytical FWHM of the Gaussian (result similar to
0106   // https://en.wikipedia.org/wiki/Full_width_at_half_maximum#Normal_distribution
0107   // but the calculation needs to be slightly modified in our case)
0108   double correctFWHM =
0109       2. * std::sqrt(2 * std::numbers::ln2 * subCovMat.determinant() /
0110                      subCovMat(0, 0));
0111 
0112   // Estimate maximum z position and seed width
0113   auto res = grid.getMaxZTPositionAndWidth(mainDensityMap);
0114   BOOST_CHECK(res.ok());
0115 
0116   // Extract variables for better readability...
0117   double maxZ = res.value().first.first;
0118   double fwhm = res.value().second * 2.355;
0119 
0120   // ... and check if they are correct (note: the optimization is not as exact
0121   // as the density values).
0122   double relTolOptimization = 1e-3;
0123   CHECK_CLOSE_REL(correctMaxZ, maxZ, relTolOptimization);
0124   CHECK_CLOSE_REL(correctFWHM, fwhm, relTolOptimization);
0125 }
0126 
0127 BOOST_AUTO_TEST_CASE(
0128     compare_to_analytical_solution_for_single_track_with_time) {
0129   // Number of bins in z- and t-direction
0130   const std::uint32_t spatialTrkGridSize = 401;
0131   const std::uint32_t temporalTrkGridSize = 401;
0132   // Bin extents
0133   const double spatialBinExtent = 3.1e-3;
0134   const double temporalBinExtent = 3.1e-3;
0135   // Arbitrary impact parameters
0136   const double d0 = -0.1;
0137   const double z0 = -0.2;
0138   const double t0 = 0.1;
0139   Vector3 impactParameters{d0, z0, t0};
0140 
0141   // symmetric covariance matrix
0142   Covariance covMat = makeRandomCovariance();
0143 
0144   BoundVector paramVec;
0145   paramVec << d0, z0, 0, 0, 0, t0;
0146   // Create perigee surface
0147   std::shared_ptr<PerigeeSurface> perigeeSurface =
0148       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0149 
0150   BoundTrackParameters params(perigeeSurface, paramVec, covMat,
0151                               ParticleHypothesis::pion());
0152 
0153   ActsSquareMatrix<3> ipCov = params.impactParameterCovariance().value();
0154 
0155   AdaptiveGridTrackDensity::Config cfg;
0156   // force track to have exactly spatialTrkGridSize spatial bins for testing
0157   // purposes
0158   cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0159   cfg.spatialBinExtent = spatialBinExtent;
0160   // force track to have exactly temporalTrkGridSize temporal bins for testing
0161   // purposes
0162   cfg.temporalTrkGridSizeRange = {temporalTrkGridSize, temporalTrkGridSize};
0163   cfg.temporalBinExtent = temporalBinExtent;
0164   cfg.useTime = true;
0165   AdaptiveGridTrackDensity grid(cfg);
0166 
0167   // Empty map
0168   AdaptiveGridTrackDensity::DensityMap mainDensityMap;
0169 
0170   // Add track
0171   auto trackDensityMap = grid.addTrack(params, mainDensityMap);
0172 
0173   double relTol = 1e-5;
0174   double small = 1e-5;
0175 
0176   auto gaussian3D = [&](const Vector3& args, const Vector3& mus,
0177                         const SquareMatrix3& sigmas) {
0178     Vector3 diffs = args - mus;
0179     double coef = 1 / std::sqrt(sigmas.determinant());
0180     double expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs);
0181     return coef * std::exp(expo);
0182   };
0183 
0184   for (const auto& [bin, density] : mainDensityMap) {
0185     // Extract variables for better readability
0186     double z = grid.getBinCenter(bin.first, spatialBinExtent);
0187     double t = grid.getBinCenter(bin.second, temporalBinExtent);
0188     // Argument for 3D gaussian
0189     Vector3 dztVec{0., z, t};
0190 
0191     // Compute correct density...
0192     float correctDensity = gaussian3D(dztVec, impactParameters, ipCov);
0193 
0194     // ... and check if our result is equivalent
0195     CHECK_CLOSE_OR_SMALL(correctDensity, density, relTol, small);
0196   }
0197 
0198   // The analytical calculations of the following can be found here:
0199   // https://acts.readthedocs.io/en/latest/white_papers/gaussian-track-densities.html
0200   // Analytical maximum of the Gaussian
0201   ActsSquareMatrix<3> ipWeights = ipCov.inverse();
0202   double denom =
0203       ipWeights(1, 1) * ipWeights(2, 2) - ipWeights(1, 2) * ipWeights(1, 2);
0204 
0205   double zNom =
0206       ipWeights(0, 1) * ipWeights(2, 2) - ipWeights(0, 2) * ipWeights(1, 2);
0207   double correctMaxZ = zNom / denom * d0 + z0;
0208 
0209   double tNom =
0210       ipWeights(0, 2) * ipWeights(1, 1) - ipWeights(0, 1) * ipWeights(1, 2);
0211   double correctMaxT = tNom / denom * d0 + t0;
0212 
0213   // Analytical FWHM of the Gaussian
0214   double correctFWHM = 2. * std::sqrt(2 * std::numbers::ln2 / ipWeights(1, 1));
0215 
0216   // Estimate maximum z position and seed width
0217   auto res = grid.getMaxZTPositionAndWidth(mainDensityMap);
0218   BOOST_CHECK(res.ok());
0219 
0220   // Extract variables for better readability...
0221   double maxZ = res.value().first.first;
0222   double maxT = res.value().first.second;
0223   double fwhm = res.value().second * 2.355;
0224 
0225   // ... and check if they are correct (note: the optimization is not as exact
0226   // as the density values).
0227   double relTolOptimization = 1e-1;
0228   CHECK_CLOSE_REL(correctMaxZ, maxZ, relTolOptimization);
0229   CHECK_CLOSE_REL(correctMaxT, maxT, relTolOptimization);
0230   CHECK_CLOSE_REL(correctFWHM, fwhm, relTolOptimization);
0231 }
0232 
0233 BOOST_AUTO_TEST_CASE(seed_width_estimation) {
0234   // Dummy track grid size (not needed for this unit test)
0235   const std::uint32_t spatialTrkGridSize = 1;
0236   double binExtent = 2.;
0237   AdaptiveGridTrackDensity::Config cfg;
0238   cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0239   cfg.spatialBinExtent = binExtent;
0240   AdaptiveGridTrackDensity grid(cfg);
0241 
0242   // Empty map
0243   AdaptiveGridTrackDensity::DensityMap mainDensityMap;
0244 
0245   // z-position of the maximum density
0246   double correctMaxZ = -2.;
0247 
0248   // Fill map with a triangular track density.
0249   // We use an isoscele triangle with a maximum density value of 1 and a width
0250   // of 20 mm. The linear approximation we use during the seed width estimation
0251   // should be exact in this case.
0252   for (std::int32_t i = -6; i <= 4; i++) {
0253     mainDensityMap[{i, 0}] =
0254         1.0 - 0.1 * std::abs(correctMaxZ - grid.getBinCenter(i, binExtent));
0255   }
0256 
0257   // Get maximum z position and corresponding seed width
0258   auto res = grid.getMaxZTPositionAndWidth(mainDensityMap);
0259   BOOST_CHECK(res.ok());
0260 
0261   // Check if we found the correct maximum
0262   double maxZ = res.value().first.first;
0263   BOOST_CHECK_EQUAL(correctMaxZ, maxZ);
0264 
0265   // Calculate full width at half maximum from seed width and check if it's
0266   // correct
0267   double fwhm = res.value().second * 2.355;
0268   CHECK_CLOSE_REL(10., fwhm, 1e-5);
0269 }
0270 
0271 BOOST_AUTO_TEST_CASE(track_adding) {
0272   const std::uint32_t spatialTrkGridSize = 15;
0273 
0274   double binExtent = 0.1;  // mm
0275 
0276   AdaptiveGridTrackDensity::Config cfg;
0277   cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0278   cfg.spatialBinExtent = binExtent;
0279   AdaptiveGridTrackDensity grid(cfg);
0280 
0281   // Create some test tracks in such a way that some tracks
0282   //  e.g. overlap and that certain tracks need to be inserted
0283   // between two other tracks
0284   Covariance covMat(Covariance::Identity());
0285 
0286   BoundVector paramVec0;
0287   paramVec0 << 100.0, -0.4, 0, 0, 0, 0;
0288   BoundVector paramVec1;
0289   paramVec1 << 0.01, -0.4, 0, 0, 0, 0;
0290   BoundVector paramVec2;
0291   paramVec2 << 0.01, 10.9, 0, 0, 0, 0;
0292   BoundVector paramVec3;
0293   paramVec3 << 0.01, 0.9, 0, 0, 0, 0;
0294 
0295   // Create perigee surface
0296   std::shared_ptr<PerigeeSurface> perigeeSurface =
0297       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0298 
0299   BoundTrackParameters params0(perigeeSurface, paramVec0, covMat,
0300                                ParticleHypothesis::pion());
0301   BoundTrackParameters params1(perigeeSurface, paramVec1, covMat,
0302                                ParticleHypothesis::pion());
0303   BoundTrackParameters params2(perigeeSurface, paramVec2, covMat,
0304                                ParticleHypothesis::pion());
0305   BoundTrackParameters params3(perigeeSurface, paramVec3, covMat,
0306                                ParticleHypothesis::pion());
0307 
0308   // Empty map
0309   AdaptiveGridTrackDensity::DensityMap mainDensityMap;
0310 
0311   // Track is too far away from z axis and was not added
0312   auto trackDensityMap = grid.addTrack(params0, mainDensityMap);
0313   BOOST_CHECK(mainDensityMap.empty());
0314 
0315   // Track should have been entirely added to both grids
0316   trackDensityMap = grid.addTrack(params1, mainDensityMap);
0317   BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap.size());
0318 
0319   // Track should have been entirely added to both grids
0320   trackDensityMap = grid.addTrack(params2, mainDensityMap);
0321   BOOST_CHECK_EQUAL(2 * spatialTrkGridSize, mainDensityMap.size());
0322 
0323   // Track 3 has overlap of 2 bins with track 1
0324   trackDensityMap = grid.addTrack(params3, mainDensityMap);
0325   BOOST_CHECK_EQUAL(3 * spatialTrkGridSize - 2, mainDensityMap.size());
0326 
0327   // Add first track again, should *not* introduce new z entries
0328   trackDensityMap = grid.addTrack(params1, mainDensityMap);
0329   BOOST_CHECK_EQUAL(3 * spatialTrkGridSize - 2, mainDensityMap.size());
0330 }
0331 
0332 BOOST_AUTO_TEST_CASE(max_z_t_and_width) {
0333   const std::uint32_t spatialTrkGridSize = 29;
0334   const std::uint32_t temporalTrkGridSize = 29;
0335 
0336   // spatial and temporal bin extent
0337   double binExtent = 0.05;
0338 
0339   // 1D grid of z values
0340   AdaptiveGridTrackDensity::Config cfg1D;
0341   cfg1D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0342   cfg1D.spatialBinExtent = binExtent;
0343   AdaptiveGridTrackDensity grid1D(cfg1D);
0344 
0345   // 2D grid of z and t values
0346   AdaptiveGridTrackDensity::Config cfg2D;
0347   cfg2D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0348   cfg2D.spatialBinExtent = binExtent;
0349   cfg2D.temporalTrkGridSizeRange = {temporalTrkGridSize, temporalTrkGridSize};
0350   cfg2D.temporalBinExtent = binExtent;
0351   cfg2D.useTime = true;
0352   AdaptiveGridTrackDensity grid2D(cfg2D);
0353 
0354   // Create some test tracks
0355   Covariance covMat(Covariance::Identity() * 0.005);
0356 
0357   double z0Trk1 = 0.25;
0358   double t0Trk1 = 0.05;
0359   double z0Trk2 = -10.95;
0360   double t0Trk2 = 0.1;
0361   BoundVector paramVec1;
0362   paramVec1 << 0.02, z0Trk1, 0, 0, 0, t0Trk1;
0363   BoundVector paramVec2;
0364   paramVec2 << 0.01, z0Trk2, 0, 0, 0, t0Trk2;
0365 
0366   // Create perigee surface
0367   std::shared_ptr<PerigeeSurface> perigeeSurface =
0368       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0369 
0370   BoundTrackParameters params1(perigeeSurface, paramVec1, covMat,
0371                                ParticleHypothesis::pion());
0372   BoundTrackParameters params2(perigeeSurface, paramVec2, covMat,
0373                                ParticleHypothesis::pion());
0374 
0375   // Empty maps
0376   AdaptiveGridTrackDensity::DensityMap mainDensityMap1D;
0377   AdaptiveGridTrackDensity::DensityMap mainDensityMap2D;
0378 
0379   // Add first track to spatial grid
0380   auto trackDensityMap = grid1D.addTrack(params1, mainDensityMap1D);
0381   auto firstRes = grid1D.getMaxZTPosition(mainDensityMap1D);
0382   BOOST_CHECK(firstRes.ok());
0383   // Maximum should be at z0Trk1 position ...
0384   BOOST_CHECK_EQUAL(z0Trk1, (*firstRes).first);
0385   // ... and the corresponding time should be set to 0
0386   BOOST_CHECK_EQUAL(0., (*firstRes).second);
0387 
0388   // Add first track to 2D grid
0389   auto trackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D);
0390   auto firstRes2D = grid2D.getMaxZTPosition(mainDensityMap2D);
0391   BOOST_CHECK(firstRes2D.ok());
0392   // Maximum should be at z0Trk1 position ...
0393   BOOST_CHECK_EQUAL(z0Trk1, (*firstRes2D).first);
0394   // ... and the corresponding time should be at t0Trk1
0395   BOOST_CHECK_EQUAL(t0Trk1, (*firstRes2D).second);
0396 
0397   // Add second track to spatial grid
0398   trackDensityMap = grid1D.addTrack(params2, mainDensityMap1D);
0399   // Calculate maximum and the corresponding width
0400   auto secondRes = grid1D.getMaxZTPositionAndWidth(mainDensityMap1D);
0401   BOOST_CHECK(secondRes.ok());
0402   // Trk 2 is closer to z-axis and should thus yield higher density values.
0403   // Therefore, the new maximum is at z0Trk2 ...
0404   CHECK_CLOSE_OR_SMALL(z0Trk2, (*secondRes).first.first, 1e-5, 1e-5);
0405   // ... the corresponding time should be set to 0...
0406   BOOST_CHECK_EQUAL(0., (*secondRes).first.second);
0407   // ... and it should have a positive width
0408   BOOST_CHECK_LT(0, (*secondRes).second);
0409 
0410   // Add second track to 2D grid
0411   trackDensityMap = grid2D.addTrack(params2, mainDensityMap2D);
0412   // Calculate maximum and the corresponding width
0413   auto secondRes2D = grid2D.getMaxZTPositionAndWidth(mainDensityMap2D);
0414   BOOST_CHECK(secondRes2D.ok());
0415   // Trk 2 is closer to z-axis and should thus yield higher density values.
0416   // Therefore, the new maximum is at z0Trk2 ...
0417   CHECK_CLOSE_OR_SMALL(z0Trk2, (*secondRes2D).first.first, 1e-5, 1e-5);
0418   // ... the corresponding time should be at t0Trk2 ...
0419   BOOST_CHECK_EQUAL(t0Trk2, (*secondRes2D).first.second);
0420   // ... and it should have approximately the same width in z direction
0421   CHECK_CLOSE_OR_SMALL((*secondRes2D).second, (*secondRes).second, 1e-5, 1e-5);
0422 }
0423 
0424 BOOST_AUTO_TEST_CASE(highest_density_sum) {
0425   const std::uint32_t spatialTrkGridSize = 29;
0426 
0427   double binExtent = 0.05;  // mm
0428 
0429   AdaptiveGridTrackDensity::Config cfg;
0430   cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0431   cfg.spatialBinExtent = binExtent;
0432   cfg.useHighestSumZPosition = true;
0433 
0434   AdaptiveGridTrackDensity grid(cfg);
0435 
0436   // Create some test tracks
0437   Covariance covMat(Covariance::Identity() * 0.005);
0438 
0439   double z0Trk1 = 0.25;
0440   double z0Trk2 = -10.95;
0441   BoundVector paramVec1;
0442   paramVec1 << 0.01, z0Trk1, 0, 0, 0, 0;
0443   BoundVector paramVec2;
0444   paramVec2 << 0.0095, z0Trk2, 0, 0, 0, 0;
0445 
0446   // Create perigee surface
0447   std::shared_ptr<PerigeeSurface> perigeeSurface =
0448       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0449 
0450   BoundTrackParameters params1(perigeeSurface, paramVec1, covMat,
0451                                ParticleHypothesis::pion());
0452   BoundTrackParameters params2(perigeeSurface, paramVec2, covMat,
0453                                ParticleHypothesis::pion());
0454 
0455   // Empty map
0456   AdaptiveGridTrackDensity::DensityMap mainDensityMap;
0457 
0458   // Fill grid with track densities
0459   auto trackDensityMap = grid.addTrack(params1, mainDensityMap);
0460 
0461   auto res1 = grid.getMaxZTPosition(mainDensityMap);
0462   BOOST_CHECK(res1.ok());
0463   // Maximum should be at z0Trk1 position
0464   BOOST_CHECK_EQUAL(z0Trk1, (*res1).first);
0465 
0466   // Add second track
0467   trackDensityMap = grid.addTrack(params2, mainDensityMap);
0468   auto res2 = grid.getMaxZTPosition(mainDensityMap);
0469   BOOST_CHECK(res2.ok());
0470   // Trk 2 is closer to z-axis and should yield higher density values
0471   // New maximum is therefore at z0Trk2
0472   CHECK_CLOSE_REL(z0Trk2, (*res2).first, 1e-5);
0473 
0474   // Add small density values around the maximum of track 1
0475   mainDensityMap[{4, 0}] += 0.5;
0476   mainDensityMap[{6, 0}] += 0.5;
0477 
0478   auto res3 = grid.getMaxZTPosition(mainDensityMap);
0479   BOOST_CHECK(res3.ok());
0480   // Trk 2 still has the highest peak density value, however, the small
0481   // added densities for track 1 around its maximum should now lead to
0482   // a prediction at z0Trk1 again
0483   BOOST_CHECK_EQUAL(z0Trk1, (*res3).first);
0484 }
0485 
0486 BOOST_AUTO_TEST_CASE(track_removing) {
0487   const std::uint32_t spatialTrkGridSize = 29;
0488   const std::uint32_t temporalTrkGridSize = 29;
0489 
0490   std::uint32_t trkGridSize = spatialTrkGridSize * temporalTrkGridSize;
0491 
0492   // bin extent in z and t direction
0493   double binExtent = 0.05;
0494 
0495   // 1D grid
0496   AdaptiveGridTrackDensity::Config cfg1D;
0497   cfg1D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0498   cfg1D.spatialBinExtent = binExtent;
0499   AdaptiveGridTrackDensity grid1D(cfg1D);
0500 
0501   // 2D grid
0502   AdaptiveGridTrackDensity::Config cfg2D;
0503   cfg2D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize};
0504   cfg2D.spatialBinExtent = binExtent;
0505   cfg2D.temporalTrkGridSizeRange = {temporalTrkGridSize, temporalTrkGridSize};
0506   cfg2D.temporalBinExtent = binExtent;
0507   cfg2D.useTime = true;
0508   AdaptiveGridTrackDensity grid2D(cfg2D);
0509 
0510   // Create some test tracks
0511   Covariance covMat = makeRandomCovariance();
0512 
0513   // Define z0 values for test tracks
0514   double z0Trk1 = -0.45;
0515   double t0Trk1 = -0.15;
0516   double z0Trk2 = -0.25;
0517   double t0Trk2 = t0Trk1;
0518 
0519   BoundVector paramVec0;
0520   paramVec0 << 0.1, z0Trk1, 0, 0, 0, t0Trk1;
0521   BoundVector paramVec1;
0522   paramVec1 << 0.1, z0Trk2, 0, 0, 0, t0Trk2;
0523 
0524   // Create perigee surface
0525   std::shared_ptr<PerigeeSurface> perigeeSurface =
0526       Surface::makeShared<PerigeeSurface>(Vector3(0., 0., 0.));
0527 
0528   BoundTrackParameters params0(perigeeSurface, paramVec0, covMat,
0529                                ParticleHypothesis::pion());
0530   BoundTrackParameters params1(perigeeSurface, paramVec1, covMat,
0531                                ParticleHypothesis::pion());
0532 
0533   // Empty maps
0534   AdaptiveGridTrackDensity::DensityMap mainDensityMap1D;
0535   AdaptiveGridTrackDensity::DensityMap mainDensityMap2D;
0536 
0537   // Lambda for calculating total density
0538   auto densitySum = [](const auto& densityMap) {
0539     double sum = 0.;
0540     for (const auto& [_, density] : densityMap) {
0541       sum += density;
0542     }
0543     return sum;
0544   };
0545 
0546   // Add track 0 to 1D grid
0547   auto firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D);
0548   BOOST_CHECK(!mainDensityMap1D.empty());
0549   // Grid size should match spatialTrkGridSize
0550   BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap1D.size());
0551   double firstDensitySum1D = densitySum(mainDensityMap1D);
0552 
0553   // Add track 0 to 2D grid
0554   auto firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D);
0555   BOOST_CHECK(!mainDensityMap2D.empty());
0556   // Grid size should match spatialTrkGridSize
0557   BOOST_CHECK_EQUAL(trkGridSize, mainDensityMap2D.size());
0558   double firstDensitySum2D = densitySum(mainDensityMap2D);
0559 
0560   // Add track 0 again to 1D grid
0561   firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D);
0562   BOOST_CHECK(!mainDensityMap1D.empty());
0563   // Grid size should still match spatialTrkGridSize
0564   BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap1D.size());
0565   // Calculate new total density ...
0566   double secondDensitySum1D = densitySum(mainDensityMap1D);
0567   // ... and check that it's twice as large as before
0568   BOOST_CHECK_EQUAL(2 * firstDensitySum1D, secondDensitySum1D);
0569 
0570   // Add track 0 again to 2D grid
0571   firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D);
0572   BOOST_CHECK(!mainDensityMap2D.empty());
0573   // Grid size should still match trkGridSize
0574   BOOST_CHECK_EQUAL(trkGridSize, mainDensityMap2D.size());
0575   // Calculate new total density ...
0576   double secondDensitySum2D = densitySum(mainDensityMap2D);
0577   // ... and check that it's twice as large as before
0578   BOOST_CHECK_EQUAL(2 * firstDensitySum2D, secondDensitySum2D);
0579 
0580   // Remove track 0 from 1D grid
0581   grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D);
0582   // Calculate new total density
0583   double thirdDensitySum1D = densitySum(mainDensityMap1D);
0584   // Density should be old one again
0585   BOOST_CHECK_EQUAL(firstDensitySum1D, thirdDensitySum1D);
0586   // Grid size should still match spatialTrkGridSize (removal does not change
0587   // grid size)
0588   BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap1D.size());
0589 
0590   // Remove track 0 from 2D grid
0591   grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D);
0592   // Calculate new total density
0593   double thirdDensitySum2D = densitySum(mainDensityMap2D);
0594   // Density should be old one again
0595   BOOST_CHECK_EQUAL(firstDensitySum2D, thirdDensitySum2D);
0596   // Grid size should still match trkGridSize (removal does not change grid
0597   // size)
0598   BOOST_CHECK_EQUAL(trkGridSize, mainDensityMap2D.size());
0599 
0600   // Add track 1 to 1D grid (overlaps with track 0!)
0601   auto secondTrackDensityMap1D = grid1D.addTrack(params1, mainDensityMap1D);
0602   std::uint32_t nNonOverlappingBins1D =
0603       static_cast<std::uint32_t>(std::abs(z0Trk1 - z0Trk2) / binExtent);
0604   BOOST_CHECK_EQUAL(spatialTrkGridSize + nNonOverlappingBins1D,
0605                     mainDensityMap1D.size());
0606   double fourthDensitySum1D = densitySum(mainDensityMap1D);
0607 
0608   // Add track 1 to 2D grid (overlaps with track 0!)
0609   auto secondTrackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D);
0610   std::uint32_t nNonOverlappingBins2D =
0611       nNonOverlappingBins1D * temporalTrkGridSize;
0612   BOOST_CHECK_EQUAL(trkGridSize + nNonOverlappingBins2D,
0613                     mainDensityMap2D.size());
0614   double fourthDensitySum2D = densitySum(mainDensityMap2D);
0615 
0616   // Remove second track 0 from 1D grid
0617   grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D);
0618   double fifthDensitySum1D = densitySum(mainDensityMap1D);
0619   // Density should match differences of removed tracks
0620   CHECK_CLOSE_REL(fifthDensitySum1D, fourthDensitySum1D - firstDensitySum1D,
0621                   1e-5);
0622 
0623   // Remove second track 0 from 2D grid
0624   grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D);
0625   double fifthDensitySum2D = densitySum(mainDensityMap2D);
0626   // Density should match differences of removed tracks
0627   CHECK_CLOSE_REL(fifthDensitySum2D, fourthDensitySum2D - firstDensitySum2D,
0628                   1e-5);
0629 
0630   // Remove track 1 from 1D grid
0631   grid1D.subtractTrack(secondTrackDensityMap1D, mainDensityMap1D);
0632   // Size should not have changed
0633   BOOST_CHECK_EQUAL(spatialTrkGridSize + nNonOverlappingBins1D,
0634                     mainDensityMap1D.size());
0635   double sixthDensitySum1D = densitySum(mainDensityMap1D);
0636   // 1D grid is now empty after all tracks were removed
0637   CHECK_CLOSE_ABS(0., sixthDensitySum1D, 1e-4);
0638 
0639   // Remove track 1 from 2D grid
0640   grid2D.subtractTrack(secondTrackDensityMap2D, mainDensityMap2D);
0641   // Size should not have changed
0642   BOOST_CHECK_EQUAL(trkGridSize + nNonOverlappingBins2D,
0643                     mainDensityMap2D.size());
0644   double sixthDensitySum2D = densitySum(mainDensityMap2D);
0645   // 2D grid is now empty after all tracks were removed
0646   CHECK_CLOSE_ABS(0., sixthDensitySum2D, 1e-4);
0647 }
0648 
0649 }  // namespace Acts::Test