00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef ThePEG_LorentzVector_H
00010 #define ThePEG_LorentzVector_H
00011
00019 #include "LorentzVector.fh"
00020 #include "Transverse.h"
00021 #include "ThePEG/Utilities/Direction.h"
00022 #include "ThePEG/Utilities/UnitIO.h"
00023 #include "LorentzRotation.h"
00024 #include "ThreeVector.h"
00025
00026 namespace {
00028 #ifdef NDEBUG
00029 inline void errorIf(bool, const std::string &) {}
00030 #else
00031 inline void errorIf(bool condition, const std::string & message) {
00032 if ( condition )
00033 throw ThePEG::Exception(message, ThePEG::Exception::eventerror);
00034 }
00035 #endif
00036 }
00037
00038 namespace ThePEG {
00039
00040 template <typename Value> class LorentzVector;
00041
00048 template <typename Value> class LorentzVector
00049 {
00050 private:
00052 typedef typename BinaryOpTraits<Value,Value>::MulT Value2;
00053
00054 public:
00057 LorentzVector()
00058 : theX(), theY(), theZ(), theT() {}
00059
00060 LorentzVector(Value x, Value y, Value z, Value t)
00061 : theX(x), theY(y), theZ(z), theT(t) {}
00062
00063 LorentzVector(const ThreeVector<Value> & v, Value t)
00064 : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {}
00065
00066 template<typename U>
00067 LorentzVector(const LorentzVector<U> & v)
00068 : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {}
00070
00072 template <typename ValueB>
00073 LorentzVector<Value> & operator=(const LorentzVector<ValueB> & b) {
00074 setX(b.x());
00075 setY(b.y());
00076 setZ(b.z());
00077 setT(b.t());
00078 return *this;
00079 }
00080
00081 public:
00083
00084 Value x() const { return theX; }
00085 Value y() const { return theY; }
00086 Value z() const { return theZ; }
00087 Value t() const { return theT; }
00088 Value e() const { return t(); }
00090
00092
00093 void setX(Value x) { theX = x; }
00094 void setY(Value y) { theY = y; }
00095 void setZ(Value z) { theZ = z; }
00096 void setT(Value t) { theT = t; }
00097 void setE(Value e) { setT(e); }
00099
00100 public:
00102 ThreeVector<Value> vect() const {
00103 return ThreeVector<Value>(x(),y(),z());
00104 }
00105
00107 operator ThreeVector<Value>() const { return vect(); }
00108
00110 void setVect(const ThreeVector<Value> & p) {
00111 theX = p.x();
00112 theY = p.y();
00113 theZ = p.z();
00114 }
00115
00116 public:
00118 LorentzVector<Value> conjugate() const
00119 {
00120 return LorentzVector<Value>(conj(x()),conj(y()),conj(z()),conj(t()));
00121 }
00122
00124 Value2 m2() const
00125 {
00126 return (t()-z())*(t()+z()) - sqr(x()) - sqr(y());
00127 }
00128
00130 Value2 m2(const LorentzVector<Value> & a) const {
00131 Value tt(a.t()+t()),zz(a.z()+z());
00132 return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y());
00133 }
00134
00136 Value m() const
00137 {
00138 Value2 tmp = m2();
00139 return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
00140 }
00141
00143 Value2 mt2() const { return (t()-z())*(t()+z()); }
00144
00146 Value mt() const
00147 {
00148 Value2 tmp = mt2();
00149 return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
00150 }
00151
00153 Value2 perp2() const { return sqr(x()) + sqr(y()); }
00154
00156 Value perp() const { return sqrt(perp2()); }
00157
00162 template <typename U>
00163 Value2 perp2(const ThreeVector<U> & p) const
00164 {
00165 return vect().perp2(p);
00166 }
00167
00172 template <typename U>
00173 Value perp(const ThreeVector<U> & p) const
00174 {
00175 return vect().perp(p);
00176 }
00177
00179 Value2 et2() const
00180 {
00181 Value2 pt2 = vect().perp2();
00182 return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z());
00183 }
00184
00186 Value et() const
00187 {
00188 Value2 etet = et2();
00189 return e() < Value() ? -sqrt(etet) : sqrt(etet);
00190 }
00191
00193 Value2 et2(const ThreeVector<double> & v) const
00194 {
00195 Value2 pt2 = vect().perp2(v);
00196 Value pv = vect().dot(v.unit());
00197 return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv);
00198 }
00199
00201 Value et(const ThreeVector<double> & v) const
00202 {
00203 Value2 etet = et2(v);
00204 return e() < Value() ? -sqrt(etet) : sqrt(etet);
00205 }
00206
00208
00209
00210 Value2 rho2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
00211
00213 Value rho() const { return sqrt(rho2()); }
00214
00216 void setRho(Value newRho)
00217 {
00218 Value oldRho = rho();
00219 if (oldRho == Value())
00220 return;
00221 double factor = newRho / oldRho;
00222 setX(x()*factor);
00223 setY(y()*factor);
00224 setZ(z()*factor);
00225 }
00226
00228 double theta() const
00229 {
00230 assert(!(x() == Value() && y() == Value() && z() == Value()));
00231 return atan2(perp(),z());
00232 }
00233
00235 double cosTheta() const
00236 {
00237 Value ptot = rho();
00238 assert( ptot > Value() );
00239 return z() / ptot;
00240 }
00241
00243 double phi() const
00244 {
00245 if ( x() == Value() && y() == Value() )
00246 return 0.;
00247
00248 return atan2(y(),x()) ;
00249 }
00251
00253 double eta() const {
00254 static const string
00255 em("Pseudorapidity for 3-vector along z-axis undefined.");
00256 Value m = rho();
00257 if ( m == Value() ) return 0.0;
00258 errorIf(m == z() || m == -z(),em);
00259 return 0.5 * log( (m+z()) / (m-z()) );
00260 }
00261
00263 double angle(const LorentzVector<Value> & w) const
00264 {
00265 return vect().angle(w.vect());
00266 }
00267
00269 double rapidity() const {
00270 static const string em1("rapidity for 4-vector with |E| = |Pz| -- infinite result");
00271 static const string em2("rapidity for spacelike 4-vector with |E| < |Pz| -- undefined");
00272 errorIf(abs(t()) == abs(z()),em1);
00273 errorIf(abs(t()) < abs(z()) ,em2);
00274 double q = (t() + z()) / (t() - z());
00275 return 0.5 * log(q);
00276 }
00277
00279 double rapidity(const Axis & ref) const {
00280 static const string em1("A zero vector used as reference to LorentzVector rapidity");
00281 static const string em2("rapidity for 4-vector with |E| = |Pu| -- infinite result");
00282 static const string em3("rapidity for spacelike 4-vector with |E|<|P*ref| undefined");
00283 double r = ref.mag2();
00284 errorIf(r == 0,em1);
00285 Value vdotu = vect().dot(ref)/sqrt(r);
00286 errorIf(abs(t()) == abs(vdotu),em2);
00287 errorIf(abs(t()) < abs(vdotu),em3);
00288 double q = (t() + vdotu) / (t() - vdotu);
00289 return 0.5 * log(q);
00290 }
00291
00296 Boost boostVector() const {
00297 static const string em1("boostVector computed for LorentzVector with t=0"
00298 " -- infinite result");
00299 static const string em2("boostVector computed for a non-timelike LorentzVector");
00300 if (t() == Value()) {
00301 if (rho2() == Value2())
00302 return Boost();
00303 else
00304 errorIf(true,em1);
00305 }
00306
00307 errorIf(m2() <= Value2(),em2);
00308 return vect() * (1./t());
00309 }
00310
00315 Boost findBoostToCM() const
00316 {
00317 return -boostVector();
00318 }
00319
00321 Value plus() const { return t() + z(); }
00323 Value minus() const { return t() - z(); }
00324
00326 bool isNear(const LorentzVector<Value> & w, double epsilon) const
00327 {
00328 Value2 limit = abs(vect().dot(w.vect()));
00329 limit += 0.25 * sqr( t() + w.t() );
00330 limit *= sqr(epsilon);
00331 Value2 delta = (vect() - w.vect()).mag2();
00332 delta += sqr( t() - w.t() );
00333 return (delta <= limit);
00334 }
00335
00337 LorentzVector<Value> & transform(const SpinOneLorentzRotation & m)
00338 {
00339 return *this = m.operator*(*this);
00340 }
00341
00343 LorentzVector<Value> & operator*=(const SpinOneLorentzRotation & m)
00344 {
00345 return transform(m);
00346 }
00347
00349 template <typename U>
00350 typename BinaryOpTraits<Value,U>::MulT
00351 dot(const LorentzVector<U> & a) const
00352 {
00353 return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() );
00354 }
00355
00356
00357 public:
00358
00370 LorentzVector<Value> &
00371 boost(double bx, double by, double bz, double gamma=-1.)
00372 {
00373 double b2 = bx*bx + by*by + bz*bz;
00374 if(gamma < 0.)
00375 gamma = 1.0 / sqrt(1.0 - b2);
00376 Value bp = bx*x() + by*y() + bz*z();
00377 double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;
00378
00379 setX(x() + gamma2*bp*bx + gamma*bx*t());
00380 setY(y() + gamma2*bp*by + gamma*by*t());
00381 setZ(z() + gamma2*bp*bz + gamma*bz*t());
00382 setT(gamma*(t() + bp));
00383 return *this;
00384 }
00385
00396 LorentzVector<Value> & boost(Boost b, double gamma=-1.) {
00397 return boost(b.x(), b.y(), b.z(),gamma);
00398 }
00399
00405 LorentzVector<Value> & rotateX (double phi) {
00406 double sinphi = sin(phi);
00407 double cosphi = cos(phi);
00408 Value ty = y() * cosphi - z() * sinphi;
00409 theZ = z() * cosphi + y() * sinphi;
00410 theY = ty;
00411 return *this;
00412 }
00413
00419 LorentzVector<Value> & rotateY (double phi) {
00420 double sinphi = sin(phi);
00421 double cosphi = cos(phi);
00422 Value tz = z() * cosphi - x() * sinphi;
00423 theX = x() * cosphi + z() * sinphi;
00424 theZ = tz;
00425 return *this;
00426 }
00427
00433 LorentzVector<Value> & rotateZ (double phi) {
00434 double sinphi = sin(phi);
00435 double cosphi = cos(phi);
00436 Value tx = x() * cosphi - y() * sinphi;
00437 theY = y() * cosphi + x() * sinphi;
00438 theX = tx;
00439 return *this;
00440 }
00441
00445 LorentzVector<Value> & rotateUz (const Axis & axis) {
00446 Axis ax = axis.unit();
00447 double u1 = ax.x();
00448 double u2 = ax.y();
00449 double u3 = ax.z();
00450 double up = u1*u1 + u2*u2;
00451 if (up>0) {
00452 up = sqrt(up);
00453 Value px = x(), py = y(), pz = z();
00454 setX( (u1*u3*px - u2*py)/up + u1*pz );
00455 setY( (u2*u3*px + u1*py)/up + u2*pz );
00456 setZ( -up*px + u3*pz );
00457 }
00458 else if (u3 < 0.) {
00459 setX(-x());
00460 setZ(-z());
00461 }
00462 return *this;
00463 }
00464
00465 public:
00467
00468 template <typename ValueB>
00469 LorentzVector<Value> & operator+=(const LorentzVector<ValueB> & a) {
00470 theX += a.x();
00471 theY += a.y();
00472 theZ += a.z();
00473 theT += a.t();
00474 return *this;
00475 }
00476
00477 template <typename ValueB>
00478 LorentzVector<Value> & operator-=(const LorentzVector<ValueB> & a) {
00479 theX -= a.x();
00480 theY -= a.y();
00481 theZ -= a.z();
00482 theT -= a.t();
00483 return *this;
00484 }
00485
00486 LorentzVector<Value> & operator*=(double a) {
00487 theX *= a;
00488 theY *= a;
00489 theZ *= a;
00490 theT *= a;
00491 return *this;
00492 }
00493
00494 LorentzVector<Value> & operator/=(double a) {
00495 theX /= a;
00496 theY /= a;
00497 theZ /= a;
00498 theT /= a;
00499 return *this;
00500 }
00502
00503 private:
00505
00506 Value theX;
00507 Value theY;
00508 Value theZ;
00509 Value theT;
00511 };
00512
00514
00515 template <typename Value>
00516 inline LorentzVector<double>
00517 operator/(const LorentzVector<Value> & v, Value a) {
00518 return LorentzVector<double>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
00519 }
00520
00521 template <typename Value>
00522 inline LorentzVector<Value> operator-(const LorentzVector<Value> & v) {
00523 return LorentzVector<Value>(-v.x(),-v.y(),-v.z(),-v.t());
00524 }
00525
00526 template <typename ValueA, typename ValueB>
00527 inline LorentzVector<ValueA>
00528 operator+(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
00529 return a += b;
00530 }
00531
00532 template <typename ValueA, typename ValueB>
00533 inline LorentzVector<ValueA>
00534 operator-(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
00535 return a -= b;
00536 }
00537
00538 template <typename Value>
00539 inline LorentzVector<Value>
00540 operator*(const LorentzVector<Value> & a, double b) {
00541 return LorentzVector<Value>(a.x()*b, a.y()*b, a.z()*b, a.t()*b);
00542 }
00543
00544 template <typename Value>
00545 inline LorentzVector<Value>
00546 operator*(double b, LorentzVector<Value> a) {
00547 return a *= b;
00548 }
00549
00550 template <typename ValueA, typename ValueB>
00551 inline
00552 LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
00553 operator*(ValueB a, const LorentzVector<ValueA> & v) {
00554 typedef typename BinaryOpTraits<ValueB,ValueA>::MulT ResultT;
00555 return LorentzVector<ResultT>(a*v.x(), a*v.y(), a*v.z(), a*v.t());
00556 }
00557
00558 template <typename ValueA, typename ValueB>
00559 inline
00560 LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
00561 operator*(const LorentzVector<ValueA> & v, ValueB b) {
00562 return b*v;
00563 }
00564
00565 template <typename ValueA, typename ValueB>
00566 inline
00567 LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::DivT>
00568 operator/(const LorentzVector<ValueA> & v, ValueB b) {
00569 typedef typename BinaryOpTraits<ValueA,ValueB>::DivT ResultT;
00570 return LorentzVector<ResultT>(v.x()/b, v.y()/b, v.z()/b, v.t()/b);
00571 }
00573
00575
00576 template <typename ValueA, typename ValueB>
00577 inline typename BinaryOpTraits<ValueA,ValueB>::MulT
00578 operator*(const LorentzVector<ValueA> & a, const LorentzVector<ValueB> & b) {
00579 return a.dot(b);
00580 }
00581
00582 template <typename Value>
00583 inline typename BinaryOpTraits<Value,Value>::MulT
00584 operator*(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
00585 return a.dot(b);
00586 }
00588
00590 template <typename Value>
00591 inline bool
00592 operator==(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
00593 return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t();
00594 }
00595
00597 inline ostream & operator<< (ostream & os, const LorentzVector<double> & v) {
00598 return os << "(" << v.x() << "," << v.y() << "," << v.z()
00599 << ";" << v.t() << ")";
00600 }
00601
00604 template <typename Value>
00605 inline double dirPlus(const LorentzVector<Value> & p) {
00606 return Direction<0>::pos()? p.plus(): p.minus();
00607 }
00608
00611 template <typename Value>
00612 inline double dirMinus(const LorentzVector<Value> & p) {
00613 return Direction<0>::neg()? p.plus(): p.minus();
00614 }
00615
00618 template <typename Value>
00619 inline double dirZ(const LorentzVector<Value> & p) {
00620 return Direction<0>::dir()*p.z();
00621 }
00622
00625 template <typename Value>
00626 inline double dirTheta(const LorentzVector<Value> & p) {
00627 return Direction<0>::pos()? p.theta(): Constants::pi - p.theta();
00628 }
00629
00632 template <typename Value>
00633 inline double dirCosTheta(const LorentzVector<Value> & p) {
00634 return Direction<0>::pos()? p.cosTheta(): -p.cosTheta();
00635 }
00636
00639 template <typename Value>
00640 inline ThreeVector<Value> dirBoostVector(const LorentzVector<Value> & p) {
00641 ThreeVector<Value> b(p.boostVector());
00642 if ( Direction<0>::neg() ) b.setZ(-b.z());
00643 return b;
00644 }
00645
00648 template <typename Value>
00649 inline LorentzVector<Value>
00650 lightCone(Value plus, Value minus, Value x, Value y) {
00651 LorentzVector<Value> r(x, y, 0.5*(plus-minus), 0.5*(plus+minus));
00652 return r;
00653 }
00654
00656 template <typename Value>
00657 inline LorentzVector<Value>
00658 lightCone(Value plus, Value minus) {
00659
00660
00661 static const Value zero = Value();
00662 LorentzVector<Value> r(zero, zero,
00663 0.5*(plus-minus), 0.5*(plus+minus));
00664 return r;
00665 }
00666
00669 template <typename Value>
00670 inline LorentzVector<Value>
00671 lightCone(Value plus, Value minus, Transverse<Value> pt) {
00672 LorentzVector<Value> r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus));
00673 return r;
00674 }
00675
00679 template <typename Value>
00680 inline LorentzVector<Value>
00681 lightConeDir(Value plus, Value minus,
00682 Value x = Value(), Value y = Value()) {
00683 LorentzVector<Value> r(x, y, Direction<0>::dir()*0.5*(plus - minus),
00684 0.5*(plus + minus));
00685 return r;
00686 }
00687
00691 template <typename Value>
00692 inline LorentzVector<Value>
00693 lightConeDir(Value plus, Value minus, Transverse<Value> pt) {
00694 LorentzVector<Value> r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus),
00695 0.5*(plus + minus));
00696 return r;
00697
00698 }
00699
00701 template <typename OStream, typename UnitT, typename Value>
00702 void ounitstream(OStream & os, const LorentzVector<Value> & p, UnitT & u) {
00703 os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u)
00704 << ounit(p.e(), u);
00705 }
00706
00708 template <typename IStream, typename UnitT, typename Value>
00709 void iunitstream(IStream & is, LorentzVector<Value> & p, UnitT & u) {
00710 Value x, y, z, e;
00711 is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u);
00712 p = LorentzVector<Value>(x, y, z, e);
00713 }
00714
00715
00717
00718 template <typename T, typename U>
00719 struct BinaryOpTraits;
00723 template <typename T, typename U>
00724 struct BinaryOpTraits<LorentzVector<T>, U> {
00727 typedef LorentzVector<typename BinaryOpTraits<T,U>::MulT> MulT;
00730 typedef LorentzVector<typename BinaryOpTraits<T,U>::DivT> DivT;
00731 };
00732
00736 template <typename T, typename U>
00737 struct BinaryOpTraits<T, LorentzVector<U> > {
00740 typedef LorentzVector<typename BinaryOpTraits<T,U>::MulT> MulT;
00741 };
00743
00744 }
00745
00746
00747 #endif
