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)
35 template <
typename Value>
class LorentzVector;
47 typedef typename BinaryOpTraits<Value,Value>::MulT
Value2;
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));
144 return tmp <
Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
157 template <
typename U>
160 return vect().perp2(p);
167 template <
typename U>
170 return vect().perp(p);
177 return pt2 ==
Value2() ?
Value2() : e()*e() * pt2/(pt2+z()*z());
184 return e() < Value() ? -sqrt(etet) : sqrt(etet);
192 return pt2 ==
Value2() ?
Value2() : e()*e() * pt2/(pt2+pv*pv);
199 return e() < Value() ? -sqrt(etet) : sqrt(etet);
205 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);
314 delta += sqr( t() - w.t() );
315 return (delta <= limit);
321 return *
this = m.operator*(*this);
331 template <
typename U>
332 typename BinaryOpTraits<Value,U>::MulT
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>