334e554299c64ded821ee70ff75d090ccda751c6
[u/mrichter/AliRoot.git] / RALICE / Ali4Vector.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 Ali4Vector
20 // Handling of Lorentz 4-vectors in various reference frames.
21 //
22 // This class is meant to serve as a base class for ALICE objects
23 // that have Lorentz 4-vector characteristics.
24 // Error propagation is performed automatically.
25 //
26 // All 4-vectors are treated in the contravariant form and the convention
27 // for the metric and the 4-vector components is according to the one
28 // used in the book "Classical Electrodynamics" by J.D. Jackson.
29 //
30 // A 4-vector is said to have a scalar part and a 3-vector part,
31 // which is indicated by the notation
32 //
33 //    x^i = (x^0,x^1,x^2,x^3)
34 // 
35 // The scalar part   = x^0
36 // The 3-vector part = (x^1,x^2,x^3)
37 //
38 // In view of accuracy and the fact that e.g. particle identity (mass)
39 // is preserved in many physics processes, the Lorentz invariant 
40 // (x^i*x_i) is internally saved together with the scalar part.
41 //
42 // This allows the following two modes of functionality :
43 //
44 // Scalar mode    : The scalar part and the 3-vector part are considered
45 //                  as basic quantities and the invariant with its error
46 //                  is derived from these.  
47 // Invariant mode : The invariant and the 3-vector part are considered
48 //                  as basic quantities and the scalar with its error
49 //                  is derived from these.
50 //
51 // The philosophy followed here is the following :
52 // ===============================================
53 //
54 // 1) Invokation of SetVector() sets the scalar and 3-vector parts 
55 //    and the invariant is calculated from these.
56 //    Automatically the scalar mode is selected and invokation of
57 //    SetErrors() will calculate the error on the invariant.
58 //
59 // 2) In case the scalar part is modified via SetScalar(), scalar mode is
60 //    automatically selected and the Lorentz invariant (x^i*x_i) and its
61 //    error are updated accordingly.
62 //    The 3-vector part is NOT modified.
63 //    This situation arises when one e.g. precisely determines the time
64 //    or energy (x^0).     
65 //
66 // 3) In case the Lorentz invariant (x^i*x_i) is modified via SetInvariant(),
67 //    invariant mode is selected automatically and the scalar part and its
68 //    error are updated accordingly.
69 //    The 3-vector part is NOT modified.
70 //    This situation arises when one e.g. precisely determines the mass.     
71 //
72 // 4) In case the vector part is modified via Set3Vector(), then the 
73 //    current mode determines whether the scalar or the invariant is updated. 
74 //    Scalar mode    : The Lorentz invariant (x^i*x_i) and its error are updated;
75 //                     the scalar part and its error are NOT modified. 
76 //                     This situation arises when one e.g. improves the 3-position
77 //                     vector for a particle with a very precise timing.     
78 //    Invariant mode : The scalar part and its error are updated;
79 //                     the Lorentz invariant (x^i*x_i) and its error are NOT modified.
80 //                     This situation arises when one e.g. improves the 3-momentum
81 //                     vector for a particle with known mass.     
82 //
83 // The dotproduct is defined such that p.Dot(p) yields the Lorentz invariant
84 // scalar of the 4-vector p (i.e. m**2 in case p is a 4-momentum).   
85 //
86 // Note :
87 // ------
88 // Vectors (v), Errors (e) and reference frames (f) are specified via
89 // SetVector(Float_t* v,TString f)
90 // SetErrors(Float_t* e,TString f)
91 // under the following conventions :
92 //
93 // f="car" ==> 3-vector part of v in Cartesian coordinates   (x,y,z)
94 // f="sph" ==> 3-vector part of v in Spherical coordinates   (r,theta,phi)
95 // f="cyl" ==> 3-vector part of v in Cylindrical coordinates (rho,phi,z)
96 //
97 // All angles are in radians.
98 //
99 // Example :
100 // ---------
101 //
102 // Ali4Vector a;
103 //
104 // Float_t v[4]={25,-1,3,7};
105 // a.SetVector(v,"car");
106 //
107 // Float_t vec[4];
108 // a.GetVector(vec,"sph");
109 //
110 // Ali4Vector b;
111 // Float_t v2[4]={33,6,-18,2};
112 // b.SetVector(v2,"car");
113 //
114 // Float_t dotpro=a.Dot(b);
115 //
116 // Float_t x0=16;
117 // Ali3Vector x;
118 // Float_t vec2[3]={1,2,3};
119 // x.SetVector(vec2,"car");
120 //
121 // Ali4Vector c;
122 // c.SetVector(x0,x);
123 // c.GetVector(vec,"car");
124 // c.Data("cyl");
125 // c=a+b;
126 // c=a-b;
127 // c=a*5;
128 //
129 //--- Author: Nick van Eijndhoven 01-apr-1999 UU-SAP Utrecht
130 //- Modified: NvE $Date$ UU-SAP Utrecht
131 ///////////////////////////////////////////////////////////////////////////
132
133 #include "Ali4Vector.h"
134 #include "Riostream.h"
135  
136 ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
137  
138 Ali4Vector::Ali4Vector()
139 {
140 // Creation of a contravariant 4-vector and initialisation of parameters.
141 // All values are initialised to 0.
142 // Scalar mode is initially selected.
143  SetZero();
144  fScalar=1;
145 }
146 ///////////////////////////////////////////////////////////////////////////
147 Ali4Vector::~Ali4Vector()
148 {
149 // Destructor to delete dynamically allocated memory
150 }
151 ///////////////////////////////////////////////////////////////////////////
152 Ali4Vector::Ali4Vector(const Ali4Vector& v)
153 {
154 // Copy constructor
155  fScalar=v.fScalar;
156  fV2=v.fV2;
157  fDv2=v.fDv2;
158  fV0=v.fV0;
159  fDv0=v.fDv0;
160  fDresult=v.fDresult;
161  fV=v.fV;
162 }
163 ///////////////////////////////////////////////////////////////////////////
164 void Ali4Vector::Load(Ali4Vector& q)
165 {
166 // Load all attributes of the input Ali4Vector into this Ali4Vector object.
167  Int_t temp1=q.GetScalarFlag();
168  Double_t temp2=q.GetResultError();
169  Double_t a[4];
170  q.GetVector(a,"sph");
171  SetVector(a,"sph");
172  q.GetErrors(a,"car");
173  SetErrors(a,"car");
174  fScalar=temp1;
175  fDresult=temp2;
176 }
177 ///////////////////////////////////////////////////////////////////////////
178 void Ali4Vector::SetZero()
179 {
180 // (Re)set all attributes to zero.
181 // Note : The (de)selection of the scalar mode is not modified.
182  fV2=0;
183  fDv2=0;
184  fV0=0;
185  fDv0=0;
186  fDresult=0;
187  fV.SetZero();
188 }
189 ///////////////////////////////////////////////////////////////////////////
190 void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v)
191 {
192 // Store contravariant vector.
193 // The error on the scalar part is initialised to 0.
194 // The errors on the vector part are taken from the input Ali3Vector.
195 // Scalar mode is automatically selected. 
196 // The error on scalar result operations is reset to 0.
197  fScalar=1;
198  fV0=v0;
199  fV=v;
200  fV2=pow(v0,2)-fV.Dot(fV);
201  SetScalarError(0);
202 }
203 ///////////////////////////////////////////////////////////////////////////
204 void Ali4Vector::SetVector(Double_t* v,TString f)
205 {
206 // Store vector according to reference frame f.
207 // All errors are initialised to 0.
208 // Scalar mode is automatically selected. 
209 // The error on scalar result operations is reset to 0.
210  fScalar=1;
211  Double_t a[3];
212  for (Int_t i=0; i<3; i++)
213  {
214   a[i]=v[i+1];
215  }
216  fV0=v[0];
217  fV.SetVector(a,f);
218  fV2=pow(fV0,2)-fV.Dot(fV);
219  fDv2=0;
220  fDv0=0;
221  fDresult=0;
222 }
223 ///////////////////////////////////////////////////////////////////////////
224 void Ali4Vector::GetVector(Double_t* v,TString f)
225 {
226 // Provide 4-vector components according to reference frame f
227 // and according to the current mode.
228 // Scalar mode    : The scalar part is directly returned via v[0].
229 // Invariant mode : The scalar part is re-calculated via the value
230 //                  of the Lorentz invariant and then returned via v[0].
231  if (fScalar)
232  {
233   v[0]=fV0;
234  }
235  else
236  {
237   v[0]=sqrt(fV.Dot(fV)+fV2);
238  } 
239  Double_t a[3];
240  fV.GetVector(a,f);
241  for (Int_t i=0; i<3; i++)
242  {
243   v[i+1]=a[i];
244  }
245 }
246 ///////////////////////////////////////////////////////////////////////////
247 void Ali4Vector::SetVector(Float_t* v,TString f)
248 {
249 // Store vector according to reference frame f.
250 // All errors are initialised to 0.
251 // Scalar mode is automatically selected. 
252 // The error on scalar result operations is reset to 0.
253  Double_t vec[4];
254  for (Int_t i=0; i<4; i++)
255  {
256   vec[i]=v[i];
257  }
258  SetVector(vec,f);
259 }
260 ///////////////////////////////////////////////////////////////////////////
261 void Ali4Vector::GetVector(Float_t* v,TString f)
262 {
263 // Provide 4-vector components according to reference frame f
264 // and according to the current mode.
265 // Scalar mode    : The scalar part is directly returned via v[0].
266 // Invariant mode : The scalar part is re-calculated via the value
267 //                  of the Lorentz invariant and then returned via v[0].
268  Double_t vec[4];
269  GetVector(vec,f);
270  for (Int_t i=0; i<4; i++)
271  {
272   v[i]=vec[i];
273  }
274 }
275 ///////////////////////////////////////////////////////////////////////////
276 Double_t Ali4Vector::GetScalar()
277 {
278 // Provide the scalar part.
279 // The error on the scalar value is available via GetResultError()
280 // after invokation of GetScalar().
281  if (fScalar)
282  {
283   fDresult=fDv0;
284   return fV0;
285  }
286  else
287  {
288   Double_t dot=fV.Dot(fV);
289   Double_t ddot=fV.GetResultError();
290   Double_t v02=dot+fV2;
291   Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2));
292   Double_t v0=sqrt(fabs(v02));
293   Double_t dv0=0;
294   if (v0) dv0=dv02/(2.*v0);
295   fDresult=dv0;
296   return v0;
297  }
298 }
299 ///////////////////////////////////////////////////////////////////////////
300 Double_t Ali4Vector::GetResultError()
301 {
302 // Provide the error on the result of an operation yielding a scalar
303 // E.g. GetScalar(), GetInvariant() or Dot()
304  return fDresult;
305 }
306 ///////////////////////////////////////////////////////////////////////////
307 void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
308 {
309 // Modify the scalar part (v0) and its error (dv0).
310 // The default value for dv0 is 0.
311 // The vector part is not modified. 
312 // Scalar mode is automatically selected
313 // ==> Lorentz invariant and its error are updated. 
314 // The error on scalar result operations is reset to 0.
315  fScalar=1;
316  fV0=v0;
317  fV2=pow(v0,2)-fV.Dot(fV);
318  SetScalarError(dv0);
319 }
320 ///////////////////////////////////////////////////////////////////////////
321 void Ali4Vector::SetScalarError(Double_t dv0)
322 {
323 // Set the error on the scalar part.
324 // If in scalar mode, update error on the invariant accordingly.
325 // The error on scalar result operations is reset to 0.
326  fDv0=dv0;
327  if (fScalar)
328  {
329   Double_t norm=fV.GetNorm();
330   Double_t dnorm=fV.GetResultError();
331   fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2));
332  } 
333  fDresult=0;
334 }
335 ///////////////////////////////////////////////////////////////////////////
336 void Ali4Vector::Set3Vector(Ali3Vector& v)
337 {
338 // Set the 3-vector part, the errors are taken from the input Ali3Vector
339 // Scalar mode    : The scalar part and its error are not modified,
340 //                  the Lorentz invariantand its error are re-calculated.
341 // Invariant mode : The Lorentz invariant and its error are not modified,
342 //                  the scalar part and its error are re-calculated.
343 // The error on scalar result operations is reset to 0.
344  fV=v;
345  if (fScalar)
346  {
347   SetScalar(fV0,fDv0);
348  }
349  else
350  {
351   SetInvariant(fV2,fDv2);
352  }
353 }
354 ///////////////////////////////////////////////////////////////////////////
355 void Ali4Vector::Set3Vector(Double_t* v,TString f)
356 {
357 // Set the 3-vector part according to reference frame f
358 // The errors on the vector part are initialised to 0
359 // Scalar mode    : The scalar part and its error are not modified,
360 //                  the Lorentz invariantand its error are re-calculated.
361 // Invariant mode : The Lorentz invariant and its error are not modified,
362 //                  the scalar part and its error are re-calculated.
363 // The error on scalar result operations is reset to 0.
364  Double_t a[3];
365  for (Int_t i=0; i<3; i++)
366  {
367   a[i]=v[i];
368  }
369  fV.SetVector(a,f);
370
371  if (fScalar)
372  {
373   SetScalar(fV0,fDv0);
374  }
375  else
376  {
377   SetInvariant(fV2,fDv2);
378  }
379 }
380 ///////////////////////////////////////////////////////////////////////////
381 void Ali4Vector::Set3Vector(Float_t* v,TString f)
382 {
383 // Set the 3-vector part according to reference frame f
384 // The errors on the vector part are initialised to 0
385 // The Lorentz invariant is not modified
386 // The error on scalar result operations is reset to 0.
387  Double_t vec[3];
388  for (Int_t i=0; i<3; i++)
389  {
390   vec[i]=v[i];
391  }
392  Set3Vector(vec,f);
393 }
394 ///////////////////////////////////////////////////////////////////////////
395 void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
396 {
397 // Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2).
398 // The default value for the error dv2 is 0.
399 // The vector part is not modified.
400 // Invariant mode is automatically selected
401 // ==> the scalar part and its error are updated.
402 // The error on scalar result operations is reset to 0.
403 //
404  fScalar=0;
405  fV2=v2;
406  fDv2=dv2;
407  fV0=GetScalar();
408  fDv0=GetResultError();
409  fDresult=0;
410 }
411 ///////////////////////////////////////////////////////////////////////////
412 void Ali4Vector::SetInvariantError(Double_t dv2)
413 {
414 // Set the error on the Lorentz invariant.
415 // If in invariant mode, update error on the scalar part accordingly.
416 // The error on scalar result operations is reset to 0.
417  fDv2=dv2;
418  if (!fScalar)
419  {
420   fDv0=GetResultError();
421  }
422  fDresult=0; 
423 }
424 ///////////////////////////////////////////////////////////////////////////
425 Double_t Ali4Vector::GetInvariant()
426 {
427 // Provide the Lorentz invariant v^i*v_i.
428 // The error on the Lorentz invariant is available via GetResultError()
429 // after invokation of GetInvariant().
430  if (!fScalar)
431  {
432   fDresult=fDv2;
433   return fV2;
434  }
435  else
436  {
437   Double_t inv=Dot(*this);
438   return inv;
439  }
440 }
441 ///////////////////////////////////////////////////////////////////////////
442 Ali3Vector Ali4Vector::Get3Vector()
443 {
444 // Provide the 3-vector part
445  return fV;
446 }
447 ///////////////////////////////////////////////////////////////////////////
448 void Ali4Vector::SetErrors(Double_t* e,TString f)
449 {
450 // Store errors for vector v^i according to reference frame f
451 // If in scalar mode, update error on the invariant accordingly.
452 // The error on scalar result operations is reset to 0.
453  Double_t a[3];
454  for (Int_t i=0; i<3; i++)
455  {
456   a[i]=e[i+1];
457  }
458  SetScalarError(e[0]);
459  fV.SetErrors(a,f);
460 }
461 ///////////////////////////////////////////////////////////////////////////
462 void Ali4Vector::SetErrors(Float_t* e,TString f)
463 {
464 // Store errors for vector v^i according to reference frame f
465 // If in scalar mode, update error on the invariant accordingly.
466 // The error on scalar result operations is reset to 0.
467  Double_t a[4];
468  for (Int_t i=0; i<4; i++)
469  {
470   a[i]=e[i];
471  }
472  SetErrors(a,f);
473 }
474 ///////////////////////////////////////////////////////////////////////////
475 void Ali4Vector::GetErrors(Double_t* e,TString f)
476 {
477 // Provide errors for vector v^i according to reference frame f
478 // and according to the current mode.
479 // Scalar mode    : The error on the scalar part is directly returned via e[0].
480 // Invariant mode : The error on the scalar part is re-calculated via the error
481 //                  value on the Lorentz invariant and then returned via e[0].
482  Double_t a[3];
483  fV.GetErrors(a,f);
484
485  // Dummy invokation of GetScalar to obtain automatic proper error determination 
486  e[0]=GetScalar();
487
488  e[0]=GetResultError();
489
490  for (Int_t i=0; i<3; i++)
491  {
492   e[i+1]=a[i];
493  }
494 }
495 ///////////////////////////////////////////////////////////////////////////
496 void Ali4Vector::GetErrors(Float_t* e,TString f)
497 {
498 // Provide errors for vector v^i according to reference frame f
499 // and according to the current mode.
500 // Scalar mode    : The error on the scalar part is directly returned via e[0].
501 // Invariant mode : The error on the scalar part is re-calculated via the error
502 //                  value on the Lorentz invariant and then returned via e[0].
503  Double_t a[4];
504  GetErrors(a,f);
505  for (Int_t i=0; i<4; i++)
506  {
507   e[i]=a[i];
508  }
509 }
510 ///////////////////////////////////////////////////////////////////////////
511 void Ali4Vector::Data(TString f)
512 {
513 // Print contravariant vector components and errors according to
514 // reference frame f and according to the current mode.
515 // Scalar mode    : The scalar part and its error are directly returned.
516 // Invariant mode : The scalar part and its error are re-calculated via the
517 //                  value (and error) of the Lorentz invariant.
518
519  if (f=="car" || f=="sph" || f=="cyl")
520  {
521   Double_t vec[4],err[4];
522   GetVector(vec,f);
523   GetErrors(err,f);
524   Double_t inv=GetInvariant();
525   Double_t dinv=GetResultError();
526   cout << " Contravariant vector in " << f.Data() << " coordinates : "
527        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl; 
528   cout << " ------------- Errors in " << f.Data() << " coordinates : "
529        << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl; 
530   cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
531  }
532  else
533  {
534   cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl
535        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl; 
536  }
537 }
538 ///////////////////////////////////////////////////////////////////////////
539 Double_t Ali4Vector::Dot(Ali4Vector& q)
540 {
541 // Provide the dot product of the current vector with vector q
542  Double_t dotpro=0;
543  Double_t a0=GetScalar();
544  Double_t da0=GetResultError();
545  if ((this) == &q) // Check for special case v.Dot(v)
546  {
547   Double_t norm=fV.GetNorm();
548   Double_t dnorm=fV.GetResultError();
549   dotpro=pow(a0,2)-pow(norm,2);
550   fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2));
551  }
552  else
553  {
554   Double_t b0=q.GetScalar();
555   Double_t db0=q.GetResultError();
556   Ali3Vector b=q.Get3Vector();
557
558   Double_t dot=fV.Dot(b);
559   Double_t ddot=fV.GetResultError();
560
561   dotpro=a0*b0-dot;
562
563   fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
564  }
565
566  return dotpro;
567 }
568 ///////////////////////////////////////////////////////////////////////////
569 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
570 {
571 // Add 4-vector q to the current 4-vector
572 // Error propagation is performed automatically
573  Double_t a0=GetScalar();
574  Double_t da0=GetResultError();
575  Ali3Vector a=Get3Vector();
576  Double_t b0=q.GetScalar();
577  Double_t db0=q.GetResultError();
578  Ali3Vector b=q.Get3Vector();
579
580  Double_t c0=a0+b0;
581  Ali3Vector c=a+b;
582  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
583
584  Ali4Vector v;
585  v.SetVector(c0,c);
586  v.SetScalarError(dc0);  
587  return v;
588 }
589 ///////////////////////////////////////////////////////////////////////////
590 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
591 {
592 // Subtract 4-vector q from the current 4-vector
593 // Error propagation is performed automatically
594  Double_t a0=GetScalar();
595  Double_t da0=GetResultError();
596  Ali3Vector a=Get3Vector();
597  Double_t b0=q.GetScalar();
598  Double_t db0=q.GetResultError();
599  Ali3Vector b=q.Get3Vector();
600
601  Double_t c0=a0-b0;
602  Ali3Vector c=a-b;
603  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
604
605  Ali4Vector v;
606  v.SetVector(c0,c);
607  v.SetScalarError(dc0);  
608  return v;
609 }
610 ///////////////////////////////////////////////////////////////////////////
611 Ali4Vector Ali4Vector::operator*(Double_t s)
612 {
613 // Multiply the current 4-vector with a scalar s
614 // Error propagation is performed automatically
615  Double_t a0=GetScalar();
616  Double_t da0=GetResultError();
617  Ali3Vector a=Get3Vector();
618
619  a0*=s;
620  a*=s;
621  da0*=s;
622
623  Ali4Vector v;
624  v.SetVector(a0,a);
625  v.SetScalarError(da0);  
626   
627  return v;
628 }
629 ///////////////////////////////////////////////////////////////////////////
630 Ali4Vector Ali4Vector::operator/(Double_t s)
631 {
632 // Divide the current vector by a scalar s
633 // Error propagation is performed automatically
634
635  if (fabs(s)<1.e-20) // Protect against division by 0
636  {
637   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
638   return *this;
639  }
640  else
641  {
642   Double_t a0=GetScalar();
643   Double_t da0=GetResultError();
644   Ali3Vector a=Get3Vector();
645
646   a0/=s;
647   a/=s;
648   da0/=s;
649
650   Ali4Vector v;
651   v.SetVector(a0,a);
652   v.SetScalarError(da0);  
653   
654   return v;
655  }
656 }
657 ///////////////////////////////////////////////////////////////////////////
658 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
659 {
660 // Add 4-vector q to the current 4-vector
661 // Error propagation is performed automatically
662  Double_t a0=GetScalar();
663  Double_t da0=GetResultError();
664  Ali3Vector a=Get3Vector();
665  Double_t b0=q.GetScalar();
666  Double_t db0=q.GetResultError();
667  Ali3Vector b=q.Get3Vector();
668
669  Double_t c0=a0+b0;
670  Ali3Vector c=a+b;
671  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
672
673  SetVector(c0,c);
674  SetScalarError(dc0);  
675   
676  return *this;
677 }
678 ///////////////////////////////////////////////////////////////////////////
679 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
680 {
681 // Subtract 4-vector q from the current 4-vector
682 // Error propagation is performed automatically
683  Double_t a0=GetScalar();
684  Double_t da0=GetResultError();
685  Ali3Vector a=Get3Vector();
686  Double_t b0=q.GetScalar();
687  Double_t db0=q.GetResultError();
688  Ali3Vector b=q.Get3Vector();
689
690  Double_t c0=a0-b0;
691  Ali3Vector c=a-b;
692  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
693
694  SetVector(c0,c);
695  SetScalarError(dc0);  
696
697  return *this;
698 }
699 ///////////////////////////////////////////////////////////////////////////
700 Ali4Vector& Ali4Vector::operator*=(Double_t s)
701 {
702 // Multiply the current 4-vector with a scalar s
703 // Error propagation is performed automatically
704  Double_t a0=GetScalar();
705  Double_t da0=GetResultError();
706  Ali3Vector a=Get3Vector();
707
708  a0*=s;
709  a*=s;
710  da0*=s;
711
712  SetVector(a0,a);
713  SetScalarError(da0);  
714
715  return *this;
716 }
717 ///////////////////////////////////////////////////////////////////////////
718 Ali4Vector& Ali4Vector::operator/=(Double_t s)
719 {
720 // Divide the current vector by a scalar s
721 // Error propagation is performed automatically
722
723  if (fabs(s)<1.e-20) // Protect against division by 0
724  {
725   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
726   return *this;
727  }
728  else
729  {
730   Double_t a0=GetScalar();
731   Double_t da0=GetResultError();
732   Ali3Vector a=Get3Vector();
733
734   a0/=s;
735   a/=s;
736   da0/=s;
737
738   SetVector(a0,a);
739   SetScalarError(da0);  
740   
741   return *this;
742  }
743 }
744 ///////////////////////////////////////////////////////////////////////////
745 Int_t Ali4Vector::GetScalarFlag()
746 {
747 // Provide the value of the fScalar flag (for internal use only).
748  return fScalar;
749 }
750 ///////////////////////////////////////////////////////////////////////////
751 Ali3Vector Ali4Vector::GetVecTrans()
752 {
753 // Provide the transverse vector part w.r.t. z-axis.
754 // Error propagation is performed automatically
755   
756  return fV.GetVecTrans();
757 }
758 ///////////////////////////////////////////////////////////////////////////
759 Ali3Vector Ali4Vector::GetVecLong()
760 {
761 // Provide the longitudinal vector part w.r.t. z-axis.
762 // Error propagation is performed automatically
763   
764  return fV.GetVecLong();
765 }
766 ///////////////////////////////////////////////////////////////////////////
767 Double_t Ali4Vector::GetScaTrans()
768 {
769 // Provide the "transverse value" of the scalar part w.r.t. z-axis.
770 // This provides a basis for e.g. E_trans calculation.
771 // Note : the returned value is always positive or zero.
772 // The error on the value is available via GetResultError()
773 // after invokation of GetScaTrans().
774  Double_t a[3],ea[3];
775
776  fV.GetVector(a,"sph");
777  fV.GetErrors(ea,"sph");
778
779  Double_t s=GetScalar();
780  Double_t ds=GetResultError();
781
782  Double_t st,dst2;
783  st=s*sin(a[1]);
784  dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2);
785
786  fDresult=sqrt(dst2);  
787  return fabs(st);
788 }
789 ///////////////////////////////////////////////////////////////////////////
790 Double_t Ali4Vector::GetScaLong()
791 {
792 // Provide the "longitudinal value" of the scalar part w.r.t. z-axis.
793 // This provides a basis for e.g. E_long calculation.
794 // Note : the returned value can also be negative.
795 // The error on the value is available via GetResultError()
796 // after invokation of GetScaLong().
797  Double_t a[3],ea[3];
798
799  fV.GetVector(a,"sph");
800  fV.GetErrors(ea,"sph");
801
802  Double_t s=GetScalar();
803  Double_t ds=GetResultError();
804
805  Double_t sl,dsl2;
806  sl=s*cos(a[1]);
807  dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2);
808
809  fDresult=sqrt(dsl2);  
810  return sl;
811 }
812 ///////////////////////////////////////////////////////////////////////////
813 Double_t Ali4Vector::GetPseudoRapidity()
814 {
815 // Provide the pseudorapidity value of the vector part w.r.t. z-axis.
816 // The error on the value is available via GetResultError()
817 // after invokation of GetPseudoRapidity().
818  Double_t eta=fV.GetPseudoRapidity();
819  fDresult=fV.GetResultError();
820  return eta;
821 }
822 ///////////////////////////////////////////////////////////////////////////
823 Ali3Vector Ali4Vector::GetBetaVector()
824 {
825 // Provide the beta 3-vector corresponding to this 4-vector.
826  Ali3Vector beta;
827  if (fabs(fV0)>0.) beta=fV/fV0;
828  if (fabs(fDv0)>0.)
829  {
830   Double_t vecv[3],errv[3],errb[3];
831   Double_t sqerr=0;
832   fV.GetVector(vecv,"car");
833   fV.GetErrors(errv,"car");
834   for (Int_t i=0; i<3; i++)
835   {
836    sqerr=pow((errv[i]/fV0),2)+pow((vecv[i]*fDv0/(fV0*fV0)),2);
837    errb[i]=sqrt(sqerr);
838   }
839   beta.SetErrors(errb,"car");
840  }
841  return beta;
842 }
843 ///////////////////////////////////////////////////////////////////////////
844 Double_t Ali4Vector::GetBeta()
845 {
846 // Provide the beta value (i.e. v/c) corresponding to this 4-vector.
847  Ali3Vector beta=GetBetaVector();
848  Double_t val=beta.GetNorm();
849  fDresult=beta.GetResultError();
850  return val;
851 }
852 ///////////////////////////////////////////////////////////////////////////
853 Double_t Ali4Vector::GetGamma()
854 {
855 // Provide the Lorentz gamma factor corresponding to this 4-vector.
856 // In case the gamma factor is infinite a value of -1 is returned.
857  Double_t gamma=-1;
858  fDresult=0;
859  Double_t inv=sqrt(fabs(fV2));
860  if (inv>0.)
861  {
862   Double_t dinv=fDv2/(2.*inv);
863   gamma=fV0/inv;
864   Double_t sqerr=pow((fDv0/inv),2)+pow((fV0*dinv/fV2),2);
865   fDresult=sqrt(sqerr);
866  }
867  return gamma;
868 }
869 ///////////////////////////////////////////////////////////////////////////