RALICE/Ali4Vector.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // $Id$
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 // Class Ali4Vector
  20 // Handling of Lorentz 4-vectors in various reference frames.
  21 //
  22 // This class is meant to serve as a base class for ALICE objects
  23 // that have Lorentz 4-vector characteristics.
  24 // Error propagation is performed automatically.
  25 //
  26 // All 4-vectors are treated in the contravariant form and the convention
  27 // for the metric and the 4-vector components is according to the one
  28 // used in the book "Classical Electrodynamics" by J.D. Jackson.
  29 //
  30 // A 4-vector is said to have a scalar part and a 3-vector part,
  31 // which is indicated by the notation
  32 //
  33 //    x^i = (x^0,x^1,x^2,x^3)
  34 //
  35 // The scalar part   = x^0
  36 // The 3-vector part = (x^1,x^2,x^3)
  37 //
  38 // In view of accuracy and the fact that e.g. particle identity (mass)
  39 // is preserved in many physics processes, the Lorentz invariant
  40 // (x^i*x_i) is internally saved together with the scalar part.
  41 //
  42 // This allows the following two modes of functionality :
  43 //
  44 // Scalar mode    : The scalar part and the 3-vector part are considered
  45 //                  as basic quantities and the invariant with its error
  46 //                  is derived from these.
  47 // Invariant mode : The invariant and the 3-vector part are considered
  48 //                  as basic quantities and the scalar with its error
  49 //                  is derived from these.
  50 //
  51 // The philosophy followed here is the following :
  52 // ===============================================
  53 //
  54 // 1) Invokation of SetVector() sets the scalar and 3-vector parts
  55 //    and the invariant is calculated from these.
  56 //    Automatically the scalar mode is selected and invokation of
  57 //    SetErrors() will calculate the error on the invariant.
  58 //
  59 // 2) In case the scalar part is modified via SetScalar(), scalar mode is
  60 //    automatically selected and the Lorentz invariant (x^i*x_i) and its
  61 //    error are updated accordingly.
  62 //    The 3-vector part is NOT modified.
  63 //    This situation arises when one e.g. precisely determines the time
  64 //    or energy (x^0).
  65 //
  66 // 3) In case the Lorentz invariant (x^i*x_i) is modified via SetInvariant(),
  67 //    invariant mode is selected automatically and the scalar part and its
  68 //    error are updated accordingly.
  69 //    The 3-vector part is NOT modified.
  70 //    This situation arises when one e.g. precisely determines the mass.
  71 //
  72 // 4) In case the vector part is modified via Set3Vector(), then the
  73 //    current mode determines whether the scalar or the invariant is updated.
  74 //    Scalar mode    : The Lorentz invariant (x^i*x_i) and its error are updated;
  75 //                     the scalar part and its error are NOT modified.
  76 //                     This situation arises when one e.g. improves the 3-position
  77 //                     vector for a particle with a very precise timing.
  78 //    Invariant mode : The scalar part and its error are updated;
  79 //                     the Lorentz invariant (x^i*x_i) and its error are NOT modified.
  80 //                     This situation arises when one e.g. improves the 3-momentum
  81 //                     vector for a particle with known mass.
  82 //
  83 // The dotproduct is defined such that p.Dot(p) yields the Lorentz invariant
  84 // scalar of the 4-vector p (i.e. m**2 in case p is a 4-momentum).
  85 //
  86 // Note :
  87 // ------
  88 // Vectors (v), Errors (e) and reference frames (f) are specified via
  89 // SetVector(Float_t* v,TString f)
  90 // SetErrors(Float_t* e,TString f)
  91 // under the following conventions :
  92 //
  93 // f="car" ==> 3-vector part of v in Cartesian coordinates   (x,y,z)
  94 // f="sph" ==> 3-vector part of v in Spherical coordinates   (r,theta,phi)
  95 // f="cyl" ==> 3-vector part of v in Cylindrical coordinates (rho,phi,z)
  96 //
  97 // All angles are in radians.
  98 //
  99 // Example :
 100 // ---------
 101 //
 102 // Ali4Vector a;
 103 //
 104 // Float_t v[4]={25,-1,3,7};
 105 // a.SetVector(v,"car");
 106 //
 107 // Float_t vec[4];
 108 // a.GetVector(vec,"sph");
 109 //
 110 // Ali4Vector b;
 111 // Float_t v2[4]={33,6,-18,2};
 112 // b.SetVector(v2,"car");
 113 //
 114 // Float_t dotpro=a.Dot(b);
 115 //
 116 // Float_t x0=16;
 117 // Ali3Vector x;
 118 // Float_t vec2[3]={1,2,3};
 119 // x.SetVector(vec2,"car");
 120 //
 121 // Ali4Vector c;
 122 // c.SetVector(x0,x);
 123 // c.GetVector(vec,"car");
 124 // c.Info("cyl");
 125 // c=a+b;
 126 // c=a-b;
 127 // c=a*5;
 128 //
 129 //--- Author: Nick van Eijndhoven 01-apr-1999 UU-SAP Utrecht
 130 //- Modified: NvE $Date$ UU-SAP Utrecht
 131 ///////////////////////////////////////////////////////////////////////////
 132
 133 #include "Ali4Vector.h"
 134
 135 ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
 136
 137 Ali4Vector::Ali4Vector()
 138 {
 139 // Creation of a contravariant 4-vector and initialisation of parameters.
 140 // All values are initialised to 0.
 141 // Scalar mode is initially selected.
 142  fScalar=1;
 143  fV2=0;
 144  fDv2=0;
 145  fV0=0;
 146  fDv0=0;
 147  fDresult=0;
 148  Double_t a[3]={0,0,0};
 149  fV.SetVector(a,"sph");
 150 }
 151 ///////////////////////////////////////////////////////////////////////////
 152 Ali4Vector::~Ali4Vector()
 153 {
 154 // Destructor to delete dynamically allocated memory
 155 }
 156 ///////////////////////////////////////////////////////////////////////////
 157 void Ali4Vector::SetVector(Double_t v0,Ali3Vector v)
 158 {
 159 // Store contravariant vector.
 160 // The error on the scalar part is initialised to 0.
 161 // The errors on the vector part are taken from the input Ali3Vector.
 162 // Scalar mode is automatically selected.
 163 // The error on scalar result operations is reset to 0.
 164  fScalar=1;
 165  fV0=v0;
 166  fV=v;
 167  fV2=pow(v0,2)-fV.Dot(fV);
 168  SetScalarError(0);
 169 }
 170 ///////////////////////////////////////////////////////////////////////////
 171 void Ali4Vector::SetVector(Double_t* v,TString f)
 172 {
 173 // Store vector according to reference frame f.
 174 // All errors are initialised to 0.
 175 // Scalar mode is automatically selected.
 176 // The error on scalar result operations is reset to 0.
 177  fScalar=1;
 178  Double_t a[3];
 179  for (Int_t i=0; i<3; i++)
 180  {
 181   a[i]=v[i+1];
 182  }
 183  fV0=v[0];
 184  fV.SetVector(a,f);
 185  fV2=pow(fV0,2)-fV.Dot(fV);
 186  fDv2=0;
 187  fDv0=0;
 188  fDresult=0;
 189 }
 190 ///////////////////////////////////////////////////////////////////////////
 191 void Ali4Vector::GetVector(Double_t* v,TString f)
 192 {
 193 // Provide 4-vector components according to reference frame f
 194 // and according to the current mode.
 195 // Scalar mode    : The scalar part is directly returned via v[0].
 196 // Invariant mode : The scalar part is re-calculated via the value
 197 //                  of the Lorentz invariant and then returned via v[0].
 198  if (fScalar)
 199  {
 200   v[0]=fV0;
 201  }
 202  else
 203  {
 204   v[0]=sqrt(fV.Dot(fV)+fV2);
 205  }
 206  Double_t a[3];
 207  fV.GetVector(a,f);
 208  for (Int_t i=0; i<3; i++)
 209  {
 210   v[i+1]=a[i];
 211  }
 212 }
 213 ///////////////////////////////////////////////////////////////////////////
 214 void Ali4Vector::SetVector(Float_t* v,TString f)
 215 {
 216 // Store vector according to reference frame f.
 217 // All errors are initialised to 0.
 218 // Scalar mode is automatically selected.
 219 // The error on scalar result operations is reset to 0.
 220  Double_t vec[4];
 221  for (Int_t i=0; i<4; i++)
 222  {
 223   vec[i]=v[i];
 224  }
 225  SetVector(vec,f);
 226 }
 227 ///////////////////////////////////////////////////////////////////////////
 228 void Ali4Vector::GetVector(Float_t* v,TString f)
 229 {
 230 // Provide 4-vector components according to reference frame f
 231 // and according to the current mode.
 232 // Scalar mode    : The scalar part is directly returned via v[0].
 233 // Invariant mode : The scalar part is re-calculated via the value
 234 //                  of the Lorentz invariant and then returned via v[0].
 235  Double_t vec[4];
 236  GetVector(vec,f);
 237  for (Int_t i=0; i<4; i++)
 238  {
 239   v[i]=vec[i];
 240  }
 241 }
 242 ///////////////////////////////////////////////////////////////////////////
 243 Double_t Ali4Vector::GetScalar()
 244 {
 245 // Provide the scalar part.
 246 // The error on the scalar value is available via GetResultError()
 247 // after invokation of GetScalar().
 248  if (fScalar)
 249  {
 250   fDresult=fDv0;
 251   return fV0;
 252  }
 253  else
 254  {
 255   Double_t dot=fV.Dot(fV);
 256   Double_t ddot=fV.GetResultError();
 257   Double_t v02=dot+fV2;
 258   Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2));
 259   Double_t v0=sqrt(fabs(v02));
 260   Double_t dv0=0;
 261   if (v0) dv0=dv02/(2.*v0);
 262   fDresult=dv0;
 263   return v0;
 264  }
 265 }
 266 ///////////////////////////////////////////////////////////////////////////
 267 Double_t Ali4Vector::GetResultError()
 268 {
 269 // Provide the error on the result of an operation yielding a scalar
 270 // E.g. GetScalar(), GetInvariant() or Dot()
 271  return fDresult;
 272 }
 273 ///////////////////////////////////////////////////////////////////////////
 274 void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
 275 {
 276 // Modify the scalar part (v0) and its error (dv0).
 277 // The default value for dv0 is 0.
 278 // The vector part is not modified.
 279 // Scalar mode is automatically selected
 280 // ==> Lorentz invariant and its error are updated.
 281 // The error on scalar result operations is reset to 0.
 282  fScalar=1;
 283  fV0=v0;
 284  fV2=pow(v0,2)-fV.Dot(fV);
 285  SetScalarError(dv0);
 286 }
 287 ///////////////////////////////////////////////////////////////////////////
 288 void Ali4Vector::SetScalarError(Double_t dv0)
 289 {
 290 // Set the error on the scalar part.
 291 // If in scalar mode, update error on the invariant accordingly.
 292 // The error on scalar result operations is reset to 0.
 293  fDv0=dv0;
 294  if (fScalar)
 295  {
 296   Double_t norm=fV.GetNorm();
 297   Double_t dnorm=fV.GetResultError();
 298   fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2));
 299  }
 300  fDresult=0;
 301 }
 302 ///////////////////////////////////////////////////////////////////////////
 303 void Ali4Vector::Set3Vector(Ali3Vector v)
 304 {
 305 // Set the 3-vector part, the errors are taken from the input Ali3Vector
 306 // Scalar mode    : The scalar part and its error are not modified,
 307 //                  the Lorentz invariantand its error are re-calculated.
 308 // Invariant mode : The Lorentz invariant and its error are not modified,
 309 //                  the scalar part and its error are re-calculated.
 310 // The error on scalar result operations is reset to 0.
 311  fV=v;
 312  if (fScalar)
 313  {
 314   SetScalar(fV0,fDv0);
 315  }
 316  else
 317  {
 318   SetInvariant(fV2,fDv2);
 319  }
 320 }
 321 ///////////////////////////////////////////////////////////////////////////
 322 void Ali4Vector::Set3Vector(Double_t* v,TString f)
 323 {
 324 // Set the 3-vector part according to reference frame f
 325 // The errors on the vector part are initialised to 0
 326 // Scalar mode    : The scalar part and its error are not modified,
 327 //                  the Lorentz invariantand its error are re-calculated.
 328 // Invariant mode : The Lorentz invariant and its error are not modified,
 329 //                  the scalar part and its error are re-calculated.
 330 // The error on scalar result operations is reset to 0.
 331  Double_t a[3];
 332  for (Int_t i=0; i<3; i++)
 333  {
 334   a[i]=v[i];
 335  }
 336  fV.SetVector(a,f);
 337
 338  if (fScalar)
 339  {
 340   SetScalar(fV0,fDv0);
 341  }
 342  else
 343  {
 344   SetInvariant(fV2,fDv2);
 345  }
 346 }
 347 ///////////////////////////////////////////////////////////////////////////
 348 void Ali4Vector::Set3Vector(Float_t* v,TString f)
 349 {
 350 // Set the 3-vector part according to reference frame f
 351 // The errors on the vector part are initialised to 0
 352 // The Lorentz invariant is not modified
 353 // The error on scalar result operations is reset to 0.
 354  Double_t vec[3];
 355  for (Int_t i=0; i<3; i++)
 356  {
 357   vec[i]=v[i];
 358  }
 359  Set3Vector(vec,f);
 360 }
 361 ///////////////////////////////////////////////////////////////////////////
 362 void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
 363 {
 364 // Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2).
 365 // The default value for the error dv2 is 0.
 366 // The vector part is not modified.
 367 // Invariant mode is automatically selected
 368 // ==> the scalar part and its error are updated.
 369 // The error on scalar result operations is reset to 0.
 370 //
 371  fScalar=0;
 372  fV2=v2;
 373  fDv2=dv2;
 374  fV0=GetScalar();
 375  fDv0=GetResultError();
 376  fDresult=0;
 377 }
 378 ///////////////////////////////////////////////////////////////////////////
 379 void Ali4Vector::SetInvariantError(Double_t dv2)
 380 {
 381 // Set the error on the Lorentz invariant.
 382 // If in invariant mode, update error on the scalar part accordingly.
 383 // The error on scalar result operations is reset to 0.
 384  fDv2=dv2;
 385  if (!fScalar)
 386  {
 387   fDv0=GetResultError();
 388  }
 389  fDresult=0;
 390 }
 391 ///////////////////////////////////////////////////////////////////////////
 392 Double_t Ali4Vector::GetInvariant()
 393 {
 394 // Provide the Lorentz invariant v^i*v_i.
 395 // The error on the Lorentz invariant is available via GetResultError()
 396 // after invokation of GetInvariant().
 397  if (!fScalar)
 398  {
 399   fDresult=fDv2;
 400   return fV2;
 401  }
 402  else
 403  {
 404   Double_t inv=Dot(*this);
 405   return inv;
 406  }
 407 }
 408 ///////////////////////////////////////////////////////////////////////////
 409 Ali3Vector Ali4Vector::Get3Vector()
 410 {
 411 // Provide the 3-vector part
 412  return fV;
 413 }
 414 ///////////////////////////////////////////////////////////////////////////
 415 void Ali4Vector::SetErrors(Double_t* e,TString f)
 416 {
 417 // Store errors for vector v^i according to reference frame f
 418 // If in scalar mode, update error on the invariant accordingly.
 419 // The error on scalar result operations is reset to 0.
 420  Double_t a[3];
 421  for (Int_t i=0; i<3; i++)
 422  {
 423   a[i]=e[i+1];
 424  }
 425  SetScalarError(e[0]);
 426  fV.SetErrors(a,f);
 427 }
 428 ///////////////////////////////////////////////////////////////////////////
 429 void Ali4Vector::SetErrors(Float_t* e,TString f)
 430 {
 431 // Store errors for vector v^i according to reference frame f
 432 // If in scalar mode, update error on the invariant accordingly.
 433 // The error on scalar result operations is reset to 0.
 434  Double_t a[4];
 435  for (Int_t i=0; i<4; i++)
 436  {
 437   a[i]=e[i];
 438  }
 439  SetErrors(a,f);
 440 }
 441 ///////////////////////////////////////////////////////////////////////////
 442 void Ali4Vector::GetErrors(Double_t* e,TString f)
 443 {
 444 // Provide errors for vector v^i according to reference frame f
 445 // and according to the current mode.
 446 // Scalar mode    : The error on the scalar part is directly returned via e[0].
 447 // Invariant mode : The error on the scalar part is re-calculated via the error
 448 //                  value on the Lorentz invariant and then returned via e[0].
 449  Double_t a[3];
 450  fV.GetErrors(a,f);
 451
 452  e[0]=GetResultError();
 453  for (Int_t i=0; i<3; i++)
 454  {
 455   e[i+1]=a[i];
 456  }
 457 }
 458 ///////////////////////////////////////////////////////////////////////////
 459 void Ali4Vector::GetErrors(Float_t* e,TString f)
 460 {
 461 // Provide errors for vector v^i according to reference frame f
 462 // and according to the current mode.
 463 // Scalar mode    : The error on the scalar part is directly returned via e[0].
 464 // Invariant mode : The error on the scalar part is re-calculated via the error
 465 //                  value on the Lorentz invariant and then returned via e[0].
 466  Double_t a[4];
 467  GetErrors(a,f);
 468  for (Int_t i=0; i<4; i++)
 469  {
 470   e[i]=a[i];
 471  }
 472 }
 473 ///////////////////////////////////////////////////////////////////////////
 474 void Ali4Vector::Info(TString f)
 475 {
 476 // Print contravariant vector components and errors according to
 477 // reference frame f and according to the current mode.
 478 // Scalar mode    : The scalar part and its error are directly returned.
 479 // Invariant mode : The scalar part and its error are re-calculated via the
 480 //                  value (and error) of the Lorentz invariant.
 481
 482  if (f=="car" || f=="sph" || f=="cyl")
 483  {
 484   Double_t vec[4],err[4];
 485   GetVector(vec,f);
 486   GetErrors(err,f);
 487   Double_t inv=GetInvariant();
 488   Double_t dinv=GetResultError();
 489   cout << " Contravariant vector in " << f << " coordinates : "
 490        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl;
 491   cout << " ------------- Errors in " << f << " coordinates : "
 492        << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl;
 493   cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
 494  }
 495  else
 496  {
 497   cout << " *Ali4Vector::Info* Unsupported frame : " << f << endl
 498        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl;
 499  }
 500 }
 501 ///////////////////////////////////////////////////////////////////////////
 502 Double_t Ali4Vector::Dot(Ali4Vector& q)
 503 {
 504 // Provide the dot product of the current vector with vector q
 505  Double_t dotpro=0;
 506  Double_t a0=GetScalar();
 507  Double_t da0=GetResultError();
 508  if ((this) == &q) // Check for special case v.Dot(v)
 509  {
 510   Double_t norm=fV.GetNorm();
 511   Double_t dnorm=fV.GetResultError();
 512   dotpro=pow(a0,2)-pow(norm,2);
 513   fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2));
 514  }
 515  else
 516  {
 517   Double_t b0=q.GetScalar();
 518   Double_t db0=q.GetResultError();
 519   Ali3Vector b=q.Get3Vector();
 520
 521   Double_t dot=fV.Dot(b);
 522   Double_t ddot=fV.GetResultError();
 523
 524   dotpro=a0*b0-dot;
 525
 526   fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
 527  }
 528
 529  return dotpro;
 530 }
 531 ///////////////////////////////////////////////////////////////////////////
 532 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
 533 {
 534 // Add 4-vector q to the current 4-vector
 535 // Error propagation is performed automatically
 536  Double_t a0=GetScalar();
 537  Double_t da0=GetResultError();
 538  Ali3Vector a=Get3Vector();
 539  Double_t b0=q.GetScalar();
 540  Double_t db0=q.GetResultError();
 541  Ali3Vector b=q.Get3Vector();
 542
 543  Double_t c0=a0+b0;
 544  Ali3Vector c=a+b;
 545  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 546
 547  Ali4Vector v;
 548  v.SetVector(c0,c);
 549  v.SetScalarError(dc0);
 550  return v;
 551 }
 552 ///////////////////////////////////////////////////////////////////////////
 553 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
 554 {
 555 // Subtract 4-vector q from the current 4-vector
 556 // Error propagation is performed automatically
 557  Double_t a0=GetScalar();
 558  Double_t da0=GetResultError();
 559  Ali3Vector a=Get3Vector();
 560  Double_t b0=q.GetScalar();
 561  Double_t db0=q.GetResultError();
 562  Ali3Vector b=q.Get3Vector();
 563
 564  Double_t c0=a0-b0;
 565  Ali3Vector c=a-b;
 566  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 567
 568  Ali4Vector v;
 569  v.SetVector(c0,c);
 570  v.SetScalarError(dc0);
 571  return v;
 572 }
 573 ///////////////////////////////////////////////////////////////////////////
 574 Ali4Vector Ali4Vector::operator*(Double_t s)
 575 {
 576 // Multiply the current 4-vector with a scalar s
 577 // Error propagation is performed automatically
 578  Double_t a0=GetScalar();
 579  Double_t da0=GetResultError();
 580  Ali3Vector a=Get3Vector();
 581
 582  a0*=s;
 583  a*=s;
 584  da0*=s;
 585
 586  Ali4Vector v;
 587  v.SetVector(a0,a);
 588  v.SetScalarError(da0);
 589
 590  return v;
 591 }
 592 ///////////////////////////////////////////////////////////////////////////
 593 Ali4Vector Ali4Vector::operator/(Double_t s)
 594 {
 595 // Divide the current vector by a scalar s
 596 // Error propagation is performed automatically
 597
 598  if (fabs(s)<1.e-20) // Protect against division by 0
 599  {
 600   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
 601   return *this;
 602  }
 603  else
 604  {
 605   Double_t a0=GetScalar();
 606   Double_t da0=GetResultError();
 607   Ali3Vector a=Get3Vector();
 608
 609   a0/=s;
 610   a/=s;
 611   da0/=s;
 612
 613   Ali4Vector v;
 614   v.SetVector(a0,a);
 615   v.SetScalarError(da0);
 616
 617   return v;
 618  }
 619 }
 620 ///////////////////////////////////////////////////////////////////////////
 621 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
 622 {
 623 // Add 4-vector q to the current 4-vector
 624 // Error propagation is performed automatically
 625  Double_t a0=GetScalar();
 626  Double_t da0=GetResultError();
 627  Ali3Vector a=Get3Vector();
 628  Double_t b0=q.GetScalar();
 629  Double_t db0=q.GetResultError();
 630  Ali3Vector b=q.Get3Vector();
 631
 632  Double_t c0=a0+b0;
 633  Ali3Vector c=a+b;
 634  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 635
 636  SetVector(c0,c);
 637  SetScalarError(dc0);
 638
 639  return *this;
 640 }
 641 ///////////////////////////////////////////////////////////////////////////
 642 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
 643 {
 644 // Subtract 4-vector q from the current 4-vector
 645 // Error propagation is performed automatically
 646  Double_t a0=GetScalar();
 647  Double_t da0=GetResultError();
 648  Ali3Vector a=Get3Vector();
 649  Double_t b0=q.GetScalar();
 650  Double_t db0=q.GetResultError();
 651  Ali3Vector b=q.Get3Vector();
 652
 653  Double_t c0=a0-b0;
 654  Ali3Vector c=a-b;
 655  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 656
 657  SetVector(c0,c);
 658  SetScalarError(dc0);
 659
 660  return *this;
 661 }
 662 ///////////////////////////////////////////////////////////////////////////
 663 Ali4Vector& Ali4Vector::operator*=(Double_t s)
 664 {
 665 // Multiply the current 4-vector with a scalar s
 666 // Error propagation is performed automatically
 667  Double_t a0=GetScalar();
 668  Double_t da0=GetResultError();
 669  Ali3Vector a=Get3Vector();
 670
 671  a0*=s;
 672  a*=s;
 673  da0*=s;
 674
 675  SetVector(a0,a);
 676  SetScalarError(da0);
 677
 678  return *this;
 679 }
 680 ///////////////////////////////////////////////////////////////////////////
 681 Ali4Vector& Ali4Vector::operator/=(Double_t s)
 682 {
 683 // Divide the current vector by a scalar s
 684 // Error propagation is performed automatically
 685
 686  if (fabs(s)<1.e-20) // Protect against division by 0
 687  {
 688   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
 689   return *this;
 690  }
 691  else
 692  {
 693   Double_t a0=GetScalar();
 694   Double_t da0=GetResultError();
 695   Ali3Vector a=Get3Vector();
 696
 697   a0/=s;
 698   a/=s;
 699   da0/=s;
 700
 701   SetVector(a0,a);
 702   SetScalarError(da0);
 703
 704   return *this;
 705  }
 706 }
 707 ///////////////////////////////////////////////////////////////////////////
 708 Int_t Ali4Vector::GetScalarFlag()
 709 {
 710 // Provide the value of the fScalar flag (for internal use only).
 711  return fScalar;
 712 }
 713 ///////////////////////////////////////////////////////////////////////////
 714 Ali3Vector Ali4Vector::GetVecTrans()
 715 {
 716 // Provide the transverse vector part w.r.t. z-axis.
 717 // Error propagation is performed automatically
 718
 719  return fV.GetVecTrans();
 720 }
 721 ///////////////////////////////////////////////////////////////////////////
 722 Ali3Vector Ali4Vector::GetVecLong()
 723 {
 724 // Provide the longitudinal vector part w.r.t. z-axis.
 725 // Error propagation is performed automatically
 726
 727  return fV.GetVecLong();
 728 }
 729 ///////////////////////////////////////////////////////////////////////////
 730 Double_t Ali4Vector::GetScaTrans()
 731 {
 732 // Provide the "transverse value" of the scalar part w.r.t. z-axis.
 733 // This provides a basis for e.g. E_trans calculation.
 734 // Note : the returned value is always positive or zero.
 735 // The error on the value is available via GetResultError()
 736 // after invokation of GetScaTrans().
 737  Double_t a[3],ea[3];
 738
 739  fV.GetVector(a,"sph");
 740  fV.GetErrors(ea,"sph");
 741
 742  Double_t s=GetScalar();
 743  Double_t ds=GetResultError();
 744
 745  Double_t st,dst2;
 746  st=s*sin(a[1]);
 747  dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2);
 748
 749  fDresult=sqrt(dst2);
 750  return fabs(st);
 751 }
 752 ///////////////////////////////////////////////////////////////////////////
 753 Double_t Ali4Vector::GetScaLong()
 754 {
 755 // Provide the "longitudinal value" of the scalar part w.r.t. z-axis.
 756 // This provides a basis for e.g. E_long calculation.
 757 // Note : the returned value can also be negative.
 758 // The error on the value is available via GetResultError()
 759 // after invokation of GetScaLong().
 760  Double_t a[3],ea[3];
 761
 762  fV.GetVector(a,"sph");
 763  fV.GetErrors(ea,"sph");
 764
 765  Double_t s=GetScalar();
 766  Double_t ds=GetResultError();
 767
 768  Double_t sl,dsl2;
 769  sl=s*cos(a[1]);
 770  dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2);
 771
 772  fDresult=sqrt(dsl2);
 773  return sl;
 774 }
 775 ///////////////////////////////////////////////////////////////////////////
 776 Double_t Ali4Vector::GetPseudoRapidity()
 777 {
 778 // Provide the pseudorapidity value of the vector part w.r.t. z-axis.
 779 // The error on the value is available via GetResultError()
 780 // after invokation of GetPseudoRapidity().
 781  Double_t eta=fV.GetPseudoRapidity();
 782  fDresult=fV.GetResultError();
 783  return eta;
 784 }
 785 ///////////////////////////////////////////////////////////////////////////