Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:11:35

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