File indexing completed on 2024-04-06 12:04:16
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/SSEVec.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 mathSSE::Vec4<T> VectorType;
0034 typedef mathSSE::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() {}
0044
0045
0046 Basic3DVector(const Basic3DVector& p) : v(p.v) {}
0047
0048
0049 Basic3DVector& operator=(const Basic3DVector&) = default;
0050
0051
0052 template <class U>
0053 Basic3DVector(const Basic3DVector<U>& p) : v(p.v) {}
0054
0055
0056 Basic3DVector(const Basic2DVector<T>& p) : v(p.x(), p.y(), 0) {}
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066 template <class OtherPoint>
0067 explicit Basic3DVector(const OtherPoint& p) : v(p.x(), p.y(), p.z()) {}
0068
0069
0070 template <class U>
0071 Basic3DVector(mathSSE::Vec4<U> const& iv) : v(iv) {}
0072
0073
0074 Basic3DVector(const T& x, const T& y, const T& z, const T& w = 0) : v(x, y, z, w) {}
0075
0076
0077
0078
0079
0080 template <typename U>
0081 Basic3DVector(const Geom::Theta<U>& theta, const Geom::Phi<U>& phi, const T& r) {
0082 Polar p(theta.value(), phi.value(), r);
0083 v.o.theX = p.x();
0084 v.o.theY = p.y();
0085 v.o.theZ = p.z();
0086 }
0087
0088 MathVector const& mathVector() const { return v; }
0089 MathVector& mathVector() { return v; }
0090
0091 T operator[](int i) const { return v[i]; }
0092 T& operator[](int i) { return v[i]; }
0093
0094
0095 T x() const { return v.o.theX; }
0096
0097
0098 T y() const { return v.o.theY; }
0099
0100
0101 T z() const { return v.o.theZ; }
0102
0103 T w() const { return v.o.theW; }
0104
0105 Basic2DVector<T> xy() const { return v.xy(); }
0106
0107
0108 bool operator==(const Basic3DVector& rh) const { return v == rh.v; }
0109
0110
0111 T mag2() const { return ::dot(v, v); }
0112
0113
0114 T mag() const { return std::sqrt(mag2()); }
0115
0116
0117 T perp2() const { return ::dotxy(v, v); }
0118
0119
0120 T perp() const { return std::sqrt(perp2()); }
0121
0122
0123 T transverse() const { return perp(); }
0124
0125
0126
0127
0128
0129 T barePhi() const { return std::atan2(y(), x()); }
0130 Geom::Phi<T> phi() const { return Geom::Phi<T>(barePhi()); }
0131
0132
0133
0134
0135
0136 T bareTheta() const { return std::atan2(perp(), z()); }
0137 Geom::Theta<T> theta() const { return Geom::Theta<T>(std::atan2(perp(), z())); }
0138
0139
0140
0141
0142
0143
0144 T eta() const { return detailsBasic3DVector::eta(x(), y(), z()); }
0145
0146
0147
0148
0149 Basic3DVector unit() const {
0150 T my_mag = mag2();
0151 return (0 != my_mag) ? (*this) * (T(1) / std::sqrt(my_mag)) : *this;
0152 }
0153
0154
0155
0156 template <class U>
0157 Basic3DVector& operator+=(const Basic3DVector<U>& p) {
0158 v = v + p.v;
0159 return *this;
0160 }
0161
0162
0163
0164 template <class U>
0165 Basic3DVector& operator-=(const Basic3DVector<U>& p) {
0166 v = v - p.v;
0167 return *this;
0168 }
0169
0170
0171 Basic3DVector operator-() const { return Basic3DVector(-v); }
0172
0173
0174 Basic3DVector& operator*=(T t) {
0175 v = t * v;
0176 return *this;
0177 }
0178
0179
0180 Basic3DVector& operator/=(T t) {
0181
0182 v = v / t;
0183 return *this;
0184 }
0185
0186
0187 T dot(const Basic3DVector& rh) const { return ::dot(v, rh.v); }
0188
0189
0190
0191
0192
0193
0194 template <class U>
0195 typename PreciseFloatType<T, U>::Type dot(const Basic3DVector<U>& lh) const {
0196 return Basic3DVector<typename PreciseFloatType<T, U>::Type>(*this).dot(
0197 Basic3DVector<typename PreciseFloatType<T, U>::Type>(lh));
0198 }
0199
0200
0201 Basic3DVector cross(const Basic3DVector& lh) const { return ::cross(v, lh.v); }
0202
0203
0204
0205
0206
0207
0208 template <class U>
0209 Basic3DVector<typename PreciseFloatType<T, U>::Type> cross(const Basic3DVector<U>& lh) const {
0210 return Basic3DVector<typename PreciseFloatType<T, U>::Type>(*this).cross(
0211 Basic3DVector<typename PreciseFloatType<T, U>::Type>(lh));
0212 }
0213
0214 public:
0215 mathSSE::Vec4<T> v;
0216 } __attribute__((aligned(16)));
0217
0218 namespace geometryDetails {
0219 std::ostream& print3D(std::ostream& s, double x, double y, double z);
0220 }
0221
0222
0223 template <class T>
0224 inline std::ostream& operator<<(std::ostream& s, const Basic3DVector<T>& v) {
0225 return geometryDetails::print3D(s, v.x(), v.y(), v.z());
0226 }
0227
0228
0229 template <class T>
0230 inline Basic3DVector<T> operator+(const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
0231 return a.v + b.v;
0232 }
0233 template <class T>
0234 inline Basic3DVector<T> operator-(const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
0235 return a.v - b.v;
0236 }
0237
0238 template <class T, class U>
0239 inline Basic3DVector<typename PreciseFloatType<T, U>::Type> operator+(const Basic3DVector<T>& a,
0240 const Basic3DVector<U>& b) {
0241 typedef Basic3DVector<typename PreciseFloatType<T, U>::Type> RT;
0242 return RT(a).v + RT(b).v;
0243 }
0244
0245 template <class T, class U>
0246 inline Basic3DVector<typename PreciseFloatType<T, U>::Type> operator-(const Basic3DVector<T>& a,
0247 const Basic3DVector<U>& b) {
0248 typedef Basic3DVector<typename PreciseFloatType<T, U>::Type> RT;
0249 return RT(a).v - RT(b).v;
0250 }
0251
0252
0253 template <class T>
0254 inline T operator*(const Basic3DVector<T>& v1, const Basic3DVector<T>& v2) {
0255 return v1.dot(v2);
0256 }
0257
0258
0259 template <class T, class U>
0260 inline typename PreciseFloatType<T, U>::Type operator*(const Basic3DVector<T>& v1, const Basic3DVector<U>& v2) {
0261 return v1.dot(v2);
0262 }
0263
0264
0265
0266
0267 template <class T>
0268 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, T t) {
0269 return v.v * t;
0270 }
0271
0272
0273 template <class T>
0274 inline Basic3DVector<T> operator*(T t, const Basic3DVector<T>& v) {
0275 return v.v * t;
0276 }
0277
0278 template <class T, typename S>
0279 inline Basic3DVector<T> operator*(S t, const Basic3DVector<T>& v) {
0280 return static_cast<T>(t) * v;
0281 }
0282
0283 template <class T, typename S>
0284 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, S t) {
0285 return static_cast<T>(t) * v;
0286 }
0287
0288
0289
0290
0291 template <class T>
0292 inline Basic3DVector<T> operator/(const Basic3DVector<T>& v, T t) {
0293 return v.v / t;
0294 }
0295
0296 template <class T, typename S>
0297 inline Basic3DVector<T> operator/(const Basic3DVector<T>& v, S s) {
0298
0299 T t = s;
0300 return v / t;
0301 }
0302
0303 typedef Basic3DVector<float> Basic3DVectorF;
0304 typedef Basic3DVector<double> Basic3DVectorD;
0305
0306
0307 #include "Basic3DVectorLD.h"
0308
0309 #endif