9 #ifndef ThePEG_ThreeVector_H
10 #define ThePEG_ThreeVector_H
19 #include "ThreeVector.fh"
20 #include "ThePEG/Config/ThePEG.h"
21 #include "ThePEG/Utilities/UnitIO.h"
32 template <
typename Value>
45 : theX(), theY(), theZ() {}
48 : theX(x), theY(y), theZ(z) {}
50 template<
typename ValueB>
51 ThreeVector(
const ThreeVector<ValueB> & v)
52 : theX(v.x()), theY(v.y()), theZ(v.z()) {}
58 Value x()
const {
return theX; }
59 Value y()
const {
return theY; }
60 Value z()
const {
return theZ; }
65 void setX(Value x) { theX = x; }
66 void setY(Value y) { theY = y; }
67 void setZ(Value z) { theZ = z; }
72 Value2
mag2()
const {
return sqr(x()) + sqr(y()) + sqr(z()); }
75 Value
mag()
const {
return sqrt(
mag2()); }
78 Value2
perp2()
const {
return sqr(x()) + sqr(y()); }
87 return x()*a.x() + y()*a.y() + z()*a.z();
94 const pSqType pMag2 = p.
mag2();
95 assert( pMag2 > pSqType() );
98 Value2 ret =
mag2() - sqr(ss)/pMag2;
105 template <
typename U>
107 return sqrt(
perp2(p));
112 double theta()
const {
114 assert(!(x() == Value() && y() == Value() && z() == Value()));
115 return atan2(
perp(),z());
120 return atan2(y(),x());
127 setX(ma*sin(th)*cos(ph));
128 setY(ma*sin(th)*sin(ph));
144 Value mg = sqrt(mg2);
154 return xx < zz ? ThreeVector<Value>(Value(),z(),-y())
157 return yy < zz ? ThreeVector<Value>(-z(),Value(),x())
163 template <
typename U>
165 double dphi = v2.
phi() -
phi();
179 template <
typename U>
183 const U ll = axis.
mag();
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();
190 setX((ca+(1-ca)*dx*dx) * xx
191 +((1-ca)*dx*dy-sa*dz) * yy
192 +((1-ca)*dx*dz+sa*dy) * zz
194 setY(((1-ca)*dy*dx+sa*dz) * xx
195 +(ca+(1-ca)*dy*dy) * yy
196 +((1-ca)*dy*dz-sa*dx) * zz
198 setZ(((1-ca)*dz*dx-sa*dy) * xx
199 +((1-ca)*dz*dy+sa*dx) * yy
200 +(ca+(1-ca)*dz*dz) * zz
214 double up = u1*u1 + u2*u2;
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 );
230 template <
typename U>
234 return ResultT( y()*a.z()-z()*a.y(),
235 -x()*a.z()+z()*a.x(),
236 x()*a.y()-y()*a.x());
249 ThreeVector<Value> & operator-=(
const ThreeVector<Value> & a) {
256 ThreeVector<Value> & operator*=(
double a) {
263 ThreeVector<Value> & operator/=(
double a) {
272 template <
typename U>
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;
285 template <
typename U>
301 operator<< (ostream & os, const ThreeVector<double> & v)
303 return os <<
'(' << v.x() <<
',' << v.y() <<
',' << v.z() <<
')';
308 template <
typename Value>
309 inline ThreeVector<Value>
310 operator+(ThreeVector<Value> a,
311 const ThreeVector<Value> & b)
316 template <
typename Value>
317 inline ThreeVector<Value>
318 operator-(ThreeVector<Value> a,
319 const ThreeVector<Value> & b)
324 template <
typename Value>
325 inline ThreeVector<Value> operator-(
const ThreeVector<Value> & v) {
326 return ThreeVector<Value>(-v.x(),-v.y(),-v.z());
329 template <
typename Value>
330 inline ThreeVector<Value> operator*(ThreeVector<Value> v,
double a) {
334 template <
typename Value>
335 inline ThreeVector<Value> operator*(
double a, ThreeVector<Value> v) {
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());
346 template <
typename ValueA,
typename ValueB>
347 inline ThreeVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
348 operator*(ThreeVector<ValueA> v, ValueB a) {
354 template <
typename ValueA,
typename ValueB>
355 inline typename BinaryOpTraits<ValueA,ValueB>::MulT
363 template <
typename Value>
370 template <
typename OStream,
typename UT,
typename Value>
376 template <
typename IStream,
typename UT,
typename Value>
Value2 perp2(const ThreeVector< U > &p) const
Squared transverse component with respect to the given axis.
BinaryOpTraits< Value, Value >::MulT Value2
Value squared.
double angle(const ThreeVector< U > &v) const
Angle between two vectors.
ThreeVector< Value > & rotate(double angle, const ThreeVector< U > &axis)
Apply a rotation.
BinaryOpTraits< Value2, Value2 >::MulT Value4
Value to the 4th power.
void ounitstream(OStream &os, const vector< T, Alloc > &v, UT &u)
Ouput a vector of objects with the specified unit.
This is the main namespace within which all identifiers in ThePEG are declared.
void setPhi(double ph)
Set the azimuthal angle.
double theta() const
Polar angle.
void iunitstream(IStream &is, vector< T, Alloc > &v, UT &u)
Input a vector of objects with the specified unit.
double deltaPhi(const ThreeVector< U > &v2) const
Azimuthal angle difference, brought into the range .
double cosTheta(const ThreeVector< U > &q) const
Cosine of the azimuthal angle between two vectors.
BinaryOpTraits< Value, U >::MulT dot(const ThreeVector< U > &a) const
Dot product.
Value perp(const ThreeVector< U > &p) const
Transverse component with respect to the given axis.
ThreeVector< Value > orthogonal() const
Orthogonal vector.
const double twopi
Good old .
OUnit< T, UT > ounit(const T &t, const UT &ut)
Helper function creating a OUnit object given an object and a unit.
ThreeVector< typename BinaryOpTraits< Value, U >::MulT > cross(const ThreeVector< U > &a) const
Vector cross-product.
Value perp() const
Transverse component .
ThreeVector< Value > & rotateUz(const Axis &axis)
Rotate the reference frame to a new z-axis.
const double pi
Good old .
double phi() const
Azimuthal angle.
Value2 mag2() const
Squared magnitude .
ThreeVector< double > unit() const
Parallel vector with unit length.
Value2 perp2() const
Squared transverse component .
BinaryOpTraits should be specialized with typdefs called MulT and DivT which gives the type resulting...
void setTheta(double th)
Set the polar angle.
Value mag() const
Magnitude .
ThreeVector< double > unitVector(const ThreeVector< Value > &v)
A parallel vector with unit length.
IUnit< T, UT > iunit(T &t, const UT &ut)
Helper function creating a IUnit object given an object and a unit.