Nc4Vector.cxx

Go to the documentation of this file.

 
#include "Nc4Vector.h"
#include "Riostream.h"
 
ClassImp(Nc4Vector); // Class implementation to enable ROOT I/O
 
Nc4Vector::Nc4Vector()
{
 
 fUser=0;
 SetZero();
 fScalar=1;
}

Nc4Vector::~Nc4Vector()
{
 
 if (fUser)
 {
  delete fUser;
  fUser=0;
 }
}

Nc4Vector::Nc4Vector(const Nc4Vector& v)
{
 
 fScalar=v.fScalar;
 fV2=v.fV2;
 fDv2=v.fDv2;
 fV0=v.fV0;
 fDv0=v.fDv0;
 fDresult=v.fDresult;
 fV=v.fV;
 fUser=0;
 SetUserData(v.fUser);
}

Nc4Vector& Nc4Vector::operator=(const Nc4Vector& q)
{
 
 fScalar=q.fScalar;
 fV2=q.fV2;
 fDv2=q.fDv2;
 fV0=q.fV0;
 fDv0=q.fDv0;
 fDresult=q.fDresult;
 fV=q.fV;
 SetUserData(q.fUser);
 
 return *this;
}

void Nc4Vector::Load(Nc4Vector& q)
{
 
 fScalar=q.fScalar;
 fV2=q.fV2;
 fDv2=q.fDv2;
 fV0=q.fV0;
 fDv0=q.fDv0;
 fDresult=q.fDresult;
 fV=q.fV;
 SetUserData(q.fUser);
}

void Nc4Vector::SetZero()
{
 
 fV2=0;
 fDv2=0;
 fV0=0;
 fDv0=0;
 fDresult=0;
 fV.SetZero();
 if (fUser)
 {
  delete fUser;
  fUser=0;
 }
}

void Nc4Vector::SetVector(Double_t v0,Nc3Vector& v)
{
 
 fScalar=1;
 fV0=v0;
 fV=v;
 fV2=pow(v0,2)-fV.Dot(fV);
 SetScalarError(0);
}

void Nc4Vector::SetVector(Double_t* v,TString f,TString u)
{
 
 fScalar=1;
 Double_t a[3];
 for (Int_t i=0; i<3; i++)
 {
  a[i]=v[i+1];
 }
 fV0=v[0];
 fV.SetVector(a,f,u);
 fV2=pow(fV0,2)-fV.Dot(fV);
 fDv2=0;
 fDv0=0;
 fDresult=0;
}

void Nc4Vector::SetVector(Double_t v0,Double_t v1,Double_t v2,Double_t v3,TString f,TString u)
{
 
 Double_t vec[4];
 vec[0]=v0;
 vec[1]=v1;
 vec[2]=v2;
 vec[3]=v3;
 SetVector(vec,f,u);
}

void Nc4Vector::GetVector(Double_t* v,TString f,TString u)
{
 
 if (fScalar)
 {
  v[0]=fV0;
 }
 else
 {
  v[0]=sqrt(fV.Dot(fV)+fV2);
 } 
 Double_t a[3];
 fV.GetVector(a,f,u);
 for (Int_t i=0; i<3; i++)
 {
  v[i+1]=a[i];
 }
}

void Nc4Vector::SetVector(Float_t* v,TString f,TString u)
{
 
 Double_t vec[4];
 for (Int_t i=0; i<4; i++)
 {
  vec[i]=v[i];
 }
 SetVector(vec,f,u);
}

void Nc4Vector::GetVector(Float_t* v,TString f,TString u)
{
 
 Double_t vec[4];
 GetVector(vec,f,u);
 for (Int_t i=0; i<4; i++)
 {
  v[i]=vec[i];
 }
}

Double_t Nc4Vector::GetScalar()
{
 
 if (fScalar)
 {
  fDresult=fDv0;
  return fV0;
 }
 else
 {
  Double_t dot=fV.Dot(fV);
  Double_t ddot=fV.GetResultError();
  Double_t v02=dot+fV2;
  Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2));
  Double_t v0=sqrt(fabs(v02));
  Double_t dv0=0;
  if (v0) dv0=dv02/(2.*v0);
  fDresult=dv0;
  return v0;
 }
}

Int_t Nc4Vector::HasVector() const
{
 
 Int_t val=fV.HasVector();
 return val;
}

Int_t Nc4Vector::HasErrors() const
{
 
 Int_t val=fV.HasErrors();
 return val;
}

Double_t Nc4Vector::GetResultError() const
{
 
 return fDresult;
}

