08-mar-2003 NvE Compiler option /GR introduced for MSVC++ in mklibs.bat to explicitly...
[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// ------
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 :
32//
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)
36//
37// All angles are in radians.
38//
39// Example :
40// ---------
41//
42// Ali3Vector a;
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");
84bb7c66 47// a.Data();
959fbac5 48//
49// Float_t vec[3];
50// Float_t err[3];
51// a.GetVector(vec,"sph");
52// a.GetErrors(vec,"sph");
53//
54// Ali3Vector b;
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");
59//
60// Float_t dotpro=a.Dot(b);
61// Float_t doterror=a.GetResultError();
62//
63// Ali3Vector c=a.Cross(b);
84bb7c66 64// c.Data("sph");
959fbac5 65// c.GetVector(vec,"cyl");
66// c.GetErrors(err,"cyl");
67//
68// Float_t norm=c.GetNorm();
69// Float_t normerror=c.GetResultError();
70//
71// c=a+b;
72// c=a-b;
73// c=a*5;
74//
75//--- Author: Nick van Eijndhoven 30-mar-1999 UU-SAP Utrecht
f531a546 76//- Modified: NvE $Date$ UU-SAP Utrecht
959fbac5 77///////////////////////////////////////////////////////////////////////////
78
d88f97cc 79#include "Ali3Vector.h"
c72198f1 80#include "Riostream.h"
d88f97cc 81
82ClassImp(Ali3Vector) // Class implementation to enable ROOT I/O
83
84Ali3Vector::Ali3Vector()
85{
86// Creation of an Ali3Vector object and initialisation of parameters
959fbac5 87// All attributes initialised to 0
c72198f1 88 SetZero();
89}
90///////////////////////////////////////////////////////////////////////////
91Ali3Vector::~Ali3Vector()
92{
93// Destructor to delete dynamically allocated memory
94}
95///////////////////////////////////////////////////////////////////////////
96Ali3Vector::Ali3Vector(const Ali3Vector& v)
97{
98// Copy constructor
99 fV=v.fV;
100 fTheta=v.fTheta;
101 fPhi=v.fPhi;
102 fDx=v.fDx;
103 fDy=v.fDy;
104 fDz=v.fDz;
105 fDresult=v.fDresult;
106}
107///////////////////////////////////////////////////////////////////////////
108void Ali3Vector::SetZero()
109{
110// (Re)set all attributes to zero.
d88f97cc 111 fV=0;
112 fTheta=0;
113 fPhi=0;
959fbac5 114 fDx=0;
115 fDy=0;
116 fDz=0;
117 fDresult=0;
d88f97cc 118}
119///////////////////////////////////////////////////////////////////////////
d88f97cc 120void Ali3Vector::SetVector(Double_t* v,TString f)
121{
122// Store vector according to reference frame f
959fbac5 123// All errors will be reset to 0
124 fDx=0;
125 fDy=0;
126 fDz=0;
127 fDresult=0;
128
d88f97cc 129 Double_t pi=acos(-1.);
959fbac5 130
d88f97cc 131 Int_t frame=0;
132 if (f == "car") frame=1;
133 if (f == "sph") frame=2;
134 if (f == "cyl") frame=3;
135
136 Double_t x,y,z,rho,phi;
137
138 switch (frame)
139 {
140 case 1: // Cartesian coordinates
141 x=v[0];
142 y=v[1];
143 z=v[2];
144 fV=sqrt(x*x+y*y+z*z);
145 fTheta=0;
146 if (fV && fabs(z/fV)<=1.)
147 {
148 fTheta=acos(z/fV);
149 }
150 else
151 {
152 if (z<0.) fTheta=pi;
153 }
154 if (fTheta<0.) fTheta+=2.*pi;
155 fPhi=0;
156 if (x || y) fPhi=atan2(y,x);
157 if (fPhi<0.) fPhi+=2.*pi;
158 break;
159
160 case 2: // Spherical coordinates
161 fV=v[0];
162 fTheta=v[1];
163 fPhi=v[2];
164 break;
165
166 case 3: // Cylindrical coordinates
167 rho=v[0];
168 phi=v[1];
169 z=v[2];
170 fV=sqrt(rho*rho+z*z);
171 fPhi=phi;
172 if (fPhi<0.) fPhi+=2.*pi;
173 fTheta=0;
174 if (fV && fabs(z/fV)<=1.)
175 {
176 fTheta=acos(z/fV);
177 }
178 else
179 {
180 if (z<0.) fTheta=pi;
181 }
182 if (fTheta<0.) fTheta+=2.*pi;
183 break;
184
185 default: // Unsupported reference frame
c72198f1 186 cout << "*Ali3Vector::SetVector* Unsupported frame : " << f.Data() << endl
d88f97cc 187 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
188 fV=0;
189 fTheta=0;
190 fPhi=0;
191 break;
192 }
193}
194///////////////////////////////////////////////////////////////////////////
195void Ali3Vector::GetVector(Double_t* v,TString f)
196{
197// Provide vector according to reference frame f
198 Int_t frame=0;
199 if (f == "car") frame=1;
200 if (f == "sph") frame=2;
201 if (f == "cyl") frame=3;
202
203 switch (frame)
204 {
205 case 1: // Cartesian coordinates
206 v[0]=fV*sin(fTheta)*cos(fPhi);
207 v[1]=fV*sin(fTheta)*sin(fPhi);
208 v[2]=fV*cos(fTheta);
209 break;
210
211 case 2: // Spherical coordinates
212 v[0]=fV;
213 v[1]=fTheta;
214 v[2]=fPhi;
215 break;
216
217 case 3: // Cylindrical coordinates
218 v[0]=fV*sin(fTheta);
219 v[1]=fPhi;
220 v[2]=fV*cos(fTheta);
221 break;
222
223 default: // Unsupported reference frame
c72198f1 224 cout << "*Ali3Vector::GetVector* Unsupported frame : " << f.Data() << endl
d88f97cc 225 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
226 for (Int_t i=0; i<3; i++)
227 {
228 v[i]=0;
229 }
230 break;
231 }
232}
233///////////////////////////////////////////////////////////////////////////
234void Ali3Vector::SetVector(Float_t* v,TString f)
235{
236// Store vector according to reference frame f
959fbac5 237// All errors will be reset to 0
d88f97cc 238 Double_t vec[3];
239 for (Int_t i=0; i<3; i++)
240 {
241 vec[i]=v[i];
242 }
243 SetVector(vec,f);
244}
245///////////////////////////////////////////////////////////////////////////
246void Ali3Vector::GetVector(Float_t* v,TString f)
247{
248// Provide vector according to reference frame f
249 Double_t vec[3];
250 GetVector(vec,f);
251 for (Int_t i=0; i<3; i++)
252 {
253 v[i]=vec[i];
254 }
255}
256///////////////////////////////////////////////////////////////////////////
959fbac5 257void Ali3Vector::SetErrors(Double_t* e,TString f)
258{
259// Store errors according to reference frame f
260// The error on scalar results is reset to 0
261 fDresult=0;
262
263 Int_t frame=0;
264 if (f == "car") frame=1;
265 if (f == "sph") frame=2;
266 if (f == "cyl") frame=3;
267
268 Double_t dx2,dy2,dz2,rho;
269
270 switch (frame)
271 {
272 case 1: // Cartesian coordinates
273 fDx=fabs(e[0]);
274 fDy=fabs(e[1]);
275 fDz=fabs(e[2]);
276 break;
277
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);
284 fDx=sqrt(dx2);
285 fDy=sqrt(dy2);
286 fDz=sqrt(dz2);
287 break;
288
289 case 3: // Cylindrical coordinates
290 rho=fV*sin(fTheta);
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);
293 fDx=sqrt(dx2);
294 fDy=sqrt(dy2);
295 fDz=fabs(e[2]);
296 break;
297
298 default: // Unsupported reference frame
c72198f1 299 cout << "*Ali3Vector::SetErrors* Unsupported frame : " << f.Data() << endl
959fbac5 300 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
301 fDx=0;
302 fDy=0;
303 fDz=0;
304 break;
305 }
306}
307///////////////////////////////////////////////////////////////////////////
308void Ali3Vector::GetErrors(Double_t* e,TString f)
309{
310// Provide errors according to reference frame f
311 Int_t frame=0;
312 if (f == "car") frame=1;
313 if (f == "sph") frame=2;
314 if (f == "cyl") frame=3;
315
d071d629 316 Double_t pi=acos(-1.);
317
959fbac5 318 Double_t dr2,dtheta2,dphi2,rho,drho2;
319 Double_t v[3];
d071d629 320 Double_t rxy2; // Shorthand for (x*x+y*y)
959fbac5 321
322 switch (frame)
323 {
324 case 1: // Cartesian coordinates
325 e[0]=fDx;
326 e[1]=fDy;
327 e[2]=fDz;
328 break;
329
330 case 2: // Spherical coordinates
331 GetVector(v,"car");
d071d629 332 rxy2=pow(v[0],2)+pow(v[1],2);
333 if (sqrt(rxy2)<(fV*1e-10)) rxy2=0;
959fbac5 334 if (fV)
335 {
336 dr2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)+pow((v[2]*fDz),2))/(fV*fV);
337 }
338 else
339 {
340 dr2=0;
341 }
d071d629 342 if (fV)
959fbac5 343 {
d071d629 344 dtheta2=rxy2*pow(fDz,2)/pow(fV,4);
345 if (v[2] && rxy2)
346 {
347 dtheta2+=rxy2*pow(v[2],2)*(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)) /
81508922 348 pow(((pow(v[2],2)*rxy2)+pow(rxy2,2)),2);
d071d629 349 }
959fbac5 350 }
351 else
352 {
959fbac5 353 dtheta2=0;
354 }
d071d629 355 if (rxy2)
959fbac5 356 {
d071d629 357 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rxy2,2));
959fbac5 358 }
359 else
360 {
361 dphi2=0;
362 }
363 e[0]=sqrt(dr2);
364 e[1]=sqrt(dtheta2);
d071d629 365 if (e[1]>pi) e[1]=pi;
959fbac5 366 e[2]=sqrt(dphi2);
d071d629 367 if (e[2]>(2.*pi)) e[2]=2.*pi;
959fbac5 368 break;
369
370 case 3: // Cylindrical coordinates
371 GetVector(v,"car");
d071d629 372 rho=fabs(fV*sin(fTheta));
373 if (rho<(fV*1e-10)) rho=0;
959fbac5 374 if (rho)
375 {
376 drho2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2))/(rho*rho);
377 }
378 else
379 {
380 drho2=0;
381 }
d071d629 382 if (rho)
959fbac5 383 {
d071d629 384 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rho,4));
959fbac5 385 }
386 else
387 {
388 dphi2=0;
389 }
390 e[0]=sqrt(drho2);
391 e[1]=sqrt(dphi2);
d071d629 392 if (e[1]>(2.*pi)) e[1]=2.*pi;
959fbac5 393 e[2]=fDz;
394 break;
395
396 default: // Unsupported reference frame
c72198f1 397 cout << "*Ali3Vector::GetErrors* Unsupported frame : " << f.Data() << endl
959fbac5 398 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
399 for (Int_t i=0; i<3; i++)
400 {
401 e[i]=0;
402 }
403 break;
404 }
405}
406///////////////////////////////////////////////////////////////////////////
407void Ali3Vector::SetErrors(Float_t* e,TString f)
408{
409// Store errors according to reference frame f
410// The error on scalar results is reset to 0
411 Double_t vec[3];
412 for (Int_t i=0; i<3; i++)
413 {
414 vec[i]=e[i];
415 }
416 SetErrors(vec,f);
417}
418///////////////////////////////////////////////////////////////////////////
419void Ali3Vector::GetErrors(Float_t* e,TString f)
420{
421// Provide errors according to reference frame f
422 Double_t vec[3];
423 GetErrors(vec,f);
424 for (Int_t i=0; i<3; i++)
425 {
426 e[i]=vec[i];
427 }
428}
429///////////////////////////////////////////////////////////////////////////
84bb7c66 430void Ali3Vector::Data(TString f)
d88f97cc 431{
432// Print vector components according to reference frame f
433 if (f=="car" || f=="sph" || f=="cyl")
434 {
959fbac5 435 Double_t vec[3],err[3];
d88f97cc 436 GetVector(vec,f);
959fbac5 437 GetErrors(err,f);
c72198f1 438 cout << " Vector in " << f.Data() << " coordinates : "
d88f97cc 439 << vec[0] << " " << vec[1] << " " << vec[2] << endl;
c72198f1 440 cout << " Err. in " << f.Data() << " coordinates : "
959fbac5 441 << err[0] << " " << err[1] << " " << err[2] << endl;
d88f97cc 442 }
443 else
444 {
c72198f1 445 cout << " *Ali3Vector::Data* Unsupported frame : " << f.Data() << endl
d88f97cc 446 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
447 }
448}
449///////////////////////////////////////////////////////////////////////////
450Double_t Ali3Vector::GetNorm()
451{
959fbac5 452// Provide the norm of the current vector
453// The error on the scalar result (norm) is updated accordingly
454 Double_t e[3];
455 GetErrors(e,"sph");
456 fDresult=e[0];
d88f97cc 457 return fV;
458}
459///////////////////////////////////////////////////////////////////////////
959fbac5 460Double_t Ali3Vector::GetPseudoRapidity()
461{
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
8adaf597 465 Double_t pi=acos(-1.);
959fbac5 466 Double_t v[3];
467 GetVector(v,"sph");
468 Double_t thetahalf=v[1]/2.;
8adaf597 469 Double_t arg=0;
470 if (v[1]<pi) arg=tan(thetahalf);
471 Double_t eta=9999;
959fbac5 472 if (arg>0) eta=-log(arg);
473 Double_t e[3];
474 GetErrors(e,"sph");
475 Double_t prod=cos(thetahalf)*sin(thetahalf);
476 fDresult=0;
477 if (prod) fDresult=fabs(e[1]/2.*prod);
478 return eta;
479}
480///////////////////////////////////////////////////////////////////////////
d88f97cc 481Double_t Ali3Vector::Dot(Ali3Vector& q)
482{
483// Provide the dot product of the current vector with vector q
959fbac5 484// The error on the scalar result (dotproduct) is updated accordingly
485
d88f97cc 486 Double_t dotpro=0;
487
959fbac5 488 if ((this) == &q) // Check for special case v.Dot(v)
d88f97cc 489 {
959fbac5 490 Double_t norm=GetNorm();
491 Double_t dnorm=GetResultError();
492 dotpro=pow(norm,2);
493 fDresult=2.*norm*dnorm;
d88f97cc 494 }
959fbac5 495 else
496 {
497 Double_t a[3],b[3];
498 Double_t ea[3],eb[3];
499 Double_t d2=0;
500
501 GetVector(a,"car");
502 GetErrors(ea,"car");
503 q.GetVector(b,"car");
504 q.GetErrors(eb,"car");
505 for (Int_t i=0; i<3; i++)
506 {
507 dotpro+=a[i]*b[i];
508 d2+=pow(b[i]*ea[i],2)+pow(a[i]*eb[i],2);
509 }
510 fDresult=sqrt(d2);
511 }
512
d88f97cc 513 return dotpro;
514}
515///////////////////////////////////////////////////////////////////////////
959fbac5 516Double_t Ali3Vector::GetResultError()
517{
518// Provide the error on the result of an operation yielding a scalar
519// E.g. GetNorm() or Dot()
520 return fDresult;
521}
522///////////////////////////////////////////////////////////////////////////
d88f97cc 523Ali3Vector Ali3Vector::Cross(Ali3Vector& q)
524{
525// Provide the cross product of the current vector with vector q
959fbac5 526// Error propagation is performed automatically
d88f97cc 527 Double_t a[3],b[3],c[3];
959fbac5 528 Double_t ea[3],eb[3],ec[3],d2;
d88f97cc 529
530 GetVector(a,"car");
959fbac5 531 GetErrors(ea,"car");
d88f97cc 532 q.GetVector(b,"car");
959fbac5 533 q.GetErrors(eb,"car");
d88f97cc 534
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];
538
959fbac5 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);
541 ec[0]=sqrt(d2);
542
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);
545 ec[1]=sqrt(d2);
546
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);
549 ec[2]=sqrt(d2);
550
d88f97cc 551 Ali3Vector v;
552 v.SetVector(c,"car");
959fbac5 553 v.SetErrors(ec,"car");
d88f97cc 554
555 return v;
556}
557///////////////////////////////////////////////////////////////////////////
558Ali3Vector Ali3Vector::operator+(Ali3Vector& q)
559{
560// Add vector q to the current vector
959fbac5 561// Error propagation is performed automatically
562 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 563
564 GetVector(a,"car");
959fbac5 565 GetErrors(ea,"car");
d88f97cc 566 q.GetVector(b,"car");
959fbac5 567 q.GetErrors(eb,"car");
d88f97cc 568
569 for (Int_t i=0; i<3; i++)
570 {
571 a[i]+=b[i];
959fbac5 572 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 573 }
574
575 Ali3Vector v;
576 v.SetVector(a,"car");
959fbac5 577 v.SetErrors(ea,"car");
d88f97cc 578
579 return v;
580}
581///////////////////////////////////////////////////////////////////////////
582Ali3Vector Ali3Vector::operator-(Ali3Vector& q)
583{
584// Subtract vector q from the current vector
959fbac5 585// Error propagation is performed automatically
586 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 587
588 GetVector(a,"car");
959fbac5 589 GetErrors(ea,"car");
d88f97cc 590 q.GetVector(b,"car");
959fbac5 591 q.GetErrors(eb,"car");
d88f97cc 592
593 for (Int_t i=0; i<3; i++)
594 {
595 a[i]-=b[i];
959fbac5 596 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 597 }
598
599 Ali3Vector v;
600 v.SetVector(a,"car");
959fbac5 601 v.SetErrors(ea,"car");
d88f97cc 602
603 return v;
604}
605///////////////////////////////////////////////////////////////////////////
606Ali3Vector Ali3Vector::operator*(Double_t s)
607{
959fbac5 608// Multiply the current vector with a scalar s.
609// Error propagation is performed automatically.
610 Double_t a[3],ea[3];
d88f97cc 611
612 GetVector(a,"car");
959fbac5 613 GetErrors(ea,"car");
d88f97cc 614
615 for (Int_t i=0; i<3; i++)
616 {
617 a[i]*=s;
959fbac5 618 ea[i]*=s;
d88f97cc 619 }
620
621 Ali3Vector v;
622 v.SetVector(a,"car");
959fbac5 623 v.SetErrors(ea,"car");
d88f97cc 624
625 return v;
626}
627///////////////////////////////////////////////////////////////////////////
628Ali3Vector Ali3Vector::operator/(Double_t s)
629{
630// Divide the current vector by a scalar s
959fbac5 631// Error propagation is performed automatically
d88f97cc 632
633 if (fabs(s)<1.e-20) // Protect against division by 0
634 {
635 cout << " *Ali3Vector::/* Division by 0 detected. No action taken." << endl;
636 return *this;
637 }
638 else
639 {
959fbac5 640 Double_t a[3],ea[3];
d88f97cc 641
642 GetVector(a,"car");
959fbac5 643 GetErrors(ea,"car");
d88f97cc 644
645 for (Int_t i=0; i<3; i++)
646 {
647 a[i]/=s;
959fbac5 648 ea[i]/=s;
d88f97cc 649 }
650
651 Ali3Vector v;
652 v.SetVector(a,"car");
959fbac5 653 v.SetErrors(ea,"car");
d88f97cc 654
655 return v;
656 }
657}
658///////////////////////////////////////////////////////////////////////////
659Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
660{
661// Add vector q to the current vector
959fbac5 662// Error propagation is performed automatically
663 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 664
665 GetVector(a,"car");
959fbac5 666 GetErrors(ea,"car");
d88f97cc 667 q.GetVector(b,"car");
959fbac5 668 q.GetErrors(eb,"car");
d88f97cc 669
670 for (Int_t i=0; i<3; i++)
671 {
672 a[i]+=b[i];
959fbac5 673 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 674 }
675
676 SetVector(a,"car");
959fbac5 677 SetErrors(ea,"car");
d88f97cc 678
679 return *this;
680}
681///////////////////////////////////////////////////////////////////////////
682Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
683{
684// Subtract vector q from the current vector
959fbac5 685// Error propagation is performed automatically
686 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 687
688 GetVector(a,"car");
959fbac5 689 GetErrors(ea,"car");
d88f97cc 690 q.GetVector(b,"car");
959fbac5 691 q.GetErrors(eb,"car");
d88f97cc 692
693 for (Int_t i=0; i<3; i++)
694 {
695 a[i]-=b[i];
959fbac5 696 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 697 }
698
699 SetVector(a,"car");
959fbac5 700 SetErrors(ea,"car");
d88f97cc 701
702 return *this;
703}
704///////////////////////////////////////////////////////////////////////////
705Ali3Vector& Ali3Vector::operator*=(Double_t s)
706{
707// Multiply the current vector with a scalar s
959fbac5 708// Error propagation is performed automatically
709 Double_t a[3],ea[3];
d88f97cc 710
711 GetVector(a,"car");
959fbac5 712 GetErrors(ea,"car");
d88f97cc 713
714 for (Int_t i=0; i<3; i++)
715 {
716 a[i]*=s;
959fbac5 717 ea[i]*=s;
d88f97cc 718 }
719
720 SetVector(a,"car");
959fbac5 721 SetErrors(ea,"car");
d88f97cc 722
723 return *this;
724}
725///////////////////////////////////////////////////////////////////////////
726Ali3Vector& Ali3Vector::operator/=(Double_t s)
727{
728// Divide the current vector by a scalar s
959fbac5 729// Error propagation is performed automatically
d88f97cc 730
731 if (fabs(s)<1.e-20) // Protect against division by 0
732 {
733 cout << " *Ali3Vector::/=* Division by 0 detected. No action taken." << endl;
734 return *this;
735 }
736 else
737 {
959fbac5 738 Double_t a[3],ea[3];
d88f97cc 739
740 GetVector(a,"car");
959fbac5 741 GetErrors(ea,"car");
d88f97cc 742
743 for (Int_t i=0; i<3; i++)
744 {
745 a[i]/=s;
959fbac5 746 ea[i]/=s;
d88f97cc 747 }
748
749 SetVector(a,"car");
959fbac5 750 SetErrors(ea,"car");
d88f97cc 751
752 return *this;
753 }
754}
755///////////////////////////////////////////////////////////////////////////
d071d629 756Ali3Vector Ali3Vector::GetVecTrans()
757{
758// Provide the transverse vector w.r.t. z-axis.
759// Error propagation is performed automatically
760 Double_t pi=acos(-1.);
761 Double_t a[3],ea[3];
762
763 GetVector(a,"sph");
764 GetErrors(ea,"sph");
765
766 Double_t vt,dvt2;
767 vt=a[0]*sin(a[1]);
768 dvt2=pow((sin(a[1])*ea[0]),2)+pow((a[0]*cos(a[1])*ea[1]),2);
769
770 a[0]=fabs(vt);
771 a[1]=pi/2.;
772
773 ea[0]=sqrt(dvt2);
774 ea[1]=0;
775
776 Ali3Vector v;
777 v.SetVector(a,"sph");
778 v.SetErrors(ea,"sph");
779
780 return v;
781}
782///////////////////////////////////////////////////////////////////////////
783Ali3Vector Ali3Vector::GetVecLong()
784{
785// Provide the longitudinal vector w.r.t. z-axis.
786// Error propagation is performed automatically
787 Double_t pi=acos(-1.);
788 Double_t a[3],ea[3];
789
790 GetVector(a,"sph");
791 GetErrors(ea,"sph");
792
793 Double_t vl,dvl2;
794 vl=a[0]*cos(a[1]);
795 dvl2=pow((cos(a[1])*ea[0]),2)+pow((a[0]*sin(a[1])*ea[1]),2);
796
797 a[0]=fabs(vl);
798 a[1]=0;
799 if (vl<0) a[1]=pi;
800 a[2]=0;
801
802 ea[0]=sqrt(dvl2);
803 ea[1]=0;
804 ea[2]=0;
805
806 Ali3Vector v;
807 v.SetVector(a,"sph");
808 v.SetErrors(ea,"sph");
809
810 return v;
811}
812///////////////////////////////////////////////////////////////////////////