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 //====================================================================================================================================================