Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:53:32

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