// Functions for four-vector class
// Glen Cowan, RHUL Physics Dept., 4 December 2001.

#include <iostream.h>
#include <math.h>
#include "ThreeVector.h"
#include "FourVector.h"

FourVector::FourVector (float x, float y, float z, float t){
  xVal = x;
  yVal = y;
  zVal = z;
  tVal = t;
}

FourVector::FourVector(){
  xVal = 0.;
  yVal = 0.;
  zVal = 0.;
  tVal = 0.;
}

float FourVector::x(){ return xVal; }
float FourVector::y(){ return yVal; }
float FourVector::z(){ return zVal; }
float FourVector::t() { return tVal; }

float FourVector::m(){
  float p2 = pow(xVal,2) + pow(yVal,2) + pow(zVal,2);
  float E2 = pow(tVal,2);
  float mVal;
  if ( E2 > p2 ){
    mVal = sqrt (E2 - p2); }
  else if ( E2 == p2 ) {
    mVal = 0.; }
  else {
    //    cout << "FourVector::m() >> negative mass set to zero" << endl;
    mVal = 0.;
  }
  return mVal;
}

void FourVector::setX(float x) { xVal = x; }
void FourVector::setY(float y) { yVal = y; }
void FourVector::setZ(float z) { zVal = z; }
void FourVector::setT(float t) { tVal = t; }

ThreeVector FourVector::betaGamma(){
  float bgx, bgy, bgz;
  float mass = this->m();
  if ( mass > 0 ) {
    bgx = this->x() / mass ;
    bgy = this->y() / mass ;  
    bgz = this->z() / mass ;
  }
  else {
    cout << "FourVector::betaGamma() >> beta_gamma set to zero" << endl;
    bgx = 0.;
    bgy = 0.;
    bgz = 0.;
  }
  ThreeVector bg(bgx, bgy, bgz);
  return bg; 
}

ThreeVector FourVector::space(){
  ThreeVector A(xVal, yVal, zVal);
  return A;
}

FourVector FourVector::operator+ (FourVector v){
  float x = xVal + v.x();
  float y = yVal + v.y();
  float z = zVal + v.z();
  float t = tVal + v.t();
  FourVector A(x, y, z, t);
  return A;
}

FourVector FourVector::operator- (FourVector v){
  float x = xVal - v.x();
  float y = yVal - v.y();
  float z = zVal - v.z();
  float t = tVal - v.t();
  FourVector A(x, y, z, t);
  return A;
}

float FourVector::operator* (FourVector v){
  float scalar_prod = pow((tVal+v.t()),2) - pow((xVal+v.x()),2) - 
                      pow((yVal+v.y()),2) - pow((zVal+v.z()),2);
  return scalar_prod;
}

// Related non-member functions

FourVector boostVector(FourVector A, ThreeVector betaGamma){

  float bg = betaGamma.r();
  float beta = bg / sqrt ( 1. + pow(bg,2) );
  float gamma;
  if ( beta > 0. ) {
    gamma = bg / beta;
  }
  else {
    gamma = 1.;
  }

  float bg_dot_A = dotProd (betaGamma, A.space());
  float Ax_prime = A.x() - betaGamma.x() * (A.t() - (bg_dot_A/(1.+gamma)));
  float Ay_prime = A.y() - betaGamma.y() * (A.t() - (bg_dot_A/(1.+gamma)));
  float Az_prime = A.z() - betaGamma.z() * (A.t() - (bg_dot_A/(1.+gamma)));
  float At_prime = gamma*A.t()  - bg_dot_A;

  FourVector boostedVector (Ax_prime, Ay_prime, Az_prime, At_prime);
  return boostedVector;

}