edm4hep/utils/vector_utils.h

0001 #ifndef EDM4HEP_UTILS_VECTOR_HH
0002 #define EDM4HEP_UTILS_VECTOR_HH
0003
0004 #include <edm4hep/Vector3f.h>
0005
0006 #include <cmath>
0007
0008 namespace edm4hep {
0009
0010 template <class V>
0011 concept VectorHasX = requires(V v) { v.x; };
0012
0013 template <class V>
0014 concept VectorHasY = requires(V v) { v.y; };
0015
0016 template <class V>
0017 concept VectorHasZ = requires(V v) { v.z; };
0018
0019 template <class V>
0020 concept VectorHasT = requires(V v) { v.t; };
0021
0022 template <class V>
0023 concept VectorHasA = requires(V v) { v.a; };
0024
0025 template <class V>
0026 concept VectorHasB = requires(V v) { v.b; };
0027
0028 template <class V>
0029 concept ClassVector = requires(V v) { v.x(); };
0030
0031 template <class V>
0032 concept Vector2D_XY = VectorHasX<V> && VectorHasY<V> && !VectorHasZ<V> && !VectorHasT<V> && !ClassVector<V>;
0033
0034 template <class V>
0035 concept Vector2D_AB = VectorHasA<V> && VectorHasB<V> && !VectorHasZ<V> && !VectorHasT<V> && !ClassVector<V>;
0036
0037 template <class V>
0038 concept Vector2D = Vector2D_XY<V> || Vector2D_AB<V>;
0039
0040 template <class V>
0041 concept Vector3D = VectorHasX<V> && VectorHasY<V> && VectorHasZ<V> && !VectorHasT<V> && !ClassVector<V>;
0042
0043 template <class V>
0044 concept Vector4D = VectorHasX<V> && VectorHasY<V> && VectorHasZ<V> && VectorHasT<V> && !ClassVector<V>;
0045
0046 template <class V>
0047 concept Vector2or3D = Vector2D<V> || Vector3D<V>;
0048
0049 template <class V>
0050 concept VectorND = Vector2D<V> || Vector3D<V> || Vector4D<V>;
0051
0052 template <class V>
0053 concept VectorND_XYZ = Vector2D_XY<V> || Vector3D<V> || Vector4D<V>;
0054
0055 namespace utils {
0056   // Utility getters to accommodate different vector types
0057   template <VectorND_XYZ V>
0058   constexpr auto vector_x(const V& v) {
0059     return v.x;
0060   }
0061
0062   template <VectorND_XYZ V>
0063   constexpr auto vector_y(const V& v) {
0064     return v.y;
0065   }
0066
0067   template <Vector2D V>
0068   constexpr auto vector_z(const V&) {
0069     return 0;
0070   }
0071
0072   template <Vector3D V>
0073   constexpr auto vector_z(const V& v) {
0074     return v.z;
0075   }
0076
0077   template <Vector4D V>
0078   constexpr auto vector_z(const V& v) {
0079     return v.z;
0080   }
0081
0082   template <Vector4D V>
0083   constexpr auto vector_t(const V& v) {
0084     return v.t;
0085   }
0086
0087   template <Vector2or3D V>
0088   constexpr auto vector_t(const V&) {
0089     return 0;
0090   }
0091
0092   template <Vector2D_AB V>
0093   constexpr auto vector_x(const V& v) {
0094     return v.a;
0095   }
0096
0097   template <Vector2D_AB V>
0098   constexpr auto vector_y(const V& v) {
0099     return v.b;
0100   }
0101
0102   namespace detail {
0103     /// Helper struct to determine the underlying type of edm4hep vector types
0104     template <typename V>
0105     struct ValueTypeHelper {
0106       using type = decltype(vector_x(std::declval<V>()));
0107     };
0108   } // namespace detail
0109
0110   /// Type alias that returns the underlying type of edm4hep
0111   template <typename V>
0112   using ValueType = typename detail::ValueTypeHelper<V>::type;
0113
0114   // inline edm4hep::Vector2f VectorFromPolar(const double r, const double theta)
0115   // {
0116   //  return {r * std::sin(theta), r * std::cos(theta)};
0117   //}
0118
0119 } // namespace utils
0120
0121 /// A vector that uses floating point numbers to represent its members
0122 template <typename V>
0123 concept FloatVectorND = std::floating_point<utils::ValueType<V>> && (Vector2D<V> || Vector3D<V>);
0124
0125 namespace utils {
0126   inline double etaToAngle(const double eta) { return std::atan(std::exp(-eta)) * 2.; }
0127
0128   inline double angleToEta(const double theta) { return -std::log(std::tan(0.5 * theta)); }
0129
0130   template <Vector3D V = edm4hep::Vector3f>
0131   V sphericalToVector(const double r, const double theta, const double phi) {
0132     using FloatType = ValueType<V>;
0133     const double sth = std::sin(theta);
0134     const double cth = std::cos(theta);
0135     const double sph = std::sin(phi);
0136     const double cph = std::cos(phi);
0137     const FloatType x = r * sth * cph;
0138     const FloatType y = r * sth * sph;
0139     const FloatType z = r * cth;
0140     return {x, y, z};
0141   }
0142
0143   template <Vector3D V>
0144   double anglePolar(const V& v) {
0145     return std::atan2(std::hypot(vector_x(v), vector_y(v)), vector_z(v));
0146   }
0147
0148   template <VectorND V>
0149   double angleAzimuthal(const V& v) {
0150     return std::atan2(vector_y(v), vector_x(v));
0151   }
0152
0153   template <Vector3D V>
0154   double eta(const V& v) {
0155     return angleToEta(anglePolar(v));
0156   }
0157
0158   template <Vector2or3D V>
0159   double magnitude(const V& v) {
0160     return std::hypot(vector_x(v), vector_y(v), vector_z(v));
0161   }
0162
0163   template <Vector4D V>
0164   double magnitude(const V& v) {
0165     return std::hypot(vector_x(v), vector_y(v), vector_z(v)) - vector_t(v) * vector_t(v);
0166   }
0167
0168   template <Vector3D V>
0169   double magnitudeTransverse(const V& v) {
0170     return std::hypot(vector_x(v), vector_y(v));
0171   }
0172
0173   template <Vector3D V>
0174   double magnitudeLongitudinal(const V& v) {
0175     return vector_z(v);
0176   }
0177
0178   template <FloatVectorND V>
0179   V normalizeVector(const V& v, double norm = 1.) {
0180     const double old = magnitude(v);
0181     if (old == 0) {
0182       return v;
0183     }
0184     return (norm > 0) ? v * norm / old : v * 0;
0185   }
0186
0187   template <Vector3D V>
0188   constexpr V vectorTransverse(const V& v) {
0189     return {vector_x(v), vector_y(v), 0};
0190   }
0191
0192   template <Vector3D V>
0193   constexpr V vectorLongitudinal(const V& v) {
0194     return {0, 0, vector_z(v)};
0195   }
0196
0197   // Two vector functions
0198   template <Vector2or3D V>
0199   double angleBetween(const V& v1, const V& v2) {
0200     const double dot = v1 * v2;
0201     if (dot == 0) {
0202       return 0.;
0203     }
0204     return std::acos(dot / (magnitude(v1) * magnitude(v2)));
0205   }
0206
0207   // Project v onto v1
0208   template <Vector2or3D V>
0209   double projection(const V& v, const V& v1) {
0210     const double norm = magnitude(v1);
0211     if (norm == 0) {
0212       return magnitude(v);
0213     }
0214     return v * v1 / norm;
0215   }
0216
0217 } // namespace utils
0218
0219 template <edm4hep::Vector2D V>
0220 constexpr V operator+(const V& v1, const V& v2) {
0221   return {edm4hep::utils::vector_x(v1) + edm4hep::utils::vector_x(v2),
0222           edm4hep::utils::vector_y(v1) + edm4hep::utils::vector_y(v2)};
0223 }
0224
0225 template <edm4hep::Vector3D V>
0226 constexpr V operator+(const V& v1, const V& v2) {
0227   return {edm4hep::utils::vector_x(v1) + edm4hep::utils::vector_x(v2),
0228           edm4hep::utils::vector_y(v1) + edm4hep::utils::vector_y(v2),
0229           edm4hep::utils::vector_z(v1) + edm4hep::utils::vector_z(v2)};
0230 }
0231
0232 template <edm4hep::Vector4D V>
0233 constexpr V operator+(const V& v1, const V& v2) {
0234   return {edm4hep::utils::vector_x(v1) + edm4hep::utils::vector_x(v2),
0235           edm4hep::utils::vector_y(v1) + edm4hep::utils::vector_y(v2),
0236           edm4hep::utils::vector_z(v1) + edm4hep::utils::vector_z(v2),
0237           edm4hep::utils::vector_t(v1) + edm4hep::utils::vector_t(v2)};
0238 }
0239
0240 template <edm4hep::Vector2D V>
0241 constexpr auto operator*(const V& v1, const V& v2) {
0242   return edm4hep::utils::vector_x(v1) * edm4hep::utils::vector_x(v2) +
0243          edm4hep::utils::vector_y(v1) * edm4hep::utils::vector_y(v2);
0244 }
0245
0246 template <edm4hep::Vector3D V>
0247 constexpr auto operator*(const V& v1, const V& v2) {
0248   return edm4hep::utils::vector_x(v1) * edm4hep::utils::vector_x(v2) +
0249          edm4hep::utils::vector_y(v1) * edm4hep::utils::vector_y(v2) +
0250          edm4hep::utils::vector_z(v1) * edm4hep::utils::vector_z(v2);
0251 }
0252
0253 template <edm4hep::Vector4D V>
0254 constexpr auto operator*(const V& v1, const V& v2) {
0255   return (edm4hep::utils::vector_x(v1) * edm4hep::utils::vector_x(v2) +
0256           edm4hep::utils::vector_y(v1) * edm4hep::utils::vector_y(v2) +
0257           edm4hep::utils::vector_z(v1) * edm4hep::utils::vector_z(v2)) -
0258          edm4hep::utils::vector_t(v1) * edm4hep::utils::vector_t(v2);
0259 }
0260
0261 template <edm4hep::Vector2D V>
0262 constexpr V operator*(const double d, const V& v) {
0263   using VT = edm4hep::utils::ValueType<V>;
0264   const VT x = d * edm4hep::utils::vector_x(v);
0265   const VT y = d * edm4hep::utils::vector_y(v);
0266   return {x, y};
0267 }
0268
0269 template <edm4hep::Vector3D V>
0270 constexpr V operator*(const double d, const V& v) {
0271   using VT = edm4hep::utils::ValueType<V>;
0272   const VT x = d * edm4hep::utils::vector_x(v);
0273   const VT y = d * edm4hep::utils::vector_y(v);
0274   const VT z = d * edm4hep::utils::vector_z(v);
0275   return {x, y, z};
0276 }
0277
0278 template <edm4hep::Vector4D V>
0279 constexpr V operator*(const double d, const V& v) {
0280   using VT = edm4hep::utils::ValueType<V>;
0281   const VT x = d * edm4hep::utils::vector_x(v);
0282   const VT y = d * edm4hep::utils::vector_y(v);
0283   const VT z = d * edm4hep::utils::vector_z(v);
0284   const VT t = d * edm4hep::utils::vector_t(v);
0285   return {x, y, z, t};
0286 }
0287
0288 template <edm4hep::VectorND V>
0289 constexpr V operator*(const V& v, const double d) {
0290   return d * v;
0291 }
0292
0293 template <edm4hep::VectorND V>
0294 constexpr V operator-(const V& v1, const V& v2) {
0295   return v1 + (-1. * v2);
0296 }
0297
0298 template <edm4hep::VectorND V>
0299 constexpr V operator/(const V& v, const double d) {
0300   return (1. / d) * v;
0301 }
0302
0303 } // namespace edm4hep
0304
0305 #endif