]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RALICE/Ali3Vector.cxx
Using AliLog (F.Carminati)
[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///////////////////////////////////////////////////////////////////////////
1fbffa23 108void Ali3Vector::Load(Ali3Vector& q)
109{
110// Load all attributes of the input Ali3Vector into this Ali3Vector object.
111 Double_t temp=q.GetResultError();
112 Double_t a[3];
113 q.GetVector(a,"sph");
114 SetVector(a,"sph");
115 q.GetErrors(a,"car");
116 SetErrors(a,"car");
117 fDresult=temp;
118}
119///////////////////////////////////////////////////////////////////////////
c72198f1 120void Ali3Vector::SetZero()
121{
122// (Re)set all attributes to zero.
d88f97cc 123 fV=0;
124 fTheta=0;
125 fPhi=0;
959fbac5 126 fDx=0;
127 fDy=0;
128 fDz=0;
129 fDresult=0;
d88f97cc 130}
131///////////////////////////////////////////////////////////////////////////
d88f97cc 132void Ali3Vector::SetVector(Double_t* v,TString f)
133{
134// Store vector according to reference frame f
959fbac5 135// All errors will be reset to 0
136 fDx=0;
137 fDy=0;
138 fDz=0;
139 fDresult=0;
140
d88f97cc 141 Double_t pi=acos(-1.);
959fbac5 142
d88f97cc 143 Int_t frame=0;
144 if (f == "car") frame=1;
145 if (f == "sph") frame=2;
146 if (f == "cyl") frame=3;
147
148 Double_t x,y,z,rho,phi;
149
150 switch (frame)
151 {
152 case 1: // Cartesian coordinates
153 x=v[0];
154 y=v[1];
155 z=v[2];
156 fV=sqrt(x*x+y*y+z*z);
157 fTheta=0;
158 if (fV && fabs(z/fV)<=1.)
159 {
160 fTheta=acos(z/fV);
161 }
162 else
163 {
164 if (z<0.) fTheta=pi;
165 }
166 if (fTheta<0.) fTheta+=2.*pi;
167 fPhi=0;
168 if (x || y) fPhi=atan2(y,x);
169 if (fPhi<0.) fPhi+=2.*pi;
170 break;
171
172 case 2: // Spherical coordinates
173 fV=v[0];
174 fTheta=v[1];
175 fPhi=v[2];
176 break;
177
178 case 3: // Cylindrical coordinates
179 rho=v[0];
180 phi=v[1];
181 z=v[2];
182 fV=sqrt(rho*rho+z*z);
183 fPhi=phi;
184 if (fPhi<0.) fPhi+=2.*pi;
185 fTheta=0;
186 if (fV && fabs(z/fV)<=1.)
187 {
188 fTheta=acos(z/fV);
189 }
190 else
191 {
192 if (z<0.) fTheta=pi;
193 }
194 if (fTheta<0.) fTheta+=2.*pi;
195 break;
196
197 default: // Unsupported reference frame
c72198f1 198 cout << "*Ali3Vector::SetVector* Unsupported frame : " << f.Data() << endl
d88f97cc 199 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
200 fV=0;
201 fTheta=0;
202 fPhi=0;
203 break;
204 }
205}
206///////////////////////////////////////////////////////////////////////////
261c0caf 207void Ali3Vector::GetVector(Double_t* v,TString f) const
d88f97cc 208{
209// Provide vector according to reference frame f
210 Int_t frame=0;
211 if (f == "car") frame=1;
212 if (f == "sph") frame=2;
213 if (f == "cyl") frame=3;
214
215 switch (frame)
216 {
217 case 1: // Cartesian coordinates
218 v[0]=fV*sin(fTheta)*cos(fPhi);
219 v[1]=fV*sin(fTheta)*sin(fPhi);
220 v[2]=fV*cos(fTheta);
221 break;
222
223 case 2: // Spherical coordinates
224 v[0]=fV;
225 v[1]=fTheta;
226 v[2]=fPhi;
227 break;
228
229 case 3: // Cylindrical coordinates
230 v[0]=fV*sin(fTheta);
231 v[1]=fPhi;
232 v[2]=fV*cos(fTheta);
233 break;
234
235 default: // Unsupported reference frame
c72198f1 236 cout << "*Ali3Vector::GetVector* Unsupported frame : " << f.Data() << endl
d88f97cc 237 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
238 for (Int_t i=0; i<3; i++)
239 {
240 v[i]=0;
241 }
242 break;
243 }
244}
245///////////////////////////////////////////////////////////////////////////
246void Ali3Vector::SetVector(Float_t* v,TString f)
247{
248// Store vector according to reference frame f
959fbac5 249// All errors will be reset to 0
d88f97cc 250 Double_t vec[3];
251 for (Int_t i=0; i<3; i++)
252 {
253 vec[i]=v[i];
254 }
255 SetVector(vec,f);
256}
257///////////////////////////////////////////////////////////////////////////
261c0caf 258void Ali3Vector::GetVector(Float_t* v,TString f) const
d88f97cc 259{
260// Provide vector according to reference frame f
261 Double_t vec[3];
262 GetVector(vec,f);
263 for (Int_t i=0; i<3; i++)
264 {
265 v[i]=vec[i];
266 }
267}
268///////////////////////////////////////////////////////////////////////////
959fbac5 269void Ali3Vector::SetErrors(Double_t* e,TString f)
270{
271// Store errors according to reference frame f
272// The error on scalar results is reset to 0
273 fDresult=0;
274
275 Int_t frame=0;
276 if (f == "car") frame=1;
277 if (f == "sph") frame=2;
278 if (f == "cyl") frame=3;
279
280 Double_t dx2,dy2,dz2,rho;
281
282 switch (frame)
283 {
284 case 1: // Cartesian coordinates
285 fDx=fabs(e[0]);
286 fDy=fabs(e[1]);
287 fDz=fabs(e[2]);
288 break;
289
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);
296 fDx=sqrt(dx2);
297 fDy=sqrt(dy2);
298 fDz=sqrt(dz2);
299 break;
300
301 case 3: // Cylindrical coordinates
302 rho=fV*sin(fTheta);
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);
305 fDx=sqrt(dx2);
306 fDy=sqrt(dy2);
307 fDz=fabs(e[2]);
308 break;
309
310 default: // Unsupported reference frame
c72198f1 311 cout << "*Ali3Vector::SetErrors* Unsupported frame : " << f.Data() << endl
959fbac5 312 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
313 fDx=0;
314 fDy=0;
315 fDz=0;
316 break;
317 }
318}
319///////////////////////////////////////////////////////////////////////////
261c0caf 320void Ali3Vector::GetErrors(Double_t* e,TString f) const
959fbac5 321{
322// Provide errors according to reference frame f
323 Int_t frame=0;
324 if (f == "car") frame=1;
325 if (f == "sph") frame=2;
326 if (f == "cyl") frame=3;
327
d071d629 328 Double_t pi=acos(-1.);
329
959fbac5 330 Double_t dr2,dtheta2,dphi2,rho,drho2;
331 Double_t v[3];
d071d629 332 Double_t rxy2; // Shorthand for (x*x+y*y)
959fbac5 333
334 switch (frame)
335 {
336 case 1: // Cartesian coordinates
337 e[0]=fDx;
338 e[1]=fDy;
339 e[2]=fDz;
340 break;
341
342 case 2: // Spherical coordinates
343 GetVector(v,"car");
d071d629 344 rxy2=pow(v[0],2)+pow(v[1],2);
345 if (sqrt(rxy2)<(fV*1e-10)) rxy2=0;
959fbac5 346 if (fV)
347 {
348 dr2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)+pow((v[2]*fDz),2))/(fV*fV);
349 }
350 else
351 {
352 dr2=0;
353 }
d071d629 354 if (fV)
959fbac5 355 {
d071d629 356 dtheta2=rxy2*pow(fDz,2)/pow(fV,4);
357 if (v[2] && rxy2)
358 {
359 dtheta2+=rxy2*pow(v[2],2)*(pow((v[0]*fDx),2)+pow((v[1]*fDy),2)) /
81508922 360 pow(((pow(v[2],2)*rxy2)+pow(rxy2,2)),2);
d071d629 361 }
959fbac5 362 }
363 else
364 {
959fbac5 365 dtheta2=0;
366 }
d071d629 367 if (rxy2)
959fbac5 368 {
d071d629 369 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rxy2,2));
959fbac5 370 }
371 else
372 {
373 dphi2=0;
374 }
375 e[0]=sqrt(dr2);
376 e[1]=sqrt(dtheta2);
d071d629 377 if (e[1]>pi) e[1]=pi;
959fbac5 378 e[2]=sqrt(dphi2);
d071d629 379 if (e[2]>(2.*pi)) e[2]=2.*pi;
959fbac5 380 break;
381
382 case 3: // Cylindrical coordinates
383 GetVector(v,"car");
d071d629 384 rho=fabs(fV*sin(fTheta));
385 if (rho<(fV*1e-10)) rho=0;
959fbac5 386 if (rho)
387 {
388 drho2=(pow((v[0]*fDx),2)+pow((v[1]*fDy),2))/(rho*rho);
389 }
390 else
391 {
392 drho2=0;
393 }
d071d629 394 if (rho)
959fbac5 395 {
d071d629 396 dphi2=(pow((v[1]*fDx),2)+pow((v[0]*fDy),2))/(pow(rho,4));
959fbac5 397 }
398 else
399 {
400 dphi2=0;
401 }
402 e[0]=sqrt(drho2);
403 e[1]=sqrt(dphi2);
d071d629 404 if (e[1]>(2.*pi)) e[1]=2.*pi;
959fbac5 405 e[2]=fDz;
406 break;
407
408 default: // Unsupported reference frame
c72198f1 409 cout << "*Ali3Vector::GetErrors* Unsupported frame : " << f.Data() << endl
959fbac5 410 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
411 for (Int_t i=0; i<3; i++)
412 {
413 e[i]=0;
414 }
415 break;
416 }
417}
418///////////////////////////////////////////////////////////////////////////
419void Ali3Vector::SetErrors(Float_t* e,TString f)
420{
421// Store errors according to reference frame f
422// The error on scalar results is reset to 0
423 Double_t vec[3];
424 for (Int_t i=0; i<3; i++)
425 {
426 vec[i]=e[i];
427 }
428 SetErrors(vec,f);
429}
430///////////////////////////////////////////////////////////////////////////
261c0caf 431void Ali3Vector::GetErrors(Float_t* e,TString f) const
959fbac5 432{
433// Provide errors according to reference frame f
434 Double_t vec[3];
435 GetErrors(vec,f);
436 for (Int_t i=0; i<3; i++)
437 {
438 e[i]=vec[i];
439 }
440}
441///////////////////////////////////////////////////////////////////////////
261c0caf 442void Ali3Vector::Data(TString f) const
d88f97cc 443{
444// Print vector components according to reference frame f
445 if (f=="car" || f=="sph" || f=="cyl")
446 {
959fbac5 447 Double_t vec[3],err[3];
d88f97cc 448 GetVector(vec,f);
959fbac5 449 GetErrors(err,f);
c72198f1 450 cout << " Vector in " << f.Data() << " coordinates : "
d88f97cc 451 << vec[0] << " " << vec[1] << " " << vec[2] << endl;
c72198f1 452 cout << " Err. in " << f.Data() << " coordinates : "
959fbac5 453 << err[0] << " " << err[1] << " " << err[2] << endl;
d88f97cc 454 }
455 else
456 {
c72198f1 457 cout << " *Ali3Vector::Data* Unsupported frame : " << f.Data() << endl
d88f97cc 458 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
459 }
460}
461///////////////////////////////////////////////////////////////////////////
462Double_t Ali3Vector::GetNorm()
463{
959fbac5 464// Provide the norm of the current vector
465// The error on the scalar result (norm) is updated accordingly
466 Double_t e[3];
467 GetErrors(e,"sph");
468 fDresult=e[0];
d88f97cc 469 return fV;
470}
471///////////////////////////////////////////////////////////////////////////
959fbac5 472Double_t Ali3Vector::GetPseudoRapidity()
473{
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
8adaf597 477 Double_t pi=acos(-1.);
959fbac5 478 Double_t v[3];
479 GetVector(v,"sph");
480 Double_t thetahalf=v[1]/2.;
8adaf597 481 Double_t arg=0;
482 if (v[1]<pi) arg=tan(thetahalf);
483 Double_t eta=9999;
959fbac5 484 if (arg>0) eta=-log(arg);
485 Double_t e[3];
486 GetErrors(e,"sph");
487 Double_t prod=cos(thetahalf)*sin(thetahalf);
488 fDresult=0;
489 if (prod) fDresult=fabs(e[1]/2.*prod);
490 return eta;
491}
492///////////////////////////////////////////////////////////////////////////
d88f97cc 493Double_t Ali3Vector::Dot(Ali3Vector& q)
494{
495// Provide the dot product of the current vector with vector q
959fbac5 496// The error on the scalar result (dotproduct) is updated accordingly
497
d88f97cc 498 Double_t dotpro=0;
499
959fbac5 500 if ((this) == &q) // Check for special case v.Dot(v)
d88f97cc 501 {
959fbac5 502 Double_t norm=GetNorm();
503 Double_t dnorm=GetResultError();
504 dotpro=pow(norm,2);
505 fDresult=2.*norm*dnorm;
d88f97cc 506 }
959fbac5 507 else
508 {
509 Double_t a[3],b[3];
510 Double_t ea[3],eb[3];
511 Double_t d2=0;
512
513 GetVector(a,"car");
514 GetErrors(ea,"car");
515 q.GetVector(b,"car");
516 q.GetErrors(eb,"car");
517 for (Int_t i=0; i<3; i++)
518 {
519 dotpro+=a[i]*b[i];
520 d2+=pow(b[i]*ea[i],2)+pow(a[i]*eb[i],2);
521 }
522 fDresult=sqrt(d2);
523 }
524
d88f97cc 525 return dotpro;
526}
527///////////////////////////////////////////////////////////////////////////
261c0caf 528Double_t Ali3Vector::GetResultError() const
959fbac5 529{
530// Provide the error on the result of an operation yielding a scalar
531// E.g. GetNorm() or Dot()
532 return fDresult;
533}
534///////////////////////////////////////////////////////////////////////////
261c0caf 535Ali3Vector Ali3Vector::Cross(Ali3Vector& q) const
d88f97cc 536{
537// Provide the cross product of the current vector with vector q
959fbac5 538// Error propagation is performed automatically
d88f97cc 539 Double_t a[3],b[3],c[3];
959fbac5 540 Double_t ea[3],eb[3],ec[3],d2;
d88f97cc 541
542 GetVector(a,"car");
959fbac5 543 GetErrors(ea,"car");
d88f97cc 544 q.GetVector(b,"car");
959fbac5 545 q.GetErrors(eb,"car");
d88f97cc 546
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];
550
959fbac5 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);
553 ec[0]=sqrt(d2);
554
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);
557 ec[1]=sqrt(d2);
558
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);
561 ec[2]=sqrt(d2);
562
d88f97cc 563 Ali3Vector v;
564 v.SetVector(c,"car");
959fbac5 565 v.SetErrors(ec,"car");
d88f97cc 566
567 return v;
568}
569///////////////////////////////////////////////////////////////////////////
261c0caf 570Ali3Vector Ali3Vector::operator+(Ali3Vector& q) const
d88f97cc 571{
572// Add vector q to the current vector
959fbac5 573// Error propagation is performed automatically
574 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 575
576 GetVector(a,"car");
959fbac5 577 GetErrors(ea,"car");
d88f97cc 578 q.GetVector(b,"car");
959fbac5 579 q.GetErrors(eb,"car");
d88f97cc 580
581 for (Int_t i=0; i<3; i++)
582 {
583 a[i]+=b[i];
959fbac5 584 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 585 }
586
587 Ali3Vector v;
588 v.SetVector(a,"car");
959fbac5 589 v.SetErrors(ea,"car");
d88f97cc 590
591 return v;
592}
593///////////////////////////////////////////////////////////////////////////
261c0caf 594Ali3Vector Ali3Vector::operator-(Ali3Vector& q) const
d88f97cc 595{
596// Subtract vector q from the current vector
959fbac5 597// Error propagation is performed automatically
598 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 599
600 GetVector(a,"car");
959fbac5 601 GetErrors(ea,"car");
d88f97cc 602 q.GetVector(b,"car");
959fbac5 603 q.GetErrors(eb,"car");
d88f97cc 604
605 for (Int_t i=0; i<3; i++)
606 {
607 a[i]-=b[i];
959fbac5 608 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 609 }
610
611 Ali3Vector v;
612 v.SetVector(a,"car");
959fbac5 613 v.SetErrors(ea,"car");
d88f97cc 614
615 return v;
616}
617///////////////////////////////////////////////////////////////////////////
261c0caf 618Ali3Vector Ali3Vector::operator*(Double_t s) const
d88f97cc 619{
959fbac5 620// Multiply the current vector with a scalar s.
621// Error propagation is performed automatically.
622 Double_t a[3],ea[3];
d88f97cc 623
624 GetVector(a,"car");
959fbac5 625 GetErrors(ea,"car");
d88f97cc 626
627 for (Int_t i=0; i<3; i++)
628 {
629 a[i]*=s;
959fbac5 630 ea[i]*=s;
d88f97cc 631 }
632
633 Ali3Vector v;
634 v.SetVector(a,"car");
959fbac5 635 v.SetErrors(ea,"car");
d88f97cc 636
637 return v;
638}
639///////////////////////////////////////////////////////////////////////////
261c0caf 640Ali3Vector Ali3Vector::operator/(Double_t s) const
d88f97cc 641{
642// Divide the current vector by a scalar s
959fbac5 643// Error propagation is performed automatically
d88f97cc 644
645 if (fabs(s)<1.e-20) // Protect against division by 0
646 {
647 cout << " *Ali3Vector::/* Division by 0 detected. No action taken." << endl;
648 return *this;
649 }
650 else
651 {
959fbac5 652 Double_t a[3],ea[3];
d88f97cc 653
654 GetVector(a,"car");
959fbac5 655 GetErrors(ea,"car");
d88f97cc 656
657 for (Int_t i=0; i<3; i++)
658 {
659 a[i]/=s;
959fbac5 660 ea[i]/=s;
d88f97cc 661 }
662
663 Ali3Vector v;
664 v.SetVector(a,"car");
959fbac5 665 v.SetErrors(ea,"car");
d88f97cc 666
667 return v;
668 }
669}
670///////////////////////////////////////////////////////////////////////////
671Ali3Vector& Ali3Vector::operator+=(Ali3Vector& q)
672{
673// Add vector q to the current vector
959fbac5 674// Error propagation is performed automatically
675 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 676
677 GetVector(a,"car");
959fbac5 678 GetErrors(ea,"car");
d88f97cc 679 q.GetVector(b,"car");
959fbac5 680 q.GetErrors(eb,"car");
d88f97cc 681
682 for (Int_t i=0; i<3; i++)
683 {
684 a[i]+=b[i];
959fbac5 685 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 686 }
687
688 SetVector(a,"car");
959fbac5 689 SetErrors(ea,"car");
d88f97cc 690
691 return *this;
692}
693///////////////////////////////////////////////////////////////////////////
694Ali3Vector& Ali3Vector::operator-=(Ali3Vector& q)
695{
696// Subtract vector q from the current vector
959fbac5 697// Error propagation is performed automatically
698 Double_t a[3],b[3],ea[3],eb[3];
d88f97cc 699
700 GetVector(a,"car");
959fbac5 701 GetErrors(ea,"car");
d88f97cc 702 q.GetVector(b,"car");
959fbac5 703 q.GetErrors(eb,"car");
d88f97cc 704
705 for (Int_t i=0; i<3; i++)
706 {
707 a[i]-=b[i];
959fbac5 708 ea[i]=sqrt(pow(ea[i],2)+pow(eb[i],2));
d88f97cc 709 }
710
711 SetVector(a,"car");
959fbac5 712 SetErrors(ea,"car");
d88f97cc 713
714 return *this;
715}
716///////////////////////////////////////////////////////////////////////////
717Ali3Vector& Ali3Vector::operator*=(Double_t s)
718{
719// Multiply the current vector with a scalar s
959fbac5 720// Error propagation is performed automatically
721 Double_t a[3],ea[3];
d88f97cc 722
723 GetVector(a,"car");
959fbac5 724 GetErrors(ea,"car");
d88f97cc 725
726 for (Int_t i=0; i<3; i++)
727 {
728 a[i]*=s;
959fbac5 729 ea[i]*=s;
d88f97cc 730 }
731
732 SetVector(a,"car");
959fbac5 733 SetErrors(ea,"car");
d88f97cc 734
735 return *this;
736}
737///////////////////////////////////////////////////////////////////////////
738Ali3Vector& Ali3Vector::operator/=(Double_t s)
739{
740// Divide the current vector by a scalar s
959fbac5 741// Error propagation is performed automatically
d88f97cc 742
743 if (fabs(s)<1.e-20) // Protect against division by 0
744 {
745 cout << " *Ali3Vector::/=* Division by 0 detected. No action taken." << endl;
746 return *this;
747 }
748 else
749 {
959fbac5 750 Double_t a[3],ea[3];
d88f97cc 751
752 GetVector(a,"car");
959fbac5 753 GetErrors(ea,"car");
d88f97cc 754
755 for (Int_t i=0; i<3; i++)
756 {
757 a[i]/=s;
959fbac5 758 ea[i]/=s;
d88f97cc 759 }
760
761 SetVector(a,"car");
959fbac5 762 SetErrors(ea,"car");
d88f97cc 763
764 return *this;
765 }
766}
767///////////////////////////////////////////////////////////////////////////
261c0caf 768Ali3Vector Ali3Vector::GetVecTrans() const
d071d629 769{
770// Provide the transverse vector w.r.t. z-axis.
771// Error propagation is performed automatically
772 Double_t pi=acos(-1.);
773 Double_t a[3],ea[3];
774
775 GetVector(a,"sph");
776 GetErrors(ea,"sph");
777
778 Double_t vt,dvt2;
779 vt=a[0]*sin(a[1]);
780 dvt2=pow((sin(a[1])*ea[0]),2)+pow((a[0]*cos(a[1])*ea[1]),2);
781
782 a[0]=fabs(vt);
783 a[1]=pi/2.;
784
785 ea[0]=sqrt(dvt2);
786 ea[1]=0;
787
788 Ali3Vector v;
789 v.SetVector(a,"sph");
790 v.SetErrors(ea,"sph");
791
792 return v;
793}
794///////////////////////////////////////////////////////////////////////////
261c0caf 795Ali3Vector Ali3Vector::GetVecLong() const
d071d629 796{
797// Provide the longitudinal vector w.r.t. z-axis.
798// Error propagation is performed automatically
799 Double_t pi=acos(-1.);
800 Double_t a[3],ea[3];
801
802 GetVector(a,"sph");
803 GetErrors(ea,"sph");
804
805 Double_t vl,dvl2;
806 vl=a[0]*cos(a[1]);
807 dvl2=pow((cos(a[1])*ea[0]),2)+pow((a[0]*sin(a[1])*ea[1]),2);
808
809 a[0]=fabs(vl);
810 a[1]=0;
811 if (vl<0) a[1]=pi;
812 a[2]=0;
813
814 ea[0]=sqrt(dvl2);
815 ea[1]=0;
816 ea[2]=0;
817
818 Ali3Vector v;
819 v.SetVector(a,"sph");
820 v.SetErrors(ea,"sph");
821
822 return v;
823}
824///////////////////////////////////////////////////////////////////////////
261c0caf 825Ali3Vector Ali3Vector::GetPrimed(TRotMatrix* m) const
826{
827// Provide vector components (and errors) in a rotated frame.
828// The orientation of the rotated frame is described by the TRotMatrix
829// input argument.
830 Ali3Vector v=*this;
831 if (!m) return v;
832
833 Double_t* mat=m->GetMatrix();
834
835 Double_t a[3],aprim[3];
836
837 GetVector(a,"car");
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");
842
843 GetErrors(a,"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");
848
849 return v;
850}
851///////////////////////////////////////////////////////////////////////////
c5555bc0 852Ali3Vector Ali3Vector::GetUnprimed(TRotMatrix* m) const
853{
854// Provide original vector components (and errors) from the rotated ones.
855// The orientation of the rotated frame is described by the TRotMatrix
856// input argument.
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.
860 Ali3Vector v=*this;
861 if (!m) return v;
862
863 Double_t* mat=m->GetMatrix();
864
865 Double_t a[3],aprim[3];
866
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");
872
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");
878
879 return v;
880}
881///////////////////////////////////////////////////////////////////////////
a7dc0627 882Double_t Ali3Vector::GetX(Int_t i,TString f)
883{
884// Provide i-th vector component according to reference frame f.
885// The vector components are addressed via the generic x1,x2,x3 notation.
886// So, i=1 denotes the first vector component.
887// The error on the selected component can be obtained via the
888// usual GetResultError() facility.
889
890 fDresult=0;
891
892 if (i<1 || i>3) return 0;
893
894 Double_t vec[3];
895 Double_t err[3];
896 GetVector(vec,f);
897 GetErrors(err,f);
898
899 fDresult=err[i-1];
900
901 return vec[i-1];
902}
903///////////////////////////////////////////////////////////////////////////