src/TrackFitting/BetheHeitlerApprox.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/TrackFitting/BetheHeitlerApprox.hpp"
0010
0011 #include <algorithm>
0012 #include <fstream>
0013 #include <stdexcept>
0014 #include <tuple>
0015
0016 namespace Acts {
0017
0018 AtlasBetheHeitlerApprox AtlasBetheHeitlerApprox::loadFromFiles(
0019     const std::string &low_parameters_path,
0020     const std::string &high_parameters_path, double lowLimit, double highLimit,
0021     bool clampToRange, double noChangeLimit, double singleGaussianLimit) {
0022   const auto read_file = [](const std::string &filepath) {
0023     std::ifstream file(filepath);
0024
0025     if (!file) {
0026       throw std::invalid_argument("Could not open '" + filepath + "'");
0027     }
0028
0029     std::size_t n_cmps = 0;
0030     std::size_t degree = 0;
0031     bool transform_code = false;
0032
0033     file >> n_cmps >> degree >> transform_code;
0034
0035     Data data;
0036
0037     for (auto &cmp : data) {
0038       for (auto &coeff : cmp.weightCoeffs) {
0039         file >> coeff;
0040       }
0041       for (auto &coeff : cmp.meanCoeffs) {
0042         file >> coeff;
0043       }
0044       for (auto &coeff : cmp.varCoeffs) {
0045         file >> coeff;
0046       }
0047     }
0048
0049     return std::make_tuple(data, transform_code);
0050   };
0051
0052   const auto [lowData, lowTransform] = read_file(low_parameters_path);
0053   const auto [highData, highTransform] = read_file(high_parameters_path);
0054
0055   return {lowData,   highData,     lowTransform,  highTransform,      lowLimit,
0056           highLimit, clampToRange, noChangeLimit, singleGaussianLimit};
0057 }
0058
0059 std::span<AtlasBetheHeitlerApprox::Component> AtlasBetheHeitlerApprox::mixture(
0060     double xOverX0, const std::span<Component> mixture) const {
0061   if (m_clampToRange) {
0062     xOverX0 = std::clamp(xOverX0, 0.0, m_highLimit);
0063   }
0064
0065   // Evaluate polynomial at x
0066   const auto poly = [](const double xx, const std::span<const double> coeffs) {
0067     double sum{0.};
0068     for (const auto c : coeffs) {
0069       sum = xx * sum + c;
0070     }
0071     return sum;
0072   };
0073
0074   // Lambda which builds the components
0075   const auto make_mixture = [&](const Data &data, double xx,
0076                                 bool transform) -> std::span<Component> {
0077     // Value initialization should garanuee that all is initialized to zero
0078     double weight_sum = 0;
0079     for (std::size_t i = 0; i < data.size(); ++i) {
0080       // These transformations must be applied to the data according to ATHENA
0081       // (TrkGaussianSumFilter/src/GsfCombinedMaterialEffects.cxx:79)
0082       if (transform) {
0083         mixture[i] = detail::inverseTransformComponent(
0084             {poly(xx, data[i].weightCoeffs), poly(xx, data[i].meanCoeffs),
0085              poly(xx, data[i].varCoeffs)});
0086       } else {
0087         mixture[i].weight = poly(xx, data[i].weightCoeffs);
0088         mixture[i].mean = poly(xx, data[i].meanCoeffs);
0089         mixture[i].var = poly(xx, data[i].varCoeffs);
0090       }
0091
0092       weight_sum += mixture[i].weight;
0093     }
0094
0095     for (std::size_t i = 0; i < data.size(); ++i) {
0096       mixture[i].weight /= weight_sum;
0097     }
0098
0099     return {mixture.data(), data.size()};
0100   };
0101
0102   // Return no change
0103   if (xOverX0 < m_noChangeLimit) {
0104     mixture[0].weight = 1.0;
0105     mixture[0].mean = 1.0;  // p_initial = p_final
0106     mixture[0].var = 0.0;
0107
0108     return {mixture.data(), 1};
0109   }
0110
0111   // Return single gaussian approximation
0112   if (xOverX0 < m_singleGaussianLimit) {
0113     BetheHeitlerApproxSingleCmp().mixture(xOverX0, mixture);
0114     return {mixture.data(), 1};
0115   }
0116
0117   // Return a component representation for lower x0
0118   if (xOverX0 < m_lowLimit) {
0119     return make_mixture(m_lowData, xOverX0, m_lowTransform);
0120   }
0121
0122   // Return a component representation for higher x0
0123   // Cap the x because beyond the parameterization goes wild
0124   const auto high_x = std::min(m_highLimit, xOverX0);
0125   return make_mixture(m_highData, high_x, m_highTransform);
0126 }
0127
0128 }  // namespace Acts
0129
0130 Acts::AtlasBetheHeitlerApprox Acts::makeDefaultBetheHeitlerApprox(
0131     bool clampToRange) {
0132   // Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par
0133   // clang-format off
0134   static AtlasBetheHeitlerApprox::Data cdf_cmps6_order5_data = {{
0135       // Component #1
0136       {
0137           {{3.74397e+004,-1.95241e+004, 3.51047e+003,-2.54377e+002, 1.81080e+001,-3.57643e+000}},
0138           {{3.56728e+004,-1.78603e+004, 2.81521e+003,-8.93555e+001,-1.14015e+001, 2.55769e-001}},
0139           {{3.73938e+004,-1.92800e+004, 3.21580e+003,-1.46203e+002,-5.65392e+000,-2.78008e+000}}
0140       },
0141       // Component #2
0142       {
0143           {{-4.14035e+004, 2.31883e+004,-4.37145e+003, 2.44289e+002, 1.13098e+001,-3.21230e+000}},
0144           {{-2.06936e+003, 2.65334e+003,-1.01413e+003, 1.78338e+002,-1.85556e+001, 1.91430e+000}},
0145           {{-5.19068e+004, 2.55327e+004,-4.22147e+003, 1.90227e+002, 9.34602e+000,-4.80961e+000}}
0146       },
0147       // Component #3
0148       {
0149           {{2.52200e+003,-4.86348e+003, 2.11942e+003,-3.84534e+002, 2.94503e+001,-2.83310e+000}},
0150           {{1.80405e+003,-1.93347e+003, 6.27196e+002,-4.32429e+001,-1.43533e+001, 3.58782e+000}},
0151           {{-4.61617e+004, 1.78221e+004,-1.95746e+003,-8.80646e+001, 3.43153e+001,-7.57830e+000}}
0152       },
0153       // Component #4
0154       {
0155           {{4.94537e+003,-2.08737e+003, 1.78089e+002, 2.29879e+001,-5.52783e+000,-1.86800e+000}},
0156           {{4.60220e+003,-1.62269e+003,-1.57552e+002, 2.01796e+002,-5.01636e+001, 6.47438e+000}},
0157           {{-9.50373e+004, 4.05517e+004,-5.62596e+003, 4.58534e+001, 6.70479e+001,-1.22430e+001}}
0158       },
0159       // Component #5
0160       {
0161           {{-1.04129e+003, 1.15222e+002,-2.70356e+001, 3.18611e+001,-7.78800e+000,-1.50242e+000}},
0162           {{-2.71361e+004, 2.00625e+004,-6.19444e+003, 1.10061e+003,-1.29354e+002, 1.08289e+001}},
0163           {{3.15252e+004,-3.31508e+004, 1.20371e+004,-2.23822e+003, 2.44396e+002,-2.09130e+001}}
0164       },
0165       // Component #6
0166       {
0167           {{1.27751e+004,-6.79813e+003, 1.24650e+003,-8.20622e+001,-2.33476e+000, 2.46459e-001}},
0168           {{3.64336e+005,-2.08457e+005, 4.33028e+004,-3.67825e+003, 4.22914e+001, 1.42701e+001}},
0169           {{-1.79298e+006, 1.01843e+006,-2.10037e+005, 1.82222e+004,-4.33573e+002,-2.72725e+001}}
0170       },
0171   }};
0172   // clang-format on
0173
0174   return {cdf_cmps6_order5_data,
0175           cdf_cmps6_order5_data,
0176           true,
0177           true,
0178           0.2,
0179           0.2,
0180           clampToRange,
0181           0.0001,
0182           0.002};
0183 }