void Nc4Vector::SetScalar(Double_t v0,Double_t dv0)
{
 
 fScalar=1;
 fV0=v0;
 fV2=pow(v0,2)-fV.Dot(fV);
 SetScalarError(dv0);
}

void Nc4Vector::SetScalarError(Double_t dv0)
{
 
 fDv0=dv0;
 if (fScalar)
 {
  Double_t norm=fV.GetNorm();
  Double_t dnorm=fV.GetResultError();
  fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2));
 } 
 fDresult=0;
}

void Nc4Vector::Set3Vector(Nc3Vector& v)
{
 
 fV=v;
 if (fScalar)
 {
  SetScalar(fV0,fDv0);
 }
 else
 {
  SetInvariant(fV2,fDv2);
 }
}

void Nc4Vector::Set3Vector(Double_t* v,TString f,TString u)
{
 
 Double_t a[3];
 for (Int_t i=0; i<3; i++)
 {
  a[i]=v[i];
 }
 fV.SetVector(a,f,u);
 
 if (fScalar)
 {
  SetScalar(fV0,fDv0);
 }
 else
 {
  SetInvariant(fV2,fDv2);
 }
}

void Nc4Vector::Set3Vector(Float_t* v,TString f,TString u)
{
 
 Double_t vec[3];
 for (Int_t i=0; i<3; i++)
 {
  vec[i]=v[i];
 }
 Set3Vector(vec,f,u);
}

void Nc4Vector::Set3Vector(Double_t v1,Double_t v2,Double_t v3,TString f,TString u)
{
 
 Double_t vec[3];
 vec[0]=v1;
 vec[1]=v2;
 vec[2]=v3;
 Set3Vector(vec,f,u);
}

void Nc4Vector::SetInvariant(Double_t v2,Double_t dv2)
{
 
 fScalar=0;
 fV2=v2;
 fDv2=dv2;
 fV0=GetScalar();
 fDv0=GetResultError();
 fDresult=0;
}

void Nc4Vector::SetInvariantError(Double_t dv2)
{
 
 fDv2=dv2;
 if (!fScalar)
 {
  fDv0=GetResultError();
 }
 fDresult=0; 
}

Double_t Nc4Vector::GetInvariant()
{
 
 if (!fScalar)
 {
  fDresult=fDv2;
  return fV2;
 }
 else
 {
  Double_t inv=Dot(*this);
  return inv;
 }
}

Nc3Vector Nc4Vector::Get3Vector() const
{
 
 return fV;
}

void Nc4Vector::SetErrors(Double_t* e,TString f,TString u)
{
 
 Double_t a[3];
 for (Int_t i=0; i<3; i++)
 {
  a[i]=e[i+1];
 }
 SetScalarError(e[0]);
 fV.SetErrors(a,f,u);
}

void Nc4Vector::SetErrors(Double_t e0,Double_t e1,Double_t e2,Double_t e3,TString f,TString u)
{
 
 Double_t vec[4];
 vec[0]=e0;
 vec[1]=e1;
 vec[2]=e2;
 vec[3]=e3;
 SetErrors(vec,f,u);
}

void Nc4Vector::SetErrors(Float_t* e,TString f,TString u)
{
 
 Double_t a[4];
 for (Int_t i=0; i<4; i++)
 {
  a[i]=e[i];
 }
 SetErrors(a,f,u);
}

void Nc4Vector::GetErrors(Double_t* e,TString f,TString u)
{
 
 Double_t a[3];
 fV.GetErrors(a,f,u);
 
 // Dummy invokation of GetScalar to obtain automatic proper error determination 
 e[0]=GetScalar();
 
 e[0]=GetResultError();
 
 for (Int_t i=0; i<3; i++)
 {
  e[i+1]=a[i];
 }
}

void Nc4Vector::GetErrors(Float_t* e,TString f,TString u)
{
 
 Double_t a[4];
 GetErrors(a,f,u);
 for (Int_t i=0; i<4; i++)
 {
  e[i]=a[i];
 }
}

void Nc4Vector::Data(TString f,TString u)
{
 
 if (f=="car" || f=="sph" || f=="cyl")
 {
  Double_t vec[4],err[4];
  GetVector(vec,f,u);
  GetErrors(err,f,u);
  Double_t inv=GetInvariant();
  Double_t dinv=GetResultError();
  if (fV.HasVector())
  {
   cout << " Contravariant vector in " << f.Data() << " (" << u.Data() << ") coordinates : "
        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl;
  }
  else
  {
   cout << " 4-Vector not initialised." << endl;
  }
  if (fV.HasErrors())
  {
   cout << " ------------- Errors in " << f.Data() << " (" << u.Data() << ") coordinates : "
        << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl;
  }
  if (fV.HasVector())
  {
   cout << " --- Lorentz invariant (v^2) : " << inv << " error : " << dinv << endl;
  }
 }
 else
 {
  cout << " *Nc4Vector::Data* Unsupported frame : " << f.Data() << endl
       << "  Possible frames are 'car', 'sph' and 'cyl'." << endl; 
 }
}

