]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - RALICE/Ali4Vector.cxx
12-oct-2005 NvE CleanTasks() invoked before execution of sub-tasks in IceConvert...
[u/mrichter/AliRoot.git] / RALICE / Ali4Vector.cxx
index a62b302f07798874d39986f51a1b40df2e0f736e..0253a4d2bba94b98353a01ca7b9247dcaa3e4e70 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-*/
+// $Id$
+
+///////////////////////////////////////////////////////////////////////////
+// Class Ali4Vector
+// Handling of Lorentz 4-vectors in various reference frames.
+//
+// This class is meant to serve as a base class for ALICE objects
+// that have Lorentz 4-vector characteristics.
+// Error propagation is performed automatically.
+//
+// All 4-vectors are treated in the contravariant form and the convention
+// for the metric and the 4-vector components is according to the one
+// used in the book "Classical Electrodynamics" by J.D. Jackson.
+//
+// A 4-vector is said to have a scalar part and a 3-vector part,
+// which is indicated by the notation
+//
+//    x^i = (x^0,x^1,x^2,x^3)
+// 
+// The scalar part   = x^0
+// The 3-vector part = (x^1,x^2,x^3)
+//
+// In view of accuracy and the fact that e.g. particle identity (mass)
+// is preserved in many physics processes, the Lorentz invariant 
+// (x^i*x_i) is internally saved together with the scalar part.
+//
+// This allows the following two modes of functionality :
+//
+// Scalar mode    : The scalar part and the 3-vector part are considered
+//                  as basic quantities and the invariant with its error
+//                  is derived from these.  
+// Invariant mode : The invariant and the 3-vector part are considered
+//                  as basic quantities and the scalar with its error
+//                  is derived from these.
+//
+// The philosophy followed here is the following :
+// ===============================================
+//
+// 1) Invokation of SetVector() sets the scalar and 3-vector parts 
+//    and the invariant is calculated from these.
+//    Automatically the scalar mode is selected and invokation of
+//    SetErrors() will calculate the error on the invariant.
+//
+// 2) In case the scalar part is modified via SetScalar(), scalar mode is
+//    automatically selected and the Lorentz invariant (x^i*x_i) and its
+//    error are updated accordingly.
+//    The 3-vector part is NOT modified.
+//    This situation arises when one e.g. precisely determines the time
+//    or energy (x^0).     
+//
+// 3) In case the Lorentz invariant (x^i*x_i) is modified via SetInvariant(),
+//    invariant mode is selected automatically and the scalar part and its
+//    error are updated accordingly.
+//    The 3-vector part is NOT modified.
+//    This situation arises when one e.g. precisely determines the mass.     
+//
+// 4) In case the vector part is modified via Set3Vector(), then the 
+//    current mode determines whether the scalar or the invariant is updated. 
+//    Scalar mode    : The Lorentz invariant (x^i*x_i) and its error are updated;
+//                     the scalar part and its error are NOT modified. 
+//                     This situation arises when one e.g. improves the 3-position
+//                     vector for a particle with a very precise timing.     
+//    Invariant mode : The scalar part and its error are updated;
+//                     the Lorentz invariant (x^i*x_i) and its error are NOT modified.
+//                     This situation arises when one e.g. improves the 3-momentum
+//                     vector for a particle with known mass.     
+//
+// The dotproduct is defined such that p.Dot(p) yields the Lorentz invariant
+// scalar of the 4-vector p (i.e. m**2 in case p is a 4-momentum).   
+//
+// Note :
+// ------
+// Vectors (v), Errors (e) and reference frames (f) are specified via
+// SetVector(Float_t* v,TString f)
+// SetErrors(Float_t* e,TString f)
+// under the following conventions :
+//
+// f="car" ==> 3-vector part of v in Cartesian coordinates   (x,y,z)
+// f="sph" ==> 3-vector part of v in Spherical coordinates   (r,theta,phi)
+// f="cyl" ==> 3-vector part of v in Cylindrical coordinates (rho,phi,z)
+//
+// All angles are in radians.
+//
+// Example :
+// ---------
+//
+// Ali4Vector a;
+//
+// Float_t v[4]={25,-1,3,7};
+// a.SetVector(v,"car");
+//
+// Float_t vec[4];
+// a.GetVector(vec,"sph");
+//
+// Ali4Vector b;
+// Float_t v2[4]={33,6,-18,2};
+// b.SetVector(v2,"car");
+//
+// Float_t dotpro=a.Dot(b);
+//
+// Float_t x0=16;
+// Ali3Vector x;
+// Float_t vec2[3]={1,2,3};
+// x.SetVector(vec2,"car");
+//
+// Ali4Vector c;
+// c.SetVector(x0,x);
+// c.GetVector(vec,"car");
+// c.Data("cyl");
+// c=a+b;
+// c=a-b;
+// c=a*5;
+//
+//--- Author: Nick van Eijndhoven 01-apr-1999 UU-SAP Utrecht
+//- Modified: NvE $Date$ UU-SAP Utrecht
+///////////////////////////////////////////////////////////////////////////
 
 #include "Ali4Vector.h"
