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