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