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::Load(Ali3Vector& q)
110 // Load all attributes of the input Ali3Vector into this Ali3Vector object.
111 Double_t temp=q.GetResultError();
113 q.GetVector(a,"sph");
115 q.GetErrors(a,"car");
119 ///////////////////////////////////////////////////////////////////////////
120 void Ali3Vector::SetZero()
122 // (Re)set all attributes to zero.
131 ///////////////////////////////////////////////////////////////////////////
132 void Ali3Vector::SetVector(Double_t* v,TString f)
134 // Store vector according to reference frame f
135 // All errors will be reset to 0
141 Double_t pi=acos(-1.);
144 if (f == "car") frame=1;
145 if (f == "sph") frame=2;
146 if (f == "cyl") frame=3;
148 Double_t x,y,z,rho,phi;
152 case 1: // Cartesian coordinates
156 fV=sqrt(x*x+y*y+z*z);
158 if (fV && fabs(z/fV)<=1.)
166 if (fTheta<0.) fTheta+=2.*pi;
168 if (x || y) fPhi=atan2(y,x);
169 if (fPhi<0.) fPhi+=2.*pi;
172 case 2: // Spherical coordinates
178 case 3: // Cylindrical coordinates
182 fV=sqrt(rho*rho+z*z);
184 if (fPhi<0.) fPhi+=2.*pi;
186 if (fV && fabs(z/fV)<=1.)
194 if (fTheta<0.) fTheta+=2.*pi;
197 default: // Unsupported reference frame
198 cout << "*Ali3Vector::SetVector* Unsupported frame : " << f.Data() << endl
199 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
206 ///////////////////////////////////////////////////////////////////////////
207 void Ali3Vector::GetVector(Double_t* v,TString f) const
209 // Provide vector according to reference frame f
211 if (f == "car") frame=1;
212 if (f == "sph") frame=2;
213 if (f == "cyl") frame=3;
217 case 1: // Cartesian coordinates
218 v[0]=fV*sin(fTheta)*cos(fPhi);
219 v[1]=fV*sin(fTheta)*sin(fPhi);
223 case 2: // Spherical coordinates
229 case 3: // Cylindrical coordinates
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++)
245 ///////////////////////////////////////////////////////////////////////////
246 void Ali3Vector::SetVector(Float_t* v,TString f)
248 // Store vector according to reference frame f
249 // All errors will be reset to 0
251 for (Int_t i=0; i<3; i++)
257 ///////////////////////////////////////////////////////////////////////////
258 void Ali3Vector::GetVector(Float_t* v,TString f) const
260 // Provide vector according to reference frame f
263 for (Int_t i=0; i<3; i++)
268 ///////////////////////////////////////////////////////////////////////////
269 void Ali3Vector::SetErrors(Double_t* e,TString f)
271 // Store errors according to reference frame f
272 // The error on scalar results is reset to 0
276 if (f == "car") frame=1;
277 if (f == "sph") frame=2;
278 if (f == "cyl") frame=3;
280 Double_t dx2,dy2,dz2,rho;
284 case 1: // Cartesian coordinates
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);
301 case 3: // Cylindrical coordinates
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);
310 default: // Unsupported reference frame
311 cout << "*Ali3Vector::SetErrors* Unsupported frame : " << f.Data() << endl
312 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
319 ///////////////////////////////////////////////////////////////////////////
320 void Ali3Vector::GetErrors(Double_t* e,TString f) const
322 // Provide errors according to reference frame f
324 if (f == "car") frame=1;
325 if (f == "sph") frame=2;
326 if (f == "cyl") frame=3;
328 Double_t pi=acos(-1.);
330 Double_t dr2,dtheta2,dphi2,rho,drho2;
332 Double_t rxy2; // Shorthand for (x*x+y*y)
336 case 1: // Cartesian coordinates
342 case 2: // Spherical coordinates
344 rxy2=pow(v[0],2)+pow(v[1],2);
345 if (sqrt(rxy2)<(fV*1e-10)) rxy2=0;
348 dr2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)+pow((v[2]*fDz),2))/(fV*fV);
356 dtheta2=rxy2*pow(fDz,2)/pow(fV,4);
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);
369 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rxy2,2));
377 if (e[1]>pi) e[1]=pi;
379 if (e[2]>(2.*pi)) e[2]=2.*pi;
382 case 3: // Cylindrical coordinates
384 rho=fabs(fV*sin(fTheta));
385 if (rho<(fV*1e-10)) rho=0;
388 drho2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2))/(rho*rho);
396 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rho,4));
404 if (e[1]>(2.*pi)) e[1]=2.*pi;
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++)
418 ///////////////////////////////////////////////////////////////////////////
419 void Ali3Vector::SetErrors(Float_t* e,TString f)
421 // Store errors according to reference frame f
422 // The error on scalar results is reset to 0
424 for (Int_t i=0; i<3; i++)
430 ///////////////////////////////////////////////////////////////////////////
431 void Ali3Vector::GetErrors(Float_t* e,TString f) const
433 // Provide errors according to reference frame f
436 for (Int_t i=0; i<3; i++)
441 ///////////////////////////////////////////////////////////////////////////
442 void Ali3Vector::Data(TString f) const
444 // Print vector components according to reference frame f
445 if (f=="car" || f=="sph" || f=="cyl")
447 Double_t vec[3],err[3];
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;
457 cout << " *Ali3Vector::Data* Unsupported frame : " << f.Data() << endl
458 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
461 ///////////////////////////////////////////////////////////////////////////
462 Double_t Ali3Vector::GetNorm()
464 // Provide the norm of the current vector
465 // The error on the scalar result (norm) is updated accordingly
471 ///////////////////////////////////////////////////////////////////////////
472 Double_t Ali3Vector::GetPseudoRapidity()
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.);
480 Double_t thetahalf=v[1]/2.;
482 if (v[1]<pi) arg=tan(thetahalf);
484 if (arg>0) eta=-log(arg);
487 Double_t prod=cos(thetahalf)*sin(thetahalf);
489 if (prod) fDresult=fabs(e[1]/2.*prod);
492 ///////////////////////////////////////////////////////////////////////////
493 Double_t Ali3Vector::Dot(Ali3Vector& q)
495 // Provide the dot product of the current vector with vector q
496 // The error on the scalar result (dotproduct) is updated accordingly
500 if ((this) == &q) // Check for special case v.Dot(v)
502 Double_t norm=GetNorm();
503 Double_t dnorm=GetResultError();
505 fDresult=2.*norm*dnorm;
510 Double_t ea[3],eb[3];
515 q.GetVector(b,"car");
516 q.GetErrors(eb,"car");
517 for (Int_t i=0; i<3; i++)
520 d2+=pow(b[i]*ea[i],2)+pow(a[i]*eb[i],2);
527 ///////////////////////////////////////////////////////////////////////////
528 Double_t Ali3Vector::GetResultError() const
530 // Provide the error on the result of an operation yielding a scalar
531 // E.g. GetNorm() or Dot()
534 ///////////////////////////////////////////////////////////////////////////
535 Ali3Vector Ali3Vector::Cross(Ali3Vector& q) const
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;
544 q.GetVector(b,"car");
545 q.GetErrors(eb,"car");
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];
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);
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);
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);
564 v.SetVector(c,"car");
565 v.SetErrors(ec,"car");
569 ///////////////////////////////////////////////////////////////////////////
570 Ali3Vector Ali3Vector::operator+(Ali3Vector& q) const
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];
578 q.GetVector(b,"car");
579 q.GetErrors(eb,"car");
581 for (Int_t i=0; i<3; i++)
584 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
588 v.SetVector(a,"car");
589 v.SetErrors(ea,"car");
593 ///////////////////////////////////////////////////////////////////////////
594 Ali3Vector Ali3Vector::operator-(Ali3Vector& q) const
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];
602 q.GetVector(b,"car");
603 q.GetErrors(eb,"car");
605 for (Int_t i=0; i<3; i++)
608 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
612 v.SetVector(a,"car");
613 v.SetErrors(ea,"car");
617 ///////////////////////////////////////////////////////////////////////////
618 Ali3Vector Ali3Vector::operator*(Double_t s) const
620 // Multiply the current vector with a scalar s.
621 // Error propagation is performed automatically.
627 for (Int_t i=0; i<3; i++)
634 v.SetVector(a,"car");
635 v.SetErrors(ea,"car");
639 ///////////////////////////////////////////////////////////////////////////
640 Ali3Vector Ali3Vector::operator/(Double_t s) const
642 // Divide the current vector by a scalar s
643 // Error propagation is performed automatically
645 if (fabs(s)<1.e-20) // Protect against division by 0
647 cout << " *Ali3Vector::/* Division by 0 detected. No action taken." << endl;
657 for (Int_t i=0; i<3; i++)
664 v.SetVector(a,"car");
665 v.SetErrors(ea,"car");
670 ///////////////////////////////////////////////////////////////////////////
671 Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
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];
679 q.GetVector(b,"car");
680 q.GetErrors(eb,"car");
682 for (Int_t i=0; i<3; i++)
685 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
693 ///////////////////////////////////////////////////////////////////////////
694 Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
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];
702 q.GetVector(b,"car");
703 q.GetErrors(eb,"car");
705 for (Int_t i=0; i<3; i++)
708 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
716 ///////////////////////////////////////////////////////////////////////////
717 Ali3Vector& Ali3Vector::operator*=(Double_t s)
719 // Multiply the current vector with a scalar s
720 // Error propagation is performed automatically
726 for (Int_t i=0; i<3; i++)
737 ///////////////////////////////////////////////////////////////////////////
738 Ali3Vector& Ali3Vector::operator/=(Double_t s)
740 // Divide the current vector by a scalar s
741 // Error propagation is performed automatically
743 if (fabs(s)<1.e-20) // Protect against division by 0
745 cout << " *Ali3Vector::/=* Division by 0 detected. No action taken." << endl;
755 for (Int_t i=0; i<3; i++)
767 ///////////////////////////////////////////////////////////////////////////
768 Ali3Vector Ali3Vector::GetVecTrans() const
770 // Provide the transverse vector w.r.t. z-axis.
771 // Error propagation is performed automatically
772 Double_t pi=acos(-1.);
780 dvt2=pow((sin(a[1])*ea[0]),2)+pow((a[0]*cos(a[1])*ea[1]),2);
789 v.SetVector(a,"sph");
790 v.SetErrors(ea,"sph");
794 ///////////////////////////////////////////////////////////////////////////
795 Ali3Vector Ali3Vector::GetVecLong() const
797 // Provide the longitudinal vector w.r.t. z-axis.
798 // Error propagation is performed automatically
799 Double_t pi=acos(-1.);
807 dvl2=pow((cos(a[1])*ea[0]),2)+pow((a[0]*sin(a[1])*ea[1]),2);
819 v.SetVector(a,"sph");
820 v.SetErrors(ea,"sph");
824 ///////////////////////////////////////////////////////////////////////////
825 Ali3Vector Ali3Vector::GetPrimed(TRotMatrix* m) const
827 // Provide vector components (and errors) in a rotated frame.
828 // The orientation of the rotated frame is described by the TRotMatrix
833 Double_t* mat=m->GetMatrix();
835 Double_t a[3],aprim[3];
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");
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");
851 ///////////////////////////////////////////////////////////////////////////
852 Ali3Vector Ali3Vector::GetUnprimed(TRotMatrix* m) const
854 // Provide original vector components (and errors) from the rotated ones.
855 // The orientation of the rotated frame is described by the TRotMatrix
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.
863 Double_t* mat=m->GetMatrix();
865 Double_t a[3],aprim[3];
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");
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");
881 ///////////////////////////////////////////////////////////////////////////