include/HepMC3/FourVector.h

0001 // -*- C++ -*-
0002 //
0003 // This file is part of HepMC
0004 // Copyright (C) 2014-2023 The HepMC collaboration (see AUTHORS for details)
0005 //
0006 #ifndef HEPMC3_FOURVECTOR_H
0007 #define HEPMC3_FOURVECTOR_H
0008 /**
0009  *  @file FourVector.h
0010  *  @brief Definition of \b class FourVector
0011  */
0012 #include <cmath>
0013 #include <limits>
0014 #ifndef M_PI
0015 /** @brief Definition of PI. Needed on some platforms */
0016 #define M_PI 3.14159265358979323846264338327950288
0017 #endif
0018 namespace HepMC3 {
0019
0020
0021 /**
0022  *  @brief Generic 4-vector
0023  *
0024  *  Interpretation of its content depends on accessors used: it's much simpler to do this
0025  *  than to distinguish between space and momentum vectors via the type system (especially
0026  *  given the need for backward compatibility with HepMC2). Be sensible and don't call
0027  *  energy functions on spatial vectors! To avoid duplication, most definitions are only
0028  *  implemented on the spatial function names, with the energy-momentum functions as aliases.
0029  *
0030  *  This is @a not intended to be a fully featured 4-vector, but does contain the majority
0031  *  of common non-boosting functionality, as well as a few support operations on
0032  *  4-vectors.
0033  *
0034  *  The implementations in this class are fully inlined.
0035  */
0036 class FourVector {
0037 public:
0038
0039     /** @brief Default constructor */
0040     FourVector()
0041         : m_v1(0.0),   m_v2(0.0), m_v3(0.0),    m_v4(0.0)  {}
0042     /** @brief Sets all FourVector fields */
0043     FourVector(double xx, double yy, double zz, double ee)
0044         : m_v1(xx),     m_v2(yy),   m_v3(zz),      m_v4(ee)    {}
0045     /** @brief Copy constructor */
0046     FourVector(const FourVector &) = default;
0047     /** @brief Move constructor */
0048     FourVector(FourVector && ) = default;
0049     /** @brief = */
0050     FourVector& operator=(const FourVector&) = default;
0051     /** @brief = */
0052     FourVector& operator=(FourVector&&) = default;
0053
0054     /// @name Component accessors
0055     /// @{
0056
0057     /** @brief Set all FourVector fields, in order x,y,z,t */
0058     void set(double x1, double x2, double x3, double x4) {
0059         m_v1 = x1;
0060         m_v2 = x2;
0061         m_v3 = x3;
0062         m_v4 = x4;
0063     }
0064
0065     /// set component of position/displacement
0066     void set_component(const int i, const double x)
0067     {
0068     if (i==0) {m_v1=x; return; }
0069     if (i==1) {m_v2=x; return; }
0070     if (i==2) {m_v3=x; return; }
0071     if (i==3) {m_v4=x; return; }
0072     }
0073     /// get component of position/displacement
0074     double get_component(const int i)   const
0075     {
0076     if (i==0) return m_v1;
0077     if (i==1) return m_v2;
0078     if (i==2) return m_v3;
0079     if (i==3) return m_v4;
0080     return 0.0;
0081     }
0082
0083
0084     /// x-component of position/displacement
0085     double x()        const { return m_v1; }
0086     /// Set x-component of position/displacement
0087     void set_x(double xx)   { m_v1 = xx;    }
0088     /// @deprecated Prefer the HepMC-style set_x() function
0089     void setX(double xx)   { set_x(xx);    }
0090
0091     /// y-component of position/displacement
0092     double y()        const { return m_v2; }
0093     /// Set y-component of position/displacement
0094     void   set_y(double yy)   { m_v2 = yy;    }
0095     /// @deprecated Prefer the HepMC-style set_y() function
0096     void setY(double yy)   { set_y(yy);    }
0097
0098     /// z-component of position/displacement
0099     double z()        const { return m_v3; }
0100     /// Set z-component of position/displacement
0101     void set_z(double zz)   { m_v3 = zz;    }
0102     /// @deprecated Prefer the HepMC-style set_z() function
0103     void setZ(double zz)   { set_z(zz);    }
0104
0105     /// Time component of position/displacement
0106     double t()        const { return m_v4; }
0107     /// Set time component of position/displacement
0108     void set_t(double tt)   { m_v4 = tt;    }
0109     /// @deprecated Prefer the HepMC-style set_t() function
0110     void setT(double tt)   { set_t(tt);    }
0111
0112
0113     /// x-component of momentum
0114     double px() const { return x(); }
0115     /// Set x-component of momentum
0116     void set_px(double pxx) { set_x(pxx);   }
0117     /// @deprecated Prefer the HepMC-style set_px() function
0118     void setPx(double pxx) { set_px(pxx);    }
0119
0120     /// y-component of momentum
0121     double py() const { return y(); }
0122     /// Set y-component of momentum
0123     void set_py(double pyy) { set_y(pyy);   }
0124     /// @deprecated Prefer the HepMC-style set_py() function
0125     void setPy(double pyy) { set_py(pyy);    }
0126
0127     /// z-component of momentum
0128     double pz() const { return z(); }
0129     /// Set z-component of momentum
0130     void set_pz(double pzz) { set_z(pzz);   }
0131     /// @deprecated Prefer the HepMC-style set_pz() function
0132     void setPz(double pzz) { set_pz(pzz);    }
0133
0134     /// Energy component of momentum
0135     double e() const { return t(); }
0136     /// Set energy component of momentum
0137     void set_e(double ee ) { this->set_t(ee);    }
0138     /// @deprecated Prefer the HepMC-style set_y() function
0139     void setE(double ee) { set_e(ee);    }
0140
0141     /// @}
0142
0143
0144     /// @name Computed properties
0145     /// @{
0146
0147     /// Squared magnitude of (x, y, z) 3-vector
0148     double length2()  const { return x()*x() + y()*y() + z()*z(); }
0149     /// Magnitude of spatial (x, y, z) 3-vector
0150     double length()  const { return std::sqrt(length2()); }
0151     /// Magnitude of spatial (x, y, z) 3-vector, for HepMC2 compatibility
0152     double rho()    const { return length(); }
0153     /// Squared magnitude of (x, y) vector
0154     double perp2()      const { return  x()*x() + y()*y(); }
0155     /// Magnitude of (x, y) vector
0156     double perp()      const { return std::sqrt(perp2()); }
0157     /// Spacetime invariant interval s^2 = t^2 - x^2 - y^2 - z^2
0158     double interval() const { return t()*t() - length2(); }
0159
0160     /// Squared magnitude of p3 = (px, py, pz) vector
0161     double p3mod2()       const { return length2(); }
0162     /// Magnitude of p3 = (px, py, pz) vector
0163     double p3mod()       const { return length(); }
0164     /// Squared transverse momentum px^2 + py^2
0165     double pt2()      const { return perp2(); }
0166     /// Transverse momentum
0167     double pt()      const { return perp(); }
0168     /// Squared invariant mass m^2 = E^2 - px^2 - py^2 - pz^2
0169     double m2()       const { return interval(); }
0170     /// Invariant mass. Returns -sqrt(-m) if e^2 - P^2 is negative
0171     double m() const { return (m2() > 0.0) ? std::sqrt(m2()) : -std::sqrt(-m2()); }
0172
0173     /// Azimuthal angle
0174     double phi() const { return std::atan2( y(), x() ); }
0175     /// Polar angle w.r.t. z direction
0176     double theta() const {  return std::atan2( perp(), z() ); }
0177     /// Pseudorapidity
0178     double eta() const {
0179       if ( p3mod() == 0.0 )  return 0.0;
0180       if ( p3mod() == fabs(pz()) ) return std::copysign(HUGE_VAL, pz());
0181       return 0.5*std::log( (p3mod() + pz()) / (p3mod() - pz()) );
0182     }
0183     /// Rapidity
0184     double rap() const {
0185       if ( e() == 0.0 ) return 0.0;
0186       if ( e() == fabs(pz()) ) return std::copysign(HUGE_VAL, pz());
0187       return 0.5*std::log( (e() + pz()) / (e() - pz()) );
0188     }
0189     /// Absolute pseudorapidity
0190     double abs_eta() const { return std::abs( eta() ); }
0191     /// Absolute rapidity
0192     double abs_rap() const { return std::abs( rap() ); }
0193
0194     /// Same as eta()
0195     ///
0196     /// @deprecated Prefer 'only one way to do it', and we don't have equivalent long names for e.g. pid, phi or eta
0197     double pseudoRapidity() const { return eta(); }
0198
0199     /// @}
0200
0201
0202     /// @name Comparisons to another FourVector
0203     /// @{
0204
0205     /// Check if the length of this vertex is zero
0206     bool is_zero() const { return x() == 0 && y() == 0 && z() == 0 && t() == 0; }
0207
0208     /// Signed azimuthal angle separation in [-pi, pi]
0209     double delta_phi(const FourVector &v) const {
0210         double dphi = phi() - v.phi();
0211         if (dphi != dphi) return dphi;
0212         while (dphi >=  M_PI) dphi -= 2.*M_PI;
0213         while (dphi <  -M_PI) dphi += 2.*M_PI;
0214         return dphi;
0215     }
0216
0217     /// Pseudorapidity separation
0218     double delta_eta(const FourVector &v) const { return eta() - v.eta(); }
0219
0220     /// Rapidity separation
0221     double delta_rap(const FourVector &v) const { return rap() - v.rap(); }
0222
0223     /// R_eta^2-distance separation dR^2 = dphi^2 + deta^2
0224     double delta_r2_eta(const FourVector &v) const {
0225         return delta_phi(v)*delta_phi(v) + delta_eta(v)*delta_eta(v);
0226     }
0227
0228     /// R_eta-distance separation dR = sqrt(dphi^2 + deta^2)
0229     double delta_r_eta(const FourVector &v) const {
0230         return std::sqrt( delta_r2_eta(v) );
0231     }
0232
0233     /// R_rap^2-distance separation dR^2 = dphi^2 + drap^2
0234     double delta_r2_rap(const FourVector &v) const {
0235         return delta_phi(v)*delta_phi(v) + delta_rap(v)*delta_rap(v);
0236     }
0237
0238     /// R-rap-distance separation dR = sqrt(dphi^2 + drap^2)
0239     double delta_r_rap(const FourVector &v) const {
0240         return std::sqrt( delta_r2_rap(v) );
0241     }
0242
0243     /// @}
0244
0245
0246     /// @name Operators
0247     /// @{
0248
0249     /// Equality
0250     bool operator==(const FourVector& rhs) const {
0251         return x() == rhs.x() && y() == rhs.y() && z() == rhs.z() && t() == rhs.t();
0252     }
0253     /// Inequality
0254     bool operator!=(const FourVector& rhs) const { return !(*this == rhs); }
0255
0256     /// Arithmetic operator +
0257     FourVector  operator+ (const FourVector& rhs) const {
0258         return FourVector( x() + rhs.x(), y() + rhs.y(), z() + rhs.z(), t() + rhs.t() );
0259     }
0260     /// Arithmetic operator -
0261     FourVector  operator- (const FourVector& rhs) const {
0262         return FourVector( x() - rhs.x(), y() - rhs.y(), z() - rhs.z(), t() - rhs.t() );
0263     }
0264     /// Arithmetic operator * by scalar
0265     FourVector  operator* (const double rhs) const {
0266         return FourVector( x()*rhs, y()*rhs, z()*rhs, t()*rhs );
0267     }
0268     /// Arithmetic operator / by scalar
0269     FourVector  operator/ (const double rhs) const {
0270         return FourVector( x()/rhs, y()/rhs, z()/rhs, t()/rhs );
0271     }
0272
0273     /// Arithmetic operator +=
0274     void operator += (const FourVector& rhs) {
0275         setX(x() + rhs.x());
0276         setY(y() + rhs.y());
0277         setZ(z() + rhs.z());
0278         setT(t() + rhs.t());
0279     }
0280     /// Arithmetic operator -=
0281     void operator -= (const FourVector& rhs) {
0282         setX(x() - rhs.x());
0283         setY(y() - rhs.y());
0284         setZ(z() - rhs.z());
0285         setT(t() - rhs.t());
0286     }
0287     /// Arithmetic operator *= by scalar
0288     void operator *= (const double rhs) {
0289         setX(x()*rhs);
0290         setY(y()*rhs);
0291         setZ(z()*rhs);
0292         setT(t()*rhs);
0293     }
0294     /// Arithmetic operator /= by scalar
0295     void operator /= (const double rhs) {
0296         setX(x()/rhs);
0297         setY(y()/rhs);
0298         setZ(z()/rhs);
0299         setT(t()/rhs);
0300     }
0301
0302     /// @}
0303
0304
0305     /// Static null FourVector = (0,0,0,0)
0306     static const FourVector& ZERO_VECTOR() {
0307         static const FourVector v;
0308         return v;
0309     }
0310
0311
0312 private:
0313
0314     double m_v1; ///< px or x. Interpretation depends on accessors used
0315     double m_v2; ///< py or y. Interpretation depends on accessors used
0316     double m_v3; ///< pz or z. Interpretation depends on accessors used
0317     double m_v4; ///< e  or t. Interpretation depends on accessors used
0318
0319 };
0320
0321
0322 /// @name Unbound vector comparison functions
0323 /// @{
0324
0325 /// Signed azimuthal angle separation in [-pi, pi] between vecs @c a and @c b
0326 inline double delta_phi(const FourVector &a, const FourVector &b) { return b.delta_phi(a); }
0327
0328 /// Pseudorapidity separation between vecs @c a and @c b
0329 inline double delta_eta(const FourVector &a, const FourVector &b) { return b.delta_eta(a); }
0330
0331 /// Rapidity separation between vecs @c a and @c b
0332 inline double delta_rap(const FourVector &a, const FourVector &b) { return b.delta_rap(a); }
0333
0334 /// R_eta^2-distance separation dR^2 = dphi^2 + deta^2 between vecs @c a and @c b
0335 inline double delta_r2_eta(const FourVector &a, const FourVector &b) { return b.delta_r2_eta(a); }
0336
0337 /// R_eta-distance separation dR = sqrt(dphi^2 + deta^2) between vecs @c a and @c b
0338 inline double delta_r_eta(const FourVector &a, const FourVector &b) { return b.delta_r_eta(a); }
0339
0340 /// R_rap^2-distance separation dR^2 = dphi^2 + drap^2 between vecs @c a and @c b
0341 inline double delta_r2_rap(const FourVector &a, const FourVector &b) { return b.delta_r2_rap(a); }
0342
0343 /// R_rap-distance separation dR = sqrt(dphi^2 + drap^2) between vecs @c a and @c b
0344 inline double delta_r_rap(const FourVector &a, const FourVector &b) { return b.delta_r_rap(a); }
0345
0346 /// @}
0347
0348
0349 } // namespace HepMC3
0350 #endif