Vector4.hh

#ifndef RIVET_MATH_VECTOR4
#define RIVET_MATH_VECTOR4
 
#include "Rivet/Math/MathConstants.hh"
#include "Rivet/Math/MathUtils.hh"
#include "Rivet/Math/Vector3.hh"
#include "Rivet/Math/VectorN.hh"
#include "Rivet/Tools/TypeTraits.hh"
 
// Forward declaration
namespace fastjet {
  class PseudoJet;
}
 
namespace Rivet {
 
 
  class FourVector;
  typedef FourVector Vector4;
  typedef FourVector V4;
 
  class FourMomentum;
  typedef FourMomentum P4;
 
  class LorentzTransform;
  FourVector transform(const LorentzTransform& lt, const FourVector& v4);
 

  class FourVector : public Vector<4> {
    friend FourVector multiply(const double a, const FourVector& v);
    friend FourVector multiply(const FourVector& v, const double a);
    friend FourVector add(const FourVector& a, const FourVector& b);
    friend FourVector transform(const LorentzTransform& lt, const FourVector& v4);
 
  public:
 
    FourVector()
        : Vector<4>() { }
 
    template <typename V4TYPE, typename std::enable_if<HasXYZT<V4TYPE>::value, int>::type DUMMY = 0>
    FourVector(const V4TYPE& other) {
      this->setT(other.t());
      this->setX(other.x());
      this->setY(other.y());
      this->setZ(other.z());
    }
 
    FourVector(const Vector<4>& other)
        : Vector<4>(other) { }
 
    FourVector(const double t, const double x, const double y, const double z) {
      this->setT(t);
      this->setX(x);
      this->setY(y);
      this->setZ(z);
    }
 
    FourVector(double t, const Vector3& v3) {
      this->setT(t);
      this->setX(v3.x());
      this->setY(v3.y());
      this->setZ(v3.z());
    }
 
    virtual ~FourVector() { }

    operator fastjet::PseudoJet() const;
 
 
  public:
 
    double t() const {
      return get(0);
    }
    double t2() const {
      return sqr(t());
    }
    FourVector& setT(const double t) {
      set(0, t);
      return *this;
    }
 
    double x() const {
      return get(1);
    }
    double x2() const {
      return sqr(x());
    }
    FourVector& setX(const double x) {
      set(1, x);
      return *this;
    }
 
    double y() const {
      return get(2);
    }
    double y2() const {
      return sqr(y());
    }
    FourVector& setY(const double y) {
      set(2, y);
      return *this;
    }
 
    double z() const {
      return get(3);
    }
    double z2() const {
      return sqr(z());
    }
    FourVector& setZ(const double z) {
      set(3, z);
      return *this;
    }
 
    double invariant() const {
      // Done this way for numerical precision
      return (t() + z()) * (t() - z()) - x() * x() - y() * y();
    }
 
    bool isNull() const {
      return Rivet::isZero(invariant());
    }

    double angle(const FourVector& v) const {
      return vector3().angle(v.vector3());
    }

    double angle(const Vector3& v3) const {
      return vector3().angle(v3);
    }

    double polarRadius2() const {
      return vector3().polarRadius2();
    }

    double perp2() const {
      return vector3().perp2();
    }

    double rho2() const {
      return vector3().rho2();
    }

    double polarRadius() const {
      return vector3().polarRadius();
    }

    double perp() const {
      return vector3().perp();
    }

    double rho() const {
      return vector3().rho();
    }

    Vector3 polarVec() const {
      return vector3().polarVec();
    }

    Vector3 perpVec() const {
      return vector3().perpVec();
    }

    Vector3 rhoVec() const {
      return vector3().rhoVec();
    }

    double azimuthalAngle(const PhiMapping mapping = ZERO_2PI) const {
      return vector3().azimuthalAngle(mapping);
    }

    double phi(const PhiMapping mapping = ZERO_2PI) const {
      return vector3().phi(mapping);
    }

    double polarAngle() const {
      return vector3().polarAngle();
    }

    double theta() const {
      return vector3().theta();
    }

    double pseudorapidity() const {
      return vector3().pseudorapidity();
    }

    double eta() const {
      return vector3().eta();
    }

    double abspseudorapidity() const {
      return fabs(eta());
    }

    double abseta() const {
      return fabs(eta());
    }

    Vector3 vector3() const {
      return Vector3(get(1), get(2), get(3));
    }

    Vector3 v3() const {
      return vector3();
    }

