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