Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:11:08

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 <array>
0012 
0013 namespace Acts::Concepts {
0014 /// @brief check types for requirements needed by interpolation
0015 ///
0016 /// @tparam Point1 type for specifying geometric positions
0017 /// @tparam Point2 type for specifying geometric positions
0018 /// @tparam Point3 type for specifying geometric positions
0019 /// @tparam Value  type of values to be interpolated
0020 ///
0021 /// This helper struct provides compile-time information whether the provided
0022 /// @c Point and @c Value types can be used in the Acts::interpolate function.
0023 ///
0024 /// The following boolean variable
0025 /// @code{.cpp}
0026 /// Acts::detail::can_interpolate<Point1,Point2,Point3,Value>::value
0027 /// @endcode
0028 ///
0029 /// is @c true if all @c Point types and @c Value fulfill the type
0030 /// requirements for being used in the interpolation function, otherwise it is
0031 /// @c false.
0032 template <typename Point1, typename Point2, typename Point3, typename Value>
0033 struct can_interpolate {
0034   template <typename C>
0035   static auto value_type_test(C* c)
0036       -> decltype(static_cast<C>(std::declval<double>() * std::declval<C>() +
0037                                  std::declval<double>() * std::declval<C>()),
0038                   std::true_type());
0039   template <typename C>
0040   static std::false_type value_type_test(...);
0041 
0042   template <typename C>
0043   static auto point_type_test(C* c)
0044       -> decltype(static_cast<double>(std::declval<C>()[0]), std::true_type());
0045   template <typename C>
0046   static std::false_type point_type_test(...);
0047 
0048   static const bool value =
0049       std::is_same_v<std::true_type,
0050                      decltype(value_type_test<Value>(nullptr))> &&
0051       std::is_same_v<std::true_type,
0052                      decltype(point_type_test<Point1>(nullptr))> &&
0053       std::is_same_v<std::true_type,
0054                      decltype(point_type_test<Point2>(nullptr))> &&
0055       std::is_same_v<std::true_type,
0056                      decltype(point_type_test<Point3>(nullptr))>;
0057 };
0058 
0059 /// @brief concept equivalent to `can_interpolate`
0060 /// @todo this is a concept based on a traditional type trait, meaning it won't
0061 /// have nice errors; should be rewritten as a real concept!
0062 template <typename T, typename P1, typename P2, typename P3>
0063 concept interpolatable = can_interpolate<P1, P2, P3, T>::value;
0064 }  // namespace Acts::Concepts
0065 
0066 namespace Acts::detail {
0067 
0068 /// @brief determine number of dimension from power of 2
0069 ///
0070 /// @tparam N power of 2
0071 template <std::size_t N>
0072 struct get_dimension {
0073   /// exponent @c d such that \f$2^d = N \f$
0074   static constexpr std::size_t value = get_dimension<(N >> 1)>::value + 1u;
0075 };
0076 
0077 /// @cond
0078 template <>
0079 struct get_dimension<2u> {
0080   static constexpr std::size_t value = 1u;
0081 };
0082 /// @endcond
0083 
0084 /// @brief helper struct for performing multi-dimensional linear interpolation
0085 ///
0086 /// @tparam T      type of values to be interpolated
0087 /// @tparam Point1 type specifying geometric positions
0088 /// @tparam Point2 type specifying geometric positions
0089 /// @tparam Point3 type specifying geometric positions
0090 /// @tparam D      current dimension on which to perform reduction
0091 /// @tparam N      number of hyper box corners
0092 ///
0093 /// @note
0094 /// - Given @c U and @c V of value type @c T as well as two @c double @c a and
0095 /// @c b, then the following must be a valid expression <tt>a * U + b * V</tt>
0096 /// yielding an object which is (implicitly) convertible to @c T.
0097 /// - The @c Point types must represent d-dimensional positions and support
0098 /// coordinate access using @c operator[]. Coordinate indices must start at 0.
0099 /// - @c N is the number of hyper box corners which is \f$2^d\f$ where \f$d\f$
0100 /// is the dimensionality of the hyper box. The dimensionality must be
0101 /// consistent with the provided @c Point types.
0102 template <typename T, class Point1, class Point2, class Point3, std::size_t D,
0103           std::size_t N>
0104 struct interpolate_impl;
0105 
0106 /// @cond
0107 // recursive implementation of linear interpolation in multiple dimensions
0108 template <typename T, class Point1, class Point2, class Point3, std::size_t D,
0109           std::size_t N>
0110 struct interpolate_impl {
0111   static T run(const Point1& pos, const Point2& lowerLeft,
0112                const Point3& upperRight, const std::array<T, N>& fields) {
0113     // get distance to lower boundary relative to total bin width
0114     const double f = (pos[D] - lowerLeft[D]) / (upperRight[D] - lowerLeft[D]);
0115 
0116     std::array<T, (N >> 1)> newFields{};
0117     for (std::size_t i = 0; i < N / 2; ++i) {
0118       newFields.at(i) = (1 - f) * fields.at(2 * i) + f * fields.at(2 * i + 1);
0119     }
0120 
0121     return interpolate_impl<T, Point1, Point2, Point3, D - 1, (N >> 1)>::run(
0122         pos, lowerLeft, upperRight, newFields);
0123   }
0124 };
0125 
0126 // simple linear interpolation in 1D
0127 template <typename T, class Point1, class Point2, class Point3, std::size_t D>
0128 struct interpolate_impl<T, Point1, Point2, Point3, D, 2u> {
0129   static T run(const Point1& pos, const Point2& lowerLeft,
0130                const Point3& upperRight, const std::array<T, 2u>& fields) {
0131     // get distance to lower boundary relative to total bin width
0132     const double f = (pos[D] - lowerLeft[D]) / (upperRight[D] - lowerLeft[D]);
0133 
0134     return (1 - f) * fields.at(0) + f * fields.at(1);
0135   }
0136 };
0137 /// @endcond
0138 }  // namespace Acts::detail