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