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