Core/TrackFitting/GsfComponentMergingTests.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/EventData/TransformationHelpers.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0017 #include "Acts/Surfaces/CylinderBounds.hpp"
0018 #include "Acts/Surfaces/CylinderSurface.hpp"
0019 #include "Acts/Surfaces/DiscSurface.hpp"
0020 #include "Acts/Surfaces/PerigeeSurface.hpp"
0021 #include "Acts/Surfaces/PlaneSurface.hpp"
0022 #include "Acts/Surfaces/Surface.hpp"
0023 #include "Acts/Surfaces/SurfaceBounds.hpp"
0024 #include "Acts/TrackFitting/GsfComponent.hpp"
0025 #include "Acts/TrackFitting/detail/GsfComponentMerging.hpp"
0026 #include "Acts/Utilities/Intersection.hpp"
0027 #include "Acts/Utilities/Result.hpp"
0028 #include "Acts/Utilities/detail/periodic.hpp"
0029
0030 #include <algorithm>
0031 #include <array>
0032 #include <cmath>
0033 #include <cstddef>
0034 #include <functional>
0035 #include <initializer_list>
0036 #include <memory>
0037 #include <numbers>
0038 #include <random>
0039 #include <stdexcept>
0040 #include <tuple>
0041 #include <utility>
0042 #include <vector>
0043
0044 #include <Eigen/Eigenvalues>
0045
0046 #define CHECK_CLOSE_MATRIX(a, b, t) \
0047   BOOST_CHECK(((a - b).array().abs() < t).all())
0048
0049 using namespace Acts;
0050 using namespace Acts::UnitLiterals;
0051
0052 // A Multivariate distribution object working in the same way as the
0053 // distributions in the standard library
0054 template <typename T, int D>
0055 class MultivariateNormalDistribution {
0056  public:
0057   using Vector = Eigen::Matrix<T, D, 1>;
0058   using Matrix = Eigen::Matrix<T, D, D>;
0059
0060  private:
0061   Vector m_mean;
0062   Matrix m_transform;
0063
0064  public:
0065   MultivariateNormalDistribution(Vector const &mean, Matrix const &boundCov)
0066       : m_mean(mean) {
0067     Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver(boundCov);
0068     m_transform = eigenSolver.eigenvectors() *
0069                   eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
0070   }
0071
0072   template <typename generator_t>
0073   Vector operator()(generator_t &gen) const {
0074     std::normal_distribution<T> normal;
0075     return m_mean +
0076            m_transform * Vector{}.unaryExpr([&](auto) { return normal(gen); });
0077   }
0078 };
0079
0080 // Sample data from a multi-component multivariate distribution
0081 auto sampleFromMultivariate(const std::vector<GsfComponent> &cmps,
0082                             std::size_t n_samples, std::mt19937 &gen) {
0083   using MultiNormal = MultivariateNormalDistribution<double, eBoundSize>;
0084
0085   std::vector<MultiNormal> dists;
0086   std::vector<double> weights;
0087   for (const auto &cmp : cmps) {
0088     dists.push_back(MultiNormal(cmp.boundPars, cmp.boundCov));
0089     weights.push_back(cmp.weight);
0090   }
0091
0092   std::discrete_distribution choice(weights.begin(), weights.end());
0093
0094   auto sample = [&]() {
0095     const auto n = choice(gen);
0096     return dists[n](gen);
0097   };
0098
0099   std::vector<BoundVector> samples(n_samples);
0100   std::generate(samples.begin(), samples.end(), sample);
0101
0102   return samples;
0103 }
0104
0105 // Simple arithmetic mean computation
0106 auto mean(const std::vector<BoundVector> &samples) -> BoundVector {
0107   BoundVector mean = BoundVector::Zero();
0108
0109   for (const auto &x : samples) {
0110     mean += x;
0111   }
0112
0113   return mean / samples.size();
0114 }
0115
0116 // A method to compute the circular mean, since the normal arithmetic mean
0117 // doesn't work for angles in general
0118 auto circularMean(const std::vector<BoundVector> &samples) -> BoundVector {
0119   BoundVector x = BoundVector::Zero();
0120   BoundVector y = BoundVector::Zero();
0121
0122   for (const auto &s : samples) {
0123     for (unsigned int i = 0; i < eBoundSize; ++i) {
0124       x[i] += std::cos(s[i]);
0125       y[i] += std::sin(s[i]);
0126     }
0127   }
0128
0129   BoundVector mean = BoundVector::Zero();
0130
0131   for (unsigned int i = 0; i < eBoundSize; ++i) {
0132     mean[i] = std::atan2(y[i], x[i]);
0133   }
0134
0135   return mean;
0136 }
0137
0138 // This general boundCovariance estimator can be equipped with a custom
0139 // subtraction object to enable circular behaviour
0140 template <typename subtract_t = std::minus<BoundVector>>
0141 auto boundCov(const std::vector<BoundVector> &samples, const BoundVector &mu,
0142               const subtract_t &sub = subtract_t{}) -> BoundMatrix {
0143   BoundMatrix boundCov = BoundMatrix::Zero();
0144   for (const auto &smpl : samples) {
0145     boundCov += sub(smpl, mu) * sub(smpl, mu).transpose();
0146   }
0147
0148   return boundCov / samples.size();
0149 }
0150
0151 // This function computes the mean of a bound gaussian mixture by converting
0152 // them to cartesian coordinates, computing the mean, and converting back to
0153 // bound.
0154 BoundVector meanFromFree(std::vector<GsfComponent> &cmps,
0155                          const Surface &surface) {
0156   // Specially handle LOC0, since the free mean would not be on the surface
0157   // likely
0158   if (surface.type() == Surface::Cylinder) {
0159     auto x = 0.0, y = 0.0;
0160     const auto r = surface.bounds().values()[CylinderBounds::eR];
0161
0162     for (const auto &cmp : cmps) {
0163       x += cmp.weight * std::cos(cmp.boundPars[eBoundLoc0] / r);
0164       y += cmp.weight * std::sin(cmp.boundPars[eBoundLoc0] / r);
0165     }
0166
0167     for (auto &cmp : cmps) {
0168       cmp.boundPars[eBoundLoc0] = std::atan2(y, x) * r;
0169     }
0170   }
0171
0172   if (surface.type() == Surface::Cone) {
0173     throw std::runtime_error("Cone surface not supported");
0174   }
0175
0176   FreeVector mean = FreeVector::Zero();
0177
0178   for (const auto &cmp : cmps) {
0179     mean += cmp.weight * transformBoundToFreeParameters(
0180                              surface,
0181                              GeometryContext::dangerouslyDefaultConstruct(),
0182                              cmp.boundPars);
0183   }
0184
0185   mean.segment<3>(eFreeDir0).normalize();
0186
0187   // Project the position on the surface.
0188   // This is mainly necessary for the perigee surface, where
0189   // the mean might not fulfill the perigee condition.
0190   Vector3 position = mean.head<3>();
0191   Vector3 direction = mean.segment<3>(eFreeDir0);
0192   Intersection3D intersection =
0193       surface
0194           .intersect(GeometryContext::dangerouslyDefaultConstruct(), position,
0195                      direction, BoundaryTolerance::Infinite())
0196           .closest();
0197   mean.head<3>() = intersection.position();
0198
0199   return *transformFreeToBoundParameters(
0200       mean, surface, GeometryContext::dangerouslyDefaultConstruct());
0201 }
0202
0203 // Typedef to describe local positions of 4 components
0204 using LocPosArray = std::array<std::pair<double, double>, 4>;
0205
0206 // Test the combination for a surface type. The local positions are given from
0207 // the outside since their meaning differs between surface types
0208 void test_surface(const Surface &surface, const LocPosArray &loc_pos,
0209                   double expectedError) {
0210   for (auto phi : {-175_degree, 0_degree, 175_degree}) {
0211     for (auto theta : {5_degree, 90_degree, 175_degree}) {
0212       // Go create mixture with 4 cmps
0213       std::vector<GsfComponent> cmps;
0214
0215       auto p_it = loc_pos.begin();
0216
0217       for (auto dphi : {-10_degree, 10_degree}) {
0218         for (auto dtheta : {-5_degree, 5_degree}) {
0219           GsfComponent a;
0220           a.weight = 1. / 4.;
0221           a.boundPars = BoundVector::Ones();
0222           a.boundPars[eBoundLoc0] *= p_it->first;
0223           a.boundPars[eBoundLoc1] *= p_it->second;
0224           a.boundPars[eBoundPhi] = detail::wrap_periodic(
0225               phi + dphi, -std::numbers::pi, 2 * std::numbers::pi);
0226           a.boundPars[eBoundTheta] = theta + dtheta;
0227
0228           // We don't look at covariance in this test
0229           a.boundCov = BoundMatrix::Zero();
0230
0231           cmps.push_back(a);
0232           ++p_it;
0233         }
0234       }
0235
0236       const auto [mean_approx, cov_approx] = detail::Gsf::mergeGaussianMixture(
0237           cmps, surface, ComponentMergeMethod::eMean);
0238
0239       const auto mean_ref = meanFromFree(cmps, surface);
0240
0241       CHECK_CLOSE_MATRIX(mean_approx, mean_ref, expectedError);
0242     }
0243   }
0244 }
0245
0246 namespace ActsTests {
0247
0248 BOOST_AUTO_TEST_SUITE(TrackFittingSuite)
0249
0250 BOOST_AUTO_TEST_CASE(test_with_data) {
0251   const std::shared_ptr<PlaneSurface> surface =
0252       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0253
0254   std::mt19937 gen(42);
0255   std::vector<GsfComponent> cmps(2);
0256
0257   cmps[0].boundPars.head<2>() << 1.0, 1.0;
0258   cmps[0].boundCov.topLeftCorner<2, 2>() << 1.0, 0.0, 0.0, 1.0;
0259   cmps[0].weight = 0.5;
0260
0261   cmps[1].boundPars.head<2>() << -2.0, -2.0;
0262   cmps[1].boundCov.topLeftCorner<2, 2>() << 1.0, 1.0, 1.0, 2.0;
0263   cmps[1].weight = 0.5;
0264
0265   const auto samples = sampleFromMultivariate(cmps, 10000, gen);
0266   const auto mean_data = mean(samples);
0267   const auto boundCov_data = boundCov(samples, mean_data);
0268
0269   const auto [mean_test, boundCov_test] = detail::Gsf::mergeGaussianMixture(
0270       cmps, *surface, ComponentMergeMethod::eMean);
0271
0272   CHECK_CLOSE_MATRIX(mean_data, mean_test, 1.e-1);
0273   CHECK_CLOSE_MATRIX(boundCov_data, boundCov_test, 1.e-1);
0274 }
0275
0276 BOOST_AUTO_TEST_CASE(test_with_data_circular) {
0277   std::mt19937 gen(42);
0278   std::vector<GsfComponent> cmps(2);
0279
0280   cmps[0].boundPars.head<2>() << 175_degree, 5_degree;
0281   cmps[0].boundCov.topLeftCorner<2, 2>() << 20_degree, 0.0, 0.0, 20_degree;
0282   cmps[0].weight = 0.5;
0283
0284   cmps[1].boundPars.head<2>() << -175_degree, -5_degree;
0285   cmps[1].boundCov.topLeftCorner<2, 2>() << 20_degree, 20_degree, 20_degree,
0286       40_degree;
0287   cmps[1].weight = 0.5;
0288
0289   const auto samples = sampleFromMultivariate(cmps, 10000, gen);
0290   const auto mean_data = circularMean(samples);
0291   const auto boundCov_data = boundCov(
0292       samples, mean_data, [](const BoundVector &a, const BoundVector &b) {
0293         BoundVector res = BoundVector::Zero();
0294         for (unsigned int i = 0; i < 2; ++i) {
0295           res[i] =
0296               detail::difference_periodic(a[i], b[i], 2 * std::numbers::pi);
0297         }
0298         return res;
0299       });
0300
0301   using detail::Gsf::CyclicAngle;
0302   const auto d = std::tuple<CyclicAngle<eBoundLoc0>, CyclicAngle<eBoundLoc1>>{};
0303   const auto [mean_test, boundCov_test] = detail::Gsf::mergeGaussianMixture(
0304       cmps,
0305       [](const GsfComponent &c) {
0306         return std::tie(c.weight, c.boundPars, c.boundCov);
0307       },
0308       d, ComponentMergeMethod::eMean);
0309
0310   BOOST_CHECK(std::abs(detail::difference_periodic(
0311                   mean_data[0], mean_test[0], 2 * std::numbers::pi)) < 1.e-1);
0312   BOOST_CHECK(std::abs(detail::difference_periodic(
0313                   mean_data[1], mean_test[1], 2 * std::numbers::pi)) < 1.e-1);
0314   CHECK_CLOSE_MATRIX(boundCov_data, boundCov_test, 1.e-1);
0315 }
0316
0317 BOOST_AUTO_TEST_CASE(test_plane_surface) {
0318   const std::shared_ptr<PlaneSurface> surface =
0319       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0320
0321   const LocPosArray p{{{1, 1}, {1, -1}, {-1, 1}, {-1, -1}}};
0322
0323   test_surface(*surface, p, 1.e-2);
0324 }
0325
0326 BOOST_AUTO_TEST_CASE(test_cylinder_surface) {
0327   const Transform3 trafo = Transform3::Identity();
0328   const double r = 2;
0329   const double halfz = 100;
0330
0331   const auto surface = Surface::makeShared<CylinderSurface>(trafo, r, halfz);
0332
0333   const double z1 = -1, z2 = 1;
0334   const double phi1 = 178_degree, phi2 = -176_degree;
0335
0336   const LocPosArray p{
0337       {{r * phi1, z1}, {r * phi1, -z2}, {r * phi2, z1}, {r * phi2, z2}}};
0338
0339   auto desc = detail::Gsf::AngleDescription<Surface::Cylinder>::Desc{};
0340   std::get<0>(desc).constant = r;
0341
0342   test_surface(*surface, p, 1.e-2);
0343 }
0344
0345 BOOST_AUTO_TEST_CASE(test_disc_surface) {
0346   const Transform3 trafo = Transform3::Identity();
0347   const auto radius = 1;
0348
0349   const auto surface = Surface::makeShared<DiscSurface>(trafo, 0.0, radius);
0350
0351   const double r1 = 0.4, r2 = 0.8;
0352   const double phi1 = -178_degree, phi2 = 176_degree;
0353
0354   const LocPosArray p{{{r1, phi1}, {r2, phi2}, {r1, phi2}, {r2, phi1}}};
0355
0356   test_surface(*surface, p, 1.e-2);
0357 }
0358
0359 BOOST_AUTO_TEST_CASE(test_perigee_surface) {
0360   const auto surface = Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0});
0361
0362   const auto z = 5;
0363   const auto d = 1;
0364
0365   const LocPosArray p{{{d, z}, {d, -z}, {2 * d, z}, {2 * d, -z}}};
0366
0367   // Here we expect a very bad approximation
0368   test_surface(*surface, p, 1.1);
0369 }
0370
0371 BOOST_AUTO_TEST_SUITE_END()
0372
0373 }  // namespace ActsTests