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"
81 ClassImp(Ali3Vector) // Class implementation to enable ROOT I/O
83 Ali3Vector::Ali3Vector()
85 // Creation of an Ali3Vector object and initialisation of parameters
86 // All attributes initialised to 0
95 ///////////////////////////////////////////////////////////////////////////
96 Ali3Vector::~Ali3Vector()
98 // Destructor to delete dynamically allocated memory
100 ///////////////////////////////////////////////////////////////////////////
101 void Ali3Vector::SetVector(Double_t* v,TString f)
103 // Store vector according to reference frame f
104 // All errors will be reset to 0
110 Double_t pi=acos(-1.);
113 if (f == "car") frame=1;
114 if (f == "sph") frame=2;
115 if (f == "cyl") frame=3;
117 Double_t x,y,z,rho,phi;
121 case 1: // Cartesian coordinates
125 fV=sqrt(x*x+y*y+z*z);
127 if (fV && fabs(z/fV)<=1.)
135 if (fTheta<0.) fTheta+=2.*pi;
137 if (x || y) fPhi=atan2(y,x);
138 if (fPhi<0.) fPhi+=2.*pi;
141 case 2: // Spherical coordinates
147 case 3: // Cylindrical coordinates
151 fV=sqrt(rho*rho+z*z);
153 if (fPhi<0.) fPhi+=2.*pi;
155 if (fV && fabs(z/fV)<=1.)
163 if (fTheta<0.) fTheta+=2.*pi;
166 default: // Unsupported reference frame
167 cout << "*Ali3Vector::SetVector* Unsupported frame : " << f << endl
168 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
175 ///////////////////////////////////////////////////////////////////////////
176 void Ali3Vector::GetVector(Double_t* v,TString f)
178 // Provide vector according to reference frame f
180 if (f == "car") frame=1;
181 if (f == "sph") frame=2;
182 if (f == "cyl") frame=3;
186 case 1: // Cartesian coordinates
187 v[0]=fV*sin(fTheta)*cos(fPhi);
188 v[1]=fV*sin(fTheta)*sin(fPhi);
192 case 2: // Spherical coordinates
198 case 3: // Cylindrical coordinates
204 default: // Unsupported reference frame
205 cout << "*Ali3Vector::GetVector* Unsupported frame : " << f << endl
206 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
207 for (Int_t i=0; i<3; i++)
214 ///////////////////////////////////////////////////////////////////////////
215 void Ali3Vector::SetVector(Float_t* v,TString f)
217 // Store vector according to reference frame f
218 // All errors will be reset to 0
220 for (Int_t i=0; i<3; i++)
226 ///////////////////////////////////////////////////////////////////////////
227 void Ali3Vector::GetVector(Float_t* v,TString f)
229 // Provide vector according to reference frame f
232 for (Int_t i=0; i<3; i++)
237 ///////////////////////////////////////////////////////////////////////////
238 void Ali3Vector::SetErrors(Double_t* e,TString f)
240 // Store errors according to reference frame f
241 // The error on scalar results is reset to 0
245 if (f == "car") frame=1;
246 if (f == "sph") frame=2;
247 if (f == "cyl") frame=3;
249 Double_t dx2,dy2,dz2,rho;
253 case 1: // Cartesian coordinates
259 case 2: // Spherical coordinates
260 dx2=pow((cos(fPhi)*sin(fTheta)*e[0]),2)+pow((fV*cos(fTheta)*cos(fPhi)*e[1]),2)
261 +pow((fV*sin(fTheta)*sin(fPhi)*e[2]),2);
262 dy2=pow((sin(fPhi)*sin(fTheta)*e[0]),2)+pow((fV*cos(fTheta)*sin(fPhi)*e[1]),2)
263 +pow((fV*sin(fTheta)*cos(fPhi)*e[2]),2);
264 dz2=pow((cos(fTheta)*e[0]),2)+pow((fV*sin(fTheta)*e[1]),2);
270 case 3: // Cylindrical coordinates
272 dx2=pow((cos(fPhi)*e[0]),2)+pow((rho*sin(fPhi)*e[1]),2);
273 dy2=pow((sin(fPhi)*e[0]),2)+pow((rho*cos(fPhi)*e[1]),2);
279 default: // Unsupported reference frame
280 cout << "*Ali3Vector::SetErrors* Unsupported frame : " << f << endl
281 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
288 ///////////////////////////////////////////////////////////////////////////
289 void Ali3Vector::GetErrors(Double_t* e,TString f)
291 // Provide errors according to reference frame f
293 if (f == "car") frame=1;
294 if (f == "sph") frame=2;
295 if (f == "cyl") frame=3;
297 Double_t pi=acos(-1.);
299 Double_t dr2,dtheta2,dphi2,rho,drho2;
301 Double_t rxy2; // Shorthand for (x*x+y*y)
305 case 1: // Cartesian coordinates
311 case 2: // Spherical coordinates
313 rxy2=pow(v[0],2)+pow(v[1],2);
314 if (sqrt(rxy2)<(fV*1e-10)) rxy2=0;
317 dr2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)+pow((v[2]*fDz),2))/(fV*fV);
325 dtheta2=rxy2*pow(fDz,2)/pow(fV,4);
328 dtheta2+=rxy2*pow(v[2],2)*(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)) /
329 pow(((pow(v[3],2)*rxy2)+pow(rxy2,2)),2);
338 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rxy2,2));
346 if (e[1]>pi) e[1]=pi;
348 if (e[2]>(2.*pi)) e[2]=2.*pi;
351 case 3: // Cylindrical coordinates
353 rho=fabs(fV*sin(fTheta));
354 if (rho<(fV*1e-10)) rho=0;
357 drho2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2))/(rho*rho);
365 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rho,4));
373 if (e[1]>(2.*pi)) e[1]=2.*pi;
377 default: // Unsupported reference frame
378 cout << "*Ali3Vector::GetErrors* Unsupported frame : " << f << endl
379 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
380 for (Int_t i=0; i<3; i++)
387 ///////////////////////////////////////////////////////////////////////////
388 void Ali3Vector::SetErrors(Float_t* e,TString f)
390 // Store errors according to reference frame f
391 // The error on scalar results is reset to 0
393 for (Int_t i=0; i<3; i++)
399 ///////////////////////////////////////////////////////////////////////////
400 void Ali3Vector::GetErrors(Float_t* e,TString f)
402 // Provide errors according to reference frame f
405 for (Int_t i=0; i<3; i++)
410 ///////////////////////////////////////////////////////////////////////////
411 void Ali3Vector::Info(TString f)
413 // Print vector components according to reference frame f
414 if (f=="car" || f=="sph" || f=="cyl")
416 Double_t vec[3],err[3];
419 cout << " Vector in " << f << " coordinates : "
420 << vec[0] << " " << vec[1] << " " << vec[2] << endl;
421 cout << " Err. in " << f << " coordinates : "
422 << err[0] << " " << err[1] << " " << err[2] << endl;
426 cout << " *Ali3Vector::Info* Unsupported frame : " << f << endl
427 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
430 ///////////////////////////////////////////////////////////////////////////
431 Double_t Ali3Vector::GetNorm()
433 // Provide the norm of the current vector
434 // The error on the scalar result (norm) is updated accordingly
440 ///////////////////////////////////////////////////////////////////////////
441 Double_t Ali3Vector::GetPseudoRapidity()
443 // Provide the pseudo-rapidity w.r.t. the z-axis.
444 // In other words : eta=-log(tan(theta/2))
445 // The error on the scalar result (pseudo-rap.) is updated accordingly
446 Double_t pi=acos(-1.);
449 Double_t thetahalf=v[1]/2.;
451 if (v[1]<pi) arg=tan(thetahalf);
453 if (arg>0) eta=-log(arg);
456 Double_t prod=cos(thetahalf)*sin(thetahalf);
458 if (prod) fDresult=fabs(e[1]/2.*prod);
461 ///////////////////////////////////////////////////////////////////////////
462 Double_t Ali3Vector::Dot(Ali3Vector& q)
464 // Provide the dot product of the current vector with vector q
465 // The error on the scalar result (dotproduct) is updated accordingly
469 if ((this) == &q) // Check for special case v.Dot(v)
471 Double_t norm=GetNorm();
472 Double_t dnorm=GetResultError();
474 fDresult=2.*norm*dnorm;
479 Double_t ea[3],eb[3];
484 q.GetVector(b,"car");
485 q.GetErrors(eb,"car");
486 for (Int_t i=0; i<3; i++)
489 d2+=pow(b[i]*ea[i],2)+pow(a[i]*eb[i],2);
496 ///////////////////////////////////////////////////////////////////////////
497 Double_t Ali3Vector::GetResultError()
499 // Provide the error on the result of an operation yielding a scalar
500 // E.g. GetNorm() or Dot()
503 ///////////////////////////////////////////////////////////////////////////
504 Ali3Vector Ali3Vector::Cross(Ali3Vector& q)
506 // Provide the cross product of the current vector with vector q
507 // Error propagation is performed automatically
508 Double_t a[3],b[3],c[3];
509 Double_t ea[3],eb[3],ec[3],d2;
513 q.GetVector(b,"car");
514 q.GetErrors(eb,"car");
516 c[0]=a[1]*b[2]-a[2]*b[1];
517 c[1]=a[2]*b[0]-a[0]*b[2];
518 c[2]=a[0]*b[1]-a[1]*b[0];
520 d2=pow(b[2]*ea[1],2)+pow(a[1]*eb[2],2)
521 +pow(b[1]*ea[2],2)+pow(a[2]*eb[1],2);
524 d2=pow(b[0]*ea[2],2)+pow(a[2]*eb[0],2)
525 +pow(b[2]*ea[0],2)+pow(a[0]*eb[2],2);
528 d2=pow(b[1]*ea[0],2)+pow(a[0]*eb[1],2)
529 +pow(b[0]*ea[1],2)+pow(a[1]*eb[0],2);
533 v.SetVector(c,"car");
534 v.SetErrors(ec,"car");
538 ///////////////////////////////////////////////////////////////////////////
539 Ali3Vector Ali3Vector::operator+(Ali3Vector& q)
541 // Add vector q to the current vector
542 // Error propagation is performed automatically
543 Double_t a[3],b[3],ea[3],eb[3];
547 q.GetVector(b,"car");
548 q.GetErrors(eb,"car");
550 for (Int_t i=0; i<3; i++)
553 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
557 v.SetVector(a,"car");
558 v.SetErrors(ea,"car");
562 ///////////////////////////////////////////////////////////////////////////
563 Ali3Vector Ali3Vector::operator-(Ali3Vector& q)
565 // Subtract vector q from the current vector
566 // Error propagation is performed automatically
567 Double_t a[3],b[3],ea[3],eb[3];
571 q.GetVector(b,"car");
572 q.GetErrors(eb,"car");
574 for (Int_t i=0; i<3; i++)
577 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
581 v.SetVector(a,"car");
582 v.SetErrors(ea,"car");
586 ///////////////////////////////////////////////////////////////////////////
587 Ali3Vector Ali3Vector::operator*(Double_t s)
589 // Multiply the current vector with a scalar s.
590 // Error propagation is performed automatically.
596 for (Int_t i=0; i<3; i++)
603 v.SetVector(a,"car");
604 v.SetErrors(ea,"car");
608 ///////////////////////////////////////////////////////////////////////////
609 Ali3Vector Ali3Vector::operator/(Double_t s)
611 // Divide the current vector by a scalar s
612 // Error propagation is performed automatically
614 if (fabs(s)<1.e-20) // Protect against division by 0
616 cout << " *Ali3Vector::/* Division by 0 detected. No action taken." << endl;
626 for (Int_t i=0; i<3; i++)
633 v.SetVector(a,"car");
634 v.SetErrors(ea,"car");
639 ///////////////////////////////////////////////////////////////////////////
640 Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
642 // Add vector q to the current vector
643 // Error propagation is performed automatically
644 Double_t a[3],b[3],ea[3],eb[3];
648 q.GetVector(b,"car");
649 q.GetErrors(eb,"car");
651 for (Int_t i=0; i<3; i++)
654 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
662 ///////////////////////////////////////////////////////////////////////////
663 Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
665 // Subtract vector q from the current vector
666 // Error propagation is performed automatically
667 Double_t a[3],b[3],ea[3],eb[3];
671 q.GetVector(b,"car");
672 q.GetErrors(eb,"car");
674 for (Int_t i=0; i<3; i++)
677 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
685 ///////////////////////////////////////////////////////////////////////////
686 Ali3Vector& Ali3Vector::operator*=(Double_t s)
688 // Multiply the current vector with a scalar s
689 // Error propagation is performed automatically
695 for (Int_t i=0; i<3; i++)
706 ///////////////////////////////////////////////////////////////////////////
707 Ali3Vector& Ali3Vector::operator/=(Double_t s)
709 // Divide the current vector by a scalar s
710 // Error propagation is performed automatically
712 if (fabs(s)<1.e-20) // Protect against division by 0
714 cout << " *Ali3Vector::/=* Division by 0 detected. No action taken." << endl;
724 for (Int_t i=0; i<3; i++)
736 ///////////////////////////////////////////////////////////////////////////
737 Ali3Vector Ali3Vector::GetVecTrans()
739 // Provide the transverse vector w.r.t. z-axis.
740 // Error propagation is performed automatically
741 Double_t pi=acos(-1.);
749 dvt2=pow((sin(a[1])*ea[0]),2)+pow((a[0]*cos(a[1])*ea[1]),2);
758 v.SetVector(a,"sph");
759 v.SetErrors(ea,"sph");
763 ///////////////////////////////////////////////////////////////////////////
764 Ali3Vector Ali3Vector::GetVecLong()
766 // Provide the longitudinal vector w.r.t. z-axis.
767 // Error propagation is performed automatically
768 Double_t pi=acos(-1.);
776 dvl2=pow((cos(a[1])*ea[0]),2)+pow((a[0]*sin(a[1])*ea[1]),2);
788 v.SetVector(a,"sph");
789 v.SetErrors(ea,"sph");
793 ///////////////////////////////////////////////////////////////////////////