ThePEG/Vectors/LorentzVector.h

0001 // -*- C++ -*-
0002 //
0003 // LorentzVector.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_LorentzVector_H
0010 #define ThePEG_LorentzVector_H
0011
0012 /**
0013  * @file LorentzVector.h contains the LorentzVector class.  Lorentz
0014  * vectors can be created with any unit type as template parameter.
0015  * All basic mathematical operations are supported, as well as a
0016  * subset of the CLHEP LorentzVector functionality.
0017  */
0018
0019 #include "LorentzVector.fh"
0020 #include "ThePEG/Utilities/Direction.h"
0021 #include "ThePEG/Utilities/UnitIO.h"
0022 #include "LorentzRotation.h"
0023 #include "ThreeVector.h"
0024
0025 /// Debug helper function
0026 #ifdef NDEBUG
0027 #define ERROR_IF(condition,message) if (false) {}
0028 #else
0029 #define ERROR_IF(condition,message) \
0030   if ( condition ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror)
0031 #endif
0032
0033 namespace ThePEG {
0034
0035 template <typename Value> class LorentzVector;
0036
0037 /**
0038  * A 4-component Lorentz vector. It can be created with any unit type
0039  * as template parameter.  All basic mathematical operations are
0040  * supported, as well as a subset of the CLHEP LorentzVector
0041  * functionality.
0042  */
0043 template <typename Value> class LorentzVector
0044 {
0045 private:
0046   /// Value squared
0047   using Value2 = decltype(sqr(std::declval<Value>()));
0048
0049 public:
0050   /** @name Constructors. */
0051   //@{
0052   LorentzVector()
0053     : theX(), theY(), theZ(), theT() {}
0054
0055   LorentzVector(Value x, Value y, Value z, Value t)
0056     : theX(x), theY(y), theZ(z), theT(t) {}
0057
0058   LorentzVector(const ThreeVector<Value> & v, Value t)
0059     : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {}
0060
0061   template<typename U>
0062   LorentzVector(const LorentzVector<U> & v)
0063     : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {}
0064   //@}
0065
0066   /// Assignment operator
0067   template <typename ValueB>
0068   LorentzVector<Value> & operator=(const LorentzVector<ValueB> & b) {
0069     setX(b.x());
0070     setY(b.y());
0071     setZ(b.z());
0072     setT(b.t());
0073     return *this;
0074   }
0075
0076 public:
0077   /// @name Component access methods.
0078   //@{
0079   Value x() const { return theX; }
0080   Value y() const { return theY; }
0081   Value z() const { return theZ; }
0082   Value t() const { return theT; }
0083   Value e() const { return t();  }
0084   //@}
0085
0086   /// @name Component set methods.
0087   //@{
0088   void setX(Value x)  {  theX = x; }
0089   void setY(Value y)  {  theY = y; }
0090   void setZ(Value z)  {  theZ = z; }
0091   void setT(Value t)  {  theT = t; }
0092   void setE(Value e)  {  setT(e);  }
0093   //@}
0094
0095 public:
0096   /// Access to the 3-component part.
0097   ThreeVector<Value> vect() const {
0098     return ThreeVector<Value>(x(),y(),z());
0099   }
0100
0101   /// Cast to the 3-component part.
0102   operator ThreeVector<Value>() const { return vect(); }
0103
0104   /// Set the 3-component part.
0105   void setVect(const ThreeVector<Value> & p) {
0106     theX = p.x();
0107     theY = p.y();
0108     theZ = p.z();
0109   }
0110
0111 public:
0112   /// The complex conjugate vector.
0113   LorentzVector<Value> conjugate() const
0114   {
0115     return LorentzVector<Value>(conj(x()),conj(y()),conj(z()),conj(t()));
0116   }
0117
0118   /// Squared magnitude \f$x^\mu\,x_\mu=t^2 - \vec{x}^2\f$.
0119   Value2 m2() const
0120   {
0121     return (t()-z())*(t()+z()) - sqr(x()) - sqr(y());
0122   }
0123
0124   /// Squared magnitude with another vector
0125   Value2 m2(const LorentzVector<Value> & a) const {
0126     Value tt(a.t()+t()),zz(a.z()+z());
0127     return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y());
0128   }
0129
0130   /// Magnitude (signed) \f$\pm\sqrt{|t^2 - \vec{x}^2|}\f$.
0131   Value  m() const
0132   {
0133     Value2 tmp = m2();
0134     return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
0135   }
0136
0137   /// Transverse mass squared \f$t^2-z^2\f$.
0138   Value2 mt2()  const { return (t()-z())*(t()+z()); }
0139
0140   /// Transverse mass (signed) \f$\pm\sqrt{|t^2 - z^2|}\f$.
0141   Value  mt()  const
0142   {
0143     Value2 tmp = mt2();
0144     return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
0145   }
0146
0147   /// Squared transverse component of the spatial vector \f$x^2+y^2\f$.
0148   Value2 perp2() const { return sqr(x()) + sqr(y()); }
0149
0150   /// Transverse component of the spatial vector \f$\pm\sqrt{x^2 + y^2}\f$.
0151   Value  perp()  const { return sqrt(perp2()); }
0152
0153   /**
0154    * Squared transverse component of the spatial vector with respect to the
0155    * given axis.
0156    */
0157   template <typename U>
0158   Value2 perp2(const ThreeVector<U> & p) const
0159   {
0160     return vect().perp2(p);
0161   }
0162
0163   /**
0164    * Transverse component of the spatial vector with respect to the
0165    * given axis.
0166    */
0167   template <typename U>
0168   Value perp(const ThreeVector<U> & p) const
0169   {
0170     return vect().perp(p);
0171   }
0172
0173   /// Transverse energy squared.
0174   Value2 et2() const
0175   {
0176     Value2 pt2 = vect().perp2();
0177     return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z());
0178   }
0179
0180   /// Transverse energy (signed).
0181   Value et() const
0182   {
0183     Value2 etet = et2();
0184     return e() < Value() ? -sqrt(etet) : sqrt(etet);
0185   }
0186
0187   /// Transverse energy squared with respect to the given axis.
0188   Value2 et2(const ThreeVector<double> & v) const
0189   {
0190     Value2 pt2 = vect().perp2(v);
0191     Value pv = vect().dot(v.unit());
0192     return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv);
0193   }
0194
0195   /// Transverse energy with respect to the given axis (signed).
0196   Value et(const ThreeVector<double> & v) const
0197   {
0198     Value2 etet = et2(v);
0199     return e() < Value() ? -sqrt(etet) : sqrt(etet);
0200   }
0201
0202   /// @name Spherical coordinates for the spatial part.
0203   //@{
0204   /// Radius squared.
0205   Value2 rho2()  const { return sqr(x()) + sqr(y()) + sqr(z()); }
0206
0207   /// Radius.
0208   Value  rho()   const { return sqrt(rho2()); }
0209
0210   /// Set new radius.
0211   void setRho(Value newRho)
0212   {
0213     Value oldRho = rho();
0214     if (oldRho == Value())
0215       return;
0216     double factor = newRho / oldRho;
0217     setX(x()*factor);
0218     setY(y()*factor);
0219     setZ(z()*factor);
0220   }
0221
0222   /// Polar angle.
0223   double theta() const
0224   {
0225     assert(!(x() == Value() && y() == Value() && z() == Value()));
0226     return atan2(perp(),z());
0227   }
0228
0229   /// Cosine of the polar angle.
0230   double cosTheta() const
0231   {
0232     Value ptot = rho();
0233     assert( ptot > Value() );
0234     return z() / ptot;
0235   }
0236
0237   /// Azimuthal angle.
0238   double phi()   const {
0239     return atan2(y(),x()) ;
0240   }
0241   //@}
0242
0243   /// Pseudorapidity of spatial part.
0244   double eta() const {
0245     Value m = rho();
0246     if ( m ==  Value() ) return  0.0;
0247     Value pt = max(Constants::epsilon*m, perp());
0248     double rap = log((m + abs(z()))/pt);
0249     return z() > ZERO? rap: -rap;
0250   }
0251
0252   /// Spatial angle with another vector.
0253   double angle(const LorentzVector<Value> & w) const
0254   {
0255     return vect().angle(w.vect());
0256   }
0257
0258   /// Rapidity \f$\frac{1}{2}\ln\frac{t+z}{t-z} \f$
0259   double rapidity() const {
0260     if ( z() == ZERO ) return 0.0;
0261     ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector");
0262     Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2() + m2()));
0263     double rap = log((t() + abs(z()))/pt);
0264     return z() > ZERO? rap: -rap;
0265   }
0266
0267   /// Rapidity with respect to another vector
0268   double rapidity(const Axis & ref) const {
0269     double r = ref.mag2();
0270     ERROR_IF(r == 0,"A zero vector used as reference to LorentzVector rapidity");
0271     Value vdotu = vect().dot(ref)/sqrt(r);
0272     if ( vdotu == ZERO ) return 0.0;
0273     ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector");
0274     Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2(ref) + m2()));
0275     double rap = log((t() + abs(z()))/pt);
0276     return z() > ZERO? rap: -rap;
0277   }
0278
0279   /**
0280    * Boost from reference frame into this vector's rest
0281    * frame: \f$\frac{\vec{x}}{t}\f$.
0282    */
0283   Boost boostVector() const {
0284     if (t() == Value()) {
0285       if (rho2() == Value2())
0286     return Boost();
0287       else
0288     ERROR_IF(true,"boostVector computed for LorentzVector with t=0 -- infinite result");
0289     }
0290     // result will make analytic sense but is physically meaningless
0291     ERROR_IF(m2() <= Value2(),"boostVector computed for a non-timelike LorentzVector");
0292     return vect() * (1./t());
0293   }
0294
0295   /**
0296    * Boost from reference frame into this vector's rest
0297    * frame: \f$-\frac{\vec{x}}{t}\f$.
0298    */
0299   Boost findBoostToCM() const
0300   {
0301     return -boostVector();
0302   }
0303
0304   /// Returns the positive light-cone component \f$t + z\f$.
0305   Value plus()  const { return t() + z(); }
0306   /// Returns the negative light-cone component \f$t - z\f$.
0307   Value minus() const { return t() - z(); }
0308
0309   /// Are two vectors nearby, using Euclidean measure \f$t^2 + |\vec{x}|^2\f$?
0310   bool isNear(const LorentzVector<Value> & w, double epsilon) const
0311   {
0312     Value2 limit = abs(vect().dot(w.vect()));
0313     limit += 0.25 * sqr( t() + w.t() );
0314     limit *= sqr(epsilon);
0315     Value2 delta = (vect() - w.vect()).mag2();
0316     delta +=  sqr( t() - w.t() );
0317     return (delta <= limit);
0318   }
0319
0320   /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$.
0321   LorentzVector<Value> & transform(const SpinOneLorentzRotation & m)
0322   {
0323     return *this = m.operator*(*this);
0324   }
0325
0326   /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$.
0327   LorentzVector<Value> & operator*=(const SpinOneLorentzRotation & m)
0328   {
0329     return transform(m);
0330   }
0331
0332   /// Dot product with metric \f$(+,-,-,-)\f$
0333   template <typename U>
0334   auto dot(const LorentzVector<U> & a) const -> decltype(this->t() * a.t())
0335   {
0336     return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() );
0337   }
0338
0339
0340 public:
0341
0342   /**
0343    * Apply boost.
0344    *
0345    * @param bx Component x of the boost.
0346    * @param by Component y of the boost.
0347    * @param bz Component z of the boost.
0348    * @param gamma Optional gamma parameter for higher numerical
0349    * accuracy. The user has to ensure consistency. If not given, it
0350    * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$.
0351    *
0352    */
0353   LorentzVector<Value> &
0354   boost(double bx, double by, double bz, double gamma=-1.)
0355   {
0356     const double b2 = bx*bx + by*by + bz*bz;
0357     if ( b2 == 0.0 ) return *this;
0358     if ( gamma < 0.0 ) {
0359         gamma = 1.0 / sqrt(1.0 - b2);
0360     }
0361     const Value bp = bx*x() + by*y() + bz*z();
0362     const double gamma2 = (gamma - 1.0)/b2;
0363
0364     setX(x() + gamma2*bp*bx + gamma*bx*t());
0365     setY(y() + gamma2*bp*by + gamma*by*t());
0366     setZ(z() + gamma2*bp*bz + gamma*bz*t());
0367     setT(gamma*(t() + bp));
0368     return *this;
0369   }
0370
0371   /**
0372    * Apply boost.
0373    *
0374    * @param b Three-vector giving the boost.
0375    *
0376    * @param gamma Optional gamma parameter for higher numerical
0377    * accuracy. The user has to ensure consistency. If not given, it
0378    * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$.
0379    *
0380    */
0381   LorentzVector<Value> &  boost(Boost b, double gamma=-1.) {
0382     return boost(b.x(), b.y(), b.z(),gamma);
0383   }
0384
0385   /**
0386    * Apply rotation around the x-axis.
0387    *
0388    * @param phi Angle in radians.
0389    */
0390   LorentzVector<Value> & rotateX (double phi) {
0391     double sinphi = sin(phi);
0392     double cosphi = cos(phi);
0393     Value  ty = y() * cosphi - z() * sinphi;
0394     theZ      = z() * cosphi + y() * sinphi;
0395     theY = ty;
0396     return *this;
0397   }
0398
0399   /**
0400    * Apply rotation around the y-axis.
0401    *
0402    * @param phi Angle in radians.
0403    */
0404   LorentzVector<Value> & rotateY (double phi) {
0405     double sinphi = sin(phi);
0406     double cosphi = cos(phi);
0407     Value  tz = z() * cosphi - x() * sinphi;
0408     theX      = x() * cosphi + z() * sinphi;
0409     theZ = tz;
0410     return *this;
0411   }
0412
0413   /**
0414    * Apply rotation around the z-axis.
0415    *
0416    * @param phi Angle in radians.
0417    */
0418   LorentzVector<Value> & rotateZ (double phi) {
0419     double sinphi = sin(phi);
0420     double cosphi = cos(phi);
0421     Value  tx = x() * cosphi - y() * sinphi;
0422     theY      = y() * cosphi + x() * sinphi;
0423     theX = tx;
0424     return *this;
0425   }
0426
0427   /**
0428    * Rotate the reference frame to a new z-axis.
0429    */
0430   LorentzVector<Value> & rotateUz (const Axis & axis) {
0431     Axis ax = axis.unit();
0432     double u1 = ax.x();
0433     double u2 = ax.y();
0434     double u3 = ax.z();
0435     double up = u1*u1 + u2*u2;
0436     if (up>0) {
0437       up = sqrt(up);
0438       Value px = x(),  py = y(),  pz = z();
0439       setX( (u1*u3*px - u2*py)/up + u1*pz );
0440       setY( (u2*u3*px + u1*py)/up + u2*pz );
0441       setZ(    -up*px +             u3*pz );
0442     }
0443     else if (u3 < 0.) {
0444       setX(-x());
0445       setZ(-z());
0446     }
0447     return *this;
0448   }
0449
0450   /**
0451    * Apply a rotation.
0452    * @param angle Rotation angle in radians.
0453    * @param axis Rotation axis.
0454    */
0455   template <typename U>
0456   LorentzVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
0457     if (angle == 0.0)
0458       return *this;
0459     const U ll = axis.mag();
0460     assert( ll > U() );
0461
0462     const double sa = sin(angle), ca = cos(angle);
0463     const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
0464     const Value  xx  = x(), yy = y(), zz = z();
0465
0466     setX((ca+(1-ca)*dx*dx)     * xx
0467      +((1-ca)*dx*dy-sa*dz) * yy
0468      +((1-ca)*dx*dz+sa*dy) * zz
0469      );
0470     setY(((1-ca)*dy*dx+sa*dz)  * xx
0471      +(ca+(1-ca)*dy*dy)    * yy
0472      +((1-ca)*dy*dz-sa*dx) * zz
0473      );
0474     setZ(((1-ca)*dz*dx-sa*dy)  * xx
0475      +((1-ca)*dz*dy+sa*dx) * yy
0476      +(ca+(1-ca)*dz*dz)    * zz
0477      );
0478     return *this;
0479   }
0480
0481
0482
0483
0484 public:
0485   /// @name Mathematical assignment operators.
0486   //@{
0487   LorentzVector<Complex> & operator+=(const LorentzVector<complex<QtyDouble> > & a) {
0488     theX += a.x();
0489     theY += a.y();
0490     theZ += a.z();
0491     theT += a.t();
0492     return *this;
0493   }
0494
0495   template <typename ValueB>
0496   LorentzVector<Value> & operator+=(const LorentzVector<ValueB> & a) {
0497     theX += a.x();
0498     theY += a.y();
0499     theZ += a.z();
0500     theT += a.t();
0501     return *this;
0502   }
0503
0504   LorentzVector<Complex> & operator-=(const LorentzVector<complex<QtyDouble> > & a) {
0505     theX -= Complex(a.x());
0506     theY -= Complex(a.y());
0507     theZ -= Complex(a.z());
0508     theT -= Complex(a.t());
0509     return *this;
0510   }
0511
0512   template <typename ValueB>
0513   LorentzVector<Value> & operator-=(const LorentzVector<ValueB> & a) {
0514     theX -= a.x();
0515     theY -= a.y();
0516     theZ -= a.z();
0517     theT -= a.t();
0518     return *this;
0519   }
0520
0521   LorentzVector<Value> & operator*=(double a) {
0522     theX *= a;
0523     theY *= a;
0524     theZ *= a;
0525     theT *= a;
0526     return *this;
0527   }
0528
0529   LorentzVector<Value> & operator/=(double a) {
0530     theX /= a;
0531     theY /= a;
0532     theZ /= a;
0533     theT /= a;
0534     return *this;
0535   }
0536   //@}
0537
0538 private:
0539   /// @name Vector components
0540   //@{
0541   Value theX;
0542   Value theY;
0543   Value theZ;
0544   Value theT;
0545   //@}
0546 };
0547
0548 /// @name Basic mathematical operations
0549 //@{
0550 template <typename Value>
0551 inline LorentzVector<double>
0552 operator/(const LorentzVector<Value> & v, Value a) {
0553   return LorentzVector<double>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
0554 }
0555
0556 inline LorentzVector<Complex>
0557 operator/(const LorentzVector<Complex> & v, Complex a) {
0558   return LorentzVector<Complex>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
0559 }
0560
0561 template <typename Value>
0562 inline LorentzVector<Value> operator-(const LorentzVector<Value> & v) {
0563   return LorentzVector<Value>(-v.x(),-v.y(),-v.z(),-v.t());
0564 }
0565
0566 template <typename ValueA, typename ValueB>
0567 inline LorentzVector<ValueA>
0568 operator+(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
0569   return a += b;
0570 }
0571
0572 template <typename ValueA, typename ValueB>
0573 inline LorentzVector<ValueA>
0574 operator-(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
0575   return a -= b;
0576 }
0577
0578 template <typename Value>
0579 inline LorentzVector<Value>
0580 operator*(const LorentzVector<Value> & a, double b) {
0581   return LorentzVector<Value>(a.x()*b, a.y()*b, a.z()*b, a.t()*b);
0582 }
0583
0584 template <typename Value>
0585 inline LorentzVector<Value>
0586 operator*(double b, LorentzVector<Value> a) {
0587   return a *= b;
0588 }
0589
0590 template <typename ValueA, typename ValueB>
0591 inline auto operator*(ValueB a, const LorentzVector<ValueA> & v)
0592 -> LorentzVector<decltype(a*v.x())>
0593 {
0594   return {a*v.x(), a*v.y(), a*v.z(), a*v.t()};
0595 }
0596
0597 template <typename ValueA, typename ValueB>
0598 inline auto operator*(const LorentzVector<ValueA> & v, ValueB b)
0599 -> LorentzVector<decltype(b*v.x())>
0600 {
0601   return b*v;
0602 }
0603
0604 template <typename ValueA, typename ValueB>
0605 inline auto operator/(const LorentzVector<ValueA> & v, ValueB b)
0606 -> LorentzVector<decltype(v.x()/b)>
0607 {
0608   return {v.x()/b, v.y()/b, v.z()/b, v.t()/b};
0609 }
0610 //@}
0611
0612 /// @name Scalar product with metric \f$(+,-,-,-)\f$
0613 //@{
0614 template <typename ValueA, typename ValueB>
0615 inline auto
0616 operator*(const LorentzVector<ValueA> & a, const LorentzVector<ValueB> & b)
0617 -> decltype(a.dot(b))
0618 {
0619   return a.dot(b);
0620 }
0621
0622 //@}
0623
0624 /// Equality
0625 template <typename Value>
0626 inline bool
0627 operator==(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
0628   return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t();
0629 }
0630
0631 /// Stream output. Format \f$(x,y,z;t)\f$.
0632 inline ostream & operator<< (ostream & os, const LorentzVector<double> & v) {
0633   return os << "(" << v.x() << "," << v.y() << "," << v.z()
0634         << ";" << v.t() << ")";
0635 }
0636
0637 /** Return the positive light-cone component. Or negative if the
0638  *  current Direction<0> is reversed. */
0639 template <typename Value>
0640 inline Value dirPlus(const LorentzVector<Value> & p) {
0641   return Direction<0>::pos()? p.plus(): p.minus();
0642 }
0643
0644 /** Return the negative light-cone component. Or positive if the
0645  *  current Direction<0> is reversed. */
0646 template <typename Value>
0647 inline Value dirMinus(const LorentzVector<Value> & p) {
0648   return Direction<0>::neg()? p.plus(): p.minus();
0649 }
0650
0651 /** Return the component along the positive z-axis. Or the negative
0652  *  z-axis if the current Direction<0> is reversed. */
0653 template <typename Value>
0654 inline Value dirZ(const LorentzVector<Value> & p) {
0655   return Direction<0>::dir()*p.z();
0656 }
0657
0658 /** Return the polar angle wrt. the positive z-axis. Or the negative
0659  *  z-axis if the current Direction<0> is reversed. */
0660 template <typename Value>
0661 inline double dirTheta(const LorentzVector<Value> & p) {
0662   return Direction<0>::pos()? p.theta(): Constants::pi - p.theta();
0663 }
0664
0665 /** Return the cosine of the polar angle wrt. the positive z-axis. Or
0666  *  the negative z-axis if the current Direction<0> is reversed. */
0667 template <typename Value>
0668 inline double dirCosTheta(const LorentzVector<Value> & p) {
0669   return Direction<0>::pos()? p.cosTheta():  -p.cosTheta();
0670 }
0671
0672 /** Get the boost vector for the LorentzVector. If the current
0673  *  Direction<0> is reversed, so is the z-component. */
0674 template <typename Value>
0675 inline ThreeVector<Value> dirBoostVector(const LorentzVector<Value> & p) {
0676   ThreeVector<Value> b(p.boostVector());
0677   if ( Direction<0>::neg() ) b.setZ(-b.z());
0678   return b;
0679 }
0680
0681 /** Create a LorentzVector giving its light-cone and transverse
0682  *  components. */
0683 template <typename Value>
0684 inline LorentzVector<Value>
0685 lightCone(Value plus, Value minus, Value x, Value y) {
0686   LorentzVector<Value> r(x, y, 0.5*(plus-minus), 0.5*(plus+minus));
0687   return r;
0688 }
0689
0690 /** Create a LorentzVector giving its light-cone components. */
0691 template <typename Value>
0692 inline LorentzVector<Value>
0693 lightCone(Value plus, Value minus) {
0694   // g++-3.3 has a problem with using Value() directly
0695   // gcc-bug c++/3650, fixed in 3.4
0696   static const Value zero = Value();
0697   LorentzVector<Value> r(zero, zero,
0698              0.5*(plus-minus), 0.5*(plus+minus));
0699   return r;
0700 }
0701
0702 }
0703
0704
0705 // delayed header inclusion to break inclusion loop:
0706 // LorentzVec -> Transverse -> Lorentz5Vec -> LorentzVec
0707 #include "Transverse.h"
0708
0709
0710
0711 namespace ThePEG {
0712
0713 /** Create a LorentzVector giving its light-cone and transverse
0714  *  components. */
0715 template <typename Value>
0716 inline LorentzVector<Value>
0717 lightCone(Value plus, Value minus, Transverse<Value> pt) {
0718   LorentzVector<Value> r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus));
0719   return r;
0720 }
0721
0722 /** Create a LorentzVector giving its light-cone and transverse
0723  *  components. If the current Direction<0> is reversed, so is the
0724  *  z-component. */
0725 template <typename Value>
0726 inline LorentzVector<Value>
0727 lightConeDir(Value plus, Value minus,
0728          Value x = Value(), Value y = Value()) {
0729   LorentzVector<Value> r(x, y, Direction<0>::dir()*0.5*(plus - minus),
0730           0.5*(plus + minus));
0731   return r;
0732 }
0733
0734 /** Create a LorentzVector giving its light-cone and transverse
0735  *  components. If the current Direction<0> is reversed, so is the
0736  *  z-component. */
0737 template <typename Value>
0738 inline LorentzVector<Value>
0739 lightConeDir(Value plus, Value minus, Transverse<Value> pt) {
0740   LorentzVector<Value> r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus),
0741           0.5*(plus + minus));
0742   return r;
0743
0744 }
0745
0746 /** Output a LorentzVector with units to a stream. */
0747 template <typename OStream, typename UnitT, typename Value>
0748 void ounitstream(OStream & os, const LorentzVector<Value> & p, UnitT & u) {
0749   os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u)
0750      << ounit(p.e(), u);
0751 }
0752
0753 /** Input a LorentzVector with units from a stream. */
0754 template <typename IStream, typename UnitT, typename Value>
0755 void iunitstream(IStream & is, LorentzVector<Value> & p, UnitT & u) {
0756   Value x, y, z, e;
0757   is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u);
0758   p = LorentzVector<Value>(x, y, z, e);
0759 }
0760
0761
0762 }
0763
0764 #undef ERROR_IF
0765 #endif /* ThePEG_LorentzVector_H */