interface/private/sseBasic3DVector.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/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 }  // namespace detailsBasic3DVector
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;  // synonym
0038
0039   /** default constructor uses default constructor of T to initialize the
0040    *  components. For built-in floating-point types this means initialization
0041    * to zero??? (force init to 0)
0042    */
0043   Basic3DVector() {}
0044
0045   /// Copy constructor from same type. Should not be needed but for gcc bug 12685
0046   Basic3DVector(const Basic3DVector& p) : v(p.v) {}
0047
0048   /// Assignment operator
0049   Basic3DVector& operator=(const Basic3DVector&) = default;
0050
0051   /// Copy constructor and implicit conversion from Basic3DVector of different precision
0052   template <class U>
0053   Basic3DVector(const Basic3DVector<U>& p) : v(p.v) {}
0054
0055   /// constructor from 2D vector (X and Y from 2D vector, z set to zero)
0056   Basic3DVector(const Basic2DVector<T>& p) : v(p.x(), p.y(), 0) {}
0057
0058   /** Explicit constructor from other (possibly unrelated) vector classes
0059    *  The only constraint on the argument type is that it has methods
0060    *  x(), y() and z(), and that these methods return a type convertible to T.
0061    *  Examples of use are
0062    *   <BR> construction from a Basic3DVector with different precision
0063    *   <BR> construction from a Hep3Vector
0064    *   <BR> construction from a coordinate system converter
0065    */
0066   template <class OtherPoint>
0067   explicit Basic3DVector(const OtherPoint& p) : v(p.x(), p.y(), p.z()) {}
0068
0069   // constructor from Vec4
0070   template <class U>
0071   Basic3DVector(mathSSE::Vec4<U> const& iv) : v(iv) {}
0072
0073   /// construct from cartesian coordinates
0074   Basic3DVector(const T& x, const T& y, const T& z, const T& w = 0) : v(x, y, z, w) {}
0075
0076   /** Deprecated construct from polar coordinates, use
0077    *  <BR> Basic3DVector<T>( Basic3DVector<T>::Polar( theta, phi, r))
0078    *  instead.
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   /// Cartesian x coordinate
0095   T x() const { return v.o.theX; }
0096
0097   /// Cartesian y coordinate
0098   T y() const { return v.o.theY; }
0099
0100   /// Cartesian z coordinate
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   // equality
0108   bool operator==(const Basic3DVector& rh) const { return v == rh.v; }
0109
0110   /// The vector magnitude squared. Equivalent to vec.dot(vec)
0111   T mag2() const { return ::dot(v, v); }
0112
0113   /// The vector magnitude. Equivalent to sqrt(vec.mag2())
0114   T mag() const { return std::sqrt(mag2()); }
0115
0116   /// Squared magnitude of transverse component
0117   T perp2() const { return ::dotxy(v, v); }
0118
0119   /// Magnitude of transverse component
0120   T perp() const { return std::sqrt(perp2()); }
0121
0122   /// Another name for perp()
0123   T transverse() const { return perp(); }
0124
0125   /** Azimuthal angle. The value is returned in radians, in the range (-pi,pi].
0126    *  Same precision as the system atan2(x,y) function.
0127    *  The return type is Geom::Phi<T>, see it's documentation.
0128    */
0129   T barePhi() const { return std::atan2(y(), x()); }
0130   Geom::Phi<T> phi() const { return Geom::Phi<T>(barePhi()); }
0131
0132   /** Polar angle. The value is returned in radians, in the range [0,pi]
0133    *  Same precision as the system atan2(x,y) function.
0134    *  The return type is Geom::Phi<T>, see it's documentation.
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   /** Pseudorapidity.
0140    *  Does not check for zero transverse component; in this case the behavior
0141    *  is as for divide-by zero, i.e. system-dependent.
0142    */
0143   // T eta() const { return -log( tan( theta()/2.));}
0144   T eta() const { return detailsBasic3DVector::eta(x(), y(), z()); }  // correct
0145
0146   /** Unit vector parallel to this.
0147    *  If mag() is zero, a zero vector is returned.
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   /** Operator += with a Basic3DVector of possibly different precision.
0155    */
0156   template <class U>
0157   Basic3DVector& operator+=(const Basic3DVector<U>& p) {
0158     v = v + p.v;
0159     return *this;
0160   }
0161
0162   /** Operator -= with a Basic3DVector of possibly different precision.
0163    */
0164   template <class U>
0165   Basic3DVector& operator-=(const Basic3DVector<U>& p) {
0166     v = v - p.v;
0167     return *this;
0168   }
0169
0170   /// Unary minus, returns a vector with components (-x(),-y(),-z())
0171   Basic3DVector operator-() const { return Basic3DVector(-v); }
0172
0173   /// Scaling by a scalar value (multiplication)
0174   Basic3DVector& operator*=(T t) {
0175     v = t * v;
0176     return *this;
0177   }
0178
0179   /// Scaling by a scalar value (division)
0180   Basic3DVector& operator/=(T t) {
0181     //t = T(1)/t;
0182     v = v / t;
0183     return *this;
0184   }
0185
0186   /// Scalar product, or "dot" product, with a vector of same type.
0187   T dot(const Basic3DVector& rh) const { return ::dot(v, rh.v); }
0188
0189   /** Scalar (or dot) product with a vector of different precision.
0190    *  The product is computed without loss of precision. The type
0191    *  of the returned scalar is the more precise of the scalar types
0192    *  of the two vectors.
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   /// Vector product, or "cross" product, with a vector of same type.
0201   Basic3DVector cross(const Basic3DVector& lh) const { return ::cross(v, lh.v); }
0202
0203   /** Vector (or cross) product with a vector of different precision.
0204    *  The product is computed without loss of precision. The type
0205    *  of the returned vector is the more precise of the types
0206    *  of the two vectors.
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 /// simple text output to standard streams
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 /// vector sum and subtraction of vectors of possibly different precision
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 /// scalar product of vectors of same precision
0253 template <class T>
0254 inline T operator*(const Basic3DVector<T>& v1, const Basic3DVector<T>& v2) {
0255   return v1.dot(v2);
0256 }
0257
0258 /// scalar product of vectors of different precision
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 /** Multiplication by scalar, does not change the precision of the vector.
0265  *  The return type is the same as the type of the vector argument.
0266  */
0267 template <class T>
0268 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, T t) {
0269   return v.v * t;
0270 }
0271
0272 /// Same as operator*( Vector, Scalar)
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 /** Division by scalar, does not change the precision of the vector.
0289  *  The return type is the same as the type of the vector argument.
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   //  T t = S(1)/s; return v*t;
0299   T t = s;
0300   return v / t;
0301 }
0302
0303 typedef Basic3DVector<float> Basic3DVectorF;
0304 typedef Basic3DVector<double> Basic3DVectorD;
0305
0306 //  add long double specialization
0307 #include "Basic3DVectorLD.h"
0308
0309 #endif  // GeometryVector_Basic3DVector_h