FourVector.h

// -*- C++ -*-
//
// This file is part of HepMC
// Copyright (C) 2014-2023 The HepMC collaboration (see AUTHORS for details)
//
#ifndef HEPMC3_FOURVECTOR_H
#define HEPMC3_FOURVECTOR_H
#include <cmath>
#include <limits>
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif
namespace HepMC3 {
 

class FourVector {
public:

    FourVector()
        : m_v1(0.0),   m_v2(0.0), m_v3(0.0),    m_v4(0.0)  {}

    FourVector(double xx, double yy, double zz, double ee)
        : m_v1(xx),     m_v2(yy),   m_v3(zz),      m_v4(ee)    {}

    FourVector(const FourVector &) = default;
    FourVector(FourVector && ) = default;
    FourVector& operator=(const FourVector&) = default;
    FourVector& operator=(FourVector&&) = default;


    void set(double x1, double x2, double x3, double x4) {
        m_v1 = x1;
        m_v2 = x2;
        m_v3 = x3;
        m_v4 = x4;
    }

    void set_component(const int i, const double x)
    {
    if (i==0) {m_v1=x; return; }
    if (i==1) {m_v2=x; return; }
    if (i==2) {m_v3=x; return; }
    if (i==3) {m_v4=x; return; }
    }

    double get_component(const int i)   const
    {
    if (i==0) return m_v1;
    if (i==1) return m_v2;
    if (i==2) return m_v3;
    if (i==3) return m_v4;
    return 0.0;
    }
 

    double x()        const { return m_v1; }
    void set_x(double xx)   { m_v1 = xx;    }
    void setX(double xx)   { set_x(xx);    }

    double y()        const { return m_v2; }
    void   set_y(double yy)   { m_v2 = yy;    }
    void setY(double yy)   { set_y(yy);    }

    double z()        const { return m_v3; }
    void set_z(double zz)   { m_v3 = zz;    }
    void setZ(double zz)   { set_z(zz);    }

    double t()        const { return m_v4; }
    void set_t(double tt)   { m_v4 = tt;    }
    void setT(double tt)   { set_t(tt);    }
 

    double px() const { return x(); }
    void set_px(double pxx) { set_x(pxx);   }
    void setPx(double pxx) { set_px(pxx);    }

    double py() const { return y(); }
    void set_py(double pyy) { set_y(pyy);   }
    void setPy(double pyy) { set_py(pyy);    }

    double pz() const { return z(); }
    void set_pz(double pzz) { set_z(pzz);   }
    void setPz(double pzz) { set_pz(pzz);    }

    double e() const { return t(); }
    void set_e(double ee ) { this->set_t(ee);    }
    void setE(double ee) { set_e(ee);    }

 


    double length2()  const { return x()*x() + y()*y() + z()*z(); }
    double length()  const { return std::sqrt(length2()); }
    double rho()    const { return length(); }
    double perp2()      const { return  x()*x() + y()*y(); }
    double perp()      const { return std::sqrt(perp2()); }
    double interval() const { return t()*t() - length2(); }

    double p3mod2()       const { return length2(); }
    double p3mod()       const { return length(); }
    double pt2()      const { return perp2(); }
    double pt()      const { return perp(); }
    double m2()       const { return interval(); }
    double m() const { return (m2() > 0.0) ? std::sqrt(m2()) : -std::sqrt(-m2()); }

    double phi() const { return std::atan2( y(), x() ); }
    double theta() const {  return std::atan2( perp(), z() ); }
    double eta() const {
      if ( p3mod() == 0.0 )  return 0.0;
      if ( p3mod() == fabs(pz()) ) return std::copysign(HUGE_VAL, pz());
      return 0.5*std::log( (p3mod() + pz()) / (p3mod() - pz()) );
    }

    double rap() const {
      if ( e() == 0.0 ) return 0.0;
      if ( e() == fabs(pz()) ) return std::copysign(HUGE_VAL, pz());
      return 0.5*std::log( (e() + pz()) / (e() - pz()) );
    }

