00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef ThePEG_ThreeVector_H
00010 #define ThePEG_ThreeVector_H
00011
00019 #include "ThreeVector.fh"
00020 #include "ThePEG/Config/ThePEG.h"
00021 #include "ThePEG/Utilities/UnitIO.h"
00022 #include <cassert>
00023
00024 namespace ThePEG {
00025
00032 template <typename Value>
00033 class ThreeVector
00034 {
00035 private:
00037 typedef typename BinaryOpTraits<Value,Value>::MulT Value2;
00039 typedef typename BinaryOpTraits<Value2,Value2>::MulT Value4;
00040
00041 public:
00044 ThreeVector()
00045 : theX(), theY(), theZ() {}
00046
00047 ThreeVector(Value x, Value y, Value z)
00048 : theX(x), theY(y), theZ(z) {}
00049
00050 template<typename ValueB>
00051 ThreeVector(const ThreeVector<ValueB> & v)
00052 : theX(v.x()), theY(v.y()), theZ(v.z()) {}
00054
00055 public:
00057
00058 Value x() const { return theX; }
00059 Value y() const { return theY; }
00060 Value z() const { return theZ; }
00062
00064
00065 void setX(Value x) { theX = x; }
00066 void setY(Value y) { theY = y; }
00067 void setZ(Value z) { theZ = z; }
00069
00070 public:
00072 Value2 mag2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
00073
00075 Value mag() const { return sqrt(mag2()); }
00076
00078 Value2 perp2() const { return sqr(x()) + sqr(y()); }
00079
00081 Value perp() const { return sqrt(perp2()); }
00082
00084 template <typename U>
00085 typename BinaryOpTraits<Value,U>::MulT
00086 dot(const ThreeVector<U> & a) const {
00087 return x()*a.x() + y()*a.y() + z()*a.z();
00088 }
00089
00091 template <typename U>
00092 Value2 perp2(const ThreeVector<U> & p) const {
00093 typedef typename BinaryOpTraits<U,U>::MulT pSqType;
00094 const pSqType pMag2 = p.mag2();
00095 assert( pMag2 > pSqType() );
00096 typename BinaryOpTraits<Value,U>::MulT
00097 ss = this->dot(p);
00098 Value2 ret = mag2() - sqr(ss)/pMag2;
00099 if ( ret <= Value2() )
00100 ret = Value2();
00101 return ret;
00102 }
00103
00105 template <typename U>
00106 Value perp(const ThreeVector<U> & p) const {
00107 return sqrt(perp2(p));
00108 }
00109
00111
00112
00113 double theta() const {
00114 assert(!(x() == Value() && y() == Value() && z() == Value()));
00115 return atan2(perp(),z());
00116 }
00117
00119 double phi() const {
00120 if ( x() == Value() && y() == Value() )
00121 return 0.;
00122
00123 return atan2(y(),x());
00124 }
00125
00127 void setTheta(double th) {
00128 double ma = mag();
00129 double ph = phi();
00130 setX(ma*sin(th)*cos(ph));
00131 setY(ma*sin(th)*sin(ph));
00132 setZ(ma*cos(th));
00133 }
00134
00136 void setPhi(double ph) {
00137 double xy = perp();
00138 setX(xy*cos(ph));
00139 setY(xy*sin(ph));
00140 }
00142
00144 ThreeVector<double> unit() const {
00145 Value2 mg2 = mag2();
00146 assert(mg2 > Value2());
00147 Value mg = sqrt(mg2);
00148 return ThreeVector<double>(x()/mg, y()/mg, z()/mg);
00149 }
00150
00152 ThreeVector<Value> orthogonal() const {
00153 Value xx = abs(x());
00154 Value yy = abs(y());
00155 Value zz = abs(z());
00156 if (xx < yy) {
00157 return xx < zz ? ThreeVector<Value>(Value(),z(),-y())
00158 : ThreeVector<Value>(y(),-x(),Value());
00159 } else {
00160 return yy < zz ? ThreeVector<Value>(-z(),Value(),x())
00161 : ThreeVector<Value>(y(),-x(),Value());
00162 }
00163 }
00164
00166 template <typename U>
00167 double deltaPhi (const ThreeVector<U> & v2) const {
00168 double dphi = v2.phi() - phi();
00169 if ( dphi > Constants::pi ) {
00170 dphi -= Constants::twopi;
00171 } else if ( dphi <= -Constants::pi ) {
00172 dphi += Constants::twopi;
00173 }
00174 return dphi;
00175 }
00176
00182 ThreeVector<Value> & rotate(double angle, const ThreeVector<double> & axis) {
00183 if (angle == 0.0)
00184 return *this;
00185 const double ll = axis.mag();
00186 assert(ll > 0.0);
00187
00188 const double sa = sin(angle), ca = cos(angle);
00189 const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
00190 const double xx = x(), yy = y(), zz = z();
00191
00192 setX((ca+(1-ca)*dx*dx) * xx
00193 +((1-ca)*dx*dy-sa*dz) * yy
00194 +((1-ca)*dx*dz+sa*dy) * zz
00195 );
00196 setY(((1-ca)*dy*dx+sa*dz) * xx
00197 +(ca+(1-ca)*dy*dy) * yy
00198 +((1-ca)*dy*dz-sa*dx) * zz
00199 );
00200 setZ(((1-ca)*dz*dx-sa*dy) * xx
00201 +((1-ca)*dz*dy+sa*dx) * yy
00202 +(ca+(1-ca)*dz*dz) * zz
00203 );
00204 return *this;
00205 }
00206
00207
00211 ThreeVector<Value> & rotateUz (const Axis & axis) {
00212 Axis ax = axis.unit();
00213 double u1 = ax.x();
00214 double u2 = ax.y();
00215 double u3 = ax.z();
00216 double up = u1*u1 + u2*u2;
00217 if (up>0) {
00218 up = sqrt(up);
00219 Value px = x(), py = y(), pz = z();
00220 setX( (u1*u3*px - u2*py)/up + u1*pz );
00221 setY( (u2*u3*px + u1*py)/up + u2*pz );
00222 setZ( -up*px + u3*pz );
00223 }
00224 else if (u3 < 0.) {
00225 setX(-x());
00226 setZ(-z());
00227 }
00228 return *this;
00229 }
00230
00232 template <typename U>
00233 ThreeVector<typename BinaryOpTraits<Value,U>::MulT>
00234 cross(const ThreeVector<U> & a) const {
00235 typedef ThreeVector<typename BinaryOpTraits<Value,U>::MulT> ResultT;
00236 return ResultT( y()*a.z()-z()*a.y(),
00237 -x()*a.z()+z()*a.x(),
00238 x()*a.y()-y()*a.x());
00239 }
00240
00241 public:
00243
00244 ThreeVector<Value> & operator+=(const ThreeVector<Value> & a) {
00245 theX += a.x();
00246 theY += a.y();
00247 theZ += a.z();
00248 return *this;
00249 }
00250
00251 ThreeVector<Value> & operator-=(const ThreeVector<Value> & a) {
00252 theX -= a.x();
00253 theY -= a.y();
00254 theZ -= a.z();
00255 return *this;
00256 }
00257
00258 ThreeVector<Value> & operator*=(double a) {
00259 theX *= a;
00260 theY *= a;
00261 theZ *= a;
00262 return *this;
00263 }
00264
00265 ThreeVector<Value> & operator/=(double a) {
00266 theX /= a;
00267 theY /= a;
00268 theZ /= a;
00269 return *this;
00270 }
00272
00274 template <typename U>
00275 double cosTheta(const ThreeVector<U> & q) const {
00276 typedef typename BinaryOpTraits<Value,U>::MulT
00277 ProdType;
00278 ProdType ptot = mag()*q.mag();
00279 assert( ptot > ProdType() );
00280 double arg = dot(q)/ptot;
00281 if (arg > 1.0) arg = 1.0;
00282 else if(arg < -1.0) arg = -1.0;
00283 return arg;
00284 }
00285
00287 template <typename U>
00288 double angle(const ThreeVector<U> & v) const {
00289 return acos(cosTheta(v));
00290 }
00291
00292 private:
00294
00295 Value theX;
00296 Value theY;
00297 Value theZ;
00299 };
00300
00302 inline ostream &
00303 operator<< (ostream & os, const ThreeVector<double> & v)
00304 {
00305 return os << '(' << v.x() << ',' << v.y() << ',' << v.z() << ')';
00306 }
00307
00309
00310 template <typename Value>
00311 inline ThreeVector<Value>
00312 operator+(ThreeVector<Value> a,
00313 const ThreeVector<Value> & b)
00314 {
00315 return a += b;
00316 }
00317
00318 template <typename Value>
00319 inline ThreeVector<Value>
00320 operator-(ThreeVector<Value> a,
00321 const ThreeVector<Value> & b)
00322 {
00323 return a -= b;
00324 }
00325
00326 template <typename Value>
00327 inline ThreeVector<Value> operator-(const ThreeVector<Value> & v) {
00328 return ThreeVector<Value>(-v.x(),-v.y(),-v.z());
00329 }
00330
00331 template <typename Value>
00332 inline ThreeVector<Value> operator*(ThreeVector<Value> v, double a) {
00333 return v *= a;
00334 }
00335
00336 template <typename Value>
00337 inline ThreeVector<Value> operator*(double a, ThreeVector<Value> v) {
00338 return v *= a;
00339 }
00340
00341 template <typename ValueA, typename ValueB>
00342 inline ThreeVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
00343 operator*(ValueB a, ThreeVector<ValueA> v) {
00344 typedef typename BinaryOpTraits<ValueA,ValueB>::MulT ResultT;
00345 return ThreeVector<ResultT>(a*v.x(), a*v.y(), a*v.z());
00346 }
00347
00348 template <typename ValueA, typename ValueB>
00349 inline ThreeVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
00350 operator*(ThreeVector<ValueA> v, ValueB a) {
00351 return a*v;
00352 }
00354
00356 template <typename ValueA, typename ValueB>
00357 inline typename BinaryOpTraits<ValueA,ValueB>::MulT
00358 operator*(const ThreeVector<ValueA> & a,
00359 const ThreeVector<ValueB> & b)
00360 {
00361 return a.dot(b);
00362 }
00363
00365 template <typename Value>
00366 ThreeVector<double> unitVector(const ThreeVector<Value> & v) {
00367 return v.unit();
00368 }
00369
00370
00372 template <typename OStream, typename UT, typename Value>
00373 void ounitstream(OStream & os, const ThreeVector<Value> & p, UT & u) {
00374 os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u);
00375 }
00376
00378 template <typename IStream, typename UT, typename Value>
00379 void iunitstream(IStream & is, ThreeVector<Value> & p, UT & u) {
00380 Value x, y, z;
00381 is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u);
00382 p = ThreeVector<Value>(x, y, z);
00383 }
00384
00385 }
00386
00387 #endif
