#ifndef FOURVECTORS_HEADER_H
#define FOURVECTORS_HEADER_H
#include<math.h>
#include <iostream.h>

class fourvector {
      public:
        static const float me;
        static const float mmu;
        static const float mK;
        static const float mp;
        static const float mN;
        static const float mpi;
        static const float mpi0;
      protected:
        float E, px, py, pz, m;
      private:
        inline float mass() {
          float M  = E*E - px*px - py*py - pz*pz;
          return (M > 0) ? sqrt (M) : -sqrt (-M);
        };
      public:
        fourvector () {
          px = 0;
          py = 0;
          pz = 0;
          E  = 0;
          m  = 0;
        };
        fourvector (float Px, float Py, float Pz, float En) {
          px = Px;
          py = Py;
          pz = Pz;
          E  = En;
          m  = mass();
        };
        inline fourvector (float Px, float Py, float Pz, float, float M) {
          px = Px;
          py = Py;
          pz = Pz;
          m  = M;
          E  = sqrt (px*px + py*py + pz*pz + m*fabs (m));
        };
        fourvector (float Px, float Py, float Pz) {
          px = Px;
          py = Py;
          pz = Pz;
          m  = mpi;
          E  = sqrt (px*px + py*py + pz*pz + m*fabs (m));
        };
        inline fourvector& operator += (const fourvector& v) {
          px += v.px;
          py += v.py;
          pz += v.pz;
          E  += v.E;
          m  = mass();
          return *this;
        };
        inline fourvector& operator -= (const fourvector& v) {
          px -= v.px;
          py -= v.py;
          pz -= v.pz;
          E  -= v.E;
          m  = mass();
          return *this;
        };
        friend fourvector operator + (const fourvector&, const fourvector&); 
        friend fourvector operator - (const fourvector&, const fourvector&); 
        inline float getpx() {return px;};
        inline float getpy() {return py;};
        inline float getpz() {return pz;};
        inline float getpt() {return sqrt (px*px + py*py);};
        inline float getp () {return sqrt (px*px + py*py + pz*pz);};
        inline float getE () {return E;};
        inline float getm () {return m;};
        inline float getphi () {return atan2 (py, px);};
        inline float gettheta () {
          float pt2 = px*px + py*py;
          return atan2 (sqrt (pt2), pz);
        };
        inline float geteta () {
          float p = getp();
          return log ((p + pz)/(p - pz));
        };
        inline float getrapidity () {
          return log ((E + pz)/(E - pz));
        };
};

        inline fourvector operator + (const fourvector& v1, const fourvector& v2) {
          fourvector v = v1;
          v += v2;
          return v;
        }; 
        inline fourvector operator - (const fourvector& v1, const fourvector& v2) {
          fourvector v = v1;
          v -= v2;
          return v;
        }; 
#endif