Rivet  1.8.3
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 
176  _vec = add(*this, v)._vec;
177  return *this;
178  }
179 
182  _vec = add(*this, -v)._vec;
183  return *this;
184  }
185 
188  FourVector result;
189  result._vec = -_vec;
190  return result;
191  }
192 
193  };
194 
195 
197  inline double contract(const FourVector& a, const FourVector& b) {
198  return a.contract(b);
199  }
200 
202  inline double dot(const FourVector& a, const FourVector& b) {
203  return contract(a, b);
204  }
205 
206  inline FourVector multiply(const double a, const FourVector& v) {
207  FourVector result;
208  result._vec = a * v._vec;
209  return result;
210  }
211 
212  inline FourVector multiply(const FourVector& v, const double a) {
213  return multiply(a, v);
214  }
215 
216  inline FourVector operator*(const double a, const FourVector& v) {
217  return multiply(a, v);
218  }
219 
220  inline FourVector operator*(const FourVector& v, const double a) {
221  return multiply(a, v);
222  }
223 
224  inline FourVector operator/(const FourVector& v, const double a) {
225  return multiply(1.0/a, v);
226  }
227 
228  inline FourVector add(const FourVector& a, const FourVector& b) {
229  FourVector result;
230  result._vec = a._vec + b._vec;
231  return result;
232  }
233 
234  inline FourVector operator+(const FourVector& a, const FourVector& b) {
235  return add(a, b);
236  }
237 
238  inline FourVector operator-(const FourVector& a, const FourVector& b) {
239  return add(a, -b);
240  }
241 
244  inline double invariant(const FourVector& lv) {
245  return lv.invariant();
246  }
247 
249  inline double angle(const FourVector& a, const FourVector& b) {
250  return a.angle(b);
251  }
252 
254  inline double angle(const Vector3& a, const FourVector& b) {
255  return angle( a, b.vector3() );
256  }
257 
259  inline double angle(const FourVector& a, const Vector3& b) {
260  return a.angle(b);
261  }
262 
264  inline double polarRadius2(const FourVector& v) {
265  return v.polarRadius2();
266  }
268  inline double perp2(const FourVector& v) {
269  return v.perp2();
270  }
272  inline double rho2(const FourVector& v) {
273  return v.rho2();
274  }
275 
277  inline double polarRadius(const FourVector& v) {
278  return v.polarRadius();
279  }
281  inline double perp(const FourVector& v) {
282  return v.perp();
283  }
285  inline double rho(const FourVector& v) {
286  return v.rho();
287  }
288 
290  inline double azimuthalAngle(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
291  return v.azimuthalAngle(mapping);
292  }
294  inline double phi(const FourVector& v, const PhiMapping mapping = ZERO_2PI) {
295  return v.phi(mapping);
296  }
297 
298 
300  inline double polarAngle(const FourVector& v) {
301  return v.polarAngle();
302  }
304  inline double theta(const FourVector& v) {
305  return v.theta();
306  }
307 
309  inline double pseudorapidity(const FourVector& v) {
310  return v.pseudorapidity();
311  }
313  inline double eta(const FourVector& v) {
314  return v.eta();
315  }
316 
317 
318 
320 
321 
322 
324  class FourMomentum : public FourVector {
325  friend FourMomentum multiply(const double a, const FourMomentum& v);
326  friend FourMomentum multiply(const FourMomentum& v, const double a);
327  friend FourMomentum add(const FourMomentum& a, const FourMomentum& b);
328  friend FourMomentum transform(const LorentzTransform& lt, const FourMomentum& v4);
329 
330  public:
331  FourMomentum() { }
332 
333  template<typename V4>
334  FourMomentum(const V4& other) {
335  this->setE(other.t());
336  this->setPx(other.x());
337  this->setPy(other.y());
338  this->setPz(other.z());
339  }
340 
341  FourMomentum(const Vector<4>& other)
342  : FourVector(other) { }
343 
344  FourMomentum(const double E, const double px, const double py, const double pz) {
345  this->setE(E);
346  this->setPx(px);
347  this->setPy(py);
348  this->setPz(pz);
349  }
350 
351  ~FourMomentum() {}
352 
353  public:
355  double E() const { return t(); }
356 
358  Vector3 p() const { return vector3(); }
359 
361  double px() const { return x(); }
362 
364  double py() const { return y(); }
365 
367  double pz() const { return z(); }
368 
370  FourMomentum& setE(double E) { setT(E); return *this; }
371 
373  FourMomentum& setPx(double px) { setX(px); return *this; }
374 
376  FourMomentum& setPy(double py) { setY(py); return *this; }
377 
379  FourMomentum& setPz(double pz) { setZ(pz); return *this; }
380 
384  double mass() const {
385  // assert(Rivet::isZero(mass2()) || mass2() > 0);
386  // if (Rivet::isZero(mass2())) {
387  // return 0.0;
388  // } else {
389  // return sqrt(mass2());
390  // }
391  return sign(mass2()) * sqrt(fabs(mass2()));
392  }
393 
395  double mass2() const {
396  return invariant();
397  }
398 
400  double rapidity() const {
401  return 0.5 * std::log( (E() + pz()) / (E() - pz()) );
402  }
403 
405  double pT2() const {
406  return vector3().polarRadius2();
407  }
408 
410  double pT() const {
411  return sqrt(pT2());
412  }
413 
415  double Et2() const {
416  return Et() * Et();
417  }
418 
420  double Et() const {
421  return E() * sin(polarAngle());
422  }
423 
426  // const Vector3 p3 = vector3();
427  // const double m2 = mass2();
428  // if (Rivet::isZero(m2)) return p3.unit();
429  // else {
430  // // Could also do this via beta = tanh(rapidity), but that's
431  // // probably more messy from a numerical hygiene point of view.
432  // const double p2 = p3.mod2();
433  // const double beta = sqrt( p2 / (m2 + p2) );
434  // return beta * p3.unit();
435  // }
437  return Vector3(px()/E(), py()/E(), pz()/E());
438  }
439 
441  struct byEAscending {
442  bool operator()(const FourMomentum& left, const FourMomentum& right) const{
443  double pt2left = left.E();
444  double pt2right = right.E();
445  return pt2left < pt2right;
446  }
447 
448  bool operator()(const FourMomentum* left, const FourMomentum* right) const{
449  return (*this)(left, right);
450  }
451  };
452 
454  struct byEDescending {
455  bool operator()(const FourMomentum& left, const FourMomentum& right) const{
456  return byEAscending()(right, left);
457  }
458 
459  bool operator()(const FourMomentum* left, const FourVector* right) const{
460  return (*this)(left, right);
461  }
462  };
463 
464 
467  _vec = multiply(a, *this)._vec;
468  return *this;
469  }
470 
473  _vec = multiply(1.0/a, *this)._vec;
474  return *this;
475  }
476 
479  _vec = add(*this, v)._vec;
480  return *this;
481  }
482 
485  _vec = add(*this, -v)._vec;
486  return *this;
487  }
488 
491  FourMomentum result;
492  result._vec = -_vec;
493  return result;
494  }
495 
496 
497  };
498 
499 
500  inline FourMomentum multiply(const double a, const FourMomentum& v) {
501  FourMomentum result;
502  result._vec = a * v._vec;
503  return result;
504  }
505 
506  inline FourMomentum multiply(const FourMomentum& v, const double a) {
507  return multiply(a, v);
508  }
509 
510  inline FourMomentum operator*(const double a, const FourMomentum& v) {
511  return multiply(a, v);
512  }
513 
514  inline FourMomentum operator*(const FourMomentum& v, const double a) {
515  return multiply(a, v);
516  }
517 
518  inline FourMomentum operator/(const FourMomentum& v, const double a) {
519  return multiply(1.0/a, v);
520  }
521 
522  inline FourMomentum add(const FourMomentum& a, const FourMomentum& b) {
523  FourMomentum result;
524  result._vec = a._vec + b._vec;
525  return result;
526  }
527 
528  inline FourMomentum operator+(const FourMomentum& a, const FourMomentum& b) {
529  return add(a, b);
530  }
531 
532  inline FourMomentum operator-(const FourMomentum& a, const FourMomentum& b) {
533  return add(a, -b);
534  }
535 
536 
537 
539  inline double mass(const FourMomentum& v) {
540  return v.mass();
541  }
542 
544  inline double mass2(const FourMomentum& v) {
545  return v.mass2();
546  }
547 
549  inline double rapidity(const FourMomentum& v) {
550  return v.rapidity();
551  }
552 
554  inline double pT2(const FourMomentum& v) {
555  return v.pT2();
556  }
557 
559  inline double pT(const FourMomentum& v) {
560  return v.pT();
561  }
562 
564  inline double Et2(const FourMomentum& v) {
565  return v.Et2();}
566 
568  inline double Et(const FourMomentum& v) {
569  return v.Et();
570  }
571 
573  inline Vector3 boostVector(const FourMomentum& v) {
574  return v.boostVector();
575  }
576 
577 
579 
580 
582 
583 
591  inline double deltaR(const FourVector& a, const FourVector& b,
592  RapScheme scheme = PSEUDORAPIDITY) {
593  switch (scheme) {
594  case PSEUDORAPIDITY :
595  return deltaR(a.vector3(), b.vector3());
596  case RAPIDITY:
597  {
598  const FourMomentum* ma = dynamic_cast<const FourMomentum*>(&a);
599  const FourMomentum* mb = dynamic_cast<const FourMomentum*>(&b);
600  if (!ma || !mb) {
601  string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
602  throw std::runtime_error(err);
603  }
604  return deltaR(*ma, *mb, scheme);
605  }
606  default:
607  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
608  }
609  }
610 
611 
617  inline double deltaR(const FourVector& v,
618  double eta2, double phi2,
619  RapScheme scheme = PSEUDORAPIDITY) {
620  switch (scheme) {
621  case PSEUDORAPIDITY :
622  return deltaR(v.vector3(), eta2, phi2);
623  case RAPIDITY:
624  {
625  const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
626  if (!mv) {
627  string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
628  throw std::runtime_error(err);
629  }
630  return deltaR(*mv, eta2, phi2, scheme);
631  }
632  default:
633  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
634  }
635  }
636 
637 
643  inline double deltaR(double eta1, double phi1,
644  const FourVector& v,
645  RapScheme scheme = PSEUDORAPIDITY) {
646  switch (scheme) {
647  case PSEUDORAPIDITY :
648  return deltaR(eta1, phi1, v.vector3());
649  case RAPIDITY:
650  {
651  const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
652  if (!mv) {
653  string err = "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
654  throw std::runtime_error(err);
655  }
656  return deltaR(eta1, phi1, *mv, scheme);
657  }
658  default:
659  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
660  }
661  }
662 
663 
669  inline double deltaR(const FourMomentum& a, const FourMomentum& b,
670  RapScheme scheme = PSEUDORAPIDITY) {
671  switch (scheme) {
672  case PSEUDORAPIDITY:
673  return deltaR(a.vector3(), b.vector3());
674  case RAPIDITY:
675  return deltaR(a.rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
676  default:
677  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
678  }
679  }
680 
686  inline double deltaR(const FourMomentum& v,
687  double eta2, double phi2,
688  RapScheme scheme = PSEUDORAPIDITY) {
689  switch (scheme) {
690  case PSEUDORAPIDITY:
691  return deltaR(v.vector3(), eta2, phi2);
692  case RAPIDITY:
693  return deltaR(v.rapidity(), v.azimuthalAngle(), eta2, phi2);
694  default:
695  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
696  }
697  }
698 
699 
705  inline double deltaR(double eta1, double phi1,
706  const FourMomentum& v,
707  RapScheme scheme = PSEUDORAPIDITY) {
708  switch (scheme) {
709  case PSEUDORAPIDITY:
710  return deltaR(eta1, phi1, v.vector3());
711  case RAPIDITY:
712  return deltaR(eta1, phi1, v.rapidity(), v.azimuthalAngle());
713  default:
714  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
715  }
716  }
717 
723  inline double deltaR(const FourMomentum& a, const FourVector& b,
724  RapScheme scheme = PSEUDORAPIDITY) {
725  switch (scheme) {
726  case PSEUDORAPIDITY:
727  return deltaR(a.vector3(), b.vector3());
728  case RAPIDITY:
729  return deltaR(a.rapidity(), a.azimuthalAngle(), FourMomentum(b).rapidity(), b.azimuthalAngle());
730  default:
731  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
732  }
733  }
734 
740  inline double deltaR(const FourVector& a, const FourMomentum& b,
741  RapScheme scheme = PSEUDORAPIDITY) {
742  switch (scheme) {
743  case PSEUDORAPIDITY:
744  return deltaR(a.vector3(), b.vector3());
745  case RAPIDITY:
746  return deltaR(FourMomentum(a).rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
747  default:
748  throw std::runtime_error("The specified deltaR scheme is not yet implemented");
749  }
750  }
751 
754  inline double deltaR(const FourMomentum& a, const Vector3& b) {
755  return deltaR(a.vector3(), b);
756  }
757 
760  inline double deltaR(const Vector3& a, const FourMomentum& b) {
761  return deltaR(a, b.vector3());
762  }
763 
766  inline double deltaR(const FourVector& a, const Vector3& b) {
767  return deltaR(a.vector3(), b);
768  }
769 
772  inline double deltaR(const Vector3& a, const FourVector& b) {
773  return deltaR(a, b.vector3());
774  }
775 
777 
779 
780 
782  inline double deltaPhi(const FourMomentum& a, const FourMomentum& b) {
783  return deltaPhi(a.vector3(), b.vector3());
784  }
785 
787  inline double deltaPhi(const FourMomentum& v, double phi2) {
788  return deltaPhi(v.vector3(), phi2);
789  }
790 
792  inline double deltaPhi(double phi1, const FourMomentum& v) {
793  return deltaPhi(phi1, v.vector3());
794  }
795 
797  inline double deltaPhi(const FourVector& a, const FourVector& b) {
798  return deltaPhi(a.vector3(), b.vector3());
799  }
800 
802  inline double deltaPhi(const FourVector& v, double phi2) {
803  return deltaPhi(v.vector3(), phi2);
804  }
805 
807  inline double deltaPhi(double phi1, const FourVector& v) {
808  return deltaPhi(phi1, v.vector3());
809  }
810 
812  inline double deltaPhi(const FourVector& a, const FourMomentum& b) {
813  return deltaPhi(a.vector3(), b.vector3());
814  }
815 
817  inline double deltaPhi(const FourMomentum& a, const FourVector& b) {
818  return deltaPhi(a.vector3(), b.vector3());
819  }
820 
822  inline double deltaPhi(const FourVector& a, const Vector3& b) {
823  return deltaPhi(a.vector3(), b);
824  }
825 
827  inline double deltaPhi(const Vector3& a, const FourVector& b) {
828  return deltaPhi(a, b.vector3());
829  }
830 
832  inline double deltaPhi(const FourMomentum& a, const Vector3& b) {
833  return deltaPhi(a.vector3(), b);
834  }
835 
837  inline double deltaPhi(const Vector3& a, const FourMomentum& b) {
838  return deltaPhi(a, b.vector3());
839  }
840 
842 
844 
845 
847  inline double deltaEta(const FourMomentum& a, const FourMomentum& b) {
848  return deltaEta(a.vector3(), b.vector3());
849  }
850 
852  inline double deltaEta(const FourMomentum& v, double eta2) {
853  return deltaEta(v.vector3(), eta2);
854  }
855 
857  inline double deltaEta(double eta1, const FourMomentum& v) {
858  return deltaEta(eta1, v.vector3());
859  }
860 
862  inline double deltaEta(const FourVector& a, const FourVector& b) {
863  return deltaEta(a.vector3(), b.vector3());
864  }
865 
867  inline double deltaEta(const FourVector& v, double eta2) {
868  return deltaEta(v.vector3(), eta2);
869  }
870 
872  inline double deltaEta(double eta1, const FourVector& v) {
873  return deltaEta(eta1, v.vector3());
874  }
875 
877  inline double deltaEta(const FourVector& a, const FourMomentum& b) {
878  return deltaEta(a.vector3(), b.vector3());
879  }
880 
882  inline double deltaEta(const FourMomentum& a, const FourVector& b) {
883  return deltaEta(a.vector3(), b.vector3());
884  }
885 
887  inline double deltaEta(const FourVector& a, const Vector3& b) {
888  return deltaEta(a.vector3(), b);
889  }
890 
892  inline double deltaEta(const Vector3& a, const FourVector& b) {
893  return deltaEta(a, b.vector3());
894  }
895 
897  inline double deltaEta(const FourMomentum& a, const Vector3& b) {
898  return deltaEta(a.vector3(), b);
899  }
900 
902  inline double deltaEta(const Vector3& a, const FourMomentum& b) {
903  return deltaEta(a, b.vector3());
904  }
905 
907 
909 
910 
912 
913 
915  inline std::string toString(const FourVector& lv) {
916  ostringstream out;
917  out << "(" << (fabs(lv.t()) < 1E-30 ? 0.0 : lv.t())
918  << "; " << (fabs(lv.x()) < 1E-30 ? 0.0 : lv.x())
919  << ", " << (fabs(lv.y()) < 1E-30 ? 0.0 : lv.y())
920  << ", " << (fabs(lv.z()) < 1E-30 ? 0.0 : lv.z())
921  << ")";
922  return out.str();
923  }
924 
926  inline std::ostream& operator<<(std::ostream& out, const FourVector& lv) {
927  out << toString(lv);
928  return out;
929  }
930 
932 
933 
934 }
935 
936 #endif