include/edm4eic/vector_utils_legacy.h

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