ThePEG  1.8.0
LorentzVector.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // LorentzVector.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 2006-2011 David Grellscheid, Leif Lonnblad
5 //
6 // ThePEG is licenced under version 2 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef ThePEG_LorentzVector_H
10 #define ThePEG_LorentzVector_H
11 
19 #include "LorentzVector.fh"
20 #include "ThePEG/Utilities/Direction.h"
21 #include "ThePEG/Utilities/UnitIO.h"
22 #include "LorentzRotation.h"
23 #include "ThreeVector.h"
24 
26 #ifdef NDEBUG
27 #define ERROR_IF(condition,message) if (false) {}
28 #else
29 #define ERROR_IF(condition,message) \
30  if ( (condition) ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror)
31 #endif
32 
33 namespace ThePEG {
34 
35 template <typename Value> class LorentzVector;
36 
43 template <typename Value> class LorentzVector
44 {
45 private:
47  typedef typename BinaryOpTraits<Value,Value>::MulT Value2;
48 
49 public:
52  LorentzVector()
53  : theX(), theY(), theZ(), theT() {}
54 
55  LorentzVector(Value x, Value y, Value z, Value t)
56  : theX(x), theY(y), theZ(z), theT(t) {}
57 
58  LorentzVector(const ThreeVector<Value> & v, Value t)
59  : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {}
60 
61  template<typename U>
62  LorentzVector(const LorentzVector<U> & v)
63  : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {}
65 
67  template <typename ValueB>
69  setX(b.x());
70  setY(b.y());
71  setZ(b.z());
72  setT(b.t());
73  return *this;
74  }
75 
76 public:
78 
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(); }
85 
87 
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); }
94 
95 public:
98  return ThreeVector<Value>(x(),y(),z());
99  }
100 
102  operator ThreeVector<Value>() const { return vect(); }
103 
105  void setVect(const ThreeVector<Value> & p) {
106  theX = p.x();
107  theY = p.y();
108  theZ = p.z();
109  }
110 
111 public:
114  {
115  return LorentzVector<Value>(conj(x()),conj(y()),conj(z()),conj(t()));
116  }
117 
119  Value2 m2() const
120  {
121  return (t()-z())*(t()+z()) - sqr(x()) - sqr(y());
122  }
123 
125  Value2 m2(const LorentzVector<Value> & a) const {
126  Value tt(a.t()+t()),zz(a.z()+z());
127  return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y());
128  }
129 
131  Value m() const
132  {
133  Value2 tmp = m2();
134  return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
135  }
136 
138  Value2 mt2() const { return (t()-z())*(t()+z()); }
139 
141  Value mt() const
142  {
143  Value2 tmp = mt2();
144  return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
145  }
146 
148  Value2 perp2() const { return sqr(x()) + sqr(y()); }
149 
151  Value perp() const { return sqrt(perp2()); }
152 
157  template <typename U>
158  Value2 perp2(const ThreeVector<U> & p) const
159  {
160  return vect().perp2(p);
161  }
162 
167  template <typename U>
168  Value perp(const ThreeVector<U> & p) const
169  {
170  return vect().perp(p);
171  }
172 
174  Value2 et2() const
175  {
176  Value2 pt2 = vect().perp2();
177  return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z());
178  }
179 
181  Value et() const
182  {
183  Value2 etet = et2();
184  return e() < Value() ? -sqrt(etet) : sqrt(etet);
185  }
186 
188  Value2 et2(const ThreeVector<double> & v) const
189  {
190  Value2 pt2 = vect().perp2(v);
191  Value pv = vect().dot(v.unit());
192  return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv);
193  }
194 
196  Value et(const ThreeVector<double> & v) const
197  {
198  Value2 etet = et2(v);
199  return e() < Value() ? -sqrt(etet) : sqrt(etet);
200  }
201 
203 
204 
205  Value2 rho2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
206 
208  Value rho() const { return sqrt(rho2()); }
209 
211  void setRho(Value newRho)
212  {
213  Value oldRho = rho();
214  if (oldRho == Value())
215  return;
216  double factor = newRho / oldRho;
217  setX(x()*factor);
218  setY(y()*factor);
219  setZ(z()*factor);
220  }
221 
223  double theta() const
224  {
225  assert(!(x() == Value() && y() == Value() && z() == Value()));
226  return atan2(perp(),z());
227  }
228 
230  double cosTheta() const
231  {
232  Value ptot = rho();
233  assert( ptot > Value() );
234  return z() / ptot;
235  }
236 
238  double phi() const {
239  return atan2(y(),x()) ;
240  }
242 
244  double eta() const {
245  Value m = rho();
246  if ( m == Value() ) return 0.0;
247  ERROR_IF(m == z() || m == -z(),
248  "Pseudorapidity for 3-vector along z-axis undefined.");
249  return 0.5 * log( (m+z()) / (m-z()) );
250  }
251 
253  double angle(const LorentzVector<Value> & w) const
254  {
255  return vect().angle(w.vect());
256  }
257 
259  double rapidity() const {
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());
263  return 0.5 * log(q);
264  }
265 
267  double rapidity(const Axis & ref) const {
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);
274  return 0.5 * log(q);
275  }
276 
281  Boost boostVector() const {
282  if (t() == Value()) {
283  if (rho2() == Value2())
284  return Boost();
285  else
286  ERROR_IF(true,"boostVector computed for LorentzVector with t=0 -- infinite result");
287  }
288  // result will make analytic sense but is physically meaningless
289  ERROR_IF(m2() <= Value2(),"boostVector computed for a non-timelike LorentzVector");
290  return vect() * (1./t());
291  }
292 
298  {
299  return -boostVector();
300  }
301 
303  Value plus() const { return t() + z(); }
305  Value minus() const { return t() - z(); }
306 
308  bool isNear(const LorentzVector<Value> & w, double epsilon) const
309  {
310  Value2 limit = abs(vect().dot(w.vect()));
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);
316  }
317 
320  {
321  return *this = m.operator*(*this);
322  }
323 
326  {
327  return transform(m);
328  }
329 
331  template <typename U>
332  typename BinaryOpTraits<Value,U>::MulT
333  dot(const LorentzVector<U> & a) const
334  {
335  return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() );
336  }
337 
338 
339 public:
340 
353  boost(double bx, double by, double bz, double gamma=-1.)
354  {
355  double b2 = bx*bx + by*by + bz*bz;
356  if(gamma < 0.)
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;
360 
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));
365  return *this;
366  }
367 
378  LorentzVector<Value> & boost(Boost b, double gamma=-1.) {
379  return boost(b.x(), b.y(), b.z(),gamma);
380  }
381 
388  double sinphi = sin(phi);
389  double cosphi = cos(phi);
390  Value ty = y() * cosphi - z() * sinphi;
391  theZ = z() * cosphi + y() * sinphi;
392  theY = ty;
393  return *this;
394  }
395 
402  double sinphi = sin(phi);
403  double cosphi = cos(phi);
404  Value tz = z() * cosphi - x() * sinphi;
405  theX = x() * cosphi + z() * sinphi;
406  theZ = tz;
407  return *this;
408  }
409 
416  double sinphi = sin(phi);
417  double cosphi = cos(phi);
418  Value tx = x() * cosphi - y() * sinphi;
419  theY = y() * cosphi + x() * sinphi;
420  theX = tx;
421  return *this;
422  }
423 
427  LorentzVector<Value> & rotateUz (const Axis & axis) {
428  Axis ax = axis.unit();
429  double u1 = ax.x();
430  double u2 = ax.y();
431  double u3 = ax.z();
432  double up = u1*u1 + u2*u2;
433  if (up>0) {
434  up = sqrt(up);
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 );
439  }
440  else if (u3 < 0.) {
441  setX(-x());
442  setZ(-z());
443  }
444  return *this;
445  }
446 
452  template <typename U>
454  if (angle == 0.0)
455  return *this;
456  const U ll = axis.mag();
457  assert( ll > U() );
458 
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();
462 
463  setX((ca+(1-ca)*dx*dx) * xx
464  +((1-ca)*dx*dy-sa*dz) * yy
465  +((1-ca)*dx*dz+sa*dy) * zz
466  );
467  setY(((1-ca)*dy*dx+sa*dz) * xx
468  +(ca+(1-ca)*dy*dy) * yy
469  +((1-ca)*dy*dz-sa*dx) * zz
470  );
471  setZ(((1-ca)*dz*dx-sa*dy) * xx
472  +((1-ca)*dz*dy+sa*dx) * yy
473  +(ca+(1-ca)*dz*dz) * zz
474  );
475  return *this;
476  }
477 
478 
479 
480 
481 public:
483 
484  template <typename ValueB>
485  LorentzVector<Value> & operator+=(const LorentzVector<ValueB> & a) {
486  theX += a.x();
487  theY += a.y();
488  theZ += a.z();
489  theT += a.t();
490  return *this;
491  }
492 
493  template <typename ValueB>
494  LorentzVector<Value> & operator-=(const LorentzVector<ValueB> & a) {
495  theX -= a.x();
496  theY -= a.y();
497  theZ -= a.z();
498  theT -= a.t();
499  return *this;
500  }
501 
502  LorentzVector<Value> & operator*=(double a) {
503  theX *= a;
504  theY *= a;
505  theZ *= a;
506  theT *= a;
507  return *this;
508  }
509 
510  LorentzVector<Value> & operator/=(double a) {
511  theX /= a;
512  theY /= a;
513  theZ /= a;
514  theT /= a;
515  return *this;
516  }
518 
519 private:
521 
522  Value theX;
523  Value theY;
524  Value theZ;
525  Value theT;
527 };
528 
530 
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);
535 }
536 
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);
540 }
541 
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());
545 }
546 
547 template <typename ValueA, typename ValueB>
548 inline LorentzVector<ValueA>
549 operator+(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
550  return a += b;
551 }
552 
553 template <typename ValueA, typename ValueB>
554 inline LorentzVector<ValueA>
555 operator-(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
556  return a -= b;
557 }
558 
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);
563 }
564 
565 template <typename Value>
566 inline LorentzVector<Value>
567 operator*(double b, LorentzVector<Value> a) {
568  return a *= b;
569 }
570 
571 template <typename ValueA, typename ValueB>
572 inline
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());
577 }
578 
579 template <typename ValueA, typename ValueB>
580 inline
581 LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
582 operator*(const LorentzVector<ValueA> & v, ValueB b) {
583  return b*v;
584 }
585 
586 template <typename ValueA, typename ValueB>
587 inline
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);
592 }
594 
596 
597 template <typename ValueA, typename ValueB>
598 inline typename BinaryOpTraits<ValueA,ValueB>::MulT
599 operator*(const LorentzVector<ValueA> & a, const LorentzVector<ValueB> & b) {
600  return a.dot(b);
601 }
602 
603 template <typename Value>
604 inline typename BinaryOpTraits<Value,Value>::MulT
605 operator*(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
606  return a.dot(b);
607 }
609 
611 template <typename Value>
612 inline bool
613 operator==(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
614  return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t();
615 }
616 
618 inline ostream & operator<< (ostream & os, const LorentzVector<double> & v) {
619  return os << "(" << v.x() << "," << v.y() << "," << v.z()
620  << ";" << v.t() << ")";
621 }
622 
625 template <typename Value>
626 inline Value dirPlus(const LorentzVector<Value> & p) {
627  return Direction<0>::pos()? p.plus(): p.minus();
628 }
629 
632 template <typename Value>
633 inline Value dirMinus(const LorentzVector<Value> & p) {
634  return Direction<0>::neg()? p.plus(): p.minus();
635 }
636 
639 template <typename Value>
640 inline Value dirZ(const LorentzVector<Value> & p) {
641  return Direction<0>::dir()*p.z();
642 }
643 
646 template <typename Value>
647 inline double dirTheta(const LorentzVector<Value> & p) {
648  return Direction<0>::pos()? p.theta(): Constants::pi - p.theta();
649 }
650 
653 template <typename Value>
654 inline double dirCosTheta(const LorentzVector<Value> & p) {
655  return Direction<0>::pos()? p.cosTheta(): -p.cosTheta();
656 }
657 
660 template <typename Value>
663  if ( Direction<0>::neg() ) b.setZ(-b.z());
664  return b;
665 }
666 
669 template <typename Value>
670 inline LorentzVector<Value>
671 lightCone(Value plus, Value minus, Value x, Value y) {
672  LorentzVector<Value> r(x, y, 0.5*(plus-minus), 0.5*(plus+minus));
673  return r;
674 }
675 
677 template <typename Value>
678 inline LorentzVector<Value>
679 lightCone(Value plus, Value minus) {
680  // g++-3.3 has a problem with using Value() directly
681  // gcc-bug c++/3650, fixed in 3.4
682  static const Value zero = Value();
683  LorentzVector<Value> r(zero, zero,
684  0.5*(plus-minus), 0.5*(plus+minus));
685  return r;
686 }
687 
688 }
689 
690 
691 // delayed header inclusion to break inclusion loop:
692 // LorentzVec -> Transverse -> Lorentz5Vec -> LorentzVec
693 #include "Transverse.h"
694 
695 
696 
697 namespace ThePEG {
698 
701 template <typename Value>
702 inline LorentzVector<Value>
703 lightCone(Value plus, Value minus, Transverse<Value> pt) {
704  LorentzVector<Value> r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus));
705  return r;
706 }
707 
711 template <typename Value>
712 inline LorentzVector<Value>
713 lightConeDir(Value plus, Value minus,
714  Value x = Value(), Value y = Value()) {
715  LorentzVector<Value> r(x, y, Direction<0>::dir()*0.5*(plus - minus),
716  0.5*(plus + minus));
717  return r;
718 }
719 
723 template <typename Value>
724 inline LorentzVector<Value>
725 lightConeDir(Value plus, Value minus, Transverse<Value> pt) {
726  LorentzVector<Value> r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus),
727  0.5*(plus + minus));
728  return r;
729 
730 }
731 
733 template <typename OStream, typename UnitT, typename Value>
734 void ounitstream(OStream & os, const LorentzVector<Value> & p, UnitT & u) {
735  os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u)
736  << ounit(p.e(), u);
737 }
738 
740 template <typename IStream, typename UnitT, typename Value>
741 void iunitstream(IStream & is, LorentzVector<Value> & p, UnitT & u) {
742  Value x, y, z, e;
743  is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u);
744  p = LorentzVector<Value>(x, y, z, e);
745 }
746 
747 
749 
750 template <typename T, typename U>
751 struct BinaryOpTraits;
755 template <typename T, typename U>
756 struct BinaryOpTraits<LorentzVector<T>, U> {
763 };
764 
768 template <typename T, typename U>
769 struct BinaryOpTraits<T, LorentzVector<U> > {
773 };
775 
776 }
777 
778 #undef ERROR_IF
779 #endif /* ThePEG_LorentzVector_H */