08-mar-2003 NvE Compiler option /GR introduced for MSVC++ in mklibs.bat to explicitly...
[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::SetZero()
165 {
166 // (Re)set all attributes to zero.
167 // Note : The (de)selection of the scalar mode is not modified.
168  fV2=0;
169  fDv2=0;
170  fV0=0;
171  fDv0=0;
172  fDresult=0;
173  fV.SetZero();
174 }
175 ///////////////////////////////////////////////////////////////////////////
176 void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v)
177 {
178 // Store contravariant vector.
179 // The error on the scalar part is initialised to 0.
180 // The errors on the vector part are taken from the input Ali3Vector.
181 // Scalar mode is automatically selected. 
182 // The error on scalar result operations is reset to 0.
183  fScalar=1;
184  fV0=v0;
185  fV=v;
186  fV2=pow(v0,2)-fV.Dot(fV);
187  SetScalarError(0);
188 }
189 ///////////////////////////////////////////////////////////////////////////
190 void Ali4Vector::SetVector(Double_t* v,TString f)
191 {
192 // Store vector according to reference frame f.
193 // All errors are initialised to 0.
194 // Scalar mode is automatically selected. 
195 // The error on scalar result operations is reset to 0.
196  fScalar=1;
197  Double_t a[3];
198  for (Int_t i=0; i<3; i++)
199  {
200   a[i]=v[i+1];
201  }
202  fV0=v[0];
203  fV.SetVector(a,f);
204  fV2=pow(fV0,2)-fV.Dot(fV);
205  fDv2=0;
206  fDv0=0;
207  fDresult=0;
208 }
209 ///////////////////////////////////////////////////////////////////////////
210 void Ali4Vector::GetVector(Double_t* v,TString f)
211 {
212 // Provide 4-vector components according to reference frame f
213 // and according to the current mode.
214 // Scalar mode    : The scalar part is directly returned via v[0].
215 // Invariant mode : The scalar part is re-calculated via the value
216 //                  of the Lorentz invariant and then returned via v[0].
217  if (fScalar)
218  {
219   v[0]=fV0;
220  }
221  else
222  {
223   v[0]=sqrt(fV.Dot(fV)+fV2);
224  } 
225  Double_t a[3];
226  fV.GetVector(a,f);
227  for (Int_t i=0; i<3; i++)
228  {
229   v[i+1]=a[i];
230  }
231 }
232 ///////////////////////////////////////////////////////////////////////////
233 void Ali4Vector::SetVector(Float_t* v,TString f)
234 {
235 // Store vector according to reference frame f.
236 // All errors are initialised to 0.
237 // Scalar mode is automatically selected. 
238 // The error on scalar result operations is reset to 0.
239  Double_t vec[4];
240  for (Int_t i=0; i<4; i++)
241  {
242   vec[i]=v[i];
243  }
244  SetVector(vec,f);
245 }
246 ///////////////////////////////////////////////////////////////////////////
247 void Ali4Vector::GetVector(Float_t* v,TString f)
248 {
249 // Provide 4-vector components according to reference frame f
250 // and according to the current mode.
251 // Scalar mode    : The scalar part is directly returned via v[0].
252 // Invariant mode : The scalar part is re-calculated via the value
253 //                  of the Lorentz invariant and then returned via v[0].
254  Double_t vec[4];
255  GetVector(vec,f);
256  for (Int_t i=0; i<4; i++)
257  {
258   v[i]=vec[i];
259  }
260 }
261 ///////////////////////////////////////////////////////////////////////////
262 Double_t Ali4Vector::GetScalar()
263 {
264 // Provide the scalar part.
265 // The error on the scalar value is available via GetResultError()
266 // after invokation of GetScalar().
267  if (fScalar)
268  {
269   fDresult=fDv0;
270   return fV0;
271  }
272  else
273  {
274   Double_t dot=fV.Dot(fV);
275   Double_t ddot=fV.GetResultError();
276   Double_t v02=dot+fV2;
277   Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2));
278   Double_t v0=sqrt(fabs(v02));
279   Double_t dv0=0;
280   if (v0) dv0=dv02/(2.*v0);
281   fDresult=dv0;
282   return v0;
283  }
284 }
285 ///////////////////////////////////////////////////////////////////////////
286 Double_t Ali4Vector::GetResultError()
287 {
288 // Provide the error on the result of an operation yielding a scalar
289 // E.g. GetScalar(), GetInvariant() or Dot()
290  return fDresult;
291 }
292 ///////////////////////////////////////////////////////////////////////////
293 void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
294 {
295 // Modify the scalar part (v0) and its error (dv0).
296 // The default value for dv0 is 0.
297 // The vector part is not modified. 
298 // Scalar mode is automatically selected
299 // ==> Lorentz invariant and its error are updated. 
300 // The error on scalar result operations is reset to 0.
301  fScalar=1;
302  fV0=v0;
303  fV2=pow(v0,2)-fV.Dot(fV);
304  SetScalarError(dv0);
305 }
306 ///////////////////////////////////////////////////////////////////////////
307 void Ali4Vector::SetScalarError(Double_t dv0)
308 {
309 // Set the error on the scalar part.
310 // If in scalar mode, update error on the invariant accordingly.
311 // The error on scalar result operations is reset to 0.
312  fDv0=dv0;
313  if (fScalar)
314  {
315   Double_t norm=fV.GetNorm();
316   Double_t dnorm=fV.GetResultError();
317   fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2));
318  } 
319  fDresult=0;
320 }
321 ///////////////////////////////////////////////////////////////////////////
322 void Ali4Vector::Set3Vector(Ali3Vector& v)
323 {
324 // Set the 3-vector part, the errors are taken from the input Ali3Vector
325 // Scalar mode    : The scalar part and its error are not modified,
326 //                  the Lorentz invariantand its error are re-calculated.
327 // Invariant mode : The Lorentz invariant and its error are not modified,
328 //                  the scalar part and its error are re-calculated.
329 // The error on scalar result operations is reset to 0.
330  fV=v;
331  if (fScalar)
332  {
333   SetScalar(fV0,fDv0);
334  }
335  else
336  {
337   SetInvariant(fV2,fDv2);
338  }
339 }
340 ///////////////////////////////////////////////////////////////////////////
341 void Ali4Vector::Set3Vector(Double_t* v,TString f)
342 {
343 // Set the 3-vector part according to reference frame f
344 // The errors on the vector part are initialised to 0
345 // Scalar mode    : The scalar part and its error are not modified,
346 //                  the Lorentz invariantand its error are re-calculated.
347 // Invariant mode : The Lorentz invariant and its error are not modified,
348 //                  the scalar part and its error are re-calculated.
349 // The error on scalar result operations is reset to 0.
350  Double_t a[3];
351  for (Int_t i=0; i<3; i++)
352  {
353   a[i]=v[i];
354  }
355  fV.SetVector(a,f);
356
357  if (fScalar)
358  {
359   SetScalar(fV0,fDv0);
360  }
361  else
362  {
363   SetInvariant(fV2,fDv2);
364  }
365 }
366 ///////////////////////////////////////////////////////////////////////////
367 void Ali4Vector::Set3Vector(Float_t* v,TString f)
368 {
369 // Set the 3-vector part according to reference frame f
370 // The errors on the vector part are initialised to 0
371 // The Lorentz invariant is not modified
372 // The error on scalar result operations is reset to 0.
373  Double_t vec[3];
374  for (Int_t i=0; i<3; i++)
375  {
376   vec[i]=v[i];
377  }
378  Set3Vector(vec,f);
379 }
380 ///////////////////////////////////////////////////////////////////////////
381 void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
382 {
383 // Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2).
384 // The default value for the error dv2 is 0.
385 // The vector part is not modified.
386 // Invariant mode is automatically selected
387 // ==> the scalar part and its error are updated.
388 // The error on scalar result operations is reset to 0.
389 //
390  fScalar=0;
391  fV2=v2;
392  fDv2=dv2;
393  fV0=GetScalar();
394  fDv0=GetResultError();
395  fDresult=0;
396 }
397 ///////////////////////////////////////////////////////////////////////////
398 void Ali4Vector::SetInvariantError(Double_t dv2)
399 {
400 // Set the error on the Lorentz invariant.
401 // If in invariant mode, update error on the scalar part accordingly.
402 // The error on scalar result operations is reset to 0.
403  fDv2=dv2;
404  if (!fScalar)
405  {
406   fDv0=GetResultError();
407  }
408  fDresult=0; 
409 }
410 ///////////////////////////////////////////////////////////////////////////
411 Double_t Ali4Vector::GetInvariant()
412 {
413 // Provide the Lorentz invariant v^i*v_i.
414 // The error on the Lorentz invariant is available via GetResultError()
415 // after invokation of GetInvariant().
416  if (!fScalar)
417  {
418   fDresult=fDv2;
419   return fV2;
420  }
421  else
422  {
423   Double_t inv=Dot(*this);
424   return inv;
425  }
426 }
427 ///////////////////////////////////////////////////////////////////////////
428 Ali3Vector Ali4Vector::Get3Vector()
429 {
430 // Provide the 3-vector part
431  return fV;
432 }
433 ///////////////////////////////////////////////////////////////////////////
434 void Ali4Vector::SetErrors(Double_t* e,TString f)
435 {
436 // Store errors for vector v^i according to reference frame f
437 // If in scalar mode, update error on the invariant accordingly.
438 // The error on scalar result operations is reset to 0.
439  Double_t a[3];
440  for (Int_t i=0; i<3; i++)
441  {
442   a[i]=e[i+1];
443  }
444  SetScalarError(e[0]);
445  fV.SetErrors(a,f);
446 }
447 ///////////////////////////////////////////////////////////////////////////
448 void Ali4Vector::SetErrors(Float_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[4];
454  for (Int_t i=0; i<4; i++)
455  {
456   a[i]=e[i];
457  }
458  SetErrors(a,f);
459 }
460 ///////////////////////////////////////////////////////////////////////////
461 void Ali4Vector::GetErrors(Double_t* e,TString f)
462 {
463 // Provide errors for vector v^i according to reference frame f
464 // and according to the current mode.
465 // Scalar mode    : The error on the scalar part is directly returned via e[0].
466 // Invariant mode : The error on the scalar part is re-calculated via the error
467 //                  value on the Lorentz invariant and then returned via e[0].
468  Double_t a[3];
469  fV.GetErrors(a,f);
470
471  e[0]=GetResultError();
472  for (Int_t i=0; i<3; i++)
473  {
474   e[i+1]=a[i];
475  }
476 }
477 ///////////////////////////////////////////////////////////////////////////
478 void Ali4Vector::GetErrors(Float_t* e,TString f)
479 {
480 // Provide errors for vector v^i according to reference frame f
481 // and according to the current mode.
482 // Scalar mode    : The error on the scalar part is directly returned via e[0].
483 // Invariant mode : The error on the scalar part is re-calculated via the error
484 //                  value on the Lorentz invariant and then returned via e[0].
485  Double_t a[4];
486  GetErrors(a,f);
487  for (Int_t i=0; i<4; i++)
488  {
489   e[i]=a[i];
490  }
491 }
492 ///////////////////////////////////////////////////////////////////////////
493 void Ali4Vector::Data(TString f)
494 {
495 // Print contravariant vector components and errors according to
496 // reference frame f and according to the current mode.
497 // Scalar mode    : The scalar part and its error are directly returned.
498 // Invariant mode : The scalar part and its error are re-calculated via the
499 //                  value (and error) of the Lorentz invariant.
500
501  if (f=="car" || f=="sph" || f=="cyl")
502  {
503   Double_t vec[4],err[4];
504   GetVector(vec,f);
505   GetErrors(err,f);
506   Double_t inv=GetInvariant();
507   Double_t dinv=GetResultError();
508   cout << " Contravariant vector in " << f.Data() << " coordinates : "
509        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl; 
510   cout << " ------------- Errors in " << f.Data() << " coordinates : "
511        << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl; 
512   cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
513  }
514  else
515  {
516   cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl
517        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl; 
518  }
519 }
520 ///////////////////////////////////////////////////////////////////////////
521 Double_t Ali4Vector::Dot(Ali4Vector& q)
522 {
523 // Provide the dot product of the current vector with vector q
524  Double_t dotpro=0;
525  Double_t a0=GetScalar();
526  Double_t da0=GetResultError();
527  if ((this) == &q) // Check for special case v.Dot(v)
528  {
529   Double_t norm=fV.GetNorm();
530   Double_t dnorm=fV.GetResultError();
531   dotpro=pow(a0,2)-pow(norm,2);
532   fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2));
533  }
534  else
535  {
536   Double_t b0=q.GetScalar();
537   Double_t db0=q.GetResultError();
538   Ali3Vector b=q.Get3Vector();
539
540   Double_t dot=fV.Dot(b);
541   Double_t ddot=fV.GetResultError();
542
543   dotpro=a0*b0-dot;
544
545   fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
546  }
547
548  return dotpro;
549 }
550 ///////////////////////////////////////////////////////////////////////////
551 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
552 {
553 // Add 4-vector q to the current 4-vector
554 // Error propagation is performed automatically
555  Double_t a0=GetScalar();
556  Double_t da0=GetResultError();
557  Ali3Vector a=Get3Vector();
558  Double_t b0=q.GetScalar();
559  Double_t db0=q.GetResultError();
560  Ali3Vector b=q.Get3Vector();
561
562  Double_t c0=a0+b0;
563  Ali3Vector c=a+b;
564  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
565
566  Ali4Vector v;
567  v.SetVector(c0,c);
568  v.SetScalarError(dc0);  
569  return v;
570 }
571 ///////////////////////////////////////////////////////////////////////////
572 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
573 {
574 // Subtract 4-vector q from the current 4-vector
575 // Error propagation is performed automatically
576  Double_t a0=GetScalar();
577  Double_t da0=GetResultError();
578  Ali3Vector a=Get3Vector();
579  Double_t b0=q.GetScalar();
580  Double_t db0=q.GetResultError();
581  Ali3Vector b=q.Get3Vector();
582
583  Double_t c0=a0-b0;
584  Ali3Vector c=a-b;
585  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
586
587  Ali4Vector v;
588  v.SetVector(c0,c);
589  v.SetScalarError(dc0);  
590  return v;
591 }
592 ///////////////////////////////////////////////////////////////////////////
593 Ali4Vector Ali4Vector::operator*(Double_t s)
594 {
595 // Multiply the current 4-vector with a scalar s
596 // Error propagation is performed automatically
597  Double_t a0=GetScalar();
598  Double_t da0=GetResultError();
599  Ali3Vector a=Get3Vector();
600
601  a0*=s;
602  a*=s;
603  da0*=s;
604
605  Ali4Vector v;
606  v.SetVector(a0,a);
607  v.SetScalarError(da0);  
608   
609  return v;
610 }
611 ///////////////////////////////////////////////////////////////////////////
612 Ali4Vector Ali4Vector::operator/(Double_t s)
613 {
614 // Divide the current vector by a scalar s
615 // Error propagation is performed automatically
616
617  if (fabs(s)<1.e-20) // Protect against division by 0
618  {
619   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
620   return *this;
621  }
622  else
623  {
624   Double_t a0=GetScalar();
625   Double_t da0=GetResultError();
626   Ali3Vector a=Get3Vector();
627
628   a0/=s;
629   a/=s;
630   da0/=s;
631
632   Ali4Vector v;
633   v.SetVector(a0,a);
634   v.SetScalarError(da0);  
635   
636   return v;
637  }
638 }
639 ///////////////////////////////////////////////////////////////////////////
640 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
641 {
642 // Add 4-vector q to the current 4-vector
643 // Error propagation is performed automatically
644  Double_t a0=GetScalar();
645  Double_t da0=GetResultError();
646  Ali3Vector a=Get3Vector();
647  Double_t b0=q.GetScalar();
648  Double_t db0=q.GetResultError();
649  Ali3Vector b=q.Get3Vector();
650
651  Double_t c0=a0+b0;
652  Ali3Vector c=a+b;
653  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
654
655  SetVector(c0,c);
656  SetScalarError(dc0);  
657   
658  return *this;
659 }
660 ///////////////////////////////////////////////////////////////////////////
661 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
662 {
663 // Subtract 4-vector q from the current 4-vector
664 // Error propagation is performed automatically
665  Double_t a0=GetScalar();
666  Double_t da0=GetResultError();
667  Ali3Vector a=Get3Vector();
668  Double_t b0=q.GetScalar();
669  Double_t db0=q.GetResultError();
670  Ali3Vector b=q.Get3Vector();
671
672  Double_t c0=a0-b0;
673  Ali3Vector c=a-b;
674  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
675
676  SetVector(c0,c);
677  SetScalarError(dc0);  
678
679  return *this;
680 }
681 ///////////////////////////////////////////////////////////////////////////
682 Ali4Vector& Ali4Vector::operator*=(Double_t s)
683 {
684 // Multiply the current 4-vector with a scalar s
685 // Error propagation is performed automatically
686  Double_t a0=GetScalar();
687  Double_t da0=GetResultError();
688  Ali3Vector a=Get3Vector();
689
690  a0*=s;
691  a*=s;
692  da0*=s;
693
694  SetVector(a0,a);
695  SetScalarError(da0);  
696
697  return *this;
698 }
699 ///////////////////////////////////////////////////////////////////////////
700 Ali4Vector& Ali4Vector::operator/=(Double_t s)
701 {
702 // Divide the current vector by a scalar s
703 // Error propagation is performed automatically
704
705  if (fabs(s)<1.e-20) // Protect against division by 0
706  {
707   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
708   return *this;
709  }
710  else
711  {
712   Double_t a0=GetScalar();
713   Double_t da0=GetResultError();
714   Ali3Vector a=Get3Vector();
715
716   a0/=s;
717   a/=s;
718   da0/=s;
719
720   SetVector(a0,a);
721   SetScalarError(da0);  
722   
723   return *this;
724  }
725 }
726 ///////////////////////////////////////////////////////////////////////////
727 Int_t Ali4Vector::GetScalarFlag()
728 {
729 // Provide the value of the fScalar flag (for internal use only).
730  return fScalar;
731 }
732 ///////////////////////////////////////////////////////////////////////////
733 Ali3Vector Ali4Vector::GetVecTrans()
734 {
735 // Provide the transverse vector part w.r.t. z-axis.
736 // Error propagation is performed automatically
737   
738  return fV.GetVecTrans();
739 }
740 ///////////////////////////////////////////////////////////////////////////
741 Ali3Vector Ali4Vector::GetVecLong()
742 {
743 // Provide the longitudinal vector part w.r.t. z-axis.
744 // Error propagation is performed automatically
745   
746  return fV.GetVecLong();
747 }
748 ///////////////////////////////////////////////////////////////////////////
749 Double_t Ali4Vector::GetScaTrans()
750 {
751 // Provide the "transverse value" of the scalar part w.r.t. z-axis.
752 // This provides a basis for e.g. E_trans calculation.
753 // Note : the returned value is always positive or zero.
754 // The error on the value is available via GetResultError()
755 // after invokation of GetScaTrans().
756  Double_t a[3],ea[3];
757
758  fV.GetVector(a,"sph");
759  fV.GetErrors(ea,"sph");
760
761  Double_t s=GetScalar();
762  Double_t ds=GetResultError();
763
764  Double_t st,dst2;
765  st=s*sin(a[1]);
766  dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2);
767
768  fDresult=sqrt(dst2);  
769  return fabs(st);
770 }
771 ///////////////////////////////////////////////////////////////////////////
772 Double_t Ali4Vector::GetScaLong()
773 {
774 // Provide the "longitudinal value" of the scalar part w.r.t. z-axis.
775 // This provides a basis for e.g. E_long calculation.
776 // Note : the returned value can also be negative.
777 // The error on the value is available via GetResultError()
778 // after invokation of GetScaLong().
779  Double_t a[3],ea[3];
780
781  fV.GetVector(a,"sph");
782  fV.GetErrors(ea,"sph");
783
784  Double_t s=GetScalar();
785  Double_t ds=GetResultError();
786
787  Double_t sl,dsl2;
788  sl=s*cos(a[1]);
789  dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2);
790
791  fDresult=sqrt(dsl2);  
792  return sl;
793 }
794 ///////////////////////////////////////////////////////////////////////////
795 Double_t Ali4Vector::GetPseudoRapidity()
796 {
797 // Provide the pseudorapidity value of the vector part w.r.t. z-axis.
798 // The error on the value is available via GetResultError()
799 // after invokation of GetPseudoRapidity().
800  Double_t eta=fV.GetPseudoRapidity();
801  fDresult=fV.GetResultError();
802  return eta;
803 }
804 ///////////////////////////////////////////////////////////////////////////