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