+#include "Riostream.h"
  
 ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
  
 Ali4Vector::Ali4Vector()
 {
-// Creation of a contravariant 4-vector and initialisation of parameters
- fV0=0;
- Double_t a[3]={0,0,0};
- fV.SetVector(a,"sph");
+// Creation of a contravariant 4-vector and initialisation of parameters.
+// All values are initialised to 0.
+// Scalar mode is initially selected.
+ SetZero();
+ fScalar=1;
 }
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector::~Ali4Vector()
@@ -34,29 +149,93 @@ Ali4Vector::~Ali4Vector()
 // Destructor to delete dynamically allocated memory
 }
 ///////////////////////////////////////////////////////////////////////////
-void Ali4Vector::SetVector(Double_t v0,Ali3Vector v)
+Ali4Vector::Ali4Vector(const Ali4Vector& v)
+{
+// Copy constructor
+ fScalar=v.fScalar;
+ fV2=v.fV2;
+ fDv2=v.fDv2;
+ fV0=v.fV0;
+ fDv0=v.fDv0;
+ fDresult=v.fDresult;
+ fV=v.fV;
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::Load(Ali4Vector& q)
+{
+// Load all attributes of the input Ali4Vector into this Ali4Vector object.
+ Int_t temp1=q.GetScalarFlag();
+ Double_t temp2=q.GetResultError();
+ Double_t a[4];
+ q.GetVector(a,"sph");
+ SetVector(a,"sph");
+ q.GetErrors(a,"car");
+ SetErrors(a,"car");
+ fScalar=temp1;
+ fDresult=temp2;
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetZero()
+{
+// (Re)set all attributes to zero.
+// Note : The (de)selection of the scalar mode is not modified.
+ fV2=0;
+ fDv2=0;
+ fV0=0;
+ fDv0=0;
+ fDresult=0;
+ fV.SetZero();
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v)
 {
-// Store contravariant vector
+// Store contravariant vector.
+// The error on the scalar part is initialised to 0.
+// The errors on the vector part are taken from the input Ali3Vector.
+// Scalar mode is automatically selected. 
+// The error on scalar result operations is reset to 0.
+ fScalar=1;
  fV0=v0;
  fV=v;
+ fV2=pow(v0,2)-fV.Dot(fV);
+ SetScalarError(0);
 }
 ///////////////////////////////////////////////////////////////////////////
 void Ali4Vector::SetVector(Double_t* v,TString f)
 {
-// Store vector according to reference frame f
- fV0=v[0];
+// Store vector according to reference frame f.
+// All errors are initialised to 0.
+// Scalar mode is automatically selected. 
+// The error on scalar result operations is reset to 0.
+ 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);
+ fV2=pow(fV0,2)-fV.Dot(fV);
+ fDv2=0;
+ fDv0=0;
+ fDresult=0;
 }
 ///////////////////////////////////////////////////////////////////////////
 void Ali4Vector::GetVector(Double_t* v,TString f)
 {
-// Provide vector according to reference frame f
- v[0]=fV0;
+// Provide 4-vector components according to reference frame f
+// and according to the current mode.
+// Scalar mode    : The scalar part is directly returned via v[0].
+// Invariant mode : The scalar part is re-calculated via the value
+//                  of the Lorentz invariant and then returned via v[0].
+ if (fScalar)
+ {
+  v[0]=fV0;
+ }
+ else
+ {
+  v[0]=sqrt(fV.Dot(fV)+fV2);
+ } 
  Double_t a[3];
  fV.GetVector(a,f);
  for (Int_t i=0; i<3; i++)
@@ -67,7 +246,10 @@ void Ali4Vector::GetVector(Double_t* v,TString f)
 ///////////////////////////////////////////////////////////////////////////
 void Ali4Vector::SetVector(Float_t* v,TString f)
 {
-// Store vector according to reference frame f
+// Store vector according to reference frame f.
+// All errors are initialised to 0.
+// Scalar mode is automatically selected. 
+// The error on scalar result operations is reset to 0.
  Double_t vec[4];
  for (Int_t i=0; i<4; i++)
  {
@@ -78,7 +260,11 @@ void Ali4Vector::SetVector(Float_t* v,TString f)
 ///////////////////////////////////////////////////////////////////////////
 void Ali4Vector::GetVector(Float_t* v,TString f)
 {
-// Provide vector according to reference frame f
+// Provide 4-vector components according to reference frame f
+// and according to the current mode.
+// Scalar mode    : The scalar part is directly returned via v[0].
+// Invariant mode : The scalar part is re-calculated via the value
+//                  of the Lorentz invariant and then returned via v[0].
  Double_t vec[4];
  GetVector(vec,f);
  for (Int_t i=0; i<4; i++)
@@ -89,29 +275,263 @@ void Ali4Vector::GetVector(Float_t* v,TString f)
 ///////////////////////////////////////////////////////////////////////////
 Double_t Ali4Vector::GetScalar()
 {
-// Provide the scalar part
- return fV0;
+// Provide the scalar part.
+// The error on the scalar value is available via GetResultError()
+// after invokation of 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;
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
-Ali3Vector Ali4Vector::Get3Vector()
+Double_t Ali4Vector::GetResultError() const
+{
+// Provide the error on the result of an operation yielding a scalar
+// E.g. GetScalar(), GetInvariant() or Dot()
+ return fDresult;
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
+{
+// Modify the scalar part (v0) and its error (dv0).
+// The default value for dv0 is 0.
+// The vector part is not modified. 
+// Scalar mode is automatically selected
+// ==> Lorentz invariant and its error are updated. 
+// The error on scalar result operations is reset to 0.
+ fScalar=1;
+ fV0=v0;
+ fV2=pow(v0,2)-fV.Dot(fV);
+ SetScalarError(dv0);
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetScalarError(Double_t dv0)
+{
+// Set the error on the scalar part.
+// If in scalar mode, update error on the invariant accordingly.
+// The error on scalar result operations is reset to 0.
+ 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 Ali4Vector::Set3Vector(Ali3Vector& v)
+{
+// Set the 3-vector part, the errors are taken from the input Ali3Vector
+// Scalar mode    : The scalar part and its error are not modified,
+//                  the Lorentz invariantand its error are re-calculated.
+// Invariant mode : The Lorentz invariant and its error are not modified,
+//                  the scalar part and its error are re-calculated.
+// The error on scalar result operations is reset to 0.
+ fV=v;
+ if (fScalar)
+ {
+  SetScalar(fV0,fDv0);
+ }
+ else
+ {
+  SetInvariant(fV2,fDv2);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::Set3Vector(Double_t* v,TString f)
+{
+// Set the 3-vector part according to reference frame f
+// The errors on the vector part are initialised to 0
+// Scalar mode    : The scalar part and its error are not modified,
+//                  the Lorentz invariantand its error are re-calculated.
+// Invariant mode : The Lorentz invariant and its error are not modified,
+//                  the scalar part and its error are re-calculated.
+// The error on scalar result operations is reset to 0.
+ Double_t a[3];
+ for (Int_t i=0; i<3; i++)
+ {
+  a[i]=v[i];
+ }
+ fV.SetVector(a,f);
+
+ if (fScalar)
+ {
+  SetScalar(fV0,fDv0);
+ }
+ else
+ {
+  SetInvariant(fV2,fDv2);
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::Set3Vector(Float_t* v,TString f)
+{
+// Set the 3-vector part according to reference frame f
+// The errors on the vector part are initialised to 0
+// The Lorentz invariant is not modified
+// The error on scalar result operations is reset to 0.
+ Double_t vec[3];
+ for (Int_t i=0; i<3; i++)
+ {
+  vec[i]=v[i];
+ }
+ Set3Vector(vec,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
+{
+// Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2).
+// The default value for the error dv2 is 0.
+// The vector part is not modified.
+// Invariant mode is automatically selected
+// ==> the scalar part and its error are updated.
+// The error on scalar result operations is reset to 0.
+//
+ fScalar=0;
+ fV2=v2;
+ fDv2=dv2;
+ fV0=GetScalar();
+ fDv0=GetResultError();
+ fDresult=0;
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetInvariantError(Double_t dv2)
+{
+// Set the error on the Lorentz invariant.
+// If in invariant mode, update error on the scalar part accordingly.
+// The error on scalar result operations is reset to 0.
+ fDv2=dv2;
+ if (!fScalar)
+ {
+  fDv0=GetResultError();
+ }
+ fDresult=0; 
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t Ali4Vector::GetInvariant()
+{
+// Provide the Lorentz invariant v^i*v_i.
+// The error on the Lorentz invariant is available via GetResultError()
+// after invokation of GetInvariant().
+ if (!fScalar)
+ {
+  fDresult=fDv2;
+  return fV2;
+ }
+ else
+ {
+  Double_t inv=Dot(*this);
+  return inv;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali4Vector::Get3Vector() const
 {
 // Provide the 3-vector part
  return fV;
 }
 ///////////////////////////////////////////////////////////////////////////
-void Ali4Vector::Info(TString f)
+void Ali4Vector::SetErrors(Double_t* e,TString f)
 {
-// Print contravariant vector components according to reference frame f
+// Store errors for vector v^i according to reference frame f
+// If in scalar mode, update error on the invariant accordingly.
+// The error on scalar result operations is reset to 0.
+ Double_t a[3];
+ for (Int_t i=0; i<3; i++)
+ {
+  a[i]=e[i+1];
+ }
+ SetScalarError(e[0]);
+ fV.SetErrors(a,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::SetErrors(Float_t* e,TString f)
+{
+// Store errors for vector v^i according to reference frame f
+// If in scalar mode, update error on the invariant accordingly.
+// The error on scalar result operations is reset to 0.
+ Double_t a[4];
+ for (Int_t i=0; i<4; i++)
+ {
+  a[i]=e[i];
+ }
+ SetErrors(a,f);
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::GetErrors(Double_t* e,TString f)
+{
+// Provide errors for vector v^i according to reference frame f
+// and according to the current mode.
+// Scalar mode    : The error on the scalar part is directly returned via e[0].
+// Invariant mode : The error on the scalar part is re-calculated via the error
+//                  value on the Lorentz invariant and then returned via e[0].
+ Double_t a[3];
+ fV.GetErrors(a,f);
+
+ // 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 Ali4Vector::GetErrors(Float_t* e,TString f)
+{
+// Provide errors for vector v^i according to reference frame f
+// and according to the current mode.
+// Scalar mode    : The error on the scalar part is directly returned via e[0].
+// Invariant mode : The error on the scalar part is re-calculated via the error
+//                  value on the Lorentz invariant and then returned via e[0].
+ Double_t a[4];
+ GetErrors(a,f);
+ for (Int_t i=0; i<4; i++)
+ {
+  e[i]=a[i];
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void Ali4Vector::Data(TString f)
+{
+// Print contravariant vector components and errors according to
+// reference frame f and according to the current mode.
+// Scalar mode    : The scalar part and its error are directly returned.
+// Invariant mode : The scalar part and its error are re-calculated via the
+//                  value (and error) of the Lorentz invariant.
+
  if (f=="car" || f=="sph" || f=="cyl")
  {
-  Double_t vec[4];
+  Double_t vec[4],err[4];
   GetVector(vec,f);
-  cout << " Contravariant vector in " << f << " coordinates : "
+  GetErrors(err,f);
+  Double_t inv=GetInvariant();
+  Double_t dinv=GetResultError();
+  cout << " Contravariant vector in " << f.Data() << " coordinates : "
        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl; 
+  cout << " ------------- Errors in " << f.Data() << " coordinates : "
+       << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl; 
+  cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
  }
  else
  {
-  cout << " *Ali4Vector::Info* Unsupported frame : " << f << endl
+  cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl
        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl; 
  }
 }
@@ -119,73 +539,90 @@ void Ali4Vector::Info(TString f)
 Double_t Ali4Vector::Dot(Ali4Vector& q)
 {
 // Provide the dot product of the current vector with vector q
- Double_t a[4],b[4];
- Double_t dotpro;
+ 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();
+  Ali3Vector b=q.Get3Vector();
 
GetVector(a,"car");
q.GetVector(b,"car");
 Double_t dot=fV.Dot(b);
 Double_t ddot=fV.GetResultError();
 
- dotpro=a[0]*b[0];
- for (Int_t i=1; i<4; i++)
- {
-  dotpro-=a[i]*b[i];
+  dotpro=a0*b0-dot;
+
+  fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
  }
+
  return dotpro;
 }
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
 {
-// Add vector q to the current vector
- Double_t a[4],b[4];
+// Add 4-vector q to the current 4-vector
+// Error propagation is performed automatically
+ Double_t a0=GetScalar();
+ Double_t da0=GetResultError();
+ Ali3Vector a=Get3Vector();
+ Double_t b0=q.GetScalar();
+ Double_t db0=q.GetResultError();
+ Ali3Vector b=q.Get3Vector();
 
- GetVector(a,"car");
- q.GetVector(b,"car");
-
- for (Int_t i=0; i<4; i++)
- {
-  a[i]+=b[i];
- }
+ Double_t c0=a0+b0;
+ Ali3Vector c=a+b;
+ Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
  Ali4Vector v;
- v.SetVector(a,"car");
-  
+ v.SetVector(c0,c);
+ v.SetScalarError(dc0);  
  return v;
 }
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
 {
-// Subtract vector q from the current vector
- Double_t a[4],b[4];
+// Subtract 4-vector q from the current 4-vector
+// Error propagation is performed automatically
+ Double_t a0=GetScalar();
+ Double_t da0=GetResultError();
+ Ali3Vector a=Get3Vector();
+ Double_t b0=q.GetScalar();
+ Double_t db0=q.GetResultError();
+ Ali3Vector b=q.Get3Vector();
 
- GetVector(a,"car");
- q.GetVector(b,"car");
-
- for (Int_t i=0; i<4; i++)
- {
-  a[i]-=b[i];
- }
+ Double_t c0=a0-b0;
+ Ali3Vector c=a-b;
+ Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
  Ali4Vector v;
- v.SetVector(a,"car");
-  
+ v.SetVector(c0,c);
+ v.SetScalarError(dc0);  
  return v;
 }
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector Ali4Vector::operator*(Double_t s)
 {
-// Multiply the current vector with a scalar s
- Double_t a[4];
-
- GetVector(a,"car");
+// Multiply the current 4-vector with a scalar s
+// Error propagation is performed automatically
+ Double_t a0=GetScalar();
+ Double_t da0=GetResultError();
+ Ali3Vector a=Get3Vector();
 
- for (Int_t i=0; i<4; i++)
- {
-  a[i]*=s;
- }
+ a0*=s;
+ a*=s;
+ da0*=s;
 
  Ali4Vector v;
- v.SetVector(a,"car");
+ v.SetVector(a0,a);
+ v.SetScalarError(da0);  
   
  return v;
 }
@@ -193,6 +630,7 @@ Ali4Vector Ali4Vector::operator*(Double_t s)
 Ali4Vector Ali4Vector::operator/(Double_t s)
 {
 // Divide the current vector by a scalar s
+// Error propagation is performed automatically
 
  if (fabs(s)<1.e-20) // Protect against division by 0
  {
@@ -201,17 +639,17 @@ Ali4Vector Ali4Vector::operator/(Double_t s)
  }
  else
  {
-  Double_t a[4];
+  Double_t a0=GetScalar();
+  Double_t da0=GetResultError();
+  Ali3Vector a=Get3Vector();
 
-  GetVector(a,"car");
-
-  for (Int_t i=0; i<4; i++)
-  {
-   a[i]/=s;
-  }
+  a0/=s;
+  a/=s;
+  da0/=s;
 
   Ali4Vector v;
-  v.SetVector(a,"car");
+  v.SetVector(a0,a);
+  v.SetScalarError(da0);  
   
   return v;
  }
@@ -219,80 +657,213 @@ Ali4Vector Ali4Vector::operator/(Double_t s)
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
 {
-// Add vector q to the current vector
- Double_t a[4],b[4];
-
- GetVector(a,"car");
- q.GetVector(b,"car");
+// Add 4-vector q to the current 4-vector
+// Error propagation is performed automatically
+ Double_t a0=GetScalar();
+ Double_t da0=GetResultError();
+ Ali3Vector a=Get3Vector();
+ Double_t b0=q.GetScalar();
+ Double_t db0=q.GetResultError();
+ Ali3Vector b=q.Get3Vector();
 
- for (Int_t i=0; i<4; i++)
- {
-  a[i]+=b[i];
- }
+ Double_t c0=a0+b0;
+ Ali3Vector c=a+b;
+ Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
- SetVector(a,"car");
+ SetVector(c0,c);
+ SetScalarError(dc0);  
   
  return *this;
 }
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
 {
-// Subtract vector q from the current vector
- Double_t a[4],b[4];
+// Subtract 4-vector q from the current 4-vector
+// Error propagation is performed automatically
+ Double_t a0=GetScalar();
+ Double_t da0=GetResultError();
+ Ali3Vector a=Get3Vector();
+ Double_t b0=q.GetScalar();
+ Double_t db0=q.GetResultError();
+ Ali3Vector b=q.Get3Vector();
 
- GetVector(a,"car");
- q.GetVector(b,"car");
+ Double_t c0=a0-b0;
+ Ali3Vector c=a-b;
+ Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 
- for (Int_t i=0; i<4; i++)
- {
-  a[i]-=b[i];
- }
+ SetVector(c0,c);
+ SetScalarError(dc0);  
 
- SetVector(a,"car");
-  
  return *this;
 }
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector& Ali4Vector::operator*=(Double_t s)
 {
-// Multiply the current vector with a scalar s
- Double_t a[4];
+// Multiply the current 4-vector with a scalar s
+// Error propagation is performed automatically
+ Double_t a0=GetScalar();
+ Double_t da0=GetResultError();
+ Ali3Vector a=Get3Vector();
 
- GetVector(a,"car");
+ a0*=s;
+ a*=s;
+ da0*=s;
 
- for (Int_t i=0; i<4; i++)
- {
-  a[i]*=s;
- }
+ SetVector(a0,a);
+ SetScalarError(da0);  
 
- SetVector(a,"car");
-  
  return *this;
 }
 ///////////////////////////////////////////////////////////////////////////
 Ali4Vector& Ali4Vector::operator/=(Double_t s)
 {
 // Divide the current vector by a scalar s
+// Error propagation is performed automatically
 
  if (fabs(s)<1.e-20) // Protect against division by 0
  {
-  cout << " *Ali4Vector::/=* Division by 0 detected. No action taken." << endl;
+  cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
   return *this;
  }
  else
  {
-  Double_t a[4];
+  Double_t a0=GetScalar();
+  Double_t da0=GetResultError();
+  Ali3Vector a=Get3Vector();
 
-  GetVector(a,"car");
-
-  for (Int_t i=0; i<4; i++)
-  {
-   a[i]/=s;
-  }
+  a0/=s;
+  a/=s;
+  da0/=s;
 
-  SetVector(a,"car");
+  SetVector(a0,a);
+  SetScalarError(da0);  
   
   return *this;
  }
 }
 ///////////////////////////////////////////////////////////////////////////
+Int_t Ali4Vector::GetScalarFlag() const
+{
+// Provide the value of the fScalar flag (for internal use only).
+ return fScalar;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali4Vector::GetVecTrans() const
+{
+// Provide the transverse vector part w.r.t. z-axis.
+// Error propagation is performed automatically
+  
+ return fV.GetVecTrans();
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali4Vector::GetVecLong() const
+{
+// Provide the longitudinal vector part w.r.t. z-axis.
+// Error propagation is performed automatically
+  
+ return fV.GetVecLong();
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t Ali4Vector::GetScaTrans()
+{
+// Provide the "transverse value" of the scalar part w.r.t. z-axis.
+// This provides a basis for e.g. E_trans calculation.
+// Note : the returned value is always positive or zero.
+// The error on the value is available via GetResultError()
+// after invokation of 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 Ali4Vector::GetScaLong()
+{
+// Provide the "longitudinal value" of the scalar part w.r.t. z-axis.
+// This provides a basis for e.g. E_long calculation.
+// Note : the returned value can also be negative.
+// The error on the value is available via GetResultError()
+// after invokation of 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 Ali4Vector::GetPseudoRapidity()
+{
+// Provide the pseudorapidity value of the vector part w.r.t. z-axis.
+// The error on the value is available via GetResultError()
+// after invokation of GetPseudoRapidity().
+ Double_t eta=fV.GetPseudoRapidity();
+ fDresult=fV.GetResultError();
+ return eta;
+}
+///////////////////////////////////////////////////////////////////////////
+Ali3Vector Ali4Vector::GetBetaVector() const
+{
+// Provide the beta 3-vector corresponding to this 4-vector.
+ Ali3Vector 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 Ali4Vector::GetBeta()
+{
+// Provide the beta value (i.e. v/c) corresponding to this 4-vector.
+ Ali3Vector beta=GetBetaVector();
+ Double_t val=beta.GetNorm();
+ fDresult=beta.GetResultError();
+ return val;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t Ali4Vector::GetGamma()
+{
+// Provide the Lorentz gamma factor corresponding to this 4-vector.
+// In case the gamma factor is infinite a value of -1 is returned.
+ 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;
+}
+///////////////////////////////////////////////////////////////////////////