00001 #ifndef RIVET_MATH_VECTOR4
00002 #define RIVET_MATH_VECTOR4
00003
00004 #include "Rivet/Math/MathHeader.hh"
00005 #include "Rivet/Math/MathUtils.hh"
00006 #include "Rivet/Math/VectorN.hh"
00007 #include "Rivet/Math/Vector3.hh"
00008
00009 namespace Rivet {
00010
00011
00012 class FourVector;
00013 class FourMomentum;
00014 class LorentzTransform;
00015 typedef FourVector Vector4;
00016 FourVector transform(const LorentzTransform& lt, const FourVector& v4);
00017
00018
00019 class FourVector : public Vector<4> {
00020 friend FourVector multiply(const double a, const FourVector& v);
00021 friend FourVector multiply(const FourVector& v, const double a);
00022 friend FourVector add(const FourVector& a, const FourVector& b);
00023 friend FourVector transform(const LorentzTransform& lt, const FourVector& v4);
00024
00025 public:
00026 FourVector() : Vector<4>() { }
00027
00028 template<typename V4>
00029 FourVector(const V4& other) {
00030 this->setT(other.t());
00031 this->setX(other.x());
00032 this->setY(other.y());
00033 this->setZ(other.z());
00034 }
00035
00036 FourVector(const Vector<4>& other)
00037 : Vector<4>(other) { }
00038
00039 FourVector(const double t, const double x, const double y, const double z) {
00040 this->setT(t);
00041 this->setX(x);
00042 this->setY(y);
00043 this->setZ(z);
00044 }
00045
00046 virtual ~FourVector() { }
00047
00048 public:
00049 double t() const { return get(0); }
00050 double x() const { return get(1); }
00051 double y() const { return get(2); }
00052 double z() const { return get(3); }
00053 FourVector& setT(const double t) { set(0, t); return *this; }
00054 FourVector& setX(const double x) { set(1, x); return *this; }
00055 FourVector& setY(const double y) { set(2, y); return *this; }
00056 FourVector& setZ(const double z) { set(3, z); return *this; }
00057
00058 double invariant() const {
00059 return t()*t() - x()*x() - y()*y() - z()*z();
00060 }
00061
00062 double angle(const FourVector& v) const {
00063 return vector3().angle( v.vector3() );
00064 }
00065
00066 double angle(const Vector3& v3) const {
00067 return vector3().angle(v3);
00068 }
00069
00070 double polarRadius2() const {
00071 return vector3().polarRadius2();
00072 }
00073
00074
00075 double perp2() const {
00076 return vector3().perp2();
00077 }
00078
00079
00080 double rho2() const {
00081 return vector3().rho2();
00082 }
00083
00084 double polarRadius() const {
00085 return vector3().polarRadius();
00086 }
00087
00088
00089 double perp() const {
00090 return vector3().perp();
00091 }
00092
00093
00094 double rho() const {
00095 return vector3().rho();
00096 }
00097
00098
00099 double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const {
00100 return vector3().azimuthalAngle(mapping);
00101 }
00102
00103
00104 double phi(const PhiMapping mapping = ZERO_2PI) const {
00105 return vector3().phi(mapping);
00106 }
00107
00108
00109 double polarAngle() const {
00110 return vector3().polarAngle();
00111 }
00112
00113
00114 double theta() const {
00115 return vector3().theta();
00116 }
00117
00118 double pseudorapidity() const {
00119 return vector3().pseudorapidity();
00120 }
00121
00122
00123 double eta() const {
00124 return vector3().eta();
00125 }
00126
00127
00128 Vector3 vector3() const {
00129 return Vector3(get(1), get(2), get(3));
00130 }
00131
00132 public:
00133
00134 double contract(const FourVector& v) const {
00135 const double result = t()*v.t() - x()*v.x() - y()*v.y() - z()*v.z();
00136 return result;
00137 }
00138
00139
00140 double dot(const FourVector& v) const {
00141 return contract(v);
00142 }
00143
00144
00145 double operator*(const FourVector& v) const {
00146 return contract(v);
00147 }
00148
00149
00150 FourVector& operator*=(double a) {
00151 _vec = multiply(a, *this)._vec;
00152 return *this;
00153 }
00154
00155
00156 FourVector& operator/=(double a) {
00157 _vec = multiply(1.0/a, *this)._vec;
00158 return *this;
00159 }
00160
00161 FourVector& operator+=(const FourVector& v) {
00162 _vec = add(*this, v)._vec;
00163 return *this;
00164 }
00165
00166 FourVector& operator-=(const FourVector& v) {
00167 _vec = add(*this, -v)._vec;
00168 return *this;
00169 }
00170
00171 FourVector operator-() const {
00172 FourVector result;
00173 result._vec = -_vec;
00174 return result;
00175 }
00176
00177 };
00178
00179
00180
00181 inline double contract(const FourVector& a, const FourVector& b) {
00182 return a.contract(b);
00183 }
00184
00185
00186 inline double dot(const FourVector& a, const FourVector& b) {
00187 return contract(a, b);
00188 }
00189
00190 inline FourVector multiply(const double a, const FourVector& v) {
00191 FourVector result;
00192 result._vec = a * v._vec;
00193 return result;
00194 }
00195
00196 inline FourVector multiply(const FourVector& v, const double a) {
00197 return multiply(a, v);
00198 }
00199
00200 inline FourVector operator*(const double a, const FourVector& v) {
00201 return multiply(a, v);
00202 }
00203
00204 inline FourVector operator*(const FourVector& v, const double a) {
00205 return multiply(a, v);
00206 }
00207
00208 inline FourVector operator/(const FourVector& v, const double a) {
00209 return multiply(1.0/a, v);
00210 }
00211
00212 inline FourVector add(const FourVector& a, const FourVector& b) {
00213 FourVector result;
00214 result._vec = a._vec + b._vec;
00215 return result;
00216 }
00217
00218 inline FourVector operator+(const FourVector& a, const FourVector& b) {
00219 return add(a, b);
00220 }
00221
00222 inline FourVector operator-(const FourVector& a, const FourVector& b) {
00223 return add(a, -b);
00224 }
00225
00226
00227
00228 inline double invariant(const FourVector& lv) {
00229 return lv.invariant();
00230 }
00231
00232
00233 inline double angle(const FourVector& a, const FourVector& b) {
00234 return a.angle(b);
00235 }
00236
00237
00238 inline double angle(const Vector3& a, const FourVector& b) {
00239 return angle( a, b.vector3() );
00240 }
00241
00242
00243 inline double angle(const FourVector& a, const Vector3& b) {
00244 return a.angle(b);
00245 }
00246
00247
00248 inline double polarRadius2(const FourVector& v) {
00249 return v.polarRadius2();
00250 }
00251
00252 inline double perp2(const FourVector& v) {
00253 return v.perp2();
00254 }
00255
00256 inline double rho2(const FourVector& v) {
00257 return v.rho2();
00258 }
00259
00260
00261 inline double polarRadius(const FourVector& v) {
00262 return v.polarRadius();
00263 }
00264
00265 inline double perp(const FourVector& v) {
00266 return v.perp();
00267 }
00268
00269 inline double rho(const FourVector& v) {
00270 return v.rho();
00271 }
00272
00273
00274 inline double azimuthalAngle(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
00275 return v.azimuthalAngle(mapping);
00276 }
00277
00278 inline double phi(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
00279 return v.phi(mapping);
00280 }
00281
00282
00283
00284 inline double polarAngle(const FourVector& v) {
00285 return v.polarAngle();
00286 }
00287
00288 inline double theta(const FourVector& v) {
00289 return v.theta();
00290 }
00291
00292
00293 inline double pseudorapidity(const FourVector& v) {
00294 return v.pseudorapidity();
00295 }
00296
00297 inline double eta(const FourVector& v) {
00298 return v.eta();
00299 }
00300
00301
00302
00303
00304
00305
00306
00307
00308 class FourMomentum : public FourVector {
00309 public:
00310 FourMomentum() { }
00311
00312 template<typename V4>
00313 FourMomentum(const V4& other) {
00314 this->setE(other.t());
00315 this->setPx(other.x());
00316 this->setPy(other.y());
00317 this->setPz(other.z());
00318 }
00319
00320 FourMomentum(const Vector<4>& other)
00321 : FourVector(other) { }
00322
00323 FourMomentum(const double E, const double px, const double py, const double pz) {
00324 this->setE(E);
00325 this->setPx(px);
00326 this->setPy(py);
00327 this->setPz(pz);
00328 }
00329
00330 ~FourMomentum() {}
00331
00332 public:
00333
00334 double E() const { return t(); }
00335
00336
00337 Vector3 p() const { return vector3(); }
00338
00339
00340 double px() const { return x(); }
00341
00342
00343 double py() const { return y(); }
00344
00345
00346 double pz() const { return z(); }
00347
00348
00349 FourMomentum& setE(double E) { setT(E); return *this; }
00350
00351
00352 FourMomentum& setPx(double px) { setX(px); return *this; }
00353
00354
00355 FourMomentum& setPy(double py) { setY(py); return *this; }
00356
00357
00358 FourMomentum& setPz(double pz) { setZ(pz); return *this; }
00359
00360
00361 double mass2() const {
00362 return invariant();
00363 }
00364
00365
00366 double mass() const {
00367 assert(Rivet::isZero(mass2()) || mass2() > 0);
00368 return sqrt(mass2());
00369 }
00370
00371
00372 double rapidity() const {
00373 return 0.5 * std::log( (E() + pz()) / (E() - pz()) );
00374 }
00375
00376
00377 double pT2() const {
00378 return vector3().polarRadius2();
00379 }
00380
00381
00382 double pT() const {
00383 return sqrt(pT2());
00384 }
00385
00386
00387 double Et2() const {
00388 return Et() * Et();
00389 }
00390
00391
00392 double Et() const {
00393 return E() * sin(polarAngle());
00394 }
00395
00396
00397 Vector3 boostVector() const {
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 return Vector3(px()/E(), py()/E(), pz()/E());
00410 }
00411
00412
00413
00414 struct byEAscending{
00415 bool operator()(const FourMomentum &left, const FourMomentum &right) const{
00416 double pt2left = left.E();
00417 double pt2right = right.E();
00418 return pt2left < pt2right;
00419 }
00420
00421 bool operator()(const FourMomentum *left, const FourMomentum *right) const{
00422 return (*this)(left, right);
00423 }
00424 };
00425
00426
00427
00428 struct byEDescending{
00429 bool operator()(const FourMomentum &left, const FourMomentum &right) const{
00430 return byEAscending()(right, left);
00431 }
00432
00433 bool operator()(const FourMomentum *left, const FourVector *right) const{
00434 return (*this)(left, right);
00435 }
00436 };
00437
00438 };
00439
00440
00441
00442 inline double mass2(const FourMomentum& v) {
00443 return v.mass2();
00444 }
00445
00446
00447 inline double mass(const FourMomentum& v) {
00448 return v.mass();
00449 }
00450
00451
00452 inline double rapidity(const FourMomentum& v) {
00453 return v.rapidity();
00454 }
00455
00456
00457 inline double pT2(const FourMomentum& v) {
00458 return v.pT2();
00459 }
00460
00461
00462 inline double pT(const FourMomentum& v) {
00463 return v.pT();
00464 }
00465
00466
00467 inline double Et2(const FourMomentum& v) {
00468 return v.Et2();}
00469
00470
00471 inline double Et(const FourMomentum& v) {
00472 return v.Et();
00473 }
00474
00475
00476 inline Vector3 boostVector(const FourMomentum& v) {
00477 return v.boostVector();
00478 }
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491 inline double deltaR(const FourVector& a, const FourVector& b,
00492 DeltaRScheme scheme = PSEUDORAPIDITY) {
00493 switch (scheme) {
00494 case PSEUDORAPIDITY :
00495 return deltaR(a.vector3(), b.vector3());
00496 case RAPIDITY:
00497 {
00498 const FourMomentum* ma = dynamic_cast<const FourMomentum*>(&a);
00499 const FourMomentum* mb = dynamic_cast<const FourMomentum*>(&b);
00500 if (!ma || !mb) {
00501 string err = "deltaR with scheme RAPIDITY, can be called with FourMomenta only";
00502 throw std::runtime_error(err);
00503 }
00504 return deltaR(*ma, *mb, scheme);
00505 }
00506 default:
00507 throw std::runtime_error("The specified deltaR scheme is not yet implemented");
00508 }
00509 }
00510
00511
00512 inline double deltaR(const FourVector& v,
00513 double eta2, double phi2,
00514 DeltaRScheme scheme = PSEUDORAPIDITY) {
00515 switch (scheme) {
00516 case PSEUDORAPIDITY :
00517 return deltaR(v.vector3(), eta2, phi2);
00518 case RAPIDITY:
00519 {
00520 const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
00521 if (!mv) {
00522 string err = "deltaR with scheme RAPIDITY, can be called with FourMomenta only";
00523 throw std::runtime_error(err);
00524 }
00525 return deltaR(*mv, eta2, phi2, scheme);
00526 }
00527 default:
00528 throw std::runtime_error("The specified deltaR scheme is not yet implemented");
00529 }
00530 }
00531
00532
00533 inline double deltaR(double eta1, double phi1,
00534 const FourVector& v,
00535 DeltaRScheme scheme = PSEUDORAPIDITY) {
00536 switch (scheme) {
00537 case PSEUDORAPIDITY :
00538 return deltaR(eta1, phi1, v.vector3());
00539 case RAPIDITY:
00540 {
00541 const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
00542 if (!mv) {
00543 string err = "deltaR with scheme RAPIDITY, can be called with FourMomenta only";
00544 throw std::runtime_error(err);
00545 }
00546 return deltaR(eta1, phi1, *mv, scheme);
00547 }
00548 default:
00549 throw std::runtime_error("The specified deltaR scheme is not yet implemented");
00550 }
00551 }
00552
00553
00554
00555
00556
00557
00558
00559 inline double deltaR(const FourMomentum& a, const FourMomentum& b,
00560 DeltaRScheme scheme = PSEUDORAPIDITY) {
00561 switch (scheme) {
00562 case PSEUDORAPIDITY:
00563 return deltaR(a.vector3(), b.vector3());
00564 case RAPIDITY:
00565 return deltaR(a.rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
00566 default:
00567 throw std::runtime_error("The specified deltaR scheme is not yet implemented");
00568 }
00569 }
00570
00571 inline double deltaR(const FourMomentum& v,
00572 double eta2, double phi2,
00573 DeltaRScheme scheme = PSEUDORAPIDITY) {
00574 switch (scheme) {
00575 case PSEUDORAPIDITY:
00576 return deltaR(v.vector3(), eta2, phi2);
00577 case RAPIDITY:
00578 return deltaR(v.rapidity(), v.azimuthalAngle(), eta2, phi2);
00579 default:
00580 throw std::runtime_error("The specified deltaR scheme is not yet implemented");
00581 }
00582 }
00583
00584
00585 inline double deltaR(double eta1, double phi1,
00586 const FourMomentum& v,
00587 DeltaRScheme scheme = PSEUDORAPIDITY) {
00588 switch (scheme) {
00589 case PSEUDORAPIDITY:
00590 return deltaR(eta1, phi1, v.vector3());
00591 case RAPIDITY:
00592 return deltaR(eta1, phi1, v.rapidity(), v.azimuthalAngle());
00593 default:
00594 throw std::runtime_error("The specified deltaR scheme is not yet implemented");
00595 }
00596 }
00597
00598
00599
00600
00601
00602
00603 inline const string toString(const FourVector& lv) {
00604 ostringstream out;
00605 out << "(" << (fabs(lv.t()) < 1E-30 ? 0.0 : lv.t())
00606 << "; " << (fabs(lv.x()) < 1E-30 ? 0.0 : lv.x())
00607 << ", " << (fabs(lv.y()) < 1E-30 ? 0.0 : lv.y())
00608 << ", " << (fabs(lv.z()) < 1E-30 ? 0.0 : lv.z())
00609 << ")";
00610 return out.str();
00611 }
00612
00613
00614 inline std::ostream& operator<<(std::ostream& out, const FourVector& lv) {
00615 out << toString(lv);
00616 return out;
00617 }
00618
00619
00620 }
00621
00622 #endif