Core/TrackFitting/GsfMixtureReductionTests.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/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0015 #include "Acts/TrackFitting/GsfMixtureReduction.hpp"
0016 #include "Acts/TrackFitting/detail/SymmetricKlDistanceMatrix.hpp"
0017
0018 #include <algorithm>
0019 #include <array>
0020 #include <cstddef>
0021 #include <memory>
0022 #include <numeric>
0023 #include <tuple>
0024 #include <utility>
0025 #include <vector>
0026
0027 using namespace Acts;
0028
0029 BOOST_AUTO_TEST_CASE(test_distance_matrix_min_distance) {
0030   std::vector<GsfComponent> cmps = {
0031       {1. / 3., BoundVector::Constant(-2.), BoundSquareMatrix::Identity()},
0032       {1. / 3., BoundVector::Constant(+0.), BoundSquareMatrix::Identity()},
0033       {1. / 3., BoundVector::Constant(+1.), BoundSquareMatrix::Identity()},
0034       {1. / 3., BoundVector::Constant(+4.), BoundSquareMatrix::Identity()}};
0035
0036   const auto proj = [](auto &a) -> decltype(auto) { return a; };
0037   detail::SymmetricKLDistanceMatrix mat(cmps, proj);
0038
0039   const auto [i, j] = mat.minDistancePair();
0040   BOOST_CHECK_EQUAL(std::min(i, j), 1);
0041   BOOST_CHECK_EQUAL(std::max(i, j), 2);
0042 }
0043
0044 BOOST_AUTO_TEST_CASE(test_distance_matrix_masking) {
0045   std::vector<GsfComponent> cmps = {
0046       {1. / 3., BoundVector::Constant(-2.), BoundSquareMatrix::Identity()},
0047       {1. / 3., BoundVector::Constant(+0.), BoundSquareMatrix::Identity()},
0048       {1. / 3., BoundVector::Constant(+1.), BoundSquareMatrix::Identity()},
0049       {1. / 3., BoundVector::Constant(+4.), BoundSquareMatrix::Identity()}};
0050
0051   const auto proj = [](auto &a) -> decltype(auto) { return a; };
0052   const std::size_t cmp_to_mask = 2;
0053
0054   detail::SymmetricKLDistanceMatrix mat_full(cmps, proj);
0055   mat_full.maskAssociatedDistances(cmp_to_mask);
0056
0057   cmps.erase(cmps.begin() + cmp_to_mask);
0058   detail::SymmetricKLDistanceMatrix mat_small(cmps, proj);
0059
0060   const auto [full_i, full_j] = mat_full.minDistancePair();
0061   const auto [small_i, small_j] = mat_small.minDistancePair();
0062
0063   BOOST_CHECK_EQUAL(std::min(full_i, full_j), 0);
0064   BOOST_CHECK_EQUAL(std::max(full_i, full_j), 1);
0065   BOOST_CHECK_EQUAL(full_i, small_i);
0066   BOOST_CHECK_EQUAL(full_j, small_j);
0067 }
0068
0069 BOOST_AUTO_TEST_CASE(test_distance_matrix_recompute_distance) {
0070   std::vector<GsfComponent> cmps = {
0071       {1. / 3., BoundVector::Constant(-2.), BoundSquareMatrix::Identity()},
0072       {1. / 3., BoundVector::Constant(+0.), BoundSquareMatrix::Identity()},
0073       {1. / 3., BoundVector::Constant(+1.), BoundSquareMatrix::Identity()},
0074       {1. / 3., BoundVector::Constant(+4.), BoundSquareMatrix::Identity()}};
0075
0076   const auto proj = [](auto &a) -> decltype(auto) { return a; };
0077   detail::SymmetricKLDistanceMatrix mat(cmps, proj);
0078
0079   {
0080     const auto [i, j] = mat.minDistancePair();
0081     BOOST_CHECK_EQUAL(std::min(i, j), 1);
0082     BOOST_CHECK_EQUAL(std::max(i, j), 2);
0083   }
0084
0085   cmps[3].boundPars = BoundVector::Constant(0.1);
0086   mat.recomputeAssociatedDistances(3, cmps, proj);
0087
0088   {
0089     const auto [i, j] = mat.minDistancePair();
0090     BOOST_CHECK_EQUAL(std::min(i, j), 1);
0091     BOOST_CHECK_EQUAL(std::max(i, j), 3);
0092   }
0093
0094   cmps[0].boundPars = BoundVector::Constant(1.01);
0095   mat.recomputeAssociatedDistances(0, cmps, proj);
0096
0097   {
0098     const auto [i, j] = mat.minDistancePair();
0099     BOOST_CHECK_EQUAL(std::min(i, j), 0);
0100     BOOST_CHECK_EQUAL(std::max(i, j), 2);
0101   }
0102 }
0103
0104 BOOST_AUTO_TEST_CASE(test_mixture_reduction) {
0105   auto meanAndSumOfWeights = [](const auto &cmps) {
0106     const auto mean = std::accumulate(
0107         cmps.begin(), cmps.end(), Acts::BoundVector::Zero().eval(),
0108         [](auto sum, const auto &cmp) -> Acts::BoundVector {
0109           return sum + cmp.weight * cmp.boundPars;
0110         });
0111
0112     const double sumOfWeights = std::accumulate(
0113         cmps.begin(), cmps.end(), 0.0,
0114         [](auto sum, const auto &cmp) { return sum + cmp.weight; });
0115
0116     return std::make_tuple(mean, sumOfWeights);
0117   };
0118
0119   // Assume that the components are on a generic plane surface
0120   auto surface = Acts::CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0})
0121                      .planeSurface();
0122   const std::size_t NComps = 4;
0123   std::vector<GsfComponent> cmps;
0124
0125   for (auto i = 0ul; i < NComps; ++i) {
0126     GsfComponent a;
0127     a.boundPars = BoundVector::Zero();
0128     a.boundCov = BoundSquareMatrix::Identity();
0129     a.weight = 1.0 / NComps;
0130     cmps.push_back(a);
0131   }
0132
0133   cmps[0].boundPars[eBoundQOverP] = 0.5_GeV;
0134   cmps[1].boundPars[eBoundQOverP] = 1.5_GeV;
0135   cmps[2].boundPars[eBoundQOverP] = 3.5_GeV;
0136   cmps[3].boundPars[eBoundQOverP] = 4.5_GeV;
0137
0138   // Check start properties
0139   const auto [mean0, sumOfWeights0] = meanAndSumOfWeights(cmps);
0140
0141   BOOST_CHECK_CLOSE(mean0[eBoundQOverP], 2.5_GeV, 1.e-8);
0142   BOOST_CHECK_CLOSE(sumOfWeights0, 1.0, 1.e-8);
0143
0144   // Reduce by factor of 2 and check if weights and QoP are correct
0145   Acts::reduceMixtureWithKLDistance(cmps, 2, *surface);
0146
0147   BOOST_CHECK_EQUAL(cmps.size(), 2);
0148
0149   std::ranges::sort(cmps, {},
0150                     [](const auto &c) { return c.boundPars[eBoundQOverP]; });
0151   BOOST_CHECK_CLOSE(cmps[0].boundPars[eBoundQOverP], 1.0_GeV, 1.e-8);
0152   BOOST_CHECK_CLOSE(cmps[1].boundPars[eBoundQOverP], 4.0_GeV, 1.e-8);
0153
0154   const auto [mean1, sumOfWeights1] = meanAndSumOfWeights(cmps);
0155
0156   BOOST_CHECK_CLOSE(mean1[eBoundQOverP], 2.5_GeV, 1.e-8);
0157   BOOST_CHECK_CLOSE(sumOfWeights1, 1.0, 1.e-8);
0158
0159   // Reduce by factor of 2 and check if weights and QoP are correct
0160   Acts::reduceMixtureWithKLDistance(cmps, 1, *surface);
0161
0162   BOOST_CHECK_EQUAL(cmps.size(), 1);
0163   BOOST_CHECK_CLOSE(cmps[0].boundPars[eBoundQOverP], 2.5_GeV, 1.e-8);
0164   BOOST_CHECK_CLOSE(cmps[0].weight, 1.0, 1.e-8);
0165 }
0166
0167 BOOST_AUTO_TEST_CASE(test_weight_cut_reduction) {
0168   auto dummy = Acts::CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0})
0169                    .planeSurface();
0170   std::vector<GsfComponent> cmps;
0171
0172   // weights do not need to be normalized for this test
0173   for (auto w : {1.0, 2.0, 3.0, 4.0}) {
0174     GsfComponent a;
0175     a.boundPars = BoundVector::Zero();
0176     a.boundCov = BoundSquareMatrix::Identity();
0177     a.weight = w;
0178     cmps.push_back(a);
0179   }
0180
0181   Acts::reduceMixtureLargestWeights(cmps, 2, *dummy);
0182
0183   BOOST_CHECK_EQUAL(cmps.size(), 2);
0184   std::ranges::sort(cmps, {}, [](const auto &c) { return c.weight; });
0185
0186   BOOST_CHECK_EQUAL(cmps[0].weight, 3.0);
0187   BOOST_CHECK_EQUAL(cmps[1].weight, 4.0);
0188 }