thepeg is hosted by Hepforge, IPPP Durham
ThePEG  2.1.5
ThreeVector.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // ThreeVector.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad
5 //
6 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef ThePEG_ThreeVector_H
10 #define ThePEG_ThreeVector_H
11 
19 #include "ThreeVector.fh"
20 #include "ThePEG/Config/ThePEG.h"
21 #include "ThePEG/Utilities/UnitIO.h"
22 #include <cassert>
23 #include <cmath>
24 
25 namespace ThePEG {
26 
33 template <typename Value>
35 {
36 private:
41 
42 public:
45  ThreeVector()
46  : theX(), theY(), theZ() {}
47 
48  ThreeVector(Value x, Value y, Value z)
49  : theX(x), theY(y), theZ(z) {}
50 
51  template<typename ValueB>
53  : theX(v.x()), theY(v.y()), theZ(v.z()) {}
55 
56 public:
58 
59  Value x() const { return theX; }
60  Value y() const { return theY; }
61  Value z() const { return theZ; }
63 
65 
66  void setX(Value x) { theX = x; }
67  void setY(Value y) { theY = y; }
68  void setZ(Value z) { theZ = z; }
70 
71 public:
73  Value2 mag2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
74 
76  Value mag() const { return sqrt(mag2()); }
77 
79  Value2 perp2() const { return sqr(x()) + sqr(y()); }
80 
82  Value perp() const { return sqrt(perp2()); }
83 
85  template <typename U>
87  dot(const ThreeVector<U> & a) const {
88  return x()*a.x() + y()*a.y() + z()*a.z();
89  }
90 
92  template <typename U>
93  Value2 perp2(const ThreeVector<U> & p) const {
94  typedef typename BinaryOpTraits<U,U>::MulT pSqType;
95  const pSqType pMag2 = p.mag2();
96  assert( pMag2 > pSqType() );
98  ss = this->dot(p);
99  Value2 ret = mag2() - sqr(ss)/pMag2;
100  if ( ret <= Value2() )
101  ret = Value2();
102  return ret;
103  }
104 
106  template <typename U>
107  Value perp(const ThreeVector<U> & p) const {
108  return sqrt(perp2(p));
109  }
110 
112 
113  double theta() const {
115  assert(!(x() == Value() && y() == Value() && z() == Value()));
116  return atan2(perp(),z());
117  }
118 
120  double phi() const {
121  return atan2(y(),x());
122  }
123 
125  void setTheta(double th) {
126  double ma = mag();
127  double ph = phi();
128  setX(ma*sin(th)*cos(ph));
129  setY(ma*sin(th)*sin(ph));
130  setZ(ma*cos(th));
131  }
132 
134  void setPhi(double ph) {
135  double xy = perp();
136  setX(xy*cos(ph));
137  setY(xy*sin(ph));
138  }
140 
143  Value2 mg2 = mag2();
144  assert(mg2 > Value2());
145  Value mg = sqrt(mg2);
146  return ThreeVector<double>(x()/mg, y()/mg, z()/mg);
147  }
148 
151  Value xx = abs(x());
152  Value yy = abs(y());
153  Value zz = abs(z());
154  if (xx < yy) {
155  return xx < zz ? ThreeVector<Value>(Value(),z(),-y())
156  : ThreeVector<Value>(y(),-x(),Value());
157  } else {
158  return yy < zz ? ThreeVector<Value>(-z(),Value(),x())
159  : ThreeVector<Value>(y(),-x(),Value());
160  }
161  }
162 
164  template <typename U>
165  double deltaPhi (const ThreeVector<U> & v2) const {
166  double dphi = v2.phi() - phi();
167  if ( dphi > Constants::pi ) {
168  dphi -= Constants::twopi;
169  } else if ( dphi <= -Constants::pi ) {
170  dphi += Constants::twopi;
171  }
172  return dphi;
173  }
174 
180  template <typename U>
181  ThreeVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
182  if (angle == 0.0)
183  return *this;
184  const U ll = axis.mag();
185  assert( ll > U() );
186 
187  const double sa = sin(angle), ca = cos(angle);
188  const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
189  const Value xx = x(), yy = y(), zz = z();
190 
191  setX((ca+(1-ca)*dx*dx) * xx
192  +((1-ca)*dx*dy-sa*dz) * yy
193  +((1-ca)*dx*dz+sa*dy) * zz
194  );
195  setY(((1-ca)*dy*dx+sa*dz) * xx
196  +(ca+(1-ca)*dy*dy) * yy
197  +((1-ca)*dy*dz-sa*dx) * zz
198  );
199  setZ(((1-ca)*dz*dx-sa*dy) * xx
200  +((1-ca)*dz*dy+sa*dx) * yy
201  +(ca+(1-ca)*dz*dz) * zz
202  );
203  return *this;
204  }
205 
206 
210  ThreeVector<Value> & rotateUz (const Axis & axis) {
211  Axis ax = axis.unit();
212  double u1 = ax.x();
213  double u2 = ax.y();
214  double u3 = ax.z();
215  double up = u1*u1 + u2*u2;
216  if (up>0) {
217  up = sqrt(up);
218  Value px = x(), py = y(), pz = z();
219  setX( (u1*u3*px - u2*py)/up + u1*pz );
220  setY( (u2*u3*px + u1*py)/up + u2*pz );
221  setZ( -up*px + u3*pz );
222  }
223  else if (u3 < 0.) {
224  setX(-x());
225  setZ(-z());
226  }
227  return *this;
228  }
229 
234  Axis ax = axis.unit();
235  double u1 = ax.x();
236  double u2 = ax.y();
237  double u3 = ax.z();
238  double up = u1*u1 + u2*u2;
239  if (up>0) {
240  up = sqrt(up);
241  Value px = x(), py = y(), pz = z();
242  setX( ( u1*u3*px + u2*u3*py)/up - up*pz );
243  setY( (-u2*px + u1*py)/up );
244  setZ( u1*px + u2*py + u3*pz );
245  }
246  else if (u3 < 0.) {
247  setX(-x());
248  setZ(-z());
249  }
250  return *this;
251  }
252 
254  template <typename U>
256  cross(const ThreeVector<U> & a) const {
258  return ResultT( y()*a.z()-z()*a.y(),
259  -x()*a.z()+z()*a.x(),
260  x()*a.y()-y()*a.x());
261  }
262 
263  public:
265 
266  bool operator==(const ThreeVector<Value> & a) const {
267  return (theX == a.x() && theY == a.y() && theZ == a.z());
268  }
269  bool operator!=(const ThreeVector<Value> & a) const {
270  return !(*this == a);
271  }
272  bool almostEqual(const ThreeVector<Value> & a, double threshold = 1e-04) const {
273  return ((std::abs(theX - a.x()) < threshold) && (std::abs(theY - a.y()) < threshold) && (std::abs(theZ - a.z()) < threshold));
274  }
275  bool almostUnequal(const ThreeVector<Value> & a, double threshold = 1e-04) const {
276  return ! this->almostEqual(a, threshold);
277  }
279 
280 public:
282 
283  ThreeVector<Value> & operator+=(const ThreeVector<Value> & a) {
284  theX += a.x();
285  theY += a.y();
286  theZ += a.z();
287  return *this;
288  }
289 
290  ThreeVector<Value> & operator-=(const ThreeVector<Value> & a) {
291  theX -= a.x();
292  theY -= a.y();
293  theZ -= a.z();
294  return *this;
295  }
296 
297  ThreeVector<Value> & operator*=(double a) {
298  theX *= a;
299  theY *= a;
300  theZ *= a;
301  return *this;
302  }
303 
304  ThreeVector<Value> & operator/=(double a) {
305  theX /= a;
306  theY /= a;
307  theZ /= a;
308  return *this;
309  }
311 
313  template <typename U>
314  double cosTheta(const ThreeVector<U> & q) const {
315  typedef typename BinaryOpTraits<Value,U>::MulT
316  ProdType;
317  ProdType ptot = mag()*q.mag();
318  assert( ptot > ProdType() );
319  double arg = dot(q)/ptot;
320  if (arg > 1.0) arg = 1.0;
321  else if(arg < -1.0) arg = -1.0;
322  return arg;
323  }
324 
326  template <typename U>
327  double angle(const ThreeVector<U> & v) const {
328  return acos(cosTheta(v));
329  }
330 
331 private:
333 
334  Value theX;
335  Value theY;
336  Value theZ;
338 };
339 
341 inline ostream &
342 operator<< (ostream & os, const ThreeVector<double> & v)
343 {
344  return os << '(' << v.x() << ',' << v.y() << ',' << v.z() << ')';
345 }
346 
348 
349 template <typename Value>
350 inline ThreeVector<Value>
351 operator+(ThreeVector<Value> a,
352  const ThreeVector<Value> & b)
353 {
354  return a += b;
355 }
356 
357 template <typename Value>
358 inline ThreeVector<Value>
359 operator-(ThreeVector<Value> a,
360  const ThreeVector<Value> & b)
361 {
362  return a -= b;
363 }
364 
365 template <typename Value>
366 inline ThreeVector<Value> operator-(const ThreeVector<Value> & v) {
367  return ThreeVector<Value>(-v.x(),-v.y(),-v.z());
368 }
369 
370 template <typename Value>
371 inline ThreeVector<Value> operator*(ThreeVector<Value> v, double a) {
372  return v *= a;
373 }
374 
375 template <typename Value>
376 inline ThreeVector<Value> operator*(double a, ThreeVector<Value> v) {
377  return v *= a;
378 }
379 
380 template <typename ValueA, typename ValueB>
382 operator*(ValueB a, ThreeVector<ValueA> v) {
383  typedef typename BinaryOpTraits<ValueA,ValueB>::MulT ResultT;
384  return ThreeVector<ResultT>(a*v.x(), a*v.y(), a*v.z());
385 }
386 
387 template <typename ValueA, typename ValueB>
389 operator*(ThreeVector<ValueA> v, ValueB a) {
390  return a*v;
391 }
393 
395 template <typename ValueA, typename ValueB>
397 operator*(const ThreeVector<ValueA> & a,
398  const ThreeVector<ValueB> & b)
399 {
400  return a.dot(b);
401 }
402 
404 template <typename Value>
406  return v.unit();
407 }
408 
409 
411 template <typename OStream, typename UT, typename Value>
412 void ounitstream(OStream & os, const ThreeVector<Value> & p, UT & u) {
413  os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u);
414 }
415 
417 template <typename IStream, typename UT, typename Value>
418 void iunitstream(IStream & is, ThreeVector<Value> & p, UT & u) {
419  Value x, y, z;
420  is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u);
421  p = ThreeVector<Value>(x, y, z);
422 }
423 
424 }
425 
426 #endif /* ThePEG_ThreeVector_H */
ThreeVector< typename BinaryOpTraits< Value, U >::MulT > cross(const ThreeVector< U > &a) const
Vector cross-product.
Definition: ThreeVector.h:256
BinaryOpTraits< Value, Value >::MulT Value2
Value squared.
Definition: ThreeVector.h:38
ThreeVector< Value > & rotateUzBack(const Axis &axis)
Rotate from a reference frame to the z-axis.
Definition: ThreeVector.h:233
Value mag() const
Magnitude .
Definition: ThreeVector.h:76
ThreeVector< Value > & rotate(double angle, const ThreeVector< U > &axis)
Apply a rotation.
Definition: ThreeVector.h:181
Value2 perp2() const
Squared transverse component .
Definition: ThreeVector.h:79
BinaryOpTraits< Value, U >::MulT dot(const ThreeVector< U > &a) const
Dot product.
Definition: ThreeVector.h:87
BinaryOpTraits< Value2, Value2 >::MulT Value4
Value to the 4th power.
Definition: ThreeVector.h:40
void ounitstream(OStream &os, const vector< T, Alloc > &v, UT &u)
Ouput a vector of objects with the specified unit.
Definition: Containers.h:275
double deltaPhi(const ThreeVector< U > &v2) const
Azimuthal angle difference, brought into the range .
Definition: ThreeVector.h:165
A 3-component vector.
Definition: ThreeVector.h:34
Value2 mag2() const
Squared magnitude .
Definition: ThreeVector.h:73
constexpr double pi
Good old .
Definition: Constants.h:54
Value perp() const
Transverse component .
Definition: ThreeVector.h:82
This is the main namespace within which all identifiers in ThePEG are declared.
Definition: FactoryBase.h:28
This is the main config header file for ThePEG.
void setPhi(double ph)
Set the azimuthal angle.
Definition: ThreeVector.h:134
void iunitstream(IStream &is, vector< T, Alloc > &v, UT &u)
Input a vector of objects with the specified unit.
Definition: Containers.h:289
double theta() const
Polar angle.
Definition: ThreeVector.h:114
double angle(const ThreeVector< U > &v) const
Angle between two vectors.
Definition: ThreeVector.h:327
constexpr double twopi
Good old .
Definition: Constants.h:57
OUnit< T, UT > ounit(const T &t, const UT &ut)
Helper function creating a OUnit object given an object and a unit.
Definition: UnitIO.h:84
double cosTheta(const ThreeVector< U > &q) const
Cosine of the azimuthal angle between two vectors.
Definition: ThreeVector.h:314
ThreeVector< double > unit() const
Parallel vector with unit length.
Definition: ThreeVector.h:142
ThreeVector< Value > orthogonal() const
Orthogonal vector.
Definition: ThreeVector.h:150
ThreeVector< Value > & rotateUz(const Axis &axis)
Rotate the reference frame to a new z-axis.
Definition: ThreeVector.h:210
Value perp(const ThreeVector< U > &p) const
Transverse component with respect to the given axis.
Definition: ThreeVector.h:107
double phi() const
Azimuthal angle.
Definition: ThreeVector.h:120
Value2 perp2(const ThreeVector< U > &p) const
Squared transverse component with respect to the given axis.
Definition: ThreeVector.h:93
BinaryOpTraits should be specialized with typdefs called MulT and DivT which gives the type resulting...
Definition: PhysicalQty.h:244
void setTheta(double th)
Set the polar angle.
Definition: ThreeVector.h:125
ThreeVector< double > unitVector(const ThreeVector< Value > &v)
A parallel vector with unit length.
Definition: ThreeVector.h:405
IUnit< T, UT > iunit(T &t, const UT &ut)
Helper function creating a IUnit object given an object and a unit.
Definition: UnitIO.h:91