Access to the number of associated clusters (M.Ivanov)
[u/mrichter/AliRoot.git] / RALICE / Ali4Vector.cxx
CommitLineData
4c039060 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
f531a546 16// $Id$
4c039060 17
959fbac5 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");
84bb7c66 124// c.Data("cyl");
959fbac5 125// c=a+b;
126// c=a-b;
127// c=a*5;
128//
129//--- Author: Nick van Eijndhoven 01-apr-1999 UU-SAP Utrecht
f531a546 130//- Modified: NvE $Date$ UU-SAP Utrecht
959fbac5 131///////////////////////////////////////////////////////////////////////////
132
d88f97cc 133#include "Ali4Vector.h"
c72198f1 134#include "Riostream.h"
d88f97cc 135
136ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O
137
138Ali4Vector::Ali4Vector()
139{
959fbac5 140// Creation of a contravariant 4-vector and initialisation of parameters.
141// All values are initialised to 0.
142// Scalar mode is initially selected.
c72198f1 143 SetZero();
959fbac5 144 fScalar=1;
d88f97cc 145}
146///////////////////////////////////////////////////////////////////////////
147Ali4Vector::~Ali4Vector()
148{
149// Destructor to delete dynamically allocated memory
150}
151///////////////////////////////////////////////////////////////////////////
c72198f1 152Ali4Vector::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///////////////////////////////////////////////////////////////////////////
1fbffa23 164void 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///////////////////////////////////////////////////////////////////////////
c72198f1 178void 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///////////////////////////////////////////////////////////////////////////
190void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v)
d88f97cc 191{
959fbac5 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;
d88f97cc 198 fV0=v0;
199 fV=v;
959fbac5 200 fV2=pow(v0,2)-fV.Dot(fV);
201 SetScalarError(0);
d88f97cc 202}
203///////////////////////////////////////////////////////////////////////////
204void Ali4Vector::SetVector(Double_t* v,TString f)
205{
959fbac5 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;
d88f97cc 211 Double_t a[3];
212 for (Int_t i=0; i<3; i++)
213 {
214 a[i]=v[i+1];
215 }
959fbac5 216 fV0=v[0];
d88f97cc 217 fV.SetVector(a,f);
959fbac5 218 fV2=pow(fV0,2)-fV.Dot(fV);
219 fDv2=0;
220 fDv0=0;
221 fDresult=0;
d88f97cc 222}
223///////////////////////////////////////////////////////////////////////////
224void Ali4Vector::GetVector(Double_t* v,TString f)
225{
959fbac5 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 }
d88f97cc 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///////////////////////////////////////////////////////////////////////////
247void Ali4Vector::SetVector(Float_t* v,TString f)
248{
959fbac5 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.
d88f97cc 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///////////////////////////////////////////////////////////////////////////
261void Ali4Vector::GetVector(Float_t* v,TString f)
262{
959fbac5 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].
d88f97cc 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///////////////////////////////////////////////////////////////////////////
276Double_t Ali4Vector::GetScalar()
277{
959fbac5 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///////////////////////////////////////////////////////////////////////////
300Double_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///////////////////////////////////////////////////////////////////////////
307void 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///////////////////////////////////////////////////////////////////////////
321void 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///////////////////////////////////////////////////////////////////////////
c72198f1 336void Ali4Vector::Set3Vector(Ali3Vector& v)
959fbac5 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///////////////////////////////////////////////////////////////////////////
355void 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///////////////////////////////////////////////////////////////////////////
381void 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///////////////////////////////////////////////////////////////////////////
395void 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///////////////////////////////////////////////////////////////////////////
412void 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///////////////////////////////////////////////////////////////////////////
425Double_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 }
d88f97cc 440}
441///////////////////////////////////////////////////////////////////////////
442Ali3Vector Ali4Vector::Get3Vector()
443{
444// Provide the 3-vector part
445 return fV;
446}
447///////////////////////////////////////////////////////////////////////////
959fbac5 448void 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///////////////////////////////////////////////////////////////////////////
462void 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///////////////////////////////////////////////////////////////////////////
475void 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
387a745b 485 // Dummy invokation of GetScalar to obtain automatic proper error determination
486 e[0]=GetScalar();
487
959fbac5 488 e[0]=GetResultError();
1c01b4f8 489
959fbac5 490 for (Int_t i=0; i<3; i++)
491 {
492 e[i+1]=a[i];
493 }
494}
495///////////////////////////////////////////////////////////////////////////
496void 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///////////////////////////////////////////////////////////////////////////
84bb7c66 511void Ali4Vector::Data(TString f)
d88f97cc 512{
959fbac5 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
d88f97cc 519 if (f=="car" || f=="sph" || f=="cyl")
520 {
959fbac5 521 Double_t vec[4],err[4];
d88f97cc 522 GetVector(vec,f);
959fbac5 523 GetErrors(err,f);
524 Double_t inv=GetInvariant();
525 Double_t dinv=GetResultError();
c72198f1 526 cout << " Contravariant vector in " << f.Data() << " coordinates : "
d88f97cc 527 << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl;
c72198f1 528 cout << " ------------- Errors in " << f.Data() << " coordinates : "
959fbac5 529 << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl;
530 cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl;
d88f97cc 531 }
532 else
533 {
c72198f1 534 cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl
d88f97cc 535 << " Possible frames are 'car', 'sph' and 'cyl'." << endl;
536 }
537}
538///////////////////////////////////////////////////////////////////////////
539Double_t Ali4Vector::Dot(Ali4Vector& q)
540{
541// Provide the dot product of the current vector with vector q
959fbac5 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();
d88f97cc 557
959fbac5 558 Double_t dot=fV.Dot(b);
559 Double_t ddot=fV.GetResultError();
d88f97cc 560
959fbac5 561 dotpro=a0*b0-dot;
562
563 fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2));
d88f97cc 564 }
959fbac5 565
d88f97cc 566 return dotpro;
567}
568///////////////////////////////////////////////////////////////////////////
569Ali4Vector Ali4Vector::operator+(Ali4Vector& q)
570{
959fbac5 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();
d88f97cc 579
959fbac5 580 Double_t c0=a0+b0;
581 Ali3Vector c=a+b;
582 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
d88f97cc 583
584 Ali4Vector v;
959fbac5 585 v.SetVector(c0,c);
586 v.SetScalarError(dc0);
d88f97cc 587 return v;
588}
589///////////////////////////////////////////////////////////////////////////
590Ali4Vector Ali4Vector::operator-(Ali4Vector& q)
591{
959fbac5 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();
d88f97cc 600
959fbac5 601 Double_t c0=a0-b0;
602 Ali3Vector c=a-b;
603 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
d88f97cc 604
605 Ali4Vector v;
959fbac5 606 v.SetVector(c0,c);
607 v.SetScalarError(dc0);
d88f97cc 608 return v;
609}
610///////////////////////////////////////////////////////////////////////////
611Ali4Vector Ali4Vector::operator*(Double_t s)
612{
959fbac5 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();
d88f97cc 618
959fbac5 619 a0*=s;
620 a*=s;
621 da0*=s;
d88f97cc 622
623 Ali4Vector v;
959fbac5 624 v.SetVector(a0,a);
625 v.SetScalarError(da0);
d88f97cc 626
627 return v;
628}
629///////////////////////////////////////////////////////////////////////////
630Ali4Vector Ali4Vector::operator/(Double_t s)
631{
632// Divide the current vector by a scalar s
959fbac5 633// Error propagation is performed automatically
d88f97cc 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 {
959fbac5 642 Double_t a0=GetScalar();
643 Double_t da0=GetResultError();
644 Ali3Vector a=Get3Vector();
d88f97cc 645
959fbac5 646 a0/=s;
647 a/=s;
648 da0/=s;
d88f97cc 649
650 Ali4Vector v;
959fbac5 651 v.SetVector(a0,a);
652 v.SetScalarError(da0);
d88f97cc 653
654 return v;
655 }
656}
657///////////////////////////////////////////////////////////////////////////
658Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q)
659{
959fbac5 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();
d88f97cc 668
959fbac5 669 Double_t c0=a0+b0;
670 Ali3Vector c=a+b;
671 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
d88f97cc 672
959fbac5 673 SetVector(c0,c);
674 SetScalarError(dc0);
d88f97cc 675
676 return *this;
677}
678///////////////////////////////////////////////////////////////////////////
679Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q)
680{
959fbac5 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();
d88f97cc 689
959fbac5 690 Double_t c0=a0-b0;
691 Ali3Vector c=a-b;
692 Double_t dc0=sqrt(pow(da0,2)+pow(db0,2));
d88f97cc 693
959fbac5 694 SetVector(c0,c);
695 SetScalarError(dc0);
d88f97cc 696
d88f97cc 697 return *this;
698}
699///////////////////////////////////////////////////////////////////////////
700Ali4Vector& Ali4Vector::operator*=(Double_t s)
701{
959fbac5 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();
d88f97cc 707
959fbac5 708 a0*=s;
709 a*=s;
710 da0*=s;
d88f97cc 711
959fbac5 712 SetVector(a0,a);
713 SetScalarError(da0);
d88f97cc 714
d88f97cc 715 return *this;
716}
717///////////////////////////////////////////////////////////////////////////
718Ali4Vector& Ali4Vector::operator/=(Double_t s)
719{
720// Divide the current vector by a scalar s
959fbac5 721// Error propagation is performed automatically
d88f97cc 722
723 if (fabs(s)<1.e-20) // Protect against division by 0
724 {
959fbac5 725 cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl;
d88f97cc 726 return *this;
727 }
728 else
729 {
959fbac5 730 Double_t a0=GetScalar();
731 Double_t da0=GetResultError();
732 Ali3Vector a=Get3Vector();
d88f97cc 733
959fbac5 734 a0/=s;
735 a/=s;
736 da0/=s;
d88f97cc 737
959fbac5 738 SetVector(a0,a);
739 SetScalarError(da0);
d88f97cc 740
741 return *this;
742 }
743}
744///////////////////////////////////////////////////////////////////////////
959fbac5 745Int_t Ali4Vector::GetScalarFlag()
746{
747// Provide the value of the fScalar flag (for internal use only).
748 return fScalar;
749}
750///////////////////////////////////////////////////////////////////////////
d071d629 751Ali3Vector 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///////////////////////////////////////////////////////////////////////////
759Ali3Vector 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///////////////////////////////////////////////////////////////////////////
767Double_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///////////////////////////////////////////////////////////////////////////
790Double_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///////////////////////////////////////////////////////////////////////////
813Double_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///////////////////////////////////////////////////////////////////////////
176f88c0 823Ali3Vector 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///////////////////////////////////////////////////////////////////////////
844Double_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///////////////////////////////////////////////////////////////////////////
853Double_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///////////////////////////////////////////////////////////////////////////