Double_t Nc4Vector::Dot(Nc4Vector& q)
{
 
 Double_t dotpro=0;
 Double_t a0=GetScalar();
 Double_t da0=GetResultError();
 if ((this) == &q) // Check for special case v.Dot(v)
 {
  Double_t norm=fV.GetNorm();
  Double_t dnorm=fV.GetResultError();
  dotpro=pow(a0,2)-pow(norm,2);
  fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2));
 }
 else
 {
  Double_t b0=q.GetScalar();
  Double_t db0=q.GetResultError();
  Nc3Vector b=q.Get3Vector();
 
  Double_t dot=fV.Dot(b);
  Double_t ddot=fV.GetResultError();
 
  dotpro=a0*b0-dot;
 
  fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
 }
 
 return dotpro;
}

Nc4Vector Nc4Vector::operator+(Nc4Vector& q)
{
 
 Double_t a0=GetScalar();
 Double_t da0=GetResultError();
 Nc3Vector a=Get3Vector();
 Double_t b0=q.GetScalar();
 Double_t db0=q.GetResultError();
 Nc3Vector b=q.Get3Vector();
 
 Double_t c0=a0+b0;
 Nc3Vector c=a+b;
 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
 Nc4Vector v;
 v.SetVector(c0,c);
 v.SetScalarError(dc0);  
 return v;
}

Nc4Vector Nc4Vector::operator-(Nc4Vector& q)
{
 
 Double_t a0=GetScalar();
 Double_t da0=GetResultError();
 Nc3Vector a=Get3Vector();
 Double_t b0=q.GetScalar();
 Double_t db0=q.GetResultError();
 Nc3Vector b=q.Get3Vector();
 
 Double_t c0=a0-b0;
 Nc3Vector c=a-b;
 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
 Nc4Vector v;
 v.SetVector(c0,c);
 v.SetScalarError(dc0);  
 return v;
}

Nc4Vector Nc4Vector::operator*(Double_t s)
{
 
 Double_t a0=GetScalar();
 Double_t da0=GetResultError();
 Nc3Vector a=Get3Vector();
 
 a0*=s;
 a*=s;
 da0*=s;
 
 Nc4Vector v;
 v.SetVector(a0,a);
 v.SetScalarError(da0);  
  
 return v;
}

Nc4Vector Nc4Vector::operator/(Double_t s)
{
 
 if (fabs(s)<1.e-20) // Protect against division by 0
 {
  cout << " *Nc4Vector::/* Division by 0 detected. No action taken." << endl;
  return *this;
 }
 else
 {
  Double_t a0=GetScalar();
  Double_t da0=GetResultError();
  Nc3Vector a=Get3Vector();
 
  a0/=s;
  a/=s;
  da0/=s;
 
  Nc4Vector v;
  v.SetVector(a0,a);
  v.SetScalarError(da0);  
  
  return v;
 }
}

Nc4Vector& Nc4Vector::operator+=(Nc4Vector& q)
{
 
 Double_t a0=GetScalar();
 Double_t da0=GetResultError();
 Nc3Vector a=Get3Vector();
 Double_t b0=q.GetScalar();
 Double_t db0=q.GetResultError();
 Nc3Vector b=q.Get3Vector();
 
 Double_t c0=a0+b0;
 Nc3Vector c=a+b;
 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
 SetVector(c0,c);
 SetScalarError(dc0);  
  
 return *this;
}

Nc4Vector& Nc4Vector::operator-=(Nc4Vector& q)
{
 
 Double_t a0=GetScalar();
 Double_t da0=GetResultError();
 Nc3Vector a=Get3Vector();
 Double_t b0=q.GetScalar();
 Double_t db0=q.GetResultError();
 Nc3Vector b=q.Get3Vector();
 
 Double_t c0=a0-b0;
 Nc3Vector c=a-b;
 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
 SetVector(c0,c);
 SetScalarError(dc0);  
 
 return *this;
}

Nc4Vector& Nc4Vector::operator*=(Double_t s)
{
 
 Double_t a0=GetScalar();
 Double_t da0=GetResultError();
 Nc3Vector a=Get3Vector();
 
 a0*=s;
 a*=s;
 da0*=s;
 
 SetVector(a0,a);
 SetScalarError(da0);  
 
 return *this;
}

