27-may-2001 NvE New class AliEvent introduced and RALICELinkDef.h & co. updated accor...
[u/mrichter/AliRoot.git] / RALICE / AliTrack.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 AliTrack
20// Handling of the attributes of a reconstructed particle track.
21//
22// Coding example :
23// ----------------
24//
25// Float_t a[4]={195.,1.2,-0.04,8.5};
26// Ali4Vector pmu;
27// pmu.SetVector(a,"car");
28// AliTrack t1;
29// t1.Set4Momentum(pmu);
30//
31// Float_t b[3]={1.2,-0.04,8.5};
32// Ali3Vector p;
33// p.SetVector(b,"car");
34// AliTrack t2;
35// t2.Set3Momentum(p);
36// t2.SetCharge(0);
37// t2.SetMass(1.115);
38//
39// t1.Info();
40// t2.Info();
41//
42// Float_t pi=acos(-1.);
43// Float_t thcms=0.2*pi; // decay theta angle in cms
44// Float_t phicms=pi/4.; // decay theta angle in cms
45// Float_t m1=0.938;
46// Float_t m2=0.140;
47// t2.Decay(m1,m2,thcms,phicms); // Track t2 decay : Lambda -> proton + pion
48//
49// t2.List();
50//
51// Int_t ndec=t2.GetNdecay();
52// AliTrack* d1=t2.GetDecayTrack(1); // Access to decay track number 1
53// AliTrack* d2=t2.GetDecayTrack(2); // Access to decay track number 2
54//
55// AliSignal s1,s2,s3,s4;
56//
57// .... // Code (e.g. detector readout) to fill AliSignal data
58//
59// AliTrack trec; // Track which will be reconstructed from signals
60// trec.AddSignal(s1);
61// trec.AddSignal(s3);
62// trec.AddSignal(s4);
63//
64// Ali3Vector P;
65// Float_t Q,M;
66//
67// ... // Code which accesses signals from trec and reconstructs
68// 3-momentum P, charge Q, mass M etc...
69//
70// trec.Set3Momentum(P);
71// trec.SetCharge(Q);
72// trec.SetMass(M);
73//
74// Float_t r1[3]={1.6,-3.8,25.7};
75// Float_t er1[3]={0.2,0.5,1.8};
76// Float_t r2[3]={8.6,23.8,-6.7};
77// Float_t er2[3]={0.93,1.78,0.8};
78// AliPosition begin,end;
79// begin.SetPosition(r1,"car");
80// begin.SetPositionErrors(er1,"car");
81// end.SetPosition(r2,"car");
82// end.SetPositionErrors(er2,"car");
83// trec.SetBeginPoint(begin);
84// trec.SetEndPoint(end);
85//
86// Note : All quantities are in GeV, GeV/c or GeV/c**2
87//
88//--- Author: Nick van Eijndhoven 10-jul-1997 UU-SAP Utrecht
f531a546 89//- Modified: NvE $Date$ UU-SAP Utrecht
959fbac5 90///////////////////////////////////////////////////////////////////////////
91
d88f97cc 92#include "AliTrack.h"
93
94ClassImp(AliTrack) // Class implementation to enable ROOT I/O
95
96AliTrack::AliTrack()
97{
98// Default constructor
99// All variables initialised to 0
100 fDecays=0;
959fbac5 101 fSignals=0;
f531a546 102 fMasses=0;
103 fDmasses=0;
104 fPmasses=0;
d88f97cc 105 Reset();
106}
107///////////////////////////////////////////////////////////////////////////
108AliTrack::~AliTrack()
109{
110// Destructor to delete memory allocated for decay tracks array
111 if (fDecays)
112 {
113 fDecays->Delete();
114 delete fDecays;
115 fDecays=0;
116 }
959fbac5 117 if (fSignals)
118 {
119 fSignals->Clear();
120 delete fSignals;
121 fSignals=0;
122 }
d88f97cc 123}
124///////////////////////////////////////////////////////////////////////////
125void AliTrack::Reset()
126{
127// Reset all variables to 0
d88f97cc 128 fQ=0;
129 fNdec=0;
959fbac5 130 fNsig=0;
f531a546 131 fNmasses=0;
d88f97cc 132 Double_t a[4]={0,0,0,0};
133 SetVector(a,"sph");
134 if (fDecays)
135 {
136 fDecays->Delete();
137 delete fDecays;
138 fDecays=0;
139 }
959fbac5 140 if (fSignals)
141 {
142 fSignals->Clear();
143 delete fSignals;
144 fSignals=0;
145 }
146 Double_t b[3]={0,0,0};
147 fBegin.SetPosition(b,"sph");
148 fEnd.SetPosition(b,"sph");
f531a546 149 if (fMasses)
150 {
151 delete fMasses;
152 fMasses=0;
153 }
154 if (fDmasses)
155 {
156 delete fDmasses;
157 fDmasses=0;
158 }
159 if (fPmasses)
160 {
161 delete fPmasses;
162 fPmasses=0;
163 }
d88f97cc 164}
165///////////////////////////////////////////////////////////////////////////
166void AliTrack::Set3Momentum(Ali3Vector& p)
167{
168// Set the track parameters according to the 3-momentum p
959fbac5 169 Set3Vector(p);
d88f97cc 170}
171///////////////////////////////////////////////////////////////////////////
172void AliTrack::Set4Momentum(Ali4Vector& p)
173{
174// Set the track parameters according to the 4-momentum p
175 Double_t E=p.GetScalar();
959fbac5 176 Double_t dE=p.GetResultError();
d88f97cc 177 Ali3Vector pv=p.Get3Vector();
178 SetVector(E,pv);
959fbac5 179 SetScalarError(dE);
d88f97cc 180}
181///////////////////////////////////////////////////////////////////////////
959fbac5 182void AliTrack::SetMass(Double_t m,Double_t dm)
d88f97cc 183{
184// Set the particle mass
959fbac5 185// The default value for the error dm is 0.
186 Double_t inv=pow(m,2);
187 Double_t dinv=fabs(2.*m*dm);
188 SetInvariant(inv,dinv);
d88f97cc 189}
190///////////////////////////////////////////////////////////////////////////
191void AliTrack::SetCharge(Float_t q)
192{
193// Set the particle charge
194 fQ=q;
195}
196///////////////////////////////////////////////////////////////////////////
197void AliTrack::Info(TString f)
198{
199// Provide track information within the coordinate frame f
959fbac5 200 Double_t m=GetMass();
201 Double_t dm=GetResultError();
202 cout << " *AliTrack::Info* Mass : " << m
203 << " error : " << dm << " Charge : " << fQ
f531a546 204 << " Momentum : " << GetMomentum() << " Nmass hyp. : " << fNmasses
205 << " Ntracks : " << fNdec << " Nsignals : " << fNsig << endl;
206 for (Int_t i=0; i<fNmasses; i++)
207 {
208 cout << " Mass hypothesis " << (i+1) << " Mass : " << fMasses->At(i)
209 << " error : " << fDmasses->At(i) << " prob. : " << fPmasses->At(i)
210 << endl;
211 }
d88f97cc 212 Ali4Vector::Info(f);
213}
214///////////////////////////////////////////////////////////////////////////
215void AliTrack::List(TString f)
216{
217// Provide current track and decay level 1 information within coordinate frame f
218
219 Info(f); // Information of the current track
220
221 // Decay products of this track
222 AliTrack* td;
223 for (Int_t id=1; id<=fNdec; id++)
224 {
225 td=GetDecayTrack(id);
226 if (td)
227 {
228 cout << " ---Level 1 sec. track no. " << id << endl;
d88f97cc 229 td->Info(f);
230 }
231 else
232 {
233 cout << " *AliTrack::List* Error : No decay track present." << endl;
234 }
235 }
236}
237///////////////////////////////////////////////////////////////////////////
238void AliTrack::ListAll(TString f)
239{
240// Provide complete track and decay information within the coordinate frame f
241
242 Info(f); // Information of the current track
959fbac5 243 cout << " Begin-point :"; fBegin.Info(f);
244 cout << " End-point :"; fEnd.Info(f);
245 for (Int_t is=1; is<=GetNsignals(); is++)
246 {
247 ((AliSignal*)GetSignal(is))->Info(f);
248 }
d88f97cc 249
250 AliTrack* t=this;
251 Dump(t,1,f); // Information of all decay products
252}
253//////////////////////////////////////////////////////////////////////////
254void AliTrack::Dump(AliTrack* t,Int_t n,TString f)
255{
256// Recursively provide the info of all decay levels of this track
257 AliTrack* td;
258 for (Int_t id=1; id<=t->GetNdecay(); id++)
259 {
260 td=t->GetDecayTrack(id);
261 if (td)
262 {
263 cout << " ---Level " << n << " sec. track no. " << id << endl;
d88f97cc 264 td->Info(f);
959fbac5 265 for (Int_t is=1; is<=td->GetNsignals(); is++)
266 {
267 ((AliSignal*)td->GetSignal(is))->Info(f);
268 }
d88f97cc 269
270 // Go for next decay level of this decay track recursively
271 Dump(td,n+1,f);
272 }
273 else
274 {
275 cout << " *AliTrack::Dump* Error : No decay track present." << endl;
276 }
277 }
278}
279//////////////////////////////////////////////////////////////////////////
280Double_t AliTrack::GetMomentum()
281{
959fbac5 282// Provide the value of the track 3-momentum.
283// The error can be obtained by invoking GetResultError() after
284// invokation of GetMomentum().
285
286// Ali3Vector p=Get3Vector();
287// return sqrt(p.Dot(p));
288 Double_t norm=fV.GetNorm();
289 return norm;
d88f97cc 290}
291///////////////////////////////////////////////////////////////////////////
292Ali3Vector AliTrack::Get3Momentum()
293{
294// Provide the track 3-momentum
295 return (Ali3Vector)Get3Vector();
296}
297///////////////////////////////////////////////////////////////////////////
298Double_t AliTrack::GetMass()
299{
959fbac5 300// Provide the particle mass.
301// The error can be obtained by invoking GetResultError() after
302// invokation of GetMass().
303 Double_t inv=GetInvariant();
304 Double_t dinv=GetResultError();
305 Double_t dm=0;
306 if (inv >= 0)
307 {
308 Double_t m=sqrt(inv);
309 if (m) dm=dinv/(2.*m);
310 fDresult=dm;
311 return m;
312 }
313 else
314 {
315 cout << "*AliTrack::GetMass* Unphysical situation m**2 = " << inv << endl;
316 cout << " Value 0 will be returned." << endl;
317 fDresult=dm;
318 return 0;
319 }
d88f97cc 320}
321///////////////////////////////////////////////////////////////////////////
322Float_t AliTrack::GetCharge()
323{
324// Provide the particle charge
325 return fQ;
326}
327///////////////////////////////////////////////////////////////////////////
328Double_t AliTrack::GetEnergy()
329{
959fbac5 330// Provide the particle's energy.
331// The error can be obtained by invoking GetResultError() after
332// invokation of GetEnergy().
333 Double_t E=GetScalar();
334 if (E>0)
335 {
336 return E;
337 }
338 else
339 {
340 cout << "*AliTrack::GetEnergy* Unphysical situation E = " << E << endl;
341 cout << " Value 0 will be returned." << endl;
342 return 0;
343 }
d88f97cc 344}
345///////////////////////////////////////////////////////////////////////////
346void AliTrack::Decay(Double_t m1,Double_t m2,Double_t thcms,Double_t phicms)
347{
348// Perform 2-body decay of current track
349// m1 : mass of decay product 1
350// m2 : mass of decay product 2
351// thcms : cms theta decay angle (in rad.) of m1
352// phicms : cms phi decay angle (in rad.) of m1
353
354 fNdec=2; // it's a 2-body decay
959fbac5 355
356 Double_t M=GetMass();
d88f97cc 357
358// Compute the 4-momenta of the decay products in the cms
359// Note : p2=p1=pnorm for a 2-body decay
959fbac5 360 Double_t e1=0;
361 if (M) e1=((M*M)+(m1*m1)-(m2*m2))/(2.*M);
362 Double_t e2=0;
363 if (M) e2=((M*M)+(m2*m2)-(m1*m1))/(2.*M);
d88f97cc 364 Double_t pnorm=(e1*e1)-(m1*m1);
365 if (pnorm>0.)
366 {
367 pnorm=sqrt(pnorm);
368 }
369 else
370 {
371 pnorm=0;
372 }
373
374 Double_t a[3];
375 a[0]=pnorm;
376 a[1]=thcms;
377 a[2]=phicms;
378 Ali3Vector p;
379 p.SetVector(a,"sph");
380
381 Ali4Vector pprim1;
382 pprim1.SetVector(e1,p);
959fbac5 383 pprim1.SetInvariant(m1*m1);
d88f97cc 384
385 Ali4Vector pprim2;
386 p*=-1;
387 pprim2.SetVector(e2,p);
959fbac5 388 pprim2.SetInvariant(m2*m2);
d88f97cc 389
390 // Determine boost parameters from the parent particle
959fbac5 391 Double_t E=GetEnergy();
d88f97cc 392 p=Get3Vector();
393 Ali4Vector pmu;
394 pmu.SetVector(E,p);
395
396 AliBoost q;
397 q.Set4Momentum(pmu);
398
399 Ali4Vector p1=q.Inverse(pprim1); // Boost decay product 1
400 Ali4Vector p2=q.Inverse(pprim2); // Boost decay product 2
401
402 // Enter the boosted data into the decay tracks array
403 if (fDecays)
404 {
405 fDecays->Delete();
406 delete fDecays;
407 }
408 fDecays=new TObjArray();
409
410 fDecays->Add(new AliTrack);
411 ((AliTrack*)fDecays->At(0))->Set4Momentum(p1);
959fbac5 412 ((AliTrack*)fDecays->At(0))->SetMass(m1);
d88f97cc 413 fDecays->Add(new AliTrack);
414 ((AliTrack*)fDecays->At(1))->Set4Momentum(p2);
d88f97cc 415 ((AliTrack*)fDecays->At(1))->SetMass(m2);
416}
417///////////////////////////////////////////////////////////////////////////
418Int_t AliTrack::GetNdecay()
419{
420// Provide the number of decay produced tracks
421 return fNdec;
422}
423///////////////////////////////////////////////////////////////////////////
424AliTrack* AliTrack::GetDecayTrack(Int_t j)
425{
426// Provide decay produced track number j
427// Note : j=1 denotes the first decay track
428 if ((j >= 1) && (j <= fNdec))
429 {
430 return (AliTrack*)fDecays->At(j-1);
431 }
432 else
433 {
434 cout << " *AliTrack* decay track number : " << j << " out of range." << endl;
435 cout << " -- Decay track number 1 (if any) returned." << endl;
436 return (AliTrack*)fDecays->At(0);
437 }
438}
439///////////////////////////////////////////////////////////////////////////
959fbac5 440void AliTrack::AddSignal(AliSignal& s)
441{
442// Relate an AliSignal object to this track.
443 if (!fSignals) fSignals=new TObjArray();
444 fNsig++;
445 fSignals->Add(&s);
446}
447///////////////////////////////////////////////////////////////////////////
448void AliTrack::RemoveSignal(AliSignal& s)
449{
450// Remove related AliSignal object to this track.
451 if (fSignals)
452 {
453 AliSignal* test=(AliSignal*)fSignals->Remove(&s);
454 if (test)
455 {
456 fNsig--;
457 fSignals->Compress();
458 }
459 }
460}
461///////////////////////////////////////////////////////////////////////////
462Int_t AliTrack::GetNsignals()
463{
464// Provide the number of related AliSignals.
465 return fNsig;
466}
467///////////////////////////////////////////////////////////////////////////
468AliSignal* AliTrack::GetSignal(Int_t j)
469{
470// Provide the related AliSignal number j.
471// Note : j=1 denotes the first signal.
472 if ((j >= 1) && (j <= fNsig))
473 {
474 return (AliSignal*)fSignals->At(j-1);
475 }
476 else
477 {
478 cout << " *AliTrack* signal number : " << j << " out of range." << endl;
479 cout << " -- Signal number 1 (if any) returned." << endl;
480 return (AliSignal*)fDecays->At(0);
481 }
482}
483///////////////////////////////////////////////////////////////////////////
484void AliTrack::SetBeginPoint(AliPosition p)
485{
486// Store the position of the track begin-point.
487 fBegin=p;
488}
489///////////////////////////////////////////////////////////////////////////
490AliPosition AliTrack::GetBeginPoint()
491{
492// Provide the position of the track begin-point.
493 return fBegin;
494}
495///////////////////////////////////////////////////////////////////////////
496void AliTrack::SetEndPoint(AliPosition p)
497{
498// Store the position of the track end-point.
499 fEnd=p;
500}
501///////////////////////////////////////////////////////////////////////////
502AliPosition AliTrack::GetEndPoint()
503{
504// Provide the position of the track end-point.
505 return fEnd;
506}
507///////////////////////////////////////////////////////////////////////////
f531a546 508void AliTrack::AddMassHypothesis(Double_t prob,Double_t m,Double_t dm)
509{
510// Add a mass hypothesis for this current track.
511// prob=probalility m=mass value dm=error on the mass value.
512// The default value for the mass error dm is 0.
513 if (!fMasses) fMasses=new TArrayD();
514 if (!fDmasses) fDmasses=new TArrayD();
515 if (!fPmasses) fPmasses=new TArrayD();
516
517 fNmasses++;
518 fMasses->Set(fNmasses);
519 fDmasses->Set(fNmasses);
520 fPmasses->Set(fNmasses);
521
522 fMasses->AddAt(m,fNmasses-1);
523 fDmasses->AddAt(dm,fNmasses-1);
524 fPmasses->AddAt(prob,fNmasses-1);
525}
526///////////////////////////////////////////////////////////////////////////
527Int_t AliTrack::GetNMassHypotheses()
528{
529// Provide the number of mass hypotheses for this track.
530 return fNmasses;
531}
532///////////////////////////////////////////////////////////////////////////
533Double_t AliTrack::GetMassHypothesis(Int_t j)
534{
535// Provide the mass of the jth hypothesis for this track.
536// Note : the first hypothesis is indicated by j=1.
537// Default : j=0 ==> Hypothesis with highest probability.
538// The error on the mass can be obtained by invoking GetResultError()
539// after invokation of GetMassHypothesis(j).
540
541 Double_t m=0,dm=0,prob=0;
542
543 // Check validity of index j
544 if (j<0 || j>fNmasses)
545 {
546 cout << " *AliTrack::GetMassHypothesis* Invalid index j : " << j
547 << " Number of mass hypotheses : " << fNmasses << endl;
548 fDresult=0;
549 return 0;
550 }
551
552 // Select mass hypothesis with highest probability
553 if (j==0)
554 {
555 if (fNmasses)
556 {
557 m=fMasses->At(0);
558 dm=fDmasses->At(0);
559 prob=fPmasses->At(0);
560 for (Int_t i=1; i<fNmasses; i++)
561 {
562 if (fPmasses->At(i)>prob)
563 {
564 m=fMasses->At(i);
565 dm=fDmasses->At(i);
566 }
567 }
568 }
569 fDresult=dm;
570 return m;
571 }
572
573 // Provide data of requested mass hypothesis
574 m=fMasses->At(j-1);
575 fDresult=fDmasses->At(j-1);
576 return m;
577}
578///////////////////////////////////////////////////////////////////////////
579Double_t AliTrack::GetMassHypothesisProb(Int_t j)
580{
581// Provide the probability of the jth hypothesis for this track.
582// Note : the first hypothesis is indicated by j=1.
583// Default : j=0 ==> Hypothesis with highest probability.
584
585 Double_t prob=0;
586
587 // Check validity of index j
588 if (j<0 || j>fNmasses)
589 {
590 cout << " *AliTrack::GetMassHypothesisProb* Invalid index j : " << j
591 << " Number of mass hypotheses : " << fNmasses << endl;
592 return 0;
593 }
594
595 // Select mass hypothesis with highest probability
596 if (j==0)
597 {
598 if (fNmasses)
599 {
600 prob=fPmasses->At(0);
601 for (Int_t i=1; i<fNmasses; i++)
602 {
603 if (fPmasses->At(i)>prob) prob=fPmasses->At(i);
604 }
605 }
606 return prob;
607 }
608
609 // Provide probability of requested mass hypothesis
610 prob=fPmasses->At(j-1);
611 return prob;
612}
613///////////////////////////////////////////////////////////////////////////
614void AliTrack::SetMass()
615{
616// Set the mass and error to the value of the hypothesis with highest prob.
617
618 Double_t m=0,dm=0,prob=0;
619
620 // Select mass hypothesis with highest probability
621 if (fNmasses)
622 {
623 m=fMasses->At(0);
624 dm=fDmasses->At(0);
625 prob=fPmasses->At(0);
626 for (Int_t i=1; i<fNmasses; i++)
627 {
628 if (fPmasses->At(i)>prob)
629 {
630 m=fMasses->At(i);
631 dm=fDmasses->At(i);
632 }
633 }
634 SetMass(m,dm);
635 }
636 else
637 {
638 cout << " *AliTrack::SetMass()* No hypothesis present => No action." << endl;
639 }
640}
641///////////////////////////////////////////////////////////////////////////
642void AliTrack::RemoveMassHypothesis(Int_t j)
643{
644// Remove the jth mass hypothesis for this track.
645// Note : the first hypothesis is indicated by j=1.
646
647 if (j<=0 || j>fNmasses) // Check validity of index j
648 {
649 cout << " *AliTrack::RemoveMassHypothesis* Invalid index j : " << j
650 << " Number of mass hypotheses : " << fNmasses << endl;
651 }
652 else
653 {
654 if (j != fNmasses)
655 {
656 fMasses->AddAt(fMasses->At(fNmasses-1),j-1);
657 fDmasses->AddAt(fDmasses->At(fNmasses-1),j-1);
658 fPmasses->AddAt(fPmasses->At(fNmasses-1),j-1);
659 }
660 fMasses->AddAt(0,fNmasses-1);
661 fDmasses->AddAt(0,fNmasses-1);
662 fPmasses->AddAt(0,fNmasses-1);
663 fNmasses--;
664 fMasses->Set(fNmasses);
665 fDmasses->Set(fNmasses);
666 fPmasses->Set(fNmasses);
667 }
668}
669///////////////////////////////////////////////////////////////////////////