1 #ifndef RIVET_MATH_VECTOR4
2 #define RIVET_MATH_VECTOR4
4 #include "Rivet/Math/MathHeader.hh"
5 #include "Rivet/Math/MathUtils.hh"
6 #include "Rivet/Math/VectorN.hh"
7 #include "Rivet/Math/Vector3.hh"
14 class LorentzTransform;
15 typedef FourVector Vector4;
16 FourVector transform(
const LorentzTransform& lt,
const FourVector& v4);
32 this->setT(other.t());
33 this->setX(other.x());
34 this->setY(other.y());
35 this->setZ(other.z());
41 FourVector(
const double t,
const double x,
const double y,
const double z) {
52 double t()
const {
return get(0); }
53 double x()
const {
return get(1); }
54 double y()
const {
return get(2); }
55 double z()
const {
return get(3); }
56 FourVector& setT(
const double t) {
set(0, t);
return *
this; }
57 FourVector& setX(
const double x) {
set(1, x);
return *
this; }
58 FourVector& setY(
const double y) {
set(2, y);
return *
this; }
59 FourVector& setZ(
const double z) {
set(3, z);
return *
this; }
61 double invariant()
const {
63 return (t() + z())*(t() - z()) - x()*x() - y()*y();
80 return vector3().polarRadius2();
140 return Vector3(
get(1),
get(2),
get(3));
148 const double result = t()*v.t() - x()*v.x() - y()*v.y() - z()*v.z();
164 _vec = multiply(a, *
this)._vec;
170 _vec = multiply(1.0/a, *
this)._vec;
175 _vec = add(*
this, v)._vec;
179 FourVector& operator-=(
const FourVector& v) {
180 _vec = add(*
this, -v)._vec;
203 inline FourVector multiply(
const double a,
const FourVector& v) {
205 result._vec = a * v._vec;
209 inline FourVector multiply(
const FourVector& v,
const double a) {
210 return multiply(a, v);
213 inline FourVector operator*(
const double a,
const FourVector& v) {
214 return multiply(a, v);
217 inline FourVector operator*(
const FourVector& v,
const double a) {
218 return multiply(a, v);
221 inline FourVector operator/(
const FourVector& v,
const double a) {
222 return multiply(1.0/a, v);
225 inline FourVector add(
const FourVector& a,
const FourVector& b) {
227 result._vec = a._vec + b._vec;
231 inline FourVector operator+(
const FourVector& a,
const FourVector& b) {
235 inline FourVector operator-(
const FourVector& a,
const FourVector& b) {
242 return lv.invariant();
292 return v.
phi(mapping);
330 template<
typename V4>
332 this->
setE(other.t());
333 this->
setPx(other.x());
334 this->
setPy(other.y());
335 this->
setPz(other.z());
352 double E()
const {
return t(); }
358 double px()
const {
return x(); }
361 double py()
const {
return y(); }
364 double pz()
const {
return z(); }
384 return sqrt(
mass2());
405 return 0.5 * std::log( (
E() +
pz()) / (
E() -
pz()) );
410 return vector3().polarRadius2();
447 double pt2left = left.
E();
448 double pt2right = right.
E();
449 return pt2left < pt2right;
453 return (*
this)(left, right);
464 return (*
this)(left, right);
471 _vec = multiply(a, *
this)._vec;
477 _vec = multiply(1.0/a, *
this)._vec;
482 _vec = add(*
this, v)._vec;
486 FourMomentum& operator-=(
const FourMomentum& v) {
487 _vec = add(*
this, -v)._vec;
501 inline FourMomentum multiply(
const double a,
const FourMomentum& v) {
503 result._vec = a * v._vec;
507 inline FourMomentum multiply(
const FourMomentum& v,
const double a) {
508 return multiply(a, v);
511 inline FourMomentum operator*(
const double a,
const FourMomentum& v) {
512 return multiply(a, v);
515 inline FourMomentum operator*(
const FourMomentum& v,
const double a) {
516 return multiply(a, v);
519 inline FourMomentum operator/(
const FourMomentum& v,
const double a) {
520 return multiply(1.0/a, v);
523 inline FourMomentum add(
const FourMomentum& a,
const FourMomentum& b) {
525 result._vec = a._vec + b._vec;
529 inline FourMomentum operator+(
const FourMomentum& a,
const FourMomentum& b) {
533 inline FourMomentum operator-(
const FourMomentum& a,
const FourMomentum& b) {
605 case PSEUDORAPIDITY :
612 string err =
"deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
613 throw std::runtime_error(err);
615 return deltaR(*ma, *mb, scheme);
618 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
629 double eta2,
double phi2,
632 case PSEUDORAPIDITY :
633 return deltaR(v.
vector3(), eta2, phi2);
638 string err =
"deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
639 throw std::runtime_error(err);
641 return deltaR(*mv, eta2, phi2, scheme);
644 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
654 inline double deltaR(
double eta1,
double phi1,
658 case PSEUDORAPIDITY :
659 return deltaR(eta1, phi1, v.
vector3());
664 string err =
"deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
665 throw std::runtime_error(err);
667 return deltaR(eta1, phi1, *mv, scheme);
670 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
688 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
698 double eta2,
double phi2,
702 return deltaR(v.
vector3(), eta2, phi2);
706 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
716 inline double deltaR(
double eta1,
double phi1,
721 return deltaR(eta1, phi1, v.
vector3());
725 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
742 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
759 throw std::runtime_error(
"The specified deltaR scheme is not yet implemented");
799 return deltaPhi(v.
vector3(), phi2);
804 return deltaPhi(phi1, v.
vector3());
814 return deltaPhi(v.
vector3(), phi2);
819 return deltaPhi(phi1, v.
vector3());
834 return deltaPhi(a.
vector3(), b);
839 return deltaPhi(a, b.
vector3());
844 return deltaPhi(a.
vector3(), b);
849 return deltaPhi(a, b.
vector3());
864 return deltaEta(v.
vector3(), eta2);
869 return deltaEta(eta1, v.
vector3());
879 return deltaEta(v.
vector3(), eta2);
884 return deltaEta(eta1, v.
vector3());
899 return deltaEta(a.
vector3(), b);
904 return deltaEta(a, b.
vector3());
909 return deltaEta(a.
vector3(), b);
914 return deltaEta(a, b.
vector3());
928 out <<
"(" << (fabs(lv.t()) < 1E-30 ? 0.0 : lv.t())
929 <<
"; " << (fabs(lv.x()) < 1E-30 ? 0.0 : lv.x())
930 <<
", " << (fabs(lv.y()) < 1E-30 ? 0.0 : lv.y())
931 <<
", " << (fabs(lv.z()) < 1E-30 ? 0.0 : lv.z())