9 #ifndef ThePEG_LorentzVector_H
10 #define ThePEG_LorentzVector_H
19 #include "LorentzVector.fh"
20 #include "ThePEG/Utilities/Direction.h"
21 #include "ThePEG/Utilities/UnitIO.h"
22 #include "LorentzRotation.h"
27 #define ERROR_IF(condition,message) if (false) {}
29 #define ERROR_IF(condition,message) \
30 if ( (condition) ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror)
53 : theX(), theY(), theZ(), theT() {}
56 : theX(x), theY(y), theZ(z), theT(t) {}
58 LorentzVector(
const ThreeVector<Value> & v, Value t)
59 : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {}
62 LorentzVector(
const LorentzVector<U> & v)
63 : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {}
67 template <
typename ValueB>
79 Value x()
const {
return theX; }
80 Value y()
const {
return theY; }
81 Value z()
const {
return theZ; }
82 Value t()
const {
return theT; }
83 Value e()
const {
return t(); }
88 void setX(Value x) { theX = x; }
89 void setY(Value y) { theY = y; }
90 void setZ(Value z) { theZ = z; }
91 void setT(Value t) { theT = t; }
92 void setE(Value e) { setT(e); }
121 return (t()-z())*(t()+z()) - sqr(x()) - sqr(y());
126 Value tt(a.t()+t()),zz(a.z()+z());
127 return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y());
134 return tmp <
Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
138 Value2
mt2()
const {
return (t()-z())*(t()+z()); }
144 return tmp <
Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
148 Value2
perp2()
const {
return sqr(x()) + sqr(y()); }
157 template <
typename U>
160 return vect().perp2(p);
167 template <
typename U>
170 return vect().perp(p);
176 Value2 pt2 =
vect().perp2();
177 return pt2 ==
Value2() ?
Value2() : e()*e() * pt2/(pt2+z()*z());
184 return e() < Value() ? -sqrt(etet) : sqrt(etet);
190 Value2 pt2 =
vect().perp2(v);
192 return pt2 ==
Value2() ?
Value2() : e()*e() * pt2/(pt2+pv*pv);
198 Value2 etet =
et2(v);
199 return e() < Value() ? -sqrt(etet) : sqrt(etet);
204 Value2
rho2()
const {
return sqr(x()) + sqr(y()) + sqr(z()); }
213 Value oldRho =
rho();
214 if (oldRho == Value())
216 double factor = newRho / oldRho;
225 assert(!(x() == Value() && y() == Value() && z() == Value()));
226 return atan2(
perp(),z());
233 assert( ptot > Value() );
239 return atan2(y(),x()) ;
246 if ( m == Value() )
return 0.0;
248 "Pseudorapidity for 3-vector along z-axis undefined.");
249 return 0.5 * log( (m+z()) / (m-z()) );
260 ERROR_IF(abs(t()) == abs(z()),
"rapidity for 4-vector with |E| = |Pz| -- infinite result");
261 ERROR_IF(abs(t()) < abs(z()) ,
"rapidity for spacelike 4-vector with |E| < |Pz| -- undefined");
262 double q = (t() + z()) / (t() - z());
268 double r = ref.
mag2();
269 ERROR_IF(r == 0,
"A zero vector used as reference to LorentzVector rapidity");
270 Value vdotu =
vect().dot(ref)/sqrt(r);
271 ERROR_IF(abs(t()) == abs(vdotu),
"rapidity for 4-vector with |E| = |Pu| -- infinite result");
272 ERROR_IF(abs(t()) < abs(vdotu),
"rapidity for spacelike 4-vector with |E|<|P*ref| undefined");
273 double q = (t() + vdotu) / (t() - vdotu);
282 if (t() == Value()) {
286 ERROR_IF(
true,
"boostVector computed for LorentzVector with t=0 -- infinite result");
289 ERROR_IF(
m2() <=
Value2(),
"boostVector computed for a non-timelike LorentzVector");
290 return vect() * (1./t());
303 Value
plus()
const {
return t() + z(); }
305 Value
minus()
const {
return t() - z(); }
311 limit += 0.25 * sqr( t() + w.t() );
312 limit *= sqr(epsilon);
313 Value2 delta = (
vect() - w.
vect()).mag2();
314 delta += sqr( t() - w.t() );
315 return (delta <= limit);
321 return *
this = m.operator*(*this);
331 template <
typename U>
335 return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() );
353 boost(
double bx,
double by,
double bz,
double gamma=-1.)
355 double b2 = bx*bx + by*by + bz*bz;
357 gamma = 1.0 / sqrt(1.0 - b2);
358 Value bp = bx*x() + by*y() + bz*z();
359 double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;
361 setX(x() + gamma2*bp*bx + gamma*bx*t());
362 setY(y() + gamma2*bp*by + gamma*by*t());
363 setZ(z() + gamma2*bp*bz + gamma*bz*t());
364 setT(gamma*(t() + bp));
379 return boost(b.x(), b.y(), b.z(),gamma);
388 double sinphi = sin(phi);
389 double cosphi = cos(phi);
390 Value ty = y() * cosphi - z() * sinphi;
391 theZ = z() * cosphi + y() * sinphi;
402 double sinphi = sin(phi);
403 double cosphi = cos(phi);
404 Value tz = z() * cosphi - x() * sinphi;
405 theX = x() * cosphi + z() * sinphi;
416 double sinphi = sin(phi);
417 double cosphi = cos(phi);
418 Value tx = x() * cosphi - y() * sinphi;
419 theY = y() * cosphi + x() * sinphi;
432 double up = u1*u1 + u2*u2;
435 Value px = x(), py = y(), pz = z();
436 setX( (u1*u3*px - u2*py)/up + u1*pz );
437 setY( (u2*u3*px + u1*py)/up + u2*pz );
438 setZ( -up*px + u3*pz );
452 template <
typename U>
456 const U ll = axis.
mag();
459 const double sa = sin(angle), ca = cos(angle);
460 const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
461 const Value xx = x(), yy = y(), zz = z();
463 setX((ca+(1-ca)*dx*dx) * xx
464 +((1-ca)*dx*dy-sa*dz) * yy
465 +((1-ca)*dx*dz+sa*dy) * zz
467 setY(((1-ca)*dy*dx+sa*dz) * xx
468 +(ca+(1-ca)*dy*dy) * yy
469 +((1-ca)*dy*dz-sa*dx) * zz
471 setZ(((1-ca)*dz*dx-sa*dy) * xx
472 +((1-ca)*dz*dy+sa*dx) * yy
473 +(ca+(1-ca)*dz*dz) * zz
484 template <
typename ValueB>
493 template <
typename ValueB>
494 LorentzVector<Value> & operator-=(
const LorentzVector<ValueB> & a) {
510 LorentzVector<Value> & operator/=(
double a) {
531 template <
typename Value>
532 inline LorentzVector<double>
533 operator/(
const LorentzVector<Value> & v, Value a) {
534 return LorentzVector<double>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
537 inline LorentzVector<Complex>
538 operator/(
const LorentzVector<Complex> & v,
Complex a) {
539 return LorentzVector<Complex>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
542 template <
typename Value>
543 inline LorentzVector<Value> operator-(
const LorentzVector<Value> & v) {
544 return LorentzVector<Value>(-v.x(),-v.y(),-v.z(),-v.t());
547 template <
typename ValueA,
typename ValueB>
548 inline LorentzVector<ValueA>
549 operator+(LorentzVector<ValueA> a,
const LorentzVector<ValueB> & b) {
553 template <
typename ValueA,
typename ValueB>
554 inline LorentzVector<ValueA>
555 operator-(LorentzVector<ValueA> a,
const LorentzVector<ValueB> & b) {
559 template <
typename Value>
560 inline LorentzVector<Value>
561 operator*(
const LorentzVector<Value> & a,
double b) {
562 return LorentzVector<Value>(a.x()*b, a.y()*b, a.z()*b, a.t()*b);
565 template <
typename Value>
566 inline LorentzVector<Value>
567 operator*(
double b, LorentzVector<Value> a) {
571 template <
typename ValueA,
typename ValueB>
573 LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
574 operator*(ValueB a,
const LorentzVector<ValueA> & v) {
575 typedef typename BinaryOpTraits<ValueB,ValueA>::MulT ResultT;
576 return LorentzVector<ResultT>(a*v.x(), a*v.y(), a*v.z(), a*v.t());
579 template <
typename ValueA,
typename ValueB>
581 LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
582 operator*(
const LorentzVector<ValueA> & v, ValueB b) {
586 template <
typename ValueA,
typename ValueB>
588 LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::DivT>
589 operator/(
const LorentzVector<ValueA> & v, ValueB b) {
590 typedef typename BinaryOpTraits<ValueA,ValueB>::DivT ResultT;
591 return LorentzVector<ResultT>(v.x()/b, v.y()/b, v.z()/b, v.t()/b);
597 template <
typename ValueA,
typename ValueB>
598 inline typename BinaryOpTraits<ValueA,ValueB>::MulT
599 operator*(
const LorentzVector<ValueA> & a,
const LorentzVector<ValueB> & b) {
603 template <
typename Value>
604 inline typename BinaryOpTraits<Value,Value>::MulT
605 operator*(
const LorentzVector<Value> & a,
const LorentzVector<Value> & b) {
611 template <
typename Value>
614 return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t();
618 inline ostream & operator<< (ostream & os, const LorentzVector<double> & v) {
619 return os <<
"(" << v.x() <<
"," << v.y() <<
"," << v.z()
620 <<
";" << v.t() <<
")";
625 template <
typename Value>
632 template <
typename Value>
639 template <
typename Value>
646 template <
typename Value>
653 template <
typename Value>
660 template <
typename Value>
669 template <
typename Value>
670 inline LorentzVector<Value>
677 template <
typename Value>
678 inline LorentzVector<Value>
682 static const Value zero = Value();
684 0.5*(plus-minus), 0.5*(plus+minus));
693 #include "Transverse.h"
701 template <
typename Value>
702 inline LorentzVector<Value>
711 template <
typename Value>
712 inline LorentzVector<Value>
714 Value x = Value(), Value y = Value()) {
723 template <
typename Value>
724 inline LorentzVector<Value>
733 template <
typename OStream,
typename UnitT,
typename Value>
740 template <
typename IStream,
typename UnitT,
typename Value>
750 template <
typename T,
typename U>
751 struct BinaryOpTraits;
755 template <
typename T,
typename U>
768 template <
typename T,
typename U>
LorentzVector< Value > & rotateUz(const Axis &axis)
Rotate the reference frame to a new z-axis.
A 4-component Lorentz vector.
double eta() const
Pseudorapidity of spatial part.
LorentzVector< Value > & boost(Boost b, double gamma=-1.)
Apply boost.
#define ERROR_IF(condition, message)
Debug helper function.
Value perp(const ThreeVector< U > &p) const
Transverse component of the spatial vector with respect to the given axis.
void setRho(Value newRho)
Set new radius.
Value et() const
Transverse energy (signed).
static Dir dir()
Return the direction.
LorentzVector< typename BinaryOpTraits< T, U >::MulT > MulT
The type resulting from multiplication of the template type with itself.
Value m() const
Magnitude (signed) .
LorentzVector< Value > lightConeDir(Value plus, Value minus, Value x=Value(), Value y=Value())
Create a LorentzVector giving its light-cone and transverse components.
LorentzVector< typename BinaryOpTraits< T, U >::MulT > MulT
The type resulting from multiplication of the template type with itself.
Value2 m2() const
Squared magnitude .
bool isNear(const LorentzVector< Value > &w, double epsilon) const
Are two vectors nearby, using Euclidean measure ?
std::complex< double > Complex
ThePEG code should use Complex for all complex scalars.
Value2 m2(const LorentzVector< Value > &a) const
Squared magnitude with another vector.
void ounitstream(OStream &os, const vector< T, Alloc > &v, UT &u)
Ouput a vector of objects with the specified unit.
LorentzVector< Value > lightCone(Value plus, Value minus, Value x, Value y)
Create a LorentzVector giving its light-cone and transverse components.
LorentzVector< Value > & operator=(const LorentzVector< ValueB > &b)
Assignment operator.
LorentzVector< Value > & operator*=(const SpinOneLorentzRotation &m)
Rotate the vector. Resets .
Value et(const ThreeVector< double > &v) const
Transverse energy with respect to the given axis (signed).
Value y() const
The y-component.
double dirTheta(const LorentzVector< Value > &p)
Return the polar angle wrt.
double cosTheta() const
Cosine of the polar angle.
Value2 perp2() const
Squared transverse component of the spatial vector .
ThreeVector< Value > vect() const
Access to the 3-component part.
This is the main namespace within which all identifiers in ThePEG are declared.
double theta() const
Polar angle.
Value2 perp2(const ThreeVector< U > &p) const
Squared transverse component of the spatial vector with respect to the given axis.
double rapidity() const
Rapidity .
Value perp() const
Transverse component of the spatial vector .
BinaryOpTraits< Value, Value >::MulT Value2
Value squared.
LorentzVector< Value > & boost(double bx, double by, double bz, double gamma=-1.)
Apply boost.
double dirCosTheta(const LorentzVector< Value > &p)
Return the cosine of the polar angle wrt.
void iunitstream(IStream &is, vector< T, Alloc > &v, UT &u)
Input a vector of objects with the specified unit.
LorentzVector< Value > & rotateY(double phi)
Apply rotation around the y-axis.
Value2 mt2() const
Transverse mass squared .
static bool pos()
Return true if the direction is positive.
Value minus() const
Returns the negative light-cone component .
LorentzVector< Value > & rotateX(double phi)
Apply rotation around the x-axis.
contains the ThreeVector class.
The SpinOneLorentzRotation class is ...
double angle(const LorentzVector< Value > &w) const
Spatial angle with another vector.
Value mt() const
Transverse mass (signed) .
OUnit< T, UT > ounit(const T &t, const UT &ut)
Helper function creating a OUnit object given an object and a unit.
static bool neg()
Return true if the direction is negative (reversed).
LorentzVector< typename BinaryOpTraits< T, U >::DivT > DivT
The type resulting from division of one template type with another.
LorentzVector< Value > & rotate(double angle, const ThreeVector< U > &axis)
Apply a rotation.
Value dirZ(const LorentzVector< Value > &p)
Return the component along the positive z-axis.
Value plus() const
Returns the positive light-cone component .
Value rho() const
Radius.
Transverse represents the transverse components of a LorentzVector.
Value2 et2() const
Transverse energy squared.
LorentzVector< Value > conjugate() const
The complex conjugate vector.
const double pi
Good old .
LorentzVector< Value > & rotateZ(double phi)
Apply rotation around the z-axis.
Value2 mag2() const
Squared magnitude .
Boost boostVector() const
Boost from reference frame into this vector's rest frame: .
void setVect(const ThreeVector< Value > &p)
Set the 3-component part.
Value dirMinus(const LorentzVector< Value > &p)
Return the negative light-cone component.
ThreeVector< Value > dirBoostVector(const LorentzVector< Value > &p)
Get the boost vector for the LorentzVector.
ThreeVector< double > unit() const
Parallel vector with unit length.
A Direction object can be used to specify that some following operations should be assumed to be perf...
double phi() const
Azimuthal angle.
double rapidity(const Axis &ref) const
Rapidity with respect to another vector.
BinaryOpTraits should be specialized with typdefs called MulT and DivT which gives the type resulting...
Value dirPlus(const LorentzVector< Value > &p)
Return the positive light-cone component.
Value mag() const
Magnitude .
BinaryOpTraits< Value, U >::MulT dot(const LorentzVector< U > &a) const
Dot product with metric .
LorentzVector< Value > & transform(const SpinOneLorentzRotation &m)
Rotate the vector. Resets .
Boost findBoostToCM() const
Boost from reference frame into this vector's rest frame: .
Value2 rho2() const
Radius squared.
Value2 et2(const ThreeVector< double > &v) const
Transverse energy squared with respect to the given axis.
Value x() const
The x-component.
IUnit< T, UT > iunit(T &t, const UT &ut)
Helper function creating a IUnit object given an object and a unit.