include/edm4eic/vector_utils.h

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Sylvester Joosten
0003
0004 #ifndef EDM4EIC_UTILS_VECTOR_HH
0005 #define EDM4EIC_UTILS_VECTOR_HH
0006
0007 // These ultilies require concepts. If not available, use the fallback
0008 // vector_utils_legacy.h instead to capture most functionality.
0009 #if !__cpp_concepts
0010 #include <edm4eic/vector_utils_legacy.h>
0011 #else
0012 #include <cmath>
0013
0014 #include <edm4hep/Vector3f.h>
0015
0016 namespace edm4eic {
0017
0018 template <class V> concept VectorHasX = requires(V v) { v.x; };
0019 template <class V> concept VectorHasY = requires(V v) { v.y; };
0020 template <class V> concept VectorHasZ = requires(V v) { v.z; };
0021 template <class V> concept VectorHasA = requires(V v) { v.a; };
0022 template <class V> concept VectorHasB = requires(V v) { v.b; };
0023 template <class V> concept ClassVector = requires(V v) { v.x(); };
0024 template <class V>
0025 concept Vector2D_XY =
0026     VectorHasX<V>&& VectorHasY<V> && !VectorHasZ<V> && !ClassVector<V>;
0027 template <class V>
0028 concept Vector2D_AB =
0029     VectorHasA<V>&& VectorHasB<V> && !VectorHasZ<V> && !ClassVector<V>;
0030 template <class V> concept Vector2D = Vector2D_XY<V> || Vector2D_AB<V>;
0031 template <class V>
0032 concept Vector3D =
0033     VectorHasX<V>&& VectorHasY<V>&& VectorHasZ<V> && !ClassVector<V>;
0034 template <class V> concept VectorND = Vector2D<V> || Vector3D<V>;
0035 template <class V> concept VectorND_XYZ = Vector2D_XY<V> || Vector3D<V>;
0036
0037 inline double etaToAngle(const double eta) {
0038   return std::atan(std::exp(-eta)) * 2.;
0039 }
0040 inline double angleToEta(const double theta) {
0041   return -std::log(std::tan(0.5 * theta));
0042 }
0043
0044 // Utility getters to accomodate different vector types
0045 template <VectorND_XYZ V> auto vector_x(const V& v) { return v.x; }
0046 template <VectorND_XYZ V> auto vector_y(const V& v) { return v.y; }
0047 template <Vector3D V> auto vector_z(const V& v) { return v.z; }
0048 template <Vector2D_AB V> auto vector_x(const V& v) { return v.a; }
0049 template <Vector2D_AB V> auto vector_y(const V& v) { return v.b; }
0050
0051 // inline edm4hep::Vector2f VectorFromPolar(const double r, const double theta)
0052 // {
0053 //  return {r * sin(theta), r * cos(theta)};
0054 //}
0055 template <Vector3D V = edm4hep::Vector3f>
0056 V sphericalToVector(const double r, const double theta, const double phi) {
0057   using FloatType = decltype(vector_x(V()));
0058   const double sth = sin(theta);
0059   const double cth = cos(theta);
0060   const double sph = sin(phi);
0061   const double cph = cos(phi);
0062   const FloatType x = r * sth * cph;
0063   const FloatType y = r * sth * sph;
0064   const FloatType z = r * cth;
0065   return {x, y, z};
0066 }
0067
0068 template <Vector3D V> double anglePolar(const V& v) {
0069   return std::atan2(std::hypot(vector_x(v), vector_y(v)), vector_z(v));
0070 }
0071 template <VectorND V> double angleAzimuthal(const V& v) {
0072   return std::atan2(vector_y(v), vector_x(v));
0073 }
0074 template <Vector3D V> double eta(const V& v) {
0075   return angleToEta(anglePolar(v));
0076 }
0077 template <Vector2D V> double magnitude(const V& v) {
0078   return std::hypot(vector_x(v), vector_y(v));
0079 }
0080 template <Vector3D V> double magnitude(const V& v) {
0081   return std::hypot(vector_x(v), vector_y(v), vector_z(v));
0082 }
0083 template <Vector3D V> double magnitudeTransverse(const V& v) {
0084   return std::hypot(vector_x(v), vector_y(v));
0085 }
0086 template <Vector3D V> double magnitudeLongitudinal(const V& v) {
0087   return vector_z(v);
0088 }
0089 template <VectorND V> V normalizeVector(const V& v, double norm = 1.) {
0090   const double old = magnitude(v);
0091   if (old == 0) {
0092     return v;
0093   }
0094   return (norm > 0) ? v * norm / old : v * 0;
0095 }
0096 template <Vector3D V> V vectorTransverse(const V& v) {
0097   return {vector_x(v), vector_y(v), 0};
0098 }
0099 template <Vector3D V> V vectorLongitudinal(const V& v) {
0100   return {0, 0, vector_z(v)};
0101 }
0102 // Two vector functions
0103 template <VectorND V> double angleBetween(const V& v1, const V& v2) {
0104   const double dot = v1 * v2;
0105   if (dot == 0) {
0106     return 0.;
0107   }
0108   return acos(dot / (magnitude(v1) * magnitude(v2)));
0109 }
0110 // Project v onto v1
0111 template <Vector3D V> double projection(const V& v, const V& v1) {
0112   const double norm = magnitude(v1);
0113   if (norm == 0) {
0114     return magnitude(v);
0115   }
0116   return v * v1 / norm;
0117 }
0118
0119 template <edm4eic::Vector2D V> V operator+(const V& v1, const V& v2) {
0120   return {edm4eic::vector_x(v1) + edm4eic::vector_x(v2),
0121           edm4eic::vector_y(v1) + edm4eic::vector_y(v2)};
0122 }
0123 template <edm4eic::Vector3D V> V operator+(const V& v1, const V& v2) {
0124   return {edm4eic::vector_x(v1) + edm4eic::vector_x(v2),
0125           edm4eic::vector_y(v1) + edm4eic::vector_y(v2),
0126           edm4eic::vector_z(v1) + edm4eic::vector_z(v2)};
0127 }
0128 template <edm4eic::Vector2D V> double operator*(const V& v1, const V& v2) {
0129   return edm4eic::vector_x(v1) * edm4eic::vector_x(v2) +
0130          edm4eic::vector_y(v1) * edm4eic::vector_y(v2);
0131 }
0132 template <edm4eic::Vector3D V> double operator*(const V& v1, const V& v2) {
0133   return edm4eic::vector_x(v1) * edm4eic::vector_x(v2) +
0134          edm4eic::vector_y(v1) * edm4eic::vector_y(v2) +
0135          edm4eic::vector_z(v1) * edm4eic::vector_z(v2);
0136 }
0137 template <edm4eic::Vector2D V> V operator*(const double d, const V& v) {
0138   return {d * edm4eic::vector_x(v), d * edm4eic::vector_y(v)};
0139 }
0140 template <edm4eic::Vector3D V> V operator*(const double d, const V& v) {
0141   return {d * edm4eic::vector_x(v), d * edm4eic::vector_y(v), d * edm4eic::vector_z(v)};
0142 }
0143 template <edm4eic::Vector2D V> V operator*(const V& v, const double d) {
0144   return {d * edm4eic::vector_x(v), d * edm4eic::vector_y(v)};
0145 }
0146 template <edm4eic::Vector3D V> V operator*(const V& v, const double d) {
0147   return {d * edm4eic::vector_x(v), d * edm4eic::vector_y(v), d * edm4eic::vector_z(v)};
0148 }
0149 template <edm4eic::VectorND V> V operator-(const V& v1, const V& v2) {
0150   return v1 + (-1. * v2);
0151 }
0152 template <edm4eic::VectorND V> V operator/(const V& v, const double d) {
0153   return (1. / d) * v;
0154 }
0155
0156 } // namespace edm4eic
0157
0158 #endif
0159 #endif