27-may-2001 NvE New class AliEvent introduced and RALICELinkDef.h & co. updated accor...
[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.Info("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  
135 ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
136  
137 Ali4Vector::Ali4Vector()
138 {
139 // Creation of a contravariant 4-vector and initialisation of parameters.
140 // All values are initialised to 0.
141 // Scalar mode is initially selected.
142  fScalar=1;
143  fV2=0;
144  fDv2=0;
145  fV0=0;
146  fDv0=0;
147  fDresult=0;
148  Double_t a[3]={0,0,0};
149  fV.SetVector(a,"sph");
150 }
151 ///////////////////////////////////////////////////////////////////////////
152 Ali4Vector::~Ali4Vector()
153 {
154 // Destructor to delete dynamically allocated memory
155 }
156 ///////////////////////////////////////////////////////////////////////////
157 void Ali4Vector::SetVector(Double_t v0,Ali3Vector v)
158 {
159 // Store contravariant vector.
160 // The error on the scalar part is initialised to 0.
161 // The errors on the vector part are taken from the input Ali3Vector.
162 // Scalar mode is automatically selected. 
163 // The error on scalar result operations is reset to 0.
164  fScalar=1;
165  fV0=v0;
166  fV=v;
167  fV2=pow(v0,2)-fV.Dot(fV);
168  SetScalarError(0);
169 }
170 ///////////////////////////////////////////////////////////////////////////
171 void Ali4Vector::SetVector(Double_t* v,TString f)
172 {
173 // Store vector according to reference frame f.
174 // All errors are initialised to 0.
175 // Scalar mode is automatically selected. 
176 // The error on scalar result operations is reset to 0.
177  fScalar=1;
178  Double_t a[3];
179  for (Int_t i=0; i<3; i++)
180  {
181   a[i]=v[i+1];
182  }
183  fV0=v[0];
184  fV.SetVector(a,f);
185  fV2=pow(fV0,2)-fV.Dot(fV);
186  fDv2=0;
187  fDv0=0;
188  fDresult=0;
189 }
190 ///////////////////////////////////////////////////////////////////////////
191 void Ali4Vector::GetVector(Double_t* v,TString f)
192 {
193 // Provide 4-vector components according to reference frame f
194 // and according to the current mode.
195 // Scalar mode    : The scalar part is directly returned via v[0].
196 // Invariant mode : The scalar part is re-calculated via the value
197 //                  of the Lorentz invariant and then returned via v[0].
198  if (fScalar)
199  {
200   v[0]=fV0;
201  }
202  else
203  {
204   v[0]=sqrt(fV.Dot(fV)+fV2);
205  } 
206  Double_t a[3];
207  fV.GetVector(a,f);
208  for (Int_t i=0; i<3; i++)
209  {
210   v[i+1]=a[i];
211  }
212 }
213 ///////////////////////////////////////////////////////////////////////////
214 void Ali4Vector::SetVector(Float_t* v,TString f)
215 {
216 // Store vector according to reference frame f.
217 // All errors are initialised to 0.
218 // Scalar mode is automatically selected. 
219 // The error on scalar result operations is reset to 0.
220  Double_t vec[4];
221  for (Int_t i=0; i<4; i++)
222  {
223   vec[i]=v[i];
224  }
225  SetVector(vec,f);
226 }
227 ///////////////////////////////////////////////////////////////////////////
228 void Ali4Vector::GetVector(Float_t* v,TString f)
229 {
230 // Provide 4-vector components according to reference frame f
231 // and according to the current mode.
232 // Scalar mode    : The scalar part is directly returned via v[0].
233 // Invariant mode : The scalar part is re-calculated via the value
234 //                  of the Lorentz invariant and then returned via v[0].
235  Double_t vec[4];
236  GetVector(vec,f);
237  for (Int_t i=0; i<4; i++)
238  {
239   v[i]=vec[i];
240  }
241 }
242 ///////////////////////////////////////////////////////////////////////////
243 Double_t Ali4Vector::GetScalar()
244 {
245 // Provide the scalar part.
246 // The error on the scalar value is available via GetResultError()
247 // after invokation of GetScalar().
248  if (fScalar)
249  {
250   fDresult=fDv0;
251   return fV0;
252  }
253  else
254  {
255   Double_t dot=fV.Dot(fV);
256   Double_t ddot=fV.GetResultError();
257   Double_t v02=dot+fV2;
258   Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2));
259   Double_t v0=sqrt(fabs(v02));
260   Double_t dv0=0;
261   if (v0) dv0=dv02/(2.*v0);
262   fDresult=dv0;
263   return v0;
264  }
265 }
266 ///////////////////////////////////////////////////////////////////////////
267 Double_t Ali4Vector::GetResultError()
268 {
269 // Provide the error on the result of an operation yielding a scalar
270 // E.g. GetScalar(), GetInvariant() or Dot()
271  return fDresult;
272 }
273 ///////////////////////////////////////////////////////////////////////////
274 void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
275 {
276 // Modify the scalar part (v0) and its error (dv0).
277 // The default value for dv0 is 0.
278 // The vector part is not modified. 
279 // Scalar mode is automatically selected
280 // ==> Lorentz invariant and its error are updated. 
281 // The error on scalar result operations is reset to 0.
282  fScalar=1;
283  fV0=v0;
284  fV2=pow(v0,2)-fV.Dot(fV);
285  SetScalarError(dv0);
286 }
287 ///////////////////////////////////////////////////////////////////////////
288 void Ali4Vector::SetScalarError(Double_t dv0)
289 {
290 // Set the error on the scalar part.
291 // If in scalar mode, update error on the invariant accordingly.
292 // The error on scalar result operations is reset to 0.
293  fDv0=dv0;
294  if (fScalar)
295  {
296   Double_t norm=fV.GetNorm();
297   Double_t dnorm=fV.GetResultError();
298   fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2));
299  } 
300  fDresult=0;
301 }
302 ///////////////////////////////////////////////////////////////////////////
303 void Ali4Vector::Set3Vector(Ali3Vector v)
304 {
305 // Set the 3-vector part, the errors are taken from the input Ali3Vector
306 // Scalar mode    : The scalar part and its error are not modified,
307 //                  the Lorentz invariantand its error are re-calculated.
308 // Invariant mode : The Lorentz invariant and its error are not modified,
309 //                  the scalar part and its error are re-calculated.
310 // The error on scalar result operations is reset to 0.
311  fV=v;
312  if (fScalar)
313  {
314   SetScalar(fV0,fDv0);
315  }
316  else
317  {
318   SetInvariant(fV2,fDv2);
319  }
320 }
321 ///////////////////////////////////////////////////////////////////////////
322 void Ali4Vector::Set3Vector(Double_t* v,TString f)
323 {
324 // Set the 3-vector part according to reference frame f
325 // The errors on the vector part are initialised to 0
326 // Scalar mode    : The scalar part and its error are not modified,
327 //                  the Lorentz invariantand its error are re-calculated.
328 // Invariant mode : The Lorentz invariant and its error are not modified,
329 //                  the scalar part and its error are re-calculated.
330 // The error on scalar result operations is reset to 0.
331  Double_t a[3];
332  for (Int_t i=0; i<3; i++)
333  {
334   a[i]=v[i];
335  }
336  fV.SetVector(a,f);
337
338  if (fScalar)
339  {
340   SetScalar(fV0,fDv0);
341  }
342  else
343  {
344   SetInvariant(fV2,fDv2);
345  }
346 }
347 ///////////////////////////////////////////////////////////////////////////
348 void Ali4Vector::Set3Vector(Float_t* v,TString f)
349 {
350 // Set the 3-vector part according to reference frame f
351 // The errors on the vector part are initialised to 0
352 // The Lorentz invariant is not modified
353 // The error on scalar result operations is reset to 0.
354  Double_t vec[3];
355  for (Int_t i=0; i<3; i++)
356  {
357   vec[i]=v[i];
358  }
359  Set3Vector(vec,f);
360 }
361 ///////////////////////////////////////////////////////////////////////////
362 void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
363 {
364 // Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2).
365 // The default value for the error dv2 is 0.
366 // The vector part is not modified.
367 // Invariant mode is automatically selected
368 // ==> the scalar part and its error are updated.
369 // The error on scalar result operations is reset to 0.
370 //
371  fScalar=0;
372  fV2=v2;
373  fDv2=dv2;
374  fV0=GetScalar();
375  fDv0=GetResultError();
376  fDresult=0;
377 }
378 ///////////////////////////////////////////////////////////////////////////
379 void Ali4Vector::SetInvariantError(Double_t dv2)
380 {
381 // Set the error on the Lorentz invariant.
382 // If in invariant mode, update error on the scalar part accordingly.
383 // The error on scalar result operations is reset to 0.
384  fDv2=dv2;
385  if (!fScalar)
386  {
387   fDv0=GetResultError();
388  }
389  fDresult=0; 
390 }
391 ///////////////////////////////////////////////////////////////////////////
392 Double_t Ali4Vector::GetInvariant()
393 {
394 // Provide the Lorentz invariant v^i*v_i.
395 // The error on the Lorentz invariant is available via GetResultError()
396 // after invokation of GetInvariant().
397  if (!fScalar)
398  {
399   fDresult=fDv2;
400   return fV2;
401  }
402  else
403  {
404   Double_t inv=Dot(*this);
405   return inv;
406  }
407 }
408 ///////////////////////////////////////////////////////////////////////////
409 Ali3Vector Ali4Vector::Get3Vector()
410 {
411 // Provide the 3-vector part
412  return fV;
413 }
414 ///////////////////////////////////////////////////////////////////////////
415 void Ali4Vector::SetErrors(Double_t* e,TString f)
416 {
417 // Store errors for vector v^i according to reference frame f
418 // If in scalar mode, update error on the invariant accordingly.
419 // The error on scalar result operations is reset to 0.
420  Double_t a[3];
421  for (Int_t i=0; i<3; i++)
422  {
423   a[i]=e[i+1];
424  }
425  SetScalarError(e[0]);
426  fV.SetErrors(a,f);
427 }
428 ///////////////////////////////////////////////////////////////////////////
429 void Ali4Vector::SetErrors(Float_t* e,TString f)
430 {
431 // Store errors for vector v^i according to reference frame f
432 // If in scalar mode, update error on the invariant accordingly.
433 // The error on scalar result operations is reset to 0.
434  Double_t a[4];
435  for (Int_t i=0; i<4; i++)
436  {
437   a[i]=e[i];
438  }
439  SetErrors(a,f);
440 }
441 ///////////////////////////////////////////////////////////////////////////
442 void Ali4Vector::GetErrors(Double_t* e,TString f)
443 {
444 // Provide errors for vector v^i according to reference frame f
445 // and according to the current mode.
446 // Scalar mode    : The error on the scalar part is directly returned via e[0].
447 // Invariant mode : The error on the scalar part is re-calculated via the error
448 //                  value on the Lorentz invariant and then returned via e[0].
449  Double_t a[3];
450  fV.GetErrors(a,f);
451
452  e[0]=GetResultError();
453  for (Int_t i=0; i<3; i++)
454  {
455   e[i+1]=a[i];
456  }
457 }
458 ///////////////////////////////////////////////////////////////////////////
459 void Ali4Vector::GetErrors(Float_t* e,TString f)
460 {
461 // Provide errors for vector v^i according to reference frame f
462 // and according to the current mode.
463 // Scalar mode    : The error on the scalar part is directly returned via e[0].
464 // Invariant mode : The error on the scalar part is re-calculated via the error
465 //                  value on the Lorentz invariant and then returned via e[0].
466  Double_t a[4];
467  GetErrors(a,f);
468  for (Int_t i=0; i<4; i++)
469  {
470   e[i]=a[i];
471  }
472 }
473 ///////////////////////////////////////////////////////////////////////////
474 void Ali4Vector::Info(TString f)
475 {
476 // Print contravariant vector components and errors according to
477 // reference frame f and according to the current mode.
478 // Scalar mode    : The scalar part and its error are directly returned.
479 // Invariant mode : The scalar part and its error are re-calculated via the
480 //                  value (and error) of the Lorentz invariant.
481
482  if (f=="car" || f=="sph" || f=="cyl")
483  {
484   Double_t vec[4],err[4];
485   GetVector(vec,f);
486   GetErrors(err,f);
487   Double_t inv=GetInvariant();
488   Double_t dinv=GetResultError();
489   cout << " Contravariant vector in " << f << " coordinates : "
490        << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl; 
491   cout << " ------------- Errors in " << f << " coordinates : "
492        << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl; 
493   cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
494  }
495  else
496  {
497   cout << " *Ali4Vector::Info* Unsupported frame : " << f << endl
498        << "  Possible frames are 'car', 'sph' and 'cyl'." << endl; 
499  }
500 }
501 ///////////////////////////////////////////////////////////////////////////
502 Double_t Ali4Vector::Dot(Ali4Vector& q)
503 {
504 // Provide the dot product of the current vector with vector q
505  Double_t dotpro=0;
506  Double_t a0=GetScalar();
507  Double_t da0=GetResultError();
508  if ((this) == &q) // Check for special case v.Dot(v)
509  {
510   Double_t norm=fV.GetNorm();
511   Double_t dnorm=fV.GetResultError();
512   dotpro=pow(a0,2)-pow(norm,2);
513   fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2));
514  }
515  else
516  {
517   Double_t b0=q.GetScalar();
518   Double_t db0=q.GetResultError();
519   Ali3Vector b=q.Get3Vector();
520
521   Double_t dot=fV.Dot(b);
522   Double_t ddot=fV.GetResultError();
523
524   dotpro=a0*b0-dot;
525
526   fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
527  }
528
529  return dotpro;
530 }
531 ///////////////////////////////////////////////////////////////////////////
532 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
533 {
534 // Add 4-vector q to the current 4-vector
535 // Error propagation is performed automatically
536  Double_t a0=GetScalar();
537  Double_t da0=GetResultError();
538  Ali3Vector a=Get3Vector();
539  Double_t b0=q.GetScalar();
540  Double_t db0=q.GetResultError();
541  Ali3Vector b=q.Get3Vector();
542
543  Double_t c0=a0+b0;
544  Ali3Vector c=a+b;
545  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
546
547  Ali4Vector v;
548  v.SetVector(c0,c);
549  v.SetScalarError(dc0);  
550  return v;
551 }
552 ///////////////////////////////////////////////////////////////////////////
553 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
554 {
555 // Subtract 4-vector q from the current 4-vector
556 // Error propagation is performed automatically
557  Double_t a0=GetScalar();
558  Double_t da0=GetResultError();
559  Ali3Vector a=Get3Vector();
560  Double_t b0=q.GetScalar();
561  Double_t db0=q.GetResultError();
562  Ali3Vector b=q.Get3Vector();
563
564  Double_t c0=a0-b0;
565  Ali3Vector c=a-b;
566  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
567
568  Ali4Vector v;
569  v.SetVector(c0,c);
570  v.SetScalarError(dc0);  
571  return v;
572 }
573 ///////////////////////////////////////////////////////////////////////////
574 Ali4Vector Ali4Vector::operator*(Double_t s)
575 {
576 // Multiply the current 4-vector with a scalar s
577 // Error propagation is performed automatically
578  Double_t a0=GetScalar();
579  Double_t da0=GetResultError();
580  Ali3Vector a=Get3Vector();
581
582  a0*=s;
583  a*=s;
584  da0*=s;
585
586  Ali4Vector v;
587  v.SetVector(a0,a);
588  v.SetScalarError(da0);  
589   
590  return v;
591 }
592 ///////////////////////////////////////////////////////////////////////////
593 Ali4Vector Ali4Vector::operator/(Double_t s)
594 {
595 // Divide the current vector by a scalar s
596 // Error propagation is performed automatically
597
598  if (fabs(s)<1.e-20) // Protect against division by 0
599  {
600   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
601   return *this;
602  }
603  else
604  {
605   Double_t a0=GetScalar();
606   Double_t da0=GetResultError();
607   Ali3Vector a=Get3Vector();
608
609   a0/=s;
610   a/=s;
611   da0/=s;
612
613   Ali4Vector v;
614   v.SetVector(a0,a);
615   v.SetScalarError(da0);  
616   
617   return v;
618  }
619 }
620 ///////////////////////////////////////////////////////////////////////////
621 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
622 {
623 // Add 4-vector q to the current 4-vector
624 // Error propagation is performed automatically
625  Double_t a0=GetScalar();
626  Double_t da0=GetResultError();
627  Ali3Vector a=Get3Vector();
628  Double_t b0=q.GetScalar();
629  Double_t db0=q.GetResultError();
630  Ali3Vector b=q.Get3Vector();
631
632  Double_t c0=a0+b0;
633  Ali3Vector c=a+b;
634  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
635
636  SetVector(c0,c);
637  SetScalarError(dc0);  
638   
639  return *this;
640 }
641 ///////////////////////////////////////////////////////////////////////////
642 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
643 {
644 // Subtract 4-vector q from the current 4-vector
645 // Error propagation is performed automatically
646  Double_t a0=GetScalar();
647  Double_t da0=GetResultError();
648  Ali3Vector a=Get3Vector();
649  Double_t b0=q.GetScalar();
650  Double_t db0=q.GetResultError();
651  Ali3Vector b=q.Get3Vector();
652
653  Double_t c0=a0-b0;
654  Ali3Vector c=a-b;
655  Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
656
657  SetVector(c0,c);
658  SetScalarError(dc0);  
659
660  return *this;
661 }
662 ///////////////////////////////////////////////////////////////////////////
663 Ali4Vector& Ali4Vector::operator*=(Double_t s)
664 {
665 // Multiply the current 4-vector with a scalar s
666 // Error propagation is performed automatically
667  Double_t a0=GetScalar();
668  Double_t da0=GetResultError();
669  Ali3Vector a=Get3Vector();
670
671  a0*=s;
672  a*=s;
673  da0*=s;
674
675  SetVector(a0,a);
676  SetScalarError(da0);  
677
678  return *this;
679 }
680 ///////////////////////////////////////////////////////////////////////////
681 Ali4Vector& Ali4Vector::operator/=(Double_t s)
682 {
683 // Divide the current vector by a scalar s
684 // Error propagation is performed automatically
685
686  if (fabs(s)<1.e-20) // Protect against division by 0
687  {
688   cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
689   return *this;
690  }
691  else
692  {
693   Double_t a0=GetScalar();
694   Double_t da0=GetResultError();
695   Ali3Vector a=Get3Vector();
696
697   a0/=s;
698   a/=s;
699   da0/=s;
700
701   SetVector(a0,a);
702   SetScalarError(da0);  
703   
704   return *this;
705  }
706 }
707 ///////////////////////////////////////////////////////////////////////////
708 Int_t Ali4Vector::GetScalarFlag()
709 {
710 // Provide the value of the fScalar flag (for internal use only).
711  return fScalar;
712 }
713 ///////////////////////////////////////////////////////////////////////////