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   Intersection3D intersection =
0203       surface
0204           .intersect(GeometryContext{}, position, direction,
0205                      BoundaryTolerance::Infinite())
0206           .closest();
0207   mean.head<3>() = intersection.position();
0208
0209   return *transformFreeToBoundParameters(mean, surface, GeometryContext{});
0210 }
0211
0212 // Typedef to describe local positions of 4 components
0213 using LocPosArray = std::array<std::pair<double, double>, 4>;
0214
0215 // Test the combination for a surface type. The local positions are given from
0216 // the outside since their meaning differs between surface types
0217 template <typename angle_description_t>
0218 void test_surface(const Surface &surface, const angle_description_t &desc,
0219                   const LocPosArray &loc_pos, double expectedError) {
0220   const auto proj = std::identity{};
0221
0222   for (auto phi : {-175_degree, 0_degree, 175_degree}) {
0223     for (auto theta : {5_degree, 90_degree, 175_degree}) {
0224       // Go create mixture with 4 cmps
0225       std::vector<DummyComponent<eBoundSize>> cmps;
0226
0227       auto p_it = loc_pos.begin();
0228
0229       for (auto dphi : {-10_degree, 10_degree}) {
0230         for (auto dtheta : {-5_degree, 5_degree}) {
0231           DummyComponent<eBoundSize> a;
0232           a.weight = 1. / 4.;
0233           a.boundPars = BoundVector::Ones();
0234           a.boundPars[eBoundLoc0] *= p_it->first;
0235           a.boundPars[eBoundLoc1] *= p_it->second;
0236           a.boundPars[eBoundPhi] = detail::wrap_periodic(
0237               phi + dphi, -std::numbers::pi, 2 * std::numbers::pi);
0238           a.boundPars[eBoundTheta] = theta + dtheta;
0239
0240           // We don't look at covariance in this test
0241           a.boundCov = BoundSquareMatrix::Zero();
0242
0243           cmps.push_back(a);
0244           ++p_it;
0245         }
0246       }
0247
0248       const auto [mean_approx, cov_approx] =
0249           detail::gaussianMixtureMeanCov(cmps, proj, desc);
0250
0251       const auto mean_ref = meanFromFree(cmps, surface);
0252
0253       CHECK_CLOSE_MATRIX(mean_approx, mean_ref, expectedError);
0254     }
0255   }
0256 }
0257
0258 namespace ActsTests {
0259
0260 BOOST_AUTO_TEST_SUITE(TrackFittingSuite)
0261
0262 BOOST_AUTO_TEST_CASE(test_with_data) {
0263   std::mt19937 gen(42);
0264   std::vector<DummyComponent<2>> cmps(2);
0265
0266   cmps[0].boundPars << 1.0, 1.0;
0267   cmps[0].boundCov << 1.0, 0.0, 0.0, 1.0;
0268   cmps[0].weight = 0.5;
0269
0270   cmps[1].boundPars << -2.0, -2.0;
0271   cmps[1].boundCov << 1.0, 1.0, 1.0, 2.0;
0272   cmps[1].weight = 0.5;
0273
0274   const auto samples = sampleFromMultivariate(cmps, 10000, gen);
0275   const auto mean_data = mean(samples);
0276   const auto boundCov_data = boundCov(samples, mean_data);
0277
0278   const auto [mean_test, boundCov_test] =
0279       detail::gaussianMixtureMeanCov(cmps, std::identity{}, std::tuple<>{});
0280
0281   CHECK_CLOSE_MATRIX(mean_data, mean_test, 1.e-1);
0282   CHECK_CLOSE_MATRIX(boundCov_data, boundCov_test, 1.e-1);
0283 }
0284
0285 BOOST_AUTO_TEST_CASE(test_with_data_circular) {
0286   std::mt19937 gen(42);
0287   std::vector<DummyComponent<2>> cmps(2);
0288
0289   cmps[0].boundPars << 175_degree, 5_degree;
0290   cmps[0].boundCov << 20_degree, 0.0, 0.0, 20_degree;
0291   cmps[0].weight = 0.5;
0292
0293   cmps[1].boundPars << -175_degree, -5_degree;
0294   cmps[1].boundCov << 20_degree, 20_degree, 20_degree, 40_degree;
0295   cmps[1].weight = 0.5;
0296
0297   const auto samples = sampleFromMultivariate(cmps, 10000, gen);
0298   const auto mean_data = circularMean(samples);
0299   const auto boundCov_data = boundCov(samples, mean_data, [](auto a, auto b) {
0300     Vector2 res = Vector2::Zero();
0301     for (int i = 0; i < 2; ++i) {
0302       res[i] = detail::difference_periodic(a[i], b[i], 2 * std::numbers::pi);
0303     }
0304     return res;
0305   });
0306
0307   using detail::CyclicAngle;
0308   const auto d = std::tuple<CyclicAngle<eBoundLoc0>, CyclicAngle<eBoundLoc1>>{};
0309   const auto [mean_test, boundCov_test] =
0310       detail::gaussianMixtureMeanCov(cmps, std::identity{}, d);
0311
0312   BOOST_CHECK(std::abs(detail::difference_periodic(
0313                   mean_data[0], mean_test[0], 2 * std::numbers::pi)) < 1.e-1);
0314   BOOST_CHECK(std::abs(detail::difference_periodic(
0315                   mean_data[1], mean_test[1], 2 * std::numbers::pi)) < 1.e-1);
0316   CHECK_CLOSE_MATRIX(boundCov_data, boundCov_test, 1.e-1);
0317 }
0318
0319 BOOST_AUTO_TEST_CASE(test_plane_surface) {
0320   const auto desc = detail::AngleDescription<Surface::Plane>::Desc{};
0321
0322   const std::shared_ptr<PlaneSurface> surface =
0323       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0324
0325   const LocPosArray p{{{1, 1}, {1, -1}, {-1, 1}, {-1, -1}}};
0326
0327   test_surface(*surface, desc, p, 1.e-2);
0328 }
0329
0330 BOOST_AUTO_TEST_CASE(test_cylinder_surface) {
0331   const Transform3 trafo = Transform3::Identity();
0332   const double r = 2;
0333   const double halfz = 100;
0334
0335   const auto surface = Surface::makeShared<CylinderSurface>(trafo, r, halfz);
0336
0337   const double z1 = -1, z2 = 1;
0338   const double phi1 = 178_degree, phi2 = -176_degree;
0339
0340   const LocPosArray p{
0341       {{r * phi1, z1}, {r * phi1, -z2}, {r * phi2, z1}, {r * phi2, z2}}};
0342
0343   auto desc = detail::AngleDescription<Surface::Cylinder>::Desc{};
0344   std::get<0>(desc).constant = r;
0345
0346   test_surface(*surface, desc, p, 1.e-2);
0347 }
0348
0349 BOOST_AUTO_TEST_CASE(test_disc_surface) {
0350   const Transform3 trafo = Transform3::Identity();
0351   const auto radius = 1;
0352
0353   const auto surface = Surface::makeShared<DiscSurface>(trafo, 0.0, radius);
0354
0355   const double r1 = 0.4, r2 = 0.8;
0356   const double phi1 = -178_degree, phi2 = 176_degree;
0357
0358   const LocPosArray p{{{r1, phi1}, {r2, phi2}, {r1, phi2}, {r2, phi1}}};
0359
0360   const auto desc = detail::AngleDescription<Surface::Disc>::Desc{};
0361
0362   test_surface(*surface, desc, p, 1.e-2);
0363 }
0364
0365 BOOST_AUTO_TEST_CASE(test_perigee_surface) {
0366   const auto desc = detail::AngleDescription<Surface::Plane>::Desc{};
0367
0368   const auto surface = Surface::makeShared<PerigeeSurface>(Vector3{0, 0, 0});
0369
0370   const auto z = 5;
0371   const auto d = 1;
0372
0373   const LocPosArray p{{{d, z}, {d, -z}, {2 * d, z}, {2 * d, -z}}};
0374
0375   // Here we expect a very bad approximation
0376   test_surface(*surface, desc, p, 1.1);
0377 }
0378
0379 BOOST_AUTO_TEST_SUITE_END()
0380
0381 }  // namespace ActsTests