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) {}
59 : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {}
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>
531 template <
typename Value>
542 template <
typename Value>
547 template <
typename ValueA,
typename ValueB>
553 template <
typename ValueA,
typename ValueB>
559 template <
typename Value>
565 template <
typename Value>
571 template <
typename ValueA,
typename ValueB>
579 template <
typename ValueA,
typename ValueB>
586 template <
typename ValueA,
typename ValueB>
597 template <
typename ValueA,
typename ValueB>
603 template <
typename Value>
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>
677 template <
typename Value>
682 static const Value zero = Value();
684 0.5*(plus-minus), 0.5*(plus+minus));
693 #include "Transverse.h" 701 template <
typename Value>
711 template <
typename Value>
714 Value x = Value(), Value y = Value()) {
723 template <
typename Value>
733 template <
typename OStream,
typename UnitT,
typename Value>
740 template <
typename IStream,
typename UnitT,
typename Value>
750 template <
typename T,
typename U>
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.
double theta() const
Polar angle.
A 4-component Lorentz vector.
LorentzVector< Value > & boost(Boost b, double gamma=-1.)
Apply boost.
BinaryOpTraits< Value, U >::MulT dot(const LorentzVector< U > &a) const
Dot product with metric .
#define ERROR_IF(condition, message)
Debug helper function.
void setRho(Value newRho)
Set new radius.
static Dir dir()
Return the direction.
LorentzVector< typename BinaryOpTraits< T, U >::MulT > MulT
The type resulting from multiplication of the template type with itself.
LorentzVector< Value > lightConeDir(Value plus, Value minus, Value x=Value(), Value y=Value())
Create a LorentzVector giving its light-cone and transverse components.
Value m() const
Magnitude (signed) .
double rapidity(const Axis &ref) const
Rapidity with respect to another vector.
LorentzVector< typename BinaryOpTraits< T, U >::MulT > MulT
The type resulting from multiplication of the template type with itself.
Value mag() const
Magnitude .
Value2 rho2() const
Radius squared.
std::complex< double > Complex
ThePEG code should use Complex for all complex scalars.
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.
double eta() const
Pseudorapidity of spatial part.
Value perp() const
Transverse component of the spatial vector .
Value2 mag2() const
Squared magnitude .
LorentzVector< Value > & operator*=(const SpinOneLorentzRotation &m)
Rotate the vector. Resets .
double dirTheta(const LorentzVector< Value > &p)
Return the polar angle wrt.
This is the main namespace within which all identifiers in ThePEG are declared.
Value2 perp2() const
Squared 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.
static bool pos()
Return true if the direction is positive.
double phi() const
Azimuthal angle.
Boost boostVector() const
Boost from reference frame into this vector's rest frame: .
LorentzVector< Value > & rotateX(double phi)
Apply rotation around the x-axis.
Value x() const
The x-component.
Value plus() const
Returns the positive light-cone component .
contains the ThreeVector class.
The SpinOneLorentzRotation class is ...
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.
Value y() const
The y-component.
ThreeVector< double > unit() const
Parallel vector with unit length.
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 mt() const
Transverse mass (signed) .
Transverse represents the transverse components of a LorentzVector.
Boost findBoostToCM() const
Boost from reference frame into this vector's rest frame: .
const double pi
Good old .
LorentzVector< Value > & rotateZ(double phi)
Apply rotation around the z-axis.
bool isNear(const LorentzVector< Value > &w, double epsilon) const
Are two vectors nearby, using Euclidean measure ?
double cosTheta() const
Cosine of the polar angle.
Value2 mt2() const
Transverse mass squared .
Value et(const ThreeVector< double > &v) const
Transverse energy with respect to the given axis (signed).
double rapidity() const
Rapidity .
ThreeVector< Value > vect() const
Access to the 3-component part.
Value2 m2(const LorentzVector< Value > &a) const
Squared magnitude with another vector.
Value2 et2() const
Transverse energy squared.
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.
A Direction object can be used to specify that some following operations should be assumed to be perf...
Value2 m2() const
Squared magnitude .
LorentzVector< Value > conjugate() const
The complex conjugate vector.
Value2 et2(const ThreeVector< double > &v) const
Transverse energy squared with respect to the given axis.
Value perp(const ThreeVector< U > &p) const
Transverse component of the spatial vector with respect to the given axis.
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.
Value2 perp2(const ThreeVector< U > &p) const
Squared transverse component of the spatial vector with respect to the given axis.
LorentzVector< Value > & transform(const SpinOneLorentzRotation &m)
Rotate the vector. Resets .
double angle(const LorentzVector< Value > &w) const
Spatial angle with another vector.
Value minus() const
Returns the negative light-cone component .
Value et() const
Transverse energy (signed).
IUnit< T, UT > iunit(T &t, const UT &ut)
Helper function creating a IUnit object given an object and a unit.