1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////
20 // Handling of 3-vectors in various reference frames.
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.
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 :
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)
37 // All angles are in radians.
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");
51 // a.GetVector(vec,"sph");
52 // a.GetErrors(vec,"sph");
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");
60 // Float_t dotpro=a.Dot(b);
61 // Float_t doterror=a.GetResultError();
63 // Ali3Vector c=a.Cross(b);
65 // c.GetVector(vec,"cyl");
66 // c.GetErrors(err,"cyl");
68 // Float_t norm=c.GetNorm();
69 // Float_t normerror=c.GetResultError();
75 //--- Author: Nick van Eijndhoven 30-mar-1999 UU-SAP Utrecht
76 //- Modified: NvE $Date$ UU-SAP Utrecht
77 ///////////////////////////////////////////////////////////////////////////
79 #include "Ali3Vector.h"
80 #include "Riostream.h"
82 ClassImp(Ali3Vector) // Class implementation to enable ROOT I/O
84 Ali3Vector::Ali3Vector()
86 // Creation of an Ali3Vector object and initialisation of parameters
87 // All attributes initialised to 0
90 ///////////////////////////////////////////////////////////////////////////
91 Ali3Vector::~Ali3Vector()
93 // Destructor to delete dynamically allocated memory
95 ///////////////////////////////////////////////////////////////////////////
96 Ali3Vector::Ali3Vector(const Ali3Vector& v)
107 ///////////////////////////////////////////////////////////////////////////
108 void Ali3Vector::SetZero()
110 // (Re)set all attributes to zero.
119 ///////////////////////////////////////////////////////////////////////////
120 void Ali3Vector::SetVector(Double_t* v,TString f)
122 // Store vector according to reference frame f
123 // All errors will be reset to 0
129 Double_t pi=acos(-1.);
132 if (f == "car") frame=1;
133 if (f == "sph") frame=2;
134 if (f == "cyl") frame=3;
136 Double_t x,y,z,rho,phi;
140 case 1: // Cartesian coordinates
144 fV=sqrt(x*x+y*y+z*z);
146 if (fV && fabs(z/fV)<=1.)
154 if (fTheta<0.) fTheta+=2.*pi;
156 if (x || y) fPhi=atan2(y,x);
157 if (fPhi<0.) fPhi+=2.*pi;
160 case 2: // Spherical coordinates
166 case 3: // Cylindrical coordinates
170 fV=sqrt(rho*rho+z*z);
172 if (fPhi<0.) fPhi+=2.*pi;
174 if (fV && fabs(z/fV)<=1.)
182 if (fTheta<0.) fTheta+=2.*pi;
185 default: // Unsupported reference frame
186 cout << "*Ali3Vector::SetVector* Unsupported frame : " << f.Data() << endl
187 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
194 ///////////////////////////////////////////////////////////////////////////
195 void Ali3Vector::GetVector(Double_t* v,TString f)
197 // Provide vector according to reference frame f
199 if (f == "car") frame=1;
200 if (f == "sph") frame=2;
201 if (f == "cyl") frame=3;
205 case 1: // Cartesian coordinates
206 v[0]=fV*sin(fTheta)*cos(fPhi);
207 v[1]=fV*sin(fTheta)*sin(fPhi);
211 case 2: // Spherical coordinates
217 case 3: // Cylindrical coordinates
223 default: // Unsupported reference frame
224 cout << "*Ali3Vector::GetVector* Unsupported frame : " << f.Data() << endl
225 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
226 for (Int_t i=0; i<3; i++)
233 ///////////////////////////////////////////////////////////////////////////
234 void Ali3Vector::SetVector(Float_t* v,TString f)
236 // Store vector according to reference frame f
237 // All errors will be reset to 0
239 for (Int_t i=0; i<3; i++)
245 ///////////////////////////////////////////////////////////////////////////
246 void Ali3Vector::GetVector(Float_t* v,TString f)
248 // Provide vector according to reference frame f
251 for (Int_t i=0; i<3; i++)
256 ///////////////////////////////////////////////////////////////////////////
257 void Ali3Vector::SetErrors(Double_t* e,TString f)
259 // Store errors according to reference frame f
260 // The error on scalar results is reset to 0
264 if (f == "car") frame=1;
265 if (f == "sph") frame=2;
266 if (f == "cyl") frame=3;
268 Double_t dx2,dy2,dz2,rho;
272 case 1: // Cartesian coordinates
278 case 2: // Spherical coordinates
279 dx2=pow((cos(fPhi)*sin(fTheta)*e[0]),2)+pow((fV*cos(fTheta)*cos(fPhi)*e[1]),2)
280 +pow((fV*sin(fTheta)*sin(fPhi)*e[2]),2);
281 dy2=pow((sin(fPhi)*sin(fTheta)*e[0]),2)+pow((fV*cos(fTheta)*sin(fPhi)*e[1]),2)
282 +pow((fV*sin(fTheta)*cos(fPhi)*e[2]),2);
283 dz2=pow((cos(fTheta)*e[0]),2)+pow((fV*sin(fTheta)*e[1]),2);
289 case 3: // Cylindrical coordinates
291 dx2=pow((cos(fPhi)*e[0]),2)+pow((rho*sin(fPhi)*e[1]),2);
292 dy2=pow((sin(fPhi)*e[0]),2)+pow((rho*cos(fPhi)*e[1]),2);
298 default: // Unsupported reference frame
299 cout << "*Ali3Vector::SetErrors* Unsupported frame : " << f.Data() << endl
300 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
307 ///////////////////////////////////////////////////////////////////////////
308 void Ali3Vector::GetErrors(Double_t* e,TString f)
310 // Provide errors according to reference frame f
312 if (f == "car") frame=1;
313 if (f == "sph") frame=2;
314 if (f == "cyl") frame=3;
316 Double_t pi=acos(-1.);
318 Double_t dr2,dtheta2,dphi2,rho,drho2;
320 Double_t rxy2; // Shorthand for (x*x+y*y)
324 case 1: // Cartesian coordinates
330 case 2: // Spherical coordinates
332 rxy2=pow(v[0],2)+pow(v[1],2);
333 if (sqrt(rxy2)<(fV*1e-10)) rxy2=0;
336 dr2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)+pow((v[2]*fDz),2))/(fV*fV);
344 dtheta2=rxy2*pow(fDz,2)/pow(fV,4);
347 dtheta2+=rxy2*pow(v[2],2)*(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)) /
348 pow(((pow(v[2],2)*rxy2)+pow(rxy2,2)),2);
357 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rxy2,2));
365 if (e[1]>pi) e[1]=pi;
367 if (e[2]>(2.*pi)) e[2]=2.*pi;
370 case 3: // Cylindrical coordinates
372 rho=fabs(fV*sin(fTheta));
373 if (rho<(fV*1e-10)) rho=0;
376 drho2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2))/(rho*rho);
384 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rho,4));
392 if (e[1]>(2.*pi)) e[1]=2.*pi;
396 default: // Unsupported reference frame
397 cout << "*Ali3Vector::GetErrors* Unsupported frame : " << f.Data() << endl
398 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
399 for (Int_t i=0; i<3; i++)
406 ///////////////////////////////////////////////////////////////////////////
407 void Ali3Vector::SetErrors(Float_t* e,TString f)
409 // Store errors according to reference frame f
410 // The error on scalar results is reset to 0
412 for (Int_t i=0; i<3; i++)
418 ///////////////////////////////////////////////////////////////////////////
419 void Ali3Vector::GetErrors(Float_t* e,TString f)
421 // Provide errors according to reference frame f
424 for (Int_t i=0; i<3; i++)
429 ///////////////////////////////////////////////////////////////////////////
430 void Ali3Vector::Data(TString f)
432 // Print vector components according to reference frame f
433 if (f=="car" || f=="sph" || f=="cyl")
435 Double_t vec[3],err[3];
438 cout << " Vector in " << f.Data() << " coordinates : "
439 << vec[0] << " " << vec[1] << " " << vec[2] << endl;
440 cout << " Err. in " << f.Data() << " coordinates : "
441 << err[0] << " " << err[1] << " " << err[2] << endl;
445 cout << " *Ali3Vector::Data* Unsupported frame : " << f.Data() << endl
446 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
449 ///////////////////////////////////////////////////////////////////////////
450 Double_t Ali3Vector::GetNorm()
452 // Provide the norm of the current vector
453 // The error on the scalar result (norm) is updated accordingly
459 ///////////////////////////////////////////////////////////////////////////
460 Double_t Ali3Vector::GetPseudoRapidity()
462 // Provide the pseudo-rapidity w.r.t. the z-axis.
463 // In other words : eta=-log(tan(theta/2))
464 // The error on the scalar result (pseudo-rap.) is updated accordingly
465 Double_t pi=acos(-1.);
468 Double_t thetahalf=v[1]/2.;
470 if (v[1]<pi) arg=tan(thetahalf);
472 if (arg>0) eta=-log(arg);
475 Double_t prod=cos(thetahalf)*sin(thetahalf);
477 if (prod) fDresult=fabs(e[1]/2.*prod);
480 ///////////////////////////////////////////////////////////////////////////
481 Double_t Ali3Vector::Dot(Ali3Vector& q)
483 // Provide the dot product of the current vector with vector q
484 // The error on the scalar result (dotproduct) is updated accordingly
488 if ((this) == &q) // Check for special case v.Dot(v)
490 Double_t norm=GetNorm();
491 Double_t dnorm=GetResultError();
493 fDresult=2.*norm*dnorm;
498 Double_t ea[3],eb[3];
503 q.GetVector(b,"car");
504 q.GetErrors(eb,"car");
505 for (Int_t i=0; i<3; i++)
508 d2+=pow(b[i]*ea[i],2)+pow(a[i]*eb[i],2);
515 ///////////////////////////////////////////////////////////////////////////
516 Double_t Ali3Vector::GetResultError()
518 // Provide the error on the result of an operation yielding a scalar
519 // E.g. GetNorm() or Dot()
522 ///////////////////////////////////////////////////////////////////////////
523 Ali3Vector Ali3Vector::Cross(Ali3Vector& q)
525 // Provide the cross product of the current vector with vector q
526 // Error propagation is performed automatically
527 Double_t a[3],b[3],c[3];
528 Double_t ea[3],eb[3],ec[3],d2;
532 q.GetVector(b,"car");
533 q.GetErrors(eb,"car");
535 c[0]=a[1]*b[2]-a[2]*b[1];
536 c[1]=a[2]*b[0]-a[0]*b[2];
537 c[2]=a[0]*b[1]-a[1]*b[0];
539 d2=pow(b[2]*ea[1],2)+pow(a[1]*eb[2],2)
540 +pow(b[1]*ea[2],2)+pow(a[2]*eb[1],2);
543 d2=pow(b[0]*ea[2],2)+pow(a[2]*eb[0],2)
544 +pow(b[2]*ea[0],2)+pow(a[0]*eb[2],2);
547 d2=pow(b[1]*ea[0],2)+pow(a[0]*eb[1],2)
548 +pow(b[0]*ea[1],2)+pow(a[1]*eb[0],2);
552 v.SetVector(c,"car");
553 v.SetErrors(ec,"car");
557 ///////////////////////////////////////////////////////////////////////////
558 Ali3Vector Ali3Vector::operator+(Ali3Vector& q)
560 // Add vector q to the current vector
561 // Error propagation is performed automatically
562 Double_t a[3],b[3],ea[3],eb[3];
566 q.GetVector(b,"car");
567 q.GetErrors(eb,"car");
569 for (Int_t i=0; i<3; i++)
572 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
576 v.SetVector(a,"car");
577 v.SetErrors(ea,"car");
581 ///////////////////////////////////////////////////////////////////////////
582 Ali3Vector Ali3Vector::operator-(Ali3Vector& q)
584 // Subtract vector q from the current vector
585 // Error propagation is performed automatically
586 Double_t a[3],b[3],ea[3],eb[3];
590 q.GetVector(b,"car");
591 q.GetErrors(eb,"car");
593 for (Int_t i=0; i<3; i++)
596 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
600 v.SetVector(a,"car");
601 v.SetErrors(ea,"car");
605 ///////////////////////////////////////////////////////////////////////////
606 Ali3Vector Ali3Vector::operator*(Double_t s)
608 // Multiply the current vector with a scalar s.
609 // Error propagation is performed automatically.
615 for (Int_t i=0; i<3; i++)
622 v.SetVector(a,"car");
623 v.SetErrors(ea,"car");
627 ///////////////////////////////////////////////////////////////////////////
628 Ali3Vector Ali3Vector::operator/(Double_t s)
630 // Divide the current vector by a scalar s
631 // Error propagation is performed automatically
633 if (fabs(s)<1.e-20) // Protect against division by 0
635 cout << " *Ali3Vector::/* Division by 0 detected. No action taken." << endl;
645 for (Int_t i=0; i<3; i++)
652 v.SetVector(a,"car");
653 v.SetErrors(ea,"car");
658 ///////////////////////////////////////////////////////////////////////////
659 Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
661 // Add vector q to the current vector
662 // Error propagation is performed automatically
663 Double_t a[3],b[3],ea[3],eb[3];
667 q.GetVector(b,"car");
668 q.GetErrors(eb,"car");
670 for (Int_t i=0; i<3; i++)
673 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
681 ///////////////////////////////////////////////////////////////////////////
682 Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
684 // Subtract vector q from the current vector
685 // Error propagation is performed automatically
686 Double_t a[3],b[3],ea[3],eb[3];
690 q.GetVector(b,"car");
691 q.GetErrors(eb,"car");
693 for (Int_t i=0; i<3; i++)
696 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
704 ///////////////////////////////////////////////////////////////////////////
705 Ali3Vector& Ali3Vector::operator*=(Double_t s)
707 // Multiply the current vector with a scalar s
708 // Error propagation is performed automatically
714 for (Int_t i=0; i<3; i++)
725 ///////////////////////////////////////////////////////////////////////////
726 Ali3Vector& Ali3Vector::operator/=(Double_t s)
728 // Divide the current vector by a scalar s
729 // Error propagation is performed automatically
731 if (fabs(s)<1.e-20) // Protect against division by 0
733 cout << " *Ali3Vector::/=* Division by 0 detected. No action taken." << endl;
743 for (Int_t i=0; i<3; i++)
755 ///////////////////////////////////////////////////////////////////////////
756 Ali3Vector Ali3Vector::GetVecTrans()
758 // Provide the transverse vector w.r.t. z-axis.
759 // Error propagation is performed automatically
760 Double_t pi=acos(-1.);
768 dvt2=pow((sin(a[1])*ea[0]),2)+pow((a[0]*cos(a[1])*ea[1]),2);
777 v.SetVector(a,"sph");
778 v.SetErrors(ea,"sph");
782 ///////////////////////////////////////////////////////////////////////////
783 Ali3Vector Ali3Vector::GetVecLong()
785 // Provide the longitudinal vector w.r.t. z-axis.
786 // Error propagation is performed automatically
787 Double_t pi=acos(-1.);
795 dvl2=pow((cos(a[1])*ea[0]),2)+pow((a[0]*sin(a[1])*ea[1]),2);
807 v.SetVector(a,"sph");
808 v.SetErrors(ea,"sph");
812 ///////////////////////////////////////////////////////////////////////////