Config.C updated
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //====================================================================================================================================================
17 //
18 //      Description of an ALICE muon forward track, combining the information of the Muon Spectrometer and the Muon Forward Tracker
19 //
20 //      Contact author: antonio.uras@cern.ch
21 //
22 //====================================================================================================================================================
23
24 #include "AliLog.h"
25 #include "AliMUONTrack.h"
26 #include "AliMFTCluster.h"
27 #include "AliMUONVCluster.h"
28 #include "AliMUONTrackParam.h"
29 #include "AliMUONTrackExtrap.h"
30 #include "TClonesArray.h"
31 #include "TMatrixD.h"
32 #include "TParticle.h"
33 #include "AliMuonForwardTrack.h"
34 #include "AliMFTConstants.h"
35 #include "TLorentzVector.h"
36 #include "TDatabasePDG.h"
37 #include "AliMUONConstants.h"
38
39 ClassImp(AliMuonForwardTrack)
40
41 //====================================================================================================================================================
42
43 AliMuonForwardTrack::AliMuonForwardTrack():
44   AliMUONTrack(),
45   fMUONTrack(0),
46   fMCTrackRef(0),
47   fMFTClusters(0),
48   fNWrongClustersMC(-1),
49   fTrackMCId(-1),
50   fKinem(0,0,0,0),
51   fParamCovMatrix(5,5)
52 {
53
54   // default constructor
55   
56   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
57   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
58     fParentMCLabel[iParent] = -1;
59     fParentPDGCode[iParent] =  0;
60   }
61   fMFTClusters = new TClonesArray("AliMFTCluster");
62   fMFTClusters -> SetOwner(kTRUE);
63
64 }
65
66 //====================================================================================================================================================
67
68 AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
69   AliMUONTrack(),
70   fMUONTrack(0),
71   fMCTrackRef(0),
72   fMFTClusters(0),
73   fNWrongClustersMC(-1),
74   fTrackMCId(-1),
75   fKinem(0,0,0,0),
76   fParamCovMatrix(5,5)
77 {
78
79   SetMUONTrack(MUONTrack);
80   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
81   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
82     fParentMCLabel[iParent] = -1;
83     fParentPDGCode[iParent] =  0;
84   }
85   fMFTClusters = new TClonesArray("AliMFTCluster");
86   fMFTClusters -> SetOwner(kTRUE);
87
88 }
89
90 //====================================================================================================================================================
91
92 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track): 
93   AliMUONTrack(track),
94   fMUONTrack(0x0),
95   fMCTrackRef(0x0),
96   fMFTClusters(0x0),
97   fNWrongClustersMC(track.fNWrongClustersMC),
98   fTrackMCId(track.fTrackMCId),
99   fKinem(track.fKinem),
100   fParamCovMatrix(track.fParamCovMatrix)
101 {
102
103   // copy constructor
104   fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
105   if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
106   fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
107   fMFTClusters->SetOwner(kTRUE);
108   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
109   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
110     fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
111     fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
112   }
113   
114 }
115
116 //====================================================================================================================================================
117
118 AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& track) {
119
120   // Asignment operator
121
122   // check assignement to self
123   if (this == &track) return *this;
124
125   // base class assignement
126   AliMUONTrack::operator=(track);
127   
128   // clear memory
129   Clear("");
130   
131   fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
132   if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
133   fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
134   fMFTClusters->SetOwner(kTRUE);
135   fNWrongClustersMC = track.fNWrongClustersMC;
136   fTrackMCId        = track.fTrackMCId;
137   fKinem            = track.fKinem;
138   fParamCovMatrix   = track.fParamCovMatrix;
139
140   for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
141   for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
142     fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
143     fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
144   }
145   
146   return *this;
147
148 }
149
150 //====================================================================================================================================================
151
152 void AliMuonForwardTrack::Clear(const Option_t* /*opt*/) {
153
154   // Clear arrays
155   fMFTClusters -> Delete(); 
156   delete fMFTClusters; fMFTClusters = 0x0;
157   delete fMUONTrack;   fMUONTrack   = 0x0;
158   delete fMCTrackRef;  fMCTrackRef  = 0x0;
159   
160 }
161
162 //====================================================================================================================================================
163
164 AliMuonForwardTrack::~AliMuonForwardTrack() {
165
166   delete fMUONTrack;
167   delete fMCTrackRef;
168   fMFTClusters -> Delete();
169   delete fMFTClusters;
170   
171 }
172
173 //====================================================================================================================================================
174
175 void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
176
177   if (fMUONTrack) {
178     AliInfo("fMUONTrack already exists, nothing will be done");
179     return;
180   }
181
182   fMUONTrack = MUONTrack;
183   SetGlobalChi2(fMUONTrack->GetGlobalChi2());
184   
185 }
186
187 //====================================================================================================================================================
188
189 void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
190
191   if (fMCTrackRef) {
192     AliInfo("fMCTrackRef already exists, nothing will be done");
193     return;
194   }
195
196   fMCTrackRef = MCTrackRef;
197   
198 }
199
200 //====================================================================================================================================================
201
202 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
203
204   AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
205   Int_t iMFTCluster = this->fMFTClusters->GetEntries();
206   AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f", 
207            mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
208   AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
209   AliDebug(1, Form("Created MUON cluster %p", muonCluster));
210   trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
211   AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
212   AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
213   // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
214   Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
215   AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
216   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
217   mftCluster.SetLocalChi2(chi2Kalman);
218   mftCluster.SetTrackChi2(newGlobalChi2);
219   new ((*(this->fMFTClusters))[iMFTCluster]) AliMFTCluster(mftCluster);
220   AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries while this->fMFTClusters=%p has %d entries",
221                    GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries(), this->fMFTClusters, this->fMFTClusters->GetEntries()));
222   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
223   SetGlobalChi2(newGlobalChi2);
224   AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
225   ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
226   for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
227     AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
228                      iPar, ((AliMUONTrackParam*)(GetTrackParamAtCluster()->At(iPar)))->GetClusterPtr()));
229   }
230
231 }
232
233 //====================================================================================================================================================
234
235 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMUONCluster(Int_t iMUONCluster) {
236
237   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
238     AliError("Invalid MUON cluster index. NULL pointer will be returned");
239     return NULL;
240   }
241
242   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) fMUONTrack->GetTrackParamAtCluster()->At(iMUONCluster); 
243
244   return trackParam;
245
246 }
247
248 //====================================================================================================================================================
249
250 AliMUONTrackParam* AliMuonForwardTrack::GetTrackParamAtMFTCluster(Int_t iMFTCluster) {
251
252   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
253     AliError("Invalid MFT cluster index. NULL pointer will be returned");
254     return NULL;
255   }
256
257   AliMUONTrackParam *trackParam = (AliMUONTrackParam*) GetTrackParamAtCluster()->At(iMFTCluster); 
258
259   return trackParam;
260
261 }
262
263 //====================================================================================================================================================
264
265 AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
266
267   if (iMUONCluster<0 || iMUONCluster>=GetNMUONClusters()) {
268     AliError("Invalid MUON cluster index. NULL pointer will be returned");
269     return NULL;
270   }
271
272   AliMUONTrackParam *trackParam = GetTrackParamAtMUONCluster(iMUONCluster);
273   AliMUONVCluster *muonCluster = trackParam->GetClusterPtr();
274
275   return muonCluster;
276
277 }
278
279 //====================================================================================================================================================
280
281 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
282
283   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
284     AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
285     return NULL;
286   }
287
288   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
289   AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
290
291   return mftCluster;
292
293 }
294
295 //====================================================================================================================================================
296
297 Double_t AliMuonForwardTrack::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster) {
298
299   AliDebug(1, Form("Running Kalman filter for parameters %p (z = %f) and cluster %p (z = %f)", 
300                    &trackParamAtCluster, trackParamAtCluster.GetZ(), trackParamAtCluster.GetClusterPtr(), trackParamAtCluster.GetClusterPtr()->GetZ()));
301   
302   // Compute new track parameters and their covariances including new cluster using kalman filter
303   // return the *additional* track chi2
304   
305   // Get actual track parameters (p)
306   TMatrixD param(trackParamAtCluster.GetParameters());
307   
308   // Get new cluster parameters (m)
309   AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
310   AliDebug(1, Form("cluster->GetX() = %f  cluster->GetY() = %f  cluster->GetErrX() = %f  cluster->GetErrY() = %f", 
311                    cluster->GetX(), cluster->GetY(), cluster->GetErrX(), cluster->GetErrY()));
312   TMatrixD clusterParam(5,1);
313   clusterParam.Zero();
314   clusterParam(0,0) = cluster->GetX();
315   clusterParam(2,0) = cluster->GetY();
316
317   // Compute the current parameter weight (W)
318   TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
319   if (paramWeight.Determinant() != 0) {
320     paramWeight.Invert();
321   } else {
322     Warning("RunKalmanFilter"," Determinant = 0");
323     return 1.e10;
324   }
325   
326   // Compute the new cluster weight (U)
327   TMatrixD clusterWeight(5,5);
328   clusterWeight.Zero();
329   clusterWeight(0,0) = 1. / cluster->GetErrX2();
330   clusterWeight(2,2) = 1. / cluster->GetErrY2();
331
332   // Compute the new parameters covariance matrix ( (W+U)^-1 )
333   TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
334   if (newParamCov.Determinant() != 0) {
335     newParamCov.Invert();
336   } 
337   else {
338     Warning("RunKalmanFilter"," Determinant = 0");
339     return 1.e10;
340   }
341   
342   // Save the new parameters covariance matrix
343   trackParamAtCluster.SetCovariances(newParamCov);
344   
345   // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
346   TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
347   TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp);    // U(m-p)
348   TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
349   newParam += param;                                   // ((W+U)^-1)U(m-p) + p
350
351   // Save the new parameters
352   trackParamAtCluster.SetParameters(newParam);
353   
354   // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
355   tmp = newParam; // p'
356   tmp -= param;   // (p'-p)
357   TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp);           // W(p'-p)
358   TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
359   tmp = newParam;      // p'
360   tmp -= clusterParam; // (p'-m)
361   TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp);            // U(p'-m)
362   addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
363   
364   AliDebug(1,Form("Adding Kalman chi2 = %f",addChi2Track(0,0)));
365
366   return addChi2Track(0,0);
367   
368 }
369
370 //====================================================================================================================================================
371
372 Double_t AliMuonForwardTrack::GetWeightedOffset(Double_t x, Double_t y, Double_t z) {
373
374   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
375   AliMUONTrackExtrap::ExtrapToZCov(param, z);
376
377   TMatrixD cov(5,5);
378   cov = param->GetCovariances();
379
380   TMatrixD covCoordinates(2,2);
381   covCoordinates(0,0) = cov(0,0);
382   covCoordinates(0,1) = cov(0,2);
383   covCoordinates(1,0) = cov(2,0);
384   covCoordinates(1,1) = cov(2,2);
385   
386   TMatrixD covCoordinatesInverse = covCoordinates.Invert();
387
388   Double_t dX = param->GetNonBendingCoor() - x;
389   Double_t dY = param->GetBendingCoor()    - y;
390   
391   Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + 
392                                              dY*dY*covCoordinatesInverse(1,1) + 
393                                              2.*dX*dY*covCoordinatesInverse(0,1)));
394
395   return weightedOffset;
396
397 }
398
399 //====================================================================================================================================================
400
401 Double_t AliMuonForwardTrack::GetOffsetX(Double_t x, Double_t z) {
402
403   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
404   AliMUONTrackExtrap::ExtrapToZCov(param, z);
405   Double_t dX = param->GetNonBendingCoor() - x;
406   return dX;
407
408 }
409
410 //====================================================================================================================================================
411
412 Double_t AliMuonForwardTrack::GetOffsetY(Double_t y, Double_t z) {
413
414   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
415   AliMUONTrackExtrap::ExtrapToZCov(param, z);
416   Double_t dY = param->GetBendingCoor() - y;
417   return dY;
418
419 }
420
421 //====================================================================================================================================================
422
423 Double_t AliMuonForwardTrack::GetOffset(Double_t x, Double_t y, Double_t z) {
424
425   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
426   AliMUONTrackExtrap::ExtrapToZCov(param, z);
427   Double_t dX = param->GetNonBendingCoor() - x;
428   Double_t dY = param->GetBendingCoor()    - y;
429   return TMath::Sqrt(dX*dX + dY*dY);
430
431 }
432
433 //====================================================================================================================================================
434
435 Double_t AliMuonForwardTrack::GetDCA(Double_t x, Double_t y, Double_t z) {
436
437   // Distance of Closest Approach, according to the standard MUON terminology. Actually, the offset of the track w.r.t. the primary vertex,
438   // where the extrapolation of the track DOES NOT include the MFT information
439
440   AliMUONTrackParam param = (*((AliMUONTrackParam*)(GetTrackParamAtMUONCluster(0))));
441
442   AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&param, z);
443   Double_t dX = param.GetNonBendingCoor() - x;
444   Double_t dY = param.GetBendingCoor()    - y;
445   return TMath::Sqrt(dX*dX + dY*dY);
446
447 }
448
449 //====================================================================================================================================================
450
451 Double_t AliMuonForwardTrack::GetMomentumSpectrometer(Double_t z) {
452
453   // Momentum of the track at the primary vertex plane, where the extrapolation of the track DOES NOT include the MFT information
454
455   AliMUONTrackParam param = (*((AliMUONTrackParam*)(GetTrackParamAtMUONCluster(0))));
456
457   AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&param, z);
458   return param.P();
459
460 }
461
462 //====================================================================================================================================================
463
464 Double_t AliMuonForwardTrack::GetThetaAbs() {
465
466   // it is the angle defined by the imaginary line goingo from the vertex to the exit point of the track at the end of the hadron absorber
467   
468   Double_t z = AliMUONConstants::AbsZEnd();
469
470   AliMUONTrackParam param = (*((AliMUONTrackParam*)(GetTrackParamAtMUONCluster(0))));
471
472   AliMUONTrackExtrap::ExtrapToZ(&param, z);
473   Double_t x = param.GetNonBendingCoor();
474   Double_t y = param.GetBendingCoor();
475
476   return 180. +TMath::ATan(TMath::Sqrt(x*x + y*y)/z)*TMath::RadToDeg();
477
478 }
479
480 //====================================================================================================================================================
481
482 Bool_t AliMuonForwardTrack::IsFromDirectResonance() {
483
484   if (!IsMuon()) return kFALSE;
485
486   for (Int_t i=GetFirstMotherID(); i>=0; i--) if (!IsPDGResonance(GetParentPDGCode(i))) return kFALSE;  // the decay chain finds a non-resonance particle
487
488   return kTRUE; 
489   
490 }
491
492 //====================================================================================================================================================
493
494 Bool_t AliMuonForwardTrack::IsFromChainResonance() {
495
496   if (!IsMuon()) return kFALSE;
497
498   if (GetFirstMotherID() == 0) return kFALSE;   // it is not a chain
499
500   if (!IsPDGResonance(GetParentPDGCode(GetFirstMotherID()))) return kFALSE;  // primordial is not a resonance
501
502   for (Int_t i=GetFirstMotherID()-1; i>=0; i--) if (!IsPDGResonance(GetParentPDGCode(i))) return kTRUE;  // the decay chain finds a non-resonance particle
503
504   return kFALSE;
505   
506 }
507
508 //====================================================================================================================================================
509
510 Bool_t AliMuonForwardTrack::IsFromDirectCharm() {
511
512   if (!IsMuon()) return kFALSE;
513
514   for (Int_t i=GetFirstMotherID(); i>=0; i--) if (!IsPDGCharm(GetParentPDGCode(i))) return kFALSE;  // the decay chain finds a non-charmed particle
515
516   return kTRUE; 
517
518 }
519
520 //====================================================================================================================================================
521
522 Bool_t AliMuonForwardTrack::IsFromChainCharm() {
523
524   if (!IsMuon()) return kFALSE;
525
526   if (GetFirstMotherID() == 0) return kFALSE;   // it is not a chain
527
528   if (!IsPDGCharm(GetParentPDGCode(GetFirstMotherID()))) return kFALSE;  // primordial is not a charmed hadron
529
530   for (Int_t i=GetFirstMotherID()-1; i>=0; i--) if (!IsPDGCharm(GetParentPDGCode(i))) return kTRUE;  // the decay chain finds a non-charmed particle
531
532   return kFALSE;
533
534 }
535
536 //====================================================================================================================================================
537
538 Bool_t AliMuonForwardTrack::IsFromDirectBeauty() {
539
540   if (!IsMuon()) return kFALSE;
541
542   for (Int_t i=GetFirstMotherID(); i>=0; i--) if (!IsPDGBeauty(GetParentPDGCode(i))) return kFALSE;  // the decay chain finds a non-beauty particle
543
544   return kTRUE; 
545
546 }
547
548 //====================================================================================================================================================
549
550 Bool_t AliMuonForwardTrack::IsFromChainBeauty() {
551
552   if (!IsMuon()) return kFALSE;
553
554   if (GetFirstMotherID() == 0) return kFALSE;   // it is not a chain
555
556   if (!IsPDGBeauty(GetParentPDGCode(GetFirstMotherID()))) return kFALSE;  // primordial is not a beauty hadron
557
558   for (Int_t i=GetFirstMotherID()-1; i>=0; i--) if (!IsPDGBeauty(GetParentPDGCode(i))) return kTRUE;  // the decay chain finds a non-beauty particle
559
560   return kFALSE;
561
562 }
563
564 //====================================================================================================================================================
565
566 Bool_t AliMuonForwardTrack::IsMuon() {
567
568   if (IsFake()) return kFALSE;
569
570   if (TMath::Abs(fMCTrackRef->GetPdgCode())==13) return kTRUE;
571
572   return kFALSE;
573
574 }
575
576 //====================================================================================================================================================
577
578 Bool_t AliMuonForwardTrack::IsFake() {
579
580   if (!fMCTrackRef) return kTRUE;
581
582   return kFALSE;
583
584 }
585
586 //====================================================================================================================================================
587
588 Bool_t AliMuonForwardTrack::IsPDGResonance(Int_t pdg) {
589
590   // if (pdg<10) return kFALSE; 
591   // Int_t id = pdg%100000; 
592   // return (!((id-id%10)%110));
593
594   Bool_t result = kFALSE;
595
596   if ( pdg ==    113 ||
597        pdg ==    221 ||
598        pdg ==    223 ||
599        pdg ==    331 ||
600        pdg ==    333 ||
601        pdg ==    443 ||
602        pdg == 100443 ||
603        pdg ==    553 ||
604        pdg == 100553 ) result = kTRUE;
605   
606   return result; 
607   
608 }
609
610 //====================================================================================================================================================
611
612 Bool_t AliMuonForwardTrack::IsPDGCharm(Int_t pdg) {
613
614   Bool_t result = kFALSE;
615
616   if ( TMath::Abs(pdg) ==   411 ||
617        TMath::Abs(pdg) ==   421 ||
618        TMath::Abs(pdg) == 10411 ||
619        TMath::Abs(pdg) == 10421 ||
620        TMath::Abs(pdg) ==   413 ||
621        TMath::Abs(pdg) ==   423 ||
622        TMath::Abs(pdg) == 10413 ||
623        TMath::Abs(pdg) == 10423 ||
624        TMath::Abs(pdg) == 20413 ||
625        TMath::Abs(pdg) == 20423 ||
626        TMath::Abs(pdg) ==   415 ||
627        TMath::Abs(pdg) ==   425 ||
628        TMath::Abs(pdg) ==   431 ||
629        TMath::Abs(pdg) == 10431 ||
630        TMath::Abs(pdg) ==   433 ||
631        TMath::Abs(pdg) == 10433 ||
632        TMath::Abs(pdg) == 20433 ||
633        TMath::Abs(pdg) ==   435 ) result = kTRUE;
634   
635   return result; 
636   
637 }
638
639 //====================================================================================================================================================
640
641 Bool_t AliMuonForwardTrack::IsPDGBeauty(Int_t pdg) {
642
643   Bool_t result = kFALSE;
644
645   if ( TMath::Abs(pdg) ==   511 ||
646        TMath::Abs(pdg) ==   521 ||
647        TMath::Abs(pdg) == 10511 ||
648        TMath::Abs(pdg) == 10521 ||
649        TMath::Abs(pdg) ==   513 ||
650        TMath::Abs(pdg) ==   523 ||
651        TMath::Abs(pdg) == 10513 ||
652        TMath::Abs(pdg) == 10523 ||
653        TMath::Abs(pdg) == 20513 ||
654        TMath::Abs(pdg) == 20523 ||
655        TMath::Abs(pdg) ==   515 ||
656        TMath::Abs(pdg) ==   525 ||
657        TMath::Abs(pdg) ==   531 ||
658        TMath::Abs(pdg) == 10531 ||
659        TMath::Abs(pdg) ==   533 ||
660        TMath::Abs(pdg) == 10533 ||
661        TMath::Abs(pdg) == 20533 ||
662        TMath::Abs(pdg) ==   535 ||
663        TMath::Abs(pdg) ==   541 ||
664        TMath::Abs(pdg) == 10541 ||
665        TMath::Abs(pdg) ==   543 ||
666        TMath::Abs(pdg) == 10543 ||
667        TMath::Abs(pdg) == 20543 ||
668        TMath::Abs(pdg) ==   545 ) result = kTRUE;
669   
670   return result; 
671   
672 }
673
674 //====================================================================================================================================================
675
676 Bool_t AliMuonForwardTrack::IsMuonFromBackground() {
677
678   Bool_t result = kFALSE;
679
680   if (!IsMuon()) return result;
681
682   if (!IsFromDirectResonance() && !IsFromChainResonance() && !IsFromDirectCharm() && !IsFromChainCharm() && !IsFromDirectBeauty() && !IsFromChainBeauty()) result = kTRUE;
683
684   if (result) AliDebug(1, Form("Muon comes from a background source %d", GetParentPDGCode(0)));
685
686   return result;
687
688 }
689
690 //====================================================================================================================================================
691
692 void AliMuonForwardTrack::EvalKinem(Double_t z) {
693
694   AliMUONTrackParam *param = GetTrackParamAtMFTCluster(0);
695   AliMUONTrackExtrap::ExtrapToZCov(param, z);
696
697   Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
698   Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
699   fKinem.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
700
701   TMatrixD cov(5,5);
702   fParamCovMatrix = param->GetCovariances();
703
704 }
705
706 //====================================================================================================================================================
707
708 Int_t AliMuonForwardTrack::GetFirstMotherID() {
709
710   Int_t motherLevel = 0; 
711   Int_t motherMCLabel = GetParentMCLabel(motherLevel); 
712   while (motherMCLabel >= 0) { 
713     motherLevel++; 
714     motherMCLabel = GetParentMCLabel(motherLevel); 
715   }
716   return motherLevel-1;
717
718 }
719
720 //====================================================================================================================================================
721
722 void AliMuonForwardTrack::PrintHistory() {
723
724   if (IsFake()) printf("Track is a fake MUON\n");
725   else {
726     TString history = "";
727     for (Int_t i=GetFirstMotherID(); i>=0; i--) {
728       history += TDatabasePDG::Instance()->GetParticle(GetParentPDGCode(i))->GetName();
729       history += " -> ";
730     }
731     history += TDatabasePDG::Instance()->GetParticle(fMCTrackRef->GetPdgCode())->GetName();
732     printf("%s\n",history.Data());
733   }
734
735 }
736
737 //====================================================================================================================================================