Nc4Vector& Nc4Vector::operator/=(Double_t s)
{
 
 if (fabs(s)<1.e-20) // Protect against division by 0
 {
  cout << " *Nc4Vector::/* Division by 0 detected. No action taken." << endl;
  return *this;
 }
 else
 {
  Double_t a0=GetScalar();
  Double_t da0=GetResultError();
  Nc3Vector a=Get3Vector();
 
  a0/=s;
  a/=s;
  da0/=s;
 
  SetVector(a0,a);
  SetScalarError(da0);  
  
  return *this;
 }
}

Int_t Nc4Vector::GetScalarFlag() const
{
 
 return fScalar;
}

Nc3Vector Nc4Vector::GetVecTrans() const
{
  
 return fV.GetVecTrans();
}

Nc3Vector Nc4Vector::GetVecLong() const
{
  
 return fV.GetVecLong();
}

Double_t Nc4Vector::GetScaTrans()
{
 
 Double_t a[3],ea[3];
 
 fV.GetVector(a,"sph");
 fV.GetErrors(ea,"sph");
 
 Double_t s=GetScalar();
 Double_t ds=GetResultError();
 
 Double_t st,dst2;
 st=s*sin(a[1]);
 dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2);
 
 fDresult=sqrt(dst2);  
 return fabs(st);
}

Double_t Nc4Vector::GetScaLong()
{
 
 Double_t a[3],ea[3];
 
 fV.GetVector(a,"sph");
 fV.GetErrors(ea,"sph");
 
 Double_t s=GetScalar();
 Double_t ds=GetResultError();
 
 Double_t sl,dsl2;
 sl=s*cos(a[1]);
 dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2);
 
 fDresult=sqrt(dsl2);  
 return sl;
}

Double_t Nc4Vector::GetPseudoRapidity()
{
 
 Double_t eta=fV.GetPseudoRapidity();
 fDresult=fV.GetResultError();
 return eta;
}

Nc3Vector Nc4Vector::GetBetaVector() const
{
 
 Nc3Vector beta;
 if (fabs(fV0)>0.) beta=fV/fV0;
 if (fabs(fDv0)>0.)
 {
  Double_t vecv[3],errv[3],errb[3];
  Double_t sqerr=0;
  fV.GetVector(vecv,"car");
  fV.GetErrors(errv,"car");
  for (Int_t i=0; i<3; i++)
  {
   sqerr=pow((errv[i]/fV0),2)+pow((vecv[i]*fDv0/(fV0*fV0)),2);
   errb[i]=sqrt(sqerr);
  }
  beta.SetErrors(errb,"car");
 }
 return beta;
}

Double_t Nc4Vector::GetBeta()
{
 
 Nc3Vector beta=GetBetaVector();
 Double_t val=beta.GetNorm();
 fDresult=beta.GetResultError();
 return val;
}

Double_t Nc4Vector::GetGamma()
{
 
 Double_t gamma=-1;
 fDresult=0;
 Double_t inv=sqrt(fabs(fV2));
 if (inv>0.)
 {
  Double_t dinv=fDv2/(2.*inv);
  gamma=fV0/inv;
  Double_t sqerr=pow((fDv0/inv),2)+pow((fV0*dinv/fV2),2);
  fDresult=sqrt(sqerr);
 }
 return gamma;
}

Double_t Nc4Vector::GetX(Int_t i,TString f,TString u)
{
 
 fDresult=0;
 
 if (i<0 || i>3) return 0;
 
 Double_t x=0;
 
 if (i==0)
 {
  x=GetScalar();
 }
 else
 {
  x=fV.GetX(i,f,u);
  fDresult=fV.GetResultError();
 }
 
 return x;
}

Double_t Nc4Vector::GetOpeningAngle(Nc4Vector& q,TString u)
{
 
 Nc3Vector v1=fV;
 Nc3Vector v2=q.Get3Vector();
 
 Double_t ang=v1.GetOpeningAngle(v2,u);
 fDresult=v1.GetResultError();
 
 return ang;
}

Double_t Nc4Vector::GetOpeningAngle(Nc3Vector& q,TString u)
{
 
 Nc3Vector v1=fV;
 
 Double_t ang=v1.GetOpeningAngle(q,u);
 fDresult=v1.GetResultError();
 
 return ang;
}

void Nc4Vector::SetUserData(NcSignal* s)
{
 
 if (fUser)
 {
  delete fUser;
  fUser=0;
 }
 
 if (s) fUser=(NcSignal*)s->Clone();
}

NcSignal* Nc4Vector::GetUserData() const
{
 
 return fUser;
}