    operator Vector3() const {
      return vector3();
    }
 
 
  public:

    double contract(const FourVector& v) const {
      const double result = t() * v.t() - x() * v.x() - y() * v.y() - z() * v.z();
      return result;
    }

    double dot(const FourVector& v) const {
      return contract(v);
    }

    double operator*(const FourVector& v) const {
      return contract(v);
    }

    FourVector& operator*=(double a) {
      _vec = multiply(a, *this)._vec;
      return *this;
    }

    FourVector& operator/=(double a) {
      _vec = multiply(1.0 / a, *this)._vec;
      return *this;
    }

    FourVector& operator+=(const FourVector& v) {
      _vec = add(*this, v)._vec;
      return *this;
    }

    FourVector& operator-=(const FourVector& v) {
      _vec = add(*this, -v)._vec;
      return *this;
    }

    FourVector operator-() const {
      FourVector result;
      result._vec = -_vec;
      return result;
    }

    FourVector reverse() const {
      FourVector result = -*this;
      result.setT(-result.t());
      return result;
    }
  };
 

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

  inline double dot(const FourVector& a, const FourVector& b) {
    return contract(a, b);
  }
 
  inline FourVector multiply(const double a, const FourVector& v) {
    FourVector result;
    result._vec = a * v._vec;
    return result;
  }
 
  inline FourVector multiply(const FourVector& v, const double a) {
    return multiply(a, v);
  }
 
  inline FourVector operator*(const double a, const FourVector& v) {
    return multiply(a, v);
  }
 
  inline FourVector operator*(const FourVector& v, const double a) {
    return multiply(a, v);
  }
 
  inline FourVector operator/(const FourVector& v, const double a) {
    return multiply(1.0 / a, v);
  }
 
  inline FourVector add(const FourVector& a, const FourVector& b) {
    FourVector result;
    result._vec = a._vec + b._vec;
    return result;
  }
 
  inline FourVector operator+(const FourVector& a, const FourVector& b) {
    return add(a, b);
  }
 
  inline FourVector operator-(const FourVector& a, const FourVector& b) {
    return add(a, -b);
  }

  inline double invariant(const FourVector& lv) {
    return lv.invariant();
  }

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

  inline double angle(const Vector3& a, const FourVector& b) {
    return angle(a, b.vector3());
  }

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

 

  class FourMomentum : public FourVector {
    friend FourMomentum multiply(const double a, const FourMomentum& v);
    friend FourMomentum multiply(const FourMomentum& v, const double a);
    friend FourMomentum add(const FourMomentum& a, const FourMomentum& b);
    friend FourMomentum transform(const LorentzTransform& lt, const FourMomentum& v4);
 
  public:
 
    FourMomentum() { }
 
    template <typename V4TYPE, typename std::enable_if<HasXYZT<V4TYPE>::value, int>::type DUMMY = 0>
    FourMomentum(const V4TYPE& other) {
      this->setE(other.t());
      this->setPx(other.x());
      this->setPy(other.y());
      this->setPz(other.z());
    }
 
    FourMomentum(const Vector<4>& other)
        : FourVector(other) { }
 
    FourMomentum(const double E, const double px, const double py, const double pz) {
      this->setE(E);
      this->setPx(px);
      this->setPy(py);
      this->setPz(pz);
    }
 
    FourMomentum(const Vector3& v3, double m) {
      const double e = sqrt(v3.mod2() + m * m);
      this->setE(e);
      this->setPx(v3.x());
      this->setPy(v3.y());
      this->setPz(v3.z());
    }
 
    ~FourMomentum() { }
 
  public:
 


    FourMomentum& setE(double E) {
      setT(E);
      return *this;
    }

    FourMomentum& setPx(double px) {
      setX(px);
      return *this;
    }

    FourMomentum& setPy(double py) {
      setY(py);
      return *this;
    }

    FourMomentum& setPz(double pz) {
      setZ(pz);
      return *this;
    }
 

    FourMomentum& setPE(double px, double py, double pz, double E) {
      if (E < 0) throw std::invalid_argument("Negative energy given as argument: " + to_str(E));
      setPx(px);
      setPy(py);
      setPz(pz);
      setE(E);
      return *this;
    }

    FourMomentum& setXYZE(double px, double py, double pz, double E) {
      return setPE(px, py, pz, E);
    }
    // /// Near-alias with switched arg order
    // FourMomentum& setEP(double E, double px, double py, double pz) {
    //   return setPE(px, py, pz, E);
    // }
    // /// Alias for setEP
    // FourMomentum& setEXYZ(double E, double px, double py, double pz) {
    //   return setEP(E, px, py, pz);
    // }
 

