edm4hep/utils/vector_utils_legacy.h

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