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