    FourMomentum& setPM(double px, double py, double pz, double mass) {
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument: " + to_str(mass));
      const double E = sqrt(sqr(mass) + sqr(px) + sqr(py) + sqr(pz));
      // setPx(px); setPy(py); setPz(pz); setE(E);
      return setPE(px, py, pz, E);
    }

    FourMomentum& setXYZM(double px, double py, double pz, double mass) {
      return setPM(px, py, pz, mass);
    }
 

    FourMomentum& setEtaPhiME(double eta, double phi, double mass, double E) {
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument");
      if (E < 0) throw std::invalid_argument("Negative energy given as argument");
      const double theta = 2 * atan(exp(-eta));
      if (theta < 0 || theta > M_PI) throw std::domain_error("Polar angle outside 0..pi in calculation");
      setThetaPhiME(theta, phi, mass, E);
      return *this;
    }

    FourMomentum& setEtaPhiMPt(double eta, double phi, double mass, double pt) {
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument");
      if (pt < 0) throw std::invalid_argument("Negative transverse momentum given as argument");
      const double theta = 2 * atan(exp(-eta));
      if (theta < 0 || theta > M_PI) throw std::domain_error("Polar angle outside 0..pi in calculation");
      const double p = pt / sin(theta);
      const double E = sqrt(sqr(p) + sqr(mass));
      setThetaPhiME(theta, phi, mass, E);
      return *this;
    }

    FourMomentum& setRapPhiME(double y, double phi, double mass, double E) {
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument");
      if (E < 0) throw std::invalid_argument("Negative energy given as argument");
      const double sqrt_pt2_m2 = E / cosh(y);
      const double pt = sqrt(sqr(sqrt_pt2_m2) - sqr(mass));
      if (pt < 0) throw std::domain_error("Negative transverse momentum in calculation");
      const double pz = sqrt_pt2_m2 * sinh(y);
      const double px = pt * cos(phi);
      const double py = pt * sin(phi);
      setPE(px, py, pz, E);
      return *this;
    }

    FourMomentum& setRapPhiMPt(double y, double phi, double mass, double pt) {
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument");
      if (pt < 0) throw std::invalid_argument("Negative transverse mass given as argument");
      const double E = sqrt(sqr(pt) + sqr(mass)) * cosh(y);
      if (E < 0) throw std::domain_error("Negative energy in calculation");
      setRapPhiME(y, phi, mass, E);
      return *this;
    }

    FourMomentum& setThetaPhiME(double theta, double phi, double mass, double E) {
      if (theta < 0 || theta > M_PI)
        throw std::invalid_argument("Polar angle outside 0..pi given as argument");
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument");
      if (E < 0) throw std::invalid_argument("Negative energy given as argument");
      const double p = sqrt(sqr(E) - sqr(mass));
      const double pz = p * cos(theta);
      const double pt = p * sin(theta);
      if (pt < 0) throw std::invalid_argument("Negative transverse momentum in calculation");
      const double px = pt * cos(phi);
      const double py = pt * sin(phi);
      setPE(px, py, pz, E);
      return *this;
    }

    FourMomentum& setThetaPhiMPt(double theta, double phi, double mass, double pt) {
      if (theta < 0 || theta > M_PI)
        throw std::invalid_argument("Polar angle outside 0..pi given as argument");
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument");
      if (pt < 0) throw std::invalid_argument("Negative transverse momentum given as argument");
      const double p = pt / sin(theta);
      const double px = pt * cos(phi);
      const double py = pt * sin(phi);
      const double pz = p * cos(theta);
      const double E = sqrt(sqr(p) + sqr(mass));
      setPE(px, py, pz, E);
      return *this;
    }

    FourMomentum& setPtPhiME(double pt, double phi, double mass, double E) {
      if (pt < 0) throw std::invalid_argument("Negative transverse momentum given as argument");
      if (mass < 0) throw std::invalid_argument("Negative mass given as argument");
      if (E < 0) throw std::invalid_argument("Negative energy given as argument");
      const double px = pt * cos(phi);
      const double py = pt * sin(phi);
      const double pz = sqrt(sqr(E) - sqr(mass) - sqr(pt));
      setPE(px, py, pz, E);
      return *this;
    }

 


    double E() const {
      return t();
    }

    double E2() const {
      return t2();
    }

