ThePEG/Vectors/ThreeVector.h

0001 // -*- C++ -*-
0002 //
0003 // ThreeVector.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 2006-2019 David Grellscheid, Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_ThreeVector_H
0010 #define ThePEG_ThreeVector_H
0011
0012 /**
0013  * @file ThreeVector.h contains the ThreeVector class. ThreeVector can be
0014  * created with any unit type as template parameter. All basic
0015  * mathematical operations are supported, as well as a subset of the
0016  * CLHEP Vector3 functionality.
0017  */
0018
0019 #include "ThreeVector.fh"
0020 #include "ThePEG/Config/ThePEG.h"
0021 #include "ThePEG/Utilities/UnitIO.h"
0022 #include <cassert>
0023 #include <cmath>
0024
0025 namespace ThePEG {
0026
0027 /**
0028  * A 3-component vector. It can be created with any unit type
0029  * as template parameter.  All basic mathematical operations are
0030  * supported, as well as a subset of the CLHEP Vector3
0031  * functionality.
0032  */
0033 template <typename Value>
0034 class ThreeVector
0035 {
0036 private:
0037   /// Value squared
0038   using Value2 = decltype(sqr(std::declval<Value>()));
0039
0040 public:
0041   /** @name Constructors. */
0042   //@{
0043   ThreeVector()
0044     : theX(), theY(), theZ() {}
0045
0046   ThreeVector(Value x, Value y, Value z)
0047     : theX(x), theY(y), theZ(z) {}
0048
0049   template<typename ValueB>
0050   ThreeVector(const ThreeVector<ValueB> & v)
0051     : theX(v.x()), theY(v.y()), theZ(v.z()) {}
0052   //@}
0053
0054 public:
0055   /// @name Component access methods.
0056   //@{
0057   Value x() const { return theX; }
0058   Value y() const { return theY; }
0059   Value z() const { return theZ; }
0060   //@}
0061
0062   /// @name Component set methods.
0063   //@{
0064   void setX(Value x)  {  theX = x; }
0065   void setY(Value y)  {  theY = y; }
0066   void setZ(Value z)  {  theZ = z; }
0067   //@}
0068
0069 public:
0070   /// Squared magnitude \f$x^2+y^2+z^2\f$.
0071   Value2 mag2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
0072
0073   /// Magnitude \f$\sqrt{x^2+y^2+z^2}\f$.
0074   Value  mag() const { return sqrt(mag2()); }
0075
0076   /// Squared transverse component \f$x^2+y^2\f$.
0077   Value2 perp2() const { return sqr(x()) + sqr(y()); }
0078
0079   /// Transverse component \f$\sqrt{x^2+y^2}\f$.
0080   Value  perp()  const { return sqrt(perp2()); }
0081
0082   /// Dot product.
0083   template <typename U>
0084   auto dot(const ThreeVector<U> & a) const
0085   -> decltype(this->x()*a.x())
0086   {
0087     return x()*a.x() + y()*a.y() + z()*a.z();
0088   }
0089
0090   /// Squared transverse component with respect to the given axis.
0091   template <typename U>
0092   Value2 perp2(const ThreeVector<U> & p) const {
0093     const auto pMag2 = p.mag2();
0094     assert( pMag2 > ZERO );
0095     auto ss = this->dot(p);
0096     Value2 ret = mag2() - sqr(ss)/pMag2;
0097     if ( ret <= ZERO )
0098       ret = ZERO;
0099     return ret;
0100   }
0101
0102   /// Transverse component with respect to the given axis.
0103   template <typename U>
0104   Value perp(const ThreeVector<U> & p) const {
0105     return sqrt(perp2(p));
0106   }
0107
0108   /// @name Spherical coordinates.
0109   //@{
0110   /// Polar angle.
0111   double theta() const {
0112     assert(!(x() == ZERO && y() == ZERO && z() == ZERO));
0113     return atan2(perp(),z());
0114   }
0115
0116   /// Azimuthal angle.
0117   double phi()   const {
0118     return atan2(y(),x());
0119   }
0120
0121   /// Set the polar angle.
0122   void setTheta(double th) {
0123     double ma  = mag();
0124     double ph  = phi();
0125     setX(ma*sin(th)*cos(ph));
0126     setY(ma*sin(th)*sin(ph));
0127     setZ(ma*cos(th));
0128   }
0129
0130   /// Set the azimuthal angle.
0131   void setPhi(double ph) {
0132     double xy = perp();
0133     setX(xy*cos(ph));
0134     setY(xy*sin(ph));
0135   }
0136   //@}
0137
0138   /// Parallel vector with unit length.
0139   ThreeVector<double> unit() const {
0140     Value2 mg2 = mag2();
0141     assert(mg2 > ZERO);
0142     Value mg = sqrt(mg2);
0143     return {x()/mg, y()/mg, z()/mg};
0144   }
0145
0146   /// Orthogonal vector.
0147   ThreeVector<Value> orthogonal() const {
0148     Value xx = abs(x());
0149     Value yy = abs(y());
0150     Value zz = abs(z());
0151     using TVec = ThreeVector<Value>;
0152     if (xx < yy) {
0153       return xx < zz ? TVec{ZERO,z(),-y()} : TVec{y(),-x(),ZERO};
0154     } else {
0155       return yy < zz ? TVec{-z(),ZERO,x()} : TVec{y(),-x(),ZERO};
0156     }
0157   }
0158
0159   /// Azimuthal angle difference, brought into the range \f$(-\pi,\pi]\f$.
0160   template <typename U>
0161   double deltaPhi  (const ThreeVector<U> & v2) const {
0162     double dphi = v2.phi() - phi();
0163     if ( dphi > Constants::pi ) {
0164       dphi -= Constants::twopi;
0165     } else if ( dphi <= -Constants::pi ) {
0166       dphi += Constants::twopi;
0167     }
0168     return dphi;
0169   }
0170
0171   /**
0172    * Apply a rotation.
0173    * @param angle Rotation angle in radians.
0174    * @param axis Rotation axis.
0175    */
0176   template <typename U>
0177   ThreeVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
0178     if (angle == 0.0)
0179       return *this;
0180     const U ll = axis.mag();
0181     assert( ll > ZERO );
0182
0183     const double sa = sin(angle), ca = cos(angle);
0184     const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
0185     const Value  xx  = x(), yy = y(), zz = z();
0186
0187     setX((ca+(1-ca)*dx*dx)     * xx
0188          +((1-ca)*dx*dy-sa*dz) * yy
0189          +((1-ca)*dx*dz+sa*dy) * zz
0190          );
0191     setY(((1-ca)*dy*dx+sa*dz)  * xx
0192          +(ca+(1-ca)*dy*dy)    * yy
0193          +((1-ca)*dy*dz-sa*dx) * zz
0194          );
0195     setZ(((1-ca)*dz*dx-sa*dy)  * xx
0196          +((1-ca)*dz*dy+sa*dx) * yy
0197          +(ca+(1-ca)*dz*dz)    * zz
0198          );
0199     return *this;
0200   }
0201
0202
0203   /**
0204    * Rotate the reference frame to a new z-axis.
0205    */
0206   ThreeVector<Value> & rotateUz (const Axis & axis) {
0207     Axis ax = axis.unit();
0208     double u1 = ax.x();
0209     double u2 = ax.y();
0210     double u3 = ax.z();
0211     double up = u1*u1 + u2*u2;
0212     if (up>0) {
0213       up = sqrt(up);
0214       Value px = x(),  py = y(),  pz = z();
0215       setX( (u1*u3*px - u2*py)/up + u1*pz );
0216       setY( (u2*u3*px + u1*py)/up + u2*pz );
0217       setZ(    -up*px +             u3*pz );
0218     }
0219     else if (u3 < 0.) {
0220       setX(-x());
0221       setZ(-z());
0222     }
0223     return *this;
0224   }
0225
0226   /**
0227    * Rotate from a reference frame to the z-axis.
0228    */
0229   ThreeVector<Value> & rotateUzBack (const Axis & axis) {
0230     Axis ax = axis.unit();
0231     double u1 = ax.x();
0232     double u2 = ax.y();
0233     double u3 = ax.z();
0234     double up = u1*u1 + u2*u2;
0235     if (up>0) {
0236       up = sqrt(up);
0237       Value px = x(),  py = y(),  pz = z();
0238       setX( ( u1*u3*px + u2*u3*py)/up - up*pz );
0239       setY( (-u2*px    + u1*py)/up );
0240       setZ(   u1*px    + u2*py        + u3*pz );
0241     }
0242     else if (u3 < 0.) {
0243       setX(-x());
0244       setZ(-z());
0245     }
0246     return *this;
0247   }
0248
0249   /// Vector cross-product
0250   template <typename U>
0251   auto cross(const ThreeVector<U> & a) const
0252   -> ThreeVector<decltype(this->y()*a.z())>
0253   {
0254     return { y()*a.z()-z()*a.y(),
0255             -x()*a.z()+z()*a.x(),
0256              x()*a.y()-y()*a.x() };
0257   }
0258
0259   public:
0260   /// @name Comparison operators.
0261   //@{
0262   bool operator==(const ThreeVector<Value> & a) const {
0263     return (theX == a.x() && theY == a.y() && theZ == a.z());
0264   }
0265   bool operator!=(const ThreeVector<Value> & a) const {
0266     return !(*this == a);
0267   }
0268   bool almostEqual(const ThreeVector<Value> & a, double threshold = 1e-04) const {
0269     return ((std::abs(theX - a.x()) < threshold) && (std::abs(theY - a.y()) < threshold) && (std::abs(theZ - a.z()) < threshold));
0270   }
0271   bool almostUnequal(const ThreeVector<Value> & a, double threshold = 1e-04) const {
0272     return ! this->almostEqual(a, threshold);
0273   }
0274      //@}
0275
0276 public:
0277   /// @name Mathematical assignment operators.
0278   //@{
0279   ThreeVector<Value> & operator+=(const ThreeVector<Value> & a) {
0280     theX += a.x();
0281     theY += a.y();
0282     theZ += a.z();
0283     return *this;
0284   }
0285
0286   ThreeVector<Value> & operator-=(const ThreeVector<Value> & a) {
0287     theX -= a.x();
0288     theY -= a.y();
0289     theZ -= a.z();
0290     return *this;
0291   }
0292
0293   ThreeVector<Value> & operator*=(double a) {
0294     theX *= a;
0295     theY *= a;
0296     theZ *= a;
0297     return *this;
0298   }
0299
0300   ThreeVector<Value> & operator/=(double a) {
0301     theX /= a;
0302     theY /= a;
0303     theZ /= a;
0304     return *this;
0305   }
0306   //@}
0307
0308   /// Cosine of the azimuthal angle between two vectors.
0309   template <typename U>
0310   double cosTheta(const ThreeVector<U> & q) const {
0311     auto ptot = mag()*q.mag();
0312     assert( ptot > ZERO );
0313     double arg = dot(q)/ptot;
0314     if     (arg >  1.0) arg =  1.0;
0315     else if(arg < -1.0) arg = -1.0;
0316     return arg;
0317   }
0318
0319   /// Angle between two vectors.
0320   template <typename U>
0321   double angle(const ThreeVector<U> & v) const {
0322     return acos(cosTheta(v));
0323   }
0324
0325 private:
0326   /// @name Vector components
0327   //@{
0328   Value theX;
0329   Value theY;
0330   Value theZ;
0331   //@}
0332 };
0333
0334 /// Stream output. Format \f$(x,y,z)\f$.
0335 inline ostream &
0336 operator<< (ostream & os, const ThreeVector<double> & v)
0337 {
0338   return os << '(' << v.x() << ',' << v.y() << ',' << v.z() << ')';
0339 }
0340
0341 /// @name Basic mathematical operations
0342 //@{
0343 template <typename Value>
0344 inline ThreeVector<Value>
0345 operator+(ThreeVector<Value> a,
0346           const ThreeVector<Value> & b)
0347 {
0348   return a += b;
0349 }
0350
0351 template <typename Value>
0352 inline ThreeVector<Value>
0353 operator-(ThreeVector<Value> a,
0354           const ThreeVector<Value> & b)
0355 {
0356   return a -= b;
0357 }
0358
0359 template <typename Value>
0360 inline ThreeVector<Value> operator-(const ThreeVector<Value> & v) {
0361   return {-v.x(),-v.y(),-v.z()};
0362 }
0363
0364 template <typename Value>
0365 inline ThreeVector<Value> operator*(ThreeVector<Value> v, double a) {
0366   return v *= a;
0367 }
0368
0369 template <typename Value>
0370 inline ThreeVector<Value> operator*(double a, ThreeVector<Value> v) {
0371   return v *= a;
0372 }
0373
0374 template <typename ValueA, typename ValueB>
0375 inline auto operator*(ValueB a, ThreeVector<ValueA> v)
0376 -> ThreeVector<decltype(a*v.x())>
0377 {
0378   return {a*v.x(), a*v.y(), a*v.z()};
0379 }
0380
0381 template <typename ValueA, typename ValueB>
0382 inline auto operator*(ThreeVector<ValueA> v, ValueB a)
0383 -> ThreeVector<decltype(v.x()*a)>
0384 {
0385   return {v.x()*a, v.y()*a, v.z()*a};
0386 }
0387 //@}
0388
0389 /// Vector dot product.
0390 template <typename ValueA, typename ValueB>
0391 inline auto operator*(const ThreeVector<ValueA> & a,
0392                       const ThreeVector<ValueB> & b)
0393 -> decltype(a.x()*b.x())
0394 {
0395   return a.dot(b);
0396 }
0397
0398 /// A parallel vector with unit length.
0399 template <typename Value>
0400 ThreeVector<double> unitVector(const ThreeVector<Value> & v) {
0401   return v.unit();
0402 }
0403
0404
0405 /** Output a ThreeVector with units to a stream. */
0406 template <typename OStream, typename UT, typename Value>
0407 void ounitstream(OStream & os, const ThreeVector<Value> & p, UT & u) {
0408   os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u);
0409 }
0410
0411 /** Input a ThreeVector with units from a stream. */
0412 template <typename IStream, typename UT, typename Value>
0413 void iunitstream(IStream & is, ThreeVector<Value> & p, UT & u) {
0414   Value x, y, z;
0415   is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u);
0416   p = ThreeVector<Value>(x, y, z);
0417 }
0418
0419 }
0420
0421 #endif /* ThePEG_ThreeVector_H */