1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////
20 // Handling of Lorentz 4-vectors in various reference frames.
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.
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.
30 // A 4-vector is said to have a scalar part and a 3-vector part,
31 // which is indicated by the notation
33 // x^i = (x^0,x^1,x^2,x^3)
35 // The scalar part = x^0
36 // The 3-vector part = (x^1,x^2,x^3)
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.
42 // This allows the following two modes of functionality :
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.
51 // The philosophy followed here is the following :
52 // ===============================================
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.
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
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.
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.
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).
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 :
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)
97 // All angles are in radians.
104 // Float_t v[4]={25,-1,3,7};
105 // a.SetVector(v,"car");
108 // a.GetVector(vec,"sph");
111 // Float_t v2[4]={33,6,-18,2};
112 // b.SetVector(v2,"car");
114 // Float_t dotpro=a.Dot(b);
118 // Float_t vec2[3]={1,2,3};
119 // x.SetVector(vec2,"car");
122 // c.SetVector(x0,x);
123 // c.GetVector(vec,"car");
129 //--- Author: Nick van Eijndhoven 01-apr-1999 UU-SAP Utrecht
130 //- Modified: NvE $Date$ UU-SAP Utrecht
131 ///////////////////////////////////////////////////////////////////////////
133 #include "Ali4Vector.h"
134 #include "Riostream.h"
136 ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
138 Ali4Vector::Ali4Vector()
140 // Creation of a contravariant 4-vector and initialisation of parameters.
141 // All values are initialised to 0.
142 // Scalar mode is initially selected.
146 ///////////////////////////////////////////////////////////////////////////
147 Ali4Vector::~Ali4Vector()
149 // Destructor to delete dynamically allocated memory
151 ///////////////////////////////////////////////////////////////////////////
152 Ali4Vector::Ali4Vector(const Ali4Vector& v)
163 ///////////////////////////////////////////////////////////////////////////
164 void Ali4Vector::Load(Ali4Vector& q)
166 // Load all attributes of the input Ali4Vector into this Ali4Vector object.
167 Int_t temp1=q.GetScalarFlag();
168 Double_t temp2=q.GetResultError();
170 q.GetVector(a,"sph");
172 q.GetErrors(a,"car");
177 ///////////////////////////////////////////////////////////////////////////
178 void Ali4Vector::SetZero()
180 // (Re)set all attributes to zero.
181 // Note : The (de)selection of the scalar mode is not modified.
189 ///////////////////////////////////////////////////////////////////////////
190 void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v)
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.
200 fV2=pow(v0,2)-fV.Dot(fV);
203 ///////////////////////////////////////////////////////////////////////////
204 void Ali4Vector::SetVector(Double_t* v,TString f)
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.
212 for (Int_t i=0; i<3; i++)
218 fV2=pow(fV0,2)-fV.Dot(fV);
223 ///////////////////////////////////////////////////////////////////////////
224 void Ali4Vector::GetVector(Double_t* v,TString f)
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].
237 v[0]=sqrt(fV.Dot(fV)+fV2);
241 for (Int_t i=0; i<3; i++)
246 ///////////////////////////////////////////////////////////////////////////
247 void Ali4Vector::SetVector(Float_t* v,TString f)
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.
254 for (Int_t i=0; i<4; i++)
260 ///////////////////////////////////////////////////////////////////////////
261 void Ali4Vector::GetVector(Float_t* v,TString f)
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].
270 for (Int_t i=0; i<4; i++)
275 ///////////////////////////////////////////////////////////////////////////
276 Double_t Ali4Vector::GetScalar()
278 // Provide the scalar part.
279 // The error on the scalar value is available via GetResultError()
280 // after invokation of GetScalar().
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));
294 if (v0) dv0=dv02/(2.*v0);
299 ///////////////////////////////////////////////////////////////////////////
300 Double_t Ali4Vector::GetResultError()
302 // Provide the error on the result of an operation yielding a scalar
303 // E.g. GetScalar(), GetInvariant() or Dot()
306 ///////////////////////////////////////////////////////////////////////////
307 void Ali4Vector::SetScalar(Double_t v0,Double_t dv0)
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.
317 fV2=pow(v0,2)-fV.Dot(fV);
320 ///////////////////////////////////////////////////////////////////////////
321 void Ali4Vector::SetScalarError(Double_t dv0)
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.
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));
335 ///////////////////////////////////////////////////////////////////////////
336 void Ali4Vector::Set3Vector(Ali3Vector& v)
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.
351 SetInvariant(fV2,fDv2);
354 ///////////////////////////////////////////////////////////////////////////
355 void Ali4Vector::Set3Vector(Double_t* v,TString f)
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.
365 for (Int_t i=0; i<3; i++)
377 SetInvariant(fV2,fDv2);
380 ///////////////////////////////////////////////////////////////////////////
381 void Ali4Vector::Set3Vector(Float_t* v,TString f)
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.
388 for (Int_t i=0; i<3; i++)
394 ///////////////////////////////////////////////////////////////////////////
395 void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2)
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.
408 fDv0=GetResultError();
411 ///////////////////////////////////////////////////////////////////////////
412 void Ali4Vector::SetInvariantError(Double_t dv2)
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.
420 fDv0=GetResultError();
424 ///////////////////////////////////////////////////////////////////////////
425 Double_t Ali4Vector::GetInvariant()
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().
437 Double_t inv=Dot(*this);
441 ///////////////////////////////////////////////////////////////////////////
442 Ali3Vector Ali4Vector::Get3Vector()
444 // Provide the 3-vector part
447 ///////////////////////////////////////////////////////////////////////////
448 void Ali4Vector::SetErrors(Double_t* e,TString f)
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.
454 for (Int_t i=0; i<3; i++)
458 SetScalarError(e[0]);
461 ///////////////////////////////////////////////////////////////////////////
462 void Ali4Vector::SetErrors(Float_t* e,TString f)
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.
468 for (Int_t i=0; i<4; i++)
474 ///////////////////////////////////////////////////////////////////////////
475 void Ali4Vector::GetErrors(Double_t* e,TString f)
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].
485 // Dummy invokation of GetScalar to obtain automatic proper error determination
488 e[0]=GetResultError();
490 for (Int_t i=0; i<3; i++)
495 ///////////////////////////////////////////////////////////////////////////
496 void Ali4Vector::GetErrors(Float_t* e,TString f)
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].
505 for (Int_t i=0; i<4; i++)
510 ///////////////////////////////////////////////////////////////////////////
511 void Ali4Vector::Data(TString f)
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.
519 if (f=="car" || f=="sph" || f=="cyl")
521 Double_t vec[4],err[4];
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;
534 cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl
535 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
538 ///////////////////////////////////////////////////////////////////////////
539 Double_t Ali4Vector::Dot(Ali4Vector& q)
541 // Provide the dot product of the current vector with vector q
543 Double_t a0=GetScalar();
544 Double_t da0=GetResultError();
545 if ((this) == &q) // Check for special case v.Dot(v)
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));
554 Double_t b0=q.GetScalar();
555 Double_t db0=q.GetResultError();
556 Ali3Vector b=q.Get3Vector();
558 Double_t dot=fV.Dot(b);
559 Double_t ddot=fV.GetResultError();
563 fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
568 ///////////////////////////////////////////////////////////////////////////
569 Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
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();
582 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
586 v.SetScalarError(dc0);
589 ///////////////////////////////////////////////////////////////////////////
590 Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
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();
603 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
607 v.SetScalarError(dc0);
610 ///////////////////////////////////////////////////////////////////////////
611 Ali4Vector Ali4Vector::operator*(Double_t s)
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();
625 v.SetScalarError(da0);
629 ///////////////////////////////////////////////////////////////////////////
630 Ali4Vector Ali4Vector::operator/(Double_t s)
632 // Divide the current vector by a scalar s
633 // Error propagation is performed automatically
635 if (fabs(s)<1.e-20) // Protect against division by 0
637 cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
642 Double_t a0=GetScalar();
643 Double_t da0=GetResultError();
644 Ali3Vector a=Get3Vector();
652 v.SetScalarError(da0);
657 ///////////////////////////////////////////////////////////////////////////
658 Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
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();
671 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
678 ///////////////////////////////////////////////////////////////////////////
679 Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
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();
692 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
699 ///////////////////////////////////////////////////////////////////////////
700 Ali4Vector& Ali4Vector::operator*=(Double_t s)
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();
717 ///////////////////////////////////////////////////////////////////////////
718 Ali4Vector& Ali4Vector::operator/=(Double_t s)
720 // Divide the current vector by a scalar s
721 // Error propagation is performed automatically
723 if (fabs(s)<1.e-20) // Protect against division by 0
725 cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
730 Double_t a0=GetScalar();
731 Double_t da0=GetResultError();
732 Ali3Vector a=Get3Vector();
744 ///////////////////////////////////////////////////////////////////////////
745 Int_t Ali4Vector::GetScalarFlag()
747 // Provide the value of the fScalar flag (for internal use only).
750 ///////////////////////////////////////////////////////////////////////////
751 Ali3Vector Ali4Vector::GetVecTrans()
753 // Provide the transverse vector part w.r.t. z-axis.
754 // Error propagation is performed automatically
756 return fV.GetVecTrans();
758 ///////////////////////////////////////////////////////////////////////////
759 Ali3Vector Ali4Vector::GetVecLong()
761 // Provide the longitudinal vector part w.r.t. z-axis.
762 // Error propagation is performed automatically
764 return fV.GetVecLong();
766 ///////////////////////////////////////////////////////////////////////////
767 Double_t Ali4Vector::GetScaTrans()
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().
776 fV.GetVector(a,"sph");
777 fV.GetErrors(ea,"sph");
779 Double_t s=GetScalar();
780 Double_t ds=GetResultError();
784 dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2);
789 ///////////////////////////////////////////////////////////////////////////
790 Double_t Ali4Vector::GetScaLong()
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().
799 fV.GetVector(a,"sph");
800 fV.GetErrors(ea,"sph");
802 Double_t s=GetScalar();
803 Double_t ds=GetResultError();
807 dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2);
812 ///////////////////////////////////////////////////////////////////////////
813 Double_t Ali4Vector::GetPseudoRapidity()
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();
822 ///////////////////////////////////////////////////////////////////////////
823 Ali3Vector Ali4Vector::GetBetaVector()
825 // Provide the beta 3-vector corresponding to this 4-vector.
827 if (fabs(fV0)>0.) beta=fV/fV0;
830 Double_t vecv[3],errv[3],errb[3];
832 fV.GetVector(vecv,"car");
833 fV.GetErrors(errv,"car");
834 for (Int_t i=0; i<3; i++)
836 sqerr=pow((errv[i]/fV0),2)+pow((vecv[i]*fDv0/(fV0*fV0)),2);
839 beta.SetErrors(errb,"car");
843 ///////////////////////////////////////////////////////////////////////////
844 Double_t Ali4Vector::GetBeta()
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();
852 ///////////////////////////////////////////////////////////////////////////
853 Double_t Ali4Vector::GetGamma()
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.
859 Double_t inv=sqrt(fabs(fV2));
862 Double_t dinv=fDv2/(2.*inv);
864 Double_t sqerr=pow((fDv0/inv),2)+pow((fV0*dinv/fV2),2);
865 fDresult=sqrt(sqerr);
869 ///////////////////////////////////////////////////////////////////////////