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.Data("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 #include "Riostream.h"
 135
 136 ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
 137
 138 Ali4Vector::Ali4Vector()
 139 {
 140 // Creation of a contravariant 4-vector and initialisation of parameters.
 141 // All values are initialised to 0.
 142 // Scalar mode is initially selected.
 143  SetZero();
 144  fScalar=1;
 145 }
 146 ///////////////////////////////////////////////////////////////////////////
 147 Ali4Vector::~Ali4Vector()
 148 {
 149 // Destructor to delete dynamically allocated memory
 150 }
 151 ///////////////////////////////////////////////////////////////////////////
 152 Ali4Vector::Ali4Vector(const Ali4Vector& v)
 153 {
 154 // Copy constructor
 155  fScalar=v.fScalar;
 156  fV2=v.fV2;
 157  fDv2=v.fDv2;
 158  fV0=v.fV0;
 159  fDv0=v.fDv0;
 160  fDresult=v.fDresult;
 161  fV=v.fV;
 162 }
 163 ///////////////////////////////////////////////////////////////////////////
 164 void Ali4Vector::Load(Ali4Vector& q)
 165 {
 166 // Load all attributes of the input Ali4Vector into this Ali4Vector object.
 167  Int_t temp1=q.GetScalarFlag();
 168  Double_t temp2=q.GetResultError();
 169  Double_t a[4];
 170  q.GetVector(a,"sph");
 171  SetVector(a,"sph");
 172  q.GetErrors(a,"car");
 173  SetErrors(a,"car");
 174  fScalar=temp1;
 175  fDresult=temp2;
 176 }
 177 ///////////////////////////////////////////////////////////////////////////
 178 void Ali4Vector::SetZero()
 179 {
 180 // (Re)set all attributes to zero.
 181 // Note : The (de)selection of the scalar mode is not modified.
 182  fV2=0;
 183  fDv2=0;
 184  fV0=0;
 185  fDv0=0;
 186  fDresult=0;
 187  fV.SetZero();
 188 }
 189 ///////////////////////////////////////////////////////////////////////////
 190 void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v)
 191 {
 192 // Store contravariant vector.
 193 // The error on the scalar part is initialised to 0.
 194 // The errors on the vector part are taken from the input Ali3Vector.
 195 // Scalar mode is automatically selected.
 196 // The error on scalar result operations is reset to 0.
 197  fScalar=1;
 198  fV0=v0;
 199  fV=v;
 200  fV2=pow(v0,2)-fV.Dot(fV);
 201  SetScalarError(0);
 202 }
 203 ///////////////////////////////////////////////////////////////////////////
 204 void Ali4Vector::SetVector(Double_t* v,TString f)
 205 {
 206 // Store vector according to reference frame f.
 207 // All errors are initialised to 0.
 208 // Scalar mode is automatically selected.
 209 // The error on scalar result operations is reset to 0.
 210  fScalar=1;
 211  Double_t a[3];
 212  for (Int_t i=0; i<3; i++)
 213  {
 214   a[i]=v[i+1];
 215  }
 216  fV0=v[0];
 217  fV.SetVector(a,f);
 218  fV2=pow(fV0,2)-fV.Dot(fV);
 219  fDv2=0;
 220  fDv0=0;
 221  fDresult=0;
 222 }
 223 ///////////////////////////////////////////////////////////////////////////
 224 void Ali4Vector::GetVector(Double_t* v,TString f)
 225 {
 226 // Provide 4-vector components according to reference frame f
 227 // and according to the current mode.
 228 // Scalar mode    : The scalar part is directly returned via v[0].
 229 // Invariant mode : The scalar part is re-calculated via the value
 230 //                  of the Lorentz invariant and then returned via v[0].
 231  if (fScalar)
 232  {
 233   v[0]=fV0;
 234  }
 235  else
 236  {
 237   v[0]=sqrt(fV.Dot(fV)+fV2);
 238  }
 239  Double_t a[3];
 240  fV.GetVector(a,f);
 241  for (Int_t i=0; i<3; i++)
 242  {
 243   v[i+1]=a[i];
 244  }
 245 }
 246 ///////////////////////////////////////////////////////////////////////////
 247 void Ali4Vector::SetVector(Float_t* v,TString f)
 248 {
 249 // Store vector according to reference frame f.
 250 // All errors are initialised to 0.
 251 // Scalar mode is automatically selected.
 252 // The error on scalar result operations is reset to 0.
 253  Double_t vec[4];
 254  for (Int_t i=0; i<4; i++)
 255  {
 256   vec[i]=v[i];
 257  }
 258  SetVector(vec,f);
 259 }
 260 ///////////////////////////////////////////////////////////////////////////
 261 void Ali4Vector::GetVector(Float_t* v,TString f)
 262 {
 263 // Provide 4-vector components according to reference frame f
 264 // and according to the current mode.
 265 // Scalar mode    : The scalar part is directly returned via v[0].
 266 // Invariant mode : The scalar part is re-calculated via the value
 267 //                  of the Lorentz invariant and then returned via v[0].
 268  Double_t vec[4];
 269  GetVector(vec,f);
 270  for (Int_t i=0; i<4; i++)
 271  {
 272   v[i]=vec[i];
 273  }
 274 }
 275 ///////////////////////////////////////////////////////////////////////////
 276 Double_t Ali4Vector::GetScalar()
 277 {
 278 // Provide the scalar part.
 279 // The error on the scalar value is available via GetResultError()
 280 // after invokation of GetScalar().
 281  if (fScalar)
 282  {
 283   fDresult=fDv0;
 284   return fV0;
 285  }
 286  else
 287  {
 288   Double_t dot=fV.Dot(fV);
 289   Double_t ddot=fV.GetResultError();
 290   Double_t v02=dot+fV2;
 291   Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2));
 292   Double_t v0=sqrt(fabs(v02));
 293   Double_t dv0=0;
 294   if (v0) dv0=dv02/(2.*v0);
 295   fDresult=dv0;
 296   return v0;
 297  }
 298 }
 299 ///////////////////////////////////////////////////////////////////////////
 300 Double_t Ali4Vector::GetResultError() const
 301 {
 302 // Provide the error on the result of an operation yielding a scalar
 303 // E.g. GetScalar(), GetInvariant() or Dot()
 304  return fDresult;
 305 }
 306 ///////////////////////////////////////////////////////////////////////////
 307 void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
 308 {
 309 // Modify the scalar part (v0) and its error (dv0).
 310 // The default value for dv0 is 0.
 311 // The vector part is not modified.
 312 // Scalar mode is automatically selected
 313 // ==> Lorentz invariant and its error are updated.
 314 // The error on scalar result operations is reset to 0.
 315  fScalar=1;
 316  fV0=v0;
 317  fV2=pow(v0,2)-fV.Dot(fV);
 318  SetScalarError(dv0);
 319 }
 320 ///////////////////////////////////////////////////////////////////////////
 321 void Ali4Vector::SetScalarError(Double_t dv0)
 322 {
 323 // Set the error on the scalar part.
 324 // If in scalar mode, update error on the invariant accordingly.
 325 // The error on scalar result operations is reset to 0.
 326  fDv0=dv0;
 327  if (fScalar)
 328  {
 329   Double_t norm=fV.GetNorm();
 330   Double_t dnorm=fV.GetResultError();
 331   fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2));
 332  }
 333  fDresult=0;
 334 }
 335 ///////////////////////////////////////////////////////////////////////////
 336 void Ali4Vector::Set3Vector(Ali3Vector& v)
 337 {
 338 // Set the 3-vector part, the errors are taken from the input Ali3Vector
 339 // Scalar mode    : The scalar part and its error are not modified,
 340 //                  the Lorentz invariantand its error are re-calculated.
 341 // Invariant mode : The Lorentz invariant and its error are not modified,
 342 //                  the scalar part and its error are re-calculated.
 343 // The error on scalar result operations is reset to 0.
 344  fV=v;
 345  if (fScalar)
 346  {
 347   SetScalar(fV0,fDv0);
 348  }
 349  else
 350  {
 351   SetInvariant(fV2,fDv2);
 352  }
 353 }
 354 ///////////////////////////////////////////////////////////////////////////
 355 void Ali4Vector::Set3Vector(Double_t* v,TString f)
 356 {
 357 // Set the 3-vector part according to reference frame f
 358 // The errors on the vector part are initialised to 0
 359 // Scalar mode    : The scalar part and its error are not modified,
 360 //                  the Lorentz invariantand its error are re-calculated.
 361 // Invariant mode : The Lorentz invariant and its error are not modified,
 362 //                  the scalar part and its error are re-calculated.
 363 // The error on scalar result operations is reset to 0.
 364  Double_t a[3];
 365  for (Int_t i=0; i<3; i++)
 366  {
 367   a[i]=v[i];
 368  }
 369  fV.SetVector(a,f);
 370
 371  if (fScalar)
 372  {
 373   SetScalar(fV0,fDv0);
 374  }
 375  else
 376  {
 377   SetInvariant(fV2,fDv2);
 378  }
 379 }
 380 ///////////////////////////////////////////////////////////////////////////
 381 void Ali4Vector::Set3Vector(Float_t* v,TString f)
 382 {
 383 // Set the 3-vector part according to reference frame f
 384 // The errors on the vector part are initialised to 0
 385 // The Lorentz invariant is not modified
 386 // The error on scalar result operations is reset to 0.
 387  Double_t vec[3];
 388  for (Int_t i=0; i<3; i++)
 389  {
 390   vec[i]=v[i];
 391  }
 392  Set3Vector(vec,f);
 393 }
 394 ///////////////////////////////////////////////////////////////////////////
 395 void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
 396 {
 397 // Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2).
 398 // The default value for the error dv2 is 0.
 399 // The vector part is not modified.
 400 // Invariant mode is automatically selected
 401 // ==> the scalar part and its error are updated.
 402 // The error on scalar result operations is reset to 0.
 403 //
 404  fScalar=0;
 405  fV2=v2;
 406  fDv2=dv2;
 407  fV0=GetScalar();
 408  fDv0=GetResultError();
 409  fDresult=0;
 410 }
 411 ///////////////////////////////////////////////////////////////////////////
 412 void Ali4Vector::SetInvariantError(Double_t dv2)
 413 {
 414 // Set the error on the Lorentz invariant.
 415 // If in invariant mode, update error on the scalar part accordingly.
 416 // The error on scalar result operations is reset to 0.
 417  fDv2=dv2;
 418  if (!fScalar)
 419  {
 420   fDv0=GetResultError();
 421  }
 422  fDresult=0;
 423 }
 424 ///////////////////////////////////////////////////////////////////////////
 425 Double_t Ali4Vector::GetInvariant()
 426 {
 427 // Provide the Lorentz invariant v^i*v_i.
 428 // The error on the Lorentz invariant is available via GetResultError()
 429 // after invokation of GetInvariant().
 430  if (!fScalar)
 431  {
 432   fDresult=fDv2;
 433   return fV2;
 434  }
 435  else
 436  {
 437   Double_t inv=Dot(*this);
 438   return inv;
 439  }
 440 }
 441 ///////////////////////////////////////////////////////////////////////////
 442 Ali3Vector Ali4Vector::Get3Vector() const
 443 {
 444 // Provide the 3-vector part
 445  return fV;
 446 }
 447 ///////////////////////////////////////////////////////////////////////////
 448 void Ali4Vector::SetErrors(Double_t* e,TString f)
 449 {
 450 // Store errors for vector v^i according to reference frame f
 451 // If in scalar mode, update error on the invariant accordingly.
 452 // The error on scalar result operations is reset to 0.
 453  Double_t a[3];
 454  for (Int_t i=0; i<3; i++)
 455  {
 456   a[i]=e[i+1];
 457  }
 458  SetScalarError(e[0]);
 459  fV.SetErrors(a,f);
 460 }
 461 ///////////////////////////////////////////////////////////////////////////
 462 void Ali4Vector::SetErrors(Float_t* e,TString f)
 463 {
 464 // Store errors for vector v^i according to reference frame f
 465 // If in scalar mode, update error on the invariant accordingly.
 466 // The error on scalar result operations is reset to 0.
 467  Double_t a[4];
 468  for (Int_t i=0; i<4; i++)
 469  {
 470   a[i]=e[i];
 471  }
 472  SetErrors(a,f);
 473 }
 474 ///////////////////////////////////////////////////////////////////////////
 475 void Ali4Vector::GetErrors(Double_t* e,TString f)
 476 {
 477 // Provide errors for vector v^i according to reference frame f
 478 // and according to the current mode.
 479 // Scalar mode    : The error on the scalar part is directly returned via e[0].
 480 // Invariant mode : The error on the scalar part is re-calculated via the error
 481 //                  value on the Lorentz invariant and then returned via e[0].
 482  Double_t a[3];
 483  fV.GetErrors(a,f);
 484
 485  // Dummy invokation of GetScalar to obtain automatic proper error determination
 486  e[0]=GetScalar();
 487
 488  e[0]=GetResultError();
 489
 490  for (Int_t i=0; i<3; i++)
 491  {
 492   e[i+1]=a[i];
 493  }
 494 }
 495 ///////////////////////////////////////////////////////////////////////////
 496 void Ali4Vector::GetErrors(Float_t* e,TString f)
 497 {
 498 // Provide errors for vector v^i according to reference frame f
 499 // and according to the current mode.
 500 // Scalar mode    : The error on the scalar part is directly returned via e[0].
 501 // Invariant mode : The error on the scalar part is re-calculated via the error
 502 //                  value on the Lorentz invariant and then returned via e[0].
 503  Double_t a[4];
 504  GetErrors(a,f);
 505  for (Int_t i=0; i<4; i++)
 506  {
 507   e[i]=a[i];
 508  }
 509 }
 510 ///////////////////////////////////////////////////////////////////////////
 511 void Ali4Vector::Data(TString f)
 512 {
 513 // Print contravariant vector components and errors according to
 514 // reference frame f and according to the current mode.
 515 // Scalar mode    : The scalar part and its error are directly returned.
 516 // Invariant mode : The scalar part and its error are re-calculated via the
 517 //                  value (and error) of the Lorentz invariant.
 518
 519  if (f=="car" || f=="sph" || f=="cyl")
 520  {
 521   Double_t vec[4],err[4];
 522   GetVector(vec,f);
 523   GetErrors(err,f);
 524   Double_t inv=GetInvariant();
 525   Double_t dinv=GetResultError();
 526   cout << " Contravariant vector in " << f.Data() << " coordinates : "
 527        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl;
 528   cout << " ------------- Errors in " << f.Data() << " coordinates : "
 529        << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl;
 530   cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
 531  }
 532  else
 533  {
 534   cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl
 535        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl;
 536  }
 537 }
 538 ///////////////////////////////////////////////////////////////////////////
 539 Double_t Ali4Vector::Dot(Ali4Vector& q)
 540 {
 541 // Provide the dot product of the current vector with vector q
 542  Double_t dotpro=0;
 543  Double_t a0=GetScalar();
 544  Double_t da0=GetResultError();
 545  if ((this) == &q) // Check for special case v.Dot(v)
 546  {
 547   Double_t norm=fV.GetNorm();
 548   Double_t dnorm=fV.GetResultError();
 549   dotpro=pow(a0,2)-pow(norm,2);
 550   fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2));
 551  }
 552  else
 553  {
 554   Double_t b0=q.GetScalar();
 555   Double_t db0=q.GetResultError();
 556   Ali3Vector b=q.Get3Vector();
 557
 558   Double_t dot=fV.Dot(b);
 559   Double_t ddot=fV.GetResultError();
 560
 561   dotpro=a0*b0-dot;
 562
 563   fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
 564  }
 565
 566  return dotpro;
 567 }
 568 ///////////////////////////////////////////////////////////////////////////
 569 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
 570 {
 571 // Add 4-vector q to the current 4-vector
 572 // Error propagation is performed automatically
 573  Double_t a0=GetScalar();
 574  Double_t da0=GetResultError();
 575  Ali3Vector a=Get3Vector();
 576  Double_t b0=q.GetScalar();
 577  Double_t db0=q.GetResultError();
 578  Ali3Vector b=q.Get3Vector();
 579
 580  Double_t c0=a0+b0;
 581  Ali3Vector c=a+b;
 582  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 583
 584  Ali4Vector v;
 585  v.SetVector(c0,c);
 586  v.SetScalarError(dc0);
 587  return v;
 588 }
 589 ///////////////////////////////////////////////////////////////////////////
 590 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
 591 {
 592 // Subtract 4-vector q from the current 4-vector
 593 // Error propagation is performed automatically
 594  Double_t a0=GetScalar();
 595  Double_t da0=GetResultError();
 596  Ali3Vector a=Get3Vector();
 597  Double_t b0=q.GetScalar();
 598  Double_t db0=q.GetResultError();
 599  Ali3Vector b=q.Get3Vector();
 600
 601  Double_t c0=a0-b0;
 602  Ali3Vector c=a-b;
 603  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 604
 605  Ali4Vector v;
 606  v.SetVector(c0,c);
 607  v.SetScalarError(dc0);
 608  return v;
 609 }
 610 ///////////////////////////////////////////////////////////////////////////
 611 Ali4Vector Ali4Vector::operator*(Double_t s)
 612 {
 613 // Multiply the current 4-vector with a scalar s
 614 // Error propagation is performed automatically
 615  Double_t a0=GetScalar();
 616  Double_t da0=GetResultError();
 617  Ali3Vector a=Get3Vector();
 618
 619  a0*=s;
 620  a*=s;
 621  da0*=s;
 622
 623  Ali4Vector v;
 624  v.SetVector(a0,a);
 625  v.SetScalarError(da0);
 626
 627  return v;
 628 }
 629 ///////////////////////////////////////////////////////////////////////////
 630 Ali4Vector Ali4Vector::operator/(Double_t s)
 631 {
 632 // Divide the current vector by a scalar s
 633 // Error propagation is performed automatically
 634
 635  if (fabs(s)<1.e-20) // Protect against division by 0
 636  {
 637   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
 638   return *this;
 639  }
 640  else
 641  {
 642   Double_t a0=GetScalar();
 643   Double_t da0=GetResultError();
 644   Ali3Vector a=Get3Vector();
 645
 646   a0/=s;
 647   a/=s;
 648   da0/=s;
 649
 650   Ali4Vector v;
 651   v.SetVector(a0,a);
 652   v.SetScalarError(da0);
 653
 654   return v;
 655  }
 656 }
 657 ///////////////////////////////////////////////////////////////////////////
 658 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
 659 {
 660 // Add 4-vector q to the current 4-vector
 661 // Error propagation is performed automatically
 662  Double_t a0=GetScalar();
 663  Double_t da0=GetResultError();
 664  Ali3Vector a=Get3Vector();
 665  Double_t b0=q.GetScalar();
 666  Double_t db0=q.GetResultError();
 667  Ali3Vector b=q.Get3Vector();
 668
 669  Double_t c0=a0+b0;
 670  Ali3Vector c=a+b;
 671  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 672
 673  SetVector(c0,c);
 674  SetScalarError(dc0);
 675
 676  return *this;
 677 }
 678 ///////////////////////////////////////////////////////////////////////////
 679 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
 680 {
 681 // Subtract 4-vector q from the current 4-vector
 682 // Error propagation is performed automatically
 683  Double_t a0=GetScalar();
 684  Double_t da0=GetResultError();
 685  Ali3Vector a=Get3Vector();
 686  Double_t b0=q.GetScalar();
 687  Double_t db0=q.GetResultError();
 688  Ali3Vector b=q.Get3Vector();
 689
 690  Double_t c0=a0-b0;
 691  Ali3Vector c=a-b;
 692  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 693
 694  SetVector(c0,c);
 695  SetScalarError(dc0);
 696
 697  return *this;
 698 }
 699 ///////////////////////////////////////////////////////////////////////////
 700 Ali4Vector& Ali4Vector::operator*=(Double_t s)
 701 {
 702 // Multiply the current 4-vector with a scalar s
 703 // Error propagation is performed automatically
 704  Double_t a0=GetScalar();
 705  Double_t da0=GetResultError();
 706  Ali3Vector a=Get3Vector();
 707
 708  a0*=s;
 709  a*=s;
 710  da0*=s;
 711
 712  SetVector(a0,a);
 713  SetScalarError(da0);
 714
 715  return *this;
 716 }
 717 ///////////////////////////////////////////////////////////////////////////
 718 Ali4Vector& Ali4Vector::operator/=(Double_t s)
 719 {
 720 // Divide the current vector by a scalar s
 721 // Error propagation is performed automatically
 722
 723  if (fabs(s)<1.e-20) // Protect against division by 0
 724  {
 725   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
 726   return *this;
 727  }
 728  else
 729  {
 730   Double_t a0=GetScalar();
 731   Double_t da0=GetResultError();
 732   Ali3Vector a=Get3Vector();
 733
 734   a0/=s;
 735   a/=s;
 736   da0/=s;
 737
 738   SetVector(a0,a);
 739   SetScalarError(da0);
 740
 741   return *this;
 742  }
 743 }
 744 ///////////////////////////////////////////////////////////////////////////
 745 Int_t Ali4Vector::GetScalarFlag() const
 746 {
 747 // Provide the value of the fScalar flag (for internal use only).
 748  return fScalar;
 749 }
 750 ///////////////////////////////////////////////////////////////////////////
 751 Ali3Vector Ali4Vector::GetVecTrans() const
 752 {
 753 // Provide the transverse vector part w.r.t. z-axis.
 754 // Error propagation is performed automatically
 755
 756  return fV.GetVecTrans();
 757 }
 758 ///////////////////////////////////////////////////////////////////////////
 759 Ali3Vector Ali4Vector::GetVecLong() const
 760 {
 761 // Provide the longitudinal vector part w.r.t. z-axis.
 762 // Error propagation is performed automatically
 763
 764  return fV.GetVecLong();
 765 }
 766 ///////////////////////////////////////////////////////////////////////////
 767 Double_t Ali4Vector::GetScaTrans()
 768 {
 769 // Provide the "transverse value" of the scalar part w.r.t. z-axis.
 770 // This provides a basis for e.g. E_trans calculation.
 771 // Note : the returned value is always positive or zero.
 772 // The error on the value is available via GetResultError()
 773 // after invokation of GetScaTrans().
 774  Double_t a[3],ea[3];
 775
 776  fV.GetVector(a,"sph");
 777  fV.GetErrors(ea,"sph");
 778
 779  Double_t s=GetScalar();
 780  Double_t ds=GetResultError();
 781
 782  Double_t st,dst2;
 783  st=s*sin(a[1]);
 784  dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2);
 785
 786  fDresult=sqrt(dst2);
 787  return fabs(st);
 788 }
 789 ///////////////////////////////////////////////////////////////////////////
 790 Double_t Ali4Vector::GetScaLong()
 791 {
 792 // Provide the "longitudinal value" of the scalar part w.r.t. z-axis.
 793 // This provides a basis for e.g. E_long calculation.
 794 // Note : the returned value can also be negative.
 795 // The error on the value is available via GetResultError()
 796 // after invokation of GetScaLong().
 797  Double_t a[3],ea[3];
 798
 799  fV.GetVector(a,"sph");
 800  fV.GetErrors(ea,"sph");
 801
 802  Double_t s=GetScalar();
 803  Double_t ds=GetResultError();
 804
 805  Double_t sl,dsl2;
 806  sl=s*cos(a[1]);
 807  dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2);
 808
 809  fDresult=sqrt(dsl2);
 810  return sl;
 811 }
 812 ///////////////////////////////////////////////////////////////////////////
 813 Double_t Ali4Vector::GetPseudoRapidity()
 814 {
 815 // Provide the pseudorapidity value of the vector part w.r.t. z-axis.
 816 // The error on the value is available via GetResultError()
 817 // after invokation of GetPseudoRapidity().
 818  Double_t eta=fV.GetPseudoRapidity();
 819  fDresult=fV.GetResultError();
 820  return eta;
 821 }
 822 ///////////////////////////////////////////////////////////////////////////
 823 Ali3Vector Ali4Vector::GetBetaVector() const
 824 {
 825 // Provide the beta 3-vector corresponding to this 4-vector.
 826  Ali3Vector beta;
 827  if (fabs(fV0)>0.) beta=fV/fV0;
 828  if (fabs(fDv0)>0.)
 829  {
 830   Double_t vecv[3],errv[3],errb[3];
 831   Double_t sqerr=0;
 832   fV.GetVector(vecv,"car");
 833   fV.GetErrors(errv,"car");
 834   for (Int_t i=0; i<3; i++)
 835   {
 836    sqerr=pow((errv[i]/fV0),2)+pow((vecv[i]*fDv0/(fV0*fV0)),2);
 837    errb[i]=sqrt(sqerr);
 838   }
 839   beta.SetErrors(errb,"car");
 840  }
 841  return beta;
 842 }
 843 ///////////////////////////////////////////////////////////////////////////
 844 Double_t Ali4Vector::GetBeta()
 845 {
 846 // Provide the beta value (i.e. v/c) corresponding to this 4-vector.
 847  Ali3Vector beta=GetBetaVector();
 848  Double_t val=beta.GetNorm();
 849  fDresult=beta.GetResultError();
 850  return val;
 851 }
 852 ///////////////////////////////////////////////////////////////////////////
 853 Double_t Ali4Vector::GetGamma()
 854 {
 855 // Provide the Lorentz gamma factor corresponding to this 4-vector.
 856 // In case the gamma factor is infinite a value of -1 is returned.
 857  Double_t gamma=-1;
 858  fDresult=0;
 859  Double_t inv=sqrt(fabs(fV2));
 860  if (inv>0.)
 861  {
 862   Double_t dinv=fDv2/(2.*inv);
 863   gamma=fV0/inv;
 864   Double_t sqerr=pow((fDv0/inv),2)+pow((fV0*dinv/fV2),2);
 865   fDresult=sqrt(sqerr);
 866  }
 867  return gamma;
 868 }
 869 ///////////////////////////////////////////////////////////////////////////