interface/private/oldBasic3DVector.h

0001 #ifndef GeometryVector_oldBasic3DVector_h
0002 #define GeometryVector_oldBasic3DVector_h
0003 #if ( defined(__CLING__) || defined(__CINT__) )  && !defined(__REFLEX__)
0004 #define __REFLEX__
0005 #endif
0006 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
0007 #include "DataFormats/GeometryVector/interface/Theta.h"
0008 #include "DataFormats/GeometryVector/interface/Phi.h"
0009 #include "DataFormats/GeometryVector/interface/PreciseFloatType.h"
0010 #include "DataFormats/GeometryVector/interface/CoordinateSets.h"
0011 #ifndef __REFLEX__
0012 #include "DataFormats/Math/interface/SIMDVec.h"
0013 #endif
0014 #include <iosfwd>
0015 #include <cmath>
0016
0017
0018 namespace detailsBasic3DVector {
0019   inline float __attribute__((always_inline)) __attribute__ ((pure))
0020   eta(float x, float y, float z) { float t(z/std::sqrt(x*x+y*y)); return ::asinhf(t);}
0021   inline double __attribute__((always_inline)) __attribute__ ((pure))
0022   eta(double x, double y, double z) { double t(z/std::sqrt(x*x+y*y)); return ::asinh(t);}
0023   inline long double __attribute__((always_inline)) __attribute__ ((pure))
0024   eta(long double x, long double y, long double z) { long double t(z/std::sqrt(x*x+y*y)); return ::asinhl(t);}
0025 }
0026
0027
0028 template < typename T>
0029 class Basic3DVector {
0030 public:
0031   typedef  Basic3DVector<T> MathVector;
0032
0033
0034   typedef T                                   ScalarType;
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() : theX(0), theY(0), theZ(0), theW(0) {}
0044
0045   /// Copy constructor from same type. Should not be needed but for gcc bug 12685
0046   Basic3DVector( const Basic3DVector & p) :
0047     theX(p.x()), theY(p.y()), theZ(p.z()), theW(p.w()) {}
0048
0049   /// Copy constructor and implicit conversion from Basic3DVector of different precision
0050   template <class U>
0051   Basic3DVector( const Basic3DVector<U> & p) :
0052     theX(p.x()), theY(p.y()), theZ(p.z()), theW(p.w()) {}
0053
0054   /// constructor from 2D vector (X and Y from 2D vector, z set to zero)
0055   Basic3DVector( const Basic2DVector<T> & p) :
0056     theX(p.x()), theY(p.y()), theZ(0), theW(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) :
0068     theX(p.x()), theY(p.y()), theZ(p.z()), theW(0) {}
0069
0070
0071 #if  defined(USE_EXTVECT)
0072   template<typename U>
0073   Basic3DVector(Vec4<U> const& iv) :
0074     theX(iv[0]), theY(iv[1]), theZ(iv[2]), theW(0) {}
0075 #elif  defined(USE_SSEVECT)
0076   // constructor from Vec4
0077   template<typename U>
0078   Basic3DVector(mathSSE::Vec4<U> const& iv) :
0079     theX(iv.arr[0]), theY(iv.arr[1]), theZ(iv.arr[2]), theW(0) {}
0080 #endif
0081
0082 #ifndef __REFLEX__
0083   /// construct from cartesian coordinates
0084   Basic3DVector( const T& x, const T& y, const T& z, const T& w=0) :
0085     theX(x), theY(y), theZ(z), theW(w) {}
0086 #else
0087   /// construct from cartesian coordinates
0088   Basic3DVector( const T& x, const T& y, const T& z) :
0089     theX(x), theY(y), theZ(z), theW(0) {}
0090   Basic3DVector( const T& x, const T& y, const T& z, const T& w) :
0091     theX(x), theY(y), theZ(z), theW(w) {}
0092 #endif
0093
0094
0095
0096   /** Deprecated construct from polar coordinates, use
0097    *  <BR> Basic3DVector<T>( Basic3DVector<T>::Polar( theta, phi, r))
0098    *  instead.
0099    */
0100   template <typename U>
0101   Basic3DVector( const Geom::Theta<U>& theta,
0102          const Geom::Phi<U>& phi, const T& r) {
0103     Polar p( theta.value(), phi.value(), r);
0104     theX = p.x(); theY = p.y(); theZ = p.z();
0105   }
0106
0107
0108   T operator[](int i) const { return *((&theX)+i)  ;}
0109   T & operator[](int i) { return *((&theX)+i);}
0110
0111
0112
0113   /// Cartesian x coordinate
0114   T x() const { return theX;}
0115
0116   /// Cartesian y coordinate
0117   T y() const { return theY;}
0118
0119   /// Cartesian z coordinate
0120   T z() const { return theZ;}
0121
0122   T w() const { return theW;}
0123
0124
0125   Basic2DVector<T> xy() const { return  Basic2DVector<T>(theX,theY);}
0126
0127
0128   // equality
0129   bool operator==(const Basic3DVector& rh) const {
0130     return x()==rh.x() && y()==rh.y() && z()==rh.z();
0131   }
0132
0133   /// The vector magnitude squared. Equivalent to vec.dot(vec)
0134   T mag2() const { return  x()*x() + y()*y()+z()*z();}
0135
0136   /// The vector magnitude. Equivalent to sqrt(vec.mag2())
0137   T mag() const  { return std::sqrt( mag2());}
0138
0139   /// Squared magnitude of transverse component
0140   T perp2() const { return x()*x() + y()*y();}
0141
0142   /// Magnitude of transverse component
0143   T perp() const { return std::sqrt( perp2());}
0144
0145   /// Another name for perp()
0146   T transverse() const { return perp();}
0147
0148   /** Azimuthal angle. The value is returned in radians, in the range (-pi,pi].
0149    *  Same precision as the system atan2(x,y) function.
0150    *  The return type is Geom::Phi<T>, see it's documentation.
0151    */
0152   T barePhi() const {return std::atan2(y(),x());}
0153   Geom::Phi<T> phi() const {return Geom::Phi<T>(barePhi());}
0154
0155   /** Polar angle. The value is returned in radians, in the range [0,pi]
0156    *  Same precision as the system atan2(x,y) function.
0157    *  The return type is Geom::Phi<T>, see it's documentation.
0158    */
0159   T bareTheta() const {return std::atan2(perp(),z());}
0160   Geom::Theta<T> theta() const {return Geom::Theta<T>(std::atan2(perp(),z()));}
0161
0162   /** Pseudorapidity.
0163    *  Does not check for zero transverse component; in this case the behavior
0164    *  is as for divide-by zero, i.e. system-dependent.
0165    */
0166   // T eta() const { return -log( tan( theta()/2.));}
0167   T eta() const { return detailsBasic3DVector::eta(x(),y(),z());} // correct
0168
0169   /** Unit vector parallel to this.
0170    *  If mag() is zero, a zero vector is returned.
0171    */
0172   Basic3DVector unit() const {
0173     T my_mag = mag2();
0174     if (my_mag==0) return *this;
0175     my_mag = T(1)/std::sqrt(my_mag);
0176     return *this * my_mag;
0177   }
0178
0179   /** Operator += with a Basic3DVector of possibly different precision.
0180    */
0181   template <class U>
0182   Basic3DVector& operator+= ( const Basic3DVector<U>& p) {
0183     theX += p.x();
0184     theY += p.y();
0185     theZ += p.z();
0186     theW += p.w();
0187     return *this;
0188   }
0189
0190   /** Operator -= with a Basic3DVector of possibly different precision.
0191    */
0192   template <class U>
0193   Basic3DVector& operator-= ( const Basic3DVector<U>& p) {
0194     theX -= p.x();
0195     theY -= p.y();
0196     theZ -= p.z();
0197     theW -= p.w();
0198     return *this;
0199   }
0200
0201   /// Unary minus, returns a vector with components (-x(),-y(),-z())
0202   Basic3DVector operator-() const { return Basic3DVector(-x(),-y(),-z());}
0203
0204   /// Scaling by a scalar value (multiplication)
0205   Basic3DVector& operator*= ( T t) {
0206     theX *= t;
0207     theY *= t;
0208     theZ *= t;
0209     theW *= t;;
0210     return *this;
0211   }
0212
0213   /// Scaling by a scalar value (division)
0214   Basic3DVector& operator/= ( T t) {
0215     t = T(1)/t;
0216     theX *= t;
0217     theY *= t;
0218     theZ *= t;
0219     theW *= t;;
0220     return *this;
0221   }
0222
0223   /// Scalar product, or "dot" product, with a vector of same type.
0224   T dot( const Basic3DVector& v) const {
0225     return x()*v.x() + y()*v.y() + z()*v.z();
0226   }
0227
0228   /** Scalar (or dot) product with a vector of different precision.
0229    *  The product is computed without loss of precision. The type
0230    *  of the returned scalar is the more precise of the scalar types
0231    *  of the two vectors.
0232    */
0233   template <class U>
0234   typename PreciseFloatType<T,U>::Type dot( const Basic3DVector<U>& v) const {
0235     return x()*v.x() + y()*v.y() + z()*v.z();
0236   }
0237
0238   /// Vector product, or "cross" product, with a vector of same type.
0239   Basic3DVector cross( const Basic3DVector& v) const {
0240     return Basic3DVector( y()*v.z() - v.y()*z(),
0241               z()*v.x() - v.z()*x(),
0242               x()*v.y() - v.x()*y());
0243   }
0244
0245
0246   /** Vector (or cross) product with a vector of different precision.
0247    *  The product is computed without loss of precision. The type
0248    *  of the returned vector is the more precise of the types
0249    *  of the two vectors.
0250    */
0251   template <class U>
0252   Basic3DVector<typename PreciseFloatType<T,U>::Type>
0253   cross( const Basic3DVector<U>& v) const {
0254     return Basic3DVector<typename PreciseFloatType<T,U>::Type>( y()*v.z() - v.y()*z(),
0255                                 z()*v.x() - v.z()*x(),
0256                                 x()*v.y() - v.x()*y());
0257   }
0258
0259 private:
0260   T theX;
0261   T theY;
0262   T theZ;
0263   T theW;
0264 }
0265 #ifndef __CINT__
0266 __attribute__ ((aligned (16)))
0267 #endif
0268 ;
0269
0270
0271 namespace geometryDetails {
0272   std::ostream & print3D(std::ostream& s, double x, double y, double z);
0273 }
0274
0275 /// simple text output to standard streams
0276 template <class T>
0277 inline std::ostream & operator<<( std::ostream& s, const Basic3DVector<T>& v) {
0278   return geometryDetails::print3D(s, v.x(),v.y(), v.z());
0279 }
0280
0281
0282 /// vector sum and subtraction of vectors of possibly different precision
0283 template <class T, class U>
0284 inline Basic3DVector<typename PreciseFloatType<T,U>::Type>
0285 operator+( const Basic3DVector<T>& a, const Basic3DVector<U>& b) {
0286   typedef Basic3DVector<typename PreciseFloatType<T,U>::Type> RT;
0287   return RT(a.x()+b.x(), a.y()+b.y(), a.z()+b.z(), a.w()+b.w());
0288 }
0289
0290 template <class T, class U>
0291 inline Basic3DVector<typename PreciseFloatType<T,U>::Type>
0292 operator-( const Basic3DVector<T>& a, const Basic3DVector<U>& b) {
0293   typedef Basic3DVector<typename PreciseFloatType<T,U>::Type> RT;
0294   return RT(a.x()-b.x(), a.y()-b.y(), a.z()-b.z(), a.w()-b.w());
0295 }
0296
0297 /// scalar product of vectors of same precision
0298 template <class T>
0299 inline T operator*( const Basic3DVector<T>& v1, const Basic3DVector<T>& v2) {
0300   return v1.dot(v2);
0301 }
0302
0303 /// scalar product of vectors of different precision
0304 template <class T, class U>
0305 inline typename PreciseFloatType<T,U>::Type operator*( const Basic3DVector<T>& v1,
0306                                const Basic3DVector<U>& v2) {
0307   return v1.x()*v2.x() + v1.y()*v2.y() + v1.z()*v2.z();
0308 }
0309
0310 /** Multiplication by scalar, does not change the precision of the vector.
0311  *  The return type is the same as the type of the vector argument.
0312  */
0313 template <class T>
0314 inline Basic3DVector<T> operator*( const Basic3DVector<T>& v, T t) {
0315   return Basic3DVector<T>(v.x()*t, v.y()*t, v.z()*t, v.w()*t);
0316 }
0317
0318 /// Same as operator*( Vector, Scalar)
0319 template <class T>
0320 inline Basic3DVector<T> operator*(T t, const Basic3DVector<T>& v) {
0321   return Basic3DVector<T>(v.x()*t, v.y()*t, v.z()*t, v.w()*t);
0322 }
0323
0324 template <class T, typename S>
0325 inline Basic3DVector<T> operator*(S t,  const Basic3DVector<T>& v) {
0326   return static_cast<T>(t)*v;
0327 }
0328
0329 template <class T, typename S>
0330 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, S t) {
0331   return static_cast<T>(t)*v;
0332 }
0333
0334
0335 /** Division by scalar, does not change the precision of the vector.
0336  *  The return type is the same as the type of the vector argument.
0337  */
0338 template <class T, typename S>
0339 inline Basic3DVector<T> operator/( const Basic3DVector<T>& v, S s) {
0340   T t = T(1)/s;
0341   return v*t;
0342 }
0343
0344
0345 typedef Basic3DVector<float> Basic3DVectorF;
0346 typedef Basic3DVector<double> Basic3DVectorD;
0347 typedef Basic3DVector<long double> Basic3DVectorLD;
0348
0349
0350 #endif // GeometryVector_Basic3DVector_h
0351
0352