]>
Commit | Line | Data |
---|---|---|
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 | // ------ | |
1f241680 | 88 | // Vectors (v), Errors (e), reference frames (f) and angular units (u) |
89 | // are specified via | |
90 | // SetVector(Float_t* v,TString f,TString u) | |
91 | // SetErrors(Float_t* e,TString f,TString u) | |
959fbac5 | 92 | // under the following conventions : |
93 | // | |
1f241680 | 94 | // f="car" ==> v in Cartesian coordinates (x,y,z) |
95 | // f="sph" ==> v in Spherical coordinates (r,theta,phi) | |
96 | // f="cyl" ==> v in Cylindrical coordinates (rho,phi,z) | |
959fbac5 | 97 | // |
1f241680 | 98 | // u="rad" ==> angles in radians |
99 | // u="deg" ==> angles in degrees | |
100 | // | |
101 | // The "f" and "u" facilities only serve as a convenient user interface. | |
102 | // Internally the actual storage of the various components is performed | |
103 | // in a unique way. This allows setting/retrieval of vector components in a | |
104 | // user selected frame/unit convention at any time. | |
959fbac5 | 105 | // |
106 | // Example : | |
107 | // --------- | |
108 | // | |
109 | // Ali4Vector a; | |
110 | // | |
111 | // Float_t v[4]={25,-1,3,7}; | |
112 | // a.SetVector(v,"car"); | |
113 | // | |
114 | // Float_t vec[4]; | |
1f241680 | 115 | // a.GetVector(vec,"sph","deg"); |
959fbac5 | 116 | // |
117 | // Ali4Vector b; | |
118 | // Float_t v2[4]={33,6,-18,2}; | |
119 | // b.SetVector(v2,"car"); | |
120 | // | |
121 | // Float_t dotpro=a.Dot(b); | |
122 | // | |
123 | // Float_t x0=16; | |
124 | // Ali3Vector x; | |
125 | // Float_t vec2[3]={1,2,3}; | |
126 | // x.SetVector(vec2,"car"); | |
127 | // | |
128 | // Ali4Vector c; | |
129 | // c.SetVector(x0,x); | |
130 | // c.GetVector(vec,"car"); | |
84bb7c66 | 131 | // c.Data("cyl"); |
959fbac5 | 132 | // c=a+b; |
133 | // c=a-b; | |
134 | // c=a*5; | |
135 | // | |
136 | //--- Author: Nick van Eijndhoven 01-apr-1999 UU-SAP Utrecht | |
f531a546 | 137 | //- Modified: NvE $Date$ UU-SAP Utrecht |
959fbac5 | 138 | /////////////////////////////////////////////////////////////////////////// |
139 | ||
d88f97cc | 140 | #include "Ali4Vector.h" |
c72198f1 | 141 | #include "Riostream.h" |
d88f97cc | 142 | |
143 | ClassImp(Ali4Vector) // Class implementation to enable ROOT I/O | |
144 | ||
145 | Ali4Vector::Ali4Vector() | |
146 | { | |
959fbac5 | 147 | // Creation of a contravariant 4-vector and initialisation of parameters. |
148 | // All values are initialised to 0. | |
149 | // Scalar mode is initially selected. | |
c72198f1 | 150 | SetZero(); |
959fbac5 | 151 | fScalar=1; |
d88f97cc | 152 | } |
153 | /////////////////////////////////////////////////////////////////////////// | |
154 | Ali4Vector::~Ali4Vector() | |
155 | { | |
156 | // Destructor to delete dynamically allocated memory | |
157 | } | |
158 | /////////////////////////////////////////////////////////////////////////// | |
c72198f1 | 159 | Ali4Vector::Ali4Vector(const Ali4Vector& v) |
160 | { | |
161 | // Copy constructor | |
162 | fScalar=v.fScalar; | |
163 | fV2=v.fV2; | |
164 | fDv2=v.fDv2; | |
165 | fV0=v.fV0; | |
166 | fDv0=v.fDv0; | |
167 | fDresult=v.fDresult; | |
168 | fV=v.fV; | |
169 | } | |
170 | /////////////////////////////////////////////////////////////////////////// | |
1fbffa23 | 171 | void Ali4Vector::Load(Ali4Vector& q) |
172 | { | |
173 | // Load all attributes of the input Ali4Vector into this Ali4Vector object. | |
174 | Int_t temp1=q.GetScalarFlag(); | |
175 | Double_t temp2=q.GetResultError(); | |
176 | Double_t a[4]; | |
177 | q.GetVector(a,"sph"); | |
178 | SetVector(a,"sph"); | |
179 | q.GetErrors(a,"car"); | |
180 | SetErrors(a,"car"); | |
181 | fScalar=temp1; | |
182 | fDresult=temp2; | |
183 | } | |
184 | /////////////////////////////////////////////////////////////////////////// | |
c72198f1 | 185 | void Ali4Vector::SetZero() |
186 | { | |
187 | // (Re)set all attributes to zero. | |
188 | // Note : The (de)selection of the scalar mode is not modified. | |
189 | fV2=0; | |
190 | fDv2=0; | |
191 | fV0=0; | |
192 | fDv0=0; | |
193 | fDresult=0; | |
194 | fV.SetZero(); | |
195 | } | |
196 | /////////////////////////////////////////////////////////////////////////// | |
197 | void Ali4Vector::SetVector(Double_t v0,Ali3Vector& v) | |
d88f97cc | 198 | { |
959fbac5 | 199 | // Store contravariant vector. |
200 | // The error on the scalar part is initialised to 0. | |
201 | // The errors on the vector part are taken from the input Ali3Vector. | |
202 | // Scalar mode is automatically selected. | |
203 | // The error on scalar result operations is reset to 0. | |
204 | fScalar=1; | |
d88f97cc | 205 | fV0=v0; |
206 | fV=v; | |
959fbac5 | 207 | fV2=pow(v0,2)-fV.Dot(fV); |
208 | SetScalarError(0); | |
d88f97cc | 209 | } |
210 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 211 | void Ali4Vector::SetVector(Double_t* v,TString f,TString u) |
d88f97cc | 212 | { |
959fbac5 | 213 | // Store vector according to reference frame f. |
1f241680 | 214 | // |
215 | // The string argument "u" allows to choose between different angular units | |
216 | // in case e.g. a spherical frame is selected. | |
217 | // u = "rad" : angles provided in radians | |
218 | // "deg" : angles provided in degrees | |
219 | // | |
220 | // The default is u="rad". | |
221 | // | |
959fbac5 | 222 | // All errors are initialised to 0. |
1f241680 | 223 | // |
959fbac5 | 224 | // Scalar mode is automatically selected. |
225 | // The error on scalar result operations is reset to 0. | |
1f241680 | 226 | |
959fbac5 | 227 | fScalar=1; |
d88f97cc | 228 | Double_t a[3]; |
229 | for (Int_t i=0; i<3; i++) | |
230 | { | |
231 | a[i]=v[i+1]; | |
232 | } | |
959fbac5 | 233 | fV0=v[0]; |
1f241680 | 234 | fV.SetVector(a,f,u); |
959fbac5 | 235 | fV2=pow(fV0,2)-fV.Dot(fV); |
236 | fDv2=0; | |
237 | fDv0=0; | |
238 | fDresult=0; | |
d88f97cc | 239 | } |
240 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 241 | void Ali4Vector::GetVector(Double_t* v,TString f,TString u) |
d88f97cc | 242 | { |
959fbac5 | 243 | // Provide 4-vector components according to reference frame f |
244 | // and according to the current mode. | |
245 | // Scalar mode : The scalar part is directly returned via v[0]. | |
246 | // Invariant mode : The scalar part is re-calculated via the value | |
247 | // of the Lorentz invariant and then returned via v[0]. | |
1f241680 | 248 | // |
249 | // The string argument "u" allows to choose between different angular units | |
250 | // in case e.g. a spherical frame is selected. | |
251 | // u = "rad" : angles provided in radians | |
252 | // "deg" : angles provided in degrees | |
253 | // | |
254 | // The default is u="rad". | |
255 | ||
959fbac5 | 256 | if (fScalar) |
257 | { | |
258 | v[0]=fV0; | |
259 | } | |
260 | else | |
261 | { | |
262 | v[0]=sqrt(fV.Dot(fV)+fV2); | |
263 | } | |
d88f97cc | 264 | Double_t a[3]; |
1f241680 | 265 | fV.GetVector(a,f,u); |
d88f97cc | 266 | for (Int_t i=0; i<3; i++) |
267 | { | |
268 | v[i+1]=a[i]; | |
269 | } | |
270 | } | |
271 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 272 | void Ali4Vector::SetVector(Float_t* v,TString f,TString u) |
d88f97cc | 273 | { |
959fbac5 | 274 | // Store vector according to reference frame f. |
1f241680 | 275 | // |
276 | // The string argument "u" allows to choose between different angular units | |
277 | // in case e.g. a spherical frame is selected. | |
278 | // u = "rad" : angles provided in radians | |
279 | // "deg" : angles provided in degrees | |
280 | // | |
281 | // The default is u="rad". | |
282 | // | |
959fbac5 | 283 | // All errors are initialised to 0. |
1f241680 | 284 | // |
959fbac5 | 285 | // Scalar mode is automatically selected. |
286 | // The error on scalar result operations is reset to 0. | |
1f241680 | 287 | |
d88f97cc | 288 | Double_t vec[4]; |
289 | for (Int_t i=0; i<4; i++) | |
290 | { | |
291 | vec[i]=v[i]; | |
292 | } | |
1f241680 | 293 | SetVector(vec,f,u); |
d88f97cc | 294 | } |
295 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 296 | void Ali4Vector::GetVector(Float_t* v,TString f,TString u) |
d88f97cc | 297 | { |
959fbac5 | 298 | // Provide 4-vector components according to reference frame f |
299 | // and according to the current mode. | |
300 | // Scalar mode : The scalar part is directly returned via v[0]. | |
301 | // Invariant mode : The scalar part is re-calculated via the value | |
302 | // of the Lorentz invariant and then returned via v[0]. | |
1f241680 | 303 | // |
304 | // The string argument "u" allows to choose between different angular units | |
305 | // in case e.g. a spherical frame is selected. | |
306 | // u = "rad" : angles provided in radians | |
307 | // "deg" : angles provided in degrees | |
308 | // | |
309 | // The default is u="rad". | |
310 | ||
d88f97cc | 311 | Double_t vec[4]; |
1f241680 | 312 | GetVector(vec,f,u); |
d88f97cc | 313 | for (Int_t i=0; i<4; i++) |
314 | { | |
315 | v[i]=vec[i]; | |
316 | } | |
317 | } | |
318 | /////////////////////////////////////////////////////////////////////////// | |
319 | Double_t Ali4Vector::GetScalar() | |
320 | { | |
959fbac5 | 321 | // Provide the scalar part. |
322 | // The error on the scalar value is available via GetResultError() | |
323 | // after invokation of GetScalar(). | |
324 | if (fScalar) | |
325 | { | |
326 | fDresult=fDv0; | |
327 | return fV0; | |
328 | } | |
329 | else | |
330 | { | |
331 | Double_t dot=fV.Dot(fV); | |
332 | Double_t ddot=fV.GetResultError(); | |
333 | Double_t v02=dot+fV2; | |
334 | Double_t dv02=sqrt(pow(ddot,2)+pow(fDv2,2)); | |
335 | Double_t v0=sqrt(fabs(v02)); | |
336 | Double_t dv0=0; | |
337 | if (v0) dv0=dv02/(2.*v0); | |
338 | fDresult=dv0; | |
339 | return v0; | |
340 | } | |
341 | } | |
342 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 343 | Double_t Ali4Vector::GetResultError() const |
959fbac5 | 344 | { |
345 | // Provide the error on the result of an operation yielding a scalar | |
346 | // E.g. GetScalar(), GetInvariant() or Dot() | |
347 | return fDresult; | |
348 | } | |
349 | /////////////////////////////////////////////////////////////////////////// | |
350 | void Ali4Vector::SetScalar(Double_t v0,Double_t dv0) | |
351 | { | |
352 | // Modify the scalar part (v0) and its error (dv0). | |
353 | // The default value for dv0 is 0. | |
354 | // The vector part is not modified. | |
355 | // Scalar mode is automatically selected | |
356 | // ==> Lorentz invariant and its error are updated. | |
357 | // The error on scalar result operations is reset to 0. | |
358 | fScalar=1; | |
359 | fV0=v0; | |
360 | fV2=pow(v0,2)-fV.Dot(fV); | |
361 | SetScalarError(dv0); | |
362 | } | |
363 | /////////////////////////////////////////////////////////////////////////// | |
364 | void Ali4Vector::SetScalarError(Double_t dv0) | |
365 | { | |
366 | // Set the error on the scalar part. | |
367 | // If in scalar mode, update error on the invariant accordingly. | |
368 | // The error on scalar result operations is reset to 0. | |
369 | fDv0=dv0; | |
370 | if (fScalar) | |
371 | { | |
372 | Double_t norm=fV.GetNorm(); | |
373 | Double_t dnorm=fV.GetResultError(); | |
374 | fDv2=sqrt(pow(2.*fV0*fDv0,2)+pow(2.*norm*dnorm,2)); | |
375 | } | |
376 | fDresult=0; | |
377 | } | |
378 | /////////////////////////////////////////////////////////////////////////// | |
c72198f1 | 379 | void Ali4Vector::Set3Vector(Ali3Vector& v) |
959fbac5 | 380 | { |
381 | // Set the 3-vector part, the errors are taken from the input Ali3Vector | |
382 | // Scalar mode : The scalar part and its error are not modified, | |
383 | // the Lorentz invariantand its error are re-calculated. | |
384 | // Invariant mode : The Lorentz invariant and its error are not modified, | |
385 | // the scalar part and its error are re-calculated. | |
386 | // The error on scalar result operations is reset to 0. | |
387 | fV=v; | |
388 | if (fScalar) | |
389 | { | |
390 | SetScalar(fV0,fDv0); | |
391 | } | |
392 | else | |
393 | { | |
394 | SetInvariant(fV2,fDv2); | |
395 | } | |
396 | } | |
397 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 398 | void Ali4Vector::Set3Vector(Double_t* v,TString f,TString u) |
959fbac5 | 399 | { |
400 | // Set the 3-vector part according to reference frame f | |
1f241680 | 401 | // |
402 | // The string argument "u" allows to choose between different angular units | |
403 | // in case e.g. a spherical frame is selected. | |
404 | // u = "rad" : angles provided in radians | |
405 | // "deg" : angles provided in degrees | |
406 | // | |
407 | // The default is u="rad". | |
408 | // | |
959fbac5 | 409 | // The errors on the vector part are initialised to 0 |
1f241680 | 410 | // |
959fbac5 | 411 | // Scalar mode : The scalar part and its error are not modified, |
412 | // the Lorentz invariantand its error are re-calculated. | |
413 | // Invariant mode : The Lorentz invariant and its error are not modified, | |
414 | // the scalar part and its error are re-calculated. | |
415 | // The error on scalar result operations is reset to 0. | |
1f241680 | 416 | |
959fbac5 | 417 | Double_t a[3]; |
418 | for (Int_t i=0; i<3; i++) | |
419 | { | |
420 | a[i]=v[i]; | |
421 | } | |
1f241680 | 422 | fV.SetVector(a,f,u); |
959fbac5 | 423 | |
424 | if (fScalar) | |
425 | { | |
426 | SetScalar(fV0,fDv0); | |
427 | } | |
428 | else | |
429 | { | |
430 | SetInvariant(fV2,fDv2); | |
431 | } | |
432 | } | |
433 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 434 | void Ali4Vector::Set3Vector(Float_t* v,TString f,TString u) |
959fbac5 | 435 | { |
436 | // Set the 3-vector part according to reference frame f | |
1f241680 | 437 | // |
438 | // The string argument "u" allows to choose between different angular units | |
439 | // in case e.g. a spherical frame is selected. | |
440 | // u = "rad" : angles provided in radians | |
441 | // "deg" : angles provided in degrees | |
442 | // | |
443 | // The default is u="rad". | |
444 | // | |
959fbac5 | 445 | // The errors on the vector part are initialised to 0 |
446 | // The Lorentz invariant is not modified | |
447 | // The error on scalar result operations is reset to 0. | |
1f241680 | 448 | |
959fbac5 | 449 | Double_t vec[3]; |
450 | for (Int_t i=0; i<3; i++) | |
451 | { | |
452 | vec[i]=v[i]; | |
453 | } | |
1f241680 | 454 | Set3Vector(vec,f,u); |
959fbac5 | 455 | } |
456 | /////////////////////////////////////////////////////////////////////////// | |
457 | void Ali4Vector::SetInvariant(Double_t v2,Double_t dv2) | |
458 | { | |
459 | // Modify the Lorentz invariant (v2) quantity v^i*v_i and its error (dv2). | |
460 | // The default value for the error dv2 is 0. | |
461 | // The vector part is not modified. | |
462 | // Invariant mode is automatically selected | |
463 | // ==> the scalar part and its error are updated. | |
464 | // The error on scalar result operations is reset to 0. | |
465 | // | |
466 | fScalar=0; | |
467 | fV2=v2; | |
468 | fDv2=dv2; | |
469 | fV0=GetScalar(); | |
470 | fDv0=GetResultError(); | |
471 | fDresult=0; | |
472 | } | |
473 | /////////////////////////////////////////////////////////////////////////// | |
474 | void Ali4Vector::SetInvariantError(Double_t dv2) | |
475 | { | |
476 | // Set the error on the Lorentz invariant. | |
477 | // If in invariant mode, update error on the scalar part accordingly. | |
478 | // The error on scalar result operations is reset to 0. | |
479 | fDv2=dv2; | |
480 | if (!fScalar) | |
481 | { | |
482 | fDv0=GetResultError(); | |
483 | } | |
484 | fDresult=0; | |
485 | } | |
486 | /////////////////////////////////////////////////////////////////////////// | |
487 | Double_t Ali4Vector::GetInvariant() | |
488 | { | |
489 | // Provide the Lorentz invariant v^i*v_i. | |
490 | // The error on the Lorentz invariant is available via GetResultError() | |
491 | // after invokation of GetInvariant(). | |
492 | if (!fScalar) | |
493 | { | |
494 | fDresult=fDv2; | |
495 | return fV2; | |
496 | } | |
497 | else | |
498 | { | |
499 | Double_t inv=Dot(*this); | |
500 | return inv; | |
501 | } | |
d88f97cc | 502 | } |
503 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 504 | Ali3Vector Ali4Vector::Get3Vector() const |
d88f97cc | 505 | { |
506 | // Provide the 3-vector part | |
507 | return fV; | |
508 | } | |
509 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 510 | void Ali4Vector::SetErrors(Double_t* e,TString f,TString u) |
959fbac5 | 511 | { |
512 | // Store errors for vector v^i according to reference frame f | |
1f241680 | 513 | // |
514 | // The string argument "u" allows to choose between different angular units | |
515 | // in case e.g. a spherical frame is selected. | |
516 | // u = "rad" : angles provided in radians | |
517 | // "deg" : angles provided in degrees | |
518 | // | |
519 | // The default is u="rad". | |
520 | // | |
959fbac5 | 521 | // If in scalar mode, update error on the invariant accordingly. |
522 | // The error on scalar result operations is reset to 0. | |
1f241680 | 523 | |
959fbac5 | 524 | Double_t a[3]; |
525 | for (Int_t i=0; i<3; i++) | |
526 | { | |
527 | a[i]=e[i+1]; | |
528 | } | |
529 | SetScalarError(e[0]); | |
1f241680 | 530 | fV.SetErrors(a,f,u); |
959fbac5 | 531 | } |
532 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 533 | void Ali4Vector::SetErrors(Float_t* e,TString f,TString u) |
959fbac5 | 534 | { |
535 | // Store errors for vector v^i according to reference frame f | |
1f241680 | 536 | // |
537 | // The string argument "u" allows to choose between different angular units | |
538 | // in case e.g. a spherical frame is selected. | |
539 | // u = "rad" : angles provided in radians | |
540 | // "deg" : angles provided in degrees | |
541 | // | |
542 | // The default is u="rad". | |
543 | // | |
959fbac5 | 544 | // If in scalar mode, update error on the invariant accordingly. |
545 | // The error on scalar result operations is reset to 0. | |
1f241680 | 546 | |
959fbac5 | 547 | Double_t a[4]; |
548 | for (Int_t i=0; i<4; i++) | |
549 | { | |
550 | a[i]=e[i]; | |
551 | } | |
1f241680 | 552 | SetErrors(a,f,u); |
959fbac5 | 553 | } |
554 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 555 | void Ali4Vector::GetErrors(Double_t* e,TString f,TString u) |
959fbac5 | 556 | { |
557 | // Provide errors for vector v^i according to reference frame f | |
558 | // and according to the current mode. | |
559 | // Scalar mode : The error on the scalar part is directly returned via e[0]. | |
560 | // Invariant mode : The error on the scalar part is re-calculated via the error | |
561 | // value on the Lorentz invariant and then returned via e[0]. | |
1f241680 | 562 | // |
563 | // The string argument "u" allows to choose between different angular units | |
564 | // in case e.g. a spherical frame is selected. | |
565 | // u = "rad" : angles provided in radians | |
566 | // "deg" : angles provided in degrees | |
567 | // | |
568 | // The default is u="rad". | |
569 | ||
959fbac5 | 570 | Double_t a[3]; |
1f241680 | 571 | fV.GetErrors(a,f,u); |
959fbac5 | 572 | |
387a745b | 573 | // Dummy invokation of GetScalar to obtain automatic proper error determination |
574 | e[0]=GetScalar(); | |
575 | ||
959fbac5 | 576 | e[0]=GetResultError(); |
1c01b4f8 | 577 | |
959fbac5 | 578 | for (Int_t i=0; i<3; i++) |
579 | { | |
580 | e[i+1]=a[i]; | |
581 | } | |
582 | } | |
583 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 584 | void Ali4Vector::GetErrors(Float_t* e,TString f,TString u) |
959fbac5 | 585 | { |
586 | // Provide errors for vector v^i according to reference frame f | |
587 | // and according to the current mode. | |
588 | // Scalar mode : The error on the scalar part is directly returned via e[0]. | |
589 | // Invariant mode : The error on the scalar part is re-calculated via the error | |
590 | // value on the Lorentz invariant and then returned via e[0]. | |
1f241680 | 591 | // |
592 | // The string argument "u" allows to choose between different angular units | |
593 | // in case e.g. a spherical frame is selected. | |
594 | // u = "rad" : angles provided in radians | |
595 | // "deg" : angles provided in degrees | |
596 | // | |
597 | // The default is u="rad". | |
598 | ||
959fbac5 | 599 | Double_t a[4]; |
1f241680 | 600 | GetErrors(a,f,u); |
959fbac5 | 601 | for (Int_t i=0; i<4; i++) |
602 | { | |
603 | e[i]=a[i]; | |
604 | } | |
605 | } | |
606 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 607 | void Ali4Vector::Data(TString f,TString u) |
d88f97cc | 608 | { |
959fbac5 | 609 | // Print contravariant vector components and errors according to |
610 | // reference frame f and according to the current mode. | |
611 | // Scalar mode : The scalar part and its error are directly returned. | |
612 | // Invariant mode : The scalar part and its error are re-calculated via the | |
613 | // value (and error) of the Lorentz invariant. | |
1f241680 | 614 | // |
615 | // The string argument "u" allows to choose between different angular units | |
616 | // in case e.g. a spherical frame is selected. | |
617 | // u = "rad" : angles provided in radians | |
618 | // "deg" : angles provided in degrees | |
619 | // | |
620 | // The defaults are f="car" and u="rad". | |
959fbac5 | 621 | |
d88f97cc | 622 | if (f=="car" || f=="sph" || f=="cyl") |
623 | { | |
959fbac5 | 624 | Double_t vec[4],err[4]; |
1f241680 | 625 | GetVector(vec,f,u); |
626 | GetErrors(err,f,u); | |
959fbac5 | 627 | Double_t inv=GetInvariant(); |
628 | Double_t dinv=GetResultError(); | |
1f241680 | 629 | cout << " Contravariant vector in " << f.Data() << " (" << u.Data() << ") coordinates : " |
d88f97cc | 630 | << vec[0] << " " << vec[1] << " " << vec[2] << " " << vec[3] << endl; |
1f241680 | 631 | cout << " ------------- Errors in " << f.Data() << " (" << u.Data() << ") coordinates : " |
959fbac5 | 632 | << err[0] << " " << err[1] << " " << err[2] << " " << err[3] << endl; |
633 | cout << " --- Lorentz invariant (v^i*v_i) : " << inv << " error : " << dinv << endl; | |
d88f97cc | 634 | } |
635 | else | |
636 | { | |
c72198f1 | 637 | cout << " *Ali4Vector::Data* Unsupported frame : " << f.Data() << endl |
d88f97cc | 638 | << " Possible frames are 'car', 'sph' and 'cyl'." << endl; |
639 | } | |
640 | } | |
641 | /////////////////////////////////////////////////////////////////////////// | |
642 | Double_t Ali4Vector::Dot(Ali4Vector& q) | |
643 | { | |
644 | // Provide the dot product of the current vector with vector q | |
959fbac5 | 645 | Double_t dotpro=0; |
646 | Double_t a0=GetScalar(); | |
647 | Double_t da0=GetResultError(); | |
648 | if ((this) == &q) // Check for special case v.Dot(v) | |
649 | { | |
650 | Double_t norm=fV.GetNorm(); | |
651 | Double_t dnorm=fV.GetResultError(); | |
652 | dotpro=pow(a0,2)-pow(norm,2); | |
653 | fDresult=sqrt(pow(2.*a0*da0,2)+pow(2.*norm*dnorm,2)); | |
654 | } | |
655 | else | |
656 | { | |
657 | Double_t b0=q.GetScalar(); | |
658 | Double_t db0=q.GetResultError(); | |
659 | Ali3Vector b=q.Get3Vector(); | |
d88f97cc | 660 | |
959fbac5 | 661 | Double_t dot=fV.Dot(b); |
662 | Double_t ddot=fV.GetResultError(); | |
d88f97cc | 663 | |
959fbac5 | 664 | dotpro=a0*b0-dot; |
665 | ||
666 | fDresult=sqrt(pow(b0*da0,2)+pow(a0*db0,2)+pow(ddot,2)); | |
d88f97cc | 667 | } |
959fbac5 | 668 | |
d88f97cc | 669 | return dotpro; |
670 | } | |
671 | /////////////////////////////////////////////////////////////////////////// | |
672 | Ali4Vector Ali4Vector::operator+(Ali4Vector& q) | |
673 | { | |
959fbac5 | 674 | // Add 4-vector q to the current 4-vector |
675 | // Error propagation is performed automatically | |
676 | Double_t a0=GetScalar(); | |
677 | Double_t da0=GetResultError(); | |
678 | Ali3Vector a=Get3Vector(); | |
679 | Double_t b0=q.GetScalar(); | |
680 | Double_t db0=q.GetResultError(); | |
681 | Ali3Vector b=q.Get3Vector(); | |
d88f97cc | 682 | |
959fbac5 | 683 | Double_t c0=a0+b0; |
684 | Ali3Vector c=a+b; | |
685 | Double_t dc0=sqrt(pow(da0,2)+pow(db0,2)); | |
d88f97cc | 686 | |
687 | Ali4Vector v; | |
959fbac5 | 688 | v.SetVector(c0,c); |
689 | v.SetScalarError(dc0); | |
d88f97cc | 690 | return v; |
691 | } | |
692 | /////////////////////////////////////////////////////////////////////////// | |
693 | Ali4Vector Ali4Vector::operator-(Ali4Vector& q) | |
694 | { | |
959fbac5 | 695 | // Subtract 4-vector q from the current 4-vector |
696 | // Error propagation is performed automatically | |
697 | Double_t a0=GetScalar(); | |
698 | Double_t da0=GetResultError(); | |
699 | Ali3Vector a=Get3Vector(); | |
700 | Double_t b0=q.GetScalar(); | |
701 | Double_t db0=q.GetResultError(); | |
702 | Ali3Vector b=q.Get3Vector(); | |
d88f97cc | 703 | |
959fbac5 | 704 | Double_t c0=a0-b0; |
705 | Ali3Vector c=a-b; | |
706 | Double_t dc0=sqrt(pow(da0,2)+pow(db0,2)); | |
d88f97cc | 707 | |
708 | Ali4Vector v; | |
959fbac5 | 709 | v.SetVector(c0,c); |
710 | v.SetScalarError(dc0); | |
d88f97cc | 711 | return v; |
712 | } | |
713 | /////////////////////////////////////////////////////////////////////////// | |
714 | Ali4Vector Ali4Vector::operator*(Double_t s) | |
715 | { | |
959fbac5 | 716 | // Multiply the current 4-vector with a scalar s |
717 | // Error propagation is performed automatically | |
718 | Double_t a0=GetScalar(); | |
719 | Double_t da0=GetResultError(); | |
720 | Ali3Vector a=Get3Vector(); | |
d88f97cc | 721 | |
959fbac5 | 722 | a0*=s; |
723 | a*=s; | |
724 | da0*=s; | |
d88f97cc | 725 | |
726 | Ali4Vector v; | |
959fbac5 | 727 | v.SetVector(a0,a); |
728 | v.SetScalarError(da0); | |
d88f97cc | 729 | |
730 | return v; | |
731 | } | |
732 | /////////////////////////////////////////////////////////////////////////// | |
733 | Ali4Vector Ali4Vector::operator/(Double_t s) | |
734 | { | |
735 | // Divide the current vector by a scalar s | |
959fbac5 | 736 | // Error propagation is performed automatically |
d88f97cc | 737 | |
738 | if (fabs(s)<1.e-20) // Protect against division by 0 | |
739 | { | |
740 | cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl; | |
741 | return *this; | |
742 | } | |
743 | else | |
744 | { | |
959fbac5 | 745 | Double_t a0=GetScalar(); |
746 | Double_t da0=GetResultError(); | |
747 | Ali3Vector a=Get3Vector(); | |
d88f97cc | 748 | |
959fbac5 | 749 | a0/=s; |
750 | a/=s; | |
751 | da0/=s; | |
d88f97cc | 752 | |
753 | Ali4Vector v; | |
959fbac5 | 754 | v.SetVector(a0,a); |
755 | v.SetScalarError(da0); | |
d88f97cc | 756 | |
757 | return v; | |
758 | } | |
759 | } | |
760 | /////////////////////////////////////////////////////////////////////////// | |
761 | Ali4Vector& Ali4Vector::operator+=(Ali4Vector& q) | |
762 | { | |
959fbac5 | 763 | // Add 4-vector q to the current 4-vector |
764 | // Error propagation is performed automatically | |
765 | Double_t a0=GetScalar(); | |
766 | Double_t da0=GetResultError(); | |
767 | Ali3Vector a=Get3Vector(); | |
768 | Double_t b0=q.GetScalar(); | |
769 | Double_t db0=q.GetResultError(); | |
770 | Ali3Vector b=q.Get3Vector(); | |
d88f97cc | 771 | |
959fbac5 | 772 | Double_t c0=a0+b0; |
773 | Ali3Vector c=a+b; | |
774 | Double_t dc0=sqrt(pow(da0,2)+pow(db0,2)); | |
d88f97cc | 775 | |
959fbac5 | 776 | SetVector(c0,c); |
777 | SetScalarError(dc0); | |
d88f97cc | 778 | |
779 | return *this; | |
780 | } | |
781 | /////////////////////////////////////////////////////////////////////////// | |
782 | Ali4Vector& Ali4Vector::operator-=(Ali4Vector& q) | |
783 | { | |
959fbac5 | 784 | // Subtract 4-vector q from the current 4-vector |
785 | // Error propagation is performed automatically | |
786 | Double_t a0=GetScalar(); | |
787 | Double_t da0=GetResultError(); | |
788 | Ali3Vector a=Get3Vector(); | |
789 | Double_t b0=q.GetScalar(); | |
790 | Double_t db0=q.GetResultError(); | |
791 | Ali3Vector b=q.Get3Vector(); | |
d88f97cc | 792 | |
959fbac5 | 793 | Double_t c0=a0-b0; |
794 | Ali3Vector c=a-b; | |
795 | Double_t dc0=sqrt(pow(da0,2)+pow(db0,2)); | |
d88f97cc | 796 | |
959fbac5 | 797 | SetVector(c0,c); |
798 | SetScalarError(dc0); | |
d88f97cc | 799 | |
d88f97cc | 800 | return *this; |
801 | } | |
802 | /////////////////////////////////////////////////////////////////////////// | |
803 | Ali4Vector& Ali4Vector::operator*=(Double_t s) | |
804 | { | |
959fbac5 | 805 | // Multiply the current 4-vector with a scalar s |
806 | // Error propagation is performed automatically | |
807 | Double_t a0=GetScalar(); | |
808 | Double_t da0=GetResultError(); | |
809 | Ali3Vector a=Get3Vector(); | |
d88f97cc | 810 | |
959fbac5 | 811 | a0*=s; |
812 | a*=s; | |
813 | da0*=s; | |
d88f97cc | 814 | |
959fbac5 | 815 | SetVector(a0,a); |
816 | SetScalarError(da0); | |
d88f97cc | 817 | |
d88f97cc | 818 | return *this; |
819 | } | |
820 | /////////////////////////////////////////////////////////////////////////// | |
821 | Ali4Vector& Ali4Vector::operator/=(Double_t s) | |
822 | { | |
823 | // Divide the current vector by a scalar s | |
959fbac5 | 824 | // Error propagation is performed automatically |
d88f97cc | 825 | |
826 | if (fabs(s)<1.e-20) // Protect against division by 0 | |
827 | { | |
959fbac5 | 828 | cout << " *Ali4Vector::/* Division by 0 detected. No action taken." << endl; |
d88f97cc | 829 | return *this; |
830 | } | |
831 | else | |
832 | { | |
959fbac5 | 833 | Double_t a0=GetScalar(); |
834 | Double_t da0=GetResultError(); | |
835 | Ali3Vector a=Get3Vector(); | |
d88f97cc | 836 | |
959fbac5 | 837 | a0/=s; |
838 | a/=s; | |
839 | da0/=s; | |
d88f97cc | 840 | |
959fbac5 | 841 | SetVector(a0,a); |
842 | SetScalarError(da0); | |
d88f97cc | 843 | |
844 | return *this; | |
845 | } | |
846 | } | |
847 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 848 | Int_t Ali4Vector::GetScalarFlag() const |
959fbac5 | 849 | { |
850 | // Provide the value of the fScalar flag (for internal use only). | |
851 | return fScalar; | |
852 | } | |
853 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 854 | Ali3Vector Ali4Vector::GetVecTrans() const |
d071d629 | 855 | { |
856 | // Provide the transverse vector part w.r.t. z-axis. | |
857 | // Error propagation is performed automatically | |
858 | ||
859 | return fV.GetVecTrans(); | |
860 | } | |
861 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 862 | Ali3Vector Ali4Vector::GetVecLong() const |
d071d629 | 863 | { |
864 | // Provide the longitudinal vector part w.r.t. z-axis. | |
865 | // Error propagation is performed automatically | |
866 | ||
867 | return fV.GetVecLong(); | |
868 | } | |
869 | /////////////////////////////////////////////////////////////////////////// | |
870 | Double_t Ali4Vector::GetScaTrans() | |
871 | { | |
872 | // Provide the "transverse value" of the scalar part w.r.t. z-axis. | |
873 | // This provides a basis for e.g. E_trans calculation. | |
874 | // Note : the returned value is always positive or zero. | |
875 | // The error on the value is available via GetResultError() | |
876 | // after invokation of GetScaTrans(). | |
877 | Double_t a[3],ea[3]; | |
878 | ||
879 | fV.GetVector(a,"sph"); | |
880 | fV.GetErrors(ea,"sph"); | |
881 | ||
882 | Double_t s=GetScalar(); | |
883 | Double_t ds=GetResultError(); | |
884 | ||
885 | Double_t st,dst2; | |
886 | st=s*sin(a[1]); | |
887 | dst2=pow((sin(a[1])*ds),2)+pow((s*cos(a[1])*ea[1]),2); | |
888 | ||
889 | fDresult=sqrt(dst2); | |
890 | return fabs(st); | |
891 | } | |
892 | /////////////////////////////////////////////////////////////////////////// | |
893 | Double_t Ali4Vector::GetScaLong() | |
894 | { | |
895 | // Provide the "longitudinal value" of the scalar part w.r.t. z-axis. | |
896 | // This provides a basis for e.g. E_long calculation. | |
897 | // Note : the returned value can also be negative. | |
898 | // The error on the value is available via GetResultError() | |
899 | // after invokation of GetScaLong(). | |
900 | Double_t a[3],ea[3]; | |
901 | ||
902 | fV.GetVector(a,"sph"); | |
903 | fV.GetErrors(ea,"sph"); | |
904 | ||
905 | Double_t s=GetScalar(); | |
906 | Double_t ds=GetResultError(); | |
907 | ||
908 | Double_t sl,dsl2; | |
909 | sl=s*cos(a[1]); | |
910 | dsl2=pow((cos(a[1])*ds),2)+pow((s*sin(a[1])*ea[1]),2); | |
911 | ||
912 | fDresult=sqrt(dsl2); | |
913 | return sl; | |
914 | } | |
915 | /////////////////////////////////////////////////////////////////////////// | |
916 | Double_t Ali4Vector::GetPseudoRapidity() | |
917 | { | |
918 | // Provide the pseudorapidity value of the vector part w.r.t. z-axis. | |
919 | // The error on the value is available via GetResultError() | |
920 | // after invokation of GetPseudoRapidity(). | |
921 | Double_t eta=fV.GetPseudoRapidity(); | |
922 | fDresult=fV.GetResultError(); | |
923 | return eta; | |
924 | } | |
925 | /////////////////////////////////////////////////////////////////////////// | |
261c0caf | 926 | Ali3Vector Ali4Vector::GetBetaVector() const |
176f88c0 | 927 | { |
928 | // Provide the beta 3-vector corresponding to this 4-vector. | |
929 | Ali3Vector beta; | |
930 | if (fabs(fV0)>0.) beta=fV/fV0; | |
931 | if (fabs(fDv0)>0.) | |
932 | { | |
933 | Double_t vecv[3],errv[3],errb[3]; | |
934 | Double_t sqerr=0; | |
935 | fV.GetVector(vecv,"car"); | |
936 | fV.GetErrors(errv,"car"); | |
937 | for (Int_t i=0; i<3; i++) | |
938 | { | |
939 | sqerr=pow((errv[i]/fV0),2)+pow((vecv[i]*fDv0/(fV0*fV0)),2); | |
940 | errb[i]=sqrt(sqerr); | |
941 | } | |
942 | beta.SetErrors(errb,"car"); | |
943 | } | |
944 | return beta; | |
945 | } | |
946 | /////////////////////////////////////////////////////////////////////////// | |
947 | Double_t Ali4Vector::GetBeta() | |
948 | { | |
949 | // Provide the beta value (i.e. v/c) corresponding to this 4-vector. | |
950 | Ali3Vector beta=GetBetaVector(); | |
951 | Double_t val=beta.GetNorm(); | |
952 | fDresult=beta.GetResultError(); | |
953 | return val; | |
954 | } | |
955 | /////////////////////////////////////////////////////////////////////////// | |
956 | Double_t Ali4Vector::GetGamma() | |
957 | { | |
958 | // Provide the Lorentz gamma factor corresponding to this 4-vector. | |
959 | // In case the gamma factor is infinite a value of -1 is returned. | |
960 | Double_t gamma=-1; | |
961 | fDresult=0; | |
962 | Double_t inv=sqrt(fabs(fV2)); | |
963 | if (inv>0.) | |
964 | { | |
965 | Double_t dinv=fDv2/(2.*inv); | |
966 | gamma=fV0/inv; | |
967 | Double_t sqerr=pow((fDv0/inv),2)+pow((fV0*dinv/fV2),2); | |
968 | fDresult=sqrt(sqerr); | |
969 | } | |
970 | return gamma; | |
971 | } | |
972 | /////////////////////////////////////////////////////////////////////////// | |
4962c850 | 973 | Double_t Ali4Vector::GetX(Int_t i,TString f,TString u) |
974 | { | |
975 | // Provide i-th vector component according to reference frame f. | |
976 | // | |
977 | // The string argument "u" allows to choose between different angular units | |
978 | // in case e.g. a spherical frame is selected. | |
979 | // u = "rad" : angles provided in radians | |
980 | // "deg" : angles provided in degrees | |
981 | // | |
982 | // The default is u="rad". | |
983 | // | |
984 | // The vector components are addressed via the generic x0,x1,x2,x3 notation. | |
985 | // So, i=0 denotes the scalar component and i=1 denotes the first 3-vector component. | |
986 | // The error on the selected component can be obtained via the | |
987 | // usual GetResultError() facility. | |
988 | ||
989 | fDresult=0; | |
990 | ||
991 | if (i<0 || i>3) return 0; | |
992 | ||
993 | Double_t x=0; | |
994 | ||
995 | if (i==0) | |
996 | { | |
997 | x=GetScalar(); | |
998 | } | |
999 | else | |
1000 | { | |
1001 | x=fV.GetX(i,f,u); | |
1002 | fDresult=fV.GetResultError(); | |
1003 | } | |
1004 | ||
1005 | return x; | |
1006 | } | |
1007 | /////////////////////////////////////////////////////////////////////////// | |
1f241680 | 1008 | Double_t Ali4Vector::GetOpeningAngle(Ali4Vector& q,TString u) |
1009 | { | |
1010 | // Provide the opening angle between 3-vector parts with 4-vector q. | |
1011 | // The string argument "u" allows to choose between different output units. | |
1012 | // u = "rad" : opening angle provided in radians | |
1013 | // "deg" : opening angle provided in degrees | |
1014 | // | |
1015 | // The default is u="rad". | |
1016 | ||
1017 | Ali3Vector v1=fV; | |
1018 | Ali3Vector v2=q.Get3Vector(); | |
1019 | ||
1020 | Double_t ang=v1.GetOpeningAngle(v2,u); | |
1021 | fDresult=v1.GetResultError(); | |
1022 | ||
1023 | return ang; | |
1024 | } | |
1025 | /////////////////////////////////////////////////////////////////////////// | |
413d0114 | 1026 | Double_t Ali4Vector::GetOpeningAngle(Ali3Vector& q,TString u) |
1027 | { | |
1028 | // Provide the opening angle between the 3-vector part and 3-vector q. | |
1029 | // The string argument "u" allows to choose between different output units. | |
1030 | // u = "rad" : opening angle provided in radians | |
1031 | // "deg" : opening angle provided in degrees | |
1032 | // | |
1033 | // The default is u="rad". | |
1034 | ||
1035 | Ali3Vector v1=fV; | |
1036 | ||
1037 | Double_t ang=v1.GetOpeningAngle(q,u); | |
1038 | fDresult=v1.GetResultError(); | |
1039 | ||
1040 | return ang; | |
1041 | } | |
1042 | /////////////////////////////////////////////////////////////////////////// |