From a4f567b950c2aa178d0cd4ebcaaed4e72db48fef Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Wed, 26 Oct 2022 10:18:40 +0200 Subject: [PATCH 01/23] checkout relevant files --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 765 ++++++++++++++++++ .../detail/BetheHeitlerApprox.hpp | 221 ----- .../TrackFitting/src/GsfFitterFunction.cpp | 78 +- 3 files changed, 839 insertions(+), 225 deletions(-) create mode 100644 Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp delete mode 100644 Core/include/Acts/TrackFitting/detail/BetheHeitlerApprox.hpp diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp new file mode 100644 index 00000000000..200362b7140 --- /dev/null +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -0,0 +1,765 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/TrackFitting/detail/GsfUtils.hpp" + +#include +#include +#include +#include + +#include +#include + +namespace Acts { + +namespace detail { + +struct GaussianComponent { + ActsScalar weight, mean, var; +}; + +/// Transform a gaussian component to a space where all values are defined from +/// [-inf, inf] +void transformComponent(const GaussianComponent &cmp, + double &transformed_weight, double &transformed_mean, + double &transformed_var) { + const auto &[weight, mean, var] = cmp; + + transformed_weight = std::log(weight) - std::log(1 - weight); + transformed_mean = std::log(mean) - std::log(1 - mean); + transformed_var = std::log(var); +} + +/// Transform a gaussian component back from the [-inf, inf]-space to the usual +/// space +auto inverseTransformComponent(double transformed_weight, + double transformed_mean, + double transformed_var) { + GaussianComponent cmp; + cmp.weight = 1. / (1 + std::exp(-transformed_weight)); + cmp.mean = 1. / (1 + std::exp(-transformed_mean)); + cmp.var = std::exp(transformed_var); + + return cmp; +} + +template +class GaussianMixtureModelPDF; + +template +class GaussianMixtureModelCDF { + std::array m_components; + + public: + GaussianMixtureModelCDF( + const std::array &components) + : m_components(components) {} + + double operator()(double x) { + auto cdf = [&](double mu, double sigma_squared) -> double { + return 0.5 * (1.0 + std::erf((x - mu) / (std::sqrt(2 * sigma_squared)))); + }; + + double sum = 0; + + for (const auto [weight, mean, sigma_squared] : m_components) { + sum += weight * cdf(mean, sigma_squared); + } + + return sum; + } + + auto pdf() const { GaussianMixtureModelPDF{m_components}; } +}; + +/// Compute the value of the gaussian mixture at x +template +class GaussianMixtureModelPDF { + std::array m_components; + + public: + GaussianMixtureModelPDF( + const std::array &components) + : m_components(components) {} + + double operator()(double x) { + auto gaussian = [&](double mu, double sigma_squared) -> double { + return (1. / std::sqrt(sigma_squared * 2 * M_PI)) * + std::exp((-(x - mu) * (x - mu)) / (2 * sigma_squared)); + }; + + double sum = 0; + + for (const auto [weight, mean, sigma_squared] : m_components) { + sum += weight * gaussian(mean, sigma_squared); + } + + return sum; + } + + auto cdf() const { return GaussianMixtureModelCDF{m_components}; } +}; + +class BetheHeitlerCDF { + double m_thickness; + + public: + BetheHeitlerCDF(double thicknessInX0) : m_thickness(thicknessInX0) {} + + double operator()(double x) { + if (x <= 0.0) { + return 0.0; + } + if (x >= 1.0) { + return 1.0; + } + + auto c = m_thickness / std::log(2); + return 1. - boost::math::gamma_p(c, -std::log(x)); + } +}; + +/// Compute the Bethe-Heitler Distribution on a value x +class BetheHeitlerPDF { + double m_thickness; + + public: + BetheHeitlerPDF(double thicknessInX0) : m_thickness(thicknessInX0) {} + + double operator()(double x) { + if (x <= 0.0 || x >= 1.0) { + return 0.0; + } + + auto c = m_thickness / std::log(2); + return std::pow(-std::log(x), c - 1) / std::tgamma(c); + } + + BetheHeitlerCDF cdf() const { return BetheHeitlerCDF{m_thickness}; } +}; + +/// Integrand for the CDF distance +template +class CDFIntegrant { + GaussianMixtureModelCDF m_mixture; + BetheHeitlerCDF m_distribution; + + public: + CDFIntegrant(const GaussianMixtureModelCDF &mixture, + const BetheHeitlerCDF &dist) + : m_mixture(mixture), m_distribution(dist) {} + + double operator()(double x) { + return std::abs(m_mixture(x) - m_distribution(x)); + } +}; + +template +class KLIntegrant { + GaussianMixtureModelPDF m_mixture; + BetheHeitlerPDF m_distribution; + + public: + KLIntegrant(const GaussianMixtureModelPDF &mixture, + const BetheHeitlerPDF &dist) + : m_mixture(mixture), m_distribution(dist) {} + + double operator()(double x) { + auto dist_x = m_distribution(x) == 0. ? std::numeric_limits::min() + : m_distribution(x); + return std::log(dist_x / m_mixture(x)) * dist_x; + } +}; + +template +class FlatCache { + using Value = + std::pair>; + std::vector m_cache; + const double m_tolerance; + + public: + FlatCache(double tol) : m_tolerance(tol) {} + + void insert(const Value &val) { + m_cache.push_back(val); + std::sort(m_cache.begin(), m_cache.end(), + [](const auto &a, const auto &b) { return a.first < b.first; }); + } + + std::optional findApprox(double x) const { +#if 1 + auto begin = std::lower_bound( + m_cache.begin(), m_cache.end(), x - m_tolerance, + [](const auto &p, const auto &v) { return p.first < v; }); + auto end = std::upper_bound( + m_cache.begin(), m_cache.end(), x + m_tolerance, + [](const auto &v, const auto &p) { return v < p.first; }); + + if (begin == m_cache.end() or begin == end) { + return std::nullopt; + } else { + auto ret = + std::min_element(begin, end, [&](const auto &a, const auto &b) { + return std::abs(a.first - m_tolerance) < + std::abs(b.first - m_tolerance); + }); + return ret->second; + } +#else + auto cached = std::min_element( + m_cache.begin(), m_cache.end(), [&](const auto &pa, const auto &pb) { + return std::abs(x - pa.first) < std::abs(x - pb.first); + }); + + if (cached != m_cache.end() and std::abs(cached->first - x) < m_tolerance) { + return cached->second; + } else { + return std::nullopt; + } +#endif + } +}; +} // namespace detail + +/// This class approximates the Bethe-Heitler with only one component. The +struct BetheHeitlerApproxSingleCmp { + /// Returns the number of components the returned mixture will have + constexpr auto numComponents() const { return 1; } + + /// Checks if an input is valid for the parameterization. Since this is for + /// debugging, it always returns false + /// + /// @param x input in terms of x/x0 + constexpr bool validXOverX0(ActsScalar) const { return false; } + + /// Returns array with length 1 containing a 1-component-representation of the + /// Bethe-Heitler-Distribution + static auto mixture(const ActsScalar x) { + std::array ret{}; + + ret[0].weight = 1.0; + + const double c = x / std::log(2); + ret[0].mean = std::pow(2, -c); + ret[0].var = std::pow(3, -c) - std::pow(4, -c); + + // ret[0].mean = std::exp(-1. * x); + // ret[0].var = + // std::exp(-1. * x * std::log(3.) / std::log(2.)) - std::exp(-2. * x); + + return ret; + } +}; + +/// This class approximates the Bethe-Heitler distribution as a gaussian +/// mixture. To enable an approximation for continous input variables, the +/// weights, means and variances are internally parametrized as a Nth order +/// polynomial. +template +class AtlasBetheHeitlerApprox { + static_assert(NComponents > 0); + static_assert(PolyDegree > 0); + + public: + struct PolyData { + std::array weightCoeffs, meanCoeffs, varCoeffs; + }; + + using Data = std::array; + + constexpr static double noChangeLimit = 0.0001; + constexpr static double singleGaussianLimit = 0.002; + constexpr static double lowerLimit = 0.10; + constexpr static double higherLimit = 0.20; + + private: + Data m_low_data; + Data m_high_data; + bool m_low_transform; + bool m_high_transform; + + public: + /// Construct the Bethe-Heitler approximation description + /// + /// @param low_data data for the lower x/x0 range + /// @param high_data data for the higher x/x0 range + /// @param transform wether the data need to be transformed (see Atlas code) + constexpr AtlasBetheHeitlerApprox(const Data &low_data, const Data &high_data, + bool low_transform, bool high_transform) + : m_low_data(low_data), + m_high_data(high_data), + m_low_transform(low_transform), + m_high_transform(high_transform) {} + + /// Returns the number of components the returned mixture will have + constexpr auto numComponents() const { return NComponents; } + + /// Checks if an input is valid for the parameterization + /// + /// @param x input in terms of x/x0 + constexpr bool validXOverX0(ActsScalar x) const { return x < higherLimit; } + + /// Generates the mixture from the polynomials and reweights them, so + /// that the sum of all weights is 1 + /// + /// @param x The input in terms of x/x0 (pathlength in terms of radiation length) + auto mixture(ActsScalar x) const { + // Build a polynom + auto poly = [](ActsScalar xx, + const std::array &coeffs) { + ActsScalar sum{0.}; + for (const auto c : coeffs) { + sum = xx * sum + c; + } + throw_assert(std::isfinite(sum), "polynom result not finite"); + return sum; + }; + + // Lambda which builds the components + auto make_mixture = [&](const Data &data, double xx, bool transform) { + // Value initialization should garanuee that all is initialized to zero + std::array ret{}; + ActsScalar weight_sum = 0; + for (int i = 0; i < NComponents; ++i) { + // These transformations must be applied to the data according to ATHENA + // (TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx:79) + if (transform) { + ret[i] = detail::inverseTransformComponent( + poly(xx, data[i].weightCoeffs), poly(xx, data[i].meanCoeffs), + poly(xx, data[i].varCoeffs)); + } else { + ret[i].weight = poly(xx, data[i].weightCoeffs); + ret[i].mean = poly(xx, data[i].meanCoeffs); + ret[i].var = poly(xx, data[i].varCoeffs); + } + + weight_sum += ret[i].weight; + } + + for (int i = 0; i < NComponents; ++i) { + ret[i].weight /= weight_sum; + } + + return ret; + }; + + // Return no change + if (x < noChangeLimit) { + std::array ret{}; + + ret[0].weight = 1.0; + ret[0].mean = 1.0; // p_initial = p_final + ret[0].var = 0.0; + + return ret; + } + // Return single gaussian approximation + if (x < singleGaussianLimit) { + std::array ret{}; + ret[0] = BetheHeitlerApproxSingleCmp::mixture(x)[0]; + return ret; + } + // Return a component representation for lower x0 + if (x < lowerLimit) { + return make_mixture(m_low_data, x, m_low_transform); + } + // Return a component representation for higher x0 + else { + // Cap the x because beyond the parameterization goes wild + const auto high_x = std::min(higherLimit, x); + return make_mixture(m_high_data, high_x, m_high_transform); + } + } + + /// Loads a parameterization from a file according to the Atlas file + /// description + /// + /// @param low_parameters_path Path to the foo.par file that stores + /// the parameterization for low x/x0 + /// @param high_parameters_path Path to the foo.par file that stores + /// the parameterization for high x/x0 + static auto loadFromFile(const std::string &low_parameters_path, + const std::string &high_parameters_path) { + auto read_file = [](const std::string &filepath) { + std::ifstream file(filepath); + + if (!file) { + throw std::invalid_argument("Could not open '" + filepath + "'"); + } + + std::size_t n_cmps, degree; + bool transform_code; + + file >> n_cmps >> degree >> transform_code; + + if (NComponents != n_cmps) { + throw std::invalid_argument("Wrong number of components in '" + + filepath + "'"); + } + + if (PolyDegree != degree) { + throw std::invalid_argument("Wrong wrong polynom order in '" + + filepath + "'"); + } + + if (!transform_code) { + throw std::invalid_argument("Transform-code is required in '" + + filepath + "'"); + } + + Data data; + + for (auto &cmp : data) { + for (auto &coeff : cmp.weightCoeffs) { + file >> coeff; + } + for (auto &coeff : cmp.meanCoeffs) { + file >> coeff; + } + for (auto &coeff : cmp.varCoeffs) { + file >> coeff; + } + } + + return std::make_tuple(data, transform_code); + }; + + const auto [low_data, low_transform] = read_file(low_parameters_path); + const auto [high_data, high_transform] = read_file(high_parameters_path); + + return AtlasBetheHeitlerApprox(low_data, high_data, low_transform, + high_transform); + } +}; + +template +struct DefaultNext { + auto operator()(std::array ps, + std::mt19937 &gen) const { + auto val_dist = std::uniform_real_distribution{-0.5, 0.5}; + for (auto &p : ps) { + p += val_dist(gen); + } + + return ps; + } +}; + +/// This class does the approximation by minimizing the CDF distance +/// individually for each point. Probably very slow, but good vor validating. +template > +class BetheHeitlerSimulatedAnnealingMinimizer { + std::vector m_temperatures; + std::array m_startValue; + std::shared_ptr m_gen; + next_t m_next; + + // save time by caching + mutable detail::FlatCache m_cache; + + public: + BetheHeitlerSimulatedAnnealingMinimizer( + const std::vector &temperatures, + const std::array &startValue, + std::shared_ptr gen, const next_t &next = next_t{}) + : m_temperatures(temperatures), + m_startValue(startValue), + m_gen(gen), + m_next(next), + m_cache(0.001) {} + + // Make a start value that is reasonable by distributing the means with to + // geometric series + static auto makeStartValue() { + std::array mixture{}; + double m = 0.; + for (auto i = 0ul; i < mixture.size(); ++i) { + m += std::pow(0.5, i + 1); + + mixture[i].weight = 1. / (mixture.size() - i); + mixture[i].mean = m; + mixture[i].var = 0.005 / (i + 1); + } + + detail::normalizeWeights(mixture, + [](auto &a) -> double & { return a.weight; }); + return mixture; + } + + /// Returns the number of components the returned mixture will have + constexpr auto numComponents() const { return NComponents; } + + /// Checks if an input is valid for the parameterization. Since we do the fit + /// on-the-fly, always return true (but of course the Bethe-Heitler model does + /// not apply to arbitrary thick materials) + /// + /// @param x input in terms of x/x0 + constexpr bool validXOverX0(ActsScalar) const { return true; } + + /// Performes a simulated-annealing minimization of the CDF distance at the + /// given x/x0. + /// + /// @param x The input in terms of x/x0 (pathlength in terms of radiation length) + auto mixture(const ActsScalar x, + std::vector *history = nullptr) const { + if( history ) { + history->reserve(m_temperatures.size()); + } + + const auto singleCmpApprox = BetheHeitlerApproxSingleCmp::mixture(x)[0]; + const double &E1 = singleCmpApprox.mean; + const double &E2 = singleCmpApprox.var; + + if (auto cached = m_cache.findApprox(x); cached) { + return *cached; + } + + // Helper function to do the integration + auto integrate = [&](const auto &mixture) { + return boost::math::quadrature::trapezoidal( +#if 1 + detail::CDFIntegrant{ + detail::GaussianMixtureModelCDF{mixture}, + detail::BetheHeitlerCDF{x}}, +#else + detail::KLIntegrant{ + detail::GaussianMixtureModelPDF{mixture}, + detail::BetheHeitlerPDF{x}}, + +#endif + -0.5, 1.5); + }; + + // Transform to coordinates defined in [-inf, inf] + auto transform = + [](const std::array &cmps) { + auto ret = std::array{}; + + auto it = ret.begin(); + for (const auto &cmp : cmps) { + detail::transformComponent(cmp, *it, *(it + 1), *(it + 2)); + it += 3; + } + + return ret; + }; + + // Transform from coordinates defined in [-inf, inf] + auto inv_transform = [](const std::array &cs) { + auto ret = std::array{}; + + auto it = cs.begin(); + for (auto &cmp : ret) { + cmp = detail::inverseTransformComponent(*it, *(it + 1), *(it + 2)); + it += 3; + } + + return ret; + }; + + // The annealing function + auto minimize = [&](const auto &start_value) { + // Initialize state + auto current_distance = integrate(start_value); + auto current_params = transform(start_value); + + // seperately keep track of best solution + auto best_distance = current_distance; + auto best_params = start_value; + + for (auto T : m_temperatures) { + if (history) { + history->push_back(best_distance); + } + + const auto new_params = m_next(current_params, *m_gen); + auto trafo_params = inv_transform(new_params); + detail::normalizeWeights(trafo_params, + [](auto &a) -> double & { return a.weight; }); + + std::sort(best_params.begin(), best_params.end(), + [](const auto &a, const auto &b) { return a.mean < b.mean; }); + + const auto sum1 = std::accumulate( + std::next(trafo_params.begin()), trafo_params.end(), 0.0, + [](const auto &s, const auto &c) { return s + c.weight * c.mean; }); + + const auto sum2 = + std::accumulate(std::next(trafo_params.begin()), trafo_params.end(), + 0.0, [](const auto &s, const auto &c) { + return s + c.weight * (c.mean * c.mean + c.var); + }); + + trafo_params[0].mean = (1. / trafo_params[0].weight) * (E1 - sum1); + + const auto weightMeanSquared = trafo_params[0].weight * + trafo_params[0].mean * + trafo_params[0].mean; + trafo_params[0].mean = + 1. / trafo_params[0].weight * (E2 - weightMeanSquared - sum2); + + const double new_distance = integrate(trafo_params); + + if (not std::isfinite(new_distance)) { + continue; + } + + const double p = std::exp(-(new_distance - current_distance) / T); + + if (new_distance < best_distance) { + best_distance = new_distance; + best_params = trafo_params; + } + + if (new_distance < current_distance or + p < std::uniform_real_distribution{0., 1.0}(*m_gen)) { + current_distance = new_distance; + current_params = new_params; + } + } + return std::make_tuple(best_distance, best_params); + }; + + // Correct mean & var +#if 0 + const double E1 = std::exp(-1. * x); + const double E2 = + std::exp(-1. * x * std::log(3.) / std::log(2.)) - std::exp(-2. * x); +#endif + + const double threshold = 0.0025; + + auto [best_distance, best_params] = minimize(m_startValue); +#if 0 + // We want to modify the component with the smallest mean + std::sort(best_params.begin(), best_params.end(), + [](const auto &a, const auto &b) { return a.mean < b.mean; }); + + std::cout << "best distance " << best_distance << " at x/x0 = " << x + << std::endl; + + std::cout << "mean result: " << std::accumulate(best_params.begin(), best_params.end(), 0.0, [](const auto &s, const auto &c){ return s + c.weight * c.mean; }) << "\n\tmeans: "; + std::transform(best_params.begin(), best_params.end(), std::ostream_iterator(std::cout, " "), [](const auto &c){ return c.mean; }); + std::cout << "\n\tweights: "; + std::transform(best_params.begin(), best_params.end(), std::ostream_iterator(std::cout, " "), [](const auto &c){ return c.weight; }); + std::cout << "\nmean true: " << singleCmpApprox.mean << "\n"; + + // Evaluate these summes before modification + const auto sum1 = std::accumulate( + std::next(best_params.begin()), best_params.end(), 0.0, + [](const auto &s, const auto &c) { return s + c.weight * c.mean; }); + + const auto sum2 = + std::accumulate(std::next(best_params.begin()), best_params.end(), 0.0, + [](const auto &s, const auto &c) { + return s + c.weight * (c.mean * c.mean + c.var); + }); + + const auto new_mean = (1. / best_params[0].weight) * (E1 - sum1); + + const auto weightMeanSquared = best_params[0].weight * new_mean * new_mean; + const auto new_var = + 1. / best_params[0].weight * (E2 - weightMeanSquared - sum2); + + auto valid = [](const detail::GaussianComponent &cmp) { + return cmp.var > 0. && cmp.mean < 1. && cmp.var > 0.0 && + std::isfinite(cmp.weight); + }; + + // TODO this can not be the case if the fit is bad for very low thicknesses + if (valid({best_params[0].weight, new_mean, new_var})) { + best_params[0].mean = new_mean; + best_params[0].var = new_var; + std::cout << "successfull corrected moments\n"; + } else { + std::cout << "could not correct moments: " << best_params[0].mean << "+-" << best_params[0].var << " -> " << new_mean << "+-" << new_var << "\n"; + auto foo = best_params; + foo[0].mean = new_mean; + foo[0].var = new_var; + std::cout << "corrected mean result: " << std::accumulate(foo.begin(), foo.end(), 0.0, [](const auto &s, const auto &c){ return s + c.weight * c.mean; }) << "\n"; + } + std::cout << "-----------------------\n"; + + detail::normalizeWeights(best_params, + [](auto &c) -> double & { return c.weight; }); + + throw_assert(std::all_of(best_params.begin(), best_params.end(), valid), + "bal"); + throw_assert(detail::weightsAreNormalized( + best_params, [](const auto &c) { return c.weight; }), + "not normalized"); +#endif + // Store in cache if result is good + if (best_distance < threshold) { + m_cache.insert({x, best_params}); + } + return best_params; + } +}; + +// template +// std::map> +// BetheHeitlerSimulatedAnnealingMinimizer::m_cache = {}; +// +// template +// std::mutex +// BetheHeitlerSimulatedAnnealingMinimizer::m_cache_mutex = {}; + +/// These data are from ATLAS and allow using the GSF without loading files. +/// However, this might not be the optimal parameterization. These data come +/// this file in Athena: +/// Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par +/// These data must be transformed, so construct the AtlasBetheHeitlerApprox +/// with transforms = true +// clang-format off +constexpr static AtlasBetheHeitlerApprox<6, 5>::Data bh_cdf_cmps6_order5_data = {{ + // Component #1 + { + {{3.74397e+004,-1.95241e+004, 3.51047e+003,-2.54377e+002, 1.81080e+001,-3.57643e+000}}, + {{3.56728e+004,-1.78603e+004, 2.81521e+003,-8.93555e+001,-1.14015e+001, 2.55769e-001}}, + {{3.73938e+004,-1.92800e+004, 3.21580e+003,-1.46203e+002,-5.65392e+000,-2.78008e+000}} + }, + // Component #2 + { + {{-4.14035e+004, 2.31883e+004,-4.37145e+003, 2.44289e+002, 1.13098e+001,-3.21230e+000}}, + {{-2.06936e+003, 2.65334e+003,-1.01413e+003, 1.78338e+002,-1.85556e+001, 1.91430e+000}}, + {{-5.19068e+004, 2.55327e+004,-4.22147e+003, 1.90227e+002, 9.34602e+000,-4.80961e+000}} + }, + // Component #3 + { + {{2.52200e+003,-4.86348e+003, 2.11942e+003,-3.84534e+002, 2.94503e+001,-2.83310e+000}}, + {{1.80405e+003,-1.93347e+003, 6.27196e+002,-4.32429e+001,-1.43533e+001, 3.58782e+000}}, + {{-4.61617e+004, 1.78221e+004,-1.95746e+003,-8.80646e+001, 3.43153e+001,-7.57830e+000}} + }, + // Component #4 + { + {{4.94537e+003,-2.08737e+003, 1.78089e+002, 2.29879e+001,-5.52783e+000,-1.86800e+000}}, + {{4.60220e+003,-1.62269e+003,-1.57552e+002, 2.01796e+002,-5.01636e+001, 6.47438e+000}}, + {{-9.50373e+004, 4.05517e+004,-5.62596e+003, 4.58534e+001, 6.70479e+001,-1.22430e+001}} + }, + // Component #5 + { + {{-1.04129e+003, 1.15222e+002,-2.70356e+001, 3.18611e+001,-7.78800e+000,-1.50242e+000}}, + {{-2.71361e+004, 2.00625e+004,-6.19444e+003, 1.10061e+003,-1.29354e+002, 1.08289e+001}}, + {{3.15252e+004,-3.31508e+004, 1.20371e+004,-2.23822e+003, 2.44396e+002,-2.09130e+001}} + }, + // Component #6 + { + {{1.27751e+004,-6.79813e+003, 1.24650e+003,-8.20622e+001,-2.33476e+000, 2.46459e-001}}, + {{3.64336e+005,-2.08457e+005, 4.33028e+004,-3.67825e+003, 4.22914e+001, 1.42701e+001}}, + {{-1.79298e+006, 1.01843e+006,-2.10037e+005, 1.82222e+004,-4.33573e+002,-2.72725e+001}} + }, +}}; +// clang-format on + +} // namespace Acts diff --git a/Core/include/Acts/TrackFitting/detail/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/detail/BetheHeitlerApprox.hpp deleted file mode 100644 index 67b08fa2e95..00000000000 --- a/Core/include/Acts/TrackFitting/detail/BetheHeitlerApprox.hpp +++ /dev/null @@ -1,221 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2021 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/Definitions/Algebra.hpp" - -#include -#include - -namespace Acts { -namespace detail { - -inline ActsScalar logistic_sigmoid(ActsScalar x) { - return 1. / (1 + std::exp(-x)); -} - -struct GaussianMixture { - ActsScalar mean, var, weight; -}; - -/// This class approximates the Bethe-Heitler distribution as a gaussian -/// mixture. To enable an approximation for continous input variables, the -/// weights, means and variances are internally parametrized as a Nth order -/// polynomial. -template -class BetheHeitlerApprox { - static_assert(NComponents > 0); - static_assert(PolyDegree > 0); - - public: - struct PolyData { - std::array weightCoeffs, meanCoeffs, varCoeffs; - }; - - using Data = std::array; - - private: - Data m_low_data; - Data m_high_data; - - ActsScalar poly(ActsScalar x, - const std::array &coeffs) const { - ActsScalar sum(0.); - for (const auto c : coeffs) { - sum = x * sum + c; - } - throw_assert(std::isfinite(sum), "polynom result not finite"); - return sum; - } - - public: - constexpr BetheHeitlerApprox(const Data &data) - : m_low_data(data), m_high_data(data) {} - constexpr BetheHeitlerApprox(const Data &low_data, const Data &high_data) - : m_low_data(low_data), m_high_data(high_data) {} - - /// @brief Returns the number of components the returned mixture will have - constexpr auto numComponents() const { return NComponents; } - - /// @brief Generates the mixture from the polynomials and reweights them, so - /// that the sum of all weights is 1 - auto mixture(const ActsScalar x) const { - // Value initialization should garanuee that all is initialized to zero - - // Some constants - constexpr double singleGaussianRange = 0.0001; - constexpr double lowerRange = 0.002; - constexpr double higherRange = 0.10; - constexpr double maxX0 = 0.20; - - // Lambda which builds the components - auto make_mixture = [&](const Data &data, double xx) { - std::array ret{}; - ActsScalar weight_sum = 0; - for (int i = 0; i < NComponents; ++i) { - // These transformations must be applied to the data according to ATHENA - // (TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx:79) - ret[i].weight = logistic_sigmoid(poly(xx, data[i].weightCoeffs)); - ret[i].mean = logistic_sigmoid(poly(xx, data[i].meanCoeffs)); - ret[i].var = std::exp(poly(xx, data[i].varCoeffs)); - - weight_sum += ret[i].weight; - } - - for (int i = 0; i < NComponents; ++i) { - ret[i].weight /= weight_sum; - } - - return ret; - }; - - // Return no change - if (x < singleGaussianRange) { - std::array ret{}; - - ret[0].weight = 1.0; - ret[0].mean = 1.0; // p_initial = p_final - ret[0].var = 0.0; - - return ret; - } - // Return single gaussian approximation - if (x < lowerRange) { - std::array ret{}; - - ret[0].weight = 1.0; - ret[0].mean = std::exp(-1. * x); - ret[0].var = - std::exp(-1. * x * std::log(3.) / std::log(2.)) - std::exp(-2. * x); - - return ret; - } - // Return a component representation for lower x0 - if (x < higherRange) { - return make_mixture(m_low_data, x); - } - // Return a component representation for higher x0 - else { - const auto high_x = x > maxX0 ? maxX0 : x; - return make_mixture(m_high_data, high_x); - } - } -}; - -template -auto load_bethe_heitler_data(const std::string &filepath) -> - typename BetheHeitlerApprox::Data { - std::ifstream file(filepath); - - if (!file) { - throw std::invalid_argument("Could not open '" + filepath + "'"); - } - - std::size_t n_cmps, order; - bool transform_code; - - file >> n_cmps >> order >> transform_code; - - if (NCmps != n_cmps) { - throw std::invalid_argument("Wrong number of components in '" + filepath + - "'"); - } - - if (Order != order) { - throw std::invalid_argument("Wrong wrong polynom order in '" + filepath + - "'"); - } - - if (!transform_code) { - throw std::invalid_argument("Transform-code is required in '" + filepath + - "'"); - } - - typename BetheHeitlerApprox::Data data; - - for (auto &cmp : data) { - for (auto &coeff : cmp.weightCoeffs) { - file >> coeff; - } - for (auto &coeff : cmp.meanCoeffs) { - file >> coeff; - } - for (auto &coeff : cmp.varCoeffs) { - file >> coeff; - } - } - - return data; -} - -/// These data are from ATLAS -/// (TrkGaussianSumFilter/Data/BetheHeitler_cdf_nC6_O5.par) -// clang-format off -constexpr static BetheHeitlerApprox<6, 5>::Data bh_cdf_cmps6_order5_data = {{ - // Component #1 - { - {{3.74397e+004,-1.95241e+004, 3.51047e+003,-2.54377e+002, 1.81080e+001,-3.57643e+000}}, - {{3.56728e+004,-1.78603e+004, 2.81521e+003,-8.93555e+001,-1.14015e+001, 2.55769e-001}}, - {{3.73938e+004,-1.92800e+004, 3.21580e+003,-1.46203e+002,-5.65392e+000,-2.78008e+000}} - }, - // Component #2 - { - {{-4.14035e+004, 2.31883e+004,-4.37145e+003, 2.44289e+002, 1.13098e+001,-3.21230e+000}}, - {{-2.06936e+003, 2.65334e+003,-1.01413e+003, 1.78338e+002,-1.85556e+001, 1.91430e+000}}, - {{-5.19068e+004, 2.55327e+004,-4.22147e+003, 1.90227e+002, 9.34602e+000,-4.80961e+000}} - }, - // Component #3 - { - {{2.52200e+003,-4.86348e+003, 2.11942e+003,-3.84534e+002, 2.94503e+001,-2.83310e+000}}, - {{1.80405e+003,-1.93347e+003, 6.27196e+002,-4.32429e+001,-1.43533e+001, 3.58782e+000}}, - {{-4.61617e+004, 1.78221e+004,-1.95746e+003,-8.80646e+001, 3.43153e+001,-7.57830e+000}} - }, - // Component #4 - { - {{4.94537e+003,-2.08737e+003, 1.78089e+002, 2.29879e+001,-5.52783e+000,-1.86800e+000}}, - {{4.60220e+003,-1.62269e+003,-1.57552e+002, 2.01796e+002,-5.01636e+001, 6.47438e+000}}, - {{-9.50373e+004, 4.05517e+004,-5.62596e+003, 4.58534e+001, 6.70479e+001,-1.22430e+001}} - }, - // Component #5 - { - {{-1.04129e+003, 1.15222e+002,-2.70356e+001, 3.18611e+001,-7.78800e+000,-1.50242e+000}}, - {{-2.71361e+004, 2.00625e+004,-6.19444e+003, 1.10061e+003,-1.29354e+002, 1.08289e+001}}, - {{3.15252e+004,-3.31508e+004, 1.20371e+004,-2.23822e+003, 2.44396e+002,-2.09130e+001}} - }, - // Component #6 - { - {{1.27751e+004,-6.79813e+003, 1.24650e+003,-8.20622e+001,-2.33476e+000, 2.46459e-001}}, - {{3.64336e+005,-2.08457e+005, 4.33028e+004,-3.67825e+003, 4.22914e+001, 1.42701e+001}}, - {{-1.79298e+006, 1.01843e+006,-2.10037e+005, 1.82222e+004,-4.33573e+002,-2.72725e+001}} - }, -}}; -// clang-format on - -} // namespace detail -} // namespace Acts diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index e1537bcc671..6e438a1fe6a 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -18,19 +18,29 @@ #include "Acts/TrackFitting/GainMatrixSmoother.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" #include "Acts/TrackFitting/GaussianSumFitter.hpp" +#include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "Acts/Utilities/Helpers.hpp" #include "ActsExamples/MagneticField/MagneticField.hpp" +#include + using namespace ActsExamples; namespace { using MultiStepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator; -using Fitter = Acts::GaussianSumFitter; using DirectPropagator = Acts::Propagator; + +#if USE_CUSTOM_BETHE_HEITLER +using BHApprox = Acts::BetheHeitlerSimulatedAnnealingMinimizer +#else +using BHApprox = Acts::AtlasBetheHeitlerApprox<6, 5>; +#endif + +using Fitter = Acts::GaussianSumFitter; using DirectFitter = - Acts::GaussianSumFitter; + Acts::GaussianSumFitter; struct GsfFitterFunctionImpl : public ActsExamples::TrackFittingAlgorithm::TrackFitterFunction { @@ -103,10 +113,70 @@ std::shared_ptr ActsExamples::makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, + std::string lowBetheHeitlerPath, std::string highBetheHeitlerPath, std::size_t maxComponents, bool abortOnError, bool disableAllMaterialHandling) { MultiStepper stepper(std::move(magneticField)); +#if USE_CUSTOM_BETHE_HEITLER + constexpr std::size_t NComponents = 12; + + auto gen = std::make_shared(23465u); + + const auto iterations = 3000; + const auto temperatures = []() { + std::vector t; + for (int i = 0; i < iterations; ++i) { + t.push_back(1 * std::exp(-0.0001 * i)); + } + return t; + }(); + + auto next = [](std::array ps, + std::mt19937& generator) { + auto val_dist = std::uniform_real_distribution{-0.5, 0.5}; + + auto idx = std::uniform_int_distribution(0ul, ps.size())(generator); + ps[idx] += ps[idx] * val_dist(generator); + + return ps; + }; + + const auto startValue = std::array{{ + {1, 0.5, 1.e-5}, + {2, 0.99, 1.e-5}, + {2, 0.99, 1.e-5}, + {2, 0.99, 1.e-5}, + {2,.99, 1.e-5}, + {2,.99, 1.e-5}, + {2, 0.99, 1.e-5}, + {2,.99, 1.e-5}, + {2,.99, 1.e-5}, + {2, 0.99, 1.e-5}, + {2,.99, 1.e-5}, + {2,.99, 1.e-5}, + }}; + + auto bhapp = Acts::BetheHeitlerSimulatedAnnealingMinimizer( + temperatures, startValue, gen, next); +#else + auto makeBehteHeitlerApprox = [&]() { + if (std::filesystem::exists(lowBetheHeitlerPath) && + std::filesystem::exists(highBetheHeitlerPath)) { + return Acts::AtlasBetheHeitlerApprox<6, 5>::loadFromFile( + lowBetheHeitlerPath, highBetheHeitlerPath); + } else { + std::cout + << "WARNING: Could not find files, use standard configuration\n"; + return Acts::AtlasBetheHeitlerApprox<6, 5>(Acts::bh_cdf_cmps6_order5_data, + Acts::bh_cdf_cmps6_order5_data, + true, true); + } + }; + + auto bhapp = makeBehteHeitlerApprox(); +#endif + // Standard fitter Acts::Navigator::Config cfg{trackingGeometry}; cfg.resolvePassive = false; @@ -114,12 +184,12 @@ ActsExamples::makeGsfFitterFunction( cfg.resolveSensitive = true; Acts::Navigator navigator(cfg); Propagator propagator(std::move(stepper), std::move(navigator)); - Fitter trackFitter(std::move(propagator)); + Fitter trackFitter(std::move(propagator), std::move(bhapp)); // Direct fitter Acts::DirectNavigator directNavigator; DirectPropagator directPropagator(stepper, directNavigator); - DirectFitter directTrackFitter(std::move(directPropagator)); + DirectFitter directTrackFitter(std::move(directPropagator), std::move(bhapp)); // build the fitter functions. owns the fitter object. auto fitterFunction = std::make_shared( From 31c7a15749285eb39cc5f6588e583fa70e6561d8 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Wed, 26 Oct 2022 11:31:25 +0200 Subject: [PATCH 02/23] removing unused stuff and adapting --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 453 +----------------- .../Acts/TrackFitting/GaussianSumFitter.hpp | 13 +- Core/include/Acts/TrackFitting/GsfError.hpp | 6 +- Core/include/Acts/TrackFitting/GsfOptions.hpp | 2 + .../Acts/TrackFitting/detail/GsfActor.hpp | 7 +- .../Acts/TrackFitting/detail/GsfSmoothing.hpp | 2 +- Core/src/TrackFitting/GsfError.cpp | 5 +- .../TrackFitting/GsfFitterFunction.hpp | 1 + .../TrackFitting/src/GsfFitterFunction.cpp | 86 +--- Examples/Python/src/TrackFitting.cpp | 2 + 10 files changed, 47 insertions(+), 530 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 200362b7140..e8d1ec441b2 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -52,185 +52,10 @@ auto inverseTransformComponent(double transformed_weight, return cmp; } -template -class GaussianMixtureModelPDF; - -template -class GaussianMixtureModelCDF { - std::array m_components; - - public: - GaussianMixtureModelCDF( - const std::array &components) - : m_components(components) {} - - double operator()(double x) { - auto cdf = [&](double mu, double sigma_squared) -> double { - return 0.5 * (1.0 + std::erf((x - mu) / (std::sqrt(2 * sigma_squared)))); - }; - - double sum = 0; - - for (const auto [weight, mean, sigma_squared] : m_components) { - sum += weight * cdf(mean, sigma_squared); - } - - return sum; - } - - auto pdf() const { GaussianMixtureModelPDF{m_components}; } -}; - -/// Compute the value of the gaussian mixture at x -template -class GaussianMixtureModelPDF { - std::array m_components; - - public: - GaussianMixtureModelPDF( - const std::array &components) - : m_components(components) {} - - double operator()(double x) { - auto gaussian = [&](double mu, double sigma_squared) -> double { - return (1. / std::sqrt(sigma_squared * 2 * M_PI)) * - std::exp((-(x - mu) * (x - mu)) / (2 * sigma_squared)); - }; - - double sum = 0; - - for (const auto [weight, mean, sigma_squared] : m_components) { - sum += weight * gaussian(mean, sigma_squared); - } - - return sum; - } - - auto cdf() const { return GaussianMixtureModelCDF{m_components}; } -}; - -class BetheHeitlerCDF { - double m_thickness; - - public: - BetheHeitlerCDF(double thicknessInX0) : m_thickness(thicknessInX0) {} - - double operator()(double x) { - if (x <= 0.0) { - return 0.0; - } - if (x >= 1.0) { - return 1.0; - } - - auto c = m_thickness / std::log(2); - return 1. - boost::math::gamma_p(c, -std::log(x)); - } -}; - -/// Compute the Bethe-Heitler Distribution on a value x -class BetheHeitlerPDF { - double m_thickness; - - public: - BetheHeitlerPDF(double thicknessInX0) : m_thickness(thicknessInX0) {} - - double operator()(double x) { - if (x <= 0.0 || x >= 1.0) { - return 0.0; - } - - auto c = m_thickness / std::log(2); - return std::pow(-std::log(x), c - 1) / std::tgamma(c); - } - - BetheHeitlerCDF cdf() const { return BetheHeitlerCDF{m_thickness}; } -}; - -/// Integrand for the CDF distance -template -class CDFIntegrant { - GaussianMixtureModelCDF m_mixture; - BetheHeitlerCDF m_distribution; - - public: - CDFIntegrant(const GaussianMixtureModelCDF &mixture, - const BetheHeitlerCDF &dist) - : m_mixture(mixture), m_distribution(dist) {} - - double operator()(double x) { - return std::abs(m_mixture(x) - m_distribution(x)); - } -}; - -template -class KLIntegrant { - GaussianMixtureModelPDF m_mixture; - BetheHeitlerPDF m_distribution; - - public: - KLIntegrant(const GaussianMixtureModelPDF &mixture, - const BetheHeitlerPDF &dist) - : m_mixture(mixture), m_distribution(dist) {} - - double operator()(double x) { - auto dist_x = m_distribution(x) == 0. ? std::numeric_limits::min() - : m_distribution(x); - return std::log(dist_x / m_mixture(x)) * dist_x; - } -}; - -template -class FlatCache { - using Value = - std::pair>; - std::vector m_cache; - const double m_tolerance; - - public: - FlatCache(double tol) : m_tolerance(tol) {} - - void insert(const Value &val) { - m_cache.push_back(val); - std::sort(m_cache.begin(), m_cache.end(), - [](const auto &a, const auto &b) { return a.first < b.first; }); - } - - std::optional findApprox(double x) const { -#if 1 - auto begin = std::lower_bound( - m_cache.begin(), m_cache.end(), x - m_tolerance, - [](const auto &p, const auto &v) { return p.first < v; }); - auto end = std::upper_bound( - m_cache.begin(), m_cache.end(), x + m_tolerance, - [](const auto &v, const auto &p) { return v < p.first; }); - - if (begin == m_cache.end() or begin == end) { - return std::nullopt; - } else { - auto ret = - std::min_element(begin, end, [&](const auto &a, const auto &b) { - return std::abs(a.first - m_tolerance) < - std::abs(b.first - m_tolerance); - }); - return ret->second; - } -#else - auto cached = std::min_element( - m_cache.begin(), m_cache.end(), [&](const auto &pa, const auto &pb) { - return std::abs(x - pa.first) < std::abs(x - pb.first); - }); - - if (cached != m_cache.end() and std::abs(cached->first - x) < m_tolerance) { - return cached->second; - } else { - return std::nullopt; - } -#endif - } -}; } // namespace detail +namespace Experimental { + /// This class approximates the Bethe-Heitler with only one component. The struct BetheHeitlerApproxSingleCmp { /// Returns the number of components the returned mixture will have @@ -442,279 +267,6 @@ class AtlasBetheHeitlerApprox { } }; -template -struct DefaultNext { - auto operator()(std::array ps, - std::mt19937 &gen) const { - auto val_dist = std::uniform_real_distribution{-0.5, 0.5}; - for (auto &p : ps) { - p += val_dist(gen); - } - - return ps; - } -}; - -/// This class does the approximation by minimizing the CDF distance -/// individually for each point. Probably very slow, but good vor validating. -template > -class BetheHeitlerSimulatedAnnealingMinimizer { - std::vector m_temperatures; - std::array m_startValue; - std::shared_ptr m_gen; - next_t m_next; - - // save time by caching - mutable detail::FlatCache m_cache; - - public: - BetheHeitlerSimulatedAnnealingMinimizer( - const std::vector &temperatures, - const std::array &startValue, - std::shared_ptr gen, const next_t &next = next_t{}) - : m_temperatures(temperatures), - m_startValue(startValue), - m_gen(gen), - m_next(next), - m_cache(0.001) {} - - // Make a start value that is reasonable by distributing the means with to - // geometric series - static auto makeStartValue() { - std::array mixture{}; - double m = 0.; - for (auto i = 0ul; i < mixture.size(); ++i) { - m += std::pow(0.5, i + 1); - - mixture[i].weight = 1. / (mixture.size() - i); - mixture[i].mean = m; - mixture[i].var = 0.005 / (i + 1); - } - - detail::normalizeWeights(mixture, - [](auto &a) -> double & { return a.weight; }); - return mixture; - } - - /// Returns the number of components the returned mixture will have - constexpr auto numComponents() const { return NComponents; } - - /// Checks if an input is valid for the parameterization. Since we do the fit - /// on-the-fly, always return true (but of course the Bethe-Heitler model does - /// not apply to arbitrary thick materials) - /// - /// @param x input in terms of x/x0 - constexpr bool validXOverX0(ActsScalar) const { return true; } - - /// Performes a simulated-annealing minimization of the CDF distance at the - /// given x/x0. - /// - /// @param x The input in terms of x/x0 (pathlength in terms of radiation length) - auto mixture(const ActsScalar x, - std::vector *history = nullptr) const { - if( history ) { - history->reserve(m_temperatures.size()); - } - - const auto singleCmpApprox = BetheHeitlerApproxSingleCmp::mixture(x)[0]; - const double &E1 = singleCmpApprox.mean; - const double &E2 = singleCmpApprox.var; - - if (auto cached = m_cache.findApprox(x); cached) { - return *cached; - } - - // Helper function to do the integration - auto integrate = [&](const auto &mixture) { - return boost::math::quadrature::trapezoidal( -#if 1 - detail::CDFIntegrant{ - detail::GaussianMixtureModelCDF{mixture}, - detail::BetheHeitlerCDF{x}}, -#else - detail::KLIntegrant{ - detail::GaussianMixtureModelPDF{mixture}, - detail::BetheHeitlerPDF{x}}, - -#endif - -0.5, 1.5); - }; - - // Transform to coordinates defined in [-inf, inf] - auto transform = - [](const std::array &cmps) { - auto ret = std::array{}; - - auto it = ret.begin(); - for (const auto &cmp : cmps) { - detail::transformComponent(cmp, *it, *(it + 1), *(it + 2)); - it += 3; - } - - return ret; - }; - - // Transform from coordinates defined in [-inf, inf] - auto inv_transform = [](const std::array &cs) { - auto ret = std::array{}; - - auto it = cs.begin(); - for (auto &cmp : ret) { - cmp = detail::inverseTransformComponent(*it, *(it + 1), *(it + 2)); - it += 3; - } - - return ret; - }; - - // The annealing function - auto minimize = [&](const auto &start_value) { - // Initialize state - auto current_distance = integrate(start_value); - auto current_params = transform(start_value); - - // seperately keep track of best solution - auto best_distance = current_distance; - auto best_params = start_value; - - for (auto T : m_temperatures) { - if (history) { - history->push_back(best_distance); - } - - const auto new_params = m_next(current_params, *m_gen); - auto trafo_params = inv_transform(new_params); - detail::normalizeWeights(trafo_params, - [](auto &a) -> double & { return a.weight; }); - - std::sort(best_params.begin(), best_params.end(), - [](const auto &a, const auto &b) { return a.mean < b.mean; }); - - const auto sum1 = std::accumulate( - std::next(trafo_params.begin()), trafo_params.end(), 0.0, - [](const auto &s, const auto &c) { return s + c.weight * c.mean; }); - - const auto sum2 = - std::accumulate(std::next(trafo_params.begin()), trafo_params.end(), - 0.0, [](const auto &s, const auto &c) { - return s + c.weight * (c.mean * c.mean + c.var); - }); - - trafo_params[0].mean = (1. / trafo_params[0].weight) * (E1 - sum1); - - const auto weightMeanSquared = trafo_params[0].weight * - trafo_params[0].mean * - trafo_params[0].mean; - trafo_params[0].mean = - 1. / trafo_params[0].weight * (E2 - weightMeanSquared - sum2); - - const double new_distance = integrate(trafo_params); - - if (not std::isfinite(new_distance)) { - continue; - } - - const double p = std::exp(-(new_distance - current_distance) / T); - - if (new_distance < best_distance) { - best_distance = new_distance; - best_params = trafo_params; - } - - if (new_distance < current_distance or - p < std::uniform_real_distribution{0., 1.0}(*m_gen)) { - current_distance = new_distance; - current_params = new_params; - } - } - return std::make_tuple(best_distance, best_params); - }; - - // Correct mean & var -#if 0 - const double E1 = std::exp(-1. * x); - const double E2 = - std::exp(-1. * x * std::log(3.) / std::log(2.)) - std::exp(-2. * x); -#endif - - const double threshold = 0.0025; - - auto [best_distance, best_params] = minimize(m_startValue); -#if 0 - // We want to modify the component with the smallest mean - std::sort(best_params.begin(), best_params.end(), - [](const auto &a, const auto &b) { return a.mean < b.mean; }); - - std::cout << "best distance " << best_distance << " at x/x0 = " << x - << std::endl; - - std::cout << "mean result: " << std::accumulate(best_params.begin(), best_params.end(), 0.0, [](const auto &s, const auto &c){ return s + c.weight * c.mean; }) << "\n\tmeans: "; - std::transform(best_params.begin(), best_params.end(), std::ostream_iterator(std::cout, " "), [](const auto &c){ return c.mean; }); - std::cout << "\n\tweights: "; - std::transform(best_params.begin(), best_params.end(), std::ostream_iterator(std::cout, " "), [](const auto &c){ return c.weight; }); - std::cout << "\nmean true: " << singleCmpApprox.mean << "\n"; - - // Evaluate these summes before modification - const auto sum1 = std::accumulate( - std::next(best_params.begin()), best_params.end(), 0.0, - [](const auto &s, const auto &c) { return s + c.weight * c.mean; }); - - const auto sum2 = - std::accumulate(std::next(best_params.begin()), best_params.end(), 0.0, - [](const auto &s, const auto &c) { - return s + c.weight * (c.mean * c.mean + c.var); - }); - - const auto new_mean = (1. / best_params[0].weight) * (E1 - sum1); - - const auto weightMeanSquared = best_params[0].weight * new_mean * new_mean; - const auto new_var = - 1. / best_params[0].weight * (E2 - weightMeanSquared - sum2); - - auto valid = [](const detail::GaussianComponent &cmp) { - return cmp.var > 0. && cmp.mean < 1. && cmp.var > 0.0 && - std::isfinite(cmp.weight); - }; - - // TODO this can not be the case if the fit is bad for very low thicknesses - if (valid({best_params[0].weight, new_mean, new_var})) { - best_params[0].mean = new_mean; - best_params[0].var = new_var; - std::cout << "successfull corrected moments\n"; - } else { - std::cout << "could not correct moments: " << best_params[0].mean << "+-" << best_params[0].var << " -> " << new_mean << "+-" << new_var << "\n"; - auto foo = best_params; - foo[0].mean = new_mean; - foo[0].var = new_var; - std::cout << "corrected mean result: " << std::accumulate(foo.begin(), foo.end(), 0.0, [](const auto &s, const auto &c){ return s + c.weight * c.mean; }) << "\n"; - } - std::cout << "-----------------------\n"; - - detail::normalizeWeights(best_params, - [](auto &c) -> double & { return c.weight; }); - - throw_assert(std::all_of(best_params.begin(), best_params.end(), valid), - "bal"); - throw_assert(detail::weightsAreNormalized( - best_params, [](const auto &c) { return c.weight; }), - "not normalized"); -#endif - // Store in cache if result is good - if (best_distance < threshold) { - m_cache.insert({x, best_params}); - } - return best_params; - } -}; - -// template -// std::map> -// BetheHeitlerSimulatedAnnealingMinimizer::m_cache = {}; -// -// template -// std::mutex -// BetheHeitlerSimulatedAnnealingMinimizer::m_cache_mutex = {}; - /// These data are from ATLAS and allow using the GSF without loading files. /// However, this might not be the optimal parameterization. These data come /// this file in Athena: @@ -762,4 +314,5 @@ constexpr static AtlasBetheHeitlerApprox<6, 5>::Data bh_cdf_cmps6_order5_data = }}; // clang-format on +} // namespace Experimental } // namespace Acts diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index 6f2151af731..f2b6f37adaa 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -14,7 +14,6 @@ #include "Acts/Propagator/StandardAborters.hpp" #include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/TrackFitting/KalmanFitter.hpp" -#include "Acts/TrackFitting/detail/BetheHeitlerApprox.hpp" #include "Acts/TrackFitting/detail/GsfActor.hpp" #include @@ -42,9 +41,12 @@ struct IsMultiComponentBoundParameters> } // namespace detail +namespace Experimental { + /// Gaussian Sum Fitter implementation. /// @tparam propagator_t The propagator type on which the algorithm is built on /// @tparam bethe_heitler_approx_t The type of the Bethe-Heitler-Approximation +/// @tparam traj_t The MultiTrajectory type (backend) /// /// @note This GSF implementation tries to be as compatible to the KalmanFitter /// as possible. However, there are certain differences at the moment: @@ -52,13 +54,11 @@ struct IsMultiComponentBoundParameters> /// * There are only measurement states in the result /// * Passed-again-surfaces is always empty at the moment /// * Probably some more differences which I don't think of at the moment. -template > +template struct GaussianSumFitter { GaussianSumFitter(propagator_t&& propagator, - bethe_heitler_approx_t&& bha = bethe_heitler_approx_t( - detail::bh_cdf_cmps6_order5_data)) - : m_propagator(std::move(propagator)), m_bethe_heitler_approx(bha) {} + bethe_heitler_approx_t&& bha) + : m_propagator(std::move(propagator)), m_bethe_heitler_approx(std::move(bha)) {} /// The propagator instance used by the fit function propagator_t m_propagator; @@ -534,4 +534,5 @@ struct GaussianSumFitter { } }; +} // namespace Experimental } // namespace Acts diff --git a/Core/include/Acts/TrackFitting/GsfError.hpp b/Core/include/Acts/TrackFitting/GsfError.hpp index 01a436fef5e..66b649b4690 100644 --- a/Core/include/Acts/TrackFitting/GsfError.hpp +++ b/Core/include/Acts/TrackFitting/GsfError.hpp @@ -11,6 +11,7 @@ #include namespace Acts { +namespace Experimental { enum class GsfError { // ensure all values are non-zero @@ -25,12 +26,13 @@ enum class GsfError { SmoothingFailed }; -std::error_code make_error_code(Acts::GsfError e); +std::error_code make_error_code(GsfError e); +} // namespace Experimental } // namespace Acts // register with STL namespace std { template <> -struct is_error_code_enum : std::true_type {}; +struct is_error_code_enum : std::true_type {}; } // namespace std diff --git a/Core/include/Acts/TrackFitting/GsfOptions.hpp b/Core/include/Acts/TrackFitting/GsfOptions.hpp index e10f485a3f4..3c85343a4db 100644 --- a/Core/include/Acts/TrackFitting/GsfOptions.hpp +++ b/Core/include/Acts/TrackFitting/GsfOptions.hpp @@ -18,6 +18,7 @@ #include "Acts/Utilities/Logger.hpp" namespace Acts { +namespace Experimental { /// The extensions needed for the GSF template @@ -74,4 +75,5 @@ struct GsfOptions { bool disableAllMaterialHandling = false; }; +} // namespace Experimental } // namespace Acts diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index a32cdf3690c..b8e5d761353 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -18,7 +18,7 @@ #include "Acts/TrackFitting/GsfError.hpp" #include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/TrackFitting/KalmanFitter.hpp" -#include "Acts/TrackFitting/detail/BetheHeitlerApprox.hpp" +#include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "Acts/TrackFitting/detail/GsfSmoothing.hpp" #include "Acts/TrackFitting/detail/GsfUtils.hpp" #include "Acts/TrackFitting/detail/KLMixtureReduction.hpp" @@ -107,7 +107,7 @@ struct GsfActor { std::optional numberMeasurements; /// The extensions - GsfExtensions extensions; + Experimental::GsfExtensions extensions; } m_cfg; /// Stores meta information about the components @@ -189,7 +189,8 @@ struct GsfActor { << result.parentTips.size() << " vs " << stepper.numberComponents(state.stepping)); - return set_error_or_abort(GsfError::ComponentNumberMismatch); + return set_error_or_abort( + Experimental::GsfError::ComponentNumberMismatch); } // There seem to be cases where this is not always after initializing the diff --git a/Core/include/Acts/TrackFitting/detail/GsfSmoothing.hpp b/Core/include/Acts/TrackFitting/detail/GsfSmoothing.hpp index 7335659f689..0f7962b2860 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfSmoothing.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfSmoothing.hpp @@ -68,7 +68,7 @@ auto bayesianSmoothing(component_iterator_t fwdBegin, } if (smoothedState.empty()) { - return ResType(GsfError::SmoothingFailed); + return ResType(Experimental::GsfError::SmoothingFailed); } normalizeWeights(smoothedState, [](auto &tuple) -> decltype(auto) { diff --git a/Core/src/TrackFitting/GsfError.cpp b/Core/src/TrackFitting/GsfError.cpp index 525425e3165..7020c62fe8e 100644 --- a/Core/src/TrackFitting/GsfError.cpp +++ b/Core/src/TrackFitting/GsfError.cpp @@ -17,7 +17,7 @@ class GsfErrorCategory : public std::error_category { // Return what each enum means in text. std::string message(int c) const final { - using Acts::GsfError; + using Acts::Experimental::GsfError; switch (static_cast(c)) { case GsfError::NavigationFailed: @@ -47,7 +47,8 @@ class GsfErrorCategory : public std::error_category { } // namespace -std::error_code Acts::make_error_code(Acts::GsfError e) { +std::error_code Acts::Experimental::make_error_code( + Acts::Experimental::GsfError e) { static GsfErrorCategory c; return {static_cast(e), c}; } diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp index 9ebc73a8810..65e67e0e5b4 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp @@ -16,6 +16,7 @@ std::shared_ptr makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, + std::string lowBetheHeitlerPath, std::string highBetheHeitlerPath, std::size_t maxComponents, bool abortOnError, bool disableAllMaterialHandling); diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index 6e438a1fe6a..e90fe626348 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -15,10 +15,10 @@ #include "Acts/Propagator/Navigator.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/Surfaces/Surface.hpp" +#include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "Acts/TrackFitting/GainMatrixSmoother.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" #include "Acts/TrackFitting/GaussianSumFitter.hpp" -#include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "Acts/Utilities/Helpers.hpp" #include "ActsExamples/MagneticField/MagneticField.hpp" @@ -32,15 +32,14 @@ using MultiStepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator; using DirectPropagator = Acts::Propagator; -#if USE_CUSTOM_BETHE_HEITLER -using BHApprox = Acts::BetheHeitlerSimulatedAnnealingMinimizer -#else -using BHApprox = Acts::AtlasBetheHeitlerApprox<6, 5>; -#endif -using Fitter = Acts::GaussianSumFitter; +using BHApprox = Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>; +using Fitter = + Acts::Experimental::GaussianSumFitter; using DirectFitter = - Acts::GaussianSumFitter; + Acts::Experimental::GaussianSumFitter; struct GsfFitterFunctionImpl : public ActsExamples::TrackFittingAlgorithm::TrackFitterFunction { @@ -59,12 +58,12 @@ struct GsfFitterFunctionImpl auto makeGsfOptions( const ActsExamples::TrackFittingAlgorithm::GeneralFitterOptions& options) const { - Acts::GsfExtensions extensions; + Acts::Experimental::GsfExtensions extensions; extensions.updater.connect< &Acts::GainMatrixUpdater::operator()>( &updater); - Acts::GsfOptions gsfOptions{ + Acts::Experimental::GsfOptions gsfOptions{ options.geoContext, options.magFieldContext, options.calibrationContext, @@ -118,64 +117,19 @@ ActsExamples::makeGsfFitterFunction( bool disableAllMaterialHandling) { MultiStepper stepper(std::move(magneticField)); -#if USE_CUSTOM_BETHE_HEITLER - constexpr std::size_t NComponents = 12; - - auto gen = std::make_shared(23465u); - - const auto iterations = 3000; - const auto temperatures = []() { - std::vector t; - for (int i = 0; i < iterations; ++i) { - t.push_back(1 * std::exp(-0.0001 * i)); - } - return t; - }(); - - auto next = [](std::array ps, - std::mt19937& generator) { - auto val_dist = std::uniform_real_distribution{-0.5, 0.5}; - - auto idx = std::uniform_int_distribution(0ul, ps.size())(generator); - ps[idx] += ps[idx] * val_dist(generator); - - return ps; - }; - - const auto startValue = std::array{{ - {1, 0.5, 1.e-5}, - {2, 0.99, 1.e-5}, - {2, 0.99, 1.e-5}, - {2, 0.99, 1.e-5}, - {2,.99, 1.e-5}, - {2,.99, 1.e-5}, - {2, 0.99, 1.e-5}, - {2,.99, 1.e-5}, - {2,.99, 1.e-5}, - {2, 0.99, 1.e-5}, - {2,.99, 1.e-5}, - {2,.99, 1.e-5}, - }}; - - auto bhapp = Acts::BetheHeitlerSimulatedAnnealingMinimizer( - temperatures, startValue, gen, next); -#else - auto makeBehteHeitlerApprox = [&]() { - if (std::filesystem::exists(lowBetheHeitlerPath) && + const auto bhapp = [&](){ + if( lowBetheHeitlerPath.empty() && highBetheHeitlerPath.empty() ) { + return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>( + Acts::Experimental::bh_cdf_cmps6_order5_data, + Acts::Experimental::bh_cdf_cmps6_order5_data, true, true); + } else if (std::filesystem::exists(lowBetheHeitlerPath) && std::filesystem::exists(highBetheHeitlerPath)) { - return Acts::AtlasBetheHeitlerApprox<6, 5>::loadFromFile( + return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>::loadFromFile( lowBetheHeitlerPath, highBetheHeitlerPath); } else { - std::cout - << "WARNING: Could not find files, use standard configuration\n"; - return Acts::AtlasBetheHeitlerApprox<6, 5>(Acts::bh_cdf_cmps6_order5_data, - Acts::bh_cdf_cmps6_order5_data, - true, true); + throw std::invalid_argument("Paths to bethe heitler parameterization do not exist. Pass empty strings to load a default parameterization"); } - }; - - auto bhapp = makeBehteHeitlerApprox(); -#endif + }(); // Standard fitter Acts::Navigator::Config cfg{trackingGeometry}; @@ -184,12 +138,12 @@ ActsExamples::makeGsfFitterFunction( cfg.resolveSensitive = true; Acts::Navigator navigator(cfg); Propagator propagator(std::move(stepper), std::move(navigator)); - Fitter trackFitter(std::move(propagator), std::move(bhapp)); + Fitter trackFitter(std::move(propagator), BHApprox(bhapp)); // Direct fitter Acts::DirectNavigator directNavigator; DirectPropagator directPropagator(stepper, directNavigator); - DirectFitter directTrackFitter(std::move(directPropagator), std::move(bhapp)); + DirectFitter directTrackFitter(std::move(directPropagator), BHApprox(bhapp)); // build the fitter functions. owns the fitter object. auto fitterFunction = std::make_shared( diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index 7d2af1590f9..08378c2e1bd 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -76,9 +76,11 @@ void addTrackFitting(Context& ctx) { "makeGsfFitterFunction", py::overload_cast, std::shared_ptr, + std::string, std::string, std::size_t, bool, bool>( &ActsExamples::makeGsfFitterFunction), py::arg("trackingGeometry"), py::arg("magneticField"), + py::arg("lowBetheHeitlerPath"), py::arg("highBetheHeitlerPath"), py::arg("maxComponents"), py::arg("abortOnError"), py::arg("disableAllMaterialHandling")); } From f9e7e1f97f0ede2bdd6c302369049d1a207ffe4b Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Wed, 26 Oct 2022 13:44:17 +0200 Subject: [PATCH 03/23] update --- .../Acts/TrackFitting/GaussianSumFitter.hpp | 9 +++++---- .../Acts/TrackFitting/detail/GsfActor.hpp | 2 +- .../TrackFitting/GsfFitterFunction.hpp | 14 ++++++++++++++ .../TrackFitting/src/GsfFitterFunction.cpp | 19 ++++++++++--------- Examples/Python/src/TrackFitting.cpp | 3 +-- 5 files changed, 31 insertions(+), 16 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp index f2b6f37adaa..2c6ce4720fc 100644 --- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp +++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp @@ -54,11 +54,12 @@ namespace Experimental { /// * There are only measurement states in the result /// * Passed-again-surfaces is always empty at the moment /// * Probably some more differences which I don't think of at the moment. -template +template struct GaussianSumFitter { - GaussianSumFitter(propagator_t&& propagator, - bethe_heitler_approx_t&& bha) - : m_propagator(std::move(propagator)), m_bethe_heitler_approx(std::move(bha)) {} + GaussianSumFitter(propagator_t&& propagator, bethe_heitler_approx_t&& bha) + : m_propagator(std::move(propagator)), + m_bethe_heitler_approx(std::move(bha)) {} /// The propagator instance used by the fit function propagator_t m_propagator; diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index b8e5d761353..5e35d615a45 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -15,10 +15,10 @@ #include "Acts/Material/ISurfaceMaterial.hpp" #include "Acts/Surfaces/CylinderSurface.hpp" #include "Acts/Surfaces/Surface.hpp" +#include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "Acts/TrackFitting/GsfError.hpp" #include "Acts/TrackFitting/GsfOptions.hpp" #include "Acts/TrackFitting/KalmanFitter.hpp" -#include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "Acts/TrackFitting/detail/GsfSmoothing.hpp" #include "Acts/TrackFitting/detail/GsfUtils.hpp" #include "Acts/TrackFitting/detail/KLMixtureReduction.hpp" diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp index 65e67e0e5b4..8d071b6fc19 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp @@ -12,6 +12,20 @@ namespace ActsExamples { +/// Makes a fitter function object for the GSF +/// +/// @param trackingGeometry the trackingGeometry for the propagator +/// @param magneticField the magnetic field for the propagator +/// @param lowBetheHeitlerPath path to a parameterization of the bethe heitler +/// distribution suitable for low thicknesses (x/x0 < 0.1). Pass empty string +/// to load a default. +/// @param highBetheHeitlerPath path to a parameterization of the bethe heitler +/// distribution suitable for higher thicknesses (0.1 < x/x0 < 0.2). Pass empty +/// string to load a default. +/// @param maxComponents number of maximum components in the track state +/// @param abortOnError wether to call std::abort on an error +/// @param disableAllMaterialHandling run the GSF like a KF (no energy loss, +/// always 1 component, ...) std::shared_ptr makeGsfFitterFunction( std::shared_ptr trackingGeometry, diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index e90fe626348..2c544ccaca6 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -32,11 +32,10 @@ using MultiStepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator; using DirectPropagator = Acts::Propagator; - using BHApprox = Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>; using Fitter = - Acts::Experimental::GaussianSumFitter; + Acts::Experimental::GaussianSumFitter; using DirectFitter = Acts::Experimental::GaussianSumFitter; @@ -117,17 +116,19 @@ ActsExamples::makeGsfFitterFunction( bool disableAllMaterialHandling) { MultiStepper stepper(std::move(magneticField)); - const auto bhapp = [&](){ - if( lowBetheHeitlerPath.empty() && highBetheHeitlerPath.empty() ) { - return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>( + const auto bhapp = [&]() { + if (lowBetheHeitlerPath.empty() && highBetheHeitlerPath.empty()) { + return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>( Acts::Experimental::bh_cdf_cmps6_order5_data, Acts::Experimental::bh_cdf_cmps6_order5_data, true, true); - } else if (std::filesystem::exists(lowBetheHeitlerPath) && - std::filesystem::exists(highBetheHeitlerPath)) { + } else if (std::filesystem::exists(lowBetheHeitlerPath) && + std::filesystem::exists(highBetheHeitlerPath)) { return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>::loadFromFile( lowBetheHeitlerPath, highBetheHeitlerPath); } else { - throw std::invalid_argument("Paths to bethe heitler parameterization do not exist. Pass empty strings to load a default parameterization"); + throw std::invalid_argument( + "Paths to bethe heitler parameterization do not exist. Pass empty " + "strings to load a default parameterization"); } }(); diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index 08378c2e1bd..444c035dbff 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -76,8 +76,7 @@ void addTrackFitting(Context& ctx) { "makeGsfFitterFunction", py::overload_cast, std::shared_ptr, - std::string, std::string, - std::size_t, bool, bool>( + std::string, std::string, std::size_t, bool, bool>( &ActsExamples::makeGsfFitterFunction), py::arg("trackingGeometry"), py::arg("magneticField"), py::arg("lowBetheHeitlerPath"), py::arg("highBetheHeitlerPath"), From b2a49a56991e1eff0825938f9a6496c6f70c2f35 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 27 Oct 2022 11:20:07 +0200 Subject: [PATCH 04/23] refactor default approximation --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 95 ++++++++++--------- .../TrackFitting/src/GsfFitterFunction.cpp | 4 +- docs/core/track_fitting.md | 2 + 3 files changed, 52 insertions(+), 49 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index e8d1ec441b2..0683b1e5f08 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -267,52 +267,55 @@ class AtlasBetheHeitlerApprox { } }; -/// These data are from ATLAS and allow using the GSF without loading files. -/// However, this might not be the optimal parameterization. These data come -/// this file in Athena: -/// Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par -/// These data must be transformed, so construct the AtlasBetheHeitlerApprox -/// with transforms = true -// clang-format off -constexpr static AtlasBetheHeitlerApprox<6, 5>::Data bh_cdf_cmps6_order5_data = {{ - // Component #1 - { - {{3.74397e+004,-1.95241e+004, 3.51047e+003,-2.54377e+002, 1.81080e+001,-3.57643e+000}}, - {{3.56728e+004,-1.78603e+004, 2.81521e+003,-8.93555e+001,-1.14015e+001, 2.55769e-001}}, - {{3.73938e+004,-1.92800e+004, 3.21580e+003,-1.46203e+002,-5.65392e+000,-2.78008e+000}} - }, - // Component #2 - { - {{-4.14035e+004, 2.31883e+004,-4.37145e+003, 2.44289e+002, 1.13098e+001,-3.21230e+000}}, - {{-2.06936e+003, 2.65334e+003,-1.01413e+003, 1.78338e+002,-1.85556e+001, 1.91430e+000}}, - {{-5.19068e+004, 2.55327e+004,-4.22147e+003, 1.90227e+002, 9.34602e+000,-4.80961e+000}} - }, - // Component #3 - { - {{2.52200e+003,-4.86348e+003, 2.11942e+003,-3.84534e+002, 2.94503e+001,-2.83310e+000}}, - {{1.80405e+003,-1.93347e+003, 6.27196e+002,-4.32429e+001,-1.43533e+001, 3.58782e+000}}, - {{-4.61617e+004, 1.78221e+004,-1.95746e+003,-8.80646e+001, 3.43153e+001,-7.57830e+000}} - }, - // Component #4 - { - {{4.94537e+003,-2.08737e+003, 1.78089e+002, 2.29879e+001,-5.52783e+000,-1.86800e+000}}, - {{4.60220e+003,-1.62269e+003,-1.57552e+002, 2.01796e+002,-5.01636e+001, 6.47438e+000}}, - {{-9.50373e+004, 4.05517e+004,-5.62596e+003, 4.58534e+001, 6.70479e+001,-1.22430e+001}} - }, - // Component #5 - { - {{-1.04129e+003, 1.15222e+002,-2.70356e+001, 3.18611e+001,-7.78800e+000,-1.50242e+000}}, - {{-2.71361e+004, 2.00625e+004,-6.19444e+003, 1.10061e+003,-1.29354e+002, 1.08289e+001}}, - {{3.15252e+004,-3.31508e+004, 1.20371e+004,-2.23822e+003, 2.44396e+002,-2.09130e+001}} - }, - // Component #6 - { - {{1.27751e+004,-6.79813e+003, 1.24650e+003,-8.20622e+001,-2.33476e+000, 2.46459e-001}}, - {{3.64336e+005,-2.08457e+005, 4.33028e+004,-3.67825e+003, 4.22914e+001, 1.42701e+001}}, - {{-1.79298e+006, 1.01843e+006,-2.10037e+005, 1.82222e+004,-4.33573e+002,-2.72725e+001}} - }, -}}; -// clang-format on +/// Creates a @ref AtlasBetheHeitlerApprox object based on a ATLAS +/// configuration, that are stored as static data in the source code. +/// This may not be an optimal configuration, but should allow running +/// the GSF without the need to load files +auto makeDefaultBetheHeitlerApprox() { + // Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par + // clang-format off + constexpr static AtlasBetheHeitlerApprox<6, 5>::Data cdf_cmps6_order5_data = {{ + // Component #1 + { + {{3.74397e+004,-1.95241e+004, 3.51047e+003,-2.54377e+002, 1.81080e+001,-3.57643e+000}}, + {{3.56728e+004,-1.78603e+004, 2.81521e+003,-8.93555e+001,-1.14015e+001, 2.55769e-001}}, + {{3.73938e+004,-1.92800e+004, 3.21580e+003,-1.46203e+002,-5.65392e+000,-2.78008e+000}} + }, + // Component #2 + { + {{-4.14035e+004, 2.31883e+004,-4.37145e+003, 2.44289e+002, 1.13098e+001,-3.21230e+000}}, + {{-2.06936e+003, 2.65334e+003,-1.01413e+003, 1.78338e+002,-1.85556e+001, 1.91430e+000}}, + {{-5.19068e+004, 2.55327e+004,-4.22147e+003, 1.90227e+002, 9.34602e+000,-4.80961e+000}} + }, + // Component #3 + { + {{2.52200e+003,-4.86348e+003, 2.11942e+003,-3.84534e+002, 2.94503e+001,-2.83310e+000}}, + {{1.80405e+003,-1.93347e+003, 6.27196e+002,-4.32429e+001,-1.43533e+001, 3.58782e+000}}, + {{-4.61617e+004, 1.78221e+004,-1.95746e+003,-8.80646e+001, 3.43153e+001,-7.57830e+000}} + }, + // Component #4 + { + {{4.94537e+003,-2.08737e+003, 1.78089e+002, 2.29879e+001,-5.52783e+000,-1.86800e+000}}, + {{4.60220e+003,-1.62269e+003,-1.57552e+002, 2.01796e+002,-5.01636e+001, 6.47438e+000}}, + {{-9.50373e+004, 4.05517e+004,-5.62596e+003, 4.58534e+001, 6.70479e+001,-1.22430e+001}} + }, + // Component #5 + { + {{-1.04129e+003, 1.15222e+002,-2.70356e+001, 3.18611e+001,-7.78800e+000,-1.50242e+000}}, + {{-2.71361e+004, 2.00625e+004,-6.19444e+003, 1.10061e+003,-1.29354e+002, 1.08289e+001}}, + {{3.15252e+004,-3.31508e+004, 1.20371e+004,-2.23822e+003, 2.44396e+002,-2.09130e+001}} + }, + // Component #6 + { + {{1.27751e+004,-6.79813e+003, 1.24650e+003,-8.20622e+001,-2.33476e+000, 2.46459e-001}}, + {{3.64336e+005,-2.08457e+005, 4.33028e+004,-3.67825e+003, 4.22914e+001, 1.42701e+001}}, + {{-1.79298e+006, 1.01843e+006,-2.10037e+005, 1.82222e+004,-4.33573e+002,-2.72725e+001}} + }, + }}; + // clang-format on + + return AtlasBetheHeitlerApprox<6,5>(cdf_cmps6_order5_data, cdf_cmps6_order5_data, true, true); +} } // namespace Experimental } // namespace Acts diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index 2c544ccaca6..99860719e5b 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -118,9 +118,7 @@ ActsExamples::makeGsfFitterFunction( const auto bhapp = [&]() { if (lowBetheHeitlerPath.empty() && highBetheHeitlerPath.empty()) { - return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>( - Acts::Experimental::bh_cdf_cmps6_order5_data, - Acts::Experimental::bh_cdf_cmps6_order5_data, true, true); + return Acts::Experimental::makeDefaultBetheHeitlerApprox(); } else if (std::filesystem::exists(lowBetheHeitlerPath) && std::filesystem::exists(highBetheHeitlerPath)) { return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>::loadFromFile( diff --git a/docs/core/track_fitting.md b/docs/core/track_fitting.md index 9c30c8e54ff..3cd1a72fd6c 100644 --- a/docs/core/track_fitting.md +++ b/docs/core/track_fitting.md @@ -75,6 +75,8 @@ This approximation of the Bethe-Heitler distribution is described in {class}`Act For small $x/x_0$ the {class}`Acts::detail::BetheHeitlerApprox` only returns a one-component mixture or no change at all. When loading a custom parametrization, it is possible to specify different parameterizations for high and for low $x/x_0$. The thresholds are currently not configurable. +A default parameterization can be created without files by {func}`Acts::Experimental::makeDefaultBetheHeitlerApprox`. + ### Further reading * *Thomas Atkinson*, Electron reconstruction with the ATLAS inner detector, 2006, see [here](https://cds.cern.ch/record/1448253) * *R Frühwirth*, Track fitting with non-Gaussian noise, 1997, see [here](https://doi.org/10.1016/S0010-4655(96)00155-5) From 1978f8433971598a8be03fc0315806d9ed2efb1c Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 27 Oct 2022 11:20:37 +0200 Subject: [PATCH 05/23] clang-format --- Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 0683b1e5f08..e7e1e15bfa4 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -314,7 +314,8 @@ auto makeDefaultBetheHeitlerApprox() { }}; // clang-format on - return AtlasBetheHeitlerApprox<6,5>(cdf_cmps6_order5_data, cdf_cmps6_order5_data, true, true); + return AtlasBetheHeitlerApprox<6, 5>(cdf_cmps6_order5_data, + cdf_cmps6_order5_data, true, true); } } // namespace Experimental From 7231618c9af03a053961d18d9da673625ef56d79 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 13:03:00 +0200 Subject: [PATCH 06/23] expose bethe heitler to python --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 4 ++-- .../TrackFitting/GsfFitterFunction.hpp | 13 ++++++------ .../TrackFitting/src/GsfFitterFunction.cpp | 21 ++++--------------- .../python/acts/examples/reconstruction.py | 1 + Examples/Python/src/TrackFitting.cpp | 16 ++++++++++---- 5 files changed, 25 insertions(+), 30 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index e7e1e15bfa4..15a3f2d89c5 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -213,8 +213,8 @@ class AtlasBetheHeitlerApprox { /// the parameterization for low x/x0 /// @param high_parameters_path Path to the foo.par file that stores /// the parameterization for high x/x0 - static auto loadFromFile(const std::string &low_parameters_path, - const std::string &high_parameters_path) { + static auto loadFromFiles(const std::string &low_parameters_path, + const std::string &high_parameters_path) { auto read_file = [](const std::string &filepath) { std::ifstream file(filepath); diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp index 8d071b6fc19..274a6fa4ab8 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp @@ -8,6 +8,7 @@ #pragma once +#include "Acts/TrackFitting/BetheHeitlerApprox.hpp" #include "ActsExamples/TrackFitting/TrackFittingAlgorithm.hpp" namespace ActsExamples { @@ -16,12 +17,10 @@ namespace ActsExamples { /// /// @param trackingGeometry the trackingGeometry for the propagator /// @param magneticField the magnetic field for the propagator -/// @param lowBetheHeitlerPath path to a parameterization of the bethe heitler -/// distribution suitable for low thicknesses (x/x0 < 0.1). Pass empty string -/// to load a default. -/// @param highBetheHeitlerPath path to a parameterization of the bethe heitler -/// distribution suitable for higher thicknesses (0.1 < x/x0 < 0.2). Pass empty -/// string to load a default. +/// @param betheHeitlerApprox The object that encapsulates the approximation. +/// For the examples framework this is fixed to a ATLAS-like approximation with +/// 6 components interpolated by a 5th order polynomial, but in general this can +/// be customized. /// @param maxComponents number of maximum components in the track state /// @param abortOnError wether to call std::abort on an error /// @param disableAllMaterialHandling run the GSF like a KF (no energy loss, @@ -30,7 +29,7 @@ std::shared_ptr makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, - std::string lowBetheHeitlerPath, std::string highBetheHeitlerPath, + Acts::Experimental::AtlasBetheHeitlerApprox<6, 5> betheHeitlerApprox, std::size_t maxComponents, bool abortOnError, bool disableAllMaterialHandling); diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index 99860719e5b..8404c8f970d 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -111,25 +111,11 @@ std::shared_ptr ActsExamples::makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, - std::string lowBetheHeitlerPath, std::string highBetheHeitlerPath, + Acts::Experimental::AtlasBetheHeitlerApprox<6, 5> betheHeitlerApprox, std::size_t maxComponents, bool abortOnError, bool disableAllMaterialHandling) { MultiStepper stepper(std::move(magneticField)); - const auto bhapp = [&]() { - if (lowBetheHeitlerPath.empty() && highBetheHeitlerPath.empty()) { - return Acts::Experimental::makeDefaultBetheHeitlerApprox(); - } else if (std::filesystem::exists(lowBetheHeitlerPath) && - std::filesystem::exists(highBetheHeitlerPath)) { - return Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>::loadFromFile( - lowBetheHeitlerPath, highBetheHeitlerPath); - } else { - throw std::invalid_argument( - "Paths to bethe heitler parameterization do not exist. Pass empty " - "strings to load a default parameterization"); - } - }(); - // Standard fitter Acts::Navigator::Config cfg{trackingGeometry}; cfg.resolvePassive = false; @@ -137,12 +123,13 @@ ActsExamples::makeGsfFitterFunction( cfg.resolveSensitive = true; Acts::Navigator navigator(cfg); Propagator propagator(std::move(stepper), std::move(navigator)); - Fitter trackFitter(std::move(propagator), BHApprox(bhapp)); + Fitter trackFitter(std::move(propagator), BHApprox(betheHeitlerApprox)); // Direct fitter Acts::DirectNavigator directNavigator; DirectPropagator directPropagator(stepper, directNavigator); - DirectFitter directTrackFitter(std::move(directPropagator), BHApprox(bhapp)); + DirectFitter directTrackFitter(std::move(directPropagator), + BHApprox(betheHeitlerApprox)); // build the fitter functions. owns the fitter object. auto fitterFunction = std::make_shared( diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index e986a553580..a2038ce36e0 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -637,6 +637,7 @@ def addTruthTrackingGsf( customLogLevel = acts.examples.defaultLogging(s, logLevel) gsfOptions = { + "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(), "maxComponents": 12, "abortOnError": False, "disableAllMaterialHandling": False, diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index 444c035dbff..f9c0e3be1b1 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -72,16 +72,24 @@ void addTrackFitting(Context& ctx) { py::arg("reverseFilteringMomThreshold"), py::arg("freeToBoundCorrection")); + using BetheHeitlerApprox = + Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>; + py::class_(mex, "AtlasBetheHeitlerApprox") + .def_static("loadFromFiles", &BetheHeitlerApprox::loadFromFiles, + py::arg("lowParametersPath"), py::arg("lowParametersPath")) + .def_static("makeDefault", []() { + return Acts::Experimental::makeDefaultBetheHeitlerApprox(); + }); + mex.def( "makeGsfFitterFunction", py::overload_cast, std::shared_ptr, - std::string, std::string, std::size_t, bool, bool>( + BetheHeitlerApprox, std::size_t, bool, bool>( &ActsExamples::makeGsfFitterFunction), py::arg("trackingGeometry"), py::arg("magneticField"), - py::arg("lowBetheHeitlerPath"), py::arg("highBetheHeitlerPath"), - py::arg("maxComponents"), py::arg("abortOnError"), - py::arg("disableAllMaterialHandling")); + py::arg("betheHeitlerApprox"), py::arg("maxComponents"), + py::arg("abortOnError"), py::arg("disableAllMaterialHandling")); } { From 00a948d1f27a9b535ec728f3dac16b2d7ace93e6 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 13:13:15 +0200 Subject: [PATCH 07/23] change to boost::static_vector and allow non-transforming files --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 21 ++++++++++++------- Examples/Python/src/TrackFitting.cpp | 2 +- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 15a3f2d89c5..9115795e038 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -16,15 +16,16 @@ #include #include -#include -#include +#include namespace Acts { namespace detail { struct GaussianComponent { - ActsScalar weight, mean, var; + ActsScalar weight = 0.0; + ActsScalar mean = 0.0; + ActsScalar var = 0.0; }; /// Transform a gaussian component to a space where all values are defined from @@ -97,7 +98,9 @@ class AtlasBetheHeitlerApprox { public: struct PolyData { - std::array weightCoeffs, meanCoeffs, varCoeffs; + std::array weightCoeffs; + std::array meanCoeffs; + std::array varCoeffs; }; using Data = std::array; @@ -139,6 +142,8 @@ class AtlasBetheHeitlerApprox { /// /// @param x The input in terms of x/x0 (pathlength in terms of radiation length) auto mixture(ActsScalar x) const { + using Array = + boost::container::static_vector; // Build a polynom auto poly = [](ActsScalar xx, const std::array &coeffs) { @@ -146,14 +151,14 @@ class AtlasBetheHeitlerApprox { for (const auto c : coeffs) { sum = xx * sum + c; } - throw_assert(std::isfinite(sum), "polynom result not finite"); + assert((std::isfinite(sum) && "polynom result not finite")); return sum; }; // Lambda which builds the components auto make_mixture = [&](const Data &data, double xx, bool transform) { // Value initialization should garanuee that all is initialized to zero - std::array ret{}; + Array ret(NComponents); ActsScalar weight_sum = 0; for (int i = 0; i < NComponents; ++i) { // These transformations must be applied to the data according to ATHENA @@ -180,7 +185,7 @@ class AtlasBetheHeitlerApprox { // Return no change if (x < noChangeLimit) { - std::array ret{}; + Array ret(1); ret[0].weight = 1.0; ret[0].mean = 1.0; // p_initial = p_final @@ -190,7 +195,7 @@ class AtlasBetheHeitlerApprox { } // Return single gaussian approximation if (x < singleGaussianLimit) { - std::array ret{}; + Array ret(1); ret[0] = BetheHeitlerApproxSingleCmp::mixture(x)[0]; return ret; } diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index f9c0e3be1b1..2b59ae91bda 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -77,7 +77,7 @@ void addTrackFitting(Context& ctx) { py::class_(mex, "AtlasBetheHeitlerApprox") .def_static("loadFromFiles", &BetheHeitlerApprox::loadFromFiles, py::arg("lowParametersPath"), py::arg("lowParametersPath")) - .def_static("makeDefault", []() { + .def_static("makeDefault ", []() { return Acts::Experimental::makeDefaultBetheHeitlerApprox(); }); From 100785b460169ee5046fccbdb7e083ca6c581406 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 13:15:30 +0200 Subject: [PATCH 08/23] remove dead code and minor refactorings --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 9115795e038..9356de735e0 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -79,10 +79,6 @@ struct BetheHeitlerApproxSingleCmp { ret[0].mean = std::pow(2, -c); ret[0].var = std::pow(3, -c) - std::pow(4, -c); - // ret[0].mean = std::exp(-1. * x); - // ret[0].var = - // std::exp(-1. * x * std::log(3.) / std::log(2.)) - std::exp(-2. * x); - return ret; } }; @@ -204,11 +200,9 @@ class AtlasBetheHeitlerApprox { return make_mixture(m_low_data, x, m_low_transform); } // Return a component representation for higher x0 - else { - // Cap the x because beyond the parameterization goes wild - const auto high_x = std::min(higherLimit, x); - return make_mixture(m_high_data, high_x, m_high_transform); - } + // Cap the x because beyond the parameterization goes wild + const auto high_x = std::min(higherLimit, x); + return make_mixture(m_high_data, high_x, m_high_transform); } /// Loads a parameterization from a file according to the Atlas file @@ -242,11 +236,6 @@ class AtlasBetheHeitlerApprox { filepath + "'"); } - if (!transform_code) { - throw std::invalid_argument("Transform-code is required in '" + - filepath + "'"); - } - Data data; for (auto &cmp : data) { From 5c78e5ebe84ea0e387efcd32e051170d89cb4f6c Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 13:36:08 +0200 Subject: [PATCH 09/23] add warning when encountering invalid x/x0 and improve docs --- Core/include/Acts/TrackFitting/detail/GsfActor.hpp | 5 +++++ docs/core/track_fitting.md | 13 ++++++++----- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 5e35d615a45..036170c42fb 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -348,6 +348,11 @@ struct GsfActor { old_bound.unitDirection()); slab.scaleThickness(pathCorrection); + // Emit a warning if the approximation is not valid for this x/x0 + if( not m_cfg.bethe_heitler_approx->validXOverX0(slab.thicknessInX0()) ) { + ACTS_WARNING("Bethe-Heitler approximation encountered invalid value for x/x0 " << slab.thicknessInX0()); + } + // Get the mixture const auto mixture = m_cfg.bethe_heitler_approx->mixture(slab.thicknessInX0()); diff --git a/docs/core/track_fitting.md b/docs/core/track_fitting.md index 3cd1a72fd6c..64907524a99 100644 --- a/docs/core/track_fitting.md +++ b/docs/core/track_fitting.md @@ -69,17 +69,20 @@ A GSF example can be found in the Acts Examples Framework [here](https://github. ### Customizing the Bethe-Heitler approximation -To be able to evaluate the approximation of the Bethe-Heitler distribution for different materials and thicknesses, the individual gaussian components (weight, mean, variance of the ratio $E_f/E_i$) are parameterized as polynomials in $x/x_0$. The default parametrization uses 6 components and 5th order polynomials. +The GSF needs an approximation of the Bethe-Heitler distribution as a Gaussian mixture on each material interaction (see above). This task is delegated to a separate class, that can be provided by a template parameter to {class}`Acts::Experimental::GaussianSumFitter`, so in principle it can be implemented in different ways. -This approximation of the Bethe-Heitler distribution is described in {class}`Acts::detail::BetheHeitlerApprox`. The class is templated on the number of components and the degree of the polynomial, and is designed to be used with the [parameterization files from ATLAS](https://gitlab.cern.ch/atlas/athena/-/tree/master/Tracking/TrkFitter/TrkGaussianSumFilter/Data). However, in principle the GSF could be constructed with custom classes with the same interface as {class}`Acts::detail::BetheHeitlerApprox`. +However, ACTS ships with the class {class}`Acts::Experimental::AtlasBetheHeitlerApprox` that implements the ATLAS strategy for this task: To be able to evaluate the approximation of the Bethe-Heitler distribution for different materials and thicknesses, the individual Gaussian components (weight, mean, variance of the ratio $E_f/E_i$) are parametrised as polynomials in $x/x_0$. This class can load files in the ATLAS format that can be found [here](https://gitlab.cern.ch/atlas/athena/-/tree/master/Tracking/TrkFitter/TrkGaussianSumFilter/Data). A default parameterization can be created with {func}`Acts::Experimental::makeDefaultBetheHeitlerApprox`. -For small $x/x_0$ the {class}`Acts::detail::BetheHeitlerApprox` only returns a one-component mixture or no change at all. When loading a custom parametrization, it is possible to specify different parameterizations for high and for low $x/x_0$. The thresholds are currently not configurable. - -A default parameterization can be created without files by {func}`Acts::Experimental::makeDefaultBetheHeitlerApprox`. +The {class}`Acts::Experimental::AtlasBetheHeitlerApprox` is constructed with two parameterizations, allowing to use different parameterizations for different $x/x_0$. In particular, it has this behaviour: +* $x/x_0 < 0.0001$: Return no change +* $x/x_0 < 0.002$: Return a single gaussian approximation +* $x/x_0 < 0.1$: Return the approximation for low $x/x_0$. +* $x/x_0 \geq 0.1$: Return the approximation for high $x/x_0$. The maximum possible value is $x/x_0 = 0.2$, for higher values it is clipped to 0.2 and the GSF emits a warning. ### Further reading * *Thomas Atkinson*, Electron reconstruction with the ATLAS inner detector, 2006, see [here](https://cds.cern.ch/record/1448253) * *R Frühwirth*, Track fitting with non-Gaussian noise, 1997, see [here](https://doi.org/10.1016/S0010-4655(96)00155-5) +* *R Frühwirth*, A Gaussian-mixture approximation of the Bethe–Heitler model of electron energy loss by bremsstrahlung, 2003, see [here](https://doi.org/10.1016/S0010-4655(03)00292-3) (gx2f_core)= ## Global Chi-Square Fitter (GX2F) [wip] From 53f393d2d7ae2ce1ad2827986e1d0eeac2455d5c Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 13:37:24 +0200 Subject: [PATCH 10/23] clang-format --- Core/include/Acts/TrackFitting/detail/GsfActor.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 036170c42fb..93ff7a3011d 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -349,8 +349,10 @@ struct GsfActor { slab.scaleThickness(pathCorrection); // Emit a warning if the approximation is not valid for this x/x0 - if( not m_cfg.bethe_heitler_approx->validXOverX0(slab.thicknessInX0()) ) { - ACTS_WARNING("Bethe-Heitler approximation encountered invalid value for x/x0 " << slab.thicknessInX0()); + if (not m_cfg.bethe_heitler_approx->validXOverX0(slab.thicknessInX0())) { + ACTS_WARNING( + "Bethe-Heitler approximation encountered invalid value for x/x0 " + << slab.thicknessInX0()); } // Get the mixture From 630e60bcf80cc15eb84be657a0aa3e5b5b190d29 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 14:04:22 +0200 Subject: [PATCH 11/23] changes & clang-format --- .../TrackFitting/GsfFitterFunction.hpp | 12 ++++++------ .../TrackFitting/src/GsfFitterFunction.cpp | 15 +++++++-------- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp index 274a6fa4ab8..4b043ea9fca 100644 --- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp +++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/GsfFitterFunction.hpp @@ -13,14 +13,15 @@ namespace ActsExamples { +/// This type is used in the Examples framework for the Bethe-Heitler +/// approximation +using BetheHeitlerApprox = Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>; + /// Makes a fitter function object for the GSF /// /// @param trackingGeometry the trackingGeometry for the propagator /// @param magneticField the magnetic field for the propagator /// @param betheHeitlerApprox The object that encapsulates the approximation. -/// For the examples framework this is fixed to a ATLAS-like approximation with -/// 6 components interpolated by a 5th order polynomial, but in general this can -/// be customized. /// @param maxComponents number of maximum components in the track state /// @param abortOnError wether to call std::abort on an error /// @param disableAllMaterialHandling run the GSF like a KF (no energy loss, @@ -29,8 +30,7 @@ std::shared_ptr makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, - Acts::Experimental::AtlasBetheHeitlerApprox<6, 5> betheHeitlerApprox, - std::size_t maxComponents, bool abortOnError, - bool disableAllMaterialHandling); + BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, + bool abortOnError, bool disableAllMaterialHandling); } // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp index 8404c8f970d..e12e09a4c50 100644 --- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp +++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp @@ -32,12 +32,11 @@ using MultiStepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator; using DirectPropagator = Acts::Propagator; -using BHApprox = Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>; using Fitter = - Acts::Experimental::GaussianSumFitter; using DirectFitter = - Acts::Experimental::GaussianSumFitter; struct GsfFitterFunctionImpl @@ -111,9 +110,8 @@ std::shared_ptr ActsExamples::makeGsfFitterFunction( std::shared_ptr trackingGeometry, std::shared_ptr magneticField, - Acts::Experimental::AtlasBetheHeitlerApprox<6, 5> betheHeitlerApprox, - std::size_t maxComponents, bool abortOnError, - bool disableAllMaterialHandling) { + BetheHeitlerApprox betheHeitlerApprox, std::size_t maxComponents, + bool abortOnError, bool disableAllMaterialHandling) { MultiStepper stepper(std::move(magneticField)); // Standard fitter @@ -123,13 +121,14 @@ ActsExamples::makeGsfFitterFunction( cfg.resolveSensitive = true; Acts::Navigator navigator(cfg); Propagator propagator(std::move(stepper), std::move(navigator)); - Fitter trackFitter(std::move(propagator), BHApprox(betheHeitlerApprox)); + Fitter trackFitter(std::move(propagator), + BetheHeitlerApprox(betheHeitlerApprox)); // Direct fitter Acts::DirectNavigator directNavigator; DirectPropagator directPropagator(stepper, directNavigator); DirectFitter directTrackFitter(std::move(directPropagator), - BHApprox(betheHeitlerApprox)); + BetheHeitlerApprox(betheHeitlerApprox)); // build the fitter functions. owns the fitter object. auto fitterFunction = std::make_shared( From 1c05ce32267014ef78924e8d40e32a12fe9669d8 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 14:05:51 +0200 Subject: [PATCH 12/23] use only one typedef --- Examples/Python/src/TrackFitting.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index 2b59ae91bda..52cfea55ef9 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -72,10 +72,9 @@ void addTrackFitting(Context& ctx) { py::arg("reverseFilteringMomThreshold"), py::arg("freeToBoundCorrection")); - using BetheHeitlerApprox = - Acts::Experimental::AtlasBetheHeitlerApprox<6, 5>; - py::class_(mex, "AtlasBetheHeitlerApprox") - .def_static("loadFromFiles", &BetheHeitlerApprox::loadFromFiles, + py::class_(mex, "AtlasBetheHeitlerApprox") + .def_static("loadFromFiles", + &ActsExamples::BetheHeitlerApprox::loadFromFiles, py::arg("lowParametersPath"), py::arg("lowParametersPath")) .def_static("makeDefault ", []() { return Acts::Experimental::makeDefaultBetheHeitlerApprox(); From 63cc3aa0d0131a47554489ce690520cfdf28429b Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Fri, 28 Oct 2022 14:19:05 +0200 Subject: [PATCH 13/23] fix unit test --- Tests/UnitTests/Core/TrackFitting/GsfTests.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp b/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp index 492d6c6ec8a..c128bcd0012 100644 --- a/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp @@ -50,8 +50,10 @@ using Propagator = Acts::Propagator; auto gsfZeroPropagator = makeConstantFieldPropagator(tester.geometry, 0_T); -const GaussianSumFitter gsfZero( - std::move(gsfZeroPropagator)); +auto betheHeitlerApprox = Acts::Experimental::makeDefaultBetheHeitlerApprox(); +const GaussianSumFitter + gsfZero(std::move(gsfZeroPropagator), std::move(betheHeitlerApprox)); std::default_random_engine rng(42); From 49f5185a44432678d3b0ff6ce282577b5a8f1ff2 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Mon, 31 Oct 2022 09:31:18 +0100 Subject: [PATCH 14/23] fix python function binding typo --- Examples/Python/src/TrackFitting.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp index 52cfea55ef9..07d99ca65e7 100644 --- a/Examples/Python/src/TrackFitting.cpp +++ b/Examples/Python/src/TrackFitting.cpp @@ -76,7 +76,7 @@ void addTrackFitting(Context& ctx) { .def_static("loadFromFiles", &ActsExamples::BetheHeitlerApprox::loadFromFiles, py::arg("lowParametersPath"), py::arg("lowParametersPath")) - .def_static("makeDefault ", []() { + .def_static("makeDefault", []() { return Acts::Experimental::makeDefaultBetheHeitlerApprox(); }); From cfa3b8ece636809c3b51a02ade1f0f29fb225270 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Tue, 1 Nov 2022 15:01:25 +0100 Subject: [PATCH 15/23] fix docs --- .../Acts/TrackFitting/BetheHeitlerApprox.hpp | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 9356de735e0..0342482cfa7 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -57,19 +57,24 @@ auto inverseTransformComponent(double transformed_weight, namespace Experimental { -/// This class approximates the Bethe-Heitler with only one component. The +/// This class approximates the Bethe-Heitler with only one component. This is +/// mainly inside @ref AtlasBetheHeitlerApprox, but can also be used as the +/// only component approximation (then probably for debugging) struct BetheHeitlerApproxSingleCmp { /// Returns the number of components the returned mixture will have constexpr auto numComponents() const { return 1; } - /// Checks if an input is valid for the parameterization. Since this is for - /// debugging, it always returns false - /// - /// @param x input in terms of x/x0 - constexpr bool validXOverX0(ActsScalar) const { return false; } + /// Checks if an input is valid for the parameterization. The threshold for + /// x/x0 is 0.002 and orientates on the values used in ATLAS + constexpr bool validXOverX0(ActsScalar x) const { + return x < 0.002; + ; + } /// Returns array with length 1 containing a 1-component-representation of the /// Bethe-Heitler-Distribution + /// + /// @param x pathlength in terms of the radiation length static auto mixture(const ActsScalar x) { std::array ret{}; @@ -130,13 +135,13 @@ class AtlasBetheHeitlerApprox { /// Checks if an input is valid for the parameterization /// - /// @param x input in terms of x/x0 + /// @param x pathlength in terms of the radiation length constexpr bool validXOverX0(ActsScalar x) const { return x < higherLimit; } /// Generates the mixture from the polynomials and reweights them, so /// that the sum of all weights is 1 /// - /// @param x The input in terms of x/x0 (pathlength in terms of radiation length) + /// @param x pathlength in terms of the radiation length auto mixture(ActsScalar x) const { using Array = boost::container::static_vector; From 8f8df582592a2f4fff4bca7ec1b0d6e184a523b1 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Tue, 1 Nov 2022 15:03:03 +0100 Subject: [PATCH 16/23] improve invalid warning --- Core/include/Acts/TrackFitting/detail/GsfActor.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 93ff7a3011d..58594a7935f 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -351,8 +351,8 @@ struct GsfActor { // Emit a warning if the approximation is not valid for this x/x0 if (not m_cfg.bethe_heitler_approx->validXOverX0(slab.thicknessInX0())) { ACTS_WARNING( - "Bethe-Heitler approximation encountered invalid value for x/x0 " - << slab.thicknessInX0()); + "Bethe-Heitler approximation encountered invalid value for x/x0=" + << slab.thicknessInX0() << " at surface " << surface.geometryId()); } // Get the mixture From 158ded00dca8b73bc912b0dd7b5fbb3eef138cc9 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Thu, 3 Nov 2022 10:24:23 +0100 Subject: [PATCH 17/23] fix docs --- Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 0342482cfa7..4f0d4159515 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -118,11 +118,14 @@ class AtlasBetheHeitlerApprox { bool m_high_transform; public: - /// Construct the Bethe-Heitler approximation description + /// Construct the Bethe-Heitler approximation description. Additional to the + /// coefficients of the polynomials, the information wether these values need + /// to be transformed beforehand must be given (see ATLAS code). /// /// @param low_data data for the lower x/x0 range /// @param high_data data for the higher x/x0 range - /// @param transform wether the data need to be transformed (see Atlas code) + /// @param low_transform wether the low data need to be transformed + /// @param high_transform wether the high data need to be transformed constexpr AtlasBetheHeitlerApprox(const Data &low_data, const Data &high_data, bool low_transform, bool high_transform) : m_low_data(low_data), From f07036ecdee92b7d2637cfbfacf07f05dcc638c8 Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Thu, 3 Nov 2022 17:08:47 +0100 Subject: [PATCH 18/23] Update root_file_hashes.txt --- Examples/Python/tests/root_file_hashes.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index b057dbbf655..6b87dc2364b 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -75,7 +75,7 @@ test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: c7b3ff9d8d3c19 test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: 76a990d595b6e097da2bed447783bd63044956e5649a5dd6fd7a6a3434786877 test_truth_tracking_gsf[generic]__trackstates_gsf.root: c9f3d30d13ee8eac8c69a01600ea89a8bd80cf89154ecc69a48d60a6afd26405 -test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 14fd6b28b8ac3722af4fde22fedd48ec9bfed6075109a8e6f327d551acaf44a8 +test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 7199b3912e089b93cc010eee74356d4f23edaae2eff7c586b836911e65ae6e9a test_truth_tracking_gsf[odd]__trackstates_gsf.root: 97ebc8f5ab5d5242596eaca2dfbed864a65580cc18c02cad634c6357cf449ab7 test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 928277c50d249abd548185f8536e1614e9ab958e0e7a9521ef0966acdcccb9d2 From 50a58d7aa0a3b7a95e6e1c77e9bb7125116d0849 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Mon, 7 Nov 2022 09:24:46 +0100 Subject: [PATCH 19/23] make clang-tidy happy --- Tests/UnitTests/Core/TrackFitting/GsfTests.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp b/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp index c128bcd0012..0fde60392dd 100644 --- a/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp @@ -47,13 +47,11 @@ const auto logger = getDefaultLogger("GSF", Logging::INFO); using Stepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator; +using GSF = GaussianSumFitter; -auto gsfZeroPropagator = - makeConstantFieldPropagator(tester.geometry, 0_T); -auto betheHeitlerApprox = Acts::Experimental::makeDefaultBetheHeitlerApprox(); -const GaussianSumFitter - gsfZero(std::move(gsfZeroPropagator), std::move(betheHeitlerApprox)); +const GSF gsfZero(makeConstantFieldPropagator(tester.geometry, 0_T), + makeDefaultBetheHeitlerApprox()); std::default_random_engine rng(42); From 83cf551a93aedeb37029742956f9f19b98b69f98 Mon Sep 17 00:00:00 2001 From: Benjamin Huth Date: Mon, 7 Nov 2022 09:39:35 +0100 Subject: [PATCH 20/23] fix --- Tests/UnitTests/Core/TrackFitting/GsfTests.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp b/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp index 0fde60392dd..7a5f468b262 100644 --- a/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp +++ b/Tests/UnitTests/Core/TrackFitting/GsfTests.cpp @@ -47,8 +47,9 @@ const auto logger = getDefaultLogger("GSF", Logging::INFO); using Stepper = Acts::MultiEigenStepperLoop<>; using Propagator = Acts::Propagator; -using GSF = GaussianSumFitter; +using BetheHeitlerApprox = AtlasBetheHeitlerApprox<6, 5>; +using GSF = + GaussianSumFitter; const GSF gsfZero(makeConstantFieldPropagator(tester.geometry, 0_T), makeDefaultBetheHeitlerApprox()); From 1b86c85cb13a2095dba7d2084174f5600004dd84 Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Mon, 7 Nov 2022 09:48:24 +0100 Subject: [PATCH 21/23] Update Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp Co-authored-by: Alexander J. Pfleger <70842573+AJPfleger@users.noreply.github.com> --- Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 4f0d4159515..f54273d8d4f 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -119,7 +119,7 @@ class AtlasBetheHeitlerApprox { public: /// Construct the Bethe-Heitler approximation description. Additional to the - /// coefficients of the polynomials, the information wether these values need + /// coefficients of the polynomials, the information whether these values need /// to be transformed beforehand must be given (see ATLAS code). /// /// @param low_data data for the lower x/x0 range From bbb157a3adebf445ddb3464bbe38fb833baf4493 Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Mon, 7 Nov 2022 09:53:25 +0100 Subject: [PATCH 22/23] Update BetheHeitlerApprox.hpp --- Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index f54273d8d4f..5e4c67338c3 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -89,7 +89,7 @@ struct BetheHeitlerApproxSingleCmp { }; /// This class approximates the Bethe-Heitler distribution as a gaussian -/// mixture. To enable an approximation for continous input variables, the +/// mixture. To enable an approximation for continuous input variables, the /// weights, means and variances are internally parametrized as a Nth order /// polynomial. template From 9ff7c34d2b470a72422a2fe1611bfca1f6cfb27a Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Mon, 7 Nov 2022 09:53:59 +0100 Subject: [PATCH 23/23] Update Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp Co-authored-by: Alexander J. Pfleger <70842573+AJPfleger@users.noreply.github.com> --- Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp index 5e4c67338c3..0be17d4e635 100644 --- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp +++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp @@ -269,9 +269,9 @@ class AtlasBetheHeitlerApprox { } }; -/// Creates a @ref AtlasBetheHeitlerApprox object based on a ATLAS +/// Creates a @ref AtlasBetheHeitlerApprox object based on an ATLAS /// configuration, that are stored as static data in the source code. -/// This may not be an optimal configuration, but should allow running +/// This may not be an optimal configuration, but should allow to run /// the GSF without the need to load files auto makeDefaultBetheHeitlerApprox() { // Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par