ThePEG  1.8.0
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-2011 David Grellscheid, Leif Lonnblad
5 //
6 // ThePEG is licenced under version 2 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 
24 namespace ThePEG {
25 
32 template <typename Value>
34 {
35 private:
40 
41 public:
44  ThreeVector()
45  : theX(), theY(), theZ() {}
46 
47  ThreeVector(Value x, Value y, Value z)
48  : theX(x), theY(y), theZ(z) {}
49 
50  template<typename ValueB>
51  ThreeVector(const ThreeVector<ValueB> & v)
52  : theX(v.x()), theY(v.y()), theZ(v.z()) {}
54 
55 public:
57 
58  Value x() const { return theX; }
59  Value y() const { return theY; }
60  Value z() const { return theZ; }
62 
64 
65  void setX(Value x) { theX = x; }
66  void setY(Value y) { theY = y; }
67  void setZ(Value z) { theZ = z; }
69 
70 public:
72  Value2 mag2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
73 
75  Value mag() const { return sqrt(mag2()); }
76 
78  Value2 perp2() const { return sqr(x()) + sqr(y()); }
79 
81  Value perp() const { return sqrt(perp2()); }
82 
84  template <typename U>
86  dot(const ThreeVector<U> & a) const {
87  return x()*a.x() + y()*a.y() + z()*a.z();
88  }
89 
91  template <typename U>
92  Value2 perp2(const ThreeVector<U> & p) const {
93  typedef typename BinaryOpTraits<U,U>::MulT pSqType;
94  const pSqType pMag2 = p.mag2();
95  assert( pMag2 > pSqType() );
97  ss = this->dot(p);
98  Value2 ret = mag2() - sqr(ss)/pMag2;
99  if ( ret <= Value2() )
100  ret = Value2();
101  return ret;
102  }
103 
105  template <typename U>
106  Value perp(const ThreeVector<U> & p) const {
107  return sqrt(perp2(p));
108  }
109 
111 
112  double theta() const {
114  assert(!(x() == Value() && y() == Value() && z() == Value()));
115  return atan2(perp(),z());
116  }
117 
119  double phi() const {
120  return atan2(y(),x());
121  }
122 
124  void setTheta(double th) {
125  double ma = mag();
126  double ph = phi();
127  setX(ma*sin(th)*cos(ph));
128  setY(ma*sin(th)*sin(ph));
129  setZ(ma*cos(th));
130  }
131 
133  void setPhi(double ph) {
134  double xy = perp();
135  setX(xy*cos(ph));
136  setY(xy*sin(ph));
137  }
139 
142  Value2 mg2 = mag2();
143  assert(mg2 > Value2());
144  Value mg = sqrt(mg2);
145  return ThreeVector<double>(x()/mg, y()/mg, z()/mg);
146  }
147 
150  Value xx = abs(x());
151  Value yy = abs(y());
152  Value zz = abs(z());
153  if (xx < yy) {
154  return xx < zz ? ThreeVector<Value>(Value(),z(),-y())
155  : ThreeVector<Value>(y(),-x(),Value());
156  } else {
157  return yy < zz ? ThreeVector<Value>(-z(),Value(),x())
158  : ThreeVector<Value>(y(),-x(),Value());
159  }
160  }
161 
163  template <typename U>
164  double deltaPhi (const ThreeVector<U> & v2) const {
165  double dphi = v2.phi() - phi();
166  if ( dphi > Constants::pi ) {
167  dphi -= Constants::twopi;
168  } else if ( dphi <= -Constants::pi ) {
169  dphi += Constants::twopi;
170  }
171  return dphi;
172  }
173 
179  template <typename U>
180  ThreeVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
181  if (angle == 0.0)
182  return *this;
183  const U ll = axis.mag();
184  assert( ll > U() );
185 
186  const double sa = sin(angle), ca = cos(angle);
187  const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
188  const Value xx = x(), yy = y(), zz = z();
189 
190  setX((ca+(1-ca)*dx*dx) * xx
191  +((1-ca)*dx*dy-sa*dz) * yy
192  +((1-ca)*dx*dz+sa*dy) * zz
193  );
194  setY(((1-ca)*dy*dx+sa*dz) * xx
195  +(ca+(1-ca)*dy*dy) * yy
196  +((1-ca)*dy*dz-sa*dx) * zz
197  );
198  setZ(((1-ca)*dz*dx-sa*dy) * xx
199  +((1-ca)*dz*dy+sa*dx) * yy
200  +(ca+(1-ca)*dz*dz) * zz
201  );
202  return *this;
203  }
204 
205 
209  ThreeVector<Value> & rotateUz (const Axis & axis) {
210  Axis ax = axis.unit();
211  double u1 = ax.x();
212  double u2 = ax.y();
213  double u3 = ax.z();
214  double up = u1*u1 + u2*u2;
215  if (up>0) {
216  up = sqrt(up);
217  Value px = x(), py = y(), pz = z();
218  setX( (u1*u3*px - u2*py)/up + u1*pz );
219  setY( (u2*u3*px + u1*py)/up + u2*pz );
220  setZ( -up*px + u3*pz );
221  }
222  else if (u3 < 0.) {
223  setX(-x());
224  setZ(-z());
225  }
226  return *this;
227  }
228 
230  template <typename U>
232  cross(const ThreeVector<U> & a) const {
234  return ResultT( y()*a.z()-z()*a.y(),
235  -x()*a.z()+z()*a.x(),
236  x()*a.y()-y()*a.x());
237  }
238 
239 public:
241 
242  ThreeVector<Value> & operator+=(const ThreeVector<Value> & a) {
243  theX += a.x();
244  theY += a.y();
245  theZ += a.z();
246  return *this;
247  }
248 
249  ThreeVector<Value> & operator-=(const ThreeVector<Value> & a) {
250  theX -= a.x();
251  theY -= a.y();
252  theZ -= a.z();
253  return *this;
254  }
255 
256  ThreeVector<Value> & operator*=(double a) {
257  theX *= a;
258  theY *= a;
259  theZ *= a;
260  return *this;
261  }
262 
263  ThreeVector<Value> & operator/=(double a) {
264  theX /= a;
265  theY /= a;
266  theZ /= a;
267  return *this;
268  }
270 
272  template <typename U>
273  double cosTheta(const ThreeVector<U> & q) const {
274  typedef typename BinaryOpTraits<Value,U>::MulT
275  ProdType;
276  ProdType ptot = mag()*q.mag();
277  assert( ptot > ProdType() );
278  double arg = dot(q)/ptot;
279  if (arg > 1.0) arg = 1.0;
280  else if(arg < -1.0) arg = -1.0;
281  return arg;
282  }
283 
285  template <typename U>
286  double angle(const ThreeVector<U> & v) const {
287  return acos(cosTheta(v));
288  }
289 
290 private:
292 
293  Value theX;
294  Value theY;
295  Value theZ;
297 };
298 
300 inline ostream &
301 operator<< (ostream & os, const ThreeVector<double> & v)
302 {
303  return os << '(' << v.x() << ',' << v.y() << ',' << v.z() << ')';
304 }
305 
307 
308 template <typename Value>
309 inline ThreeVector<Value>
310 operator+(ThreeVector<Value> a,
311  const ThreeVector<Value> & b)
312 {
313  return a += b;
314 }
315 
316 template <typename Value>
317 inline ThreeVector<Value>
318 operator-(ThreeVector<Value> a,
319  const ThreeVector<Value> & b)
320 {
321  return a -= b;
322 }
323 
324 template <typename Value>
325 inline ThreeVector<Value> operator-(const ThreeVector<Value> & v) {
326  return ThreeVector<Value>(-v.x(),-v.y(),-v.z());
327 }
328 
329 template <typename Value>
330 inline ThreeVector<Value> operator*(ThreeVector<Value> v, double a) {
331  return v *= a;
332 }
333 
334 template <typename Value>
335 inline ThreeVector<Value> operator*(double a, ThreeVector<Value> v) {
336  return v *= a;
337 }
338 
339 template <typename ValueA, typename ValueB>
340 inline ThreeVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
341 operator*(ValueB a, ThreeVector<ValueA> v) {
342  typedef typename BinaryOpTraits<ValueA,ValueB>::MulT ResultT;
343  return ThreeVector<ResultT>(a*v.x(), a*v.y(), a*v.z());
344 }
345 
346 template <typename ValueA, typename ValueB>
347 inline ThreeVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
348 operator*(ThreeVector<ValueA> v, ValueB a) {
349  return a*v;
350 }
352 
354 template <typename ValueA, typename ValueB>
355 inline typename BinaryOpTraits<ValueA,ValueB>::MulT
356 operator*(const ThreeVector<ValueA> & a,
357  const ThreeVector<ValueB> & b)
358 {
359  return a.dot(b);
360 }
361 
363 template <typename Value>
365  return v.unit();
366 }
367 
368 
370 template <typename OStream, typename UT, typename Value>
371 void ounitstream(OStream & os, const ThreeVector<Value> & p, UT & u) {
372  os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u);
373 }
374 
376 template <typename IStream, typename UT, typename Value>
377 void iunitstream(IStream & is, ThreeVector<Value> & p, UT & u) {
378  Value x, y, z;
379  is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u);
380  p = ThreeVector<Value>(x, y, z);
381 }
382 
383 }
384 
385 #endif /* ThePEG_ThreeVector_H */
Value2 perp2(const ThreeVector< U > &p) const
Squared transverse component with respect to the given axis.
Definition: ThreeVector.h:92
BinaryOpTraits< Value, Value >::MulT Value2
Value squared.
Definition: ThreeVector.h:37
double angle(const ThreeVector< U > &v) const
Angle between two vectors.
Definition: ThreeVector.h:286
ThreeVector< Value > & rotate(double angle, const ThreeVector< U > &axis)
Apply a rotation.
Definition: ThreeVector.h:180
BinaryOpTraits< Value2, Value2 >::MulT Value4
Value to the 4th power.
Definition: ThreeVector.h:39
void ounitstream(OStream &os, const vector< T, Alloc > &v, UT &u)
Ouput a vector of objects with the specified unit.
Definition: Containers.h:275
A 3-component vector.
Definition: ThreeVector.h:33
This is the main namespace within which all identifiers in ThePEG are declared.
Definition: FactoryBase.h:28
void setPhi(double ph)
Set the azimuthal angle.
Definition: ThreeVector.h:133
double theta() const
Polar angle.
Definition: ThreeVector.h:113
void iunitstream(IStream &is, vector< T, Alloc > &v, UT &u)
Input a vector of objects with the specified unit.
Definition: Containers.h:289
double deltaPhi(const ThreeVector< U > &v2) const
Azimuthal angle difference, brought into the range .
Definition: ThreeVector.h:164
double cosTheta(const ThreeVector< U > &q) const
Cosine of the azimuthal angle between two vectors.
Definition: ThreeVector.h:273
BinaryOpTraits< Value, U >::MulT dot(const ThreeVector< U > &a) const
Dot product.
Definition: ThreeVector.h:86
Value perp(const ThreeVector< U > &p) const
Transverse component with respect to the given axis.
Definition: ThreeVector.h:106
ThreeVector< Value > orthogonal() const
Orthogonal vector.
Definition: ThreeVector.h:149
const 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:90
ThreeVector< typename BinaryOpTraits< Value, U >::MulT > cross(const ThreeVector< U > &a) const
Vector cross-product.
Definition: ThreeVector.h:232
Value perp() const
Transverse component .
Definition: ThreeVector.h:81
ThreeVector< Value > & rotateUz(const Axis &axis)
Rotate the reference frame to a new z-axis.
Definition: ThreeVector.h:209
const double pi
Good old .
Definition: Constants.h:54
double phi() const
Azimuthal angle.
Definition: ThreeVector.h:119
Value2 mag2() const
Squared magnitude .
Definition: ThreeVector.h:72
ThreeVector< double > unit() const
Parallel vector with unit length.
Definition: ThreeVector.h:141
Value2 perp2() const
Squared transverse component .
Definition: ThreeVector.h:78
BinaryOpTraits should be specialized with typdefs called MulT and DivT which gives the type resulting...
Definition: PhysicalQty.h:239
void setTheta(double th)
Set the polar angle.
Definition: ThreeVector.h:124
Value mag() const
Magnitude .
Definition: ThreeVector.h:75
ThreeVector< double > unitVector(const ThreeVector< Value > &v)
A parallel vector with unit length.
Definition: ThreeVector.h:364
IUnit< T, UT > iunit(T &t, const UT &ut)
Helper function creating a IUnit object given an object and a unit.
Definition: UnitIO.h:97