PWG3/muon/AliMCMuonTrack.cxx

   1 #include <TParticle.h>
   2 #include <TLorentzVector.h>
   3 #include <TClonesArray.h>
   4
   5 #include "AliMCEventHandler.h"
   6 #include "AliMCEvent.h"
   7 #include "AliStack.h"
   8 #include "AliAODMCParticle.h"
   9 #include "AliESDEvent.h"
  10 #include "AliESDMuonTrack.h"
  11 #include "AliESDMuonCluster.h"
  12 #include "AliAODTrack.h"
  13 #include "AliAODMuonTrack.h"
  14 #include "AliMCMuonTrack.h"
  15
  16 ClassImp(AliMCMuonTrack)
  17
  18 const Double_t AliMCMuonTrack::fgkSigmaCut = 10.;
  19
  20 //-----------------------------------------------------------------------------
  21 AliMCMuonTrack::AliMCMuonTrack() :
  22 AliAODMuonTrack(),
  23 fIsFull(kFALSE),
  24 fPGen(),
  25 fTrackIndex(-1),
  26 fTrackPDGCode(0),
  27 fSource(-1),
  28 fNParents(0),
  29 fOscillation(kFALSE),
  30 fWeight(0.)
  31 {
  32   //
  33   // default constructor
  34   //
  35   for (Int_t i=0; i<fgkNParentsMax; i++) {
  36     fParentIndex[i] = -1;
  37     fParentPDGCode[i] = 0;
  38   }
  39   for (Int_t i=0; i<4; i++) {
  40     fQuarkIndex[i] = -1;
  41     fQuarkPDGCode[i] = 0;
  42   }
  43 }
  44
  45 //-----------------------------------------------------------------------------
  46 AliMCMuonTrack::AliMCMuonTrack(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full) :
  47 AliAODMuonTrack(trkAOD),
  48 fIsFull(full),
  49 fPGen(),
  50 fTrackIndex(-1),
  51 fTrackPDGCode(0),
  52 fSource(-1),
  53 fNParents(0),
  54 fOscillation(kFALSE),
  55 fWeight(0.)
  56 {
  57   //
  58   // default constructor
  59   //
  60   AliAODMCParticle *pMC = this->FindTrackRef(trkAOD, mcClArr);
  61   if (pMC) this->SetMCInfo(pMC, mcClArr);
  62 }
  63
  64 //-----------------------------------------------------------------------------
  65 AliMCMuonTrack::AliMCMuonTrack(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full) :
  66 AliAODMuonTrack(trkESD),
  67 fIsFull(full),
  68 fPGen(),
  69 fTrackIndex(-1),
  70 fTrackPDGCode(0),
  71 fSource(-1),
  72 fNParents(0),
  73 fOscillation(kFALSE),
  74 fWeight(0.)
  75 {
  76   //
  77   // default constructor
  78   //
  79   TParticle *pMC = this->FindTrackRef(trkESD, esd, mcH);  // do nothing when running on official train
  80   if (pMC) this->SetMCInfo(pMC, mcH);
  81   // MC infomation is not implement when running on the official train for they depend on the MUON module.
  82   // In the case of private running, one can get the information from ESD by uncommenting the definitions
  83   // of FindTrackRef() and CovESDtoMuonTrack()
  84   // and including the following head files in both AliMCMuonTrack.h and this file,
  85   // #include "AliMUONRecoCheck.h"
  86   // #include "AliMUONVTrackStore.h"
  87   // #include "AliMUONVClusterStore.h"
  88   // #include "AliMUONTrackParam.h"
  89   // #include "AliMUONTrack.h"
  90   // #include "AliMUONVCluster.h"
  91   // #include "AliMUONESDInterface.h"
  92 }
  93
  94
  95
  96 //-----------------------------------------------------------------------------
  97 AliMCMuonTrack::~AliMCMuonTrack()
  98 {
  99   //
 100   // destructor
 101   //
 102 }
 103
 104 //-----------------------------------------------------------------------------
 105 AliAODMCParticle* AliMCMuonTrack::FindTrackRef(AliAODTrack *trkAOD, TClonesArray *mcClArr)
 106 {
 107   AliAODMCParticle *pMC = 0;
 108   fTrackIndex = trkAOD->GetLabel();
 109   if (fTrackIndex>=0)
 110     pMC = (AliAODMCParticle*)mcClArr->At(fTrackIndex);
 111   return pMC;
 112 }
 113
 114 //-----------------------------------------------------------------------------
 115 TParticle* AliMCMuonTrack::FindTrackRef(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH)
 116 {
 117   TParticle *pMCRef = 0;
 118
 119   /*AliMUONTrack *trkMuon = CovESDtoMuonTrack(*trkESD);
 120   AliMUONRecoCheck rc(esd,mcH);
 121   AliMUONVTrackStore* trkRefArr = rc.TrackRefs(-1);
 122   TIter next(trkRefArr->CreateIterator());
 123   AliMUONTrack* trkRef = 0;
 124   while ((trkRef=static_cast<AliMUONTrack*>(next()))) {
 125     Bool_t trkCompArr[10];
 126     Int_t nMatchClusters = trkMuon->CompatibleTrack(trkRef, fgkSigmaCut, trkCompArr);
 127     Double_t matchClusterFrac = ((Double_t)nMatchClusters) / ((Double_t)trkMuon->GetNClusters());
 128     if ((trkCompArr[0] || trkCompArr[1] || trkCompArr[2] || trkCompArr[3]) &&
 129         (trkCompArr[6] || trkCompArr[7] || trkCompArr[8] || trkCompArr[9]) && matchClusterFrac>0.5) {
 130       fTrackIndex = trkRef->GetUniqueID();
 131       trkRefArr->Remove(*trkRef);
 132       break;
 133     }
 134   }
 135   if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);*/
 136
 137   return pMCRef;
 138 }
 139
 140 //-----------------------------------------------------------------------------
 141 void AliMCMuonTrack::SetMCInfo(AliAODMCParticle *pMC, TClonesArray *mcClArr)
 142 {
 143   fPGen.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
 144   fTrackPDGCode = pMC->GetPdgCode();
 145   if (TMath::Abs(fTrackPDGCode)!=13) {
 146     fSource = 4;
 147     return;
 148   }
 149
 150   Int_t lineM = pMC->GetMother();
 151   if (lineM<0) {
 152     fSource = 2;
 153     return;
 154   }
 155
 156   Bool_t isPrimary =
 157     ((AliAODMCParticle*)mcClArr->At(lineM))->IsPrimary();
 158   if (!isPrimary) {
 159     fSource = 3;
 160     return;
 161   }
 162
 163   this->FillHistoryParents(lineM, mcClArr);
 164   fSource = this->SelectHFMuon();
 165   return;
 166 }
 167
 168 //-----------------------------------------------------------------------------
 169 void AliMCMuonTrack::SetMCInfo(TParticle *pMC, AliMCEventHandler *mcH)
 170 {
 171   fPGen.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
 172   fTrackPDGCode = pMC->GetPdgCode();
 173   if (TMath::Abs(fTrackPDGCode)!=13) {
 174     fSource = 4;
 175     return;
 176   }
 177
 178   Int_t lineM = pMC->GetFirstMother();
 179   if (lineM<0) {
 180     fSource = 2;
 181     return;
 182   }
 183
 184   AliStack *stack = mcH->MCEvent()->Stack();
 185   if (lineM>=stack->GetNprimary()) {
 186     fSource = 3;
 187     return;
 188   }
 189
 190   this->FillHistoryParents(lineM, stack);
 191   fSource = this->SelectHFMuon();
 192   return;
 193 }
 194
 195 //-----------------------------------------------------------------------------
 196 void AliMCMuonTrack::FillHistoryParents(Int_t lineM, TClonesArray *mcClArr)
 197 {
 198   Int_t countP=-1, pdg=0;
 199   Int_t parents[10], parLine[10];
 200   AliAODMCParticle *mother = 0;
 201   while(lineM>=0){
 202     mother = (AliAODMCParticle*)mcClArr->At(lineM);
 203     pdg = mother->GetPdgCode();
 204     if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
 205     parents[++countP] = pdg;
 206     parLine[countP] = lineM;
 207     lineM = mother->GetMother();
 208   }
 209   for(Int_t i=0; i<=countP; i++){
 210     fParentIndex[i] = parLine[countP-i];
 211     fParentPDGCode[i] = parents[countP-i];
 212   }
 213   fNParents = countP + 1;
 214
 215   if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, mcClArr);
 216   return;
 217 }
 218
 219 //-----------------------------------------------------------------------------
 220 void AliMCMuonTrack::FillHistoryParents(Int_t lineM, AliStack *stack)
 221 {
 222   Int_t countP=-1, pdg=0;
 223   Int_t parents[10], parLine[10];
 224   TParticle *mother = 0;
 225   while(lineM>=0){
 226     mother = stack->Particle(lineM);
 227     pdg = mother->GetPdgCode();
 228     if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
 229     parents[++countP] = pdg;
 230     parLine[countP] = lineM;
 231     lineM = mother->GetFirstMother();
 232   }
 233   for(Int_t i=0; i<=countP; i++){
 234     fParentIndex[i] = parLine[countP-i];
 235     fParentPDGCode[i] = parents[countP-i];
 236   }
 237   fNParents = countP + 1;
 238
 239   if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, stack);
 240   return;
 241 }
 242
 243 //-----------------------------------------------------------------------------
 244 void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, TClonesArray *mcClArr)
 245 {
 246   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 247   if (lineM<0) return;
 248   Int_t countP=-1, pdg=0;
 249   AliAODMCParticle *mother = 0;
 250   while(lineM>=0){
 251     mother = (AliAODMCParticle*)mcClArr->At(lineM);
 252     pdg = mother->GetPdgCode();
 253     fQuarkIndex[++countP] = lineM;
 254     fQuarkPDGCode[countP] = pdg;
 255     lineM = mother->GetMother();
 256   }
 257
 258   // for PYTHIA checking
 259   countP = 1;
 260   for(Int_t par=0; par<4; par++) {
 261     if(TMath::Abs(this->GetQuarkPDGCode(par))<6) { countP=par; break; }
 262   }
 263   if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)) {
 264     if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
 265       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 266                  this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));
 267
 268       pdg = this->GetQuarkPDGCode(countP);
 269       Int_t line = this->GetQuarkIndex(countP);
 270       this->ResetQuarkInfo();
 271       while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
 272         pdg = ((AliAODMCParticle*)mcClArr->At(++line))->GetPdgCode();
 273       }
 274       while(line>=0){
 275         mother = (AliAODMCParticle*)mcClArr->At(line);
 276         pdg = mother->GetPdgCode();
 277         fQuarkIndex[countP] = line;
 278         fQuarkPDGCode[countP++] = pdg;
 279         line = mother->GetMother();
 280       }
 281     }
 282   }
 283   return;
 284 }
 285
 286 //-----------------------------------------------------------------------------
 287 void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, AliStack *stack)
 288 {
 289   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 290   if (lineM<0) return;
 291   Int_t countP=-1, pdg=0;
 292   TParticle *mother = 0;
 293   while(lineM>=0){
 294     mother = stack->Particle(lineM);
 295     pdg = mother->GetPdgCode();
 296     fQuarkIndex[++countP] = lineM;
 297     fQuarkPDGCode[countP] = pdg;
 298     lineM = mother->GetFirstMother();
 299   }
 300
 301   // for PYTHIA checking
 302   countP = 1;
 303   for(Int_t par=0; par<4; par++) {
 304     if(TMath::Abs(this->GetQuarkPDGCode(par))<6) { countP=par; break; }
 305   }
 306   if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)) {
 307     if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
 308       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 309                  this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));
 310
 311       pdg = this->GetQuarkPDGCode(countP);
 312       Int_t line = this->GetQuarkIndex(countP);
 313       this->ResetQuarkInfo();
 314       while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
 315         pdg = stack->Particle(++line)->GetPdgCode();
 316       }
 317       while(line>=0){
 318         mother = stack->Particle(line);
 319         pdg = mother->GetPdgCode();
 320         fQuarkIndex[countP] = line;
 321         fQuarkPDGCode[countP++] = pdg;
 322         line = mother->GetFirstMother();
 323       }
 324     }
 325   }
 326   return;
 327 }
 328
 329 //-----------------------------------------------------------------------------
 330 /*AliMUONTrack* AliMCMuonTrack::CovESDtoMuonTrack(AliESDMuonTrack &trkESD)
 331 {
 332   // method in $ALICE_ROOT/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx
 333   AliMUONTrack *trkMuon = new AliMUONTrack();
 334   if(!trkESD.ClustersStored()) return trkMuon;
 335
 336   AliMUONTrackParam param;
 337   AliMUONESDInterface::GetParamAtFirstCluster(trkESD, param);
 338   AliMUONESDInterface::GetParamCov(trkESD, param);
 339   AliMUONVClusterStore *cStore = AliMUONESDInterface::NewClusterStore();
 340   AliMUONVCluster *cluster = cStore->CreateCluster(0,0,0);
 341
 342   AliESDMuonCluster *esdCluster = (AliESDMuonCluster*)trkESD.GetClusters().First();
 343   while (esdCluster) {
 344     AliMUONESDInterface::ESDToMUON(*esdCluster, *cluster);
 345     param.SetZ(cluster->GetZ());
 346     trkMuon->AddTrackParamAtCluster(param, *cluster, kTRUE);
 347     esdCluster = (AliESDMuonCluster*)trkESD.GetClusters().After(esdCluster);
 348   }
 349
 350   delete cluster; cluster = 0;
 351   delete cStore;  cStore = 0;
 352   return trkMuon;
 353 }*/
 354
 355 //-----------------------------------------------------------------------------
 356 Int_t AliMCMuonTrack::SelectHFMuon()
 357 {
 358   Int_t flv = GetParentFlavour(0);
 359   if (flv!=4 && flv!=5) return 2;
 360
 361   Bool_t isRes = kFALSE;
 362   Int_t i=0, nparents=this->GetNParents();
 363   while (i<nparents && !isRes) {
 364     isRes = IsMotherAResonance(i++);
 365   }
 366
 367   if (isRes) return 2;
 368   if (flv==5) return 0;
 369   else return 1;
 370 }
 371
 372 //-----------------------------------------------------------------------------
 373 Bool_t AliMCMuonTrack::IsDiquark(Int_t pdg)
 374 {
 375   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 376   pdg = TMath::Abs(pdg);
 377   if(pdg>1000 && (pdg%100)<10) return kTRUE;
 378   else return kFALSE;
 379 }
 380
 381 //-----------------------------------------------------------------------------
 382 void AliMCMuonTrack::ResetQuarkInfo()
 383 {
 384   for(Int_t pos=1; pos<4; pos++) {
 385     fQuarkIndex[pos] = -1;
 386     fQuarkPDGCode[pos] = 0;
 387   }
 388   return;
 389 }
 390
 391 //-----------------------------------------------------------------------------
 392 Int_t AliMCMuonTrack::GetParentFlavour(Int_t i) const
 393 {
 394   Int_t pdg = GetParentPDGCode(i);
 395   pdg = TMath::Abs(pdg/100);
 396   if(pdg>9) pdg /= 10;
 397   return pdg;
 398 }
 399
 400 //-----------------------------------------------------------------------------
 401 Bool_t AliMCMuonTrack::IsMotherAResonance(Int_t i) const
 402 {
 403   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 404   Int_t pdg = GetParentPDGCode(i);
 405   Int_t id=pdg%100000;
 406   return (!((id-id%10)%110));
 407 }