PWG3/hfe/AliHFEmcQA.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 // QA class of Heavy Flavor quark and fragmeted/decayed particles
  17 // -Check kinematics of Heavy Quarks/hadrons, and decayed leptons
  18 //    pT, rapidity
  19 //    decay lepton kinematics w/wo acceptance
  20 //    heavy hadron decay length, electron pT fraction carried from decay
  21 // -Check yield of Heavy Quarks/hadrons
  22 //    Number of produced heavy quark
  23 //    Number of produced hadron of given pdg code
  24 //
  25 //
  26 // Authors:
  27 //   MinJung Kweon <minjung@physi.uni-heidelberg.de>
  28 //
  29
  30 #include <TH2F.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TParticle.h>
  34
  35 #include <AliLog.h>
  36 #include <AliStack.h>
  37 #include <AliGenEventHeader.h>
  38 #include <AliAODMCParticle.h>
  39
  40 #include "AliHFEmcQA.h"
  41
  42 const Int_t AliHFEmcQA::fgkGluon=21;
  43 const Int_t AliHFEmcQA::fgkMaxGener=10;
  44 const Int_t AliHFEmcQA::fgkMaxIter=100;
  45 const Int_t AliHFEmcQA::fgkqType = 7;    // number of species waiting for QA done
  46
  47 ClassImp(AliHFEmcQA)
  48
  49 //_______________________________________________________________________________________________
  50 AliHFEmcQA::AliHFEmcQA() :
  51         fStack(0x0)
  52         ,fMCHeader(0x0)
  53         ,fMCArray(0x0)
  54         ,fNparents(0)
  55 {
  56         // Default constructor
  57
  58 }
  59
  60 //_______________________________________________________________________________________________
  61 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
  62         TObject(p)
  63         ,fStack(0x0)
  64         ,fMCHeader(0x0)
  65         ,fMCArray(0x0)
  66         ,fNparents(p.fNparents)
  67 {
  68         // Copy constructor
  69 }
  70
  71 //_______________________________________________________________________________________________
  72 AliHFEmcQA&
  73 AliHFEmcQA::operator=(const AliHFEmcQA &)
  74 {
  75   // Assignment operator
  76
  77   AliInfo("Not yet implemented.");
  78   return *this;
  79 }
  80
  81 //_______________________________________________________________________________________________
  82 AliHFEmcQA::~AliHFEmcQA()
  83 {
  84         // Destructor
  85
  86         AliInfo("Analysis Done.");
  87 }
  88
  89 //_______________________________________________________________________________________________
  90 void AliHFEmcQA::PostAnalyze() const
  91 {
  92 }
  93
  94 //__________________________________________
  95 void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
  96 {
  97   // create histograms
  98
  99   if (kquark != kCharm && kquark != kBeauty) {
 100     AliDebug(1, "This task is only for heavy quark QA, return\n");
 101     return;
 102   }
 103   Int_t iq = kquark - kCharm;
 104
 105   TString kqTypeLabel[fgkqType];
 106   if (kquark == kCharm){
 107     kqTypeLabel[kQuark]="c";
 108     kqTypeLabel[kantiQuark]="cbar";
 109     kqTypeLabel[kHadron]="cHadron";
 110     kqTypeLabel[keHadron]="ceHadron";
 111     kqTypeLabel[kDeHadron]="nullHadron";
 112     kqTypeLabel[kElectron]="ce";
 113     kqTypeLabel[kElectron2nd]="nulle";
 114   } else if (kquark == kBeauty){
 115     kqTypeLabel[kQuark]="b";
 116     kqTypeLabel[kantiQuark]="bbar";
 117     kqTypeLabel[kHadron]="bHadron";
 118     kqTypeLabel[keHadron]="beHadron";
 119     kqTypeLabel[kDeHadron]="bDeHadron";
 120     kqTypeLabel[kElectron]="be";
 121     kqTypeLabel[kElectron2nd]="bce";
 122   }
 123
 124
 125   TString hname;
 126   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
 127      if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
 128      hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
 129      fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
 130      hname = hnopt+"Pt_"+kqTypeLabel[iqType];
 131      fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",250,0,50);
 132      hname = hnopt+"Y_"+kqTypeLabel[iqType];
 133      fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
 134      hname = hnopt+"Eta_"+kqTypeLabel[iqType];
 135      fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
 136   }
 137
 138   if (icut == 0){
 139     hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
 140     fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
 141     hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
 142     fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
 143   }
 144   hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
 145   fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
 146   hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
 147   fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
 148   hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
 149   fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
 150   hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
 151   fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
 152
 153 }
 154
 155 //__________________________________________
 156 void AliHFEmcQA::Init()
 157 {
 158   // called at begining every event
 159
 160   for (Int_t i=0; i<2; i++){
 161      fIsHeavy[i] = 0;
 162   }
 163
 164   fNparents = 7;
 165
 166   fParentSelect[0][0] =  411;
 167   fParentSelect[0][1] =  421;
 168   fParentSelect[0][2] =  431;
 169   fParentSelect[0][3] = 4122;
 170   fParentSelect[0][4] = 4132;
 171   fParentSelect[0][5] = 4232;
 172   fParentSelect[0][6] = 4332;
 173
 174   fParentSelect[1][0] =  511;
 175   fParentSelect[1][1] =  521;
 176   fParentSelect[1][2] =  531;
 177   fParentSelect[1][3] = 5122;
 178   fParentSelect[1][4] = 5132;
 179   fParentSelect[1][5] = 5232;
 180   fParentSelect[1][6] = 5332;
 181
 182 }
 183
 184 //__________________________________________
 185 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
 186 {
 187   // get heavy quark kinematics
 188
 189     if (kquark != kCharm && kquark != kBeauty) {
 190       AliDebug(1, "This task is only for heavy quark QA, return\n");
 191       return;
 192     }
 193     Int_t iq = kquark - kCharm;
 194
 195     if (iTrack < 0) {
 196       AliDebug(1, "Stack label is negative, return\n");
 197       return;
 198     }
 199
 200     //TParticle *part = fStack->Particle(iTrack);
 201     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 202
 203     // select heavy hadron or not fragmented heavy quark
 204     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
 205
 206       TParticle *partMother;
 207       Int_t iLabel;
 208
 209       if (partPdgcode == kquark){ // in case of not fragmented heavy quark
 210         partMother = part;
 211         iLabel = iTrack;
 212       } else{ // in case of heavy hadron, start to search for mother heavy parton
 213         iLabel = part->GetFirstMother();
 214         if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
 215         else {
 216           AliDebug(1, "Stack label is negative, return\n");
 217           return;
 218         }
 219       }
 220
 221       // heavy parton selection as a mother of heavy hadron
 222       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
 223       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
 224       // should I make a condition that partMother should be quark or diquark? -> not necessary
 225       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
 226       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
 227
 228         if ( abs(partMother->GetPdgCode()) != kquark ){
 229           // search fragmented partons in the same string
 230           Bool_t isSameString = kTRUE;
 231           for (Int_t i=1; i<fgkMaxIter; i++){
 232              iLabel = iLabel - 1;
 233              if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
 234              else {
 235                AliDebug(1, "Stack label is negative, return\n");
 236                return;
 237              }
 238              if ( abs(partMother->GetPdgCode()) == kquark ) break;
 239              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
 240              if (!isSameString) return;
 241           }
 242         }
 243         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
 244         if (abs(partMother->GetPdgCode()) != kquark) return;
 245
 246         if (fIsHeavy[iq] >= 50) return;
 247         fHeavyQuark[fIsHeavy[iq]] = partMother;
 248         fIsHeavy[iq]++;
 249
 250         // fill kinematics for heavy parton
 251         if (partMother->GetPdgCode() > 0) { // quark
 252           fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 253           fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
 254           fHist[iq][kQuark][0].fY->Fill(GetRapidity(partMother));
 255           fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
 256         } else{ // antiquark
 257           fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 258           fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
 259           fHist[iq][kantiQuark][0].fY->Fill(GetRapidity(partMother));
 260           fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
 261         }
 262
 263       } // end of heavy parton slection loop
 264
 265     } // end of heavy hadron or quark selection
 266
 267 }
 268
 269 //__________________________________________
 270 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
 271 {
 272   // end of event analysis
 273
 274   if (kquark != kCharm && kquark != kBeauty) {
 275     AliDebug(1, "This task is only for heavy quark QA, return\n");
 276     return;
 277   }
 278   Int_t iq = kquark - kCharm;
 279
 280
 281   // # of heavy quark per event
 282   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
 283   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
 284
 285   Int_t motherID[fgkMaxGener];
 286   Int_t motherType[fgkMaxGener];
 287   Int_t motherLabel[fgkMaxGener];
 288   Int_t ancestorPdg[fgkMaxGener];
 289   Int_t ancestorLabel[fgkMaxGener];
 290
 291   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
 292      motherID[i] = 0;
 293      motherType[i] = 0;
 294      motherLabel[i] = 0;
 295      ancestorPdg[i] = 0;
 296      ancestorLabel[i] = 0;
 297   }
 298
 299   // check history of found heavy quarks
 300   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 301
 302      ancestorLabel[0] = i;
 303      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
 304      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
 305
 306      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
 307      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
 308
 309      Int_t ig = 1;
 310      while (ancestorLabel[ig] != -1){
 311           // in case there is mother, get mother's pdg code and grandmother's label
 312           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
 313           // if mother is still heavy, find again mother's ancestor
 314           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
 315             ig++;
 316             continue; // if it is from same heavy
 317           }
 318           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
 319           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 320           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 321           // if it is not the above case, something is strange
 322           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
 323           break;
 324      }
 325      if (ancestorLabel[ig] == -1){ // from hard scattering
 326        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
 327      }
 328
 329   } // end of found heavy quark loop
 330
 331
 332   // check process type
 333   Int_t processID = 0;
 334   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 335      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
 336   }
 337
 338
 339   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
 340   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
 341
 342      Int_t id1 = ipair*2;
 343      Int_t id2 = ipair*2 + 1;
 344
 345      if (motherType[id1] == 2 && motherType[id2] == 2){
 346        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
 347        else processID = -9;
 348      }
 349      else if (motherType[id1] == -1 && motherType[id2] == -1) {
 350        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
 351          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
 352          else processID = kPairCreationFromq; // q-qbar pair creation
 353        }
 354        else processID = -8;
 355      }
 356      else if (motherType[id1] == -1 || motherType[id2] == -1) {
 357        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
 358          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
 359          else processID = kLightQuarkShower;
 360        }
 361        else processID = -7;
 362      }
 363      else if (motherType[id1] == -2 || motherType[id2] == -2) {
 364        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
 365        else processID = -6;
 366
 367      }
 368      else processID = -5;
 369
 370      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
 371      else fHistComm[iq][0].fProcessID->Fill(processID);
 372      AliDebug(1,Form("Process ID = %d\n",processID));
 373   } // end of # heavy quark pair loop
 374
 375 }
 376
 377 //__________________________________________
 378 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
 379 {
 380     // decay electron kinematics
 381
 382     if (kquark != kCharm && kquark != kBeauty) {
 383       AliDebug(1, "This task is only for heavy quark QA, return\n");
 384       return;
 385     }
 386     Int_t iq = kquark - kCharm;
 387
 388     //TParticle* mcpart = fStack->Particle(iTrack);
 389
 390     Int_t iLabel = mcpart->GetFirstMother();
 391     if (iLabel<0){
 392       AliDebug(1, "Stack label is negative, return\n");
 393       return;
 394     }
 395
 396     TParticle *partCopy = mcpart;
 397     Int_t pdgcode = mcpart->GetPdgCode();
 398     Int_t pdgcodeCopy = pdgcode;
 399
 400     // if the mother is charmed hadron
 401     Bool_t isDirectCharm = kFALSE;
 402     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
 403
 404           // iterate until you find B hadron as a mother or become top ancester
 405           for (Int_t i=1; i<fgkMaxIter; i++){
 406
 407              Int_t jLabel = mcpart->GetFirstMother();
 408              if (jLabel == -1){
 409                isDirectCharm = kTRUE;
 410                break; // if there is no ancester
 411              }
 412              if (jLabel < 0){ // safety protection
 413                AliDebug(1, "Stack label is negative, return\n");
 414                return;
 415              }
 416              // if there is an ancester
 417              TParticle* mother = fStack->Particle(jLabel);
 418              Int_t motherPDG = mother->GetPdgCode();
 419
 420              for (Int_t j=0; j<fNparents; j++){
 421                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
 422              }
 423
 424              mcpart = mother;
 425           } // end of iteration
 426     } // end of if
 427     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 428          for (Int_t i=0; i<fNparents; i++){
 429             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
 430
 431               // fill hadron kinematics
 432               fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
 433               fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
 434               fHist[iq][kHadron][0].fY->Fill(GetRapidity(partCopy));
 435               fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
 436
 437             }
 438          }
 439     } // end of if
 440 }
 441
 442 //__________________________________________
 443 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
 444 {
 445     // decay electron kinematics
 446
 447     if (kquark != kCharm && kquark != kBeauty) {
 448       AliDebug(1, "This task is only for heavy quark QA, return\n");
 449       return;
 450     }
 451     Int_t iq = kquark - kCharm;
 452     Bool_t isFinalOpenCharm = kFALSE;
 453
 454 /*
 455     if (iTrack < 0) {
 456       AliDebug(1, "Stack label is negative, return\n");
 457       return;
 458     }
 459     */
 460
 461     //TParticle* mcpart = fStack->Particle(iTrack);
 462
 463     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 464
 465     Int_t iLabel = mcpart->GetFirstMother();
 466     if (iLabel<0){
 467       AliDebug(1, "Stack label is negative, return\n");
 468       return;
 469     }
 470
 471     TParticle *partMother = fStack->Particle(iLabel);
 472     TParticle *partMotherCopy = partMother;
 473     Int_t maPdgcode = partMother->GetPdgCode();
 474     Int_t maPdgcodeCopy = maPdgcode;
 475
 476     // get mc primary vertex
 477     TArrayF mcPrimVtx;
 478     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
 479
 480     // get electron production vertex
 481     TLorentzVector ePoint;
 482     mcpart->ProductionVertex(ePoint);
 483
 484     // calculated production vertex to primary vertex (in xy plane)
 485     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
 486
 487     // if the mother is charmed hadron
 488     Bool_t isMotherDirectCharm = kFALSE;
 489     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 490
 491          for (Int_t i=0; i<fNparents; i++){
 492             if (abs(maPdgcode)==fParentSelect[0][i]){
 493               isFinalOpenCharm = kTRUE;
 494             }
 495          }
 496          if (!isFinalOpenCharm) return ;
 497
 498           // iterate until you find B hadron as a mother or become top ancester
 499           for (Int_t i=1; i<fgkMaxIter; i++){
 500
 501              Int_t jLabel = partMother->GetFirstMother();
 502              if (jLabel == -1){
 503                isMotherDirectCharm = kTRUE;
 504                break; // if there is no ancester
 505              }
 506              if (jLabel < 0){ // safety protection
 507                AliDebug(1, "Stack label is negative, return\n");
 508                return;
 509              }
 510
 511              // if there is an ancester
 512              TParticle* grandMa = fStack->Particle(jLabel);
 513              Int_t grandMaPDG = grandMa->GetPdgCode();
 514
 515              for (Int_t j=0; j<fNparents; j++){
 516                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 517
 518                   if (kquark == kCharm) return;
 519                   // fill electron kinematics
 520                   fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 521                   fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
 522                   fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart));
 523                   fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
 524
 525                   // fill mother hadron kinematics
 526                   fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
 527                   fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
 528                   fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa));
 529                   fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
 530
 531                   // ratio between pT of electron and pT of mother B hadron
 532                   if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 533
 534                   // distance between electron production point and primary vertex
 535                   fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
 536                   return;
 537                 }
 538              }
 539
 540              partMother = grandMa;
 541           } // end of iteration
 542     } // end of if
 543     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 544          for (Int_t i=0; i<fNparents; i++){
 545             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 546
 547               // fill electron kinematics
 548               fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 549               fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 550               fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart));
 551               fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 552
 553               // fill mother hadron kinematics
 554               fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
 555               fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
 556               fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy));
 557               fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
 558
 559               // ratio between pT of electron and pT of mother B hadron
 560               if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
 561
 562               // distance between electron production point and primary vertex
 563               fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
 564             }
 565          }
 566     } // end of if
 567 }
 568
 569 //____________________________________________________________________
 570 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
 571 {
 572   // decay electron kinematics
 573
 574   if (kquark != kCharm && kquark != kBeauty) {
 575     AliDebug(1, "This task is only for heavy quark QA, return\n");
 576     return;
 577   }
 578
 579   Int_t iq = kquark - kCharm;
 580   Bool_t isFinalOpenCharm = kFALSE;
 581
 582   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 583
 584   // mother
 585   Int_t iLabel = mcpart->GetMother();
 586   if (iLabel<0){
 587     AliDebug(1, "Stack label is negative, return\n");
 588     return;
 589   }
 590
 591   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
 592   AliAODMCParticle *partMotherCopy = partMother;
 593   Int_t maPdgcode = partMother->GetPdgCode();
 594   Int_t maPdgcodeCopy = maPdgcode;
 595
 596   Bool_t isMotherDirectCharm = kFALSE;
 597   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 598
 599     for (Int_t i=0; i<fNparents; i++){
 600        if (abs(maPdgcode)==fParentSelect[0][i]){
 601          isFinalOpenCharm = kTRUE;
 602        }
 603     }
 604     if (!isFinalOpenCharm) return;
 605
 606     for (Int_t i=1; i<fgkMaxIter; i++){
 607
 608        Int_t jLabel = partMother->GetMother();
 609        if (jLabel == -1){
 610          isMotherDirectCharm = kTRUE;
 611          break; // if there is no ancester
 612        }
 613        if (jLabel < 0){ // safety protection
 614          AliDebug(1, "Stack label is negative, return\n");
 615          return;
 616        }
 617
 618        // if there is an ancester
 619        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
 620        Int_t grandMaPDG = grandMa->GetPdgCode();
 621
 622        for (Int_t j=0; j<fNparents; j++){
 623           if (abs(grandMaPDG)==fParentSelect[1][j]){
 624
 625             if (kquark == kCharm) return;
 626             // fill electron kinematics
 627             fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 628             fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
 629             fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart));
 630             fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
 631
 632             // fill mother hadron kinematics
 633             fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
 634             fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
 635             fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa));
 636             fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
 637
 638             return;
 639           }
 640        }
 641
 642        partMother = grandMa;
 643     } // end of iteration
 644   } // end of if
 645   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 646     for (Int_t i=0; i<fNparents; i++){
 647        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 648
 649          // fill electron kinematics
 650          fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 651          fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 652          fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart));
 653          fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 654
 655          // fill mother hadron kinematics
 656          fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
 657          fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
 658          fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy));
 659          fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
 660
 661        }
 662     }
 663   } // end of if
 664
 665 }
 666
 667 //__________________________________________
 668 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
 669 {
 670        // find mother pdg code and label
 671
 672        if (motherlabel < 0) {
 673          AliDebug(1, "Stack label is negative, return\n");
 674          return;
 675        }
 676        TParticle *heavysMother = fStack->Particle(motherlabel);
 677        motherpdg = heavysMother->GetPdgCode();
 678        grandmotherlabel = heavysMother->GetFirstMother();
 679        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
 680 }
 681
 682 //__________________________________________
 683 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 684 {
 685        // mothertype -1 means this heavy quark coming from hard vertex
 686
 687        TParticle *afterinitialrad1  = fStack->Particle(4);
 688        TParticle *afterinitialrad2  = fStack->Particle(5);
 689
 690        motherlabel = -1;
 691
 692        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
 693          AliDebug(1,"heavy from gluon gluon pair creation!\n");
 694          mothertype = -1;
 695          motherID = fgkGluon;
 696        }
 697        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
 698          AliDebug(1,"heavy from flavor exitation!\n");
 699          mothertype = -1;
 700          motherID = kquark;
 701        }
 702        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
 703          AliDebug(1,"heavy from q-qbar pair creation!\n");
 704          mothertype = -1;
 705          motherID = 1;
 706        }
 707        else {
 708          AliDebug(1,"something strange!\n");
 709          mothertype = -999;
 710          motherlabel = -999;
 711          motherID = -999;
 712        }
 713 }
 714
 715 //__________________________________________
 716 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 717 {
 718        // mothertype -2 means this heavy quark coming from initial state
 719
 720        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
 721          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
 722          motherID = heavysMother->GetPdgCode();
 723          mothertype = -2; // there is mother before initial state radiation
 724          motherlabel = inputmotherlabel;
 725          AliDebug(1,"initial parton shower! \n");
 726
 727          return kTRUE;
 728        }
 729
 730        return kFALSE;
 731 }
 732
 733 //__________________________________________
 734 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 735 {
 736        // mothertype 2 means this heavy quark coming from final state
 737
 738        if (inputmotherlabel > 5){ // mother exist after hard scattering
 739          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
 740          motherID = heavysMother->GetPdgCode();
 741          mothertype = 2; //
 742          motherlabel = inputmotherlabel;
 743          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
 744
 745          return kTRUE;
 746        }
 747        return kFALSE;
 748 }
 749
 750 //__________________________________________
 751 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 752 {
 753       // mark strange behavior
 754
 755        mothertype = -888;
 756        motherlabel = -888;
 757        motherID = -888;
 758        AliDebug(1,"something strange!\n");
 759 }
 760
 761 //__________________________________________
 762 Float_t AliHFEmcQA::GetRapidity(TParticle *part) const
 763 {
 764       // return rapidity
 765
 766        Float_t rapidity;
 767        if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
 768        else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 769        return rapidity;
 770 }
 771
 772 //__________________________________________
 773 Float_t AliHFEmcQA::GetRapidity(AliAODMCParticle *part) const
 774 {
 775       // return rapidity
 776
 777        Float_t rapidity;
 778        if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999;
 779        else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz())));
 780        return rapidity;
 781 }
 782
 783 //__________________________________________
 784 Int_t AliHFEmcQA::GetSource(AliAODMCParticle * const mcpart)
 785 {
 786   // decay particle's origin
 787
 788   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
 789
 790   Int_t origin = -1;
 791   Bool_t isFinalOpenCharm = kFALSE;
 792
 793   // mother
 794   Int_t iLabel = mcpart->GetMother();
 795   if (iLabel<0){
 796     AliDebug(1, "Stack label is negative, return\n");
 797     return -1;
 798   }
 799
 800   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
 801   Int_t maPdgcode = partMother->GetPdgCode();
 802
 803   // if the mother is charmed hadron
 804   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 805
 806     for (Int_t i=0; i<fNparents; i++){
 807        if (abs(maPdgcode)==fParentSelect[0][i]){
 808          isFinalOpenCharm = kTRUE;
 809        }
 810     }
 811     if (!isFinalOpenCharm) return -1;
 812
 813     // iterate until you find B hadron as a mother or become top ancester
 814     for (Int_t i=1; i<fgkMaxIter; i++){
 815
 816        Int_t jLabel = partMother->GetMother();
 817        if (jLabel == -1){
 818          origin = kDirectCharm;
 819          return origin;
 820        }
 821        if (jLabel < 0){ // safety protection
 822          AliDebug(1, "Stack label is negative, return\n");
 823          return -1;
 824        }
 825
 826        // if there is an ancester
 827        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
 828        Int_t grandMaPDG = grandMa->GetPdgCode();
 829
 830        for (Int_t j=0; j<fNparents; j++){
 831           if (abs(grandMaPDG)==fParentSelect[1][j]){
 832             origin = kBeautyCharm;
 833             return origin;
 834           }
 835        }
 836
 837        partMother = grandMa;
 838     } // end of iteration
 839   } // end of if
 840   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
 841     for (Int_t i=0; i<fNparents; i++){
 842        if (abs(maPdgcode)==fParentSelect[1][i]){
 843          origin = kDirectBeauty;
 844          return origin;
 845        }
 846     }
 847   } // end of if
 848   else if ( abs(maPdgcode) == 22 ) {
 849     origin = kGamma;
 850     return origin;
 851   } // end of if
 852   else if ( abs(maPdgcode) == 111 ) {
 853     origin = kPi0;
 854     return origin;
 855   } // end of if
 856
 857   return origin;
 858 }
 859
 860 //__________________________________________
 861 Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
 862 {
 863   // decay particle's origin
 864
 865   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
 866
 867   Int_t origin = -1;
 868   Bool_t isFinalOpenCharm = kFALSE;
 869
 870   Int_t iLabel = mcpart->GetFirstMother();
 871   if (iLabel<0){
 872     AliDebug(1, "Stack label is negative, return\n");
 873     return -1;
 874   }
 875
 876   TParticle *partMother = fStack->Particle(iLabel);
 877   Int_t maPdgcode = partMother->GetPdgCode();
 878
 879    // if the mother is charmed hadron
 880    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 881
 882      for (Int_t i=0; i<fNparents; i++){
 883         if (abs(maPdgcode)==fParentSelect[0][i]){
 884           isFinalOpenCharm = kTRUE;
 885         }
 886      }
 887      if (!isFinalOpenCharm) return -1;
 888
 889      // iterate until you find B hadron as a mother or become top ancester
 890      for (Int_t i=1; i<fgkMaxIter; i++){
 891
 892         Int_t jLabel = partMother->GetFirstMother();
 893         if (jLabel == -1){
 894           origin = kDirectCharm;
 895           return origin;
 896         }
 897         if (jLabel < 0){ // safety protection
 898           AliDebug(1, "Stack label is negative, return\n");
 899           return -1;
 900         }
 901
 902         // if there is an ancester
 903         TParticle* grandMa = fStack->Particle(jLabel);
 904         Int_t grandMaPDG = grandMa->GetPdgCode();
 905
 906         for (Int_t j=0; j<fNparents; j++){
 907            if (abs(grandMaPDG)==fParentSelect[1][j]){
 908              origin = kBeautyCharm;
 909              return origin;
 910            }
 911         }
 912
 913         partMother = grandMa;
 914      } // end of iteration
 915    } // end of if
 916    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
 917      for (Int_t i=0; i<fNparents; i++){
 918         if (abs(maPdgcode)==fParentSelect[1][i]){
 919           origin = kDirectBeauty;
 920           return origin;
 921         }
 922      }
 923    } // end of if
 924    else if ( abs(maPdgcode) == 22 ) {
 925      origin = kGamma;
 926      return origin;
 927    } // end of if
 928    else if ( abs(maPdgcode) == 111 ) {
 929      origin = kPi0;
 930      return origin;
 931    } // end of if
 932    else origin = kElse;
 933
 934    return origin;
 935 }