Rivet  1.8.0
Vector3.hh
1 #ifndef RIVET_MATH_VECTOR3
2 #define RIVET_MATH_VECTOR3
3 
4 #include "Rivet/Math/MathHeader.hh"
5 #include "Rivet/Math/MathUtils.hh"
6 #include "Rivet/Math/VectorN.hh"
7 
8 namespace Rivet {
9 
10 
11  class Vector3;
12  typedef Vector3 ThreeVector;
13  class Matrix3;
14 
15  Vector3 multiply(const double, const Vector3&);
16  Vector3 multiply(const Vector3&, const double);
17  Vector3 add(const Vector3&, const Vector3&);
18  Vector3 operator*(const double, const Vector3&);
19  Vector3 operator*(const Vector3&, const double);
20  Vector3 operator/(const Vector3&, const double);
21  Vector3 operator+(const Vector3&, const Vector3&);
22  Vector3 operator-(const Vector3&, const Vector3&);
23 
24 
26  class Vector3 : public Vector<3> {
27 
28  friend class Matrix3;
29  friend Vector3 multiply(const double, const Vector3&);
30  friend Vector3 multiply(const Vector3&, const double);
31  friend Vector3 add(const Vector3&, const Vector3&);
32  friend Vector3 subtract(const Vector3&, const Vector3&);
33 
34  public:
35  Vector3() : Vector<3>() { }
36 
37  template<typename V3>
38  Vector3(const V3& other) {
39  this->setX(other.x());
40  this->setY(other.y());
41  this->setZ(other.z());
42  }
43 
44  Vector3(const Vector<3>& other) {
45  this->setX(other.get(0));
46  this->setY(other.get(1));
47  this->setZ(other.get(2));
48  }
49 
50  Vector3(double x, double y, double z) {
51  this->setX(x);
52  this->setY(y);
53  this->setZ(z);
54  }
55 
56  ~Vector3() { }
57 
58  public:
59  static Vector3 mkX() { return Vector3(1,0,0); }
60  static Vector3 mkY() { return Vector3(0,1,0); }
61  static Vector3 mkZ() { return Vector3(0,0,1); }
62 
63  public:
64  double x() const { return get(0); }
65  double y() const { return get(1); }
66  double z() const { return get(2); }
67  Vector3& setX(double x) { set(0, x); return *this; }
68  Vector3& setY(double y) { set(1, y); return *this; }
69  Vector3& setZ(double z) { set(2, z); return *this; }
70 
71  double dot(const Vector3& v) const {
72  return _vec.dot(v._vec);
73  }
74 
75  Vector3 cross(const Vector3& v) const {
76  Vector3 result;
77  result._vec = _vec.cross(v._vec);
78  return result;
79  }
80 
81  double angle(const Vector3& v) const {
82  double localDotOther = unit().dot(v.unit());
83  if(Rivet::isZero(localDotOther - 1.0)) return 0.0;
84  return acos( localDotOther );
85  }
86 
87  Vector3 unit() const {
89  if (isZero()) return *this;
90  else return *this * 1.0/this->mod();
91  }
92 
93  double polarRadius2() const {
94  return x()*x() + y()*y();
95  }
96 
98  double perp2() const {
99  return polarRadius2();
100  }
101 
103  double rho2() const {
104  return polarRadius2();
105  }
106 
107  double polarRadius() const {
108  return sqrt(polarRadius2());
109  }
110 
112  double perp() const {
113  return polarRadius();
114  }
115 
117  double rho() const {
118  return polarRadius();
119  }
120 
122  double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const {
123  // If this is a null vector, return zero rather than let atan2 set an error state
124  if (Rivet::isZero(mod2())) return 0.0;
125 
126  // Calculate the arctan and correct for numerical boundary cases
127  double value = atan2( y(), x() );
128  if (value > 2*PI || value < -2*PI){
129  value = fmod(value, 2*PI);
130  }
131  if (value <= -PI) value += 2*PI;
132  if (value > PI) value -= 2*PI;
133 
134  // Return in the requested range
135  switch (mapping) {
136  case MINUSPI_PLUSPI:
137  assert(value > -PI && value <= PI);
138  return value;
139  case ZERO_2PI:
140  if (value >= 0) {
141  assert(value >= 0 && value < 2*PI);
142  return value;
143  } else if (Rivet::isZero(value)) {
144  value = 0.0;
145  return value;
146  } else {
147  value = 2*PI + value;
148  assert(value >= 0 && value < 2*PI);
149  return value;
150  }
151  default:
152  throw std::runtime_error("The specified phi mapping scheme is not yet implemented");
153  }
154  }
155 
157  double phi(const PhiMapping mapping = ZERO_2PI) const {
158  return azimuthalAngle(mapping);
159  }
160 
162  double polarAngle() const {
163  // Get number beween [0,PI]
164  double polarangle = atan2(polarRadius(), z());
165  assert(polarangle >= -PI && polarangle <= PI);
166  return polarangle;
167  }
168 
170  double theta() const {
171  return polarAngle();
172  }
173 
176  double pseudorapidity() const {
177  return -std::log(tan( 0.5 * polarAngle() ));
178  }
179 
181  double eta() const {
182  return pseudorapidity();
183  }
184 
185  public:
186  Vector3& operator*=(const double a) {
187  _vec = multiply(a, *this)._vec;
188  return *this;
189  }
190 
191  Vector3& operator/=(const double a) {
192  _vec = multiply(1.0/a, *this)._vec;
193  return *this;
194  }
195 
196  Vector3& operator+=(const Vector3& v) {
197  _vec = add(*this, v)._vec;
198  return *this;
199  }
200 
201  Vector3& operator-=(const Vector3& v) {
202  _vec = subtract(*this, v)._vec;
203  return *this;
204  }
205 
206  Vector3 operator-() const {
207  Vector3 rtn;
208  rtn._vec = -_vec;
209  return rtn;
210  }
211 
212  };
213 
214 
215 
216  inline double dot(const Vector3& a, const Vector3& b) {
217  return a.dot(b);
218  }
219 
220  inline Vector3 cross(const Vector3& a, const Vector3& b) {
221  return a.cross(b);
222  }
223 
224  inline Vector3 multiply(const double a, const Vector3& v) {
225  Vector3 result;
226  result._vec = a * v._vec;
227  return result;
228  }
229 
230  inline Vector3 multiply(const Vector3& v, const double a) {
231  return multiply(a, v);
232  }
233 
234  inline Vector3 operator*(const double a, const Vector3& v) {
235  return multiply(a, v);
236  }
237 
238  inline Vector3 operator*(const Vector3& v, const double a) {
239  return multiply(a, v);
240  }
241 
242  inline Vector3 operator/(const Vector3& v, const double a) {
243  return multiply(1.0/a, v);
244  }
245 
246  inline Vector3 add(const Vector3& a, const Vector3& b) {
247  Vector3 result;
248  result._vec = a._vec + b._vec;
249  return result;
250  }
251 
252  inline Vector3 subtract(const Vector3& a, const Vector3& b) {
253  Vector3 result;
254  result._vec = a._vec - b._vec;
255  return result;
256  }
257 
258  inline Vector3 operator+(const Vector3& a, const Vector3& b) {
259  return add(a, b);
260  }
261 
262  inline Vector3 operator-(const Vector3& a, const Vector3& b) {
263  return subtract(a, b);
264  }
265 
266  // More physicsy coordinates etc.
267 
269  inline double angle(const Vector3& a, const Vector3& b) {
270  return a.angle(b);
271  }
272 
274  inline double polarRadius2(const Vector3& v) {
275  return v.polarRadius2();
276  }
278  inline double perp2(const Vector3& v) {
279  return v.perp2();
280  }
282  inline double rho2(const Vector3& v) {
283  return v.rho2();
284  }
285 
287  inline double polarRadius(const Vector3& v) {
288  return v.polarRadius();
289  }
291  inline double perp(const Vector3& v) {
292  return v.perp();
293  }
295  inline double rho(const Vector3& v) {
296  return v.rho();
297  }
298 
299 
302  inline double azimuthalAngle(const Vector3& v, const PhiMapping mapping = ZERO_2PI) {
303  return v.azimuthalAngle(mapping);
304  }
306  inline double phi(const Vector3& v, const PhiMapping mapping = ZERO_2PI) {
307  return v.phi(mapping);
308  }
309 
311  inline double polarAngle(const Vector3& v) {
312  return v.polarAngle();
313  }
315  inline double theta(const Vector3& v) {
316  return v.theta();
317  }
318 
320  inline double pseudorapidity(const Vector3& v) {
321  return v.pseudorapidity();
322  }
324  inline double eta(const Vector3& v) {
325  return v.eta();
326  }
327 
328 
330 
332 
333 
335  inline double deltaEta(const Vector3& a, const Vector3& b) {
336  return deltaEta(a.pseudorapidity(), b.pseudorapidity());
337  }
338 
340  inline double deltaEta(const Vector3& v, double eta2) {
341  return deltaEta(v.pseudorapidity(), eta2);
342  }
343 
345  inline double deltaEta(double eta1, const Vector3& v) {
346  return deltaEta(eta1, v.pseudorapidity());
347  }
348 
350 
351 
353 
354 
356  inline double deltaPhi(const Vector3& a, const Vector3& b) {
357  return deltaPhi(a.azimuthalAngle(), b.azimuthalAngle());
358  }
359 
361  inline double deltaPhi(const Vector3& v, double phi2) {
362  return deltaPhi(v.azimuthalAngle(), phi2);
363  }
364 
366  inline double deltaPhi(double phi1, const Vector3& v) {
367  return deltaPhi(phi1, v.azimuthalAngle());
368  }
369 
371 
372 
374 
375 
377  inline double deltaR(const Vector3& a, const Vector3& b) {
378  return deltaR(a.pseudorapidity(), a.azimuthalAngle(),
379  b.pseudorapidity(), b.azimuthalAngle());
380  }
381 
383  inline double deltaR(const Vector3& v, double eta2, double phi2) {
384  return deltaR(v.pseudorapidity(), v.azimuthalAngle(), eta2, phi2);
385  }
386 
388  inline double deltaR(double eta1, double phi1, const Vector3& v) {
389  return deltaR(eta1, phi1, v.pseudorapidity(), v.azimuthalAngle());
390  }
391 
393 
394 }
395 
396 #endif