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