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  *                                                                        *
  17  * QA class of Heavy Flavor quark and fragmeted/decayed particles         *
  18  * -Check kinematics of Heavy Quarks/hadrons, and decayed leptons         *
  19  *    pT, rapidity                                                        *
  20  *    decay lepton kinematics w/wo acceptance                             *
  21  *    heavy hadron decay length, electron pT fraction carried from decay  *
  22  * -Check yield of Heavy Quarks/hadrons                                   *
  23  *    Number of produced heavy quark                                      *
  24  *    Number of produced hadron of given pdg code                         *
  25  *                                                                        *
  26  *                                                                        *
  27  * Authors:                                                               *
  28  *   MinJung Kweon <minjung@physi.uni-heidelberg.de>                      *
  29  *                                                                        *
  30  **************************************************************************/
  31
  32
  33 #include <iostream>
  34 #include <TH2F.h>
  35 #include <TCanvas.h>
  36 #include <TFile.h>
  37 #include <TH1F.h>
  38 #include <TH2F.h>
  39 #include <TCut.h>
  40 #include <TRandom.h>
  41 #include <TDatabasePDG.h>
  42 #include <TROOT.h>
  43 #include <TParticle.h>
  44
  45 #include <AliLog.h>
  46 #include <AliESDEvent.h>
  47 #include <AliESDtrack.h>
  48 #include <AliESDInputHandler.h>
  49 #include <AliESDVertex.h>
  50 #include <AliStack.h>
  51
  52 #include "AliHFEmcQA.h"
  53
  54 const Int_t AliHFEmcQA::fgkGluon=21;
  55 const Int_t AliHFEmcQA::fgkMaxGener=10;
  56 const Int_t AliHFEmcQA::fgkMaxIter=100;
  57 const Int_t AliHFEmcQA::fgkqType = 7;    // number of species waiting for QA done
  58
  59 ClassImp(AliHFEmcQA)
  60
  61 //_______________________________________________________________________________________________
  62 AliHFEmcQA::AliHFEmcQA() :
  63         fStack(0x0)
  64         ,fNparents(0)
  65 {
  66         // Default constructor
  67
  68 }
  69
  70 //_______________________________________________________________________________________________
  71 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
  72         TObject(p)
  73         ,fStack(0x0)
  74         ,fNparents(p.fNparents)
  75 {
  76         // Copy constructor
  77 }
  78
  79 //_______________________________________________________________________________________________
  80 AliHFEmcQA&
  81 AliHFEmcQA::operator=(const AliHFEmcQA &)
  82 {
  83   // Assignment operator
  84
  85   AliInfo("Not yet implemented.");
  86   return *this;
  87 }
  88
  89 //_______________________________________________________________________________________________
  90 AliHFEmcQA::~AliHFEmcQA()
  91 {
  92         // Destructor
  93
  94         AliInfo("Analysis Done.");
  95 }
  96
  97 //_______________________________________________________________________________________________
  98 void AliHFEmcQA::PostAnalyze()
  99 {
 100 }
 101
 102 //__________________________________________
 103 void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
 104 {
 105   // create histograms
 106
 107   if (kquark != kCharm && kquark != kBeauty) {
 108     AliDebug(1, "This task is only for heavy quark QA, return\n");
 109     return;
 110   }
 111   Int_t iq = kquark - kCharm;
 112
 113   TString kqTypeLabel[fgkqType];
 114   if (kquark == kCharm){
 115     kqTypeLabel[kQuark]="c";
 116     kqTypeLabel[kantiQuark]="cbar";
 117     kqTypeLabel[kHadron]="cHadron";
 118     kqTypeLabel[keHadron]="ceHadron";
 119     kqTypeLabel[kDeHadron]="nullHadron";
 120     kqTypeLabel[kElectron]="ce";
 121     kqTypeLabel[kElectron2nd]="nulle";
 122   } else if (kquark == kBeauty){
 123     kqTypeLabel[kQuark]="b";
 124     kqTypeLabel[kantiQuark]="bbar";
 125     kqTypeLabel[kHadron]="bHadron";
 126     kqTypeLabel[keHadron]="beHadron";
 127     kqTypeLabel[kDeHadron]="bDeHadron";
 128     kqTypeLabel[kElectron]="be";
 129     kqTypeLabel[kElectron2nd]="bce";
 130   }
 131
 132
 133   TString hname;
 134   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
 135      if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
 136      hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
 137      fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
 138      hname = hnopt+"Pt_"+kqTypeLabel[iqType];
 139      fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",250,0,50);
 140      hname = hnopt+"Y_"+kqTypeLabel[iqType];
 141      fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
 142      hname = hnopt+"Eta_"+kqTypeLabel[iqType];
 143      fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
 144   }
 145
 146   if (icut == 0){
 147     hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
 148     fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
 149     hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
 150     fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
 151   }
 152   hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
 153   fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
 154   hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
 155   fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
 156   hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
 157   fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
 158   hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
 159   fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
 160
 161 }
 162
 163 //__________________________________________
 164 void AliHFEmcQA::Init()
 165 {
 166   // called at begining every event
 167
 168   for (Int_t i=0; i<2; i++){
 169      fIsHeavy[i] = 0;
 170   }
 171
 172   fNparents = 7;
 173
 174   fParentSelect[0][0] =  411;
 175   fParentSelect[0][1] =  421;
 176   fParentSelect[0][2] =  431;
 177   fParentSelect[0][3] = 4122;
 178   fParentSelect[0][4] = 4132;
 179   fParentSelect[0][5] = 4232;
 180   fParentSelect[0][6] = 4332;
 181
 182   fParentSelect[1][0] =  511;
 183   fParentSelect[1][1] =  521;
 184   fParentSelect[1][2] =  531;
 185   fParentSelect[1][3] = 5122;
 186   fParentSelect[1][4] = 5132;
 187   fParentSelect[1][5] = 5232;
 188   fParentSelect[1][6] = 5332;
 189
 190 }
 191
 192 //__________________________________________
 193 void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark)
 194 {
 195   // get heavy quark kinematics
 196
 197     if (kquark != kCharm && kquark != kBeauty) {
 198       AliDebug(1, "This task is only for heavy quark QA, return\n");
 199       return;
 200     }
 201     Int_t iq = kquark - kCharm;
 202
 203
 204     if (iTrack < 0) {
 205       AliDebug(1, "Stack label is negative, return\n");
 206       return;
 207     }
 208
 209     TParticle *part = fStack->Particle(iTrack);
 210     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 211
 212     // select heavy hadron or not fragmented heavy quark
 213     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
 214
 215       TParticle *partMother;
 216       Int_t iLabel;
 217
 218       if (partPdgcode == kquark){ // in case of not fragmented heavy quark
 219         partMother = part;
 220         iLabel = iTrack;
 221       } else{ // in case of heavy hadron, start to search for mother heavy parton
 222         iLabel = part->GetFirstMother();
 223         if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
 224         else {
 225           AliDebug(1, "Stack label is negative, return\n");
 226           return;
 227         }
 228       }
 229
 230       // heavy parton selection as a mother of heavy hadron
 231       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
 232       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
 233       // should I make a condition that partMother should be quark or diquark? -> not necessary
 234       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
 235       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
 236
 237         if ( abs(partMother->GetPdgCode()) != kquark ){
 238           // search fragmented partons in the same string
 239           Bool_t isSameString = kTRUE;
 240           for (Int_t i=1; i<fgkMaxIter; i++){
 241              iLabel = iLabel - 1;
 242              if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
 243              else {
 244                AliDebug(1, "Stack label is negative, return\n");
 245                return;
 246              }
 247              if ( abs(partMother->GetPdgCode()) == kquark ) break;
 248              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
 249              if (!isSameString) return;
 250           }
 251         }
 252         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
 253         if (abs(partMother->GetPdgCode()) != kquark) return;
 254
 255         if (fIsHeavy[iq] >= 50) return;
 256         fHeavyQuark[fIsHeavy[iq]] = partMother;
 257         fIsHeavy[iq]++;
 258
 259         // fill kinematics for heavy parton
 260         if (partMother->GetPdgCode() > 0) { // quark
 261           fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 262           fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
 263           fHist[iq][kQuark][0].fY->Fill(GetRapidity(partMother));
 264           fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
 265         } else{ // antiquark
 266           fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 267           fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
 268           fHist[iq][kantiQuark][0].fY->Fill(GetRapidity(partMother));
 269           fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
 270         }
 271
 272       } // end of heavy parton slection loop
 273
 274     } // end of heavy hadron or quark selection
 275
 276 }
 277
 278 //__________________________________________
 279 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
 280 {
 281   // end of event analysis
 282
 283   if (kquark != kCharm && kquark != kBeauty) {
 284     AliDebug(1, "This task is only for heavy quark QA, return\n");
 285     return;
 286   }
 287   Int_t iq = kquark - kCharm;
 288
 289
 290   // # of heavy quark per event
 291   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
 292   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
 293
 294   Int_t motherID[fgkMaxGener];
 295   Int_t motherType[fgkMaxGener];
 296   Int_t motherLabel[fgkMaxGener];
 297   Int_t ancestorPdg[fgkMaxGener];
 298   Int_t ancestorLabel[fgkMaxGener];
 299
 300   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
 301      motherID[i] = 0;
 302      motherType[i] = 0;
 303      motherLabel[i] = 0;
 304      ancestorPdg[i] = 0;
 305      ancestorLabel[i] = 0;
 306   }
 307
 308   // check history of found heavy quarks
 309   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 310
 311      ancestorLabel[0] = i;
 312      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
 313      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
 314
 315      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
 316      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
 317
 318      Int_t ig = 1;
 319      while (ancestorLabel[ig] != -1){
 320           // in case there is mother, get mother's pdg code and grandmother's label
 321           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
 322           // if mother is still heavy, find again mother's ancestor
 323           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
 324             ig++;
 325             continue; // if it is from same heavy
 326           }
 327           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
 328           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 329           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 330           // if it is not the above case, something is strange
 331           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
 332           break;
 333      }
 334      if (ancestorLabel[ig] == -1){ // from hard scattering
 335        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
 336      }
 337
 338   } // end of found heavy quark loop
 339
 340
 341   // check process type
 342   Int_t processID = 0;
 343   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 344      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
 345   }
 346
 347
 348   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
 349   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
 350
 351      Int_t id1 = ipair*2;
 352      Int_t id2 = ipair*2 + 1;
 353
 354      if (motherType[id1] == 2 && motherType[id2] == 2){
 355        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
 356        else processID = -9;
 357      }
 358      else if (motherType[id1] == -1 && motherType[id2] == -1) {
 359        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
 360          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
 361          else processID = kPairCreationFromq; // q-qbar pair creation
 362        }
 363        else processID = -8;
 364      }
 365      else if (motherType[id1] == -1 || motherType[id2] == -1) {
 366        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
 367          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
 368          else processID = kLightQuarkShower;
 369        }
 370        else processID = -7;
 371      }
 372      else if (motherType[id1] == -2 || motherType[id2] == -2) {
 373        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
 374        else processID = -6;
 375
 376      }
 377      else processID = -5;
 378
 379      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
 380      else fHistComm[iq][0].fProcessID->Fill(processID);
 381      AliDebug(1,Form("Process ID = %d\n",processID));
 382   } // end of # heavy quark pair loop
 383
 384 }
 385
 386 //__________________________________________
 387 void AliHFEmcQA::GetHadronKine(Int_t iTrack, const Int_t kquark)
 388 {
 389     // decay electron kinematics
 390
 391     if (kquark != kCharm && kquark != kBeauty) {
 392       AliDebug(1, "This task is only for heavy quark QA, return\n");
 393       return;
 394     }
 395     Int_t iq = kquark - kCharm;
 396
 397     if (iTrack < 0) {
 398       AliDebug(1, "Stack label is negative, return\n");
 399       return;
 400     }
 401
 402     TParticle* mcpart = fStack->Particle(iTrack);
 403
 404     Int_t iLabel = mcpart->GetFirstMother();
 405     if (iLabel<0){
 406       AliDebug(1, "Stack label is negative, return\n");
 407       return;
 408     }
 409
 410     TParticle *partCopy = mcpart;
 411     Int_t pdgcode = mcpart->GetPdgCode();
 412     Int_t pdgcodeCopy = pdgcode;
 413
 414     // if the mother is charmed hadron
 415     Bool_t isDirectCharm = kFALSE;
 416     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
 417
 418           // iterate until you find B hadron as a mother or become top ancester
 419           for (Int_t i=1; i<fgkMaxIter; i++){
 420
 421              Int_t jLabel = mcpart->GetFirstMother();
 422              if (jLabel == -1){
 423                isDirectCharm = kTRUE;
 424                break; // if there is no ancester
 425              }
 426              if (jLabel < 0){ // safety protection
 427                AliDebug(1, "Stack label is negative, return\n");
 428                return;
 429              }
 430              // if there is an ancester
 431              TParticle* mother = fStack->Particle(jLabel);
 432              Int_t motherPDG = mother->GetPdgCode();
 433
 434              for (Int_t j=0; j<fNparents; j++){
 435                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
 436              }
 437
 438              mcpart = mother;
 439           } // end of iteration
 440     } // end of if
 441     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 442          for (Int_t i=0; i<fNparents; i++){
 443             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
 444
 445               // fill hadron kinematics
 446               fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
 447               fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
 448               fHist[iq][kHadron][0].fY->Fill(GetRapidity(partCopy));
 449               fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
 450
 451             }
 452          }
 453     } // end of if
 454 }
 455
 456 //__________________________________________
 457 void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed, Int_t icut, Bool_t isbarrel)
 458 {
 459     // decay electron kinematics
 460
 461     if (kquark != kCharm && kquark != kBeauty) {
 462       AliDebug(1, "This task is only for heavy quark QA, return\n");
 463       return;
 464     }
 465     Int_t iq = kquark - kCharm;
 466
 467     if (iTrack < 0) {
 468       AliDebug(1, "Stack label is negative, return\n");
 469       return;
 470     }
 471
 472     TParticle* mcpart = fStack->Particle(iTrack);
 473
 474     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 475     if ( isbarrel && TMath::Abs(mcpart->Eta()) > 0.9 ) return;
 476
 477     Int_t iLabel = mcpart->GetFirstMother();
 478     if (iLabel<0){
 479       AliDebug(1, "Stack label is negative, return\n");
 480       return;
 481     }
 482
 483     TParticle *partMother = fStack->Particle(iLabel);
 484     TParticle *partMotherCopy = partMother;
 485     Int_t maPdgcode = partMother->GetPdgCode();
 486     Int_t maPdgcodeCopy = maPdgcode;
 487
 488     // get electron production vertex
 489     TLorentzVector ePoint;
 490     mcpart->ProductionVertex(ePoint);
 491
 492     // if the mother is charmed hadron
 493     Bool_t isMotherDirectCharm = kFALSE;
 494     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 495
 496           // iterate until you find B hadron as a mother or become top ancester
 497           for (Int_t i=1; i<fgkMaxIter; i++){
 498
 499              Int_t jLabel = partMother->GetFirstMother();
 500              if (jLabel == -1){
 501                isMotherDirectCharm = kTRUE;
 502                break; // if there is no ancester
 503              }
 504              if (jLabel < 0){ // safety protection
 505                AliDebug(1, "Stack label is negative, return\n");
 506                return;
 507              }
 508
 509              // if there is an ancester
 510              TParticle* grandMa = fStack->Particle(jLabel);
 511              Int_t grandMaPDG = grandMa->GetPdgCode();
 512
 513              for (Int_t j=0; j<fNparents; j++){
 514                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 515
 516                   if (kquark == kCharm) return;
 517                   // fill electron kinematics
 518                   fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 519                   fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
 520                   fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart));
 521                   fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
 522
 523                   // fill mother hadron kinematics
 524                   fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
 525                   fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
 526                   fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa));
 527                   fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
 528
 529                   // ratio between pT of electron and pT of mother B hadron
 530                   if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 531
 532                   // distance between electron production point and mother hadron production point
 533                   TLorentzVector debPoint;
 534                   grandMa->ProductionVertex(debPoint);
 535                   TLorentzVector dedistance = ePoint - debPoint;
 536                   fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),dedistance.P());
 537                   return;
 538                 }
 539              }
 540
 541              partMother = grandMa;
 542           } // end of iteration
 543     } // end of if
 544     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 545          for (Int_t i=0; i<fNparents; i++){
 546             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 547
 548               // fill electron kinematics
 549               fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 550               fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 551               fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart));
 552               fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 553
 554               // fill mother hadron kinematics
 555               fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
 556               fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
 557               fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy));
 558               fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
 559
 560               // ratio between pT of electron and pT of mother B hadron
 561               if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
 562
 563               // distance between electron production point and mother hadron production point
 564               TLorentzVector ebPoint;
 565               partMotherCopy->ProductionVertex(ebPoint);
 566               TLorentzVector edistance = ePoint - ebPoint;
 567               fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),edistance.P());
 568             }
 569          }
 570     } // end of if
 571 }
 572
 573
 574 //__________________________________________
 575 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
 576 {
 577        // find mother pdg code and label
 578
 579        if (motherlabel < 0) {
 580          AliDebug(1, "Stack label is negative, return\n");
 581          return;
 582        }
 583        TParticle *heavysMother = fStack->Particle(motherlabel);
 584        motherpdg = heavysMother->GetPdgCode();
 585        grandmotherlabel = heavysMother->GetFirstMother();
 586        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
 587 }
 588
 589 //__________________________________________
 590 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 591 {
 592        // mothertype -1 means this heavy quark coming from hard vertex
 593
 594        TParticle *afterinitialrad1  = fStack->Particle(4);
 595        TParticle *afterinitialrad2  = fStack->Particle(5);
 596
 597        motherlabel = -1;
 598
 599        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
 600          AliDebug(1,"heavy from gluon gluon pair creation!\n");
 601          mothertype = -1;
 602          motherID = fgkGluon;
 603        }
 604        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
 605          AliDebug(1,"heavy from flavor exitation!\n");
 606          mothertype = -1;
 607          motherID = kquark;
 608        }
 609        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
 610          AliDebug(1,"heavy from q-qbar pair creation!\n");
 611          mothertype = -1;
 612          motherID = 1;
 613        }
 614        else {
 615          AliDebug(1,"something strange!\n");
 616          mothertype = -999;
 617          motherlabel = -999;
 618          motherID = -999;
 619        }
 620 }
 621
 622 //__________________________________________
 623 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 624 {
 625        // mothertype -2 means this heavy quark coming from initial state
 626
 627        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
 628          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
 629          motherID = heavysMother->GetPdgCode();
 630          mothertype = -2; // there is mother before initial state radiation
 631          motherlabel = inputmotherlabel;
 632          AliDebug(1,"initial parton shower! \n");
 633
 634          return kTRUE;
 635        }
 636
 637        return kFALSE;
 638 }
 639
 640 //__________________________________________
 641 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 642 {
 643        // mothertype 2 means this heavy quark coming from final state
 644
 645        if (inputmotherlabel > 5){ // mother exist after hard scattering
 646          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
 647          motherID = heavysMother->GetPdgCode();
 648          mothertype = 2; //
 649          motherlabel = inputmotherlabel;
 650          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
 651
 652          return kTRUE;
 653        }
 654        return kFALSE;
 655 }
 656
 657 //__________________________________________
 658 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 659 {
 660       // mark strange behavior
 661
 662        mothertype = -888;
 663        motherlabel = -888;
 664        motherID = -888;
 665        AliDebug(1,"something strange!\n");
 666 }
 667
 668 //__________________________________________
 669 Float_t AliHFEmcQA::GetRapidity(TParticle *part)
 670 {
 671       // return rapidity
 672
 673        Float_t rapidity;
 674        if(part->Energy() - part->Pz() == 0 || part->Energy() + part->Pz() == 0) rapidity=-999;
 675        else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
 676        return rapidity;
 677 }