    double px() const {
      return x();
    }

    double px2() const {
      return x2();
    }

    double py() const {
      return y();
    }

    double py2() const {
      return y2();
    }

    double pz() const {
      return z();
    }

    double pz2() const {
      return z2();
    }
 

    double mass() const {
      // assert(Rivet::isZero(mass2()) || mass2() > 0);
      // if (Rivet::isZero(mass2())) {
      //   return 0.0;
      // } else {
      //   return sqrt(mass2());
      // }
      return sign(mass2()) * sqrt(fabs(mass2()));
    }

    double mass2() const {
      return invariant();
    }
 

    Vector3 p3() const {
      return vector3();
    }

    double p() const {
      return p3().mod();
    }

    double p2() const {
      return p3().mod2();
    }
 

    double rapidity() const {
      if (E() == 0.0) return 0.0;                             
      if (E() == fabs(pz())) return std::copysign(INF, pz()); 
      return 0.5 * std::log((E() + pz()) / (E() - pz()));
    }

    double rap() const {
      return rapidity();
    }

    double absrapidity() const {
      return fabs(rapidity());
    }

    double absrap() const {
      return fabs(rap());
    }

    Vector3 pTvec() const {
      return p3().polarVec();
    }

    Vector3 ptvec() const {
      return pTvec();
    }

    double pT2() const {
      return vector3().polarRadius2();
    }

    double pt2() const {
      return vector3().polarRadius2();
    }

    double pT() const {
      return sqrt(pT2());
    }

    double pt() const {
      return sqrt(pT2());
    }

    double Et2() const {
      return Et() * Et();
    }

    double Et() const {
      return E() * sin(polarAngle());
    }

 


    double gamma() const {
      return sqrt(E2() / mass2());
    }

    Vector3 gammaVec() const {
      return gamma() * p3().unit();
    }

    double beta() const {
      return p() / E();
    }

    Vector3 betaVec() const {
      // return Vector3(px()/E(), py()/E(), pz()/E());
      return p3() / E();
    }

 



    FourMomentum& operator*=(double a) {
      _vec = multiply(a, *this)._vec;
      return *this;
    }

    FourMomentum& operator/=(double a) {
      _vec = multiply(1.0 / a, *this)._vec;
      return *this;
    }

    FourMomentum& operator+=(const FourMomentum& v) {
      _vec = add(*this, v)._vec;
      return *this;
    }

    FourMomentum& operator-=(const FourMomentum& v) {
      _vec = add(*this, -v)._vec;
      return *this;
    }

    FourMomentum operator-() const {
      FourMomentum result;
      result._vec = -_vec;
      return result;
    }

    FourMomentum reverse() const {
      FourMomentum result = -*this;
      result.setE(-result.E());
      return result;
    }

 

 


    static FourMomentum mkXYZE(double px, double py, double pz, double E) {
      return FourMomentum().setPE(px, py, pz, E);
    }

    static FourMomentum mkXYZM(double px, double py, double pz, double mass) {
      return FourMomentum().setPM(px, py, pz, mass);
    }

    static FourMomentum mkEtaPhiME(double eta, double phi, double mass, double E) {
      return FourMomentum().setEtaPhiME(eta, phi, mass, E);
    }

    static FourMomentum mkEtaPhiMPt(double eta, double phi, double mass, double pt) {
      return FourMomentum().setEtaPhiMPt(eta, phi, mass, pt);
    }

    static FourMomentum mkRapPhiME(double y, double phi, double mass, double E) {
      return FourMomentum().setRapPhiME(y, phi, mass, E);
    }

    static FourMomentum mkRapPhiMPt(double y, double phi, double mass, double pt) {
      return FourMomentum().setRapPhiMPt(y, phi, mass, pt);
    }

    static FourMomentum mkThetaPhiME(double theta, double phi, double mass, double E) {
      return FourMomentum().setThetaPhiME(theta, phi, mass, E);
    }

    static FourMomentum mkThetaPhiMPt(double theta, double phi, double mass, double pt) {
      return FourMomentum().setThetaPhiMPt(theta, phi, mass, pt);
    }

    static FourMomentum mkPtPhiME(double pt, double phi, double mass, double E) {
      return FourMomentum().setPtPhiME(pt, phi, mass, E);
    }

  };
 
 
  inline FourMomentum multiply(const double a, const FourMomentum& v) {
    FourMomentum result;
    result._vec = a * v._vec;
    return result;
  }
 