    double abs_eta() const { return std::abs( eta() ); }
    double abs_rap() const { return std::abs( rap() ); }

    double pseudoRapidity() const { return eta(); }

 


    bool is_zero() const { return x() == 0 && y() == 0 && z() == 0 && t() == 0; }

    double delta_phi(const FourVector &v) const {
        double dphi = phi() - v.phi();
        if (dphi != dphi) return dphi;
        while (dphi >=  M_PI) dphi -= 2.*M_PI;
        while (dphi <  -M_PI) dphi += 2.*M_PI;
        return dphi;
    }

    double delta_eta(const FourVector &v) const { return eta() - v.eta(); }

    double delta_rap(const FourVector &v) const { return rap() - v.rap(); }

    double delta_r2_eta(const FourVector &v) const {
        return delta_phi(v)*delta_phi(v) + delta_eta(v)*delta_eta(v);
    }

    double delta_r_eta(const FourVector &v) const {
        return std::sqrt( delta_r2_eta(v) );
    }

    double delta_r2_rap(const FourVector &v) const {
        return delta_phi(v)*delta_phi(v) + delta_rap(v)*delta_rap(v);
    }

    double delta_r_rap(const FourVector &v) const {
        return std::sqrt( delta_r2_rap(v) );
    }

 


    bool operator==(const FourVector& rhs) const {
        return x() == rhs.x() && y() == rhs.y() && z() == rhs.z() && t() == rhs.t();
    }

    bool operator!=(const FourVector& rhs) const { return !(*this == rhs); }

    FourVector  operator+ (const FourVector& rhs) const {
        return FourVector( x() + rhs.x(), y() + rhs.y(), z() + rhs.z(), t() + rhs.t() );
    }

    FourVector  operator- (const FourVector& rhs) const {
        return FourVector( x() - rhs.x(), y() - rhs.y(), z() - rhs.z(), t() - rhs.t() );
    }

    FourVector  operator* (const double rhs) const {
        return FourVector( x()*rhs, y()*rhs, z()*rhs, t()*rhs );
    }

    FourVector  operator/ (const double rhs) const {
        return FourVector( x()/rhs, y()/rhs, z()/rhs, t()/rhs );
    }

    void operator += (const FourVector& rhs) {
        setX(x() + rhs.x());
        setY(y() + rhs.y());
        setZ(z() + rhs.z());
        setT(t() + rhs.t());
    }

    void operator -= (const FourVector& rhs) {
        setX(x() - rhs.x());
        setY(y() - rhs.y());
        setZ(z() - rhs.z());
        setT(t() - rhs.t());
    }

    void operator *= (const double rhs) {
        setX(x()*rhs);
        setY(y()*rhs);
        setZ(z()*rhs);
        setT(t()*rhs);
    }

    void operator /= (const double rhs) {
        setX(x()/rhs);
        setY(y()/rhs);
        setZ(z()/rhs);
        setT(t()/rhs);
    }

 

    static const FourVector& ZERO_VECTOR() {
        static const FourVector v;
        return v;
    }
 
 
private:
 
    double m_v1; 
    double m_v2; 
    double m_v3; 
    double m_v4; 
 
};
 


inline double delta_phi(const FourVector &a, const FourVector &b) { return b.delta_phi(a); }

inline double delta_eta(const FourVector &a, const FourVector &b) { return b.delta_eta(a); }

inline double delta_rap(const FourVector &a, const FourVector &b) { return b.delta_rap(a); }

inline double delta_r2_eta(const FourVector &a, const FourVector &b) { return b.delta_r2_eta(a); }

inline double delta_r_eta(const FourVector &a, const FourVector &b) { return b.delta_r_eta(a); }

inline double delta_r2_rap(const FourVector &a, const FourVector &b) { return b.delta_r2_rap(a); }

inline double delta_r_rap(const FourVector &a, const FourVector &b) { return b.delta_r_rap(a); }

 
 
} // namespace HepMC3
#endif