interface/private/extBasic3DVector.h

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 "FWCore/Utilities/interface/Likely.h"
0011 #include <iosfwd>
0012 #include <cmath>
0013
0014 namespace detailsBasic3DVector {
0015   inline float __attribute__((always_inline)) __attribute__((pure)) eta(float x, float y, float z) {
0016     float t(z / std::sqrt(x * x + y * y));
0017     return ::asinhf(t);
0018   }
0019   inline double __attribute__((always_inline)) __attribute__((pure)) eta(double x, double y, double z) {
0020     double t(z / std::sqrt(x * x + y * y));
0021     return ::asinh(t);
0022   }
0023   inline long double __attribute__((always_inline)) __attribute__((pure)) eta(long double x,
0024                                                                               long double y,
0025                                                                               long double z) {
0026     long double t(z / std::sqrt(x * x + y * y));
0027     return ::asinhl(t);
0028   }
0029 }  // namespace detailsBasic3DVector
0030
0031 template <typename T>
0032 class Basic3DVector {
0033 public:
0034   typedef T ScalarType;
0035   typedef Vec4<T> VectorType;
0036   typedef Vec4<T> MathVector;
0037   typedef Geom::Cylindrical2Cartesian<T> Cylindrical;
0038   typedef Geom::Spherical2Cartesian<T> Spherical;
0039   typedef Spherical Polar;  // synonym
0040
0041   /** default constructor uses default constructor of T to initialize the
0042    *  components. For built-in floating-point types this means initialization
0043    * to zero??? (force init to 0)
0044    */
0045   Basic3DVector() : v{0, 0, 0, 0} {}
0046
0047   /// Copy constructor from same type. Should not be needed but for gcc bug 12685
0048   Basic3DVector(const Basic3DVector& p) : v(p.v) {}
0049
0050   /// Assignment operator
0051   Basic3DVector& operator=(const Basic3DVector&) = default;
0052
0053   /// Copy constructor and implicit conversion from Basic3DVector of different precision
0054   template <class U>
0055   Basic3DVector(const Basic3DVector<U>& p) : v{T(p.v[0]), T(p.v[1]), T(p.v[2]), T(p.v[3])} {}
0056
0057   /// constructor from 2D vector (X and Y from 2D vector, z set to zero)
0058   Basic3DVector(const Basic2DVector<T>& p) : v{p.x(), p.y(), 0} {}
0059
0060   /** Explicit constructor from other (possibly unrelated) vector classes
0061    *  The only constraint on the argument type is that it has methods
0062    *  x(), y() and z(), and that these methods return a type convertible to T.
0063    *  Examples of use are
0064    *   <BR> construction from a Basic3DVector with different precision
0065    *   <BR> construction from a Hep3Vector
0066    *   <BR> construction from a coordinate system converter
0067    */
0068   template <class OtherPoint>
0069   explicit Basic3DVector(const OtherPoint& p) : v{T(p.x()), T(p.y()), T(p.z())} {}
0070
0071   // constructor from Vec4
0072   Basic3DVector(MathVector const& iv) : v(iv) {}
0073
0074   template <class U>
0075   Basic3DVector(Vec4<U> const& iv) : v{T(iv[0]), T(iv[1]), T(iv[2]), T(iv[3])} {}
0076
0077   /// construct from cartesian coordinates
0078   Basic3DVector(const T& x, const T& y, const T& z, const T& w = 0) : v{x, y, z, w} {}
0079
0080   /** Deprecated construct from polar coordinates, use
0081    *  <BR> Basic3DVector<T>( Basic3DVector<T>::Polar( theta, phi, r))
0082    *  instead.
0083    */
0084   template <typename U>
0085   Basic3DVector(const Geom::Theta<U>& theta, const Geom::Phi<U>& phi, const T& r) {
0086     Polar p(theta.value(), phi.value(), r);
0087     v[0] = p.x();
0088     v[1] = p.y();
0089     v[2] = p.z();
0090   }
0091
0092   MathVector const& mathVector() const { return v; }
0093   MathVector& mathVector() { return v; }
0094
0095   T operator[](int i) const { return v[i]; }
0096
0097   /// Cartesian x coordinate
0098   T x() const { return v[0]; }
0099
0100   /// Cartesian y coordinate
0101   T y() const { return v[1]; }
0102
0103   /// Cartesian z coordinate
0104   T z() const { return v[2]; }
0105
0106   T w() const { return v[3]; }
0107
0108   Basic2DVector<T> xy() const { return ::xy(v); }
0109
0110   // equality
0111   bool operator==(const Basic3DVector& rh) const {
0112     auto res = v == rh.v;
0113     return res[0] & res[1] & res[2] & res[3];
0114   }
0115
0116   /// The vector magnitude squared. Equivalent to vec.dot(vec)
0117   T mag2() const { return ::dot(v, v); }
0118
0119   /// The vector magnitude. Equivalent to sqrt(vec.mag2())
0120   T mag() const { return std::sqrt(mag2()); }
0121
0122   /// Squared magnitude of transverse component
0123   T perp2() const { return ::dot2(v, v); }
0124
0125   /// Magnitude of transverse component
0126   T perp() const { return std::sqrt(perp2()); }
0127
0128   /// Another name for perp()
0129   T transverse() const { return perp(); }
0130
0131   /** Azimuthal angle. The value is returned in radians, in the range (-pi,pi].
0132    *  Same precision as the system atan2(x,y) function.
0133    *  The return type is Geom::Phi<T>, see it's documentation.
0134    */
0135   T barePhi() const { return std::atan2(y(), x()); }
0136   Geom::Phi<T> phi() const { return Geom::Phi<T>(barePhi()); }
0137
0138   /** Polar angle. The value is returned in radians, in the range [0,pi]
0139    *  Same precision as the system atan2(x,y) function.
0140    *  The return type is Geom::Phi<T>, see it's documentation.
0141    */
0142   T bareTheta() const { return std::atan2(perp(), z()); }
0143   Geom::Theta<T> theta() const { return Geom::Theta<T>(std::atan2(perp(), z())); }
0144
0145   /** Pseudorapidity.
0146    *  Does not check for zero transverse component; in this case the behavior
0147    *  is as for divide-by zero, i.e. system-dependent.
0148    */
0149   // T eta() const { return -log( tan( theta()/2.));}
0150   T eta() const { return detailsBasic3DVector::eta(x(), y(), z()); }  // correct
0151
0152   /** Unit vector parallel to this.
0153    *  If mag() is zero, a zero vector is returned.
0154    */
0155   Basic3DVector unit() const {
0156     T my_mag = mag2();
0157     return LIKELY(0 != my_mag) ? (*this) * (T(1) / std::sqrt(my_mag)) : *this;
0158   }
0159
0160   /** Operator += with a Basic3DVector of possibly different precision.
0161    */
0162   template <class U>
0163   Basic3DVector& operator+=(const Basic3DVector<U>& p) {
0164     v = v + p.v;
0165     return *this;
0166   }
0167
0168   /** Operator -= with a Basic3DVector of possibly different precision.
0169    */
0170   template <class U>
0171   Basic3DVector& operator-=(const Basic3DVector<U>& p) {
0172     v = v - p.v;
0173     return *this;
0174   }
0175
0176   /// Unary minus, returns a vector with components (-x(),-y(),-z())
0177   Basic3DVector operator-() const { return Basic3DVector(-v); }
0178
0179   /// Scaling by a scalar value (multiplication)
0180   Basic3DVector& operator*=(T t) {
0181     v = t * v;
0182     return *this;
0183   }
0184
0185   /// Scaling by a scalar value (division)
0186   Basic3DVector& operator/=(T t) {
0187     //t = T(1)/t;
0188     v = v / t;
0189     return *this;
0190   }
0191
0192   /// Scalar product, or "dot" product, with a vector of same type.
0193   T dot(const Basic3DVector& rh) const { return ::dot(v, rh.v); }
0194
0195   /** Scalar (or dot) product with a vector of different precision.
0196    *  The product is computed without loss of precision. The type
0197    *  of the returned scalar is the more precise of the scalar types
0198    *  of the two vectors.
0199    */
0200   template <class U>
0201   typename PreciseFloatType<T, U>::Type dot(const Basic3DVector<U>& lh) const {
0202     return Basic3DVector<typename PreciseFloatType<T, U>::Type>(*this).dot(
0203         Basic3DVector<typename PreciseFloatType<T, U>::Type>(lh));
0204   }
0205
0206   /// Vector product, or "cross" product, with a vector of same type.
0207   Basic3DVector cross(const Basic3DVector& lh) const { return ::cross3(v, lh.v); }
0208
0209   /** Vector (or cross) product with a vector of different precision.
0210    *  The product is computed without loss of precision. The type
0211    *  of the returned vector is the more precise of the types
0212    *  of the two vectors.
0213    */
0214   template <class U>
0215   Basic3DVector<typename PreciseFloatType<T, U>::Type> cross(const Basic3DVector<U>& lh) const {
0216     return Basic3DVector<typename PreciseFloatType<T, U>::Type>(*this).cross(
0217         Basic3DVector<typename PreciseFloatType<T, U>::Type>(lh));
0218   }
0219
0220 public:
0221   Vec4<T> v;
0222 } __attribute__((aligned(16)));
0223
0224 namespace geometryDetails {
0225   std::ostream& print3D(std::ostream& s, double x, double y, double z);
0226 }
0227
0228 /// simple text output to standard streams
0229 template <class T>
0230 inline std::ostream& operator<<(std::ostream& s, const Basic3DVector<T>& v) {
0231   return geometryDetails::print3D(s, v.x(), v.y(), v.z());
0232 }
0233
0234 /// vector sum and subtraction of vectors of possibly different precision
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 template <class T>
0240 inline Basic3DVector<T> operator-(const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
0241   return a.v - b.v;
0242 }
0243
0244 template <class T, class U>
0245 inline Basic3DVector<typename PreciseFloatType<T, U>::Type> operator+(const Basic3DVector<T>& a,
0246                                                                       const Basic3DVector<U>& b) {
0247   typedef Basic3DVector<typename PreciseFloatType<T, U>::Type> RT;
0248   return RT(a).v + RT(b).v;
0249 }
0250
0251 template <class T, class U>
0252 inline Basic3DVector<typename PreciseFloatType<T, U>::Type> operator-(const Basic3DVector<T>& a,
0253                                                                       const Basic3DVector<U>& b) {
0254   typedef Basic3DVector<typename PreciseFloatType<T, U>::Type> RT;
0255   return RT(a).v - RT(b).v;
0256 }
0257
0258 /// scalar product of vectors of same precision
0259 template <class T>
0260 inline T operator*(const Basic3DVector<T>& v1, const Basic3DVector<T>& v2) {
0261   return v1.dot(v2);
0262 }
0263
0264 /// scalar product of vectors of different precision
0265 template <class T, class U>
0266 inline typename PreciseFloatType<T, U>::Type operator*(const Basic3DVector<T>& v1, const Basic3DVector<U>& v2) {
0267   return v1.dot(v2);
0268 }
0269
0270 /** Multiplication by scalar, does not change the precision of the vector.
0271  *  The return type is the same as the type of the vector argument.
0272  */
0273 template <class T>
0274 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, T t) {
0275   return v.v * t;
0276 }
0277
0278 /// Same as operator*( Vector, Scalar)
0279 template <class T>
0280 inline Basic3DVector<T> operator*(T t, const Basic3DVector<T>& v) {
0281   return v.v * t;
0282 }
0283
0284 template <class T, typename S>
0285 inline Basic3DVector<T> operator*(S t, const Basic3DVector<T>& v) {
0286   return static_cast<T>(t) * v;
0287 }
0288
0289 template <class T, typename S>
0290 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, S t) {
0291   return static_cast<T>(t) * v;
0292 }
0293
0294 /** Division by scalar, does not change the precision of the vector.
0295  *  The return type is the same as the type of the vector argument.
0296  */
0297 template <class T>
0298 inline Basic3DVector<T> operator/(const Basic3DVector<T>& v, T t) {
0299   return v.v / t;
0300 }
0301
0302 template <class T, typename S>
0303 inline Basic3DVector<T> operator/(const Basic3DVector<T>& v, S s) {
0304   //  T t = S(1)/s; return v*t;
0305   T t = s;
0306   return v / t;
0307 }
0308
0309 typedef Basic3DVector<float> Basic3DVectorF;
0310 typedef Basic3DVector<double> Basic3DVectorD;
0311
0312 //  add long double specialization
0313 #include "Basic3DVectorLD.h"
0314
0315 #endif  // GeometryVector_Basic3DVector_h