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