Core/TrackFitting/GsfMixtureReductionTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2022-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/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/TrackFitting/GsfMixtureReduction.hpp"
0015 #include "Acts/TrackFitting/detail/SymmetricKlDistanceMatrix.hpp"
0016
0017 #include <algorithm>
0018 #include <array>
0019 #include <cstddef>
0020 #include <memory>
0021 #include <numeric>
0022 #include <tuple>
0023 #include <utility>
0024 #include <vector>
0025
0026 using namespace Acts;
0027
0028 BOOST_AUTO_TEST_CASE(test_distance_matrix_min_distance) {
0029   std::vector<GsfComponent> cmps = {
0030       {1. / 3., BoundVector::Constant(-2.), BoundSquareMatrix::Identity()},
0031       {1. / 3., BoundVector::Constant(+0.), BoundSquareMatrix::Identity()},
0032       {1. / 3., BoundVector::Constant(+1.), BoundSquareMatrix::Identity()},
0033       {1. / 3., BoundVector::Constant(+4.), BoundSquareMatrix::Identity()}};
0034
0035   const auto proj = [](auto &a) -> decltype(auto) { return a; };
0036   detail::SymmetricKLDistanceMatrix mat(cmps, proj);
0037
0038   const auto [i, j] = mat.minDistancePair();
0039   BOOST_CHECK_EQUAL(std::min(i, j), 1);
0040   BOOST_CHECK_EQUAL(std::max(i, j), 2);
0041 }
0042
0043 BOOST_AUTO_TEST_CASE(test_distance_matrix_masking) {
0044   std::vector<GsfComponent> cmps = {
0045       {1. / 3., BoundVector::Constant(-2.), BoundSquareMatrix::Identity()},
0046       {1. / 3., BoundVector::Constant(+0.), BoundSquareMatrix::Identity()},
0047       {1. / 3., BoundVector::Constant(+1.), BoundSquareMatrix::Identity()},
0048       {1. / 3., BoundVector::Constant(+4.), BoundSquareMatrix::Identity()}};
0049
0050   const auto proj = [](auto &a) -> decltype(auto) { return a; };
0051   const std::size_t cmp_to_mask = 2;
0052
0053   detail::SymmetricKLDistanceMatrix mat_full(cmps, proj);
0054   mat_full.maskAssociatedDistances(cmp_to_mask);
0055
0056   cmps.erase(cmps.begin() + cmp_to_mask);
0057   detail::SymmetricKLDistanceMatrix mat_small(cmps, proj);
0058
0059   const auto [full_i, full_j] = mat_full.minDistancePair();
0060   const auto [small_i, small_j] = mat_small.minDistancePair();
0061
0062   BOOST_CHECK_EQUAL(std::min(full_i, full_j), 0);
0063   BOOST_CHECK_EQUAL(std::max(full_i, full_j), 1);
0064   BOOST_CHECK_EQUAL(full_i, small_i);
0065   BOOST_CHECK_EQUAL(full_j, small_j);
0066 }
0067
0068 BOOST_AUTO_TEST_CASE(test_distance_matrix_recompute_distance) {
0069   std::vector<GsfComponent> cmps = {
0070       {1. / 3., BoundVector::Constant(-2.), BoundSquareMatrix::Identity()},
0071       {1. / 3., BoundVector::Constant(+0.), BoundSquareMatrix::Identity()},
0072       {1. / 3., BoundVector::Constant(+1.), BoundSquareMatrix::Identity()},
0073       {1. / 3., BoundVector::Constant(+4.), BoundSquareMatrix::Identity()}};
0074
0075   const auto proj = [](auto &a) -> decltype(auto) { return a; };
0076   detail::SymmetricKLDistanceMatrix mat(cmps, proj);
0077
0078   {
0079     const auto [i, j] = mat.minDistancePair();
0080     BOOST_CHECK_EQUAL(std::min(i, j), 1);
0081     BOOST_CHECK_EQUAL(std::max(i, j), 2);
0082   }
0083
0084   cmps[3].boundPars = BoundVector::Constant(0.1);
0085   mat.recomputeAssociatedDistances(3, cmps, proj);
0086
0087   {
0088     const auto [i, j] = mat.minDistancePair();
0089     BOOST_CHECK_EQUAL(std::min(i, j), 1);
0090     BOOST_CHECK_EQUAL(std::max(i, j), 3);
0091   }
0092
0093   cmps[0].boundPars = BoundVector::Constant(1.01);
0094   mat.recomputeAssociatedDistances(0, cmps, proj);
0095
0096   {
0097     const auto [i, j] = mat.minDistancePair();
0098     BOOST_CHECK_EQUAL(std::min(i, j), 0);
0099     BOOST_CHECK_EQUAL(std::max(i, j), 2);
0100   }
0101 }
0102
0103 BOOST_AUTO_TEST_CASE(test_mixture_reduction) {
0104   auto meanAndSumOfWeights = [](const auto &cmps) {
0105     const auto mean = std::accumulate(
0106         cmps.begin(), cmps.end(), Acts::BoundVector::Zero().eval(),
0107         [](auto sum, const auto &cmp) -> Acts::BoundVector {
0108           return sum + cmp.weight * cmp.boundPars;
0109         });
0110
0111     const double sumOfWeights = std::accumulate(
0112         cmps.begin(), cmps.end(), 0.0,
0113         [](auto sum, const auto &cmp) { return sum + cmp.weight; });
0114
0115     return std::make_tuple(mean, sumOfWeights);
0116   };
0117
0118   // Assume that the components are on a generic plane surface
0119   auto surface = Acts::Surface::makeShared<PlaneSurface>(Vector3{0, 0, 0},
0120                                                          Vector3{1, 0, 0});
0121   const std::size_t NComps = 4;
0122   std::vector<GsfComponent> cmps;
0123
0124   for (auto i = 0ul; i < NComps; ++i) {
0125     GsfComponent a;
0126     a.boundPars = BoundVector::Zero();
0127     a.boundCov = BoundSquareMatrix::Identity();
0128     a.weight = 1.0 / NComps;
0129     cmps.push_back(a);
0130   }
0131
0132   cmps[0].boundPars[eBoundQOverP] = 0.5_GeV;
0133   cmps[1].boundPars[eBoundQOverP] = 1.5_GeV;
0134   cmps[2].boundPars[eBoundQOverP] = 3.5_GeV;
0135   cmps[3].boundPars[eBoundQOverP] = 4.5_GeV;
0136
0137   // Check start properties
0138   const auto [mean0, sumOfWeights0] = meanAndSumOfWeights(cmps);
0139
0140   BOOST_CHECK_CLOSE(mean0[eBoundQOverP], 2.5_GeV, 1.e-8);
0141   BOOST_CHECK_CLOSE(sumOfWeights0, 1.0, 1.e-8);
0142
0143   // Reduce by factor of 2 and check if weights and QoP are correct
0144   Acts::reduceMixtureWithKLDistance(cmps, 2, *surface);
0145
0146   BOOST_CHECK_EQUAL(cmps.size(), 2);
0147
0148   std::sort(cmps.begin(), cmps.end(), [](const auto &a, const auto &b) {
0149     return a.boundPars[eBoundQOverP] < b.boundPars[eBoundQOverP];
0150   });
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::Surface::makeShared<PlaneSurface>(Vector3{0, 0, 0},
0169                                                        Vector3{1, 0, 0});
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::sort(cmps.begin(), cmps.end(),
0185             [](const auto &a, const auto &b) { return a.weight < b.weight; });
0186
0187   BOOST_CHECK_EQUAL(cmps[0].weight, 3.0);
0188   BOOST_CHECK_EQUAL(cmps[1].weight, 4.0);
0189 }