  inline FourMomentum multiply(const FourMomentum& v, const double a) {
    return multiply(a, v);
  }
 
  inline FourMomentum operator*(const double a, const FourMomentum& v) {
    return multiply(a, v);
  }
 
  inline FourMomentum operator*(const FourMomentum& v, const double a) {
    return multiply(a, v);
  }
 
  inline FourMomentum operator/(const FourMomentum& v, const double a) {
    return multiply(1.0 / a, v);
  }
 
  inline FourMomentum add(const FourMomentum& a, const FourMomentum& b) {
    FourMomentum result;
    result._vec = a._vec + b._vec;
    return result;
  }
 
  inline FourMomentum operator+(const FourMomentum& a, const FourMomentum& b) {
    return add(a, b);
  }
 
  inline FourMomentum operator-(const FourMomentum& a, const FourMomentum& b) {
    return add(a, -b);
  }
 

 


  inline double deltaR2(const FourVector& a, const FourVector& b, RapScheme scheme = PSEUDORAPIDITY) {
    switch (scheme) {
      case PSEUDORAPIDITY: return deltaR2(a.vector3(), b.vector3());
      case RAPIDITY: {
        const FourMomentum* ma = dynamic_cast<const FourMomentum*>(&a);
        const FourMomentum* mb = dynamic_cast<const FourMomentum*>(&b);
        if (!ma || !mb) {
          string err =
              "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
          throw std::runtime_error(err);
        }
        return deltaR2(*ma, *mb, scheme);
      }
      default: throw std::runtime_error("The specified deltaR scheme is not yet implemented");
    }
  }

  inline double deltaR(const FourVector& a, const FourVector& b, RapScheme scheme = PSEUDORAPIDITY) {
    return sqrt(deltaR2(a, b, scheme));
  }
 

  inline double deltaR2(const FourVector& v, double eta2, double phi2, RapScheme scheme = PSEUDORAPIDITY) {
    switch (scheme) {
      case PSEUDORAPIDITY: return deltaR2(v.vector3(), eta2, phi2);
      case RAPIDITY: {
        const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
        if (!mv) {
          string err =
              "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
          throw std::runtime_error(err);
        }
        return deltaR2(*mv, eta2, phi2, scheme);
      }
      default: throw std::runtime_error("The specified deltaR scheme is not yet implemented");
    }
  }

  inline double deltaR(const FourVector& v, double eta2, double phi2, RapScheme scheme = PSEUDORAPIDITY) {
    return sqrt(deltaR2(v, eta2, phi2, scheme));
  }
 

  inline double deltaR2(double eta1, double phi1, const FourVector& v, RapScheme scheme = PSEUDORAPIDITY) {
    switch (scheme) {
      case PSEUDORAPIDITY: return deltaR2(eta1, phi1, v.vector3());
      case RAPIDITY: {
        const FourMomentum* mv = dynamic_cast<const FourMomentum*>(&v);
        if (!mv) {
          string err =
              "deltaR with scheme RAPIDITY can only be called with FourMomentum objects, not FourVectors";
          throw std::runtime_error(err);
        }
        return deltaR2(eta1, phi1, *mv, scheme);
      }
      default: throw std::runtime_error("The specified deltaR scheme is not yet implemented");
    }
  }

  inline double deltaR(double eta1, double phi1, const FourVector& v, RapScheme scheme = PSEUDORAPIDITY) {
    return sqrt(deltaR2(eta1, phi1, v, scheme));
  }
 

  inline double deltaR2(const FourMomentum& a, const FourMomentum& b, RapScheme scheme = PSEUDORAPIDITY) {
    switch (scheme) {
      case PSEUDORAPIDITY: return deltaR2(a.vector3(), b.vector3());
      case RAPIDITY: return deltaR2(a.rapidity(), a.azimuthalAngle(), b.rapidity(), b.azimuthalAngle());
      default: throw std::runtime_error("The specified deltaR scheme is not yet implemented");
    }
  }

  inline double deltaR(const FourMomentum& a, const FourMomentum& b, RapScheme scheme = PSEUDORAPIDITY) {
    return sqrt(deltaR2(a, b, scheme));
  }
 

  inline double deltaR2(const FourMomentum& v, double eta2, double phi2, RapScheme scheme = PSEUDORAPIDITY) {
    switch (scheme) {
      case PSEUDORAPIDITY: return deltaR2(v.vector3(), eta2, phi2);
      case RAPIDITY: return deltaR2(v.rapidity(), v.azimuthalAngle(), eta2, phi2);
      default: throw std::runtime_error("The specified deltaR scheme is not yet implemented");
    }
  }

