rivet is hosted by Hepforge, IPPP Durham
 The Rivet MC analysis system  2.5.2
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
```