File indexing completed on 2023-10-25 09:39:05
0001 #ifndef GeometryVector_newBasic3DVector_h
0002 #define GeometryVector_newBasic3DVector_h
0003
0004 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
0005 #include "DataFormats/GeometryVector/interface/Theta.h"
0006 #include "DataFormats/GeometryVector/interface/Phi.h"
0007 #include "DataFormats/GeometryVector/interface/PreciseFloatType.h"
0008 #include "DataFormats/GeometryVector/interface/CoordinateSets.h"
0009 #include "DataFormats/Math/interface/ExtVec.h"
0010 #include <iosfwd>
0011 #include <cmath>
0012
0013 namespace detailsBasic3DVector {
0014 inline float __attribute__((always_inline)) __attribute__((pure)) eta(float x, float y, float z) {
0015 float t(z / std::sqrt(x * x + y * y));
0016 return ::asinhf(t);
0017 }
0018 inline double __attribute__((always_inline)) __attribute__((pure)) eta(double x, double y, double z) {
0019 double t(z / std::sqrt(x * x + y * y));
0020 return ::asinh(t);
0021 }
0022 inline long double __attribute__((always_inline)) __attribute__((pure))
0023 eta(long double x, long double y, long double z) {
0024 long double t(z / std::sqrt(x * x + y * y));
0025 return ::asinhl(t);
0026 }
0027 }
0028
0029 template <typename T>
0030 class Basic3DVector {
0031 public:
0032 typedef T ScalarType;
0033 typedef Vec4<T> VectorType;
0034 typedef Vec4<T> MathVector;
0035 typedef Geom::Cylindrical2Cartesian<T> Cylindrical;
0036 typedef Geom::Spherical2Cartesian<T> Spherical;
0037 typedef Spherical Polar;
0038
0039
0040
0041
0042
0043 Basic3DVector() : v{0, 0, 0, 0} {}
0044
0045
0046 Basic3DVector(const Basic3DVector& p) : v(p.v) {}
0047
0048
0049 template <class U>
0050 Basic3DVector(const Basic3DVector<U>& p) : v{T(p.v[0]), T(p.v[1]), T(p.v[2]), T(p.v[3])} {}
0051
0052
0053 Basic3DVector(const Basic2DVector<T>& p) : v{p.x(), p.y(), 0} {}
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063 template <class OtherPoint>
0064 explicit Basic3DVector(const OtherPoint& p) : v{T(p.x()), T(p.y()), T(p.z())} {}
0065
0066
0067 Basic3DVector(MathVector const& iv) : v(iv) {}
0068
0069 template <class U>
0070 Basic3DVector(Vec4<U> const& iv) : v{T(iv[0]), T(iv[1]), T(iv[2]), T(iv[3])} {}
0071
0072
0073 Basic3DVector(const T& x, const T& y, const T& z, const T& w = 0) : v{x, y, z, w} {}
0074
0075
0076
0077
0078
0079 template <typename U>
0080 Basic3DVector(const Geom::Theta<U>& theta, const Geom::Phi<U>& phi, const T& r) {
0081 Polar p(theta.value(), phi.value(), r);
0082 v[0] = p.x();
0083 v[1] = p.y();
0084 v[2] = p.z();
0085 }
0086
0087 MathVector const& mathVector() const { return v; }
0088 MathVector& mathVector() { return v; }
0089
0090 T operator[](int i) const { return v[i]; }
0091 T& operator[](int i) { return v[i]; }
0092
0093
0094 T x() const { return v[0]; }
0095
0096
0097 T y() const { return v[1]; }
0098
0099
0100 T z() const { return v[2]; }
0101
0102 T w() const { return v[3]; }
0103
0104 Basic2DVector<T> xy() const { return ::xy(v); }
0105
0106
0107 bool operator==(const Basic3DVector& rh) const {
0108 auto res = v == rh.v;
0109 return res[0] & res[1] & res[2] & res[3];
0110 }
0111
0112
0113 T mag2() const { return ::dot(v, v); }
0114
0115
0116 T mag() const { return std::sqrt(mag2()); }
0117
0118
0119 T perp2() const { return ::dot2(v, v); }
0120
0121
0122 T perp() const { return std::sqrt(perp2()); }
0123
0124
0125 T transverse() const { return perp(); }
0126
0127
0128
0129
0130
0131 T barePhi() const { return std::atan2(y(), x()); }
0132 Geom::Phi<T> phi() const { return Geom::Phi<T>(barePhi()); }
0133
0134
0135
0136
0137
0138 T bareTheta() const { return std::atan2(perp(), z()); }
0139 Geom::Theta<T> theta() const { return Geom::Theta<T>(std::atan2(perp(), z())); }
0140
0141
0142
0143
0144
0145
0146 T eta() const { return detailsBasic3DVector::eta(x(), y(), z()); }
0147
0148
0149
0150
0151 Basic3DVector unit() const {
0152 T my_mag = mag2();
0153 return (0 != my_mag) ? (*this) * (T(1) / std::sqrt(my_mag)) : *this;
0154 }
0155
0156
0157
0158 template <class U>
0159 Basic3DVector& operator+=(const Basic3DVector<U>& p) {
0160 v = v + p.v;
0161 return *this;
0162 }
0163
0164
0165
0166 template <class U>
0167 Basic3DVector& operator-=(const Basic3DVector<U>& p) {
0168 v = v - p.v;
0169 return *this;
0170 }
0171
0172
0173 Basic3DVector operator-() const { return Basic3DVector(-v); }
0174
0175
0176 Basic3DVector& operator*=(T t) {
0177 v = t * v;
0178 return *this;
0179 }
0180
0181
0182 Basic3DVector& operator/=(T t) {
0183
0184 v = v / t;
0185 return *this;
0186 }
0187
0188
0189 T dot(const Basic3DVector& rh) const { return ::dot(v, rh.v); }
0190
0191
0192
0193
0194
0195
0196 template <class U>
0197 typename PreciseFloatType<T, U>::Type dot(const Basic3DVector<U>& lh) const {
0198 return Basic3DVector<typename PreciseFloatType<T, U>::Type>(*this).dot(
0199 Basic3DVector<typename PreciseFloatType<T, U>::Type>(lh));
0200 }
0201
0202
0203 Basic3DVector cross(const Basic3DVector& lh) const { return ::cross3(v, lh.v); }
0204
0205
0206
0207
0208
0209
0210 template <class U>
0211 Basic3DVector<typename PreciseFloatType<T, U>::Type> cross(const Basic3DVector<U>& lh) const {
0212 return Basic3DVector<typename PreciseFloatType<T, U>::Type>(*this).cross(
0213 Basic3DVector<typename PreciseFloatType<T, U>::Type>(lh));
0214 }
0215
0216 public:
0217 Vec4<T> v;
0218 } __attribute__((aligned(16)));
0219
0220 namespace geometryDetails {
0221 std::ostream& print3D(std::ostream& s, double x, double y, double z);
0222 }
0223
0224
0225 template <class T>
0226 inline std::ostream& operator<<(std::ostream& s, const Basic3DVector<T>& v) {
0227 return geometryDetails::print3D(s, v.x(), v.y(), v.z());
0228 }
0229
0230
0231 template <class T>
0232 inline Basic3DVector<T> operator+(const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
0233 return a.v + b.v;
0234 }
0235 template <class T>
0236 inline Basic3DVector<T> operator-(const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
0237 return a.v - b.v;
0238 }
0239
0240 template <class T, class U>
0241 inline Basic3DVector<typename PreciseFloatType<T, U>::Type> operator+(const Basic3DVector<T>& a,
0242 const Basic3DVector<U>& b) {
0243 typedef Basic3DVector<typename PreciseFloatType<T, U>::Type> RT;
0244 return RT(a).v + RT(b).v;
0245 }
0246
0247 template <class T, class U>
0248 inline Basic3DVector<typename PreciseFloatType<T, U>::Type> operator-(const Basic3DVector<T>& a,
0249 const Basic3DVector<U>& b) {
0250 typedef Basic3DVector<typename PreciseFloatType<T, U>::Type> RT;
0251 return RT(a).v - RT(b).v;
0252 }
0253
0254
0255 template <class T>
0256 inline T operator*(const Basic3DVector<T>& v1, const Basic3DVector<T>& v2) {
0257 return v1.dot(v2);
0258 }
0259
0260
0261 template <class T, class U>
0262 inline typename PreciseFloatType<T, U>::Type operator*(const Basic3DVector<T>& v1, const Basic3DVector<U>& v2) {
0263 return v1.dot(v2);
0264 }
0265
0266
0267
0268
0269 template <class T>
0270 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, T t) {
0271 return v.v * t;
0272 }
0273
0274
0275 template <class T>
0276 inline Basic3DVector<T> operator*(T t, const Basic3DVector<T>& v) {
0277 return v.v * t;
0278 }
0279
0280 template <class T, typename S>
0281 inline Basic3DVector<T> operator*(S t, const Basic3DVector<T>& v) {
0282 return static_cast<T>(t) * v;
0283 }
0284
0285 template <class T, typename S>
0286 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, S t) {
0287 return static_cast<T>(t) * v;
0288 }
0289
0290
0291
0292
0293 template <class T>
0294 inline Basic3DVector<T> operator/(const Basic3DVector<T>& v, T t) {
0295 return v.v / t;
0296 }
0297
0298 template <class T, typename S>
0299 inline Basic3DVector<T> operator/(const Basic3DVector<T>& v, S s) {
0300
0301 T t = s;
0302 return v / t;
0303 }
0304
0305 typedef Basic3DVector<float> Basic3DVectorF;
0306 typedef Basic3DVector<double> Basic3DVectorD;
0307
0308
0309 #include "Basic3DVectorLD.h"
0310
0311 #endif