  inline double deltaR(const FourMomentum& v, double eta2, double phi2, RapScheme scheme = PSEUDORAPIDITY) {
    return sqrt(deltaR2(v, eta2, phi2, scheme));
  }
 

  inline double deltaR2(double eta1, double phi1, const FourMomentum& v, RapScheme scheme = PSEUDORAPIDITY) {
    switch (scheme) {
      case PSEUDORAPIDITY: return deltaR2(eta1, phi1, v.vector3());
      case RAPIDITY: return deltaR2(eta1, phi1, v.rapidity(), v.azimuthalAngle());
      default: throw std::runtime_error("The specified deltaR scheme is not yet implemented");
    }
  }

  inline double deltaR(double eta1, double phi1, const FourMomentum& v, RapScheme scheme = PSEUDORAPIDITY) {
    return sqrt(deltaR2(eta1, phi1, v, scheme));
  }
 

  inline double deltaR2(const FourMomentum& a, const FourVector& b, RapScheme scheme = PSEUDORAPIDITY) {
    switch (scheme) {
      case PSEUDORAPIDITY: return deltaR2(a.vector3(), b.vector3());
      case RAPIDITY:
        return deltaR2(a.rapidity(), a.azimuthalAngle(), FourMomentum(b).rapidity(), b.azimuthalAngle());
      default: throw std::runtime_error("The specified deltaR scheme is not yet implemented");
    }
  }

  inline double deltaR(const FourMomentum& a, const FourVector& b, RapScheme scheme = PSEUDORAPIDITY) {
    return sqrt(deltaR2(a, b, scheme));
  }
 

  inline double deltaR2(const FourVector& a, const FourMomentum& b, RapScheme scheme = PSEUDORAPIDITY) {
    return deltaR2(b, a, scheme); //< note reversed args
  }

  inline double deltaR(const FourVector& a, const FourMomentum& b, RapScheme scheme = PSEUDORAPIDITY) {
    return deltaR(b, a, scheme); //< note reversed args
  }
 

  inline double deltaR2(const FourMomentum& a, const Vector3& b) {
    return deltaR2(a.vector3(), b);
  }

  inline double deltaR(const FourMomentum& a, const Vector3& b) {
    return deltaR(a.vector3(), b);
  }

  inline double deltaR2(const Vector3& a, const FourMomentum& b) {
    return deltaR2(a, b.vector3());
  }

  inline double deltaR(const Vector3& a, const FourMomentum& b) {
    return deltaR(a, b.vector3());
  }

  inline double deltaR2(const FourVector& a, const Vector3& b) {
    return deltaR2(a.vector3(), b);
  }

  inline double deltaR(const FourVector& a, const Vector3& b) {
    return deltaR(a.vector3(), b);
  }

  inline double deltaR2(const Vector3& a, const FourVector& b) {
    return deltaR2(a, b.vector3());
  }

  inline double deltaR(const Vector3& a, const FourVector& b) {
    return deltaR(a, b.vector3());
  }

 

 


  inline double deltaPhi(const FourMomentum& a, const FourMomentum& b, bool sign = false) {
    return deltaPhi(a.vector3(), b.vector3(), sign);
  }

  inline double deltaPhi(const FourMomentum& v, double phi2, bool sign = false) {
    return deltaPhi(v.vector3(), phi2, sign);
  }

  inline double deltaPhi(double phi1, const FourMomentum& v, bool sign = false) {
    return deltaPhi(phi1, v.vector3(), sign);
  }

  inline double deltaPhi(const FourVector& a, const FourVector& b, bool sign = false) {
    return deltaPhi(a.vector3(), b.vector3(), sign);
  }

  inline double deltaPhi(const FourVector& v, double phi2, bool sign = false) {
    return deltaPhi(v.vector3(), phi2, sign);
  }

  inline double deltaPhi(double phi1, const FourVector& v, bool sign = false) {
    return deltaPhi(phi1, v.vector3(), sign);
  }

  inline double deltaPhi(const FourVector& a, const FourMomentum& b, bool sign = false) {
    return deltaPhi(a.vector3(), b.vector3(), sign);
  }

  inline double deltaPhi(const FourMomentum& a, const FourVector& b, bool sign = false) {
    return deltaPhi(a.vector3(), b.vector3(), sign);
  }

