rivet is hosted by Hepforge, IPPP Durham
Vector3.hh
Go to the documentation of this file.
00001 #ifndef RIVET_MATH_VECTOR3
00002 #define RIVET_MATH_VECTOR3
00003 
00004 #include "Rivet/Math/MathHeader.hh"
00005 #include "Rivet/Math/MathUtils.hh"
00006 #include "Rivet/Math/VectorN.hh"
00007 
00008 namespace Rivet {
00009 
00010 
00011   class Vector3;
00012   typedef Vector3 ThreeVector;
00013   class Matrix3;
00014 
00015   Vector3 multiply(const double, const Vector3&);
00016   Vector3 multiply(const Vector3&, const double);
00017   Vector3 add(const Vector3&, const Vector3&);
00018   Vector3 operator*(const double, const Vector3&);
00019   Vector3 operator*(const Vector3&, const double);
00020   Vector3 operator/(const Vector3&, const double);
00021   Vector3 operator+(const Vector3&, const Vector3&);
00022   Vector3 operator-(const Vector3&, const Vector3&);
00023 
00024 
00025   /// @brief Three-dimensional specialisation of Vector.
00026   class Vector3 : public Vector<3> {
00027 
00028     friend class Matrix3;
00029     friend Vector3 multiply(const double, const Vector3&);
00030     friend Vector3 multiply(const Vector3&, const double);
00031     friend Vector3 add(const Vector3&, const Vector3&);
00032     friend Vector3 subtract(const Vector3&, const Vector3&);
00033 
00034   public:
00035     Vector3() : Vector<3>() { }
00036 
00037     template<typename V3>
00038     Vector3(const V3& other) {
00039       this->setX(other.x());
00040       this->setY(other.y());
00041       this->setZ(other.z());
00042     }
00043 
00044     Vector3(const Vector<3>& other) {
00045       this->setX(other.get(0));
00046       this->setY(other.get(1));
00047       this->setZ(other.get(2));
00048     }
00049 
00050     Vector3(double x, double y, double z) {
00051       this->setX(x);
00052       this->setY(y);
00053       this->setZ(z);
00054     }
00055 
00056     ~Vector3() { }
00057 
00058   public:
00059     static Vector3 mkX() { return Vector3(1,0,0); }
00060     static Vector3 mkY() { return Vector3(0,1,0); }
00061     static Vector3 mkZ() { return Vector3(0,0,1); }
00062 
00063   public:
00064     double x() const { return get(0); }
00065     double y() const { return get(1); }
00066     double z() const { return get(2); }
00067     Vector3& setX(double x) { set(0, x); return *this; }
00068     Vector3& setY(double y) { set(1, y); return *this; }
00069     Vector3& setZ(double z) { set(2, z); return *this; }
00070 
00071     double dot(const Vector3& v) const {
00072       return _vec.dot(v._vec);
00073     }
00074 
00075     Vector3 cross(const Vector3& v) const {
00076       Vector3 result;
00077       result._vec = _vec.cross(v._vec);
00078       return result;
00079     }
00080 
00081     double angle(const Vector3& v) const {
00082       const double localDotOther = unit().dot(v.unit());
00083       if (localDotOther > 1.0) return 0.0;
00084       if (localDotOther < -1.0) return M_PI;
00085       return acos(localDotOther);
00086     }
00087 
00088     Vector3 unit() const {
00089       /// @todo What to do in this situation?
00090       if (isZero()) return *this;
00091       else return *this * 1.0/this->mod();
00092     }
00093 
00094     double polarRadius2() const {
00095       return x()*x() + y()*y();
00096     }
00097 
00098     /// Synonym for polarRadius2
00099     double perp2() const {
00100       return polarRadius2();
00101     }
00102 
00103     /// Synonym for polarRadius2
00104     double rho2() const {
00105       return polarRadius2();
00106     }
00107 
00108     double polarRadius() const {
00109       return sqrt(polarRadius2());
00110     }
00111 
00112     /// Synonym for polarRadius
00113     double perp() const {
00114       return polarRadius();
00115     }
00116 
00117     /// Synonym for polarRadius
00118     double rho() const {
00119       return polarRadius();
00120     }
00121 
00122     /// Angle subtended by the vector's projection in x-y and the x-axis.
00123     double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const {
00124       // If this is a null vector, return zero rather than let atan2 set an error state
00125       if (Rivet::isZero(mod2())) return 0.0;
00126 
00127       // Calculate the arctan and return in the requested range
00128       const double value = atan2( y(), x() );
00129       return mapAngle(value, mapping);
00130     }
00131 
00132     /// Synonym for azimuthalAngle.
00133     double phi(const PhiMapping mapping = ZERO_2PI) const {
00134       return azimuthalAngle(mapping);
00135     }
00136 
00137     /// Angle subtended by the vector and the z-axis.
00138     double polarAngle() const {
00139       // Get number beween [0,PI]
00140       const double polarangle = atan2(polarRadius(), z());
00141       return mapAngle0ToPi(polarangle);
00142     }
00143 
00144     /// Synonym for polarAngle
00145     double theta() const {
00146       return polarAngle();
00147     }
00148 
00149     /// Purely geometric approximation to rapidity; exact for massless particles
00150     /// and in the central region.
00151     // cut-off such that |eta| < log(2/DBL_EPSILON)
00152     double pseudorapidity() const {
00153       const double epsilon = DBL_EPSILON;
00154       double m = mod();
00155       if ( m == 0.0 ) return  0.0;
00156       double pt = max(epsilon*m, perp());
00157       double rap = std::log((m + fabs(z()))/pt);
00158       return z() > 0.0 ? rap: -rap;
00159     }
00160 
00161     /// Synonym for pseudorapidity.
00162     double eta() const {
00163       return pseudorapidity();
00164     }
00165 
00166   public:
00167     Vector3& operator*=(const double a) {
00168       _vec = multiply(a, *this)._vec;
00169       return *this;
00170     }
00171 
00172     Vector3& operator/=(const double a) {
00173       _vec = multiply(1.0/a, *this)._vec;
00174       return *this;
00175     }
00176 
00177     Vector3& operator+=(const Vector3& v) {
00178       _vec = add(*this, v)._vec;
00179       return *this;
00180     }
00181 
00182     Vector3& operator-=(const Vector3& v) {
00183       _vec = subtract(*this, v)._vec;
00184       return *this;
00185     }
00186 
00187     Vector3 operator-() const {
00188       Vector3 rtn;
00189       rtn._vec = -_vec;
00190       return rtn;
00191     }
00192 
00193   };
00194 
00195 
00196 
00197   inline double dot(const Vector3& a, const Vector3& b) {
00198     return a.dot(b);
00199   }
00200 
00201   inline Vector3 cross(const Vector3& a, const Vector3& b) {
00202     return a.cross(b);
00203   }
00204 
00205   inline Vector3 multiply(const double a, const Vector3& v) {
00206     Vector3 result;
00207     result._vec = a * v._vec;
00208     return result;
00209   }
00210 
00211   inline Vector3 multiply(const Vector3& v, const double a) {
00212     return multiply(a, v);
00213   }
00214 
00215   inline Vector3 operator*(const double a, const Vector3& v) {
00216     return multiply(a, v);
00217   }
00218 
00219   inline Vector3 operator*(const Vector3& v, const double a) {
00220     return multiply(a, v);
00221   }
00222 
00223   inline Vector3 operator/(const Vector3& v, const double a) {
00224     return multiply(1.0/a, v);
00225   }
00226 
00227   inline Vector3 add(const Vector3& a, const Vector3& b) {
00228     Vector3 result;
00229     result._vec = a._vec + b._vec;
00230     return result;
00231   }
00232 
00233   inline Vector3 subtract(const Vector3& a, const Vector3& b) {
00234     Vector3 result;
00235     result._vec = a._vec - b._vec;
00236     return result;
00237   }
00238 
00239   inline Vector3 operator+(const Vector3& a, const Vector3& b) {
00240     return add(a, b);
00241   }
00242 
00243   inline Vector3 operator-(const Vector3& a, const Vector3& b) {
00244     return subtract(a, b);
00245   }
00246 
00247   // More physicsy coordinates etc.
00248 
00249   /// Angle (in radians) between two 3-vectors.
00250   inline double angle(const Vector3& a, const Vector3& b) {
00251     return a.angle(b);
00252   }
00253 
00254   /////////////////////////////////////////////////////
00255 
00256   /// @name \f$ |\Delta eta| \f$ calculations from 3-vectors
00257   //@{
00258 
00259   /// Calculate the difference in pseudorapidity between two spatial vectors.
00260   inline double deltaEta(const Vector3& a, const Vector3& b) {
00261     return deltaEta(a.pseudorapidity(), b.pseudorapidity());
00262   }
00263 
00264   /// Calculate the difference in pseudorapidity between two spatial vectors.
00265   inline double deltaEta(const Vector3& v, double eta2) {
00266     return deltaEta(v.pseudorapidity(), eta2);
00267   }
00268 
00269   /// Calculate the difference in pseudorapidity between two spatial vectors.
00270   inline double deltaEta(double eta1, const Vector3& v) {
00271     return deltaEta(eta1, v.pseudorapidity());
00272   }
00273 
00274   //@}
00275 
00276 
00277   /// @name \f$ \Delta phi \f$ calculations from 3-vectors
00278   //@{
00279 
00280   /// Calculate the difference in azimuthal angle between two spatial vectors.
00281   inline double deltaPhi(const Vector3& a, const Vector3& b) {
00282     return deltaPhi(a.azimuthalAngle(), b.azimuthalAngle());
00283   }
00284 
00285   /// Calculate the difference in azimuthal angle between two spatial vectors.
00286   inline double deltaPhi(const Vector3& v, double phi2) {
00287     return deltaPhi(v.azimuthalAngle(), phi2);
00288   }
00289 
00290   /// Calculate the difference in azimuthal angle between two spatial vectors.
00291   inline double deltaPhi(double phi1, const Vector3& v) {
00292     return deltaPhi(phi1, v.azimuthalAngle());
00293   }
00294 
00295   //@}
00296 
00297 
00298   /// @name \f$ \Delta R \f$ calculations from 3-vectors
00299   //@{
00300 
00301   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
00302   inline double deltaR(const Vector3& a, const Vector3& b) {
00303     return deltaR(a.pseudorapidity(), a.azimuthalAngle(),
00304                   b.pseudorapidity(), b.azimuthalAngle());
00305   }
00306 
00307   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
00308   inline double deltaR(const Vector3& v, double eta2, double phi2) {
00309     return deltaR(v.pseudorapidity(), v.azimuthalAngle(), eta2, phi2);
00310   }
00311 
00312   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
00313   inline double deltaR(double eta1, double phi1, const Vector3& v) {
00314     return deltaR(eta1, phi1, v.pseudorapidity(), v.azimuthalAngle());
00315   }
00316 
00317   //@}
00318 
00319 
00320   /// @name Typedefs of vector types to short names
00321   /// @todo Switch canonical and alias names
00322   //@{
00323   //typedef Vector3 V3; //< generic
00324   typedef Vector3 X3; //< spatial
00325   //@}
00326 
00327 
00328 }
00329 
00330 #endif