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