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), reference frames (f) and angular units (u)
  89 // are specified via
  90 // SetVector(Float_t* v,TString f,TString u)
  91 // SetErrors(Float_t* e,TString f,TString u)
  92 // under the following conventions :
  93 //
  94 // f="car" ==> v in Cartesian coordinates   (x,y,z)
  95 // f="sph" ==> v in Spherical coordinates   (r,theta,phi)
  96 // f="cyl" ==> v in Cylindrical coordinates (rho,phi,z)
  97 //
  98 // u="rad" ==> angles in radians
  99 // u="deg" ==> angles in degrees
 100 //
 101 // The "f" and "u" facilities only serve as a convenient user interface.
 102 // Internally the actual storage of the various components is performed
 103 // in a unique way. This allows setting/retrieval of vector components in a
 104 // user selected frame/unit convention at any time.
 105 //
 106 // Example :
 107 // ---------
 108 //
 109 // Ali4Vector a;
 110 //
 111 // Float_t v[4]={25,-1,3,7};
 112 // a.SetVector(v,"car");
 113 //
 114 // Float_t vec[4];
 115 // a.GetVector(vec,"sph","deg");
 116 //
 117 // Ali4Vector b;
 118 // Float_t v2[4]={33,6,-18,2};
 119 // b.SetVector(v2,"car");
 120 //
 121 // Float_t dotpro=a.Dot(b);
 122 //
 123 // Float_t x0=16;
 124 // Ali3Vector x;
 125 // Float_t vec2[3]={1,2,3};
 126 // x.SetVector(vec2,"car");
 127 //
 128 // Ali4Vector c;
 129 // c.SetVector(x0,x);
 130 // c.GetVector(vec,"car");
 131 // c.Data("cyl");
 132 // c=a+b;
 133 // c=a-b;
 134 // c=a*5;
 135 //
 136 //--- Author: Nick van Eijndhoven 01-apr-1999 UU-SAP Utrecht
 137 //- Modified: NvE $Date$ UU-SAP Utrecht
 138 ///////////////////////////////////////////////////////////////////////////
 139
 140 #include "Ali4Vector.h"
 141 #include "Riostream.h"
 142
 143 ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
 144
 145 Ali4Vector::Ali4Vector()
 146 {
 147 // Creation of a contravariant 4-vector and initialisation of parameters.
 148 // All values are initialised to 0.
 149 // Scalar mode is initially selected.
 150  SetZero();
 151  fScalar=1;
 152 }
 153 ///////////////////////////////////////////////////////////////////////////
 154 Ali4Vector::~Ali4Vector()
 155 {
 156 // Destructor to delete dynamically allocated memory
 157 }
 158 ///////////////////////////////////////////////////////////////////////////
 159 Ali4Vector::Ali4Vector(const Ali4Vector& v)
 160 {
 161 // Copy constructor
 162  fScalar=v.fScalar;
 163  fV2=v.fV2;
 164  fDv2=v.fDv2;
 165  fV0=v.fV0;
 166  fDv0=v.fDv0;
 167  fDresult=v.fDresult;
 168  fV=v.fV;
 169 }
 170 ///////////////////////////////////////////////////////////////////////////
 171 void Ali4Vector::Load(Ali4Vector& q)
 172 {
 173 // Load all attributes of the input Ali4Vector into this Ali4Vector object.
 174  Int_t temp1=q.GetScalarFlag();
 175  Double_t temp2=q.GetResultError();
 176  Double_t a[4];
 177  q.GetVector(a,"sph");
 178  SetVector(a,"sph");
 179  q.GetErrors(a,"car");
 180  SetErrors(a,"car");
 181  fScalar=temp1;
 182  fDresult=temp2;
 183 }
 184 ///////////////////////////////////////////////////////////////////////////
 185 void Ali4Vector::SetZero()
 186 {
 187 // (Re)set all attributes to zero.
 188 // Note : The (de)selection of the scalar mode is not modified.
 189  fV2=0;
 190  fDv2=0;
 191  fV0=0;
 192  fDv0=0;
 193  fDresult=0;
 194  fV.SetZero();
 195 }
 196 ///////////////////////////////////////////////////////////////////////////
 197 void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v)
 198 {
 199 // Store contravariant vector.
 200 // The error on the scalar part is initialised to 0.
 201 // The errors on the vector part are taken from the input Ali3Vector.
 202 // Scalar mode is automatically selected.
 203 // The error on scalar result operations is reset to 0.
 204  fScalar=1;
 205  fV0=v0;
 206  fV=v;
 207  fV2=pow(v0,2)-fV.Dot(fV);
 208  SetScalarError(0);
 209 }
 210 ///////////////////////////////////////////////////////////////////////////
 211 void Ali4Vector::SetVector(Double_t* v,TString f,TString u)
 212 {
 213 // Store vector according to reference frame f.
 214 //
 215 // The string argument "u" allows to choose between different angular units
 216 // in case e.g. a spherical frame is selected.
 217 // u = "rad" : angles provided in radians
 218 //     "deg" : angles provided in degrees
 219 //
 220 // The default is u="rad".
 221 //
 222 // All errors are initialised to 0.
 223 //
 224 // Scalar mode is automatically selected.
 225 // The error on scalar result operations is reset to 0.
 226
 227  fScalar=1;
 228  Double_t a[3];
 229  for (Int_t i=0; i<3; i++)
 230  {
 231   a[i]=v[i+1];
 232  }
 233  fV0=v[0];
 234  fV.SetVector(a,f,u);
 235  fV2=pow(fV0,2)-fV.Dot(fV);
 236  fDv2=0;
 237  fDv0=0;
 238  fDresult=0;
 239 }
 240 ///////////////////////////////////////////////////////////////////////////
 241 void Ali4Vector::GetVector(Double_t* v,TString f,TString u)
 242 {
 243 // Provide 4-vector components according to reference frame f
 244 // and according to the current mode.
 245 // Scalar mode    : The scalar part is directly returned via v[0].
 246 // Invariant mode : The scalar part is re-calculated via the value
 247 //                  of the Lorentz invariant and then returned via v[0].
 248 //
 249 // The string argument "u" allows to choose between different angular units
 250 // in case e.g. a spherical frame is selected.
 251 // u = "rad" : angles provided in radians
 252 //     "deg" : angles provided in degrees
 253 //
 254 // The default is u="rad".
 255
 256  if (fScalar)
 257  {
 258   v[0]=fV0;
 259  }
 260  else
 261  {
 262   v[0]=sqrt(fV.Dot(fV)+fV2);
 263  }
 264  Double_t a[3];
 265  fV.GetVector(a,f,u);
 266  for (Int_t i=0; i<3; i++)
 267  {
 268   v[i+1]=a[i];
 269  }
 270 }
 271 ///////////////////////////////////////////////////////////////////////////
 272 void Ali4Vector::SetVector(Float_t* v,TString f,TString u)
 273 {
 274 // Store vector according to reference frame f.
 275 //
 276 // The string argument "u" allows to choose between different angular units
 277 // in case e.g. a spherical frame is selected.
 278 // u = "rad" : angles provided in radians
 279 //     "deg" : angles provided in degrees
 280 //
 281 // The default is u="rad".
 282 //
 283 // All errors are initialised to 0.
 284 //
 285 // Scalar mode is automatically selected.
 286 // The error on scalar result operations is reset to 0.
 287
 288  Double_t vec[4];
 289  for (Int_t i=0; i<4; i++)
 290  {
 291   vec[i]=v[i];
 292  }
 293  SetVector(vec,f,u);
 294 }
 295 ///////////////////////////////////////////////////////////////////////////
 296 void Ali4Vector::GetVector(Float_t* v,TString f,TString u)
 297 {
 298 // Provide 4-vector components according to reference frame f
 299 // and according to the current mode.
 300 // Scalar mode    : The scalar part is directly returned via v[0].
 301 // Invariant mode : The scalar part is re-calculated via the value
 302 //                  of the Lorentz invariant and then returned via v[0].
 303 //
 304 // The string argument "u" allows to choose between different angular units
 305 // in case e.g. a spherical frame is selected.
 306 // u = "rad" : angles provided in radians
 307 //     "deg" : angles provided in degrees
 308 //
 309 // The default is u="rad".
 310
 311  Double_t vec[4];
 312  GetVector(vec,f,u);
 313  for (Int_t i=0; i<4; i++)
 314  {
 315   v[i]=vec[i];
 316  }
 317 }
 318 ///////////////////////////////////////////////////////////////////////////
 319 Double_t Ali4Vector::GetScalar()
 320 {
 321 // Provide the scalar part.
 322 // The error on the scalar value is available via GetResultError()
 323 // after invokation of GetScalar().
 324  if (fScalar)
 325  {
 326   fDresult=fDv0;
 327   return fV0;
 328  }
 329  else
 330  {
 331   Double_t dot=fV.Dot(fV);
 332   Double_t ddot=fV.GetResultError();
 333   Double_t v02=dot+fV2;
 334   Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2));
 335   Double_t v0=sqrt(fabs(v02));
 336   Double_t dv0=0;
 337   if (v0) dv0=dv02/(2.*v0);
 338   fDresult=dv0;
 339   return v0;
 340  }
 341 }
 342 ///////////////////////////////////////////////////////////////////////////
 343 Double_t Ali4Vector::GetResultError() const
 344 {
 345 // Provide the error on the result of an operation yielding a scalar
 346 // E.g. GetScalar(), GetInvariant() or Dot()
 347  return fDresult;
 348 }
 349 ///////////////////////////////////////////////////////////////////////////
 350 void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
 351 {
 352 // Modify the scalar part (v0) and its error (dv0).
 353 // The default value for dv0 is 0.
 354 // The vector part is not modified.
 355 // Scalar mode is automatically selected
 356 // ==> Lorentz invariant and its error are updated.
 357 // The error on scalar result operations is reset to 0.
 358  fScalar=1;
 359  fV0=v0;
 360  fV2=pow(v0,2)-fV.Dot(fV);
 361  SetScalarError(dv0);
 362 }
 363 ///////////////////////////////////////////////////////////////////////////
 364 void Ali4Vector::SetScalarError(Double_t dv0)
 365 {
 366 // Set the error on the scalar part.
 367 // If in scalar mode, update error on the invariant accordingly.
 368 // The error on scalar result operations is reset to 0.
 369  fDv0=dv0;
 370  if (fScalar)
 371  {
 372   Double_t norm=fV.GetNorm();
 373   Double_t dnorm=fV.GetResultError();
 374   fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2));
 375  }
 376  fDresult=0;
 377 }
 378 ///////////////////////////////////////////////////////////////////////////
 379 void Ali4Vector::Set3Vector(Ali3Vector& v)
 380 {
 381 // Set the 3-vector part, the errors are taken from the input Ali3Vector
 382 // Scalar mode    : The scalar part and its error are not modified,
 383 //                  the Lorentz invariantand its error are re-calculated.
 384 // Invariant mode : The Lorentz invariant and its error are not modified,
 385 //                  the scalar part and its error are re-calculated.
 386 // The error on scalar result operations is reset to 0.
 387  fV=v;
 388  if (fScalar)
 389  {
 390   SetScalar(fV0,fDv0);
 391  }
 392  else
 393  {
 394   SetInvariant(fV2,fDv2);
 395  }
 396 }
 397 ///////////////////////////////////////////////////////////////////////////
 398 void Ali4Vector::Set3Vector(Double_t* v,TString f,TString u)
 399 {
 400 // Set the 3-vector part according to reference frame f
 401 //
 402 // The string argument "u" allows to choose between different angular units
 403 // in case e.g. a spherical frame is selected.
 404 // u = "rad" : angles provided in radians
 405 //     "deg" : angles provided in degrees
 406 //
 407 // The default is u="rad".
 408 //
 409 // The errors on the vector part are initialised to 0
 410 //
 411 // Scalar mode    : The scalar part and its error are not modified,
 412 //                  the Lorentz invariantand its error are re-calculated.
 413 // Invariant mode : The Lorentz invariant and its error are not modified,
 414 //                  the scalar part and its error are re-calculated.
 415 // The error on scalar result operations is reset to 0.
 416
 417  Double_t a[3];
 418  for (Int_t i=0; i<3; i++)
 419  {
 420   a[i]=v[i];
 421  }
 422  fV.SetVector(a,f,u);
 423
 424  if (fScalar)
 425  {
 426   SetScalar(fV0,fDv0);
 427  }
 428  else
 429  {
 430   SetInvariant(fV2,fDv2);
 431  }
 432 }
 433 ///////////////////////////////////////////////////////////////////////////
 434 void Ali4Vector::Set3Vector(Float_t* v,TString f,TString u)
 435 {
 436 // Set the 3-vector part according to reference frame f
 437 //
 438 // The string argument "u" allows to choose between different angular units
 439 // in case e.g. a spherical frame is selected.
 440 // u = "rad" : angles provided in radians
 441 //     "deg" : angles provided in degrees
 442 //
 443 // The default is u="rad".
 444 //
 445 // The errors on the vector part are initialised to 0
 446 // The Lorentz invariant is not modified
 447 // The error on scalar result operations is reset to 0.
 448
 449  Double_t vec[3];
 450  for (Int_t i=0; i<3; i++)
 451  {
 452   vec[i]=v[i];
 453  }
 454  Set3Vector(vec,f,u);
 455 }
 456 ///////////////////////////////////////////////////////////////////////////
 457 void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
 458 {
 459 // Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2).
 460 // The default value for the error dv2 is 0.
 461 // The vector part is not modified.
 462 // Invariant mode is automatically selected
 463 // ==> the scalar part and its error are updated.
 464 // The error on scalar result operations is reset to 0.
 465 //
 466  fScalar=0;
 467  fV2=v2;
 468  fDv2=dv2;
 469  fV0=GetScalar();
 470  fDv0=GetResultError();
 471  fDresult=0;
 472 }
 473 ///////////////////////////////////////////////////////////////////////////
 474 void Ali4Vector::SetInvariantError(Double_t dv2)
 475 {
 476 // Set the error on the Lorentz invariant.
 477 // If in invariant mode, update error on the scalar part accordingly.
 478 // The error on scalar result operations is reset to 0.
 479  fDv2=dv2;
 480  if (!fScalar)
 481  {
 482   fDv0=GetResultError();
 483  }
 484  fDresult=0;
 485 }
 486 ///////////////////////////////////////////////////////////////////////////
 487 Double_t Ali4Vector::GetInvariant()
 488 {
 489 // Provide the Lorentz invariant v^i*v_i.
 490 // The error on the Lorentz invariant is available via GetResultError()
 491 // after invokation of GetInvariant().
 492  if (!fScalar)
 493  {
 494   fDresult=fDv2;
 495   return fV2;
 496  }
 497  else
 498  {
 499   Double_t inv=Dot(*this);
 500   return inv;
 501  }
 502 }
 503 ///////////////////////////////////////////////////////////////////////////
 504 Ali3Vector Ali4Vector::Get3Vector() const
 505 {
 506 // Provide the 3-vector part
 507  return fV;
 508 }
 509 ///////////////////////////////////////////////////////////////////////////
 510 void Ali4Vector::SetErrors(Double_t* e,TString f,TString u)
 511 {
 512 // Store errors for vector v^i according to reference frame f
 513 //
 514 // The string argument "u" allows to choose between different angular units
 515 // in case e.g. a spherical frame is selected.
 516 // u = "rad" : angles provided in radians
 517 //     "deg" : angles provided in degrees
 518 //
 519 // The default is u="rad".
 520 //
 521 // If in scalar mode, update error on the invariant accordingly.
 522 // The error on scalar result operations is reset to 0.
 523
 524  Double_t a[3];
 525  for (Int_t i=0; i<3; i++)
 526  {
 527   a[i]=e[i+1];
 528  }
 529  SetScalarError(e[0]);
 530  fV.SetErrors(a,f,u);
 531 }
 532 ///////////////////////////////////////////////////////////////////////////
 533 void Ali4Vector::SetErrors(Float_t* e,TString f,TString u)
 534 {
 535 // Store errors for vector v^i according to reference frame f
 536 //
 537 // The string argument "u" allows to choose between different angular units
 538 // in case e.g. a spherical frame is selected.
 539 // u = "rad" : angles provided in radians
 540 //     "deg" : angles provided in degrees
 541 //
 542 // The default is u="rad".
 543 //
 544 // If in scalar mode, update error on the invariant accordingly.
 545 // The error on scalar result operations is reset to 0.
 546
 547  Double_t a[4];
 548  for (Int_t i=0; i<4; i++)
 549  {
 550   a[i]=e[i];
 551  }
 552  SetErrors(a,f,u);
 553 }
 554 ///////////////////////////////////////////////////////////////////////////
 555 void Ali4Vector::GetErrors(Double_t* e,TString f,TString u)
 556 {
 557 // Provide errors for vector v^i according to reference frame f
 558 // and according to the current mode.
 559 // Scalar mode    : The error on the scalar part is directly returned via e[0].
 560 // Invariant mode : The error on the scalar part is re-calculated via the error
 561 //                  value on the Lorentz invariant and then returned via e[0].
 562 //
 563 // The string argument "u" allows to choose between different angular units
 564 // in case e.g. a spherical frame is selected.
 565 // u = "rad" : angles provided in radians
 566 //     "deg" : angles provided in degrees
 567 //
 568 // The default is u="rad".
 569
 570  Double_t a[3];
 571  fV.GetErrors(a,f,u);
 572
 573  // Dummy invokation of GetScalar to obtain automatic proper error determination
 574  e[0]=GetScalar();
 575
 576  e[0]=GetResultError();
 577
 578  for (Int_t i=0; i<3; i++)
 579  {
 580   e[i+1]=a[i];
 581  }
 582 }
 583 ///////////////////////////////////////////////////////////////////////////
 584 void Ali4Vector::GetErrors(Float_t* e,TString f,TString u)
 585 {
 586 // Provide errors for vector v^i according to reference frame f
 587 // and according to the current mode.
 588 // Scalar mode    : The error on the scalar part is directly returned via e[0].
 589 // Invariant mode : The error on the scalar part is re-calculated via the error
 590 //                  value on the Lorentz invariant and then returned via e[0].
 591 //
 592 // The string argument "u" allows to choose between different angular units
 593 // in case e.g. a spherical frame is selected.
 594 // u = "rad" : angles provided in radians
 595 //     "deg" : angles provided in degrees
 596 //
 597 // The default is u="rad".
 598
 599  Double_t a[4];
 600  GetErrors(a,f,u);
 601  for (Int_t i=0; i<4; i++)
 602  {
 603   e[i]=a[i];
 604  }
 605 }
 606 ///////////////////////////////////////////////////////////////////////////
 607 void Ali4Vector::Data(TString f,TString u)
 608 {
 609 // Print contravariant vector components and errors according to
 610 // reference frame f and according to the current mode.
 611 // Scalar mode    : The scalar part and its error are directly returned.
 612 // Invariant mode : The scalar part and its error are re-calculated via the
 613 //                  value (and error) of the Lorentz invariant.
 614 //
 615 // The string argument "u" allows to choose between different angular units
 616 // in case e.g. a spherical frame is selected.
 617 // u = "rad" : angles provided in radians
 618 //     "deg" : angles provided in degrees
 619 //
 620 // The defaults are f="car" and u="rad".
 621
 622  if (f=="car" || f=="sph" || f=="cyl")
 623  {
 624   Double_t vec[4],err[4];
 625   GetVector(vec,f,u);
 626   GetErrors(err,f,u);
 627   Double_t inv=GetInvariant();
 628   Double_t dinv=GetResultError();
 629   cout << " Contravariant vector in " << f.Data() << " (" << u.Data() << ") coordinates : "
 630        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl;
 631   cout << " ------------- Errors in " << f.Data() << " (" << u.Data() << ") coordinates : "
 632        << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl;
 633   cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
 634  }
 635  else
 636  {
 637   cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl
 638        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl;
 639  }
 640 }
 641 ///////////////////////////////////////////////////////////////////////////
 642 Double_t Ali4Vector::Dot(Ali4Vector& q)
 643 {
 644 // Provide the dot product of the current vector with vector q
 645  Double_t dotpro=0;
 646  Double_t a0=GetScalar();
 647  Double_t da0=GetResultError();
 648  if ((this) == &q) // Check for special case v.Dot(v)
 649  {
 650   Double_t norm=fV.GetNorm();
 651   Double_t dnorm=fV.GetResultError();
 652   dotpro=pow(a0,2)-pow(norm,2);
 653   fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2));
 654  }
 655  else
 656  {
 657   Double_t b0=q.GetScalar();
 658   Double_t db0=q.GetResultError();
 659   Ali3Vector b=q.Get3Vector();
 660
 661   Double_t dot=fV.Dot(b);
 662   Double_t ddot=fV.GetResultError();
 663
 664   dotpro=a0*b0-dot;
 665
 666   fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
 667  }
 668
 669  return dotpro;
 670 }
 671 ///////////////////////////////////////////////////////////////////////////
 672 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
 673 {
 674 // Add 4-vector q to the current 4-vector
 675 // Error propagation is performed automatically
 676  Double_t a0=GetScalar();
 677  Double_t da0=GetResultError();
 678  Ali3Vector a=Get3Vector();
 679  Double_t b0=q.GetScalar();
 680  Double_t db0=q.GetResultError();
 681  Ali3Vector b=q.Get3Vector();
 682
 683  Double_t c0=a0+b0;
 684  Ali3Vector c=a+b;
 685  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 686
 687  Ali4Vector v;
 688  v.SetVector(c0,c);
 689  v.SetScalarError(dc0);
 690  return v;
 691 }
 692 ///////////////////////////////////////////////////////////////////////////
 693 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
 694 {
 695 // Subtract 4-vector q from the current 4-vector
 696 // Error propagation is performed automatically
 697  Double_t a0=GetScalar();
 698  Double_t da0=GetResultError();
 699  Ali3Vector a=Get3Vector();
 700  Double_t b0=q.GetScalar();
 701  Double_t db0=q.GetResultError();
 702  Ali3Vector b=q.Get3Vector();
 703
 704  Double_t c0=a0-b0;
 705  Ali3Vector c=a-b;
 706  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 707
 708  Ali4Vector v;
 709  v.SetVector(c0,c);
 710  v.SetScalarError(dc0);
 711  return v;
 712 }
 713 ///////////////////////////////////////////////////////////////////////////
 714 Ali4Vector Ali4Vector::operator*(Double_t s)
 715 {
 716 // Multiply the current 4-vector with a scalar s
 717 // Error propagation is performed automatically
 718  Double_t a0=GetScalar();
 719  Double_t da0=GetResultError();
 720  Ali3Vector a=Get3Vector();
 721
 722  a0*=s;
 723  a*=s;
 724  da0*=s;
 725
 726  Ali4Vector v;
 727  v.SetVector(a0,a);
 728  v.SetScalarError(da0);
 729
 730  return v;
 731 }
 732 ///////////////////////////////////////////////////////////////////////////
 733 Ali4Vector Ali4Vector::operator/(Double_t s)
 734 {
 735 // Divide the current vector by a scalar s
 736 // Error propagation is performed automatically
 737
 738  if (fabs(s)<1.e-20) // Protect against division by 0
 739  {
 740   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
 741   return *this;
 742  }
 743  else
 744  {
 745   Double_t a0=GetScalar();
 746   Double_t da0=GetResultError();
 747   Ali3Vector a=Get3Vector();
 748
 749   a0/=s;
 750   a/=s;
 751   da0/=s;
 752
 753   Ali4Vector v;
 754   v.SetVector(a0,a);
 755   v.SetScalarError(da0);
 756
 757   return v;
 758  }
 759 }
 760 ///////////////////////////////////////////////////////////////////////////
 761 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
 762 {
 763 // Add 4-vector q to the current 4-vector
 764 // Error propagation is performed automatically
 765  Double_t a0=GetScalar();
 766  Double_t da0=GetResultError();
 767  Ali3Vector a=Get3Vector();
 768  Double_t b0=q.GetScalar();
 769  Double_t db0=q.GetResultError();
 770  Ali3Vector b=q.Get3Vector();
 771
 772  Double_t c0=a0+b0;
 773  Ali3Vector c=a+b;
 774  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 775
 776  SetVector(c0,c);
 777  SetScalarError(dc0);
 778
 779  return *this;
 780 }
 781 ///////////////////////////////////////////////////////////////////////////
 782 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
 783 {
 784 // Subtract 4-vector q from the current 4-vector
 785 // Error propagation is performed automatically
 786  Double_t a0=GetScalar();
 787  Double_t da0=GetResultError();
 788  Ali3Vector a=Get3Vector();
 789  Double_t b0=q.GetScalar();
 790  Double_t db0=q.GetResultError();
 791  Ali3Vector b=q.Get3Vector();
 792
 793  Double_t c0=a0-b0;
 794  Ali3Vector c=a-b;
 795  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
 796
 797  SetVector(c0,c);
 798  SetScalarError(dc0);
 799
 800  return *this;
 801 }
 802 ///////////////////////////////////////////////////////////////////////////
 803 Ali4Vector& Ali4Vector::operator*=(Double_t s)
 804 {
 805 // Multiply the current 4-vector with a scalar s
 806 // Error propagation is performed automatically
 807  Double_t a0=GetScalar();
 808  Double_t da0=GetResultError();
 809  Ali3Vector a=Get3Vector();
 810
 811  a0*=s;
 812  a*=s;
 813  da0*=s;
 814
 815  SetVector(a0,a);
 816  SetScalarError(da0);
 817
 818  return *this;
 819 }
 820 ///////////////////////////////////////////////////////////////////////////
 821 Ali4Vector& Ali4Vector::operator/=(Double_t s)
 822 {
 823 // Divide the current vector by a scalar s
 824 // Error propagation is performed automatically
 825
 826  if (fabs(s)<1.e-20) // Protect against division by 0
 827  {
 828   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
 829   return *this;
 830  }
 831  else
 832  {
 833   Double_t a0=GetScalar();
 834   Double_t da0=GetResultError();
 835   Ali3Vector a=Get3Vector();
 836
 837   a0/=s;
 838   a/=s;
 839   da0/=s;
 840
 841   SetVector(a0,a);
 842   SetScalarError(da0);
 843
 844   return *this;
 845  }
 846 }
 847 ///////////////////////////////////////////////////////////////////////////
 848 Int_t Ali4Vector::GetScalarFlag() const
 849 {
 850 // Provide the value of the fScalar flag (for internal use only).
 851  return fScalar;
 852 }
 853 ///////////////////////////////////////////////////////////////////////////
 854 Ali3Vector Ali4Vector::GetVecTrans() const
 855 {
 856 // Provide the transverse vector part w.r.t. z-axis.
 857 // Error propagation is performed automatically
 858
 859  return fV.GetVecTrans();
 860 }
 861 ///////////////////////////////////////////////////////////////////////////
 862 Ali3Vector Ali4Vector::GetVecLong() const
 863 {
 864 // Provide the longitudinal vector part w.r.t. z-axis.
 865 // Error propagation is performed automatically
 866
 867  return fV.GetVecLong();
 868 }
 869 ///////////////////////////////////////////////////////////////////////////
 870 Double_t Ali4Vector::GetScaTrans()
 871 {
 872 // Provide the "transverse value" of the scalar part w.r.t. z-axis.
 873 // This provides a basis for e.g. E_trans calculation.
 874 // Note : the returned value is always positive or zero.
 875 // The error on the value is available via GetResultError()
 876 // after invokation of GetScaTrans().
 877  Double_t a[3],ea[3];
 878
 879  fV.GetVector(a,"sph");
 880  fV.GetErrors(ea,"sph");
 881
 882  Double_t s=GetScalar();
 883  Double_t ds=GetResultError();
 884
 885  Double_t st,dst2;
 886  st=s*sin(a[1]);
 887  dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2);
 888
 889  fDresult=sqrt(dst2);
 890  return fabs(st);
 891 }
 892 ///////////////////////////////////////////////////////////////////////////
 893 Double_t Ali4Vector::GetScaLong()
 894 {
 895 // Provide the "longitudinal value" of the scalar part w.r.t. z-axis.
 896 // This provides a basis for e.g. E_long calculation.
 897 // Note : the returned value can also be negative.
 898 // The error on the value is available via GetResultError()
 899 // after invokation of GetScaLong().
 900  Double_t a[3],ea[3];
 901
 902  fV.GetVector(a,"sph");
 903  fV.GetErrors(ea,"sph");
 904
 905  Double_t s=GetScalar();
 906  Double_t ds=GetResultError();
 907
 908  Double_t sl,dsl2;
 909  sl=s*cos(a[1]);
 910  dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2);
 911
 912  fDresult=sqrt(dsl2);
 913  return sl;
 914 }
 915 ///////////////////////////////////////////////////////////////////////////
 916 Double_t Ali4Vector::GetPseudoRapidity()
 917 {
 918 // Provide the pseudorapidity value of the vector part w.r.t. z-axis.
 919 // The error on the value is available via GetResultError()
 920 // after invokation of GetPseudoRapidity().
 921  Double_t eta=fV.GetPseudoRapidity();
 922  fDresult=fV.GetResultError();
 923  return eta;
 924 }
 925 ///////////////////////////////////////////////////////////////////////////
 926 Ali3Vector Ali4Vector::GetBetaVector() const
 927 {
 928 // Provide the beta 3-vector corresponding to this 4-vector.
 929  Ali3Vector beta;
 930  if (fabs(fV0)>0.) beta=fV/fV0;
 931  if (fabs(fDv0)>0.)
 932  {
 933   Double_t vecv[3],errv[3],errb[3];
 934   Double_t sqerr=0;
 935   fV.GetVector(vecv,"car");
 936   fV.GetErrors(errv,"car");
 937   for (Int_t i=0; i<3; i++)
 938   {
 939    sqerr=pow((errv[i]/fV0),2)+pow((vecv[i]*fDv0/(fV0*fV0)),2);
 940    errb[i]=sqrt(sqerr);
 941   }
 942   beta.SetErrors(errb,"car");
 943  }
 944  return beta;
 945 }
 946 ///////////////////////////////////////////////////////////////////////////
 947 Double_t Ali4Vector::GetBeta()
 948 {
 949 // Provide the beta value (i.e. v/c) corresponding to this 4-vector.
 950  Ali3Vector beta=GetBetaVector();
 951  Double_t val=beta.GetNorm();
 952  fDresult=beta.GetResultError();
 953  return val;
 954 }
 955 ///////////////////////////////////////////////////////////////////////////
 956 Double_t Ali4Vector::GetGamma()
 957 {
 958 // Provide the Lorentz gamma factor corresponding to this 4-vector.
 959 // In case the gamma factor is infinite a value of -1 is returned.
 960  Double_t gamma=-1;
 961  fDresult=0;
 962  Double_t inv=sqrt(fabs(fV2));
 963  if (inv>0.)
 964  {
 965   Double_t dinv=fDv2/(2.*inv);
 966   gamma=fV0/inv;
 967   Double_t sqerr=pow((fDv0/inv),2)+pow((fV0*dinv/fV2),2);
 968   fDresult=sqrt(sqerr);
 969  }
 970  return gamma;
 971 }
 972 ///////////////////////////////////////////////////////////////////////////
 973 Double_t Ali4Vector::GetX(Int_t i,TString f,TString u)
 974 {
 975 // Provide i-th vector component according to reference frame f.
 976 //
 977 // The string argument "u" allows to choose between different angular units
 978 // in case e.g. a spherical frame is selected.
 979 // u = "rad" : angles provided in radians
 980 //     "deg" : angles provided in degrees
 981 //
 982 // The default is u="rad".
 983 //
 984 // The vector components are addressed via the generic x0,x1,x2,x3 notation.
 985 // So, i=0 denotes the scalar component and i=1 denotes the first 3-vector component.
 986 // The error on the selected component can be obtained via the
 987 // usual GetResultError() facility.
 988
 989  fDresult=0;
 990
 991  if (i<0 || i>3) return 0;
 992
 993  Double_t x=0;
 994
 995  if (i==0)
 996  {
 997   x=GetScalar();
 998  }
 999  else
1000  {
1001   x=fV.GetX(i,f,u);
1002   fDresult=fV.GetResultError();
1003  }
1004
1005  return x;
1006 }
1007 ///////////////////////////////////////////////////////////////////////////
1008 Double_t Ali4Vector::GetOpeningAngle(Ali4Vector& q,TString u)
1009 {
1010 // Provide the opening angle between 3-vector parts with 4-vector q.
1011 // The string argument "u" allows to choose between different output units.
1012 // u = "rad" : opening angle provided in radians
1013 //     "deg" : opening angle provided in degrees
1014 //
1015 // The default is u="rad".
1016
1017  Ali3Vector v1=fV;
1018  Ali3Vector v2=q.Get3Vector();
1019
1020  Double_t ang=v1.GetOpeningAngle(v2,u);
1021  fDresult=v1.GetResultError();
1022
1023  return ang;
1024 }
1025 ///////////////////////////////////////////////////////////////////////////
1026 Double_t Ali4Vector::GetOpeningAngle(Ali3Vector& q,TString u)
1027 {
1028 // Provide the opening angle between the 3-vector part and 3-vector q.
1029 // The string argument "u" allows to choose between different output units.
1030 // u = "rad" : opening angle provided in radians
1031 //     "deg" : opening angle provided in degrees
1032 //
1033 // The default is u="rad".
1034
1035  Ali3Vector v1=fV;
1036
1037  Double_t ang=v1.GetOpeningAngle(q,u);
1038  fDresult=v1.GetResultError();
1039
1040  return ang;
1041 }
1042 ///////////////////////////////////////////////////////////////////////////