RALICE/Ali3Vector.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 Ali3Vector
  20 // Handling of 3-vectors in various reference frames.
  21 //
  22 // This class is meant to serve as a base class for ALICE objects
  23 // that have 3-dimensional vector characteristics.
  24 // Error propagation is performed automatically.
  25 //
  26 // Note :
  27 // ------
  28 // Vectors (v), Errors (e) and reference frames (f) are specified via
  29 // SetVector(Float_t* v,TString f)
  30 // SetErrors(Float_t* e,TString f)
  31 // under the following conventions :
  32 //
  33 // f="car" ==> v in Cartesian coordinates   (x,y,z)
  34 // f="sph" ==> v in Spherical coordinates   (r,theta,phi)
  35 // f="cyl" ==> v in Cylindrical coordinates (rho,phi,z)
  36 //
  37 // All angles are in radians.
  38 //
  39 // Example :
  40 // ---------
  41 //
  42 // Ali3Vector a;
  43 // Float_t v[3]={-1,25,7};
  44 // Float_t e[3]={0.03,0.5,0.21};
  45 // a.SetVector(v,"car");
  46 // a.SetErrors(e,"car");
  47 // a.Data();
  48 //
  49 // Float_t vec[3];
  50 // Float_t err[3];
  51 // a.GetVector(vec,"sph");
  52 // a.GetErrors(vec,"sph");
  53 //
  54 // Ali3Vector b;
  55 // Float_t v2[3]={6,-18,33};
  56 // Float_t e2[3]={0.19,0.45,0.93};
  57 // b.SetVector(v2,"car");
  58 // b.SetErrors(e2,"car");
  59 //
  60 // Float_t dotpro=a.Dot(b);
  61 // Float_t doterror=a.GetResultError();
  62 //
  63 // Ali3Vector c=a.Cross(b);
  64 // c.Data("sph");
  65 // c.GetVector(vec,"cyl");
  66 // c.GetErrors(err,"cyl");
  67 //
  68 // Float_t norm=c.GetNorm();
  69 // Float_t normerror=c.GetResultError();
  70 //
  71 // c=a+b;
  72 // c=a-b;
  73 // c=a*5;
  74 //
  75 //--- Author: Nick van Eijndhoven 30-mar-1999 UU-SAP Utrecht
  76 //- Modified: NvE $Date$ UU-SAP Utrecht
  77 ///////////////////////////////////////////////////////////////////////////
  78
  79 #include "Ali3Vector.h"
  80 #include "Riostream.h"
  81
  82 ClassImp(Ali3Vector) // Class implementation to enable ROOT I/O
  83
  84 Ali3Vector::Ali3Vector()
  85 {
  86 // Creation of an Ali3Vector object and initialisation of parameters
  87 // All attributes initialised to 0
  88  SetZero();
  89 }
  90 ///////////////////////////////////////////////////////////////////////////
  91 Ali3Vector::~Ali3Vector()
  92 {
  93 // Destructor to delete dynamically allocated memory
  94 }
  95 ///////////////////////////////////////////////////////////////////////////
  96 Ali3Vector::Ali3Vector(const Ali3Vector& v)
  97 {
  98 // Copy constructor
  99  fV=v.fV;
 100  fTheta=v.fTheta;
 101  fPhi=v.fPhi;
 102  fDx=v.fDx;
 103  fDy=v.fDy;
 104  fDz=v.fDz;
 105  fDresult=v.fDresult;
 106 }
 107 ///////////////////////////////////////////////////////////////////////////
 108 void Ali3Vector::Load(Ali3Vector& q)
 109 {
 110 // Load all attributes of the input Ali3Vector into this Ali3Vector object.
 111  Double_t temp=q.GetResultError();
 112  Double_t a[3];
 113  q.GetVector(a,"sph");
 114  SetVector(a,"sph");
 115  q.GetErrors(a,"car");
 116  SetErrors(a,"car");
 117  fDresult=temp;
 118 }
 119 ///////////////////////////////////////////////////////////////////////////
 120 void Ali3Vector::SetZero()
 121 {
 122 // (Re)set all attributes to zero.
 123  fV=0;
 124  fTheta=0;
 125  fPhi=0;
 126  fDx=0;
 127  fDy=0;
 128  fDz=0;
 129  fDresult=0;
 130 }
 131 ///////////////////////////////////////////////////////////////////////////
 132 void Ali3Vector::SetVector(Double_t* v,TString f)
 133 {
 134 // Store vector according to reference frame f
 135 // All errors will be reset to 0
 136  fDx=0;
 137  fDy=0;
 138  fDz=0;
 139  fDresult=0;
 140
 141  Double_t pi=acos(-1.);
 142
 143  Int_t frame=0;
 144  if (f == "car") frame=1;
 145  if (f == "sph") frame=2;
 146  if (f == "cyl") frame=3;
 147
 148  Double_t x,y,z,rho,phi;
 149
 150  switch (frame)
 151  {
 152   case 1: // Cartesian coordinates
 153    x=v[0];
 154    y=v[1];
 155    z=v[2];
 156    fV=sqrt(x*x+y*y+z*z);
 157    fTheta=0;
 158    if (fV && fabs(z/fV)<=1.)
 159    {
 160     fTheta=acos(z/fV);
 161    }
 162    else
 163    {
 164     if (z<0.) fTheta=pi;
 165    }
 166    if (fTheta<0.) fTheta+=2.*pi;
 167    fPhi=0;
 168    if (x || y) fPhi=atan2(y,x);
 169    if (fPhi<0.) fPhi+=2.*pi;
 170    break;
 171
 172   case 2: // Spherical coordinates
 173    fV=v[0];
 174    fTheta=v[1];
 175    fPhi=v[2];
 176    break;
 177
 178   case 3: // Cylindrical coordinates
 179    rho=v[0];
 180    phi=v[1];
 181    z=v[2];
 182    fV=sqrt(rho*rho+z*z);
 183    fPhi=phi;
 184    if (fPhi<0.) fPhi+=2.*pi;
 185    fTheta=0;
 186    if (fV && fabs(z/fV)<=1.)
 187    {
 188     fTheta=acos(z/fV);
 189    }
 190    else
 191    {
 192     if (z<0.) fTheta=pi;
 193    }
 194    if (fTheta<0.) fTheta+=2.*pi;
 195    break;
 196
 197   default: // Unsupported reference frame
 198    cout << "*Ali3Vector::SetVector* Unsupported frame : " << f.Data() << endl
 199         << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
 200    fV=0;
 201    fTheta=0;
 202    fPhi=0;
 203    break;
 204  }
 205 }
 206 ///////////////////////////////////////////////////////////////////////////
 207 void Ali3Vector::GetVector(Double_t* v,TString f) const
 208 {
 209 // Provide vector according to reference frame f
 210  Int_t frame=0;
 211  if (f == "car") frame=1;
 212  if (f == "sph") frame=2;
 213  if (f == "cyl") frame=3;
 214
 215  switch (frame)
 216  {
 217   case 1: // Cartesian coordinates
 218    v[0]=fV*sin(fTheta)*cos(fPhi);
 219    v[1]=fV*sin(fTheta)*sin(fPhi);
 220    v[2]=fV*cos(fTheta);
 221    break;
 222
 223   case 2: // Spherical coordinates
 224    v[0]=fV;
 225    v[1]=fTheta;
 226    v[2]=fPhi;
 227    break;
 228
 229   case 3: // Cylindrical coordinates
 230    v[0]=fV*sin(fTheta);
 231    v[1]=fPhi;
 232    v[2]=fV*cos(fTheta);
 233    break;
 234
 235   default: // Unsupported reference frame
 236    cout << "*Ali3Vector::GetVector* Unsupported frame : " << f.Data() << endl
 237         << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
 238    for (Int_t i=0; i<3; i++)
 239    {
 240     v[i]=0;
 241    }
 242    break;
 243  }
 244 }
 245 ///////////////////////////////////////////////////////////////////////////
 246 void Ali3Vector::SetVector(Float_t* v,TString f)
 247 {
 248 // Store vector according to reference frame f
 249 // All errors will be reset to 0
 250  Double_t vec[3];
 251  for (Int_t i=0; i<3; i++)
 252  {
 253   vec[i]=v[i];
 254  }
 255  SetVector(vec,f);
 256 }
 257 ///////////////////////////////////////////////////////////////////////////
 258 void Ali3Vector::GetVector(Float_t* v,TString f) const
 259 {
 260 // Provide vector according to reference frame f
 261  Double_t vec[3];
 262  GetVector(vec,f);
 263  for (Int_t i=0; i<3; i++)
 264  {
 265   v[i]=vec[i];
 266  }
 267 }
 268 ///////////////////////////////////////////////////////////////////////////
 269 void Ali3Vector::SetErrors(Double_t* e,TString f)
 270 {
 271 // Store errors according to reference frame f
 272 // The error on scalar results is reset to 0
 273  fDresult=0;
 274
 275  Int_t frame=0;
 276  if (f == "car") frame=1;
 277  if (f == "sph") frame=2;
 278  if (f == "cyl") frame=3;
 279
 280  Double_t dx2,dy2,dz2,rho;
 281
 282  switch (frame)
 283  {
 284   case 1: // Cartesian coordinates
 285    fDx=fabs(e[0]);
 286    fDy=fabs(e[1]);
 287    fDz=fabs(e[2]);
 288    break;
 289
 290   case 2: // Spherical coordinates
 291    dx2=pow((cos(fPhi)*sin(fTheta)*e[0]),2)+pow((fV*cos(fTheta)*cos(fPhi)*e[1]),2)
 292        +pow((fV*sin(fTheta)*sin(fPhi)*e[2]),2);
 293    dy2=pow((sin(fPhi)*sin(fTheta)*e[0]),2)+pow((fV*cos(fTheta)*sin(fPhi)*e[1]),2)
 294        +pow((fV*sin(fTheta)*cos(fPhi)*e[2]),2);
 295    dz2=pow((cos(fTheta)*e[0]),2)+pow((fV*sin(fTheta)*e[1]),2);
 296    fDx=sqrt(dx2);
 297    fDy=sqrt(dy2);
 298    fDz=sqrt(dz2);
 299    break;
 300
 301   case 3: // Cylindrical coordinates
 302    rho=fV*sin(fTheta);
 303    dx2=pow((cos(fPhi)*e[0]),2)+pow((rho*sin(fPhi)*e[1]),2);
 304    dy2=pow((sin(fPhi)*e[0]),2)+pow((rho*cos(fPhi)*e[1]),2);
 305    fDx=sqrt(dx2);
 306    fDy=sqrt(dy2);
 307    fDz=fabs(e[2]);
 308    break;
 309
 310   default: // Unsupported reference frame
 311    cout << "*Ali3Vector::SetErrors* Unsupported frame : " << f.Data() << endl
 312         << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
 313    fDx=0;
 314    fDy=0;
 315    fDz=0;
 316    break;
 317  }
 318 }
 319 ///////////////////////////////////////////////////////////////////////////
 320 void Ali3Vector::GetErrors(Double_t* e,TString f) const
 321 {
 322 // Provide errors according to reference frame f
 323  Int_t frame=0;
 324  if (f == "car") frame=1;
 325  if (f == "sph") frame=2;
 326  if (f == "cyl") frame=3;
 327
 328  Double_t pi=acos(-1.);
 329
 330  Double_t dr2,dtheta2,dphi2,rho,drho2;
 331  Double_t v[3];
 332  Double_t rxy2; // Shorthand for (x*x+y*y)
 333
 334  switch (frame)
 335  {
 336   case 1: // Cartesian coordinates
 337    e[0]=fDx;
 338    e[1]=fDy;
 339    e[2]=fDz;
 340    break;
 341
 342   case 2: // Spherical coordinates
 343    GetVector(v,"car");
 344    rxy2=pow(v[0],2)+pow(v[1],2);
 345    if (sqrt(rxy2)<(fV*1e-10)) rxy2=0;
 346    if (fV)
 347    {
 348     dr2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)+pow((v[2]*fDz),2))/(fV*fV);
 349    }
 350    else
 351    {
 352     dr2=0;
 353    }
 354    if (fV)
 355    {
 356     dtheta2=rxy2*pow(fDz,2)/pow(fV,4);
 357     if (v[2] && rxy2)
 358     {
 359      dtheta2+=rxy2*pow(v[2],2)*(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)) /
 360               pow(((pow(v[2],2)*rxy2)+pow(rxy2,2)),2);
 361     }
 362    }
 363    else
 364    {
 365     dtheta2=0;
 366    }
 367    if (rxy2)
 368    {
 369     dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rxy2,2));
 370    }
 371    else
 372    {
 373     dphi2=0;
 374    }
 375    e[0]=sqrt(dr2);
 376    e[1]=sqrt(dtheta2);
 377    if (e[1]>pi) e[1]=pi;
 378    e[2]=sqrt(dphi2);
 379    if (e[2]>(2.*pi)) e[2]=2.*pi;
 380    break;
 381
 382   case 3: // Cylindrical coordinates
 383    GetVector(v,"car");
 384    rho=fabs(fV*sin(fTheta));
 385    if (rho<(fV*1e-10)) rho=0;
 386    if (rho)
 387    {
 388     drho2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2))/(rho*rho);
 389    }
 390    else
 391    {
 392     drho2=0;
 393    }
 394    if (rho)
 395    {
 396     dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rho,4));
 397    }
 398    else
 399    {
 400     dphi2=0;
 401    }
 402    e[0]=sqrt(drho2);
 403    e[1]=sqrt(dphi2);
 404    if (e[1]>(2.*pi)) e[1]=2.*pi;
 405    e[2]=fDz;
 406    break;
 407
 408   default: // Unsupported reference frame
 409    cout << "*Ali3Vector::GetErrors* Unsupported frame : " << f.Data() << endl
 410         << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
 411    for (Int_t i=0; i<3; i++)
 412    {
 413     e[i]=0;
 414    }
 415    break;
 416  }
 417 }
 418 ///////////////////////////////////////////////////////////////////////////
 419 void Ali3Vector::SetErrors(Float_t* e,TString f)
 420 {
 421 // Store errors according to reference frame f
 422 // The error on scalar results is reset to 0
 423  Double_t vec[3];
 424  for (Int_t i=0; i<3; i++)
 425  {
 426   vec[i]=e[i];
 427  }
 428  SetErrors(vec,f);
 429 }
 430 ///////////////////////////////////////////////////////////////////////////
 431 void Ali3Vector::GetErrors(Float_t* e,TString f) const
 432 {
 433 // Provide errors according to reference frame f
 434  Double_t vec[3];
 435  GetErrors(vec,f);
 436  for (Int_t i=0; i<3; i++)
 437  {
 438   e[i]=vec[i];
 439  }
 440 }
 441 ///////////////////////////////////////////////////////////////////////////
 442 void Ali3Vector::Data(TString f) const
 443 {
 444 // Print vector components according to reference frame f
 445  if (f=="car" || f=="sph" || f=="cyl")
 446  {
 447   Double_t vec[3],err[3];
 448   GetVector(vec,f);
 449   GetErrors(err,f);
 450   cout << " Vector in " << f.Data() << " coordinates : "
 451        << vec[0] << " " << vec[1] << " " << vec[2] << endl;
 452   cout << "   Err. in " << f.Data() << " coordinates : "
 453        << err[0] << " " << err[1] << " " << err[2] << endl;
 454  }
 455  else
 456  {
 457   cout << " *Ali3Vector::Data* Unsupported frame : " << f.Data() << endl
 458        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl;
 459  }
 460 }
 461 ///////////////////////////////////////////////////////////////////////////
 462 Double_t Ali3Vector::GetNorm()
 463 {
 464 // Provide the norm of the current vector
 465 // The error on the scalar result (norm) is updated accordingly
 466  Double_t e[3];
 467  GetErrors(e,"sph");
 468  fDresult=e[0];
 469  return fV;
 470 }
 471 ///////////////////////////////////////////////////////////////////////////
 472 Double_t Ali3Vector::GetPseudoRapidity()
 473 {
 474 // Provide the pseudo-rapidity w.r.t. the z-axis.
 475 // In other words : eta=-log(tan(theta/2))
 476 // The error on the scalar result (pseudo-rap.) is updated accordingly
 477  Double_t pi=acos(-1.);
 478  Double_t v[3];
 479  GetVector(v,"sph");
 480  Double_t thetahalf=v[1]/2.;
 481  Double_t arg=0;
 482  if (v[1]<pi) arg=tan(thetahalf);
 483  Double_t eta=9999;
 484  if (arg>0) eta=-log(arg);
 485  Double_t e[3];
 486  GetErrors(e,"sph");
 487  Double_t prod=cos(thetahalf)*sin(thetahalf);
 488  fDresult=0;
 489  if (prod) fDresult=fabs(e[1]/2.*prod);
 490  return eta;
 491 }
 492 ///////////////////////////////////////////////////////////////////////////
 493 Double_t Ali3Vector::Dot(Ali3Vector& q)
 494 {
 495 // Provide the dot product of the current vector with vector q
 496 // The error on the scalar result (dotproduct) is updated accordingly
 497
 498  Double_t dotpro=0;
 499
 500  if ((this) == &q) // Check for special case v.Dot(v)
 501  {
 502   Double_t norm=GetNorm();
 503   Double_t dnorm=GetResultError();
 504   dotpro=pow(norm,2);
 505   fDresult=2.*norm*dnorm;
 506  }
 507  else
 508  {
 509   Double_t a[3],b[3];
 510   Double_t ea[3],eb[3];
 511   Double_t d2=0;
 512
 513   GetVector(a,"car");
 514   GetErrors(ea,"car");
 515   q.GetVector(b,"car");
 516   q.GetErrors(eb,"car");
 517   for (Int_t i=0; i<3; i++)
 518   {
 519    dotpro+=a[i]*b[i];
 520    d2+=pow(b[i]*ea[i],2)+pow(a[i]*eb[i],2);
 521   }
 522   fDresult=sqrt(d2);
 523  }
 524
 525  return dotpro;
 526 }
 527 ///////////////////////////////////////////////////////////////////////////
 528 Double_t Ali3Vector::GetResultError() const
 529 {
 530 // Provide the error on the result of an operation yielding a scalar
 531 // E.g. GetNorm() or Dot()
 532  return fDresult;
 533 }
 534 ///////////////////////////////////////////////////////////////////////////
 535 Ali3Vector Ali3Vector::Cross(Ali3Vector& q) const
 536 {
 537 // Provide the cross product of the current vector with vector q
 538 // Error propagation is performed automatically
 539  Double_t a[3],b[3],c[3];
 540  Double_t ea[3],eb[3],ec[3],d2;
 541
 542  GetVector(a,"car");
 543  GetErrors(ea,"car");
 544  q.GetVector(b,"car");
 545  q.GetErrors(eb,"car");
 546
 547  c[0]=a[1]*b[2]-a[2]*b[1];
 548  c[1]=a[2]*b[0]-a[0]*b[2];
 549  c[2]=a[0]*b[1]-a[1]*b[0];
 550
 551  d2=pow(b[2]*ea[1],2)+pow(a[1]*eb[2],2)
 552    +pow(b[1]*ea[2],2)+pow(a[2]*eb[1],2);
 553  ec[0]=sqrt(d2);
 554
 555  d2=pow(b[0]*ea[2],2)+pow(a[2]*eb[0],2)
 556    +pow(b[2]*ea[0],2)+pow(a[0]*eb[2],2);
 557  ec[1]=sqrt(d2);
 558
 559  d2=pow(b[1]*ea[0],2)+pow(a[0]*eb[1],2)
 560    +pow(b[0]*ea[1],2)+pow(a[1]*eb[0],2);
 561  ec[2]=sqrt(d2);
 562
 563  Ali3Vector v;
 564  v.SetVector(c,"car");
 565  v.SetErrors(ec,"car");
 566
 567  return v;
 568 }
 569 ///////////////////////////////////////////////////////////////////////////
 570 Ali3Vector Ali3Vector::operator+(Ali3Vector& q) const
 571 {
 572 // Add vector q to the current vector
 573 // Error propagation is performed automatically
 574  Double_t a[3],b[3],ea[3],eb[3];
 575
 576  GetVector(a,"car");
 577  GetErrors(ea,"car");
 578  q.GetVector(b,"car");
 579  q.GetErrors(eb,"car");
 580
 581  for (Int_t i=0; i<3; i++)
 582  {
 583   a[i]+=b[i];
 584   ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
 585  }
 586
 587  Ali3Vector v;
 588  v.SetVector(a,"car");
 589  v.SetErrors(ea,"car");
 590
 591  return v;
 592 }
 593 ///////////////////////////////////////////////////////////////////////////
 594 Ali3Vector Ali3Vector::operator-(Ali3Vector& q) const
 595 {
 596 // Subtract vector q from the current vector
 597 // Error propagation is performed automatically
 598  Double_t a[3],b[3],ea[3],eb[3];
 599
 600  GetVector(a,"car");
 601  GetErrors(ea,"car");
 602  q.GetVector(b,"car");
 603  q.GetErrors(eb,"car");
 604
 605  for (Int_t i=0; i<3; i++)
 606  {
 607   a[i]-=b[i];
 608   ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
 609  }
 610
 611  Ali3Vector v;
 612  v.SetVector(a,"car");
 613  v.SetErrors(ea,"car");
 614
 615  return v;
 616 }
 617 ///////////////////////////////////////////////////////////////////////////
 618 Ali3Vector Ali3Vector::operator*(Double_t s) const
 619 {
 620 // Multiply the current vector with a scalar s.
 621 // Error propagation is performed automatically.
 622  Double_t a[3],ea[3];
 623
 624  GetVector(a,"car");
 625  GetErrors(ea,"car");
 626
 627  for (Int_t i=0; i<3; i++)
 628  {
 629   a[i]*=s;
 630   ea[i]*=s;
 631  }
 632
 633  Ali3Vector v;
 634  v.SetVector(a,"car");
 635  v.SetErrors(ea,"car");
 636
 637  return v;
 638 }
 639 ///////////////////////////////////////////////////////////////////////////
 640 Ali3Vector Ali3Vector::operator/(Double_t s) const
 641 {
 642 // Divide the current vector by a scalar s
 643 // Error propagation is performed automatically
 644
 645  if (fabs(s)<1.e-20) // Protect against division by 0
 646  {
 647   cout << " *Ali3Vector::/* Division by 0 detected. No action taken." << endl;
 648   return *this;
 649  }
 650  else
 651  {
 652   Double_t a[3],ea[3];
 653
 654   GetVector(a,"car");
 655   GetErrors(ea,"car");
 656
 657   for (Int_t i=0; i<3; i++)
 658   {
 659    a[i]/=s;
 660    ea[i]/=s;
 661   }
 662
 663   Ali3Vector v;
 664   v.SetVector(a,"car");
 665   v.SetErrors(ea,"car");
 666
 667   return v;
 668  }
 669 }
 670 ///////////////////////////////////////////////////////////////////////////
 671 Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
 672 {
 673 // Add vector q to the current vector
 674 // Error propagation is performed automatically
 675  Double_t a[3],b[3],ea[3],eb[3];
 676
 677  GetVector(a,"car");
 678  GetErrors(ea,"car");
 679  q.GetVector(b,"car");
 680  q.GetErrors(eb,"car");
 681
 682  for (Int_t i=0; i<3; i++)
 683  {
 684   a[i]+=b[i];
 685   ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
 686  }
 687
 688  SetVector(a,"car");
 689  SetErrors(ea,"car");
 690
 691  return *this;
 692 }
 693 ///////////////////////////////////////////////////////////////////////////
 694 Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
 695 {
 696 // Subtract vector q from the current vector
 697 // Error propagation is performed automatically
 698  Double_t a[3],b[3],ea[3],eb[3];
 699
 700  GetVector(a,"car");
 701  GetErrors(ea,"car");
 702  q.GetVector(b,"car");
 703  q.GetErrors(eb,"car");
 704
 705  for (Int_t i=0; i<3; i++)
 706  {
 707   a[i]-=b[i];
 708   ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
 709  }
 710
 711  SetVector(a,"car");
 712  SetErrors(ea,"car");
 713
 714  return *this;
 715 }
 716 ///////////////////////////////////////////////////////////////////////////
 717 Ali3Vector& Ali3Vector::operator*=(Double_t s)
 718 {
 719 // Multiply the current vector with a scalar s
 720 // Error propagation is performed automatically
 721  Double_t a[3],ea[3];
 722
 723  GetVector(a,"car");
 724  GetErrors(ea,"car");
 725
 726  for (Int_t i=0; i<3; i++)
 727  {
 728   a[i]*=s;
 729   ea[i]*=s;
 730  }
 731
 732  SetVector(a,"car");
 733  SetErrors(ea,"car");
 734
 735  return *this;
 736 }
 737 ///////////////////////////////////////////////////////////////////////////
 738 Ali3Vector& Ali3Vector::operator/=(Double_t s)
 739 {
 740 // Divide the current vector by a scalar s
 741 // Error propagation is performed automatically
 742
 743  if (fabs(s)<1.e-20) // Protect against division by 0
 744  {
 745   cout << " *Ali3Vector::/=* Division by 0 detected. No action taken." << endl;
 746   return *this;
 747  }
 748  else
 749  {
 750   Double_t a[3],ea[3];
 751
 752   GetVector(a,"car");
 753   GetErrors(ea,"car");
 754
 755   for (Int_t i=0; i<3; i++)
 756   {
 757    a[i]/=s;
 758    ea[i]/=s;
 759   }
 760
 761   SetVector(a,"car");
 762   SetErrors(ea,"car");
 763
 764   return *this;
 765  }
 766 }
 767 ///////////////////////////////////////////////////////////////////////////
 768 Ali3Vector Ali3Vector::GetVecTrans() const
 769 {
 770 // Provide the transverse vector w.r.t. z-axis.
 771 // Error propagation is performed automatically
 772  Double_t pi=acos(-1.);
 773  Double_t a[3],ea[3];
 774
 775  GetVector(a,"sph");
 776  GetErrors(ea,"sph");
 777
 778  Double_t vt,dvt2;
 779  vt=a[0]*sin(a[1]);
 780  dvt2=pow((sin(a[1])*ea[0]),2)+pow((a[0]*cos(a[1])*ea[1]),2);
 781
 782  a[0]=fabs(vt);
 783  a[1]=pi/2.;
 784
 785  ea[0]=sqrt(dvt2);
 786  ea[1]=0;
 787
 788  Ali3Vector v;
 789  v.SetVector(a,"sph");
 790  v.SetErrors(ea,"sph");
 791
 792  return v;
 793 }
 794 ///////////////////////////////////////////////////////////////////////////
 795 Ali3Vector Ali3Vector::GetVecLong() const
 796 {
 797 // Provide the longitudinal vector w.r.t. z-axis.
 798 // Error propagation is performed automatically
 799  Double_t pi=acos(-1.);
 800  Double_t a[3],ea[3];
 801
 802  GetVector(a,"sph");
 803  GetErrors(ea,"sph");
 804
 805  Double_t vl,dvl2;
 806  vl=a[0]*cos(a[1]);
 807  dvl2=pow((cos(a[1])*ea[0]),2)+pow((a[0]*sin(a[1])*ea[1]),2);
 808
 809  a[0]=fabs(vl);
 810  a[1]=0;
 811  if (vl<0) a[1]=pi;
 812  a[2]=0;
 813
 814  ea[0]=sqrt(dvl2);
 815  ea[1]=0;
 816  ea[2]=0;
 817
 818  Ali3Vector v;
 819  v.SetVector(a,"sph");
 820  v.SetErrors(ea,"sph");
 821
 822  return v;
 823 }
 824 ///////////////////////////////////////////////////////////////////////////
 825 Ali3Vector Ali3Vector::GetPrimed(TRotMatrix* m) const
 826 {
 827 // Provide vector components (and errors) in a rotated frame.
 828 // The orientation of the rotated frame is described by the TRotMatrix
 829 // input argument.
 830  Ali3Vector v=*this;
 831  if (!m) return v;
 832
 833  Double_t* mat=m->GetMatrix();
 834
 835  Double_t a[3],aprim[3];
 836
 837  GetVector(a,"car");
 838  aprim[0]=a[0]*mat[0]+a[1]*mat[1]+a[2]*mat[2];
 839  aprim[1]=a[0]*mat[3]+a[1]*mat[4]+a[2]*mat[5];
 840  aprim[2]=a[0]*mat[6]+a[1]*mat[7]+a[2]*mat[8];
 841  v.SetVector(aprim,"car");
 842
 843  GetErrors(a,"car");
 844  aprim[0]=sqrt(pow(a[0]*mat[0],2)+pow(a[1]*mat[1],2)+pow(a[2]*mat[2],2));
 845  aprim[1]=sqrt(pow(a[0]*mat[3],2)+pow(a[1]*mat[4],2)+pow(a[2]*mat[5],2));
 846  aprim[2]=sqrt(pow(a[0]*mat[6],2)+pow(a[1]*mat[7],2)+pow(a[2]*mat[8],2));
 847  v.SetErrors(aprim,"car");
 848
 849  return v;
 850 }
 851 ///////////////////////////////////////////////////////////////////////////
 852 Ali3Vector Ali3Vector::GetUnprimed(TRotMatrix* m) const
 853 {
 854 // Provide original vector components (and errors) from the rotated ones.
 855 // The orientation of the rotated frame is described by the TRotMatrix
 856 // input argument.
 857 // So, this is the inverse of the GetPrimed() memberfunction.
 858 // This memberfunction makes use of the fact that the inverse of a certain
 859 // TRotMatrix is given by its transposed matrix.
 860  Ali3Vector v=*this;
 861  if (!m) return v;
 862
 863  Double_t* mat=m->GetMatrix();
 864
 865  Double_t a[3],aprim[3];
 866
 867  GetVector(aprim,"car");
 868  a[0]=aprim[0]*mat[0]+aprim[1]*mat[3]+aprim[2]*mat[6];
 869  a[1]=aprim[0]*mat[1]+aprim[1]*mat[4]+aprim[2]*mat[7];
 870  a[2]=aprim[0]*mat[2]+aprim[1]*mat[5]+aprim[2]*mat[8];
 871  v.SetVector(a,"car");
 872
 873  GetErrors(aprim,"car");
 874  a[0]=sqrt(pow(aprim[0]*mat[0],2)+pow(aprim[1]*mat[3],2)+pow(aprim[2]*mat[6],2));
 875  a[1]=sqrt(pow(aprim[0]*mat[1],2)+pow(aprim[1]*mat[4],2)+pow(aprim[2]*mat[7],2));
 876  a[2]=sqrt(pow(aprim[0]*mat[2],2)+pow(aprim[1]*mat[5],2)+pow(aprim[2]*mat[8],2));
 877  v.SetErrors(a,"car");
 878
 879  return v;
 880 }
 881 ///////////////////////////////////////////////////////////////////////////
 882 Double_t Ali3Vector::GetX(Int_t i,TString f)
 883 {
 884 // Provide i-th vector component according to reference frame f.
 885 // The vector components are addressed via the generic x1,x2,x3 notation.
 886 // So, i=1 denotes the first vector component.
 887 // The error on the selected component can be obtained via the
 888 // usual GetResultError() facility.
 889
 890  fDresult=0;
 891
 892  if (i<1 || i>3) return 0;
 893
 894  Double_t vec[3];
 895  Double_t err[3];
 896  GetVector(vec,f);
 897  GetErrors(err,f);
 898
 899  fDresult=err[i-1];
 900
 901  return vec[i-1];
 902 }
 903 ///////////////////////////////////////////////////////////////////////////