Removing warnings
[u/mrichter/AliRoot.git] / RALICE / Ali3Vector.cxx
CommitLineData
4c039060 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
f531a546 16// $Id$
4c039060 17
959fbac5 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// ------
1f241680 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)
959fbac5 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//
1f241680 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.
959fbac5 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");
84bb7c66 54// a.Data();
959fbac5 55//
56// Float_t vec[3];
57// Float_t err[3];
1f241680 58// a.GetVector(vec,"sph","deg");
59// a.GetErrors(vec,"sph","deg");
959fbac5 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);
84bb7c66 71// c.Data("sph");
959fbac5 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
f531a546 83//- Modified: NvE $Date$ UU-SAP Utrecht
959fbac5 84///////////////////////////////////////////////////////////////////////////
85
d88f97cc 86#include "Ali3Vector.h"
c72198f1 87#include "Riostream.h"
d88f97cc 88
89ClassImp(Ali3Vector) // Class implementation to enable ROOT I/O
90
91Ali3Vector::Ali3Vector()
92{
93// Creation of an Ali3Vector object and initialisation of parameters
959fbac5 94// All attributes initialised to 0
c72198f1 95 SetZero();
96}
97///////////////////////////////////////////////////////////////////////////
98Ali3Vector::~Ali3Vector()
99{
100// Destructor to delete dynamically allocated memory
101}
102///////////////////////////////////////////////////////////////////////////
103Ali3Vector::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///////////////////////////////////////////////////////////////////////////
1fbffa23 115void 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///////////////////////////////////////////////////////////////////////////
c72198f1 127void Ali3Vector::SetZero()
128{
129// (Re)set all attributes to zero.
d88f97cc 130 fV=0;
131 fTheta=0;
132 fPhi=0;
959fbac5 133 fDx=0;
134 fDy=0;
135 fDz=0;
136 fDresult=0;
d88f97cc 137}
138///////////////////////////////////////////////////////////////////////////
1f241680 139void Ali3Vector::SetVector(Double_t* v,TString f,TString u)
d88f97cc 140{
141// Store vector according to reference frame f
1f241680 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//
959fbac5 150// All errors will be reset to 0
1f241680 151
959fbac5 152 fDx=0;
153 fDy=0;
154 fDz=0;
155 fDresult=0;
156
d88f97cc 157 Double_t pi=acos(-1.);
959fbac5 158
1f241680 159 Double_t fu=1.;
160 if (u == "deg") fu=pi/180.;
161
d88f97cc 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];
1f241680 193 fTheta=v[1]*fu;
194 fPhi=v[2]*fu;
d88f97cc 195 break;
196
197 case 3: // Cylindrical coordinates
198 rho=v[0];
1f241680 199 phi=v[1]*fu;
d88f97cc 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
c72198f1 217 cout << "*Ali3Vector::SetVector* Unsupported frame : " << f.Data() << endl
d88f97cc 218 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
219 fV=0;
220 fTheta=0;
221 fPhi=0;
222 break;
223 }
224}
225///////////////////////////////////////////////////////////////////////////
1f241680 226void Ali3Vector::GetVector(Double_t* v,TString f,TString u) const
d88f97cc 227{
228// Provide vector according to reference frame f
1f241680 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
d88f97cc 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;
1f241680 257 v[1]=fTheta*fu;
258 v[2]=fPhi*fu;
d88f97cc 259 break;
260
261 case 3: // Cylindrical coordinates
262 v[0]=fV*sin(fTheta);
1f241680 263 v[1]=fPhi*fu;
d88f97cc 264 v[2]=fV*cos(fTheta);
265 break;
266
267 default: // Unsupported reference frame
c72198f1 268 cout << "*Ali3Vector::GetVector* Unsupported frame : " << f.Data() << endl
d88f97cc 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///////////////////////////////////////////////////////////////////////////
1f241680 278void Ali3Vector::SetVector(Float_t* v,TString f,TString u)
d88f97cc 279{
280// Store vector according to reference frame f
1f241680 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//
959fbac5 289// All errors will be reset to 0
1f241680 290
d88f97cc 291 Double_t vec[3];
292 for (Int_t i=0; i<3; i++)
293 {
294 vec[i]=v[i];
295 }
1f241680 296 SetVector(vec,f,u);
d88f97cc 297}
298///////////////////////////////////////////////////////////////////////////
1f241680 299void Ali3Vector::GetVector(Float_t* v,TString f,TString u) const
d88f97cc 300{
301// Provide vector according to reference frame f
1f241680 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
d88f97cc 310 Double_t vec[3];
1f241680 311 GetVector(vec,f,u);
d88f97cc 312 for (Int_t i=0; i<3; i++)
313 {
314 v[i]=vec[i];
315 }
316}
317///////////////////////////////////////////////////////////////////////////
1f241680 318void Ali3Vector::SetErrors(Double_t* e,TString f,TString u)
959fbac5 319{
320// Store errors according to reference frame f
1f241680 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//
959fbac5 329// The error on scalar results is reset to 0
1f241680 330
331 Double_t pi=acos(-1.);
332
333 Double_t fu=1.;
334 if (u == "deg") fu=pi/180.;
335
959fbac5 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
1f241680 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);
959fbac5 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);
1f241680 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);
959fbac5 368 fDx=sqrt(dx2);
369 fDy=sqrt(dy2);
370 fDz=fabs(e[2]);
371 break;
372
373 default: // Unsupported reference frame
c72198f1 374 cout << "*Ali3Vector::SetErrors* Unsupported frame : " << f.Data() << endl
959fbac5 375 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
376 fDx=0;
377 fDy=0;
378 fDz=0;
379 break;
380 }
381}
382///////////////////////////////////////////////////////////////////////////
1f241680 383void Ali3Vector::GetErrors(Double_t* e,TString f,TString u) const
959fbac5 384{
385// Provide errors according to reference frame f
1f241680 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
959fbac5 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];
d071d629 406 Double_t rxy2; // Shorthand for (x*x+y*y)
959fbac5 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");
d071d629 418 rxy2=pow(v[0],2)+pow(v[1],2);
419 if (sqrt(rxy2)<(fV*1e-10)) rxy2=0;
959fbac5 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 }
d071d629 428 if (fV)
959fbac5 429 {
d071d629 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)) /
81508922 434 pow(((pow(v[2],2)*rxy2)+pow(rxy2,2)),2);
d071d629 435 }
959fbac5 436 }
437 else
438 {
959fbac5 439 dtheta2=0;
440 }
d071d629 441 if (rxy2)
959fbac5 442 {
d071d629 443 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rxy2,2));
959fbac5 444 }
445 else
446 {
447 dphi2=0;
448 }
449 e[0]=sqrt(dr2);
450 e[1]=sqrt(dtheta2);
d071d629 451 if (e[1]>pi) e[1]=pi;
959fbac5 452 e[2]=sqrt(dphi2);
d071d629 453 if (e[2]>(2.*pi)) e[2]=2.*pi;
1f241680 454 e[1]*=fu;
455 e[2]*=fu;
959fbac5 456 break;
457
458 case 3: // Cylindrical coordinates
459 GetVector(v,"car");
d071d629 460 rho=fabs(fV*sin(fTheta));
461 if (rho<(fV*1e-10)) rho=0;
959fbac5 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 }
d071d629 470 if (rho)
959fbac5 471 {
d071d629 472 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rho,4));
959fbac5 473 }
474 else
475 {
476 dphi2=0;
477 }
478 e[0]=sqrt(drho2);
479 e[1]=sqrt(dphi2);
d071d629 480 if (e[1]>(2.*pi)) e[1]=2.*pi;
959fbac5 481 e[2]=fDz;
1f241680 482 e[1]*=fu;
959fbac5 483 break;
484
485 default: // Unsupported reference frame
c72198f1 486 cout << "*Ali3Vector::GetErrors* Unsupported frame : " << f.Data() << endl
959fbac5 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///////////////////////////////////////////////////////////////////////////
1f241680 496void Ali3Vector::SetErrors(Float_t* e,TString f,TString u)
959fbac5 497{
498// Store errors according to reference frame f
1f241680 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//
959fbac5 507// The error on scalar results is reset to 0
1f241680 508
959fbac5 509 Double_t vec[3];
510 for (Int_t i=0; i<3; i++)
511 {
512 vec[i]=e[i];
513 }
1f241680 514 SetErrors(vec,f,u);
959fbac5 515}
516///////////////////////////////////////////////////////////////////////////
1f241680 517void Ali3Vector::GetErrors(Float_t* e,TString f,TString u) const
959fbac5 518{
519// Provide errors according to reference frame f
1f241680 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
959fbac5 528 Double_t vec[3];
1f241680 529 GetErrors(vec,f,u);
959fbac5 530 for (Int_t i=0; i<3; i++)
531 {
532 e[i]=vec[i];
533 }
534}
535///////////////////////////////////////////////////////////////////////////
1f241680 536void Ali3Vector::Data(TString f,TString u) const
d88f97cc 537{
538// Print vector components according to reference frame f
1f241680 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
d88f97cc 547 if (f=="car" || f=="sph" || f=="cyl")
548 {
959fbac5 549 Double_t vec[3],err[3];
1f241680 550 GetVector(vec,f,u);
551 GetErrors(err,f,u);
552 cout << " Vector in " << f.Data() << " (" << u.Data() << ") coordinates : "
d88f97cc 553 << vec[0] << " " << vec[1] << " " << vec[2] << endl;
1f241680 554 cout << " Err. in " << f.Data() << " (" << u.Data() << ") coordinates : "
959fbac5 555 << err[0] << " " << err[1] << " " << err[2] << endl;
d88f97cc 556 }
557 else
558 {
c72198f1 559 cout << " *Ali3Vector::Data* Unsupported frame : " << f.Data() << endl
d88f97cc 560 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
561 }
562}
563///////////////////////////////////////////////////////////////////////////
564Double_t Ali3Vector::GetNorm()
565{
959fbac5 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];
d88f97cc 571 return fV;
572}
573///////////////////////////////////////////////////////////////////////////
959fbac5 574Double_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
8adaf597 579 Double_t pi=acos(-1.);
959fbac5 580 Double_t v[3];
581 GetVector(v,"sph");
582 Double_t thetahalf=v[1]/2.;
8adaf597 583 Double_t arg=0;
584 if (v[1]<pi) arg=tan(thetahalf);
585 Double_t eta=9999;
959fbac5 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///////////////////////////////////////////////////////////////////////////
d88f97cc 595Double_t Ali3Vector::Dot(Ali3Vector& q)
596{
597// Provide the dot product of the current vector with vector q
959fbac5 598// The error on the scalar result (dotproduct) is updated accordingly
599
d88f97cc 600 Double_t dotpro=0;
601
959fbac5 602 if ((this) == &q) // Check for special case v.Dot(v)
d88f97cc 603 {
959fbac5 604 Double_t norm=GetNorm();
605 Double_t dnorm=GetResultError();
606 dotpro=pow(norm,2);
607 fDresult=2.*norm*dnorm;
d88f97cc 608 }
959fbac5 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
d88f97cc 627 return dotpro;
628}
629///////////////////////////////////////////////////////////////////////////
261c0caf 630Double_t Ali3Vector::GetResultError() const
959fbac5 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///////////////////////////////////////////////////////////////////////////
261c0caf 637Ali3Vector Ali3Vector::Cross(Ali3Vector& q) const
d88f97cc 638{
639// Provide the cross product of the current vector with vector q
959fbac5 640// Error propagation is performed automatically
d88f97cc 641 Double_t a[3],b[3],c[3];
959fbac5 642 Double_t ea[3],eb[3],ec[3],d2;
d88f97cc 643
644 GetVector(a,"car");
959fbac5 645 GetErrors(ea,"car");
d88f97cc 646 q.GetVector(b,"car");
959fbac5 647 q.GetErrors(eb,"car");
d88f97cc 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
959fbac5 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
d88f97cc 665 Ali3Vector v;
666 v.SetVector(c,"car");
959fbac5 667 v.SetErrors(ec,"car");
d88f97cc 668
669 return v;
670}
671///////////////////////////////////////////////////////////////////////////
261c0caf 672Ali3Vector Ali3Vector::operator+(Ali3Vector& q) const
d88f97cc 673{
674// Add vector q to the current vector
959fbac5 675// Error propagation is performed automatically
676 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 677
678 GetVector(a,"car");
959fbac5 679 GetErrors(ea,"car");
d88f97cc 680 q.GetVector(b,"car");
959fbac5 681 q.GetErrors(eb,"car");
d88f97cc 682
683 for (Int_t i=0; i<3; i++)
684 {
685 a[i]+=b[i];
959fbac5 686 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 687 }
688
689 Ali3Vector v;
690 v.SetVector(a,"car");
959fbac5 691 v.SetErrors(ea,"car");
d88f97cc 692
693 return v;
694}
695///////////////////////////////////////////////////////////////////////////
261c0caf 696Ali3Vector Ali3Vector::operator-(Ali3Vector& q) const
d88f97cc 697{
698// Subtract vector q from the current vector
959fbac5 699// Error propagation is performed automatically
700 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 701
702 GetVector(a,"car");
959fbac5 703 GetErrors(ea,"car");
d88f97cc 704 q.GetVector(b,"car");
959fbac5 705 q.GetErrors(eb,"car");
d88f97cc 706
707 for (Int_t i=0; i<3; i++)
708 {
709 a[i]-=b[i];
959fbac5 710 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 711 }
712
713 Ali3Vector v;
714 v.SetVector(a,"car");
959fbac5 715 v.SetErrors(ea,"car");
d88f97cc 716
717 return v;
718}
719///////////////////////////////////////////////////////////////////////////
261c0caf 720Ali3Vector Ali3Vector::operator*(Double_t s) const
d88f97cc 721{
959fbac5 722// Multiply the current vector with a scalar s.
723// Error propagation is performed automatically.
724 Double_t a[3],ea[3];
d88f97cc 725
726 GetVector(a,"car");
959fbac5 727 GetErrors(ea,"car");
d88f97cc 728
729 for (Int_t i=0; i<3; i++)
730 {
731 a[i]*=s;
959fbac5 732 ea[i]*=s;
d88f97cc 733 }
734
735 Ali3Vector v;
736 v.SetVector(a,"car");
959fbac5 737 v.SetErrors(ea,"car");
d88f97cc 738
739 return v;
740}
741///////////////////////////////////////////////////////////////////////////
261c0caf 742Ali3Vector Ali3Vector::operator/(Double_t s) const
d88f97cc 743{
744// Divide the current vector by a scalar s
959fbac5 745// Error propagation is performed automatically
d88f97cc 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 {
959fbac5 754 Double_t a[3],ea[3];
d88f97cc 755
756 GetVector(a,"car");
959fbac5 757 GetErrors(ea,"car");
d88f97cc 758
759 for (Int_t i=0; i<3; i++)
760 {
761 a[i]/=s;
959fbac5 762 ea[i]/=s;
d88f97cc 763 }
764
765 Ali3Vector v;
766 v.SetVector(a,"car");
959fbac5 767 v.SetErrors(ea,"car");
d88f97cc 768
769 return v;
770 }
771}
772///////////////////////////////////////////////////////////////////////////
773Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
774{
775// Add vector q to the current vector
959fbac5 776// Error propagation is performed automatically
777 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 778
779 GetVector(a,"car");
959fbac5 780 GetErrors(ea,"car");
d88f97cc 781 q.GetVector(b,"car");
959fbac5 782 q.GetErrors(eb,"car");
d88f97cc 783
784 for (Int_t i=0; i<3; i++)
785 {
786 a[i]+=b[i];
959fbac5 787 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 788 }
789
790 SetVector(a,"car");
959fbac5 791 SetErrors(ea,"car");
d88f97cc 792
793 return *this;
794}
795///////////////////////////////////////////////////////////////////////////
796Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
797{
798// Subtract vector q from the current vector
959fbac5 799// Error propagation is performed automatically
800 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 801
802 GetVector(a,"car");
959fbac5 803 GetErrors(ea,"car");
d88f97cc 804 q.GetVector(b,"car");
959fbac5 805 q.GetErrors(eb,"car");
d88f97cc 806
807 for (Int_t i=0; i<3; i++)
808 {
809 a[i]-=b[i];
959fbac5 810 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 811 }
812
813 SetVector(a,"car");
959fbac5 814 SetErrors(ea,"car");
d88f97cc 815
816 return *this;
817}
818///////////////////////////////////////////////////////////////////////////
819Ali3Vector& Ali3Vector::operator*=(Double_t s)
820{
821// Multiply the current vector with a scalar s
959fbac5 822// Error propagation is performed automatically
823 Double_t a[3],ea[3];
d88f97cc 824
825 GetVector(a,"car");
959fbac5 826 GetErrors(ea,"car");
d88f97cc 827
828 for (Int_t i=0; i<3; i++)
829 {
830 a[i]*=s;
959fbac5 831 ea[i]*=s;
d88f97cc 832 }
833
834 SetVector(a,"car");
959fbac5 835 SetErrors(ea,"car");
d88f97cc 836
837 return *this;
838}
839///////////////////////////////////////////////////////////////////////////
840Ali3Vector& Ali3Vector::operator/=(Double_t s)
841{
842// Divide the current vector by a scalar s
959fbac5 843// Error propagation is performed automatically
d88f97cc 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 {
959fbac5 852 Double_t a[3],ea[3];
d88f97cc 853
854 GetVector(a,"car");
959fbac5 855 GetErrors(ea,"car");
d88f97cc 856
857 for (Int_t i=0; i<3; i++)
858 {
859 a[i]/=s;
959fbac5 860 ea[i]/=s;
d88f97cc 861 }
862
863 SetVector(a,"car");
959fbac5 864 SetErrors(ea,"car");
d88f97cc 865
866 return *this;
867 }
868}
869///////////////////////////////////////////////////////////////////////////
261c0caf 870Ali3Vector Ali3Vector::GetVecTrans() const
d071d629 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///////////////////////////////////////////////////////////////////////////
261c0caf 897Ali3Vector Ali3Vector::GetVecLong() const
d071d629 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///////////////////////////////////////////////////////////////////////////
261c0caf 927Ali3Vector 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///////////////////////////////////////////////////////////////////////////
c5555bc0 954Ali3Vector 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///////////////////////////////////////////////////////////////////////////
1f241680 984Double_t Ali3Vector::GetX(Int_t i,TString f,TString u)
a7dc0627 985{
986// Provide i-th vector component according to reference frame f.
1f241680 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//
a7dc0627 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];
1f241680 1006 GetVector(vec,f,u);
1007 GetErrors(err,f,u);
a7dc0627 1008
1009 fDresult=err[i-1];
1010
1011 return vec[i-1];
1012}
1013///////////////////////////////////////////////////////////////////////////
1f241680 1014Double_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;
35e721bf 1048 if (x>1.) x=1;
1049 if (x<-1.) x=-1;
1f241680 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///////////////////////////////////////////////////////////////////////////