diff --git a/Vectors/LorentzVector.h b/Vectors/LorentzVector.h --- a/Vectors/LorentzVector.h +++ b/Vectors/LorentzVector.h @@ -1,765 +1,765 @@ // -*- C++ -*- // // LorentzVector.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_LorentzVector_H #define ThePEG_LorentzVector_H /** * @file LorentzVector.h contains the LorentzVector class. Lorentz * vectors can be created with any unit type as template parameter. * All basic mathematical operations are supported, as well as a * subset of the CLHEP LorentzVector functionality. */ #include "LorentzVector.fh" #include "ThePEG/Utilities/Direction.h" #include "ThePEG/Utilities/UnitIO.h" #include "LorentzRotation.h" #include "ThreeVector.h" /// Debug helper function #ifdef NDEBUG #define ERROR_IF(condition,message) if (false) {} #else #define ERROR_IF(condition,message) \ if ( condition ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror) #endif namespace ThePEG { template class LorentzVector; /** * A 4-component Lorentz vector. It can be created with any unit type * as template parameter. All basic mathematical operations are * supported, as well as a subset of the CLHEP LorentzVector * functionality. */ template class LorentzVector { private: /// Value squared using Value2 = decltype(sqr(std::declval())); public: /** @name Constructors. */ //@{ LorentzVector() : theX(), theY(), theZ(), theT() {} LorentzVector(Value x, Value y, Value z, Value t) : theX(x), theY(y), theZ(z), theT(t) {} LorentzVector(const ThreeVector & v, Value t) : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {} template LorentzVector(const LorentzVector & v) : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {} //@} /// Assignment operator template LorentzVector & operator=(const LorentzVector & b) { setX(b.x()); setY(b.y()); setZ(b.z()); setT(b.t()); return *this; } public: /// @name Component access methods. //@{ Value x() const { return theX; } Value y() const { return theY; } Value z() const { return theZ; } Value t() const { return theT; } Value e() const { return t(); } //@} /// @name Component set methods. //@{ void setX(Value x) { theX = x; } void setY(Value y) { theY = y; } void setZ(Value z) { theZ = z; } void setT(Value t) { theT = t; } void setE(Value e) { setT(e); } //@} public: /// Access to the 3-component part. ThreeVector vect() const { return ThreeVector(x(),y(),z()); } /// Cast to the 3-component part. operator ThreeVector() const { return vect(); } /// Set the 3-component part. void setVect(const ThreeVector & p) { theX = p.x(); theY = p.y(); theZ = p.z(); } public: /// The complex conjugate vector. LorentzVector conjugate() const { return LorentzVector(conj(x()),conj(y()),conj(z()),conj(t())); } /// Squared magnitude \f$x^\mu\,x_\mu=t^2 - \vec{x}^2\f$. Value2 m2() const { return (t()-z())*(t()+z()) - sqr(x()) - sqr(y()); } /// Squared magnitude with another vector Value2 m2(const LorentzVector & a) const { Value tt(a.t()+t()),zz(a.z()+z()); return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y()); } /// Magnitude (signed) \f$\pm\sqrt{|t^2 - \vec{x}^2|}\f$. Value m() const { Value2 tmp = m2(); return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp)); } /// Transverse mass squared \f$t^2-z^2\f$. Value2 mt2() const { return (t()-z())*(t()+z()); } /// Transverse mass (signed) \f$\pm\sqrt{|t^2 - z^2|}\f$. Value mt() const { Value2 tmp = mt2(); return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp)); } /// Squared transverse component of the spatial vector \f$x^2+y^2\f$. Value2 perp2() const { return sqr(x()) + sqr(y()); } /// Transverse component of the spatial vector \f$\pm\sqrt{x^2 + y^2}\f$. Value perp() const { return sqrt(perp2()); } /** * Squared transverse component of the spatial vector with respect to the * given axis. */ template Value2 perp2(const ThreeVector & p) const { return vect().perp2(p); } /** * Transverse component of the spatial vector with respect to the * given axis. */ template Value perp(const ThreeVector & p) const { return vect().perp(p); } /// Transverse energy squared. Value2 et2() const { Value2 pt2 = vect().perp2(); return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z()); } /// Transverse energy (signed). Value et() const { Value2 etet = et2(); return e() < Value() ? -sqrt(etet) : sqrt(etet); } /// Transverse energy squared with respect to the given axis. Value2 et2(const ThreeVector & v) const { Value2 pt2 = vect().perp2(v); Value pv = vect().dot(v.unit()); return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv); } /// Transverse energy with respect to the given axis (signed). Value et(const ThreeVector & v) const { Value2 etet = et2(v); return e() < Value() ? -sqrt(etet) : sqrt(etet); } /// @name Spherical coordinates for the spatial part. //@{ /// Radius squared. Value2 rho2() const { return sqr(x()) + sqr(y()) + sqr(z()); } /// Radius. Value rho() const { return sqrt(rho2()); } /// Set new radius. void setRho(Value newRho) { Value oldRho = rho(); if (oldRho == Value()) return; double factor = newRho / oldRho; setX(x()*factor); setY(y()*factor); setZ(z()*factor); } /// Polar angle. double theta() const { assert(!(x() == Value() && y() == Value() && z() == Value())); return atan2(perp(),z()); } /// Cosine of the polar angle. double cosTheta() const { Value ptot = rho(); assert( ptot > Value() ); return z() / ptot; } /// Azimuthal angle. double phi() const { return atan2(y(),x()) ; } //@} /// Pseudorapidity of spatial part. double eta() const { Value m = rho(); if ( m == Value() ) return 0.0; Value pt = max(Constants::epsilon*m, perp()); double rap = log((m + abs(z()))/pt); return z() > ZERO? rap: -rap; } /// Spatial angle with another vector. double angle(const LorentzVector & w) const { return vect().angle(w.vect()); } /// Rapidity \f$\frac{1}{2}\ln\frac{t+z}{t-z} \f$ double rapidity() const { if ( z() == ZERO ) return 0.0; ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector"); Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2() + m2())); double rap = log((t() + abs(z()))/pt); return z() > ZERO? rap: -rap; } /// Rapidity with respect to another vector double rapidity(const Axis & ref) const { double r = ref.mag2(); ERROR_IF(r == 0,"A zero vector used as reference to LorentzVector rapidity"); Value vdotu = vect().dot(ref)/sqrt(r); if ( vdotu == ZERO ) return 0.0; ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector"); Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2(ref) + m2())); double rap = log((t() + abs(z()))/pt); return z() > ZERO? rap: -rap; } /** * Boost from reference frame into this vector's rest * frame: \f$\frac{\vec{x}}{t}\f$. */ Boost boostVector() const { if (t() == Value()) { if (rho2() == Value2()) return Boost(); else ERROR_IF(true,"boostVector computed for LorentzVector with t=0 -- infinite result"); } // result will make analytic sense but is physically meaningless ERROR_IF(m2() <= Value2(),"boostVector computed for a non-timelike LorentzVector"); return vect() * (1./t()); } /** * Boost from reference frame into this vector's rest * frame: \f$-\frac{\vec{x}}{t}\f$. */ Boost findBoostToCM() const { return -boostVector(); } /// Returns the positive light-cone component \f$t + z\f$. Value plus() const { return t() + z(); } /// Returns the negative light-cone component \f$t - z\f$. Value minus() const { return t() - z(); } /// Are two vectors nearby, using Euclidean measure \f$t^2 + |\vec{x}|^2\f$? bool isNear(const LorentzVector & w, double epsilon) const { Value2 limit = abs(vect().dot(w.vect())); limit += 0.25 * sqr( t() + w.t() ); limit *= sqr(epsilon); Value2 delta = (vect() - w.vect()).mag2(); delta += sqr( t() - w.t() ); return (delta <= limit); } /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$. LorentzVector & transform(const SpinOneLorentzRotation & m) { return *this = m.operator*(*this); } /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$. LorentzVector & operator*=(const SpinOneLorentzRotation & m) { return transform(m); } /// Dot product with metric \f$(+,-,-,-)\f$ template - auto dot(const LorentzVector & a) const -> decltype(t() * a.t()) + auto dot(const LorentzVector & a) const -> decltype(this->t() * a.t()) { return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() ); } public: /** * Apply boost. * * @param bx Component x of the boost. * @param by Component y of the boost. * @param bz Component z of the boost. * @param gamma Optional gamma parameter for higher numerical * accuracy. The user has to ensure consistency. If not given, it * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$. * */ LorentzVector & boost(double bx, double by, double bz, double gamma=-1.) { const double b2 = bx*bx + by*by + bz*bz; if ( b2 == 0.0 ) return *this; if ( gamma < 0.0 ) { gamma = 1.0 / sqrt(1.0 - b2); } const Value bp = bx*x() + by*y() + bz*z(); const double gamma2 = (gamma - 1.0)/b2; setX(x() + gamma2*bp*bx + gamma*bx*t()); setY(y() + gamma2*bp*by + gamma*by*t()); setZ(z() + gamma2*bp*bz + gamma*bz*t()); setT(gamma*(t() + bp)); return *this; } /** * Apply boost. * * @param b Three-vector giving the boost. * * @param gamma Optional gamma parameter for higher numerical * accuracy. The user has to ensure consistency. If not given, it * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$. * */ LorentzVector & boost(Boost b, double gamma=-1.) { return boost(b.x(), b.y(), b.z(),gamma); } /** * Apply rotation around the x-axis. * * @param phi Angle in radians. */ LorentzVector & rotateX (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value ty = y() * cosphi - z() * sinphi; theZ = z() * cosphi + y() * sinphi; theY = ty; return *this; } /** * Apply rotation around the y-axis. * * @param phi Angle in radians. */ LorentzVector & rotateY (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value tz = z() * cosphi - x() * sinphi; theX = x() * cosphi + z() * sinphi; theZ = tz; return *this; } /** * Apply rotation around the z-axis. * * @param phi Angle in radians. */ LorentzVector & rotateZ (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value tx = x() * cosphi - y() * sinphi; theY = y() * cosphi + x() * sinphi; theX = tx; return *this; } /** * Rotate the reference frame to a new z-axis. */ LorentzVector & rotateUz (const Axis & axis) { Axis ax = axis.unit(); double u1 = ax.x(); double u2 = ax.y(); double u3 = ax.z(); double up = u1*u1 + u2*u2; if (up>0) { up = sqrt(up); Value px = x(), py = y(), pz = z(); setX( (u1*u3*px - u2*py)/up + u1*pz ); setY( (u2*u3*px + u1*py)/up + u2*pz ); setZ( -up*px + u3*pz ); } else if (u3 < 0.) { setX(-x()); setZ(-z()); } return *this; } /** * Apply a rotation. * @param angle Rotation angle in radians. * @param axis Rotation axis. */ template LorentzVector & rotate(double angle, const ThreeVector & axis) { if (angle == 0.0) return *this; const U ll = axis.mag(); assert( ll > U() ); const double sa = sin(angle), ca = cos(angle); const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll; const Value xx = x(), yy = y(), zz = z(); setX((ca+(1-ca)*dx*dx) * xx +((1-ca)*dx*dy-sa*dz) * yy +((1-ca)*dx*dz+sa*dy) * zz ); setY(((1-ca)*dy*dx+sa*dz) * xx +(ca+(1-ca)*dy*dy) * yy +((1-ca)*dy*dz-sa*dx) * zz ); setZ(((1-ca)*dz*dx-sa*dy) * xx +((1-ca)*dz*dy+sa*dx) * yy +(ca+(1-ca)*dz*dz) * zz ); return *this; } public: /// @name Mathematical assignment operators. //@{ LorentzVector & operator+=(const LorentzVector > & a) { theX += a.x(); theY += a.y(); theZ += a.z(); theT += a.t(); return *this; } template LorentzVector & operator+=(const LorentzVector & a) { theX += a.x(); theY += a.y(); theZ += a.z(); theT += a.t(); return *this; } LorentzVector & operator-=(const LorentzVector > & a) { theX -= Complex(a.x()); theY -= Complex(a.y()); theZ -= Complex(a.z()); theT -= Complex(a.t()); return *this; } template LorentzVector & operator-=(const LorentzVector & a) { theX -= a.x(); theY -= a.y(); theZ -= a.z(); theT -= a.t(); return *this; } LorentzVector & operator*=(double a) { theX *= a; theY *= a; theZ *= a; theT *= a; return *this; } LorentzVector & operator/=(double a) { theX /= a; theY /= a; theZ /= a; theT /= a; return *this; } //@} private: /// @name Vector components //@{ Value theX; Value theY; Value theZ; Value theT; //@} }; /// @name Basic mathematical operations //@{ template inline LorentzVector operator/(const LorentzVector & v, Value a) { return LorentzVector(v.x()/a, v.y()/a, v.z()/a, v.t()/a); } inline LorentzVector operator/(const LorentzVector & v, Complex a) { return LorentzVector(v.x()/a, v.y()/a, v.z()/a, v.t()/a); } template inline LorentzVector operator-(const LorentzVector & v) { return LorentzVector(-v.x(),-v.y(),-v.z(),-v.t()); } template inline LorentzVector operator+(LorentzVector a, const LorentzVector & b) { return a += b; } template inline LorentzVector operator-(LorentzVector a, const LorentzVector & b) { return a -= b; } template inline LorentzVector operator*(const LorentzVector & a, double b) { return LorentzVector(a.x()*b, a.y()*b, a.z()*b, a.t()*b); } template inline LorentzVector operator*(double b, LorentzVector a) { return a *= b; } template inline auto operator*(ValueB a, const LorentzVector & v) -> LorentzVector { return {a*v.x(), a*v.y(), a*v.z(), a*v.t()}; } template inline auto operator*(const LorentzVector & v, ValueB b) -> LorentzVector { return b*v; } template inline auto operator/(const LorentzVector & v, ValueB b) -> LorentzVector { return {v.x()/b, v.y()/b, v.z()/b, v.t()/b}; } //@} /// @name Scalar product with metric \f$(+,-,-,-)\f$ //@{ template inline auto operator*(const LorentzVector & a, const LorentzVector & b) -> decltype(a.dot(b)) { return a.dot(b); } //@} /// Equality template inline bool operator==(const LorentzVector & a, const LorentzVector & b) { return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t(); } /// Stream output. Format \f$(x,y,z;t)\f$. inline ostream & operator<< (ostream & os, const LorentzVector & v) { return os << "(" << v.x() << "," << v.y() << "," << v.z() << ";" << v.t() << ")"; } /** Return the positive light-cone component. Or negative if the * current Direction<0> is reversed. */ template inline Value dirPlus(const LorentzVector & p) { return Direction<0>::pos()? p.plus(): p.minus(); } /** Return the negative light-cone component. Or positive if the * current Direction<0> is reversed. */ template inline Value dirMinus(const LorentzVector & p) { return Direction<0>::neg()? p.plus(): p.minus(); } /** Return the component along the positive z-axis. Or the negative * z-axis if the current Direction<0> is reversed. */ template inline Value dirZ(const LorentzVector & p) { return Direction<0>::dir()*p.z(); } /** Return the polar angle wrt. the positive z-axis. Or the negative * z-axis if the current Direction<0> is reversed. */ template inline double dirTheta(const LorentzVector & p) { return Direction<0>::pos()? p.theta(): Constants::pi - p.theta(); } /** Return the cosine of the polar angle wrt. the positive z-axis. Or * the negative z-axis if the current Direction<0> is reversed. */ template inline double dirCosTheta(const LorentzVector & p) { return Direction<0>::pos()? p.cosTheta(): -p.cosTheta(); } /** Get the boost vector for the LorentzVector. If the current * Direction<0> is reversed, so is the z-component. */ template inline ThreeVector dirBoostVector(const LorentzVector & p) { ThreeVector b(p.boostVector()); if ( Direction<0>::neg() ) b.setZ(-b.z()); return b; } /** Create a LorentzVector giving its light-cone and transverse * components. */ template inline LorentzVector lightCone(Value plus, Value minus, Value x, Value y) { LorentzVector r(x, y, 0.5*(plus-minus), 0.5*(plus+minus)); return r; } /** Create a LorentzVector giving its light-cone components. */ template inline LorentzVector lightCone(Value plus, Value minus) { // g++-3.3 has a problem with using Value() directly // gcc-bug c++/3650, fixed in 3.4 static const Value zero = Value(); LorentzVector r(zero, zero, 0.5*(plus-minus), 0.5*(plus+minus)); return r; } } // delayed header inclusion to break inclusion loop: // LorentzVec -> Transverse -> Lorentz5Vec -> LorentzVec #include "Transverse.h" namespace ThePEG { /** Create a LorentzVector giving its light-cone and transverse * components. */ template inline LorentzVector lightCone(Value plus, Value minus, Transverse pt) { LorentzVector r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus)); return r; } /** Create a LorentzVector giving its light-cone and transverse * components. If the current Direction<0> is reversed, so is the * z-component. */ template inline LorentzVector lightConeDir(Value plus, Value minus, Value x = Value(), Value y = Value()) { LorentzVector r(x, y, Direction<0>::dir()*0.5*(plus - minus), 0.5*(plus + minus)); return r; } /** Create a LorentzVector giving its light-cone and transverse * components. If the current Direction<0> is reversed, so is the * z-component. */ template inline LorentzVector lightConeDir(Value plus, Value minus, Transverse pt) { LorentzVector r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus), 0.5*(plus + minus)); return r; } /** Output a LorentzVector with units to a stream. */ template void ounitstream(OStream & os, const LorentzVector & p, UnitT & u) { os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u) << ounit(p.e(), u); } /** Input a LorentzVector with units from a stream. */ template void iunitstream(IStream & is, LorentzVector & p, UnitT & u) { Value x, y, z, e; is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u); p = LorentzVector(x, y, z, e); } } #undef ERROR_IF #endif /* ThePEG_LorentzVector_H */ diff --git a/Vectors/ThreeVector.h b/Vectors/ThreeVector.h --- a/Vectors/ThreeVector.h +++ b/Vectors/ThreeVector.h @@ -1,421 +1,421 @@ // -*- C++ -*- // // ThreeVector.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_ThreeVector_H #define ThePEG_ThreeVector_H /** * @file ThreeVector.h contains the ThreeVector class. ThreeVector can be * created with any unit type as template parameter. All basic * mathematical operations are supported, as well as a subset of the * CLHEP Vector3 functionality. */ #include "ThreeVector.fh" #include "ThePEG/Config/ThePEG.h" #include "ThePEG/Utilities/UnitIO.h" #include #include namespace ThePEG { /** * A 3-component vector. It can be created with any unit type * as template parameter. All basic mathematical operations are * supported, as well as a subset of the CLHEP Vector3 * functionality. */ template class ThreeVector { private: /// Value squared using Value2 = decltype(sqr(std::declval())); public: /** @name Constructors. */ //@{ ThreeVector() : theX(), theY(), theZ() {} ThreeVector(Value x, Value y, Value z) : theX(x), theY(y), theZ(z) {} template ThreeVector(const ThreeVector & v) : theX(v.x()), theY(v.y()), theZ(v.z()) {} //@} public: /// @name Component access methods. //@{ Value x() const { return theX; } Value y() const { return theY; } Value z() const { return theZ; } //@} /// @name Component set methods. //@{ void setX(Value x) { theX = x; } void setY(Value y) { theY = y; } void setZ(Value z) { theZ = z; } //@} public: /// Squared magnitude \f$x^2+y^2+z^2\f$. Value2 mag2() const { return sqr(x()) + sqr(y()) + sqr(z()); } /// Magnitude \f$\sqrt{x^2+y^2+z^2}\f$. Value mag() const { return sqrt(mag2()); } /// Squared transverse component \f$x^2+y^2\f$. Value2 perp2() const { return sqr(x()) + sqr(y()); } /// Transverse component \f$\sqrt{x^2+y^2}\f$. Value perp() const { return sqrt(perp2()); } /// Dot product. template auto dot(const ThreeVector & a) const - -> decltype(x()*a.x()) + -> decltype(this->x()*a.x()) { return x()*a.x() + y()*a.y() + z()*a.z(); } /// Squared transverse component with respect to the given axis. template Value2 perp2(const ThreeVector & p) const { const auto pMag2 = p.mag2(); assert( pMag2 > ZERO ); auto ss = this->dot(p); Value2 ret = mag2() - sqr(ss)/pMag2; if ( ret <= ZERO ) ret = ZERO; return ret; } /// Transverse component with respect to the given axis. template Value perp(const ThreeVector & p) const { return sqrt(perp2(p)); } /// @name Spherical coordinates. //@{ /// Polar angle. double theta() const { assert(!(x() == ZERO && y() == ZERO && z() == ZERO)); return atan2(perp(),z()); } /// Azimuthal angle. double phi() const { return atan2(y(),x()); } /// Set the polar angle. void setTheta(double th) { double ma = mag(); double ph = phi(); setX(ma*sin(th)*cos(ph)); setY(ma*sin(th)*sin(ph)); setZ(ma*cos(th)); } /// Set the azimuthal angle. void setPhi(double ph) { double xy = perp(); setX(xy*cos(ph)); setY(xy*sin(ph)); } //@} /// Parallel vector with unit length. ThreeVector unit() const { Value2 mg2 = mag2(); assert(mg2 > ZERO); Value mg = sqrt(mg2); return {x()/mg, y()/mg, z()/mg}; } /// Orthogonal vector. ThreeVector orthogonal() const { Value xx = abs(x()); Value yy = abs(y()); Value zz = abs(z()); using TVec = ThreeVector; if (xx < yy) { return xx < zz ? TVec{ZERO,z(),-y()} : TVec{y(),-x(),ZERO}; } else { return yy < zz ? TVec{-z(),ZERO,x()} : TVec{y(),-x(),ZERO}; } } /// Azimuthal angle difference, brought into the range \f$(-\pi,\pi]\f$. template double deltaPhi (const ThreeVector & v2) const { double dphi = v2.phi() - phi(); if ( dphi > Constants::pi ) { dphi -= Constants::twopi; } else if ( dphi <= -Constants::pi ) { dphi += Constants::twopi; } return dphi; } /** * Apply a rotation. * @param angle Rotation angle in radians. * @param axis Rotation axis. */ template ThreeVector & rotate(double angle, const ThreeVector & axis) { if (angle == 0.0) return *this; const U ll = axis.mag(); assert( ll > ZERO ); const double sa = sin(angle), ca = cos(angle); const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll; const Value xx = x(), yy = y(), zz = z(); setX((ca+(1-ca)*dx*dx) * xx +((1-ca)*dx*dy-sa*dz) * yy +((1-ca)*dx*dz+sa*dy) * zz ); setY(((1-ca)*dy*dx+sa*dz) * xx +(ca+(1-ca)*dy*dy) * yy +((1-ca)*dy*dz-sa*dx) * zz ); setZ(((1-ca)*dz*dx-sa*dy) * xx +((1-ca)*dz*dy+sa*dx) * yy +(ca+(1-ca)*dz*dz) * zz ); return *this; } /** * Rotate the reference frame to a new z-axis. */ ThreeVector & rotateUz (const Axis & axis) { Axis ax = axis.unit(); double u1 = ax.x(); double u2 = ax.y(); double u3 = ax.z(); double up = u1*u1 + u2*u2; if (up>0) { up = sqrt(up); Value px = x(), py = y(), pz = z(); setX( (u1*u3*px - u2*py)/up + u1*pz ); setY( (u2*u3*px + u1*py)/up + u2*pz ); setZ( -up*px + u3*pz ); } else if (u3 < 0.) { setX(-x()); setZ(-z()); } return *this; } /** * Rotate from a reference frame to the z-axis. */ ThreeVector & rotateUzBack (const Axis & axis) { Axis ax = axis.unit(); double u1 = ax.x(); double u2 = ax.y(); double u3 = ax.z(); double up = u1*u1 + u2*u2; if (up>0) { up = sqrt(up); Value px = x(), py = y(), pz = z(); setX( ( u1*u3*px + u2*u3*py)/up - up*pz ); setY( (-u2*px + u1*py)/up ); setZ( u1*px + u2*py + u3*pz ); } else if (u3 < 0.) { setX(-x()); setZ(-z()); } return *this; } /// Vector cross-product template auto cross(const ThreeVector & a) const - -> ThreeVector + -> ThreeVectory()*a.z())> { return { y()*a.z()-z()*a.y(), -x()*a.z()+z()*a.x(), x()*a.y()-y()*a.x() }; } public: /// @name Comparison operators. //@{ bool operator==(const ThreeVector & a) const { return (theX == a.x() && theY == a.y() && theZ == a.z()); } bool operator!=(const ThreeVector & a) const { return !(*this == a); } bool almostEqual(const ThreeVector & a, double threshold = 1e-04) const { return ((std::abs(theX - a.x()) < threshold) && (std::abs(theY - a.y()) < threshold) && (std::abs(theZ - a.z()) < threshold)); } bool almostUnequal(const ThreeVector & a, double threshold = 1e-04) const { return ! this->almostEqual(a, threshold); } //@} public: /// @name Mathematical assignment operators. //@{ ThreeVector & operator+=(const ThreeVector & a) { theX += a.x(); theY += a.y(); theZ += a.z(); return *this; } ThreeVector & operator-=(const ThreeVector & a) { theX -= a.x(); theY -= a.y(); theZ -= a.z(); return *this; } ThreeVector & operator*=(double a) { theX *= a; theY *= a; theZ *= a; return *this; } ThreeVector & operator/=(double a) { theX /= a; theY /= a; theZ /= a; return *this; } //@} /// Cosine of the azimuthal angle between two vectors. template double cosTheta(const ThreeVector & q) const { auto ptot = mag()*q.mag(); assert( ptot > ZERO ); double arg = dot(q)/ptot; if (arg > 1.0) arg = 1.0; else if(arg < -1.0) arg = -1.0; return arg; } /// Angle between two vectors. template double angle(const ThreeVector & v) const { return acos(cosTheta(v)); } private: /// @name Vector components //@{ Value theX; Value theY; Value theZ; //@} }; /// Stream output. Format \f$(x,y,z)\f$. inline ostream & operator<< (ostream & os, const ThreeVector & v) { return os << '(' << v.x() << ',' << v.y() << ',' << v.z() << ')'; } /// @name Basic mathematical operations //@{ template inline ThreeVector operator+(ThreeVector a, const ThreeVector & b) { return a += b; } template inline ThreeVector operator-(ThreeVector a, const ThreeVector & b) { return a -= b; } template inline ThreeVector operator-(const ThreeVector & v) { return {-v.x(),-v.y(),-v.z()}; } template inline ThreeVector operator*(ThreeVector v, double a) { return v *= a; } template inline ThreeVector operator*(double a, ThreeVector v) { return v *= a; } template inline auto operator*(ValueB a, ThreeVector v) -> ThreeVector { return {a*v.x(), a*v.y(), a*v.z()}; } template inline auto operator*(ThreeVector v, ValueB a) -> ThreeVector { return {v.x()*a, v.y()*a, v.z()*a}; } //@} /// Vector dot product. template inline auto operator*(const ThreeVector & a, const ThreeVector & b) -> decltype(a.x()*b.x()) { return a.dot(b); } /// A parallel vector with unit length. template ThreeVector unitVector(const ThreeVector & v) { return v.unit(); } /** Output a ThreeVector with units to a stream. */ template void ounitstream(OStream & os, const ThreeVector & p, UT & u) { os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u); } /** Input a ThreeVector with units from a stream. */ template void iunitstream(IStream & is, ThreeVector & p, UT & u) { Value x, y, z; is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u); p = ThreeVector(x, y, z); } } #endif /* ThePEG_ThreeVector_H */