EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-12-16 09:22:47
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 #pragma once 0010 0011 #include <cmath> 0012 #include <cstddef> 0013 #include <numbers> 0014 #include <span> 0015 #include <string> 0016 #include <vector> 0017 0018 namespace Acts { 0019 0020 namespace detail { 0021 0022 struct GaussianComponent { 0023 double weight = 0; 0024 double mean = 0; 0025 double var = 0; 0026 }; 0027 0028 /// Transform a gaussian component to a space where all values are defined from 0029 /// [-inf, inf] 0030 inline GaussianComponent transformComponent(const GaussianComponent &cmp) { 0031 GaussianComponent transformed_cmp; 0032 transformed_cmp.weight = std::log(cmp.weight) - std::log(1 - cmp.weight); 0033 transformed_cmp.mean = std::log(cmp.mean) - std::log(1 - cmp.mean); 0034 transformed_cmp.var = std::log(cmp.var); 0035 return transformed_cmp; 0036 } 0037 0038 /// Transform a gaussian component back from the [-inf, inf]-space to the usual 0039 /// space 0040 inline GaussianComponent inverseTransformComponent( 0041 const GaussianComponent &transformed_cmp) { 0042 GaussianComponent cmp; 0043 cmp.weight = 1 / (1 + std::exp(-transformed_cmp.weight)); 0044 cmp.mean = 1 / (1 + std::exp(-transformed_cmp.mean)); 0045 cmp.var = std::exp(transformed_cmp.var); 0046 return cmp; 0047 } 0048 0049 } // namespace detail 0050 0051 class BetheHeitlerApprox { 0052 public: 0053 using Component = detail::GaussianComponent; 0054 0055 virtual ~BetheHeitlerApprox() = default; 0056 0057 virtual std::size_t maxComponents() const = 0; 0058 0059 virtual bool validXOverX0(double xOverX0) const = 0; 0060 0061 virtual std::span<Component> mixture(double xOverX0, 0062 std::span<Component> mixture) const = 0; 0063 }; 0064 0065 /// This class approximates the Bethe-Heitler with only one component. This is 0066 /// mainly inside @ref AtlasBetheHeitlerApprox, but can also be used as the 0067 /// only component approximation (then probably for debugging) 0068 class BetheHeitlerApproxSingleCmp final : public BetheHeitlerApprox { 0069 public: 0070 /// Returns the number of components the returned mixture will have 0071 /// @return Number of components (always 1 for single component approximation) 0072 std::size_t maxComponents() const override { return 1; } 0073 0074 /// Checks if an input is valid for the parameterization. The threshold for 0075 /// x/x0 is 0.002 and orientates on the values used in ATLAS 0076 /// @param xOverX0 The x/x0 value to check 0077 /// @return True if x/x0 is below the threshold for single component approximation 0078 bool validXOverX0(const double xOverX0) const override { 0079 return xOverX0 < 0.002; 0080 } 0081 0082 /// Returns array with length 1 containing a 1-component-representation of the 0083 /// Bethe-Heitler-Distribution 0084 /// 0085 /// @param xOverX0 pathlength in terms of the radiation length 0086 /// @param mixture preallocated array to store the result 0087 /// @return the unmodified input span containing the single component 0088 std::span<Component> mixture( 0089 const double xOverX0, const std::span<Component> mixture) const override { 0090 mixture[0].weight = 1.0; 0091 0092 const double c = xOverX0 / std::numbers::ln2; 0093 mixture[0].mean = std::pow(2, -c); 0094 mixture[0].var = std::pow(3, -c) - std::pow(4, -c); 0095 0096 return mixture; 0097 } 0098 }; 0099 0100 /// This class approximates the Bethe-Heitler distribution as a gaussian 0101 /// mixture. To enable an approximation for continuous input variables, the 0102 /// weights, means and variances are internally parametrized as a Nth order 0103 /// polynomial. 0104 /// @todo This class is rather inflexible: It forces two data representations, 0105 /// making it a bit awkward to add a single parameterization. It would be good 0106 /// to generalize this at some point. 0107 class AtlasBetheHeitlerApprox : public BetheHeitlerApprox { 0108 public: 0109 struct PolyData { 0110 std::vector<double> weightCoeffs; 0111 std::vector<double> meanCoeffs; 0112 std::vector<double> varCoeffs; 0113 }; 0114 0115 /// Type alias for array of polynomial data for all components 0116 using Data = std::vector<PolyData>; 0117 0118 /// Loads a parameterization from a file according to the Atlas file 0119 /// description 0120 /// 0121 /// @param low_parameters_path Path to the foo.par file that stores 0122 /// the parameterization for low x/x0 0123 /// @param high_parameters_path Path to the foo.par file that stores 0124 /// the parameterization for high x/x0 0125 /// @param lowLimit the upper limit for the low x/x0-data 0126 /// @param highLimit the upper limit for the high x/x0-data 0127 /// @param clampToRange forwarded to constructor 0128 /// @param noChangeLimit forwarded to constructor 0129 /// @param singleGaussianLimit forwarded to constructor 0130 /// @return AtlasBetheHeitlerApprox instance loaded from parameter files 0131 static AtlasBetheHeitlerApprox loadFromFiles( 0132 const std::string &low_parameters_path, 0133 const std::string &high_parameters_path, double lowLimit, 0134 double highLimit, bool clampToRange, double noChangeLimit, 0135 double singleGaussianLimit); 0136 0137 /// Construct the Bethe-Heitler approximation description with two 0138 /// parameterizations, one for lower ranges, one for higher ranges. 0139 /// Is it assumed that the lower limit of the high-x/x0 data is equal 0140 /// to the upper limit of the low-x/x0 data. 0141 /// 0142 /// @param lowData data for the lower x/x0 range 0143 /// @param highData data for the higher x/x0 range 0144 /// @param lowTransform whether the low data need to be transformed 0145 /// @param highTransform whether the high data need to be transformed 0146 /// @param lowLimit the upper limit for the low data 0147 /// @param highLimit the upper limit for the high data 0148 /// @param clampToRange whether to clamp the input x/x0 to the allowed range 0149 /// @param noChangeLimit limit below which no change is applied 0150 /// @param singleGaussianLimit limit below which a single Gaussian is used 0151 AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData, 0152 bool lowTransform, bool highTransform, 0153 double lowLimit, double highLimit, bool clampToRange, 0154 double noChangeLimit, double singleGaussianLimit) 0155 : m_lowData(lowData), 0156 m_highData(highData), 0157 m_lowTransform(lowTransform), 0158 m_highTransform(highTransform), 0159 m_lowLimit(lowLimit), 0160 m_highLimit(highLimit), 0161 m_clampToRange(clampToRange), 0162 m_noChangeLimit(noChangeLimit), 0163 m_singleGaussianLimit(singleGaussianLimit) {} 0164 0165 /// Returns the number of components the returned mixture will have 0166 /// @return Number of components in the mixture 0167 std::size_t maxComponents() const override { 0168 return std::max(m_lowData.size(), m_highData.size()); 0169 } 0170 0171 /// Checks if an input is valid for the parameterization 0172 /// 0173 /// @param xOverX0 pathlength in terms of the radiation length 0174 /// @return True if x/x0 is within valid range for this parameterization 0175 bool validXOverX0(const double xOverX0) const override { 0176 if (m_clampToRange) { 0177 return true; 0178 } else { 0179 return xOverX0 < m_highLimit; 0180 } 0181 } 0182 0183 /// Generates the mixture from the polynomials and reweights them, so 0184 /// that the sum of all weights is 1 0185 /// 0186 /// @param xOverX0 pathlength in terms of the radiation length 0187 /// @param mixture preallocated array to store the result 0188 /// @return the potentially modified input span containing the mixture 0189 std::span<Component> mixture( 0190 double xOverX0, const std::span<Component> mixture) const override; 0191 0192 private: 0193 Data m_lowData; 0194 Data m_highData; 0195 bool m_lowTransform = false; 0196 bool m_highTransform = false; 0197 double m_lowLimit = 0; 0198 double m_highLimit = 0; 0199 bool m_clampToRange = false; 0200 double m_noChangeLimit = 0; 0201 double m_singleGaussianLimit = 0; 0202 }; 0203 0204 /// Creates a @ref AtlasBetheHeitlerApprox object based on an ATLAS 0205 /// configuration, that are stored as static data in the source code. 0206 /// This may not be an optimal configuration, but should allow to run 0207 /// the GSF without the need to load files 0208 /// @param clampToRange Whether to clamp values to the valid range 0209 /// @return AtlasBetheHeitlerApprox with default ATLAS configuration parameters 0210 AtlasBetheHeitlerApprox makeDefaultBetheHeitlerApprox( 0211 bool clampToRange = false); 0212 0213 } // namespace Acts
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|