Rivet
1.8.0
|
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 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 double localDotOther = unit().dot(v.unit()); 00083 if(Rivet::isZero(localDotOther - 1.0)) return 0.0; 00084 return acos( localDotOther ); 00085 } 00086 00087 Vector3 unit() const { 00089 if (isZero()) return *this; 00090 else return *this * 1.0/this->mod(); 00091 } 00092 00093 double polarRadius2() const { 00094 return x()*x() + y()*y(); 00095 } 00096 00098 double perp2() const { 00099 return polarRadius2(); 00100 } 00101 00103 double rho2() const { 00104 return polarRadius2(); 00105 } 00106 00107 double polarRadius() const { 00108 return sqrt(polarRadius2()); 00109 } 00110 00112 double perp() const { 00113 return polarRadius(); 00114 } 00115 00117 double rho() const { 00118 return polarRadius(); 00119 } 00120 00122 double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const { 00123 // If this is a null vector, return zero rather than let atan2 set an error state 00124 if (Rivet::isZero(mod2())) return 0.0; 00125 00126 // Calculate the arctan and correct for numerical boundary cases 00127 double value = atan2( y(), x() ); 00128 if (value > 2*PI || value < -2*PI){ 00129 value = fmod(value, 2*PI); 00130 } 00131 if (value <= -PI) value += 2*PI; 00132 if (value > PI) value -= 2*PI; 00133 00134 // Return in the requested range 00135 switch (mapping) { 00136 case MINUSPI_PLUSPI: 00137 assert(value > -PI && value <= PI); 00138 return value; 00139 case ZERO_2PI: 00140 if (value >= 0) { 00141 assert(value >= 0 && value < 2*PI); 00142 return value; 00143 } else if (Rivet::isZero(value)) { 00144 value = 0.0; 00145 return value; 00146 } else { 00147 value = 2*PI + value; 00148 assert(value >= 0 && value < 2*PI); 00149 return value; 00150 } 00151 default: 00152 throw std::runtime_error("The specified phi mapping scheme is not yet implemented"); 00153 } 00154 } 00155 00157 double phi(const PhiMapping mapping = ZERO_2PI) const { 00158 return azimuthalAngle(mapping); 00159 } 00160 00162 double polarAngle() const { 00163 // Get number beween [0,PI] 00164 double polarangle = atan2(polarRadius(), z()); 00165 assert(polarangle >= -PI && polarangle <= PI); 00166 return polarangle; 00167 } 00168 00170 double theta() const { 00171 return polarAngle(); 00172 } 00173 00176 double pseudorapidity() const { 00177 return -std::log(tan( 0.5 * polarAngle() )); 00178 } 00179 00181 double eta() const { 00182 return pseudorapidity(); 00183 } 00184 00185 public: 00186 Vector3& operator*=(const double a) { 00187 _vec = multiply(a, *this)._vec; 00188 return *this; 00189 } 00190 00191 Vector3& operator/=(const double a) { 00192 _vec = multiply(1.0/a, *this)._vec; 00193 return *this; 00194 } 00195 00196 Vector3& operator+=(const Vector3& v) { 00197 _vec = add(*this, v)._vec; 00198 return *this; 00199 } 00200 00201 Vector3& operator-=(const Vector3& v) { 00202 _vec = subtract(*this, v)._vec; 00203 return *this; 00204 } 00205 00206 Vector3 operator-() const { 00207 Vector3 rtn; 00208 rtn._vec = -_vec; 00209 return rtn; 00210 } 00211 00212 }; 00213 00214 00215 00216 inline double dot(const Vector3& a, const Vector3& b) { 00217 return a.dot(b); 00218 } 00219 00220 inline Vector3 cross(const Vector3& a, const Vector3& b) { 00221 return a.cross(b); 00222 } 00223 00224 inline Vector3 multiply(const double a, const Vector3& v) { 00225 Vector3 result; 00226 result._vec = a * v._vec; 00227 return result; 00228 } 00229 00230 inline Vector3 multiply(const Vector3& v, const double a) { 00231 return multiply(a, v); 00232 } 00233 00234 inline Vector3 operator*(const double a, const Vector3& v) { 00235 return multiply(a, v); 00236 } 00237 00238 inline Vector3 operator*(const Vector3& v, const double a) { 00239 return multiply(a, v); 00240 } 00241 00242 inline Vector3 operator/(const Vector3& v, const double a) { 00243 return multiply(1.0/a, v); 00244 } 00245 00246 inline Vector3 add(const Vector3& a, const Vector3& b) { 00247 Vector3 result; 00248 result._vec = a._vec + b._vec; 00249 return result; 00250 } 00251 00252 inline Vector3 subtract(const Vector3& a, const Vector3& b) { 00253 Vector3 result; 00254 result._vec = a._vec - b._vec; 00255 return result; 00256 } 00257 00258 inline Vector3 operator+(const Vector3& a, const Vector3& b) { 00259 return add(a, b); 00260 } 00261 00262 inline Vector3 operator-(const Vector3& a, const Vector3& b) { 00263 return subtract(a, b); 00264 } 00265 00266 // More physicsy coordinates etc. 00267 00269 inline double angle(const Vector3& a, const Vector3& b) { 00270 return a.angle(b); 00271 } 00272 00274 inline double polarRadius2(const Vector3& v) { 00275 return v.polarRadius2(); 00276 } 00278 inline double perp2(const Vector3& v) { 00279 return v.perp2(); 00280 } 00282 inline double rho2(const Vector3& v) { 00283 return v.rho2(); 00284 } 00285 00287 inline double polarRadius(const Vector3& v) { 00288 return v.polarRadius(); 00289 } 00291 inline double perp(const Vector3& v) { 00292 return v.perp(); 00293 } 00295 inline double rho(const Vector3& v) { 00296 return v.rho(); 00297 } 00298 00299 00302 inline double azimuthalAngle(const Vector3& v, const PhiMapping mapping = ZERO_2PI) { 00303 return v.azimuthalAngle(mapping); 00304 } 00306 inline double phi(const Vector3& v, const PhiMapping mapping = ZERO_2PI) { 00307 return v.phi(mapping); 00308 } 00309 00311 inline double polarAngle(const Vector3& v) { 00312 return v.polarAngle(); 00313 } 00315 inline double theta(const Vector3& v) { 00316 return v.theta(); 00317 } 00318 00320 inline double pseudorapidity(const Vector3& v) { 00321 return v.pseudorapidity(); 00322 } 00324 inline double eta(const Vector3& v) { 00325 return v.eta(); 00326 } 00327 00328 00330 00332 00333 00335 inline double deltaEta(const Vector3& a, const Vector3& b) { 00336 return deltaEta(a.pseudorapidity(), b.pseudorapidity()); 00337 } 00338 00340 inline double deltaEta(const Vector3& v, double eta2) { 00341 return deltaEta(v.pseudorapidity(), eta2); 00342 } 00343 00345 inline double deltaEta(double eta1, const Vector3& v) { 00346 return deltaEta(eta1, v.pseudorapidity()); 00347 } 00348 00350 00351 00353 00354 00356 inline double deltaPhi(const Vector3& a, const Vector3& b) { 00357 return deltaPhi(a.azimuthalAngle(), b.azimuthalAngle()); 00358 } 00359 00361 inline double deltaPhi(const Vector3& v, double phi2) { 00362 return deltaPhi(v.azimuthalAngle(), phi2); 00363 } 00364 00366 inline double deltaPhi(double phi1, const Vector3& v) { 00367 return deltaPhi(phi1, v.azimuthalAngle()); 00368 } 00369 00371 00372 00374 00375 00377 inline double deltaR(const Vector3& a, const Vector3& b) { 00378 return deltaR(a.pseudorapidity(), a.azimuthalAngle(), 00379 b.pseudorapidity(), b.azimuthalAngle()); 00380 } 00381 00383 inline double deltaR(const Vector3& v, double eta2, double phi2) { 00384 return deltaR(v.pseudorapidity(), v.azimuthalAngle(), eta2, phi2); 00385 } 00386 00388 inline double deltaR(double eta1, double phi1, const Vector3& v) { 00389 return deltaR(eta1, phi1, v.pseudorapidity(), v.azimuthalAngle()); 00390 } 00391 00393 00394 } 00395 00396 #endif