source: huonglan/MYSIMULATION/MyVector4.C @ 15

#include <iostream>
#include <cstdio>
#include <math.h>

#include "MyVector4.h"

//using namespace std;


MyVector4::MyVector4()
{
  MyVector3 v3(0.,0.,0.);
  this->mv3 = v3;
  this->mt = 0.;
}

MyVector4::MyVector4(MyVector3 v3, double t)
{
  this->mv3 = v3;
  this->mt = t;
}

MyVector3 MyVector4::vector3() const
{
  return this->mv3;
}

double MyVector4::perp() const
{
 return sqrt((this->mv3.X())*(this->mv3.X())+ (this->mv3.Y())*(this->mv3.Y()) + (this->mv3.Z())*(this->mv3.Z()));
}

double MyVector4::time() const
{
  return this->mt;
}

double MyVector4::E() const
{
  return this->mt;
}

double MyVector4::mass() const
{
  double mass = pow((this->mt),2)-pow((this->mv3.mag()),2);
  if (mass<0){
    //cout << "Tachyon!" << endl;
    return -sqrt(-mass);}
  else
    return sqrt(mass);
}

double MyVector4::phi() const
{
  double phi4 = atan2(this->mv3.Y(),this->mv3.X());
  if (phi4 < - M_PI && phi4 > M_PI){
    cout << " Incorrect value of phi ! " << endl;
    return 0.;
  }
  else
    return phi4;
}

double MyVector4::Pt() const
{
  return this->mv3.perp();
}

double MyVector4::eta() const
{
  return this->mv3.eta();
}

double MyVector4::deltaR(MyVector4 &vec4) const
{
  MyVector3 vec3 = this -> mv3;
  MyVector3 vec33 = vec4.vector3();
  double delta_R = vec3.deltaR(vec33);
  return delta_R;
}

void MyVector4::SetPtEtaPhiM(double Pt, double eta, double phi, double M)
{
  MyVector3 v3(Pt*cos(phi),Pt*sin(phi),Pt*sinh(eta));
  this->mv3 = v3;
  this->mt = sqrt(M*M+v3.mag()*v3.mag());  
}

//###### OPERATEURS ##############

MyVector4& MyVector4::operator=(MyVector4 &vec)
{
  MyVector3 vec3 = vec.vector3();
  this->mv3 = vec3;
  this->mt = vec.time();

  return *this;
}

//Add 10 June !!!
bool MyVector4::operator==(MyVector4 &vec)
{
  MyVector3 vec1 = this->mv3;
  MyVector3 vec2 = vec.mv3;
  double t1 = this->mt;
  double t2 = vec.mt;
  if ((vec1 == vec2) && (t1 == t2))
    return true;
  else
    return false;
}
//***************

MyVector4& MyVector4::operator+=(MyVector4 &vec)
{
  MyVector3 vec3 = vec.vector3();
  this->mv3 += vec3;
  this->mt += vec.time();

  return *this;
}

MyVector4 operator+ (MyVector4 &vec1, MyVector4 &vec2)
{
  MyVector3 vec31 = vec1.vector3();
  MyVector3 vec32 = vec2.vector3();
  MyVector3 vec3_final = vec31 + vec32;

  double t1 = vec1.time();
  double t2 = vec2.time();

  double t_final = t1 + t2;

  return MyVector4(vec3_final,t_final);
}

MyVector4 operator- (MyVector4 &vec1, MyVector4 &vec2)
{
  MyVector3 vec31 = vec1.vector3();
  MyVector3 vec32 = vec2.vector3();
  MyVector3 vec3_final = vec31 - vec32;

  double t1 = vec1.time();
  double t2 = vec2.time();

  double t_final = t1 - t2;

  return MyVector4(vec3_final,t_final);
}

//####### BOOST #########

void MyVector4::boost(MyVector3 beta_vect)
{
  double beta = beta_vect.mag();

  if (beta>1)
    {
      cout <<"Beta > 1 : Forbidden! Do nothing !!!" << endl;
    }
  
  double gamma = 1./sqrt(1. - beta*beta);
  MyVector3 unit = beta_vect*(double)(1./beta);
  MyVector3 v3 = this->mv3;
  double time = this->mt;

  double scala = (double)(v3*unit);

  MyVector3 v3_paral = unit*scala;
  MyVector3 v3_perpen = v3 - v3_paral;

  MyVector3 v3_perpen_boost = v3_perpen;
  MyVector3 v3_paral_boost = unit*((scala + beta*time)*gamma);

  MyVector3 v3_boost = v3_perpen_boost + v3_paral_boost;        
  double t_boost = (double)gamma*(time + beta*scala);

  this->mv3 = v3_boost;
  this->mt = t_boost;
}

//********************* Print operator ************************
std::ostream& operator<<(std::ostream& out, const MyVector4& vec) 
{
  out << vec.vector3()
     << " " << vec.time()
     << " " << vec.mass()
     << std::endl;

 return out;
}