Core/TrackFitting/GsfComponentMergingTests.cpp

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