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 /* $Id$ */
  17
  18 //
  19 // QA class of Heavy Flavor quark and fragmeted/decayed particles
  20 // -Check kinematics of Heavy Quarks/hadrons, and decayed leptons
  21 //    pT, rapidity
  22 //    decay lepton kinematics w/wo acceptance
  23 //    heavy hadron decay length, electron pT fraction carried from decay
  24 // -Check yield of Heavy Quarks/hadrons
  25 //    Number of produced heavy quark
  26 //    Number of produced hadron of given pdg code
  27 //
  28 //
  29 // Authors:
  30 //   MinJung Kweon <minjung@physi.uni-heidelberg.de>
  31 //
  32
  33 #include <TH2F.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TList.h>
  37 #include <TParticle.h>
  38
  39 #include <AliLog.h>
  40 #include <AliMCEvent.h>
  41 #include <AliGenEventHeader.h>
  42 #include <AliAODMCParticle.h>
  43
  44 #include "AliHFEmcQA.h"
  45 #include "AliHFEtools.h"
  46
  47 const Int_t AliHFEmcQA::fgkGluon=21;
  48 const Int_t AliHFEmcQA::fgkMaxGener=10;
  49 const Int_t AliHFEmcQA::fgkMaxIter=100;
  50 const Int_t AliHFEmcQA::fgkqType = 7;    // number of species waiting for QA done
  51
  52 ClassImp(AliHFEmcQA)
  53
  54 //_______________________________________________________________________________________________
  55 AliHFEmcQA::AliHFEmcQA() :
  56                fMCEvent(NULL)
  57         ,fMCHeader(NULL)
  58         ,fMCArray(NULL)
  59         ,fQAhistos(NULL)
  60         ,fNparents(0)
  61 {
  62         // Default constructor
  63
  64   for(Int_t mom = 0; mom < 50; mom++){
  65     fHeavyQuark[mom] = NULL;
  66   }
  67   for(Int_t mom = 0; mom < 2; mom++){
  68     fIsHeavy[mom] = 0;
  69   }
  70
  71 }
  72
  73 //_______________________________________________________________________________________________
  74 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
  75         TObject(p)
  76         ,fMCEvent(NULL)
  77         ,fMCHeader(NULL)
  78         ,fMCArray(NULL)
  79         ,fQAhistos(p.fQAhistos)
  80         ,fNparents(p.fNparents)
  81 {
  82         // Copy constructor
  83   for(Int_t mom = 0; mom < 50; mom++){
  84     fHeavyQuark[mom] = NULL;
  85   }
  86   for(Int_t mom = 0; mom < 2; mom++){
  87     fIsHeavy[mom] = 0;
  88   }
  89 }
  90
  91 //_______________________________________________________________________________________________
  92 AliHFEmcQA&
  93 AliHFEmcQA::operator=(const AliHFEmcQA &)
  94 {
  95   // Assignment operator
  96
  97   AliInfo("Not yet implemented.");
  98   return *this;
  99 }
 100
 101 //_______________________________________________________________________________________________
 102 AliHFEmcQA::~AliHFEmcQA()
 103 {
 104         // Destructor
 105
 106         AliInfo("Analysis Done.");
 107 }
 108
 109 //_______________________________________________________________________________________________
 110 void AliHFEmcQA::PostAnalyze() const
 111 {
 112         //
 113         // Post analysis
 114         //
 115 }
 116
 117 //__________________________________________
 118 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
 119 {
 120   //
 121   // make default histograms
 122   //
 123
 124   if(!qaList) return;
 125
 126   fQAhistos = qaList;
 127   fQAhistos->SetName("MCqa");
 128
 129   CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
 130   CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
 131   CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
 132   CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm
 133   CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
 134   CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
 135   CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm
 136   CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
 137   CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
 138   CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_");        // create histograms for charm
 139   CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_");       // create histograms for beauty
 140   CreateHistograms(AliHFEmcQA::kOthers,3,"mcqa_reccut_");       // create histograms for beauty
 141   CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_");     // create histograms for charm
 142   CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_");    // create histograms for beauty
 143   CreateHistograms(AliHFEmcQA::kOthers,4,"mcqa_recpidcut_");    // create histograms for beauty
 144   CreateHistograms(AliHFEmcQA::kCharm,5,"mcqa_secvtxcut_");     // create histograms for charm
 145   CreateHistograms(AliHFEmcQA::kBeauty,5,"mcqa_secvtxcut_");    // create histograms for beauty
 146   CreateHistograms(AliHFEmcQA::kOthers,5,"mcqa_secvtxcut_");    // create histograms for beauty
 147
 148 // check D spectra
 149   TString kDspecies[9];
 150   kDspecies[0]="411";
 151   kDspecies[1]="421";
 152   kDspecies[2]="431";
 153   kDspecies[3]="4122";
 154   kDspecies[4]="4132";
 155   kDspecies[5]="4232";
 156   kDspecies[6]="4332";
 157   kDspecies[7]="413";
 158   kDspecies[8]="423";
 159
 160   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 161   Int_t iBin[1];
 162   iBin[0] = 44; // bins in pt
 163   Double_t* binEdges[1];
 164   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 165
 166   // bin size is chosen to consider ALICE D measurement
 167   Double_t xbins[13]={0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,8.0,12,16,20,50};
 168   Double_t ybins[10]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5};
 169   TString hname;
 170   for (Int_t iDmeson=0; iDmeson<9; iDmeson++){
 171      hname = "Dmeson"+kDspecies[iDmeson];
 172      fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",12,xbins,9,ybins);
 173      hname = "DmesonLogbin"+kDspecies[iDmeson];
 174      fhDLogbin[iDmeson] = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
 175      if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
 176      if(fQAhistos) fQAhistos->Add(fhDLogbin[iDmeson]);
 177   }
 178
 179 }
 180
 181 //__________________________________________
 182 void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
 183 {
 184   // create histograms
 185
 186   if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
 187     AliDebug(1, "This task is only for heavy quark QA, return\n");
 188     return;
 189   }
 190   Int_t iq = kquark - kCharm;
 191
 192   TString kqTypeLabel[fgkqType];
 193   if (kquark == kCharm){
 194     kqTypeLabel[kQuark]="c";
 195     kqTypeLabel[kantiQuark]="cbar";
 196     kqTypeLabel[kHadron]="cHadron";
 197     kqTypeLabel[keHadron]="ceHadron";
 198     kqTypeLabel[kDeHadron]="nullHadron";
 199     kqTypeLabel[kElectron]="ce";
 200     kqTypeLabel[kElectron2nd]="nulle";
 201   } else if (kquark == kBeauty){
 202     kqTypeLabel[kQuark]="b";
 203     kqTypeLabel[kantiQuark]="bbar";
 204     kqTypeLabel[kHadron]="bHadron";
 205     kqTypeLabel[keHadron]="beHadron";
 206     kqTypeLabel[kDeHadron]="bDeHadron";
 207     kqTypeLabel[kElectron]="be";
 208     kqTypeLabel[kElectron2nd]="bce";
 209   } else if (kquark == kOthers){
 210     kqTypeLabel[kGamma-4]="gammae";
 211     kqTypeLabel[kPi0-4]="pi0e";
 212     kqTypeLabel[kElse-4]="elsee";
 213     kqTypeLabel[kMisID-4]="miside";
 214   }
 215
 216   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 217   Int_t iBin[1];
 218   iBin[0] = 44; // bins in pt
 219   Double_t* binEdges[1];
 220   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 221
 222   TString hname;
 223   if(kquark == kOthers){
 224     for (Int_t iqType = 0; iqType < 4; iqType++ ){
 225        hname = hnopt+"Pt_"+kqTypeLabel[iqType];
 226        //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL
 227        fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25);
 228        hname = hnopt+"Y_"+kqTypeLabel[iqType];
 229        fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
 230        hname = hnopt+"Eta_"+kqTypeLabel[iqType];
 231        fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
 232        // Fill List
 233        if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
 234     }
 235     return;
 236   }
 237   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
 238      if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
 239      hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
 240      fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
 241      hname = hnopt+"Pt_"+kqTypeLabel[iqType];
 242      //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
 243      fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL
 244      hname = hnopt+"Y_"+kqTypeLabel[iqType];
 245      fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
 246      hname = hnopt+"Eta_"+kqTypeLabel[iqType];
 247      fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
 248      // Fill List
 249      if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
 250   }
 251
 252   if (icut == 0){
 253     hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
 254     fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5);
 255     hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
 256     fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
 257   }
 258   hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
 259   fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
 260   hname = hnopt+"PtCorr_"+kqTypeLabel[kQuark];
 261   fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";p_{T} (GeV/c);p_{T} (GeV/c)",200,0,20,200,0,20);
 262   hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
 263   fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
 264   hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
 265   fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
 266   hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
 267   fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
 268   if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
 269 }
 270
 271 //__________________________________________
 272 void AliHFEmcQA::Init()
 273 {
 274   // called at begining every event
 275
 276   for (Int_t i=0; i<2; i++){
 277      fIsHeavy[i] = 0;
 278   }
 279
 280   fNparents = 7;
 281
 282   fParentSelect[0][0] =  411;
 283   fParentSelect[0][1] =  421;
 284   fParentSelect[0][2] =  431;
 285   fParentSelect[0][3] = 4122;
 286   fParentSelect[0][4] = 4132;
 287   fParentSelect[0][5] = 4232;
 288   fParentSelect[0][6] = 4332;
 289
 290   fParentSelect[1][0] =  511;
 291   fParentSelect[1][1] =  521;
 292   fParentSelect[1][2] =  531;
 293   fParentSelect[1][3] = 5122;
 294   fParentSelect[1][4] = 5132;
 295   fParentSelect[1][5] = 5232;
 296   fParentSelect[1][6] = 5332;
 297
 298 }
 299
 300 //__________________________________________
 301 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
 302 {
 303   // get heavy quark kinematics
 304
 305     if (kquark != kCharm && kquark != kBeauty) {
 306       AliDebug(1, "This task is only for heavy quark QA, return\n");
 307       return;
 308     }
 309     Int_t iq = kquark - kCharm;
 310
 311     if (iTrack < 0 || !part) {
 312       AliDebug(1, "Stack label is negative or no mcparticle, return\n");
 313       return;
 314     }
 315
 316     AliMCParticle *mctrack = NULL;
 317     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 318
 319     // select heavy hadron or not fragmented heavy quark
 320     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
 321
 322       TParticle *partMother;
 323       Int_t iLabel;
 324
 325       if (partPdgcode == kquark){ // in case of not fragmented heavy quark
 326         partMother = part;
 327         iLabel = iTrack;
 328       } else{ // in case of heavy hadron, start to search for mother heavy parton
 329         iLabel = part->GetFirstMother();
 330         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 331         if (iLabel>-1) { partMother = mctrack->Particle(); }
 332         else {
 333           AliDebug(1, "Stack label is negative, return\n");
 334           return;
 335         }
 336       }
 337
 338       // heavy parton selection as a mother of heavy hadron
 339       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
 340       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
 341       // should I make a condition that partMother should be quark or diquark? -> not necessary
 342       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
 343       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
 344
 345         if ( abs(partMother->GetPdgCode()) != kquark ){
 346           // search fragmented partons in the same string
 347           Bool_t isSameString = kTRUE;
 348           for (Int_t i=1; i<fgkMaxIter; i++){
 349              iLabel = iLabel - 1;
 350              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 351              if (iLabel>-1) { partMother = mctrack->Particle(); }
 352              else {
 353                AliDebug(1, "Stack label is negative, return\n");
 354                return;
 355              }
 356              if ( abs(partMother->GetPdgCode()) == kquark ) break;
 357              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
 358              if (!isSameString) return;
 359           }
 360         }
 361         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
 362         if (abs(partMother->GetPdgCode()) != kquark) return;
 363
 364         if (fIsHeavy[iq] >= 50) return;
 365         fHeavyQuark[fIsHeavy[iq]] = partMother;
 366         fIsHeavy[iq]++;
 367
 368         // fill kinematics for heavy parton
 369         if (partMother->GetPdgCode() > 0) { // quark
 370           fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 371           fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
 372           fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
 373           fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
 374         } else{ // antiquark
 375           fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 376           fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
 377           fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
 378           fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
 379         }
 380
 381       } // end of heavy parton slection loop
 382
 383     } // end of heavy hadron or quark selection
 384
 385 }
 386
 387 //__________________________________________
 388 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
 389 {
 390   // end of event analysis
 391
 392   if (kquark != kCharm && kquark != kBeauty) {
 393     AliDebug(1, "This task is only for heavy quark QA, return\n");
 394     return;
 395   }
 396   Int_t iq = kquark - kCharm;
 397
 398
 399   // # of heavy quark per event
 400   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
 401   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
 402
 403   Int_t motherID[fgkMaxGener];
 404   Int_t motherType[fgkMaxGener];
 405   Int_t motherLabel[fgkMaxGener];
 406   Int_t ancestorPdg[fgkMaxGener];
 407   Int_t ancestorLabel[fgkMaxGener];
 408
 409   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
 410      motherID[i] = 0;
 411      motherType[i] = 0;
 412      motherLabel[i] = 0;
 413      ancestorPdg[i] = 0;
 414      ancestorLabel[i] = 0;
 415   }
 416
 417
 418   // check history of found heavy quarks
 419   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 420
 421      if(!fHeavyQuark[i]) return;
 422
 423      ancestorLabel[0] = i;
 424      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
 425      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
 426
 427      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
 428      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
 429
 430      Int_t ig = 1;
 431      while (ancestorLabel[ig] != -1){
 432           // in case there is mother, get mother's pdg code and grandmother's label
 433           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
 434           // if mother is still heavy, find again mother's ancestor
 435           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
 436             ig++;
 437             continue; // if it is from same heavy
 438           }
 439           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
 440           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 441           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 442           // if it is not the above case, something is strange
 443           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
 444           break;
 445      }
 446      if (ancestorLabel[ig] == -1){ // from hard scattering
 447        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
 448      }
 449
 450   } // end of found heavy quark loop
 451
 452
 453   // check process type
 454   Int_t processID = 0;
 455   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 456      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
 457   }
 458
 459
 460   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
 461   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
 462
 463      Int_t id1 = ipair*2;
 464      Int_t id2 = ipair*2 + 1;
 465
 466      if (motherType[id1] == 2 && motherType[id2] == 2){
 467        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
 468        else processID = -9;
 469      }
 470      else if (motherType[id1] == -1 && motherType[id2] == -1) {
 471        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
 472          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
 473          else processID = kPairCreationFromq; // q-qbar pair creation
 474        }
 475        else processID = -8;
 476      }
 477      else if (motherType[id1] == -1 || motherType[id2] == -1) {
 478        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
 479          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
 480          else processID = kLightQuarkShower;
 481        }
 482        else processID = -7;
 483      }
 484      else if (motherType[id1] == -2 || motherType[id2] == -2) {
 485        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
 486        else processID = -6;
 487
 488      }
 489      else processID = -5;
 490
 491      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
 492      else fHistComm[iq][0].fProcessID->Fill(processID);
 493      AliDebug(1,Form("Process ID = %d\n",processID));
 494   } // end of # heavy quark pair loop
 495
 496 }
 497
 498 //__________________________________________
 499 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
 500 {
 501     // decay electron kinematics
 502
 503     if (kquark != kCharm && kquark != kBeauty) {
 504       AliDebug(1, "This task is only for heavy quark QA, return\n");
 505       return;
 506     }
 507     Int_t iq = kquark - kCharm;
 508
 509     if(!mcpart){
 510       AliDebug(1, "no mc particle, return\n");
 511       return;
 512     }
 513
 514     Int_t iLabel = mcpart->GetFirstMother();
 515     if (iLabel<0){
 516       AliDebug(1, "Stack label is negative, return\n");
 517       return;
 518     }
 519
 520     TParticle *partCopy = mcpart;
 521     Int_t pdgcode = mcpart->GetPdgCode();
 522     Int_t pdgcodeCopy = pdgcode;
 523
 524     AliMCParticle *mctrack = NULL;
 525
 526     // if the mother is charmed hadron
 527     Bool_t isDirectCharm = kFALSE;
 528     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
 529
 530           // iterate until you find B hadron as a mother or become top ancester
 531           for (Int_t i=1; i<fgkMaxIter; i++){
 532
 533              Int_t jLabel = mcpart->GetFirstMother();
 534              if (jLabel == -1){
 535                isDirectCharm = kTRUE;
 536                break; // if there is no ancester
 537              }
 538              if (jLabel < 0){ // safety protection
 539                AliDebug(1, "Stack label is negative, return\n");
 540                return;
 541              }
 542              // if there is an ancester
 543              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
 544              TParticle* mother = mctrack->Particle();
 545              Int_t motherPDG = mother->GetPdgCode();
 546
 547              for (Int_t j=0; j<fNparents; j++){
 548                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
 549              }
 550
 551              mcpart = mother;
 552           } // end of iteration
 553     } // end of if
 554     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 555          for (Int_t i=0; i<fNparents; i++){
 556             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
 557
 558               // fill hadron kinematics
 559               fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
 560               fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
 561               fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
 562               fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
 563
 564               if(iq==0) {
 565                fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 566                if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[i]->Fill(partCopy->Pt());
 567               }
 568             }
 569          }
 570          // I also want to store D* info to compare with D* measurement
 571          if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
 572                fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 573                if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[7]->Fill(partCopy->Pt());
 574          }
 575          if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
 576                fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 577                if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[8]->Fill(partCopy->Pt());
 578          }
 579     } // end of if
 580 }
 581
 582 //__________________________________________
 583 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
 584 {
 585     // decay electron kinematics
 586
 587     if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
 588       AliDebug(1, "This task is only for heavy quark QA, return\n");
 589       return;
 590     }
 591     Int_t iq = kquark - kCharm;
 592     Bool_t isFinalOpenCharm = kFALSE;
 593
 594     if(!mcpart){
 595       AliDebug(1, "no mcparticle, return\n");
 596       return;
 597     }
 598
 599     if(kquark==kOthers){
 600       Int_t esource = -1;
 601       if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
 602       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
 603       if(esource==0|| esource==1 || esource==2 || esource==3){
 604         fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
 605         fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 606         fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
 607         return;
 608       }
 609       else {
 610         AliDebug(1, "e source is out of defined ranges, return\n");
 611         return;
 612       }
 613     }
 614
 615     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 616
 617     Int_t iLabel = mcpart->GetFirstMother();
 618     if (iLabel<0){
 619       AliDebug(1, "Stack label is negative, return\n");
 620       return;
 621     }
 622
 623     AliMCParticle *mctrack = NULL;
 624     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 625     TParticle *partMother = mctrack->Particle();
 626     TParticle *partMotherCopy = partMother;
 627     Int_t maPdgcode = partMother->GetPdgCode();
 628     Int_t maPdgcodeCopy = maPdgcode;
 629
 630     // get mc primary vertex
 631     /*
 632     TArrayF mcPrimVtx;
 633     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
 634
 635     // get electron production vertex
 636     TLorentzVector ePoint;
 637     mcpart->ProductionVertex(ePoint);
 638
 639     // calculated production vertex to primary vertex (in xy plane)
 640     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
 641     */
 642     Float_t decayLxy = 0;
 643
 644     // if the mother is charmed hadron
 645     Bool_t isMotherDirectCharm = kFALSE;
 646     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 647
 648          for (Int_t i=0; i<fNparents; i++){
 649             if (abs(maPdgcode)==fParentSelect[0][i]){
 650               isFinalOpenCharm = kTRUE;
 651             }
 652          }
 653          if (!isFinalOpenCharm) return ;
 654
 655           // iterate until you find B hadron as a mother or become top ancester
 656           for (Int_t i=1; i<fgkMaxIter; i++){
 657
 658              Int_t jLabel = partMother->GetFirstMother();
 659              if (jLabel == -1){
 660                isMotherDirectCharm = kTRUE;
 661                break; // if there is no ancester
 662              }
 663              if (jLabel < 0){ // safety protection
 664                AliDebug(1, "Stack label is negative, return\n");
 665                return;
 666              }
 667
 668              // if there is an ancester
 669              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
 670              TParticle* grandMa = mctrack->Particle();
 671              Int_t grandMaPDG = grandMa->GetPdgCode();
 672
 673              for (Int_t j=0; j<fNparents; j++){
 674                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 675
 676                   if (kquark == kCharm) return;
 677                   // fill electron kinematics
 678                   fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 679                   fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
 680                   fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 681                   fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
 682
 683                   // fill mother hadron kinematics
 684                   fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
 685                   fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
 686                   fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
 687                   fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
 688
 689                   // ratio between pT of electron and pT of mother B hadron
 690                   if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 691
 692                   // distance between electron production point and primary vertex
 693                   fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
 694                   return;
 695                 }
 696              }
 697
 698              partMother = grandMa;
 699           } // end of iteration
 700     } // end of if
 701     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 702          for (Int_t i=0; i<fNparents; i++){
 703             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 704
 705 //mj weighting to consider measured spectra!!!
 706               Double_t mpt=partMotherCopy->Pt();
 707               Double_t wfactor=2*(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225)); // 2* considering in pythia I used particle + antiparticle differently from the measurement
 708               //Double_t wfactor=(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225));
 709               // fill electron kinematics
 710               if(iq==0){
 711                 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode(),wfactor);
 712                 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt(),wfactor);
 713                 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart),wfactor);
 714                 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta(),wfactor);
 715
 716                 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 717                 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 718                 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 719                 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 720               }
 721               else{
 722                 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 723                 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 724                 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 725                 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 726               }
 727 //--------------
 728               // fill mother hadron kinematics
 729               fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
 730               fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
 731               fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
 732               fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
 733
 734               // ratio between pT of electron and pT of mother B or direct D hadron
 735               if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
 736               fHistComm[iq][icut].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
 737
 738               // distance between electron production point and primary vertex
 739               fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
 740             }
 741          }
 742     } // end of if
 743 }
 744
 745 //____________________________________________________________________
 746 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
 747 {
 748   // decay electron kinematics
 749
 750   if (kquark != kCharm && kquark != kBeauty) {
 751     AliDebug(1, "This task is only for heavy quark QA, return\n");
 752     return;
 753   }
 754
 755   Int_t iq = kquark - kCharm;
 756   Bool_t isFinalOpenCharm = kFALSE;
 757
 758   if(!mcpart){
 759     AliDebug(1, "no mcparticle, return\n");
 760     return;
 761   }
 762
 763   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 764
 765   // mother
 766   Int_t iLabel = mcpart->GetMother();
 767   if (iLabel<0){
 768     AliDebug(1, "Stack label is negative, return\n");
 769     return;
 770   }
 771
 772   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
 773   AliAODMCParticle *partMotherCopy = partMother;
 774   Int_t maPdgcode = partMother->GetPdgCode();
 775   Int_t maPdgcodeCopy = maPdgcode;
 776
 777   Bool_t isMotherDirectCharm = kFALSE;
 778   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 779
 780     for (Int_t i=0; i<fNparents; i++){
 781        if (abs(maPdgcode)==fParentSelect[0][i]){
 782          isFinalOpenCharm = kTRUE;
 783        }
 784     }
 785     if (!isFinalOpenCharm) return;
 786
 787     for (Int_t i=1; i<fgkMaxIter; i++){
 788
 789        Int_t jLabel = partMother->GetMother();
 790        if (jLabel == -1){
 791          isMotherDirectCharm = kTRUE;
 792          break; // if there is no ancester
 793        }
 794        if (jLabel < 0){ // safety protection
 795          AliDebug(1, "Stack label is negative, return\n");
 796          return;
 797        }
 798
 799        // if there is an ancester
 800        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
 801        Int_t grandMaPDG = grandMa->GetPdgCode();
 802
 803        for (Int_t j=0; j<fNparents; j++){
 804           if (abs(grandMaPDG)==fParentSelect[1][j]){
 805
 806             if (kquark == kCharm) return;
 807             // fill electron kinematics
 808             fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 809             fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
 810             fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 811             fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
 812
 813             // fill mother hadron kinematics
 814             fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
 815             fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
 816             fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
 817             fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
 818
 819             return;
 820           }
 821        }
 822
 823        partMother = grandMa;
 824     } // end of iteration
 825   } // end of if
 826   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 827     for (Int_t i=0; i<fNparents; i++){
 828        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 829
 830          // fill electron kinematics
 831          fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 832          fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 833          fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 834          fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 835
 836          // fill mother hadron kinematics
 837          fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
 838          fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
 839          fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
 840          fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
 841
 842        }
 843     }
 844   } // end of if
 845
 846 }
 847
 848 //__________________________________________
 849 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
 850 {
 851        // find mother pdg code and label
 852
 853        if (motherlabel < 0) {
 854          AliDebug(1, "Stack label is negative, return\n");
 855          return;
 856        }
 857        AliMCParticle *mctrack = NULL;
 858        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
 859        TParticle *heavysMother = mctrack->Particle();
 860        motherpdg = heavysMother->GetPdgCode();
 861        grandmotherlabel = heavysMother->GetFirstMother();
 862        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
 863 }
 864
 865 //__________________________________________
 866 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 867 {
 868        // mothertype -1 means this heavy quark coming from hard vertex
 869
 870        AliMCParticle *mctrack1 = NULL;
 871        AliMCParticle *mctrack2 = NULL;
 872        if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
 873        if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
 874        TParticle *afterinitialrad1  = mctrack1->Particle();
 875        TParticle *afterinitialrad2  = mctrack2->Particle();
 876
 877        motherlabel = -1;
 878
 879        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
 880          AliDebug(1,"heavy from gluon gluon pair creation!\n");
 881          mothertype = -1;
 882          motherID = fgkGluon;
 883        }
 884        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
 885          AliDebug(1,"heavy from flavor exitation!\n");
 886          mothertype = -1;
 887          motherID = kquark;
 888        }
 889        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
 890          AliDebug(1,"heavy from q-qbar pair creation!\n");
 891          mothertype = -1;
 892          motherID = 1;
 893        }
 894        else {
 895          AliDebug(1,"something strange!\n");
 896          mothertype = -999;
 897          motherlabel = -999;
 898          motherID = -999;
 899        }
 900 }
 901
 902 //__________________________________________
 903 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 904 {
 905        // mothertype -2 means this heavy quark coming from initial state
 906
 907        AliMCParticle *mctrack = NULL;
 908        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
 909          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
 910          TParticle *heavysMother = mctrack->Particle();
 911          motherID = heavysMother->GetPdgCode();
 912          mothertype = -2; // there is mother before initial state radiation
 913          motherlabel = inputmotherlabel;
 914          AliDebug(1,"initial parton shower! \n");
 915
 916          return kTRUE;
 917        }
 918
 919        return kFALSE;
 920 }
 921
 922 //__________________________________________
 923 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 924 {
 925        // mothertype 2 means this heavy quark coming from final state
 926
 927        AliMCParticle *mctrack = NULL;
 928        if (inputmotherlabel > 5){ // mother exist after hard scattering
 929          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
 930          TParticle *heavysMother = mctrack->Particle();
 931          motherID = heavysMother->GetPdgCode();
 932          mothertype = 2; //
 933          motherlabel = inputmotherlabel;
 934          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
 935
 936          return kTRUE;
 937        }
 938        return kFALSE;
 939 }
 940
 941 //__________________________________________
 942 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 943 {
 944       // mark strange behavior
 945
 946        mothertype = -888;
 947        motherlabel = -888;
 948        motherID = -888;
 949        AliDebug(1,"something strange!\n");
 950 }
 951
 952 //__________________________________________
 953 Int_t AliHFEmcQA::GetSource(AliAODMCParticle * const mcpart)
 954 {
 955   // decay particle's origin
 956
 957   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
 958
 959   Int_t origin = -1;
 960   Bool_t isFinalOpenCharm = kFALSE;
 961
 962   if(!mcpart){
 963     AliDebug(1, "Stack label is negative or no mcparticle, return\n");
 964     return -1;
 965   }
 966
 967   // mother
 968   Int_t iLabel = mcpart->GetMother();
 969   if (iLabel<0){
 970     AliDebug(1, "Stack label is negative, return\n");
 971     return -1;
 972   }
 973
 974   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
 975   Int_t maPdgcode = partMother->GetPdgCode();
 976
 977   // if the mother is charmed hadron
 978   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 979
 980     for (Int_t i=0; i<fNparents; i++){
 981        if (abs(maPdgcode)==fParentSelect[0][i]){
 982          isFinalOpenCharm = kTRUE;
 983        }
 984     }
 985     if (!isFinalOpenCharm) return -1;
 986
 987     // iterate until you find B hadron as a mother or become top ancester
 988     for (Int_t i=1; i<fgkMaxIter; i++){
 989
 990        Int_t jLabel = partMother->GetMother();
 991        if (jLabel == -1){
 992          origin = kDirectCharm;
 993          return origin;
 994        }
 995        if (jLabel < 0){ // safety protection
 996          AliDebug(1, "Stack label is negative, return\n");
 997          return -1;
 998        }
 999
1000        // if there is an ancester
1001        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1002        Int_t grandMaPDG = grandMa->GetPdgCode();
1003
1004        for (Int_t j=0; j<fNparents; j++){
1005           if (abs(grandMaPDG)==fParentSelect[1][j]){
1006             origin = kBeautyCharm;
1007             return origin;
1008           }
1009        }
1010
1011        partMother = grandMa;
1012     } // end of iteration
1013   } // end of if
1014   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1015     for (Int_t i=0; i<fNparents; i++){
1016        if (abs(maPdgcode)==fParentSelect[1][i]){
1017          origin = kDirectBeauty;
1018          return origin;
1019        }
1020     }
1021   } // end of if
1022   else if ( abs(maPdgcode) == 22 ) {
1023     origin = kGamma;
1024     return origin;
1025   } // end of if
1026   else if ( abs(maPdgcode) == 111 ) {
1027     origin = kPi0;
1028     return origin;
1029   } // end of if
1030
1031   return origin;
1032 }
1033
1034 //__________________________________________
1035 Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
1036 {
1037   // decay particle's origin
1038
1039   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1040
1041   Int_t origin = -1;
1042   Bool_t isFinalOpenCharm = kFALSE;
1043
1044   if(!mcpart){
1045     AliDebug(1, "no mcparticle, return\n");
1046     return -1;
1047   }
1048
1049   Int_t iLabel = mcpart->GetFirstMother();
1050   if (iLabel<0){
1051     AliDebug(1, "Stack label is negative, return\n");
1052     return -1;
1053   }
1054
1055   AliMCParticle *mctrack = NULL;
1056   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1057   TParticle *partMother = mctrack->Particle();
1058   Int_t maPdgcode = partMother->GetPdgCode();
1059
1060    // if the mother is charmed hadron
1061    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1062
1063      for (Int_t i=0; i<fNparents; i++){
1064         if (abs(maPdgcode)==fParentSelect[0][i]){
1065           isFinalOpenCharm = kTRUE;
1066         }
1067      }
1068      if (!isFinalOpenCharm) return -1;
1069
1070      // iterate until you find B hadron as a mother or become top ancester
1071      for (Int_t i=1; i<fgkMaxIter; i++){
1072
1073         Int_t jLabel = partMother->GetFirstMother();
1074         if (jLabel == -1){
1075           origin = kDirectCharm;
1076           return origin;
1077         }
1078         if (jLabel < 0){ // safety protection
1079           AliDebug(1, "Stack label is negative, return\n");
1080           return -1;
1081         }
1082
1083         // if there is an ancester
1084         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1085         TParticle* grandMa = mctrack->Particle();
1086         Int_t grandMaPDG = grandMa->GetPdgCode();
1087
1088         for (Int_t j=0; j<fNparents; j++){
1089            if (abs(grandMaPDG)==fParentSelect[1][j]){
1090              origin = kBeautyCharm;
1091              return origin;
1092            }
1093         }
1094
1095         partMother = grandMa;
1096      } // end of iteration
1097    } // end of if
1098    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1099      for (Int_t i=0; i<fNparents; i++){
1100         if (abs(maPdgcode)==fParentSelect[1][i]){
1101           origin = kDirectBeauty;
1102           return origin;
1103         }
1104      }
1105    } // end of if
1106    else if ( abs(maPdgcode) == 22 ) {
1107      origin = kGamma;
1108      return origin;
1109    } // end of if
1110    else if ( abs(maPdgcode) == 111 ) {
1111      origin = kPi0;
1112      return origin;
1113    } // end of if
1114    else origin = kElse;
1115
1116    return origin;
1117 }
1118
1119 //__________________________________________
1120 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1121 {
1122   // decay particle's origin
1123
1124   if(!mcpart){
1125     AliDebug(1, "no mcparticle, return\n");
1126     return -1;
1127   }
1128
1129   if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1130
1131   Int_t origin = -1;
1132   Bool_t isFinalOpenCharm = kFALSE;
1133
1134   Int_t iLabel = mcpart->GetFirstMother();
1135   if (iLabel<0){
1136     AliDebug(1, "Stack label is negative, return\n");
1137     return -1;
1138   }
1139
1140   AliMCParticle *mctrack = NULL;
1141   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1142   TParticle *partMother = mctrack->Particle();
1143   Int_t maPdgcode = partMother->GetPdgCode();
1144
1145    // if the mother is charmed hadron
1146    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1147
1148      for (Int_t i=0; i<fNparents; i++){
1149         if (abs(maPdgcode)==fParentSelect[0][i]){
1150           isFinalOpenCharm = kTRUE;
1151         }
1152      }
1153      if (!isFinalOpenCharm) return -1;
1154
1155      // iterate until you find B hadron as a mother or become top ancester
1156      for (Int_t i=1; i<fgkMaxIter; i++){
1157
1158         Int_t jLabel = partMother->GetFirstMother();
1159         if (jLabel == -1){
1160           origin = kDirectCharm;
1161           return origin;
1162         }
1163         if (jLabel < 0){ // safety protection
1164           AliDebug(1, "Stack label is negative, return\n");
1165           return -1;
1166         }
1167
1168         // if there is an ancester
1169         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1170         TParticle* grandMa = mctrack->Particle();
1171         Int_t grandMaPDG = grandMa->GetPdgCode();
1172
1173         for (Int_t j=0; j<fNparents; j++){
1174            if (abs(grandMaPDG)==fParentSelect[1][j]){
1175              origin = kBeautyCharm;
1176              return origin;
1177            }
1178         }
1179
1180         partMother = grandMa;
1181      } // end of iteration
1182    } // end of if
1183    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1184      for (Int_t i=0; i<fNparents; i++){
1185         if (abs(maPdgcode)==fParentSelect[1][i]){
1186           origin = kDirectBeauty;
1187           return origin;
1188         }
1189      }
1190    } // end of if
1191    else if ( abs(maPdgcode) == 22 ) {
1192      origin = kGamma;
1193      return origin;
1194    } // end of if
1195    else if ( abs(maPdgcode) == 111 ) {
1196      origin = kPi0;
1197      return origin;
1198    } // end of if
1199    else if ( abs(maPdgcode) == 221 ) {
1200      origin = kEta;
1201      return origin;
1202    } // end of if
1203    else if ( abs(maPdgcode) == 223 ) {
1204      origin = kOmega;
1205      return origin;
1206    } // end of if
1207    else if ( abs(maPdgcode) == 333 ) {
1208      origin = kPhi;
1209      return origin;
1210    } // end of if
1211    else if ( abs(maPdgcode) == 331 ) {
1212      origin = kEtaPrime;
1213      return origin;
1214    } // end of if
1215    else if ( abs(maPdgcode) == 113 ) {
1216      origin = kRho0;
1217      return origin;
1218    } // end of if
1219    else origin = kElse;
1220
1221    return origin;
1222 }
1223
1224 //__________________________________________
1225 void AliHFEmcQA::AliHists::FillList(TList *l) const {
1226   //
1227   // Fill Histos into a list for output
1228   //
1229   if(fPdgCode) l->Add(fPdgCode);
1230   if(fPt) l->Add(fPt);
1231   if(fY) l->Add(fY);
1232   if(fEta) l->Add(fEta);
1233 }
1234
1235 //__________________________________________
1236 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
1237   //
1238   // Fill Histos into a list for output
1239   //
1240   if(fNq) l->Add(fNq);
1241   if(fProcessID) l->Add(fProcessID);
1242   if(fePtRatio) l->Add(fePtRatio);
1243   if(fPtCorr) l->Add(fPtCorr);
1244   if(fDePtRatio) l->Add(fDePtRatio);
1245   if(feDistance) l->Add(feDistance);
1246   if(fDeDistance) l->Add(fDeDistance);
1247 }
1248