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