Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/Rivet/Math/Vector3.hh was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 #ifndef RIVET_MATH_VECTOR3
0002 #define RIVET_MATH_VECTOR3
0003 
0004 #include "Rivet/Tools/TypeTraits.hh"
0005 #include "Rivet/Math/MathConstants.hh"
0006 #include "Rivet/Math/MathUtils.hh"
0007 #include "Rivet/Math/VectorN.hh"
0008 
0009 namespace Rivet {
0010 
0011 
0012   class Vector3;
0013   typedef Vector3 ThreeVector;
0014   typedef Vector3 V3;
0015   Vector3 multiply(const double, const Vector3&);
0016   Vector3 multiply(const Vector3&, const double);
0017   Vector3 add(const Vector3&, const Vector3&);
0018   Vector3 operator*(const double, const Vector3&);
0019   Vector3 operator*(const Vector3&, const double);
0020   Vector3 operator/(const Vector3&, const double);
0021   Vector3 operator+(const Vector3&, const Vector3&);
0022   Vector3 operator-(const Vector3&, const Vector3&);
0023 
0024   class ThreeMomentum;
0025   typedef ThreeMomentum P3;
0026   ThreeMomentum multiply(const double, const ThreeMomentum&);
0027   ThreeMomentum multiply(const ThreeMomentum&, const double);
0028   ThreeMomentum add(const ThreeMomentum&, const ThreeMomentum&);
0029   ThreeMomentum operator*(const double, const ThreeMomentum&);
0030   ThreeMomentum operator*(const ThreeMomentum&, const double);
0031   ThreeMomentum operator/(const ThreeMomentum&, const double);
0032   ThreeMomentum operator+(const ThreeMomentum&, const ThreeMomentum&);
0033   ThreeMomentum operator-(const ThreeMomentum&, const ThreeMomentum&);
0034 
0035   class Matrix3;
0036 
0037 
0038 
0039   /// @brief Three-dimensional specialisation of Vector.
0040   class Vector3 : public Vector<3> {
0041 
0042     friend class Matrix3;
0043     friend Vector3 multiply(const double, const Vector3&);
0044     friend Vector3 multiply(const Vector3&, const double);
0045     friend Vector3 add(const Vector3&, const Vector3&);
0046     friend Vector3 subtract(const Vector3&, const Vector3&);
0047 
0048   public:
0049     Vector3() : Vector<3>() { }
0050 
0051     template<typename V3TYPE>
0052     Vector3(const V3TYPE& other) {
0053       this->setX(other.x());
0054       this->setY(other.y());
0055       this->setZ(other.z());
0056     }
0057 
0058     Vector3(const Vector<3>& other) {
0059       this->setX(other.get(0));
0060       this->setY(other.get(1));
0061       this->setZ(other.get(2));
0062     }
0063 
0064     Vector3(double x, double y, double z) {
0065       this->setX(x);
0066       this->setY(y);
0067       this->setZ(z);
0068     }
0069 
0070     ~Vector3() { }
0071 
0072 
0073   public:
0074 
0075     static Vector3 mkX() { return Vector3(1,0,0); }
0076     static Vector3 mkY() { return Vector3(0,1,0); }
0077     static Vector3 mkZ() { return Vector3(0,0,1); }
0078 
0079 
0080   public:
0081 
0082     double x() const { return get(0); }
0083     double x2() const { return sqr(x()); }
0084     Vector3& setX(double x) { set(0, x); return *this; }
0085 
0086     double y() const { return get(1); }
0087     double y2() const { return sqr(y()); }
0088     Vector3& setY(double y) { set(1, y); return *this; }
0089 
0090     double z() const { return get(2); }
0091     double z2() const { return sqr(z()); }
0092     Vector3& setZ(double z) { set(2, z); return *this; }
0093 
0094 
0095     /// Dot-product with another vector
0096     double dot(const Vector3& v) const {
0097       return _vec.dot(v._vec);
0098     }
0099 
0100     /// Cross-product with another vector
0101     Vector3 cross(const Vector3& v) const {
0102       Vector3 result;
0103       result._vec = _vec.cross(v._vec);
0104       return result;
0105     }
0106 
0107     /// Angle in radians to another vector
0108     double angle(const Vector3& v) const {
0109       const double localDotOther = unit().dot(v.unit());
0110       if (localDotOther > 1.0) return 0.0;
0111       if (localDotOther < -1.0) return M_PI;
0112       return acos(localDotOther);
0113     }
0114 
0115 
0116     /// Unit-normalized version of this vector.
0117     Vector3 unitVec() const {
0118       double md = mod();
0119       if ( fuzzyLessEquals(md, 0.0) ) return Vector3();
0120       else return *this * 1.0/md;
0121     }
0122 
0123     /// Synonym for unitVec
0124     Vector3 unit() const {
0125       return unitVec();
0126     }
0127 
0128     /// Polar projection of this vector into the x-y plane
0129     Vector3 polarVec() const {
0130       Vector3 rtn = *this;
0131       rtn.setZ(0.);
0132       return rtn;
0133     }
0134     /// Synonym for polarVec
0135     Vector3 perpVec() const {
0136       return polarVec();
0137     }
0138     /// Synonym for polarVec
0139     Vector3 rhoVec() const {
0140       return polarVec();
0141     }
0142 
0143     /// Square of the polar radius (
0144     double polarRadius2() const {
0145       return x()*x() + y()*y();
0146     }
0147     /// Synonym for polarRadius2
0148     double perp2() const {
0149       return polarRadius2();
0150     }
0151     /// Synonym for polarRadius2
0152     double rho2() const {
0153       return polarRadius2();
0154     }
0155 
0156     /// Polar radius
0157     double polarRadius() const {
0158       return sqrt(polarRadius2());
0159     }
0160     /// Synonym for polarRadius
0161     double perp() const {
0162       return polarRadius();
0163     }
0164     /// Synonym for polarRadius
0165     double rho() const {
0166       return polarRadius();
0167     }
0168 
0169     /// @brief Angle subtended by the vector's projection in x-y and the x-axis.
0170     ///
0171     /// @note Returns zero in the case of a vector with null x and y components.
0172     /// @todo Would it be better to return NaN in the null-perp case? Or throw?!
0173     double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const {
0174       // If this has a null perp-vector, return zero rather than let atan2 set an error state
0175       // This isn't necessary if the implementation supports IEEE floating-point arithmetic (IEC 60559)... are we sure?
0176       if (x() == 0 && y() == 0) return 0.0; //< Or return nan / throw an exception?
0177       // Calculate the arctan and return in the requested range
0178       const double value = atan2( y(), x() );
0179       return mapAngle(value, mapping);
0180     }
0181     /// Synonym for azimuthalAngle.
0182     double phi(const PhiMapping mapping = ZERO_2PI) const {
0183       return azimuthalAngle(mapping);
0184     }
0185 
0186     /// Tangent of the polar angle
0187     double tanTheta() const {
0188       return polarRadius()/z();
0189     }
0190 
0191     /// Sine of the polar angle
0192     double sinTheta() const {
0193       return sqrt(polarRadius2()/mod2());
0194     }
0195 
0196     /// Cosine of the polar angle
0197     double cosTheta() const {
0198       return z()/mod2();
0199     }
0200 
0201     /// Angle subtended by the vector and the z-axis.
0202     double polarAngle() const {
0203       // Get number beween [0,PI]
0204       const double polarangle = atan2(polarRadius(), z());
0205       return mapAngle0ToPi(polarangle);
0206     }
0207 
0208     /// Synonym for polarAngle
0209     double theta() const {
0210       return polarAngle();
0211     }
0212 
0213     /// @brief Purely geometric approximation to rapidity
0214     ///
0215     /// eta = -ln[ tan(theta/2) ]
0216     ///
0217     /// Also invariant under z-boosts, equal to y for massless particles.
0218     ///
0219     /// Implemented using the tan half-angle formula
0220     /// tan(theta/2) = sin(theta) / [1 + cos(theta)] = pT / (p + pz)
0221     double pseudorapidity() const {
0222       if (mod() == 0.0) return 0.0; ///< @todo Add [[ unlikely ]] with C++20
0223       if (mod() == fabs(z()) ) return std::copysign(INF, z()); ///< @todo Add [[ unlikely ]] with C++20
0224       const double eta = std::log((mod() + fabs(z())) / perp());
0225       return std::copysign(eta, z());
0226     }
0227 
0228     /// Synonym for pseudorapidity
0229     double eta() const {
0230       return pseudorapidity();
0231     }
0232 
0233     /// Convenience shortcut for fabs(eta())
0234     double abseta() const {
0235       return fabs(eta());
0236     }
0237 
0238 
0239   public:
0240 
0241     /// In-place scalar multiplication operator
0242     Vector3& operator *= (const double a) {
0243       _vec = multiply(a, *this)._vec;
0244       return *this;
0245     }
0246 
0247     /// In-place scalar division operator
0248     Vector3& operator /= (const double a) {
0249       _vec = multiply(1.0/a, *this)._vec;
0250       return *this;
0251     }
0252 
0253     /// In-place addition operator
0254     Vector3& operator += (const Vector3& v) {
0255       _vec = add(*this, v)._vec;
0256       return *this;
0257     }
0258 
0259     /// In-place subtraction operator
0260     Vector3& operator -= (const Vector3& v) {
0261       _vec = subtract(*this, v)._vec;
0262       return *this;
0263     }
0264 
0265     /// In-place negation operator
0266     Vector3 operator - () const {
0267       Vector3 rtn;
0268       rtn._vec = -_vec;
0269       return rtn;
0270     }
0271 
0272   };
0273 
0274 
0275 
0276   /// Unbound dot-product function
0277   inline double dot(const Vector3& a, const Vector3& b) {
0278     return a.dot(b);
0279   }
0280 
0281   /// Unbound cross-product function
0282   inline Vector3 cross(const Vector3& a, const Vector3& b) {
0283     return a.cross(b);
0284   }
0285 
0286   /// Unbound scalar-product function
0287   inline Vector3 multiply(const double a, const Vector3& v) {
0288     Vector3 result;
0289     result._vec = a * v._vec;
0290     return result;
0291   }
0292 
0293   /// Unbound scalar-product function
0294   inline Vector3 multiply(const Vector3& v, const double a) {
0295     return multiply(a, v);
0296   }
0297 
0298   /// Unbound scalar multiplication operator
0299   inline Vector3 operator * (const double a, const Vector3& v) {
0300     return multiply(a, v);
0301   }
0302 
0303   /// Unbound scalar multiplication operator
0304   inline Vector3 operator * (const Vector3& v, const double a) {
0305     return multiply(a, v);
0306   }
0307 
0308   /// Unbound scalar division operator
0309   inline Vector3 operator / (const Vector3& v, const double a) {
0310     return multiply(1.0/a, v);
0311   }
0312 
0313   /// Unbound vector addition function
0314   inline Vector3 add(const Vector3& a, const Vector3& b) {
0315     Vector3 result;
0316     result._vec = a._vec + b._vec;
0317     return result;
0318   }
0319 
0320   /// Unbound vector subtraction function
0321   inline Vector3 subtract(const Vector3& a, const Vector3& b) {
0322     Vector3 result;
0323     result._vec = a._vec - b._vec;
0324     return result;
0325   }
0326 
0327   /// Unbound vector addition operator
0328   inline Vector3 operator + (const Vector3& a, const Vector3& b) {
0329     return add(a, b);
0330   }
0331 
0332   /// Unbound vector subtraction operator
0333   inline Vector3 operator - (const Vector3& a, const Vector3& b) {
0334     return subtract(a, b);
0335   }
0336 
0337   // More physicsy coordinates etc.
0338 
0339   /// Angle (in radians) between two 3-vectors.
0340   inline double angle(const Vector3& a, const Vector3& b) {
0341     return a.angle(b);
0342   }
0343 
0344 
0345   /////////////////////////////////////////////////////
0346 
0347 
0348   /// Specialized version of the ThreeVector with momentum functionality.
0349   class ThreeMomentum : public ThreeVector {
0350   public:
0351     ThreeMomentum() { }
0352 
0353     template<typename V3TYPE, typename std::enable_if<HasXYZ<V3TYPE>::value, int>::type DUMMY=0>
0354     ThreeMomentum(const V3TYPE& other) {
0355       this->setPx(other.x());
0356       this->setPy(other.y());
0357       this->setPz(other.z());
0358     }
0359 
0360     ThreeMomentum(const Vector<3>& other)
0361       : ThreeVector(other) { }
0362 
0363     ThreeMomentum(const double px, const double py, const double pz) {
0364       this->setPx(px);
0365       this->setPy(py);
0366       this->setPz(pz);
0367     }
0368 
0369     ~ThreeMomentum() {}
0370 
0371   public:
0372 
0373 
0374     /// @name Coordinate setters
0375     /// @{
0376 
0377     /// Set x-component of momentum \f$ p_x \f$.
0378     ThreeMomentum& setPx(double px) {
0379       setX(px);
0380       return *this;
0381     }
0382 
0383     /// Set y-component of momentum \f$ p_y \f$.
0384     ThreeMomentum& setPy(double py) {
0385       setY(py);
0386       return *this;
0387     }
0388 
0389     /// Set z-component of momentum \f$ p_z \f$.
0390     ThreeMomentum& setPz(double pz) {
0391       setZ(pz);
0392       return *this;
0393     }
0394 
0395     /// @}
0396 
0397 
0398     /// @name Accessors
0399     /// @{
0400 
0401     /// Get x-component of momentum \f$ p_x \f$.
0402     double px() const { return x(); }
0403     /// Get x-squared \f$ p_x^2 \f$.
0404     double px2() const { return x2(); }
0405 
0406     /// Get y-component of momentum \f$ p_y \f$.
0407     double py() const { return y(); }
0408     /// Get y-squared \f$ p_y^2 \f$.
0409     double py2() const { return y2(); }
0410 
0411     /// Get z-component of momentum \f$ p_z \f$.
0412     double pz() const { return z(); }
0413     /// Get z-squared \f$ p_z^2 \f$.
0414     double pz2() const { return z2(); }
0415 
0416 
0417     /// Get the modulus of the 3-momentum
0418     double p() const { return mod(); }
0419     /// Get the modulus-squared of the 3-momentum
0420     double p2() const { return mod2(); }
0421 
0422 
0423     /// Calculate the transverse momentum vector \f$ \vec{p}_T \f$.
0424     ThreeMomentum pTvec() const {
0425       return polarVec();
0426     }
0427     /// Synonym for pTvec
0428     ThreeMomentum ptvec() const {
0429       return pTvec();
0430     }
0431 
0432     /// Calculate the squared transverse momentum \f$ p_T^2 \f$.
0433     double pT2() const {
0434       return polarRadius2();
0435     }
0436     /// Calculate the squared transverse momentum \f$ p_T^2 \f$.
0437     double pt2() const {
0438       return polarRadius2();
0439     }
0440 
0441     /// Calculate the transverse momentum \f$ p_T \f$.
0442     double pT() const {
0443       return sqrt(pT2());
0444     }
0445     /// Calculate the transverse momentum \f$ p_T \f$.
0446     double pt() const {
0447       return sqrt(pT2());
0448     }
0449 
0450     /// @}
0451 
0452 
0453     ////////////////////////////////////////
0454 
0455 
0456     /// @name Arithmetic operators (needed again for covariant returns)
0457     /// @{
0458 
0459     /// Multiply by a scalar
0460     ThreeMomentum& operator *= (double a) {
0461       _vec = multiply(a, *this)._vec;
0462       return *this;
0463     }
0464 
0465     /// Divide by a scalar
0466     ThreeMomentum& operator /= (double a) {
0467       _vec = multiply(1.0/a, *this)._vec;
0468       return *this;
0469     }
0470 
0471     /// Add two 3-momenta
0472     ThreeMomentum& operator += (const ThreeMomentum& v) {
0473       _vec = add(*this, v)._vec;
0474       return *this;
0475     }
0476 
0477     /// Subtract two 3-momenta
0478     ThreeMomentum& operator -= (const ThreeMomentum& v) {
0479       _vec = add(*this, -v)._vec;
0480       return *this;
0481     }
0482 
0483     /// Multiply all components by -1.
0484     ThreeMomentum operator - () const {
0485       ThreeMomentum result;
0486       result._vec = -_vec;
0487       return result;
0488     }
0489 
0490     // /// Multiply space (i.e. all!) components by -1.
0491     // ThreeMomentum reverse() const {
0492     //   return -*this;
0493     // }
0494 
0495     /// @}
0496 
0497   };
0498 
0499 
0500   inline ThreeMomentum multiply(const double a, const ThreeMomentum& v) {
0501     ThreeMomentum result;
0502     result._vec = a * v._vec;
0503     return result;
0504   }
0505 
0506   inline ThreeMomentum multiply(const ThreeMomentum& v, const double a) {
0507     return multiply(a, v);
0508   }
0509 
0510   inline ThreeMomentum operator*(const double a, const ThreeMomentum& v) {
0511     return multiply(a, v);
0512   }
0513 
0514   inline ThreeMomentum operator*(const ThreeMomentum& v, const double a) {
0515     return multiply(a, v);
0516   }
0517 
0518   inline ThreeMomentum operator/(const ThreeMomentum& v, const double a) {
0519     return multiply(1.0/a, v);
0520   }
0521 
0522   inline ThreeMomentum add(const ThreeMomentum& a, const ThreeMomentum& b) {
0523     ThreeMomentum result;
0524     result._vec = a._vec + b._vec;
0525     return result;
0526   }
0527 
0528   inline ThreeMomentum operator+(const ThreeMomentum& a, const ThreeMomentum& b) {
0529     return add(a, b);
0530   }
0531 
0532   inline ThreeMomentum operator-(const ThreeMomentum& a, const ThreeMomentum& b) {
0533     return add(a, -b);
0534   }
0535 
0536 
0537   /// @todo Mixed-arg operators: better via SFINAE??
0538   /// @note Why *does* this actually cause compiler trouble, given (V3, V3) is a correct sig-match for (V3,P3) and (P3, P3) is not?
0539   inline Vector3 operator+(const ThreeMomentum& a, const Vector3& b) {
0540     return add(static_cast<const Vector3&>(a), b);
0541   }
0542   inline Vector3 operator+(const Vector3& a, const ThreeMomentum& b) {
0543     return add(a, static_cast<const Vector3&>(b));
0544   }
0545 
0546   inline Vector3 operator-(const ThreeMomentum& a, const Vector3& b) {
0547     return add(static_cast<const Vector3&>(a), -b);
0548   }
0549   inline Vector3 operator-(const Vector3& a, const ThreeMomentum& b) {
0550     return add(a, -static_cast<const Vector3&>(b));
0551   }
0552 
0553 
0554 
0555   /////////////////////////////////////////////////////
0556 
0557 
0558   /// @defgroup momutils_vec3_deta \f$ |\Delta eta| \f$ calculations from 3-vectors
0559   /// @{
0560 
0561   /// Calculate the difference in pseudorapidity between two spatial vectors.
0562   inline double deltaEta(const Vector3& a, const Vector3& b, bool sign=false) {
0563     return deltaEta(a.pseudorapidity(), b.pseudorapidity(), sign);
0564   }
0565 
0566   /// Calculate the difference in pseudorapidity between two spatial vectors.
0567   inline double deltaEta(const Vector3& v, double eta2, bool sign=false) {
0568     return deltaEta(v.pseudorapidity(), eta2, sign);
0569   }
0570 
0571   /// Calculate the difference in pseudorapidity between two spatial vectors.
0572   inline double deltaEta(double eta1, const Vector3& v, bool sign=false) {
0573     return deltaEta(eta1, v.pseudorapidity(), sign);
0574   }
0575 
0576   /// @}
0577 
0578 
0579   /// @defgroup momutils_vec3_dphi \f$ \Delta phi \f$ calculations from 3-vectors
0580   /// @{
0581 
0582   /// Calculate the difference in azimuthal angle between two spatial vectors.
0583   inline double deltaPhi(const Vector3& a, const Vector3& b, bool sign=false) {
0584     return deltaPhi(a.azimuthalAngle(), b.azimuthalAngle(), sign);
0585   }
0586 
0587   /// Calculate the difference in azimuthal angle between two spatial vectors.
0588   inline double deltaPhi(const Vector3& v, double phi2, bool sign=false) {
0589     return deltaPhi(v.azimuthalAngle(), phi2, sign);
0590   }
0591 
0592   /// Calculate the difference in azimuthal angle between two spatial vectors.
0593   inline double deltaPhi(double phi1, const Vector3& v, bool sign=false) {
0594     return deltaPhi(phi1, v.azimuthalAngle(), sign);
0595   }
0596 
0597   /// @}
0598 
0599 
0600   /// @defgroup momutils_vec3_dr \f$ \Delta R \f$ calculations from 3-vectors
0601   /// @{
0602 
0603   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
0604   inline double deltaR2(const Vector3& a, const Vector3& b) {
0605     return deltaR2(a.pseudorapidity(), a.azimuthalAngle(),
0606                    b.pseudorapidity(), b.azimuthalAngle());
0607   }
0608 
0609   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
0610   inline double deltaR(const Vector3& a, const Vector3& b) {
0611     return sqrt(deltaR2(a,b));
0612   }
0613 
0614   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
0615   inline double deltaR2(const Vector3& v, double eta2, double phi2) {
0616     return deltaR2(v.pseudorapidity(), v.azimuthalAngle(), eta2, phi2);
0617   }
0618 
0619   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
0620   inline double deltaR(const Vector3& v, double eta2, double phi2) {
0621     return sqrt(deltaR2(v, eta2, phi2));
0622   }
0623 
0624   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
0625   inline double deltaR2(double eta1, double phi1, const Vector3& v) {
0626     return deltaR2(eta1, phi1, v.pseudorapidity(), v.azimuthalAngle());
0627   }
0628 
0629   /// Calculate the 2D rapidity-azimuthal ("eta-phi") distance between two spatial vectors.
0630   inline double deltaR(double eta1, double phi1, const Vector3& v) {
0631     return sqrt(deltaR2(eta1, phi1, v));
0632   }
0633 
0634   /// @}
0635 
0636 
0637   /// @defgroup momutils_vec3_mt MT calculation
0638   /// @{
0639 
0640   /// Calculate transverse mass of a visible and an invisible 3-vector
0641   ///
0642   /// @note Note assumption of zero-mass particles
0643   inline double mT(const Vector3& vis, const Vector3& invis) {
0644     // return sqrt(2*vis.perp()*invis.perp() * (1 - cos(deltaPhi(vis, invis))) );
0645     return mT(vis.perp(), invis.perp(), deltaPhi(vis, invis));
0646   }
0647 
0648   /// Calculate transverse momentum of pair of 3-vectors
0649   inline double pT(const Vector3& a, const Vector3& b) {
0650     return (a+b).perp();
0651   }
0652 
0653   /// @}
0654 
0655 
0656 }
0657 
0658 #endif