  inline double deltaPhi(const FourVector& a, const Vector3& b, bool sign = false) {
    return deltaPhi(a.vector3(), b, sign);
  }

  inline double deltaPhi(const Vector3& a, const FourVector& b, bool sign = false) {
    return deltaPhi(a, b.vector3(), sign);
  }

  inline double deltaPhi(const FourMomentum& a, const Vector3& b, bool sign = false) {
    return deltaPhi(a.vector3(), b, sign);
  }

  inline double deltaPhi(const Vector3& a, const FourMomentum& b, bool sign = false) {
    return deltaPhi(a, b.vector3(), sign);
  }

 

 


  inline double deltaEta(const FourMomentum& a, const FourMomentum& b, bool sign = false) {
    return deltaEta(a.vector3(), b.vector3(), sign);
  }

  inline double deltaEta(const FourMomentum& v, double eta2, bool sign = false) {
    return deltaEta(v.vector3(), eta2, sign);
  }

  inline double deltaEta(double eta1, const FourMomentum& v, bool sign = false) {
    return deltaEta(eta1, v.vector3(), sign);
  }

  inline double deltaEta(const FourVector& a, const FourVector& b, bool sign = false) {
    return deltaEta(a.vector3(), b.vector3(), sign);
  }

  inline double deltaEta(const FourVector& v, double eta2, bool sign = false) {
    return deltaEta(v.vector3(), eta2, sign);
  }

  inline double deltaEta(double eta1, const FourVector& v, bool sign = false) {
    return deltaEta(eta1, v.vector3(), sign);
  }

  inline double deltaEta(const FourVector& a, const FourMomentum& b, bool sign = false) {
    return deltaEta(a.vector3(), b.vector3(), sign);
  }

  inline double deltaEta(const FourMomentum& a, const FourVector& b, bool sign = false) {
    return deltaEta(a.vector3(), b.vector3(), sign);
  }

  inline double deltaEta(const FourVector& a, const Vector3& b, bool sign = false) {
    return deltaEta(a.vector3(), b, sign);
  }

  inline double deltaEta(const Vector3& a, const FourVector& b, bool sign = false) {
    return deltaEta(a, b.vector3(), sign);
  }

  inline double deltaEta(const FourMomentum& a, const Vector3& b, bool sign = false) {
    return deltaEta(a.vector3(), b, sign);
  }

  inline double deltaEta(const Vector3& a, const FourMomentum& b, bool sign = false) {
    return deltaEta(a, b.vector3(), sign);
  }

 


  inline double deltaRap(const FourMomentum& a, const FourMomentum& b, bool sign = false) {
    return deltaRap(a.rapidity(), b.rapidity(), sign);
  }

  inline double deltaRap(const FourMomentum& v, double y2, bool sign = false) {
    return deltaRap(v.rapidity(), y2, sign);
  }

  inline double deltaRap(double y1, const FourMomentum& v, bool sign = false) {
    return deltaRap(y1, v.rapidity(), sign);
  }

 

 



  inline bool cmpMomByPt(const FourMomentum& a, const FourMomentum& b) {
    return a.pt() > b.pt();
  }

  inline bool cmpMomByAscPt(const FourMomentum& a, const FourMomentum& b) {
    return a.pt() < b.pt();
  }

  inline bool cmpMomByP(const FourMomentum& a, const FourMomentum& b) {
    return a.vector3().mod() > b.vector3().mod();
  }

  inline bool cmpMomByAscP(const FourMomentum& a, const FourMomentum& b) {
    return a.vector3().mod() < b.vector3().mod();
  }

  inline bool cmpMomByEt(const FourMomentum& a, const FourMomentum& b) {
    return a.Et() > b.Et();
  }

  inline bool cmpMomByAscEt(const FourMomentum& a, const FourMomentum& b) {
    return a.Et() < b.Et();
  }

  inline bool cmpMomByE(const FourMomentum& a, const FourMomentum& b) {
    return a.E() > b.E();
  }

  inline bool cmpMomByAscE(const FourMomentum& a, const FourMomentum& b) {
    return a.E() < b.E();
  }

  inline bool cmpMomByMass(const FourMomentum& a, const FourMomentum& b) {
    return a.mass() > b.mass();
  }

  inline bool cmpMomByAscMass(const FourMomentum& a, const FourMomentum& b) {
    return a.mass() < b.mass();
  }

  inline bool cmpMomByEta(const FourMomentum& a, const FourMomentum& b) {
    return a.eta() < b.eta();
  }

