PWG/muon/AliMuonInfoStoreMC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, 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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // class used to extract and store info of MC particle
  21 //
  22 // Author: X-M. Zhang, zhang@clermont.in2p3.fr
  23 //                     zhangxm@iopp.ccnu.edu.cn
  24 /////////////////////////////////////////////////////////////
  25
  26 #include <TParticle.h>
  27
  28 #include "AliMCEvent.h"
  29 #include "AliAODMCParticle.h"
  30 #include "AliESDMuonTrack.h"
  31 #include "AliAODTrack.h"
  32 #include "AliMuonInfoStoreRD.h"
  33 #include "AliMuonInfoStoreMC.h"
  34 #include "AliMCEvent.h"
  35
  36 ClassImp(AliMuonInfoStoreMC)
  37
  38 const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
  39 const Int_t   AliMuonInfoStoreMC::fgkSourcesN = 6;
  40
  41 //-----------------------------------------------------------------------------
  42 AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
  43 AliMuonInfoStoreRD(),
  44 fIsFull(kFALSE),
  45 fLorentzP(),
  46 fTrackIndex(-1),
  47 fTrackPDGCode(0),
  48 fSource(-1),
  49 fNParents(0),
  50 fOscillation(kFALSE),
  51 fWeight(0.)
  52 {
  53   //
  54   // default constructor
  55   //
  56   for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  57   for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
  58 }
  59
  60 //-----------------------------------------------------------------------------
  61 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack *trkAOD, AliMCEvent *mcEvent, UInt_t selMask, Bool_t full) :
  62 AliMuonInfoStoreRD(trkAOD,selMask),
  63 fIsFull(full),
  64 fLorentzP(),
  65 fTrackIndex(-1),
  66 fTrackPDGCode(0),
  67 fSource(-1),
  68 fNParents(0),
  69 fOscillation(kFALSE),
  70 fWeight(0.)
  71 {
  72   //
  73   // default constructor
  74   //
  75   for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  76   for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
  77
  78   this->SetMCInfoAOD(mcEvent, trkAOD->GetLabel());
  79 }
  80
  81 //-----------------------------------------------------------------------------
  82 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEvent *mcEvent, UInt_t selMask, Bool_t full) :
  83 AliMuonInfoStoreRD(trkESD,selMask),
  84 fIsFull(full),
  85 fLorentzP(),
  86 fTrackIndex(-1),
  87 fTrackPDGCode(0),
  88 fSource(-1),
  89 fNParents(0),
  90 fOscillation(kFALSE),
  91 fWeight(0.)
  92 {
  93   //
  94   // default constructor
  95   //
  96   for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  97   for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
  98
  99   this->SetMCInfoESD(mcEvent, trkESD->GetLabel());
 100 }
 101
 102 //-----------------------------------------------------------------------------
 103 AliMuonInfoStoreMC::AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src) :
 104 AliMuonInfoStoreRD(src),
 105 fIsFull(src.fIsFull),
 106 fLorentzP(src.fLorentzP),
 107 fTrackIndex(src.fTrackIndex),
 108 fTrackPDGCode(src.fTrackPDGCode),
 109 fSource(src.fSource),
 110 fNParents(src.fNParents),
 111 fOscillation(src.fOscillation),
 112 fWeight(src.fWeight)
 113 {
 114   //
 115   // copy constructor
 116   //
 117   for (Int_t i=5; i--;) {
 118     fParentIndex[i]   = src.fParentIndex[i];
 119     fParentPDGCode[i] = src.fParentPDGCode[i];
 120   }
 121   for (Int_t i=4; i--;) {
 122     fQuarkIndex[i]    = src.fQuarkIndex[i];
 123     fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
 124   }
 125 }
 126
 127 //-----------------------------------------------------------------------------
 128 AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
 129 {
 130   //
 131   // assignment constructor
 132   //
 133   if(&src==this) return *this;
 134   AliMuonInfoStoreRD::operator=(src);
 135
 136   fIsFull       = src.fIsFull;
 137   fLorentzP     = src.fLorentzP;
 138   fTrackIndex   = src.fTrackIndex;
 139   fTrackPDGCode = src.fTrackPDGCode;
 140   fSource       = src.fSource;
 141   fNParents     = src.fNParents;
 142   fOscillation  = src.fOscillation;
 143   fWeight       = src.fWeight;
 144   for (Int_t i=5; i--;) {
 145     fParentIndex[i]   = src.fParentIndex[i];
 146     fParentPDGCode[i] = src.fParentPDGCode[i];
 147   }
 148   for (Int_t i=4; i--;) {
 149     fQuarkIndex[i]    = src.fQuarkIndex[i];
 150     fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
 151   }
 152
 153   return *this;
 154 }
 155
 156 //-----------------------------------------------------------------------------
 157 AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
 158 {
 159   //
 160   // destructor
 161   //
 162 }
 163
 164 //-----------------------------------------------------------------------------
 165 void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
 166 {
 167   // fill track MC info with AOD base
 168   fTrackIndex = label;
 169   if (fTrackIndex<0) { fSource=5; return; }
 170
 171   AliAODMCParticle *pMC = (AliAODMCParticle*)mcEvent->GetTrack(fTrackIndex);
 172   fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
 173
 174   fTrackPDGCode = pMC->PdgCode();
 175   if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
 176
 177   Int_t lineM = pMC->GetMother();
 178   if (lineM<0) { fSource=2; return; }
 179
 180   Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPrimary();
 181   if (!isPrimary) { fSource=3; return; }
 182
 183   Int_t countP=-1, pdg=0;
 184   Int_t parents[10], parLine[10];
 185   AliAODMCParticle *mother = 0;
 186   while(lineM>=0){
 187     mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
 188     pdg = mother->GetPdgCode();
 189     if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
 190     parents[++countP] = pdg;
 191     parLine[countP] = lineM;
 192     lineM = mother->GetMother();
 193   }
 194   for(Int_t i=0; i<=countP; i++){
 195     fParentIndex[i] = parLine[countP-i];
 196     fParentPDGCode[i] = parents[countP-i];
 197   }
 198   fNParents = countP + 1;
 199
 200   if (fIsFull && lineM>=0) this->FillHistoryQuarksAOD(mcEvent, lineM);
 201
 202   fSource = this->SelectHFMuon();
 203   return;
 204 }
 205
 206 //-----------------------------------------------------------------------------
 207 void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
 208 {
 209   // fill track MC info with ESD base
 210   fTrackIndex = label;
 211   if (fTrackIndex<0) { fSource=5; return; }
 212
 213   TParticle *pMC = ((AliMCParticle*)mcEvent->GetTrack(fTrackIndex))->Particle();
 214   fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
 215
 216   fTrackPDGCode = pMC->GetPdgCode();
 217   if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
 218
 219   Int_t lineM = pMC->GetFirstMother();
 220   if (lineM<0) { fSource=2; return; }
 221
 222   if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource=3; return; }
 223
 224   Int_t countP=-1, pdg=0;
 225   Int_t parents[10], parLine[10];
 226   TParticle *mother = 0;
 227   while(lineM>=0){
 228     mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
 229     pdg = mother->GetPdgCode();
 230     if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
 231     parents[++countP] = pdg;
 232     parLine[countP] = lineM;
 233     lineM = mother->GetFirstMother();
 234   }
 235   for(Int_t i=0; i<=countP; i++){
 236     fParentIndex[i] = parLine[countP-i];
 237     fParentPDGCode[i] = parents[countP-i];
 238   }
 239   fNParents = countP + 1;
 240
 241   if (fIsFull && lineM>=0) this->FillHistoryQuarksESD(mcEvent, lineM);
 242
 243   fSource = this->SelectHFMuon();
 244   return;
 245 }
 246
 247 //-----------------------------------------------------------------------------
 248 void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t lineM)
 249 {
 250   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 251
 252   if (lineM<0) return;
 253   Int_t countP=-1, pdg=0;
 254   AliAODMCParticle *mother = 0;
 255   while(lineM>=0){
 256     mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
 257     pdg = mother->GetPdgCode();
 258     fQuarkIndex[++countP] = lineM;
 259     fQuarkPDGCode[countP] = pdg;
 260     lineM = mother->GetMother();
 261   }
 262
 263   // for PYTHIA checking
 264   countP = 1;
 265   for(Int_t par=0; par<4; par++) {
 266     if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
 267   }
 268   if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
 269     if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
 270       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 271                  this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
 272
 273       pdg = this->QuarkPDGCode(countP);
 274       Int_t line = this->QuarkIndex(countP);
 275       this->ResetQuarkInfo();
 276       while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
 277         pdg = ((AliAODMCParticle*)mcEvent->GetTrack(++line))->GetPdgCode();
 278       }
 279       while(line>=0){
 280         mother = (AliAODMCParticle*)mcEvent->GetTrack(line);
 281         pdg = mother->GetPdgCode();
 282         fQuarkIndex[countP] = line;
 283         fQuarkPDGCode[countP++] = pdg;
 284         line = mother->GetMother();
 285       }
 286     }
 287   }
 288   return;
 289 }
 290
 291 //-----------------------------------------------------------------------------
 292 void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t lineM)
 293 {
 294   // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 295
 296   if (lineM<0) return;
 297   Int_t countP=-1, pdg=0;
 298   TParticle *mother = 0;
 299   while(lineM>=0){
 300     mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
 301     pdg = mother->GetPdgCode();
 302     fQuarkIndex[++countP] = lineM;
 303     fQuarkPDGCode[countP] = pdg;
 304     lineM = mother->GetFirstMother();
 305   }
 306
 307   // for PYTHIA checking
 308   countP = 1;
 309   for(Int_t par=0; par<4; par++) {
 310     if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
 311   }
 312   if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
 313     if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
 314       AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
 315                  this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
 316
 317       pdg = this->QuarkPDGCode(countP);
 318       Int_t line = this->QuarkIndex(countP);
 319       this->ResetQuarkInfo();
 320       while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
 321         pdg = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle()->GetPdgCode();
 322       }
 323       while(line>=0){
 324         mother = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle();
 325         pdg = mother->GetPdgCode();
 326         fQuarkIndex[countP] = line;
 327         fQuarkPDGCode[countP++] = pdg;
 328         line = mother->GetFirstMother();
 329       }
 330     }
 331   }
 332   return;
 333 }
 334
 335 //-----------------------------------------------------------------------------
 336 Int_t AliMuonInfoStoreMC::SelectHFMuon()
 337 {
 338   // set info of muon from HF
 339
 340   Int_t flv = ParentFlavour(0);
 341   if (flv!=4 && flv!=5) return 2;
 342
 343   Bool_t isRes = kFALSE;
 344   Int_t i=0, nparents=this->ParentsN();
 345   while (i<nparents && !isRes) {
 346     isRes = IsMotherAResonance(i++);
 347   }
 348
 349   if (isRes) return 2;
 350   if (flv==5) return 0;
 351   else return 1;
 352 }
 353
 354 //-----------------------------------------------------------------------------
 355 Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
 356 {
 357   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 358   pdg = TMath::Abs(pdg);
 359   if(pdg>1000 && (pdg%100)<10) return kTRUE;
 360   else return kFALSE;
 361 }
 362
 363 //-----------------------------------------------------------------------------
 364 void AliMuonInfoStoreMC::ResetQuarkInfo()
 365 {
 366   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 367   for(Int_t pos=1; pos<4; pos++) {
 368     fQuarkIndex[pos] = -1;
 369     fQuarkPDGCode[pos] = 0;
 370   }
 371   return;
 372 }
 373
 374 //-----------------------------------------------------------------------------
 375 Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
 376 {
 377   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 378   Int_t pdg = ParentPDGCode(i);
 379   pdg = TMath::Abs(pdg/100);
 380   if(pdg>9) pdg /= 10;
 381   return pdg;
 382 }
 383
 384 //-----------------------------------------------------------------------------
 385 Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
 386 {
 387   // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
 388   Int_t pdg = ParentPDGCode(i);
 389   Int_t id=pdg%100000;
 390   return (!((id-id%10)%110));
 391 }