Rivet  1.8.0
Vector4.hh
1 #ifndef RIVET_MATH_VECTOR4
2 #define RIVET_MATH_VECTOR4
3 
4 #include "Rivet/Math/MathHeader.hh"
5 #include "Rivet/Math/MathUtils.hh"
6 #include "Rivet/Math/VectorN.hh"
7 #include "Rivet/Math/Vector3.hh"
8 
9 namespace Rivet {
10 
11 
12  class FourVector;
13  class FourMomentum;
14  class LorentzTransform;
15  typedef FourVector Vector4;
16  FourVector transform(const LorentzTransform& lt, const FourVector& v4);
17 
18 
20  class FourVector : public Vector<4> {
21  friend FourVector multiply(const double a, const FourVector& v);
22  friend FourVector multiply(const FourVector& v, const double a);
23  friend FourVector add(const FourVector& a, const FourVector& b);
24  friend FourVector transform(const LorentzTransform& lt, const FourVector& v4);
25 
26  public:
27 
28  FourVector() : Vector<4>() { }
29 
30  template<typename V4>
31  FourVector(const V4& other) {
32  this->setT(other.t());
33  this->setX(other.x());
34  this->setY(other.y());
35  this->setZ(other.z());
36  }
37 
38  FourVector(const Vector<4>& other)
39  : Vector<4>(other) { }
40 
41  FourVector(const double t, const double x, const double y, const double z) {
42  this->setT(t);
43  this->setX(x);
44  this->setY(y);
45  this->setZ(z);
46  }
47 
48  virtual ~FourVector() { }
49 
50  public:
51 
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; }
60 
61  double invariant() const {
62  // Done this way for numerical precision
63  return (t() + z())*(t() - z()) - x()*x() - y()*y();
64  }
65 
67  double angle(const FourVector& v) const {
68  return vector3().angle( v.vector3() );
69  }
70 
72  double angle(const Vector3& v3) const {
73  return vector3().angle(v3);
74  }
75 
79  double polarRadius2() const {
80  return vector3().polarRadius2();
81  }
82 
84  double perp2() const {
85  return vector3().perp2();
86  }
87 
89  double rho2() const {
90  return vector3().rho2();
91  }
92 
94  double polarRadius() const {
95  return vector3().polarRadius();
96  }
97 
99  double perp() const {
100  return vector3().perp();
101  }
102 
104  double rho() const {
105  return vector3().rho();
106  }
107 
109  double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const {
110  return vector3().azimuthalAngle(mapping);
111  }
112 
114  double phi(const PhiMapping mapping = ZERO_2PI) const {
115  return vector3().phi(mapping);
116  }
117 
119  double polarAngle() const {
120  return vector3().polarAngle();
121  }
122 
124  double theta() const {
125  return vector3().theta();
126  }
127 
129  double pseudorapidity() const {
130  return vector3().pseudorapidity();
131  }
132 
134  double eta() const {
135  return vector3().eta();
136  }
137 
139  Vector3 vector3() const {
140  return Vector3(get(1), get(2), get(3));
141  }
142 
143 
144  public:
145 
147  double contract(const FourVector& v) const {
148  const double result = t()*v.t() - x()*v.x() - y()*v.y() - z()*v.z();
149  return result;
150  }
151 
153  double dot(const FourVector& v) const {
154  return contract(v);
155  }
156 
158  double operator*(const FourVector& v) const {
159  return contract(v);
160  }
161 
163  FourVector& operator*=(double a) {
164  _vec = multiply(a, *this)._vec;
165  return *this;
166  }
167 
169  FourVector& operator/=(double a) {
170  _vec = multiply(1.0/a, *this)._vec;
171  return *this;
172  }
173 
174  FourVector& operator+=(const FourVector& v) {
175  _vec = add(*this, v)._vec;
176  return *this;
177  }
178 
179  FourVector& operator-=(const FourVector& v) {
180  _vec = add(*this, -v)._vec;
181  return *this;
182  }
183 
185  FourVector result;
186  result._vec = -_vec;
187  return result;
188  }
189 
190  };
191 
192 
194  inline double contract(const FourVector& a, const FourVector& b) {
195  return a.contract(b);
196  }
197 
199  inline double dot(const FourVector& a, const FourVector& b) {
200  return contract(a, b);
201  }
202 
203  inline FourVector multiply(const double a, const FourVector& v) {
204  FourVector result;
205  result._vec = a * v._vec;
206  return result;
207  }
208 
209  inline FourVector multiply(const FourVector& v, const double a) {
210  return multiply(a, v);
211  }
212 
213  inline FourVector operator*(const double a, const FourVector& v) {
214  return multiply(a, v);
215  }
216 
217  inline FourVector operator*(const FourVector& v, const double a) {
218  return multiply(a, v);
219  }
220 
221  inline FourVector operator/(const FourVector& v, const double a) {
222  return multiply(1.0/a, v);
223  }
224 
225  inline FourVector add(const FourVector& a, const FourVector& b) {
226  FourVector result;
227  result._vec = a._vec + b._vec;
228  return result;
229  }
230 
231  inline FourVector operator+(const FourVector& a, const FourVector& b) {
232  return add(a, b);
233  }
234 
235  inline FourVector operator-(const FourVector& a, const FourVector& b) {
236  return add(a, -b);
237  }
238 
241  inline double invariant(const FourVector& lv) {
242  return lv.invariant();
243  }
244 
246  inline double angle(const FourVector& a, const FourVector& b) {
247  return a.angle(b);
248  }
249 
251  inline double angle(const Vector3& a, const FourVector& b) {
252  return angle( a, b.vector3() );
253  }
254 
256  inline double angle(const FourVector& a, const Vector3& b) {
257  return a.angle(b);
258  }
259 
261  inline double polarRadius2(const FourVector& v) {
262  return v.polarRadius2();
263  }
265  inline double perp2(const FourVector& v) {
266  return v.perp2();
267  }
269  inline double rho2(const FourVector& v) {
270  return v.rho2();
271  }
272 
274  inline double polarRadius(const FourVector& v) {
275  return v.polarRadius();
276  }
278  inline double perp(const FourVector& v) {
279  return v.perp();
280  }
282  inline double rho(const FourVector& v) {
283  return v.rho();
284  }
285 
287  inline double azimuthalAngle(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
288  return v.azimuthalAngle(mapping);
289  }
291  inline double phi(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
292  return v.phi(mapping);
293  }
294 
295 
297  inline double polarAngle(const FourVector& v) {
298  return v.polarAngle();
299  }
301  inline double theta(const FourVector& v) {
302  return v.theta();
303  }
304 
306  inline double pseudorapidity(const FourVector& v) {
307  return v.pseudorapidity();
308  }
310  inline double eta(const FourVector& v) {
311  return v.eta();
312  }
313 
314 
315 
317 
318 
319 
321  class FourMomentum : public FourVector {
322  friend FourMomentum multiply(const double a, const FourMomentum& v);
323  friend FourMomentum multiply(const FourMomentum& v, const double a);
324  friend FourMomentum add(const FourMomentum& a, const FourMomentum& b);
325  friend FourMomentum transform(const LorentzTransform& lt, const FourMomentum& v4);
326 
327  public:
328  FourMomentum() { }
329 
330  template<typename V4>
331  FourMomentum(const V4& other) {
332  this->setE(other.t());
333  this->setPx(other.x());
334  this->setPy(other.y());
335  this->setPz(other.z());
336  }
337 
338  FourMomentum(const Vector<4>& other)
339  : FourVector(other) { }
340 
341  FourMomentum(const double E, const double px, const double py, const double pz) {
342  this->setE(E);
343  this->setPx(px);
344  this->setPy(py);
345  this->setPz(pz);
346  }
347 
348  ~FourMomentum() {}
349 
350  public:
352  double E() const { return t(); }
353 
355  Vector3 p() const { return vector3(); }
356 
358  double px() const { return x(); }
359 
361  double py() const { return y(); }
362 
364  double pz() const { return z(); }
365 
367  FourMomentum& setE(double E) { setT(E); return *this; }
368 
370  FourMomentum& setPx(double px) { setX(px); return *this; }
371 
373  FourMomentum& setPy(double py) { setY(py); return *this; }
374 
376  FourMomentum& setPz(double pz) { setZ(pz); return *this; }
377 
379  double mass() const {
380  assert(Rivet::isZero(mass2()) || mass2() > 0);
381  if (Rivet::isZero(mass2())) {
382  return 0.0;
383  } else {
384  return sqrt(mass2());
385  }
386  }
387 
389  double massT() const {
390  return mass() * sin(polarAngle());
391  }
392 
394  double mass2() const {
395  return invariant();
396  }
397 
399  double massT2() const {
400  return massT() * massT();
401  }
402 
404  double rapidity() const {
405  return 0.5 * std::log( (E() + pz()) / (E() - pz()) );
406  }
407 
409  double pT2() const {
410  return vector3().polarRadius2();
411  }
412 
414  double pT() const {
415  return sqrt(pT2());
416  }
417 
419  double Et2() const {
420  return Et() * Et();
421  }
422 
424  double Et() const {
425  return E() * sin(polarAngle());
426  }
427 
430  // const Vector3 p3 = vector3();
431  // const double m2 = mass2();
432  // if (Rivet::isZero(m2)) return p3.unit();
433  // else {
434  // // Could also do this via beta = tanh(rapidity), but that's
435  // // probably more messy from a numerical hygiene point of view.
436  // const double p2 = p3.mod2();
437  // const double beta = sqrt( p2 / (m2 + p2) );
438  // return beta * p3.unit();
439  // }
441  return Vector3(px()/E(), py()/E(), pz()/E());
442  }
443 
445  struct byEAscending {
446  bool operator()(const FourMomentum& left, const FourMomentum& right) const{
447  double pt2left = left.E();
448  double pt2right = right.E();
449  return pt2left < pt2right;
450  }
451 
452  bool operator()(const FourMomentum* left, const FourMomentum* right) const{
453  return (*this)(left, right);
454  }
455  };
456 
458  struct byEDescending {
459  bool operator()(const FourMomentum& left, const FourMomentum& right) const{
460  return byEAscending()(right, left);
461  }
462 
463  bool operator()(const FourMomentum* left, const FourVector* right) const{
464  return (*this)(left, right);
465  }
466  };
467 
468 
471  _vec = multiply(a, *this)._vec;
472  return *this;
473  }
474 
477  _vec = multiply(1.0/a, *this)._vec;
478  return *this;
479  }
480 
481  FourMomentum& operator+=(const FourMomentum& v) {
482  _vec = add(*this, v)._vec;
483  return *this;
484  }
485 
486  FourMomentum& operator-=(const FourMomentum& v) {
487  _vec = add(*this, -v)._vec;
488  return *this;
489  }
490 
492  FourMomentum result;
493  result._vec = -_vec;
494  return result;
495  }
496 
497 
498  };
499 
500 
501  inline FourMomentum multiply(const double a, const FourMomentum& v) {
502  FourMomentum result;
503  result._vec = a * v._vec;
504  return result;
505  }
506 
507  inline FourMomentum multiply(const FourMomentum& v, const double a) {
508  return multiply(a, v);
509  }
510 
511  inline FourMomentum operator*(const double a, const FourMomentum& v) {
512  return multiply(a, v);
513  }
514 
515  inline FourMomentum operator*(const FourMomentum& v, const double a) {
516  return multiply(a, v);
517  }
518 
519  inline FourMomentum operator/(const FourMomentum& v, const double a) {
520  return multiply(1.0/a, v);
521  }
522 
523  inline FourMomentum add(const FourMomentum& a, const FourMomentum& b) {
524  FourMomentum result;
525  result._vec = a._vec + b._vec;
526  return result;
527  }
528 
529  inline FourMomentum operator+(const FourMomentum& a, const FourMomentum& b) {
530  return add(a, b);
531  }
532 
533  inline FourMomentum operator-(const FourMomentum& a, const FourMomentum& b) {
534  return add(a, -b);
535  }
536 
537 
538 
540  inline double mass(const FourMomentum& v) {
541  return v.mass();
542  }
543 
545  inline double massT(const FourMomentum& v) {
546  return v.massT();
547  }
548 
550  inline double mass2(const FourMomentum& v) {
551  return v.mass2();
552  }
553 
555  inline double massT2(const FourMomentum& v) {
556  return v.massT2();
557  }
558 
560  inline double rapidity(const FourMomentum& v) {
561  return v.rapidity();
562  }
563 
565  inline double pT2(const FourMomentum& v) {
566  return v.pT2();
567  }
568 
570  inline double pT(const FourMomentum& v) {
571  return v.pT();
572  }
573 
575  inline double Et2(const FourMomentum& v) {
576  return v.Et2();}
577 
579  inline double Et(const FourMomentum& v) {
580  return v.Et();
581  }
582 
584  inline Vector3 boostVector(const FourMomentum& v) {
585  return v.boostVector();
586  }
587 
588 
590 
591 
593 
594 
602  inline double deltaR(const FourVector& a, const FourVector& b,
603  RapScheme scheme = PSEUDORAPIDITY) {
604  switch (scheme) {
605  case PSEUDORAPIDITY :
606  return deltaR(a.vector3(), b.vector3());
607  case RAPIDITY:
608  {
609  const FourMomentum* ma = dynamic_cast<const FourMomentum*>(&a);
610  const FourMomentum* mb = dynamic_cast<const FourMomentum*>(&b);
611  if (!ma || !mb) {
612  string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
613  throw std::runtime_error(err);
614  }
615  return deltaR(*ma, *mb, scheme);
616  }
617  default:
618  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
619  }
620  }
621 
622 
628  inline double deltaR(const FourVector& v,
629  double eta2, double phi2,
630  RapScheme scheme = PSEUDORAPIDITY) {
631  switch (scheme) {
632  case PSEUDORAPIDITY :
633  return deltaR(v.vector3(), eta2, phi2);
634  case RAPIDITY:
635  {
636  const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
637  if (!mv) {
638  string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
639  throw std::runtime_error(err);
640  }
641  return deltaR(*mv, eta2, phi2, scheme);
642  }
643  default:
644  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
645  }
646  }
647 
648 
654  inline double deltaR(double eta1, double phi1,
655  const FourVector& v,
656  RapScheme scheme = PSEUDORAPIDITY) {
657  switch (scheme) {
658  case PSEUDORAPIDITY :
659  return deltaR(eta1, phi1, v.vector3());
660  case RAPIDITY:
661  {
662  const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
663  if (!mv) {
664  string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
665  throw std::runtime_error(err);
666  }
667  return deltaR(eta1, phi1, *mv, scheme);
668  }
669  default:
670  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
671  }
672  }
673 
674 
680  inline double deltaR(const FourMomentum& a, const FourMomentum& b,
681  RapScheme scheme = PSEUDORAPIDITY) {
682  switch (scheme) {
683  case PSEUDORAPIDITY:
684  return deltaR(a.vector3(), b.vector3());
685  case RAPIDITY:
686  return deltaR(a.rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
687  default:
688  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
689  }
690  }
691 
697  inline double deltaR(const FourMomentum& v,
698  double eta2, double phi2,
699  RapScheme scheme = PSEUDORAPIDITY) {
700  switch (scheme) {
701  case PSEUDORAPIDITY:
702  return deltaR(v.vector3(), eta2, phi2);
703  case RAPIDITY:
704  return deltaR(v.rapidity(), v.azimuthalAngle(), eta2, phi2);
705  default:
706  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
707  }
708  }
709 
710 
716  inline double deltaR(double eta1, double phi1,
717  const FourMomentum& v,
718  RapScheme scheme = PSEUDORAPIDITY) {
719  switch (scheme) {
720  case PSEUDORAPIDITY:
721  return deltaR(eta1, phi1, v.vector3());
722  case RAPIDITY:
723  return deltaR(eta1, phi1, v.rapidity(), v.azimuthalAngle());
724  default:
725  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
726  }
727  }
728 
734  inline double deltaR(const FourMomentum& a, const FourVector& b,
735  RapScheme scheme = PSEUDORAPIDITY) {
736  switch (scheme) {
737  case PSEUDORAPIDITY:
738  return deltaR(a.vector3(), b.vector3());
739  case RAPIDITY:
740  return deltaR(a.rapidity(), a.azimuthalAngle(), FourMomentum(b).rapidity(), b.azimuthalAngle());
741  default:
742  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
743  }
744  }
745 
751  inline double deltaR(const FourVector& a, const FourMomentum& b,
752  RapScheme scheme = PSEUDORAPIDITY) {
753  switch (scheme) {
754  case PSEUDORAPIDITY:
755  return deltaR(a.vector3(), b.vector3());
756  case RAPIDITY:
757  return deltaR(FourMomentum(a).rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
758  default:
759  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
760  }
761  }
762 
765  inline double deltaR(const FourMomentum& a, const Vector3& b) {
766  return deltaR(a.vector3(), b);
767  }
768 
771  inline double deltaR(const Vector3& a, const FourMomentum& b) {
772  return deltaR(a, b.vector3());
773  }
774 
777  inline double deltaR(const FourVector& a, const Vector3& b) {
778  return deltaR(a.vector3(), b);
779  }
780 
783  inline double deltaR(const Vector3& a, const FourVector& b) {
784  return deltaR(a, b.vector3());
785  }
786 
788 
790 
791 
793  inline double deltaPhi(const FourMomentum& a, const FourMomentum& b) {
794  return deltaPhi(a.vector3(), b.vector3());
795  }
796 
798  inline double deltaPhi(const FourMomentum& v, double phi2) {
799  return deltaPhi(v.vector3(), phi2);
800  }
801 
803  inline double deltaPhi(double phi1, const FourMomentum& v) {
804  return deltaPhi(phi1, v.vector3());
805  }
806 
808  inline double deltaPhi(const FourVector& a, const FourVector& b) {
809  return deltaPhi(a.vector3(), b.vector3());
810  }
811 
813  inline double deltaPhi(const FourVector& v, double phi2) {
814  return deltaPhi(v.vector3(), phi2);
815  }
816 
818  inline double deltaPhi(double phi1, const FourVector& v) {
819  return deltaPhi(phi1, v.vector3());
820  }
821 
823  inline double deltaPhi(const FourVector& a, const FourMomentum& b) {
824  return deltaPhi(a.vector3(), b.vector3());
825  }
826 
828  inline double deltaPhi(const FourMomentum& a, const FourVector& b) {
829  return deltaPhi(a.vector3(), b.vector3());
830  }
831 
833  inline double deltaPhi(const FourVector& a, const Vector3& b) {
834  return deltaPhi(a.vector3(), b);
835  }
836 
838  inline double deltaPhi(const Vector3& a, const FourVector& b) {
839  return deltaPhi(a, b.vector3());
840  }
841 
843  inline double deltaPhi(const FourMomentum& a, const Vector3& b) {
844  return deltaPhi(a.vector3(), b);
845  }
846 
848  inline double deltaPhi(const Vector3& a, const FourMomentum& b) {
849  return deltaPhi(a, b.vector3());
850  }
851 
853 
855 
856 
858  inline double deltaEta(const FourMomentum& a, const FourMomentum& b) {
859  return deltaEta(a.vector3(), b.vector3());
860  }
861 
863  inline double deltaEta(const FourMomentum& v, double eta2) {
864  return deltaEta(v.vector3(), eta2);
865  }
866 
868  inline double deltaEta(double eta1, const FourMomentum& v) {
869  return deltaEta(eta1, v.vector3());
870  }
871 
873  inline double deltaEta(const FourVector& a, const FourVector& b) {
874  return deltaEta(a.vector3(), b.vector3());
875  }
876 
878  inline double deltaEta(const FourVector& v, double eta2) {
879  return deltaEta(v.vector3(), eta2);
880  }
881 
883  inline double deltaEta(double eta1, const FourVector& v) {
884  return deltaEta(eta1, v.vector3());
885  }
886 
888  inline double deltaEta(const FourVector& a, const FourMomentum& b) {
889  return deltaEta(a.vector3(), b.vector3());
890  }
891 
893  inline double deltaEta(const FourMomentum& a, const FourVector& b) {
894  return deltaEta(a.vector3(), b.vector3());
895  }
896 
898  inline double deltaEta(const FourVector& a, const Vector3& b) {
899  return deltaEta(a.vector3(), b);
900  }
901 
903  inline double deltaEta(const Vector3& a, const FourVector& b) {
904  return deltaEta(a, b.vector3());
905  }
906 
908  inline double deltaEta(const FourMomentum& a, const Vector3& b) {
909  return deltaEta(a.vector3(), b);
910  }
911 
913  inline double deltaEta(const Vector3& a, const FourMomentum& b) {
914  return deltaEta(a, b.vector3());
915  }
916 
918 
920 
921 
923 
924 
926  inline const string toString(const FourVector& lv) {
927  ostringstream out;
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())
932  << ")";
933  return out.str();
934  }
935 
937  inline std::ostream& operator<<(std::ostream& out, const FourVector& lv) {
938  out << toString(lv);
939  return out;
940  }
941 
943 
944 
945 }
946 
947 #endif