  inline bool cmpMomByDescEta(const FourMomentum& a, const FourMomentum& b) {
    return a.pseudorapidity() > b.pseudorapidity();
  }

  inline bool cmpMomByAbsEta(const FourMomentum& a, const FourMomentum& b) {
    return fabs(a.eta()) < fabs(b.eta());
  }

  inline bool cmpMomByDescAbsEta(const FourMomentum& a, const FourMomentum& b) {
    return fabs(a.eta()) > fabs(b.eta());
  }

  inline bool cmpMomByRap(const FourMomentum& a, const FourMomentum& b) {
    return a.rapidity() < b.rapidity();
  }

  inline bool cmpMomByDescRap(const FourMomentum& a, const FourMomentum& b) {
    return a.rapidity() > b.rapidity();
  }

  inline bool cmpMomByAbsRap(const FourMomentum& a, const FourMomentum& b) {
    return fabs(a.rapidity()) < fabs(b.rapidity());
  }

  inline bool cmpMomByDescAbsRap(const FourMomentum& a, const FourMomentum& b) {
    return fabs(a.rapidity()) > fabs(b.rapidity());
  }

 

  template <typename MOMS, typename CMP>
  inline MOMS& isortBy(MOMS& pbs, const CMP& cmp) {
    std::sort(pbs.begin(), pbs.end(), cmp);
    return pbs;
  }

  template <typename MOMS, typename CMP>
  inline MOMS sortBy(const MOMS& pbs, const CMP& cmp) {
    MOMS rtn = pbs;
    std::sort(rtn.begin(), rtn.end(), cmp);
    return rtn;
  }

  template <typename MOMS>
  inline MOMS& isortByPt(MOMS& pbs) {
    return isortBy(pbs, cmpMomByPt);
  }

  template <typename MOMS>
  inline MOMS sortByPt(const MOMS& pbs) {
    return sortBy(pbs, cmpMomByPt);
  }

  template <typename MOMS>
  inline MOMS& isortByE(MOMS& pbs) {
    return isortBy(pbs, cmpMomByE);
  }

  template <typename MOMS>
  inline MOMS sortByE(const MOMS& pbs) {
    return sortBy(pbs, cmpMomByE);
  }

  template <typename MOMS>
  inline MOMS& isortByEt(MOMS& pbs) {
    return isortBy(pbs, cmpMomByEt);
  }

  template <typename MOMS>
  inline MOMS sortByEt(const MOMS& pbs) {
    return sortBy(pbs, cmpMomByEt);
  }

 


  inline double mass(const FourMomentum& a, const FourMomentum& b) {
    return (a + b).mass();
  }

  inline double mass2(const FourMomentum& a, const FourMomentum& b) {
    return (a + b).mass2();
  }

  inline double mT(const FourMomentum& vis, const FourMomentum& invis) {
    return mT(vis.p3(), invis.p3());
  }

  inline double mT(const FourMomentum& vis, const Vector3& invis) {
    return mT(vis.p3(), invis);
  }

  inline double mT(const Vector3& vis, const FourMomentum& invis) {
    return mT(vis, invis.p3());
  }

  inline double mT(const FourMomentum& vis, const Vector4& invis) {
    return mT(vis.p3(), invis.vector3());
  }

  inline double mT(const Vector4& vis, const FourMomentum& invis) {
    return mT(vis.vector3(), invis.p3());
  }

  inline double pT(const FourMomentum& vis, const FourMomentum& invis) {
    return pT(vis.p3(), invis.p3());
  }

  inline double pT(const FourMomentum& vis, const Vector3& invis) {
    return pT(vis.p3(), invis);
  }

  inline double pT(const Vector3& vis, const FourMomentum& invis) {
    return pT(vis, invis.p3());
  }

 

 


  inline std::string toString(const FourVector& lv) {
    std::ostringstream out;
    out << "(" << (fabs(lv.t()) < 1E-30 ? 0.0 : lv.t()) << "; " << (fabs(lv.x()) < 1E-30 ? 0.0 : lv.x())
        << ", " << (fabs(lv.y()) < 1E-30 ? 0.0 : lv.y()) << ", " << (fabs(lv.z()) < 1E-30 ? 0.0 : lv.z())
        << ")";
    return out.str();
  }

  inline std::ostream& operator<<(std::ostream& out, const FourVector& lv) {
    out << toString(lv);
    return out;
  }


  typedef std::vector<FourVector> FourVectors;
  typedef std::vector<FourMomentum> FourMomenta;

 
 
}
 
#endif