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