PWGHF/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 #include "TTreeStream.h"
  36
  37 #include <AliLog.h>
  38 #include <AliMCEvent.h>
  39 #include <AliGenEventHeader.h>
  40 #include <AliAODMCParticle.h>
  41 #include <AliStack.h>
  42
  43 #include "AliHFEmcQA.h"
  44 #include "AliHFEtools.h"
  45 #include "AliHFEcollection.h"
  46
  47 ClassImp(AliHFEmcQA)
  48
  49 //_______________________________________________________________________________________________
  50 AliHFEmcQA::AliHFEmcQA() :
  51 fMCEvent(NULL)
  52   ,fMCHeader(NULL)
  53   ,fMCArray(NULL)
  54   ,fQAhistos(NULL)
  55   ,fMCQACollection(NULL)
  56   ,fNparents(0)
  57   ,fCentrality(0)
  58   ,fPerCentrality(-1)
  59   ,fIsPbPb(kFALSE)
  60   ,fIsppMultiBin(kFALSE)
  61   ,fContainerStep(0)
  62   ,fIsDebugStreamerON(kFALSE)
  63   ,fRecPt(-999)
  64   ,fRecEta(-999)
  65   ,fRecPhi(-999)
  66   ,fLyrhit(0)
  67   ,fLyrstat(0)
  68   ,fHfeImpactR(-999)
  69   ,fHfeImpactnsigmaR(-999)
  70   ,fTreeStream(NULL)
  71 {
  72   // Default constructor
  73   for(Int_t mom = 0; mom < 9; mom++){
  74     fhD[mom] = NULL;
  75   }
  76   for(Int_t mom = 0; mom < 50; mom++){
  77     fHeavyQuark[mom] = NULL;
  78   }
  79   for(Int_t mom = 0; mom < 2; mom++){
  80     fIsHeavy[mom] = 0;
  81   }
  82   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
  83   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
  84 }
  85
  86 //_______________________________________________________________________________________________
  87 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
  88   TObject(p)
  89   ,fMCEvent(NULL)
  90   ,fMCHeader(NULL)
  91   ,fMCArray(NULL)
  92   ,fQAhistos(p.fQAhistos)
  93   ,fMCQACollection(p.fMCQACollection)
  94   ,fNparents(p.fNparents)
  95   ,fCentrality(0)
  96   ,fPerCentrality(-1)
  97   ,fIsPbPb(kFALSE)
  98   ,fIsppMultiBin(kFALSE)
  99   ,fContainerStep(0)
 100   ,fIsDebugStreamerON(kFALSE)
 101   ,fRecPt(-999)
 102   ,fRecEta(-999)
 103   ,fRecPhi(-999)
 104   ,fLyrhit(0)
 105   ,fLyrstat(0)
 106   ,fHfeImpactR(0)
 107   ,fHfeImpactnsigmaR(0)
 108   ,fTreeStream(NULL)
 109 {
 110   // Copy constructor
 111   for(Int_t mom = 0; mom < 9; mom++){
 112     fhD[mom] = NULL;
 113   }
 114   for(Int_t mom = 0; mom < 50; mom++){
 115     fHeavyQuark[mom] = NULL;
 116   }
 117   for(Int_t mom = 0; mom < 2; mom++){
 118     fIsHeavy[mom] = 0;
 119   }
 120   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 121   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 122 }
 123 //_______________________________________________________________________________________________
 124 AliHFEmcQA&
 125 AliHFEmcQA::operator=(const AliHFEmcQA &)
 126 {
 127   // Assignment operator
 128
 129   AliInfo("Not yet implemented.");
 130   return *this;
 131 }
 132
 133 //_______________________________________________________________________________________________
 134 AliHFEmcQA::~AliHFEmcQA()
 135 {
 136   // Destructor
 137
 138   if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
 139   AliInfo("Analysis Done.");
 140 }
 141 //_______________________________________________________________________________________________
 142 void AliHFEmcQA::PostAnalyze() const
 143 {
 144   //
 145   // Post analysis
 146   //
 147 }
 148 //_______________________________________________________________________________________________
 149 void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
 150 {
 151   //
 152   // copy background weighting factors into data member
 153   //
 154
 155   memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 156   memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
 157 }
 158 //__________________________________________
 159 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
 160 {
 161   //
 162   // make default histograms
 163   //
 164
 165   if(!qaList) return;
 166
 167   fQAhistos = qaList;
 168   fQAhistos->SetName("MCqa");
 169
 170   CreateHistograms(AliHFEmcQA::kCharm);               // create histograms for charm
 171   CreateHistograms(AliHFEmcQA::kBeauty);              // create histograms for beauty
 172   CreateHistograms(AliHFEmcQA::kOthers);              // create histograms for beauty
 173
 174   // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
 175   const Int_t nbspecies = 9;
 176   TString kDspecies[nbspecies];
 177   kDspecies[0]="411";   //D+
 178   kDspecies[1]="421";   //D0
 179   kDspecies[2]="431";   //Ds+
 180   kDspecies[3]="4122";  //Lambdac+
 181   kDspecies[4]="4132";  //Ksic0
 182   kDspecies[5]="4232";  //Ksic+
 183   kDspecies[6]="4332";  //OmegaC0
 184   kDspecies[7]="413";   //D*(2010)+
 185   kDspecies[8]="423";   //D*(2007)0
 186
 187   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 188   Int_t iBin[2];
 189   iBin[0] = 44; // bins in pt for log binning
 190   iBin[1] = 23; // bins in pt for pi0 measurement binning
 191   Double_t* binEdges[1];
 192   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 193
 194   // bin size is chosen to consider ALICE D measurement
 195   const Int_t nptbins = 15;
 196   const Int_t nybins = 9;
 197   Double_t xbins[nptbins+1]={0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,12,16,24,32,40,50}; //pt binning for the final 7 TeV D measurement
 198   Double_t ybins[nybins+1]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5}; // y binning
 199   TString hname;
 200   for (Int_t iDmeson=0; iDmeson<nbspecies; iDmeson++){
 201      hname = "Dmeson"+kDspecies[iDmeson];
 202      fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",nptbins,xbins,nybins,ybins);
 203      if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
 204   }
 205
 206   const Double_t kPtRange[24] = {0.,0.3,0.4,0.5,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.6,2.8,3.,3.5,4.,5.,6.,7.,20.,30.}; // to cope with Ana's bin
 207
 208   Int_t kNcent;
 209   if(fIsPbPb) kNcent=11;
 210   else
 211   {
 212       if(fIsppMultiBin) kNcent=8;
 213       else kNcent = 1;
 214   }
 215
 216   fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
 217
 218   for(Int_t centbin=0; centbin<kNcent; centbin++)
 219   {
 220       fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
 221       fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
 222       fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
 223       fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
 224       fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
 225       fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
 226
 227       fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 228       fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 229       fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 230       fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 231       fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 232       fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 233       fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 234       fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 235
 236       fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 237       fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 238       fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 239       fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 240       fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 241       fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 242   }
 243
 244   fQAhistos->Add(fMCQACollection->GetList());
 245
 246   if(!fTreeStream && fIsDebugStreamerON){
 247    fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
 248   }
 249
 250 }
 251
 252 //__________________________________________
 253 void AliHFEmcQA::CreateHistograms(const Int_t kquark)
 254 {
 255   // create histograms
 256
 257   if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
 258     AliDebug(1, "This task is only for heavy quark QA, return\n");
 259     return;
 260   }
 261   Int_t iq = kquark - kCharm;
 262
 263   TString kqTypeLabel[fgkqType];
 264   if (kquark == kCharm){
 265     kqTypeLabel[kQuark]="c";
 266     kqTypeLabel[kantiQuark]="cbar";
 267     kqTypeLabel[kHadron]="cHadron";
 268     kqTypeLabel[keHadron]="ceHadron";
 269     kqTypeLabel[kDeHadron]="nullHadron";
 270     kqTypeLabel[kElectron]="ce";
 271     kqTypeLabel[kElectron2nd]="nulle";
 272   } else if (kquark == kBeauty){
 273     kqTypeLabel[kQuark]="b";
 274     kqTypeLabel[kantiQuark]="bbar";
 275     kqTypeLabel[kHadron]="bHadron";
 276     kqTypeLabel[keHadron]="beHadron";
 277     kqTypeLabel[kDeHadron]="bDeHadron";
 278     kqTypeLabel[kElectron]="be";
 279     kqTypeLabel[kElectron2nd]="bce";
 280   } else if (kquark == kOthers){
 281     kqTypeLabel[kGamma-4]="gammae";
 282     kqTypeLabel[kPi0-4]="pi0e";
 283     kqTypeLabel[kElse-4]="elsee";
 284     kqTypeLabel[kMisID-4]="miside";
 285   }
 286
 287   TString kqEtaRangeLabel[fgkEtaRanges];
 288   kqEtaRangeLabel[0] = "mcqa_";
 289   kqEtaRangeLabel[1] = "mcqa_barrel_";
 290   kqEtaRangeLabel[2] = "mcqa_unitY_";
 291
 292   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 293   const Int_t nptbinning1 = 35;
 294   Int_t iBin[2];
 295   iBin[0] = 44; // bins in pt
 296   iBin[1] = nptbinning1; // bins in pt
 297   Double_t* binEdges[1];
 298   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 299
 300   // new binning for final electron analysis
 301   const Double_t kPtbinning1[nptbinning1+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
 302
 303   const Int_t ndptbins = 500;
 304   Double_t xcorrbin[ndptbins+1];
 305   for (int icorrbin = 0; icorrbin< ndptbins+1; icorrbin++){
 306     xcorrbin[icorrbin]=icorrbin*0.1;
 307   }
 308
 309   Int_t fCentrmax = 0;
 310   if(!fIsPbPb&&!fIsppMultiBin) fCentrmax=0;
 311   else fCentrmax = kCentBins;
 312   TString hname;
 313   if(kquark == kOthers){
 314    for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
 315        for (Int_t iqType = 0; iqType < 4; iqType++ ){
 316            for(Int_t icentr = 0; icentr < (fCentrmax+1); icentr++)
 317            {
 318                hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
 319                fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),60,0.25,30.25);
 320                hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
 321                fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
 322                hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
 323                fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
 324                // Fill List
 325                if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
 326            }
 327        }
 328    }
 329    return;
 330   }
 331   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
 332    for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
 333        if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
 334         for(Int_t icentr = 0; icentr<(fCentrmax+1); icentr++)
 335            {
 336                hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
 337                fHist[iq][iqType][icut][icentr].fPdgCode = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;PdgCode",hname.Data(),icentr),20001,-10000.5,10000.5);
 338                hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
 339                fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),iBin[1],kPtbinning1); // new binning
 340                hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
 341                fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
 342                hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
 343                fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
 344                // Fill List
 345                if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
 346            }
 347    }
 348   }
 349
 350   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
 351     hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
 352     fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 353     hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
 354     fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 355     hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
 356     fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 357     hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
 358     fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 359
 360     hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
 361     fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
 362     hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
 363     fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
 364     hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
 365     fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
 366     hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
 367     fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
 368
 369     if(icut <1){
 370       hname = kqEtaRangeLabel[icut]+"PtCorrDinein_"+kqTypeLabel[kQuark];
 371       fHistComm[iq][icut].fPtCorrDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 372       hname = kqEtaRangeLabel[icut]+"PtCorrDineout_"+kqTypeLabel[kQuark];
 373       fHistComm[iq][icut].fPtCorrDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 374       hname = kqEtaRangeLabel[icut]+"PtCorrDoutein_"+kqTypeLabel[kQuark];
 375       fHistComm[iq][icut].fPtCorrDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 376       hname = kqEtaRangeLabel[icut]+"PtCorrDouteout_"+kqTypeLabel[kQuark];
 377       fHistComm[iq][icut].fPtCorrDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 378
 379       hname = kqEtaRangeLabel[icut]+"PtCorrDpDinein_"+kqTypeLabel[kQuark];
 380       fHistComm[iq][icut].fPtCorrDpDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 381       hname = kqEtaRangeLabel[icut]+"PtCorrDpDineout_"+kqTypeLabel[kQuark];
 382       fHistComm[iq][icut].fPtCorrDpDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 383       hname = kqEtaRangeLabel[icut]+"PtCorrDpDoutein_"+kqTypeLabel[kQuark];
 384       fHistComm[iq][icut].fPtCorrDpDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 385       hname = kqEtaRangeLabel[icut]+"PtCorrDpDouteout_"+kqTypeLabel[kQuark];
 386       fHistComm[iq][icut].fPtCorrDpDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 387
 388       hname = kqEtaRangeLabel[icut]+"PtCorrD0Dinein_"+kqTypeLabel[kQuark];
 389       fHistComm[iq][icut].fPtCorrD0Dinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 390       hname = kqEtaRangeLabel[icut]+"PtCorrD0Dineout_"+kqTypeLabel[kQuark];
 391       fHistComm[iq][icut].fPtCorrD0Dineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 392       hname = kqEtaRangeLabel[icut]+"PtCorrD0Doutein_"+kqTypeLabel[kQuark];
 393       fHistComm[iq][icut].fPtCorrD0Doutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 394       hname = kqEtaRangeLabel[icut]+"PtCorrD0Douteout_"+kqTypeLabel[kQuark];
 395       fHistComm[iq][icut].fPtCorrD0Douteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 396
 397       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDinein_"+kqTypeLabel[kQuark];
 398       fHistComm[iq][icut].fPtCorrDrestDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 399       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDineout_"+kqTypeLabel[kQuark];
 400       fHistComm[iq][icut].fPtCorrDrestDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 401       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDoutein_"+kqTypeLabel[kQuark];
 402       fHistComm[iq][icut].fPtCorrDrestDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 403       hname = kqEtaRangeLabel[icut]+"PtCorrDrestDouteout_"+kqTypeLabel[kQuark];
 404       fHistComm[iq][icut].fPtCorrDrestDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 405
 406       hname = kqEtaRangeLabel[icut]+"fEtaCorrD_"+kqTypeLabel[kQuark];
 407       fHistComm[iq][icut].fEtaCorrD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 408       hname = kqEtaRangeLabel[icut]+"fEtaCorrDp_"+kqTypeLabel[kQuark];
 409       fHistComm[iq][icut].fEtaCorrDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 410       hname = kqEtaRangeLabel[icut]+"fEtaCorrD0_"+kqTypeLabel[kQuark];
 411       fHistComm[iq][icut].fEtaCorrD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 412       hname = kqEtaRangeLabel[icut]+"fEtaCorrDrest_"+kqTypeLabel[kQuark];
 413       fHistComm[iq][icut].fEtaCorrDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 414
 415       hname = kqEtaRangeLabel[icut]+"fEtaCorrGD_"+kqTypeLabel[kQuark];
 416       fHistComm[iq][icut].fEtaCorrGD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 417       hname = kqEtaRangeLabel[icut]+"fEtaCorrGDp_"+kqTypeLabel[kQuark];
 418       fHistComm[iq][icut].fEtaCorrGDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 419       hname = kqEtaRangeLabel[icut]+"fEtaCorrGD0_"+kqTypeLabel[kQuark];
 420       fHistComm[iq][icut].fEtaCorrGD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 421       hname = kqEtaRangeLabel[icut]+"fEtaCorrGDrest_"+kqTypeLabel[kQuark];
 422       fHistComm[iq][icut].fEtaCorrGDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10);
 423
 424       hname = kqEtaRangeLabel[icut]+"fEtaCorrB_"+kqTypeLabel[kQuark];
 425       fHistComm[iq][icut].fEtaCorrB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10);
 426       hname = kqEtaRangeLabel[icut]+"fEtaCorrGB_"+kqTypeLabel[kQuark];
 427       fHistComm[iq][icut].fEtaCorrGB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10);
 428
 429       hname = kqEtaRangeLabel[icut]+"PtCorrBinein_"+kqTypeLabel[kQuark];
 430       fHistComm[iq][icut].fPtCorrBinein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 431       hname = kqEtaRangeLabel[icut]+"PtCorrBineout_"+kqTypeLabel[kQuark];
 432       fHistComm[iq][icut].fPtCorrBineout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 433       hname = kqEtaRangeLabel[icut]+"PtCorrBoutein_"+kqTypeLabel[kQuark];
 434       fHistComm[iq][icut].fPtCorrBoutein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 435       hname = kqEtaRangeLabel[icut]+"PtCorrBouteout_"+kqTypeLabel[kQuark];
 436       fHistComm[iq][icut].fPtCorrBouteout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
 437     }
 438     if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
 439   }
 440
 441
 442   hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
 443   fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
 444   hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
 445   fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
 446
 447 }
 448
 449 //__________________________________________
 450 void AliHFEmcQA::Init()
 451 {
 452   // called at begining every event
 453
 454   for (Int_t i=0; i<2; i++){
 455      fIsHeavy[i] = 0;
 456   }
 457
 458   fNparents = 7;
 459
 460   fParentSelect[0][0] =  411; //D+
 461   fParentSelect[0][1] =  421; //D0
 462   fParentSelect[0][2] =  431; //Ds+
 463   fParentSelect[0][3] = 4122; //Lambdac+
 464   fParentSelect[0][4] = 4132; //Ksic0
 465   fParentSelect[0][5] = 4232; //Ksic+
 466   fParentSelect[0][6] = 4332; //OmegaC0
 467
 468   fParentSelect[1][0] =  511; //B0
 469   fParentSelect[1][1] =  521; //B+
 470   fParentSelect[1][2] =  531; //Bs0
 471   fParentSelect[1][3] = 5122; //Lambdab0
 472   fParentSelect[1][4] = 5132; //Ksib-
 473   fParentSelect[1][5] = 5232; //Ksib0
 474   fParentSelect[1][6] = 5332; //Omegab-
 475
 476
 477 }
 478
 479 //__________________________________________
 480 void AliHFEmcQA::GetMesonKine()
 481 {
 482   //
 483   // get meson pt spectra
 484   //
 485
 486   AliVParticle *mctrack2 = NULL;
 487   AliMCParticle *mctrack0 = NULL;
 488   AliVParticle *mctrackdaugt= NULL;
 489   AliMCParticle *mctrackd= NULL;
 490   Int_t id1=0, id2=0;
 491
 492
 493   if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
 494   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 495
 496   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
 497      if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
 498      TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
 499      if(!mcpart0) continue;
 500      mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
 501      if(!mctrack0) continue;
 502
 503 //     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 504
 505      if(TMath::Abs(mctrack0->PdgCode()) == 111) // pi0
 506        {
 507           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 508             fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
 509             fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 510           }
 511           id1=mctrack0->GetFirstDaughter();
 512           id2=mctrack0->GetLastDaughter();
 513           if(!((id2-id1)==2)) continue;
 514           for(int idx=id1; idx<=id2; idx++){
 515             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 516             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 517             if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 518              fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
 519           }
 520        }
 521      else if(TMath::Abs(mctrack0->PdgCode()) == 221) // eta
 522        {
 523           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 524             fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
 525             fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 526           }
 527           id1=mctrack0->GetFirstDaughter();
 528           id2=mctrack0->GetLastDaughter();
 529           if(!((id2-id1)==2||(id2-id1)==3)) continue;
 530           for(int idx=id1; idx<=id2; idx++){
 531             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 532             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 533             if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 534              fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
 535           }
 536        }
 537      else if(TMath::Abs(mctrack0->PdgCode()) == 223) // omega
 538        {
 539           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 540             fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
 541             fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 542           }
 543           id1=mctrack0->GetFirstDaughter();
 544           id2=mctrack0->GetLastDaughter();
 545           if(!((id2-id1)==1||(id2-id1)==2)) continue;
 546           for(int idx=id1; idx<=id2; idx++){
 547             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 548             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 549             if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 550              fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
 551           }
 552        }
 553      else if(TMath::Abs(mctrack0->PdgCode()) == 333) // phi
 554        {
 555           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 556             fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
 557             fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 558           }
 559           id1=mctrack0->GetFirstDaughter();
 560           id2=mctrack0->GetLastDaughter();
 561           if(!((id2-id1)==1)) continue;
 562           for(int idx=id1; idx<=id2; idx++){
 563             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 564             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 565             if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 566              fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
 567           }
 568        }
 569      else if(TMath::Abs(mctrack0->PdgCode()) == 331) // eta prime
 570        {
 571           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 572             fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
 573             fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 574           }
 575           id1=mctrack0->GetFirstDaughter();
 576           id2=mctrack0->GetLastDaughter();
 577           if(!((id2-id1)==2||(id2-id1)==3)) continue;
 578           for(int idx=id1; idx<=id2; idx++){
 579             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 580             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 581             if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 582              fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
 583           }
 584        }
 585      else if(TMath::Abs(mctrack0->PdgCode()) == 113) // rho
 586        {
 587           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 588             fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
 589             fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 590           }
 591           id1=mctrack0->GetFirstDaughter();
 592           id2=mctrack0->GetLastDaughter();
 593           if(!((id2-id1)==1)) continue;
 594           for(int idx=id1; idx<=id2; idx++){
 595             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 596             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 597             if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 598              fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
 599           }
 600        }
 601      else if(TMath::Abs(mctrack0->PdgCode()) == 321) // kaon+-
 602        {
 603           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 604             fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 605           }
 606        }
 607      else if(TMath::Abs(mctrack0->PdgCode()) == 130) // k0L
 608        {
 609           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 610             fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 611           }
 612        }
 613      }
 614
 615 }
 616 //__________________________________________
 617 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
 618 {
 619   // get heavy quark kinematics
 620
 621     if (kquark != kCharm && kquark != kBeauty) {
 622       AliDebug(1, "This task is only for heavy quark QA, return\n");
 623       return;
 624     }
 625     Int_t iq = kquark - kCharm;
 626
 627     if (iTrack < 0 || !part) {
 628       AliDebug(1, "Stack label is negative or no mcparticle, return\n");
 629       return;
 630     }
 631
 632     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 633
 634     AliMCParticle *mctrack = NULL;
 635     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 636
 637     // select heavy hadron or not fragmented heavy quark
 638     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
 639
 640       TParticle *partMother;
 641       Int_t iLabel;
 642
 643       if (partPdgcode == kquark){ // in case of not fragmented heavy quark
 644         partMother = part;
 645         iLabel = iTrack;
 646       } else{ // in case of heavy hadron, start to search for mother heavy parton
 647         iLabel = part->GetFirstMother();
 648         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 649         if (iLabel>-1) { partMother = mctrack->Particle(); }
 650         else {
 651           AliDebug(1, "Stack label is negative, return\n");
 652           return;
 653         }
 654       }
 655
 656       // heavy parton selection as a mother of heavy hadron
 657       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
 658       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
 659       // should I make a condition that partMother should be quark or diquark? -> not necessary
 660       if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
 661       //if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
 662
 663         if ( TMath::Abs(partMother->GetPdgCode()) != kquark ){
 664           // search fragmented partons in the same string
 665           Bool_t isSameString = kTRUE;
 666           for (Int_t i=1; i<fgkMaxIter; i++){
 667              iLabel = iLabel - 1;
 668              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 669              if (iLabel>-1) { partMother = mctrack->Particle(); }
 670              else {
 671                AliDebug(1, "Stack label is negative, return\n");
 672                return;
 673              }
 674              if ( TMath::Abs(partMother->GetPdgCode()) == kquark ) break;
 675              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
 676              if (!isSameString) return;
 677           }
 678         }
 679         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
 680         if (TMath::Abs(partMother->GetPdgCode()) != kquark) return;
 681
 682         if (fIsHeavy[iq] >= 50) return;
 683         fHeavyQuark[fIsHeavy[iq]] = partMother;
 684         fIsHeavy[iq]++;
 685
 686         // fill kinematics for heavy parton
 687         if (partMother->GetPdgCode() > 0) { // quark
 688           fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
 689           fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
 690           fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
 691           fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
 692         } else{ // antiquark
 693           fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
 694           fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
 695           fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
 696           fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
 697         }
 698
 699       } // end of heavy parton slection loop
 700
 701     } // end of heavy hadron or quark selection
 702
 703 }
 704
 705 //__________________________________________
 706 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
 707 {
 708   // end of event analysis
 709
 710   if (kquark != kCharm && kquark != kBeauty) {
 711     AliDebug(1, "This task is only for heavy quark QA, return\n");
 712     return;
 713   }
 714   Int_t iq = kquark - kCharm;
 715
 716
 717   // # of heavy quark per event
 718   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
 719   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
 720
 721   Int_t motherID[fgkMaxGener];
 722   Int_t motherType[fgkMaxGener];
 723   Int_t motherLabel[fgkMaxGener];
 724   Int_t ancestorPdg[fgkMaxGener];
 725   Int_t ancestorLabel[fgkMaxGener];
 726
 727   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
 728      motherID[i] = 0;
 729      motherType[i] = 0;
 730      motherLabel[i] = 0;
 731      ancestorPdg[i] = 0;
 732      ancestorLabel[i] = 0;
 733   }
 734
 735
 736   // check history of found heavy quarks
 737   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 738
 739      if(!fHeavyQuark[i]) return;
 740
 741      ancestorLabel[0] = i;
 742      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
 743      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
 744
 745      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
 746      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
 747
 748      Int_t ig = 1;
 749      while (ancestorLabel[ig] != -1){
 750           // in case there is mother, get mother's pdg code and grandmother's label
 751           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
 752           // if mother is still heavy, find again mother's ancestor
 753           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
 754             ig++;
 755             continue; // if it is from same heavy
 756           }
 757           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
 758           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 759           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 760           // if it is not the above case, something is strange
 761           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
 762           break;
 763      }
 764      if (ancestorLabel[ig] == -1){ // from hard scattering
 765        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
 766      }
 767
 768   } // end of found heavy quark loop
 769
 770
 771   // check process type
 772   Int_t processID = 0;
 773   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 774      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
 775   }
 776
 777
 778   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
 779   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
 780
 781      Int_t id1 = ipair*2;
 782      Int_t id2 = ipair*2 + 1;
 783
 784      if (motherType[id1] == 2 && motherType[id2] == 2){
 785        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
 786        else processID = -9;
 787      }
 788      else if (motherType[id1] == -1 && motherType[id2] == -1) {
 789        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
 790          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
 791          else processID = kPairCreationFromq; // q-qbar pair creation
 792        }
 793        else processID = -8;
 794      }
 795      else if (motherType[id1] == -1 || motherType[id2] == -1) {
 796        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
 797          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
 798          else processID = kLightQuarkShower;
 799        }
 800        else processID = -7;
 801      }
 802      else if (motherType[id1] == -2 || motherType[id2] == -2) {
 803        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
 804        else processID = -6;
 805
 806      }
 807      else processID = -5;
 808
 809      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
 810      else fHistComm[iq][0].fProcessID->Fill(processID);
 811      AliDebug(1,Form("Process ID = %d\n",processID));
 812   } // end of # heavy quark pair loop
 813
 814 }
 815
 816 //__________________________________________
 817 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
 818 {
 819     // decay electron kinematics
 820
 821     if (kquark != kCharm && kquark != kBeauty) {
 822       AliDebug(1, "This task is only for heavy quark QA, return\n");
 823       return;
 824     }
 825     Int_t iq = kquark - kCharm;
 826
 827     if(!mcpart){
 828       AliDebug(1, "no mc particle, return\n");
 829       return;
 830     }
 831
 832     Int_t iLabel = mcpart->GetFirstMother();
 833     if (iLabel<0){
 834       AliDebug(1, "Stack label is negative, return\n");
 835       return;
 836     }
 837
 838     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 839
 840     TParticle *partCopy = mcpart;
 841     Int_t pdgcode = mcpart->GetPdgCode();
 842     Int_t pdgcodeCopy = pdgcode;
 843
 844     AliMCParticle *mctrack = NULL;
 845
 846     // if the mother is charmed hadron
 847     Bool_t isDirectCharm = kFALSE;
 848     if ( int(TMath::Abs(pdgcode)/100.) == kCharm || int(TMath::Abs(pdgcode)/1000.) == kCharm ) {
 849
 850           // iterate until you find B hadron as a mother or become top ancester
 851           for (Int_t i=1; i<fgkMaxIter; i++){
 852
 853              Int_t jLabel = mcpart->GetFirstMother();
 854              if (jLabel == -1){
 855                isDirectCharm = kTRUE;
 856                break; // if there is no ancester
 857              }
 858              if (jLabel < 0){ // safety protection
 859                AliDebug(1, "Stack label is negative, return\n");
 860                return;
 861              }
 862              // if there is an ancester
 863              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
 864              TParticle* mother = mctrack->Particle();
 865              Int_t motherPDG = mother->GetPdgCode();
 866
 867              for (Int_t j=0; j<fNparents; j++){
 868                 if (TMath::Abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
 869              }
 870
 871              mcpart = mother;
 872           } // end of iteration
 873     } // end of if
 874     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 875          for (Int_t i=0; i<fNparents; i++){
 876             if (TMath::Abs(pdgcodeCopy)==fParentSelect[iq][i]){
 877
 878               // fill hadron kinematics
 879               fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
 880               fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
 881               fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
 882               fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
 883
 884               if(iq==0) {
 885                fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 886               }
 887             }
 888          }
 889          // I also want to store D* info to compare with D* measurement
 890          if (TMath::Abs(pdgcodeCopy)==413 && iq==0) { //D*+
 891                fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 892          }
 893          if (TMath::Abs(pdgcodeCopy)==423 && iq==0) { //D*0
 894                fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 895          }
 896     } // end of if
 897 }
 898
 899 //__________________________________________
 900 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed)
 901 {
 902     // decay electron kinematics
 903
 904     if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
 905       AliDebug(1, "This task is only for heavy quark QA, return\n");
 906       return;
 907     }
 908     Int_t iq = kquark - kCharm;
 909     Bool_t isFinalOpenCharm = kFALSE;
 910
 911     if(!mcpart){
 912       AliDebug(1, "no mcparticle, return\n");
 913       return;
 914     }
 915
 916     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 917
 918     Double_t eabsEta = TMath::Abs(mcpart->Eta());
 919     Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
 920
 921     if(kquark==kOthers){
 922       Int_t esource = -1;
 923       if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
 924       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
 925       if(esource==0|| esource==1 || esource==2 || esource==3){
 926         fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
 927         fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 928         fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
 929         if(eabsEta<0.9){
 930           fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
 931           fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 932           fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
 933         }
 934         if(eabsY<0.5){
 935           fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
 936           fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 937           fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
 938         }
 939         return;
 940       }
 941       else {
 942         AliDebug(1, "e source is out of defined ranges, return\n");
 943         return;
 944       }
 945     }
 946
 947     if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
 948
 949     Int_t iLabel = mcpart->GetFirstMother();
 950     if (iLabel<0){
 951       AliDebug(1, "Stack label is negative, return\n");
 952       return;
 953     }
 954
 955     AliMCParticle *mctrack = NULL;
 956     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 957     TParticle *partMother = mctrack->Particle();
 958     TParticle *partMotherCopy = partMother;
 959     Int_t maPdgcode = partMother->GetPdgCode();
 960     Int_t maPdgcodeCopy = maPdgcode;
 961
 962     // get mc primary vertex
 963     /*
 964     TArrayF mcPrimVtx;
 965     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
 966
 967     // get electron production vertex
 968     TLorentzVector ePoint;
 969     mcpart->ProductionVertex(ePoint);
 970
 971     // calculated production vertex to primary vertex (in xy plane)
 972     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
 973     */
 974     Float_t decayLxy = 0;
 975
 976     // if the mother is charmed hadron
 977     Bool_t isMotherDirectCharm = kFALSE;
 978     if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
 979
 980          for (Int_t i=0; i<fNparents; i++){
 981             if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
 982               isFinalOpenCharm = kTRUE;
 983             }
 984          }
 985          if (!isFinalOpenCharm) return ;
 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->GetFirstMother();
 991              if (jLabel == -1){
 992                isMotherDirectCharm = kTRUE;
 993                break; // if there is no ancester
 994              }
 995              if (jLabel < 0){ // safety protection
 996                AliDebug(1, "Stack label is negative, return\n");
 997                return;
 998              }
 999
1000              // if there is an ancester
1001              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
1002              TParticle* grandMa = mctrack->Particle();
1003              Int_t grandMaPDG = grandMa->GetPdgCode();
1004
1005              for (Int_t j=0; j<fNparents; j++){
1006                 if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1007
1008                   if (kquark == kCharm) return;
1009                   // fill electron kinematics
1010                   fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1011                   fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1012                   fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1013                   fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1014
1015                   // fill mother hadron kinematics
1016                   fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1017                   fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1018                   fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1019                   fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1020
1021                   if(eabsEta<0.9){
1022                     fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1023                     fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1024                     fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1025                     fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1026
1027                     // fill mother hadron kinematics
1028                     fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1029                     fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1030                     fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1031                     fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1032                   }
1033
1034                   if(eabsY<0.5){
1035                     fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1036                     fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1037                     fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1038                     fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1039
1040                     // fill mother hadron kinematics
1041                     fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1042                     fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1043                     fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1044                     fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1045                   }
1046
1047                   //mj: to calculate B to e eta correlation to calculate total heavy quark cross section
1048                   Int_t kLabel0 = grandMa->GetFirstMother();
1049                   Bool_t isGGrandmaYes = kFALSE;
1050                   Double_t ggmrapidwstmp=0;
1051                   if (!(kLabel0 < 0)){ // safety protection
1052                     if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel0))))){
1053                       TParticle* ggrandMatmp = mctrack->Particle();
1054                       Int_t ggrandMaPDGtmp = ggrandMatmp->GetPdgCode();
1055                       if ( int(TMath::Abs(ggrandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(ggrandMaPDGtmp)/1000.) == kBeauty) isGGrandmaYes = kTRUE;
1056                       ggmrapidwstmp = AliHFEtools::GetRapidity(ggrandMatmp);
1057                     }
1058                   }
1059
1060                   Double_t gmrapidwstmp0 = AliHFEtools::GetRapidity(grandMa);
1061                   Double_t eetawstmp0 = mcpart->Eta();
1062
1063                   Double_t gmrapidtmp0 = TMath::Abs(gmrapidwstmp0);
1064                   Double_t eetatmp0 = TMath::Abs(eetawstmp0);
1065
1066                   fHistComm[iq][0].fEtaCorrB->Fill(gmrapidwstmp0,eetawstmp0);
1067                   if(isGGrandmaYes) fHistComm[iq][0].fEtaCorrGB->Fill(ggmrapidwstmp,eetawstmp0);
1068                   else fHistComm[iq][0].fEtaCorrGB->Fill(gmrapidwstmp0,eetawstmp0);
1069
1070                   if(gmrapidtmp0<0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBinein->Fill(grandMa->Pt(),mcpart->Pt());
1071                   else if(gmrapidtmp0<0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBineout->Fill(grandMa->Pt(),mcpart->Pt());
1072                   else if(gmrapidtmp0>0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBoutein->Fill(grandMa->Pt(),mcpart->Pt());
1073                   else if(gmrapidtmp0>0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBouteout->Fill(grandMa->Pt(),mcpart->Pt());
1074                   //======================================================================================
1075
1076                   // ratio between pT of electron and pT of mother B hadron
1077                   if(grandMa->Pt()) {
1078                     fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1079                     if(eabsEta<0.9){
1080                       fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1081                     }
1082                     if(eabsY<0.5){
1083                       fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
1084                     }
1085                   }
1086
1087                   // distance between electron production point and primary vertex
1088                   fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1089                   if(eabsEta<0.9){
1090                     fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1091                   }
1092                   if(eabsY<0.5){
1093                     fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
1094                   }
1095                   return;
1096                 }
1097              }
1098
1099              partMother = grandMa;
1100           } // end of iteration
1101     } // end of if
1102     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1103          for (Int_t i=0; i<fNparents; i++){
1104             if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1105
1106               fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1107               fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1108               fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1109               fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1110
1111               // fill mother hadron kinematics
1112               fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1113               fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1114               fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1115               fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1116
1117               if(eabsEta<0.9){
1118                 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1119                 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1120                 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1121                 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1122
1123                 // fill mother hadron kinematics
1124                 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1125                 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1126                 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1127                 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1128               }
1129
1130               if(eabsY<0.5){
1131                 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1132                 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1133                 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1134                 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1135
1136                 // fill mother hadron kinematics
1137                 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1138                 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1139                 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1140                 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1141               }
1142
1143               // ratio between pT of electron and pT of mother B or direct D hadron
1144               if(partMotherCopy->Pt()) {
1145                  fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1146                  if(eabsEta<0.9){
1147                    fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1148                  }
1149                  if(eabsY<0.5){
1150                    fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1151                  }
1152               }
1153               fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1154               if(eabsEta<0.9){
1155                 fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1156               }
1157               if(eabsY<0.5){
1158                 fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1159               }
1160               if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1161                 fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1162                 if(eabsEta<0.9){
1163                   fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1164                 }
1165                 if(eabsY<0.5){
1166                   fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1167                 }
1168               }
1169               else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1170                 fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1171                 if(eabsEta<0.9){
1172                   fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1173                 }
1174                 if(eabsY<0.5){
1175                   fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1176                 }
1177               }
1178               else {
1179                 fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1180                 if(eabsEta<0.9){
1181                   fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1182                 }
1183                 if(eabsY<0.5){
1184                   fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1185                 }
1186               }
1187
1188               //mj: to calculate D to e eta correlation to calculate total heavy quark cross section
1189               Int_t kLabel = partMotherCopy->GetFirstMother();
1190               Bool_t isGrandmaYes = kFALSE;
1191               Double_t gmrapidwstmp =0;
1192               if (!(kLabel < 0)){ // safety protection
1193                 if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel))))){
1194                   TParticle* grandMatmp = mctrack->Particle();
1195                   Int_t grandMaPDGtmp = grandMatmp->GetPdgCode();
1196                   if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kCharm || int(TMath::Abs(grandMaPDGtmp)/1000.) == kCharm ) isGrandmaYes = kTRUE;
1197                   if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(grandMaPDGtmp)/1000.) == kBeauty) isGrandmaYes = kTRUE;
1198                   gmrapidwstmp = AliHFEtools::GetRapidity(grandMatmp);
1199                 }
1200               }
1201
1202               Double_t mrapidwstmp = AliHFEtools::GetRapidity(partMotherCopy);
1203               Double_t eetawstmp = mcpart->Eta();
1204
1205               Double_t mrapidtmp = TMath::Abs(mrapidwstmp);
1206               Double_t eetatmp = TMath::Abs(eetawstmp);
1207
1208               fHistComm[iq][0].fEtaCorrD->Fill(mrapidwstmp,eetawstmp);
1209               if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD->Fill(gmrapidwstmp,eetawstmp);
1210               else fHistComm[iq][0].fEtaCorrGD->Fill(mrapidwstmp,eetawstmp);
1211
1212               if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1213               else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1214               else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1215               else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1216               if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1217                 fHistComm[iq][0].fEtaCorrDp->Fill(mrapidwstmp,eetawstmp);
1218                 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDp->Fill(gmrapidwstmp,eetawstmp);
1219                 else fHistComm[iq][0].fEtaCorrGDp->Fill(mrapidwstmp,eetawstmp);
1220                 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1221                 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1222                 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1223                 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1224               }
1225               else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1226                 fHistComm[iq][0].fEtaCorrD0->Fill(mrapidwstmp,eetawstmp);
1227                 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD0->Fill(gmrapidwstmp,eetawstmp);
1228                 else fHistComm[iq][0].fEtaCorrGD0->Fill(mrapidwstmp,eetawstmp);
1229                 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Dinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1230                 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Dineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1231                 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Doutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1232                 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Douteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1233               }
1234               else {
1235                 fHistComm[iq][0].fEtaCorrDrest->Fill(mrapidwstmp,eetawstmp);
1236                 if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDrest->Fill(gmrapidwstmp,eetawstmp);
1237                 else fHistComm[iq][0].fEtaCorrGDrest->Fill(mrapidwstmp,eetawstmp);
1238                 if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1239                 else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1240                 else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
1241                 else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
1242               }
1243
1244               // distance between electron production point and primary vertex
1245               fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1246               if(eabsEta<0.9){
1247                 fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1248               }
1249               if(eabsY<0.5){
1250                 fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1251               }
1252             }
1253          }
1254     } // end of if
1255 }
1256
1257 //____________________________________________________________________
1258 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
1259 {
1260   // decay electron kinematics
1261
1262   if (kquark != kCharm && kquark != kBeauty) {
1263     AliDebug(1, "This task is only for heavy quark QA, return\n");
1264     return;
1265   }
1266
1267   Int_t iq = kquark - kCharm;
1268   Bool_t isFinalOpenCharm = kFALSE;
1269
1270   if(!mcpart){
1271     AliDebug(1, "no mcparticle, return\n");
1272     return;
1273   }
1274
1275   if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
1276
1277   Double_t eabsEta = TMath::Abs(mcpart->Eta());
1278   Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
1279
1280   // mother
1281   Int_t iLabel = mcpart->GetMother();
1282   if (iLabel<0){
1283     AliDebug(1, "Stack label is negative, return\n");
1284     return;
1285   }
1286
1287   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
1288
1289   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1290   AliAODMCParticle *partMotherCopy = partMother;
1291   Int_t maPdgcode = partMother->GetPdgCode();
1292   Int_t maPdgcodeCopy = maPdgcode;
1293
1294   Bool_t isMotherDirectCharm = kFALSE;
1295   if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1296
1297     for (Int_t i=0; i<fNparents; i++){
1298        if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1299          isFinalOpenCharm = kTRUE;
1300        }
1301     }
1302     if (!isFinalOpenCharm) return;
1303
1304     for (Int_t i=1; i<fgkMaxIter; i++){
1305
1306        Int_t jLabel = partMother->GetMother();
1307        if (jLabel == -1){
1308          isMotherDirectCharm = kTRUE;
1309          break; // if there is no ancester
1310        }
1311        if (jLabel < 0){ // safety protection
1312          AliDebug(1, "Stack label is negative, return\n");
1313          return;
1314        }
1315
1316        // if there is an ancester
1317        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1318        Int_t grandMaPDG = grandMa->GetPdgCode();
1319
1320        for (Int_t j=0; j<fNparents; j++){
1321           if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1322
1323             if (kquark == kCharm) return;
1324             // fill electron kinematics
1325             fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1326             fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1327             fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1328             fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1329
1330             // fill mother hadron kinematics
1331             fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1332             fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1333             fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1334             fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1335
1336             if(eabsEta<0.9){
1337               // fill electron kinematics
1338               fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1339               fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1340               fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1341               fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1342
1343               // fill mother hadron kinematics
1344               fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1345               fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1346               fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1347               fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1348             }
1349             if(eabsY<0.5){
1350               // fill electron kinematics
1351               fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1352               fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1353               fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1354               fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1355
1356               // fill mother hadron kinematics
1357               fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1358               fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1359               fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1360               fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1361             }
1362
1363             return;
1364           }
1365        }
1366
1367        partMother = grandMa;
1368     } // end of iteration
1369   } // end of if
1370   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1371     for (Int_t i=0; i<fNparents; i++){
1372        if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1373
1374          // fill electron kinematics
1375          fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1376          fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1377          fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1378          fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1379
1380          // fill mother hadron kinematics
1381          fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1382          fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1383          fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1384          fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1385
1386          if(eabsEta<0.9){
1387            // fill electron kinematics
1388            fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1389            fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1390            fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1391            fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1392
1393            // fill mother hadron kinematics
1394            fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1395            fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1396            fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1397            fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1398          }
1399          if(eabsY<0.5){
1400            // fill electron kinematics
1401            fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1402            fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1403            fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1404            fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1405
1406            // fill mother hadron kinematics
1407            fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1408            fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1409            fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1410            fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1411          }
1412
1413        }
1414     }
1415   } // end of if
1416
1417 }
1418
1419 //__________________________________________
1420 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
1421 {
1422        // find mother pdg code and label
1423
1424        if (motherlabel < 0) {
1425          AliDebug(1, "Stack label is negative, return\n");
1426          return;
1427        }
1428        AliMCParticle *mctrack = NULL;
1429        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1430        TParticle *heavysMother = mctrack->Particle();
1431        motherpdg = heavysMother->GetPdgCode();
1432        grandmotherlabel = heavysMother->GetFirstMother();
1433        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
1434 }
1435
1436 //__________________________________________
1437 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1438 {
1439        // mothertype -1 means this heavy quark coming from hard vertex
1440
1441        AliMCParticle *mctrack1 = NULL;
1442        AliMCParticle *mctrack2 = NULL;
1443        if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1444        if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1445        TParticle *afterinitialrad1  = mctrack1->Particle();
1446        TParticle *afterinitialrad2  = mctrack2->Particle();
1447
1448        motherlabel = -1;
1449
1450        if (TMath::Abs(afterinitialrad1->GetPdgCode()) == fgkGluon && TMath::Abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1451          AliDebug(1,"heavy from gluon gluon pair creation!\n");
1452          mothertype = -1;
1453          motherID = fgkGluon;
1454        }
1455        else if (TMath::Abs(afterinitialrad1->GetPdgCode()) == kquark || TMath::Abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1456          AliDebug(1,"heavy from flavor exitation!\n");
1457          mothertype = -1;
1458          motherID = kquark;
1459        }
1460        else if  (TMath::Abs(afterinitialrad1->GetPdgCode()) == TMath::Abs(afterinitialrad2->GetPdgCode())){
1461          AliDebug(1,"heavy from q-qbar pair creation!\n");
1462          mothertype = -1;
1463          motherID = 1;
1464        }
1465        else {
1466          AliDebug(1,"something strange!\n");
1467          mothertype = -999;
1468          motherlabel = -999;
1469          motherID = -999;
1470        }
1471 }
1472
1473 //__________________________________________
1474 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1475 {
1476        // mothertype -2 means this heavy quark coming from initial state
1477
1478        AliMCParticle *mctrack = NULL;
1479        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
1480          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1481          TParticle *heavysMother = mctrack->Particle();
1482          motherID = heavysMother->GetPdgCode();
1483          mothertype = -2; // there is mother before initial state radiation
1484          motherlabel = inputmotherlabel;
1485          AliDebug(1,"initial parton shower! \n");
1486
1487          return kTRUE;
1488        }
1489
1490        return kFALSE;
1491 }
1492
1493 //__________________________________________
1494 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1495 {
1496        // mothertype 2 means this heavy quark coming from final state
1497
1498        AliMCParticle *mctrack = NULL;
1499        if (inputmotherlabel > 5){ // mother exist after hard scattering
1500          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1501          TParticle *heavysMother = mctrack->Particle();
1502          motherID = heavysMother->GetPdgCode();
1503          mothertype = 2; //
1504          motherlabel = inputmotherlabel;
1505          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1506
1507          return kTRUE;
1508        }
1509        return kFALSE;
1510 }
1511
1512 //__________________________________________
1513 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1514 {
1515       // mark strange behavior
1516
1517        mothertype = -888;
1518        motherlabel = -888;
1519        motherID = -888;
1520        AliDebug(1,"something strange!\n");
1521 }
1522
1523 //__________________________________________
1524 Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart)
1525 {
1526   // decay particle's origin
1527
1528   //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1529
1530   Int_t origin = -1;
1531   Bool_t isFinalOpenCharm = kFALSE;
1532
1533   if(!mcpart){
1534     AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1535     return -1;
1536   }
1537
1538   // mother
1539   // Information not in the base class, cast necessary
1540   Int_t iLabel = GetMother(mcpart);
1541   if (iLabel<0){
1542     AliDebug(1, "Stack label is negative, return\n");
1543     return -1;
1544   }
1545
1546   const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
1547   Int_t maPdgcode = partMother->PdgCode();
1548
1549   // if the mother is charmed hadron
1550   if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1551
1552     for (Int_t i=0; i<fNparents; i++){
1553        if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1554          isFinalOpenCharm = kTRUE;
1555        }
1556     }
1557     if (!isFinalOpenCharm) return -1;
1558
1559     // iterate until you find B hadron as a mother or become top ancester
1560     for (Int_t i=1; i<fgkMaxIter; i++){
1561
1562        Int_t jLabel = GetMother(partMother);
1563        if (jLabel == -1){
1564          origin = kDirectCharm;
1565          return origin;
1566        }
1567        if (jLabel < 0){ // safety protection
1568          AliDebug(1, "Stack label is negative, return\n");
1569          return -1;
1570        }
1571
1572        // if there is an ancester
1573        const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
1574        Int_t grandMaPDG = grandMa->PdgCode();
1575
1576        for (Int_t j=0; j<fNparents; j++){
1577           if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1578             origin = kBeautyCharm;
1579             return origin;
1580           }
1581        }
1582
1583        partMother = grandMa;
1584     } // end of iteration
1585   } // end of if
1586   else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
1587     for (Int_t i=0; i<fNparents; i++){
1588        if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1589          origin = kDirectBeauty;
1590          return origin;
1591        }
1592     }
1593   } // end of if
1594   else if ( TMath::Abs(maPdgcode) == 22 ) {
1595     origin = kGamma;
1596     return origin;
1597   } // end of if
1598   else if ( TMath::Abs(maPdgcode) == 111 ) {
1599     origin = kPi0;
1600     return origin;
1601   } // end of if
1602
1603   return origin;
1604 }
1605
1606 //__________________________________________
1607 Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
1608 {
1609   // decay particle's origin
1610
1611   //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1612
1613   Int_t origin = -1;
1614   Bool_t isFinalOpenCharm = kFALSE;
1615
1616   if(!mcpart){
1617     AliDebug(1, "no mcparticle, return\n");
1618     return -1;
1619   }
1620
1621   Int_t iLabel = mcpart->GetFirstMother();
1622   if (iLabel<0){
1623     AliDebug(1, "Stack label is negative, return\n");
1624     return -1;
1625   }
1626
1627   AliMCParticle *mctrack = NULL;
1628   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1629   TParticle *partMother = mctrack->Particle();
1630   Int_t maPdgcode = partMother->GetPdgCode();
1631
1632    // if the mother is charmed hadron
1633    if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
1634
1635      for (Int_t i=0; i<fNparents; i++){
1636         if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1637           isFinalOpenCharm = kTRUE;
1638         }
1639      }
1640      if (!isFinalOpenCharm) return -1;
1641
1642      // iterate until you find B hadron as a mother or become top ancester
1643      for (Int_t i=1; i<fgkMaxIter; i++){
1644
1645         Int_t jLabel = partMother->GetFirstMother();
1646         if (jLabel == -1){
1647           origin = kDirectCharm;
1648           return origin;
1649         }
1650         if (jLabel < 0){ // safety protection
1651           AliDebug(1, "Stack label is negative, return\n");
1652           return -1;
1653         }
1654
1655         // if there is an ancester
1656         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1657         TParticle* grandMa = mctrack->Particle();
1658         Int_t grandMaPDG = grandMa->GetPdgCode();
1659
1660         for (Int_t j=0; j<fNparents; j++){
1661            if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1662              origin = kBeautyCharm;
1663              return origin;
1664            }
1665         }
1666
1667         partMother = grandMa;
1668      } // end of iteration
1669    } // end of if
1670    else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
1671      for (Int_t i=0; i<fNparents; i++){
1672         if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1673           origin = kDirectBeauty;
1674           return origin;
1675         }
1676      }
1677    } // end of if
1678    else if ( TMath::Abs(maPdgcode) == 22 ) {
1679      origin = kGamma;
1680      return origin;
1681    } // end of if
1682    else if ( TMath::Abs(maPdgcode) == 111 ) {
1683      origin = kPi0;
1684      return origin;
1685    } // end of if
1686    else origin = kElse;
1687
1688    return origin;
1689 }
1690
1691 //__________________________________________
1692 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1693 {
1694   // decay particle's origin
1695
1696   if(!mcpart){
1697     AliDebug(1, "no mcparticle, return\n");
1698     return -1;
1699   }
1700
1701   if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1702
1703   Int_t origin = -1;
1704   Bool_t isFinalOpenCharm = kFALSE;
1705
1706   Int_t iLabel = mcpart->GetFirstMother();
1707   if (iLabel<0){
1708     AliDebug(1, "Stack label is negative, return\n");
1709     return -1;
1710   }
1711
1712   AliMCParticle *mctrack = NULL;
1713   Int_t tmpMomLabel=0;
1714   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1715   TParticle *partMother = mctrack->Particle();
1716   TParticle *partMotherCopy = mctrack->Particle();
1717   Int_t maPdgcode = partMother->GetPdgCode();
1718   Int_t grmaPdgcode;
1719   Int_t ggrmaPdgcode;
1720
1721    // if the mother is charmed hadron
1722
1723    if(TMath::Abs(maPdgcode)==443){ // J/spi
1724       Int_t jLabel = partMother->GetFirstMother();
1725       if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))){
1726         TParticle* grandMa = mctrack->Particle();
1727         Int_t grandMaPDG = grandMa->GetPdgCode();
1728         if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
1729       }
1730       return kJpsi;
1731    }
1732    else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
1733
1734      for (Int_t i=0; i<fNparents; i++){
1735         if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1736           isFinalOpenCharm = kTRUE;
1737         }
1738      }
1739      if (!isFinalOpenCharm) return -1;
1740
1741      // iterate until you find B hadron as a mother or become top ancester
1742      for (Int_t i=1; i<fgkMaxIter; i++){
1743
1744         Int_t jLabel = partMother->GetFirstMother();
1745         if (jLabel == -1){
1746           return kDirectCharm;
1747         }
1748         if (jLabel < 0){ // safety protection
1749           AliDebug(1, "Stack label is negative, return\n");
1750           return -1;
1751         }
1752
1753         // if there is an ancester
1754         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1755         TParticle* grandMa = mctrack->Particle();
1756         Int_t grandMaPDG = grandMa->GetPdgCode();
1757
1758         for (Int_t j=0; j<fNparents; j++){
1759            if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1760              return kBeautyCharm;
1761            }
1762         }
1763
1764         partMother = grandMa;
1765      } // end of iteration
1766    } // end of if
1767    else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
1768      for (Int_t i=0; i<fNparents; i++){
1769         if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1770           return kDirectBeauty;
1771         }
1772      }
1773    } // end of if
1774    else if ( TMath::Abs(maPdgcode) == 22 ) { //conversion
1775
1776      tmpMomLabel = partMotherCopy->GetFirstMother();
1777      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1778      partMother = mctrack->Particle();
1779      maPdgcode = partMother->GetPdgCode();
1780
1781      // check if the ligth meson is the decay product of heavy mesons
1782      tmpMomLabel = partMother->GetFirstMother();
1783      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1784       partMother = mctrack->Particle();
1785       grmaPdgcode = partMother->GetPdgCode();
1786
1787       if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1788       if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1789       if ( TMath::Abs(grmaPdgcode) == 221 ) return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
1790
1791       tmpMomLabel = partMother->GetFirstMother();
1792       if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1793        partMother = mctrack->Particle();
1794        ggrmaPdgcode = partMother->GetPdgCode();
1795
1796        if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1797        if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1798       }
1799      }
1800
1801      if ( TMath::Abs(maPdgcode) == 111 ) {
1802        return kGammaPi0;
1803      }
1804      else if ( TMath::Abs(maPdgcode) == 221 ) {
1805        return kGammaEta;
1806      }
1807      else if ( TMath::Abs(maPdgcode) == 223 ) {
1808        return kGammaOmega;
1809      }
1810      else if ( TMath::Abs(maPdgcode) == 333 ) {
1811        return kGammaPhi;
1812      }
1813      else if ( TMath::Abs(maPdgcode) == 331 ) {
1814        return kGammaEtaPrime;
1815      }
1816      else if ( TMath::Abs(maPdgcode) == 113 ) {
1817        return kGammaRho0;
1818      }
1819      else origin = kElse;
1820      return origin;
1821
1822    }
1823    else {
1824
1825      // check if the ligth meson is the decay product of heavy mesons
1826      tmpMomLabel = partMotherCopy->GetFirstMother();
1827      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1828       partMother = mctrack->Particle();
1829       grmaPdgcode = partMother->GetPdgCode();
1830
1831       if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1832       if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1833       if ( TMath::Abs(grmaPdgcode) == 221 ) return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
1834
1835       tmpMomLabel = partMother->GetFirstMother();
1836       if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
1837        partMother = mctrack->Particle();
1838        ggrmaPdgcode = partMother->GetPdgCode();
1839
1840        if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
1841        if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
1842       }
1843      }
1844
1845      if ( TMath::Abs(maPdgcode) == 111 ) {
1846        return kPi0;
1847      }
1848      else if ( TMath::Abs(maPdgcode) == 221 ) {
1849        return kEta;
1850      }
1851      else if ( TMath::Abs(maPdgcode) == 223 ) {
1852        return kOmega;
1853      }
1854      else if ( TMath::Abs(maPdgcode) == 333 ) {
1855        return kPhi;
1856      }
1857      else if ( TMath::Abs(maPdgcode) == 331 ) {
1858        return kEtaPrime;
1859      }
1860      else if ( TMath::Abs(maPdgcode) == 113 ) {
1861        return kRho0;
1862      }
1863      else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
1864        return kKe3;
1865      }
1866      else{
1867       origin = kElse;
1868      }
1869
1870    }
1871    return origin;
1872 }
1873
1874 //__________________________________________
1875 Int_t AliHFEmcQA::GetElecSource(const AliVParticle * const mctrack)
1876 {
1877   //
1878   // decay particle's origin
1879   //
1880
1881   if(!mctrack){
1882     AliDebug(1, "no mcparticle, return\n");
1883     return -1;
1884   }
1885
1886   TClass *type = mctrack->IsA();
1887
1888   if(type == AliMCParticle::Class()) {
1889     const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
1890     if(esdmc){
1891       TParticle *mcpart =  esdmc->Particle();
1892       return GetElecSource(mcpart);
1893     }
1894     else return -1;
1895   }
1896   if(type == AliAODMCParticle::Class()) {
1897     const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(mctrack);
1898     if(aodmc){
1899       return GetElecSource(aodmc);
1900     }
1901     else return -1;
1902   }
1903   return -1;
1904 }
1905 //__________________________________________
1906 Int_t AliHFEmcQA::GetElecSource(const AliAODMCParticle * const mcpart)
1907 {
1908   //
1909   // Function for AliAODMCParticle
1910   //
1911
1912   if (!mcpart) return -1;
1913   if (!fMCArray) return -1;
1914
1915   if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1916
1917   Int_t origin = -1;
1918   Bool_t isFinalOpenCharm = kFALSE;
1919
1920   Int_t iLabel = mcpart->GetMother();
1921   if ((iLabel<0) || (iLabel>=fMCArray->GetEntriesFast())){
1922     AliDebug(1, "label is out of range, return\n");
1923     return -1;
1924   }
1925
1926   AliAODMCParticle *mctrack = NULL; // will change all the time
1927   Int_t tmpMomLabel=0;
1928   if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return -1;
1929   AliAODMCParticle *partMother = mctrack;
1930   AliAODMCParticle *partMotherCopy = mctrack;
1931   Int_t maPdgcode = mctrack->GetPdgCode();
1932   Int_t grmaPdgcode;
1933   Int_t ggrmaPdgcode;
1934
1935   // if the mother is charmed hadron
1936
1937   if(TMath::Abs(maPdgcode)==443){
1938     //
1939     // J/spi
1940     //
1941     Int_t jLabel = partMother->GetMother();
1942     if ((jLabel>=0) && (jLabel<fMCArray->GetEntriesFast())){
1943       if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))){
1944         Int_t grandMaPDG = mctrack->GetPdgCode();
1945         if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) {
1946           return kB2Jpsi;
1947         }
1948       }
1949     }
1950     return kJpsi;
1951   }
1952   else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
1953     //
1954     // charm
1955     //
1956     for (Int_t i=0; i<fNparents; i++){
1957       if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
1958         isFinalOpenCharm = kTRUE;
1959       }
1960     }
1961     if (!isFinalOpenCharm) {
1962       return -1;
1963     }
1964
1965     // iterate until you find B hadron as a mother or become top ancester
1966     for (Int_t i=1; i<fgkMaxIter; i++){
1967
1968       Int_t jLabel = partMother->GetMother();
1969       if (jLabel == -1){
1970         return kDirectCharm;
1971       }
1972       if ((jLabel<0) || (jLabel>=fMCArray->GetEntriesFast())){
1973         AliDebug(1, "Stack label is negative, return\n");
1974         return -1;
1975       }
1976
1977       // if there is an ancester
1978       if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))) {
1979         return -1;
1980       }
1981       Int_t grandMaPDG = mctrack->GetPdgCode();
1982       for (Int_t j=0; j<fNparents; j++){
1983         if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
1984           return kBeautyCharm;
1985         }
1986       }
1987       partMother = mctrack;
1988     } // end of iteration
1989
1990   } // end of if
1991   else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
1992     //
1993     // beauty
1994     //
1995     for (Int_t i=0; i<fNparents; i++){
1996       if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
1997         return kDirectBeauty;
1998       }
1999     }
2000   } // end of if
2001   else if ( TMath::Abs(maPdgcode) == 22 ) {
2002     //
2003     //conversion
2004     //
2005     tmpMomLabel = partMotherCopy->GetMother();
2006     if((tmpMomLabel<0) || (tmpMomLabel>=fMCArray->GetEntriesFast())) {
2007       return -1;
2008     }
2009     if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2010       return -1;
2011     }
2012     partMother = mctrack;
2013     maPdgcode = partMother->GetPdgCode();
2014
2015     // check if the ligth meson is the decay product of heavy mesons
2016     tmpMomLabel = partMother->GetMother();
2017     if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2018       if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2019         partMother = mctrack;
2020         grmaPdgcode = partMother->GetPdgCode();
2021
2022         if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {
2023           return kGammaB2M;
2024         }
2025         if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {
2026           return kGammaD2M;
2027         }
2028         if ( TMath::Abs(grmaPdgcode) == 221 ) {
2029           return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
2030         }
2031
2032         tmpMomLabel = partMother->GetMother();
2033         if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2034           if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2035             partMother = mctrack;
2036             ggrmaPdgcode = partMother->GetPdgCode();
2037
2038             if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
2039               return kGammaB2M;
2040             }
2041             if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
2042               return kGammaD2M;
2043             }
2044            }
2045         }
2046       }
2047     }
2048
2049     if ( TMath::Abs(maPdgcode) == 111 ) {
2050       return kGammaPi0;
2051     }
2052     else if ( TMath::Abs(maPdgcode) == 221 ) {
2053       return kGammaEta;
2054     }
2055     else if ( TMath::Abs(maPdgcode) == 223 ) {
2056       return kGammaOmega;
2057     }
2058     else if ( TMath::Abs(maPdgcode) == 333 ) {
2059       return kGammaPhi;
2060     }
2061     else if ( TMath::Abs(maPdgcode) == 331 ) {
2062       return kGammaEtaPrime;
2063     }
2064     else if ( TMath::Abs(maPdgcode) == 113 ) {
2065       return kGammaRho0;
2066     }
2067     else origin = kElse;
2068
2069     return origin;
2070
2071   }
2072    else {
2073      //
2074      // check if the ligth meson is the decay product of heavy mesons
2075      //
2076      tmpMomLabel = partMotherCopy->GetMother();
2077      if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2078        if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2079          partMother = mctrack;
2080          grmaPdgcode = partMother->GetPdgCode();
2081
2082          if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {
2083            return kGammaB2M;
2084          }
2085          if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {
2086            return kGammaD2M;
2087          }
2088          if ( TMath::Abs(grmaPdgcode) == 221 ) {
2089            return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
2090          }
2091
2092          tmpMomLabel = partMother->GetMother();
2093          if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
2094            if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
2095              partMother = mctrack;
2096              ggrmaPdgcode = partMother->GetPdgCode();
2097
2098              if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
2099                return kGammaB2M;
2100              }
2101              if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
2102                return kGammaD2M;
2103              }
2104            }
2105          }
2106        }
2107
2108        if ( TMath::Abs(maPdgcode) == 111 ) {
2109          return kPi0;
2110        }
2111        else if ( TMath::Abs(maPdgcode) == 221 ) {
2112          return kEta;
2113        }
2114        else if ( TMath::Abs(maPdgcode) == 223 ) {
2115          return kOmega;
2116        }
2117        else if ( TMath::Abs(maPdgcode) == 333 ) {
2118          return kPhi;
2119        }
2120        else if ( TMath::Abs(maPdgcode) == 331 ) {
2121          return kEtaPrime;
2122        }
2123        else if ( TMath::Abs(maPdgcode) == 113 ) {
2124          return kRho0;
2125        }
2126        else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
2127          return kKe3;
2128        }
2129        else{
2130          origin = kElse;
2131        }
2132      }
2133    }
2134   return origin;
2135 }
2136 //__________________________________________
2137 Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
2138   //
2139   // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
2140   //
2141   AliMCParticle *mctrackmother = NULL;
2142   Double_t weightElecBg = 0.;
2143   Double_t mesonPt = 0.;
2144   Double_t bgcategory = 0.;
2145   Int_t mArr = -1;
2146   Int_t mesonID = GetElecSource(mctrack->Particle());
2147   if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0;                //pion
2148   else if(mesonID==kGammaEta || mesonID==kEta) mArr=1;           //eta
2149   else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2;       //omega
2150   else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3;           //phi
2151   else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
2152   else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
2153
2154   Double_t datamc[25]={-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999, -999, -999, -999, -999, -999, -999, -999, -999, -999};
2155   Double_t xr[3]={-999,-999,-999};
2156   datamc[0] = mesonID;
2157   datamc[17] = mctrack->Pt(); //electron pt
2158   datamc[18] = mctrack->Eta(); //electron eta
2159
2160   mctrack->XvYvZv(xr);
2161   datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2162   datamc[10] = xr[2];
2163
2164   TParticle *mcpart = mctrack->Particle();
2165   if(mcpart){
2166     datamc[14] = mcpart->GetUniqueID();
2167   }
2168
2169   if(!(mArr<0)){
2170      if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering
2171         Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
2172         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2173           glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
2174           if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2175             mesonPt = mctrackmother->Pt(); //meson pt
2176             bgcategory = 1.;
2177             datamc[1] = bgcategory;
2178             datamc[2] = mesonPt;
2179             mctrackmother->XvYvZv(xr);
2180             datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2181             datamc[12] = xr[2];
2182
2183             mcpart = mctrackmother->Particle();
2184             if(mcpart){
2185               datamc[15] = mcpart->GetUniqueID();
2186             }
2187             if(glabel>fMCEvent->GetNumberOfPrimaries()) {
2188               bgcategory = 2.;
2189               datamc[1] = bgcategory;
2190               //printf("I should be gamma meson = %d  mesonlabel= %d  NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries());
2191               glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2192               if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2193                 datamc[3]=mctrackmother->PdgCode();
2194                 datamc[4]=mctrackmother->Pt();
2195                 if(TMath::Abs(mctrackmother->PdgCode())==310){
2196                   bgcategory = 3.;
2197                   glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
2198                   if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2199                     datamc[5]=mctrackmother->PdgCode();
2200                     glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2201                     if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2202                        datamc[6]=mctrackmother->PdgCode();
2203                     }
2204                   }
2205                 }
2206               }
2207             }
2208           }
2209         }
2210      }
2211      else{ // nonHFE except for the conversion electron
2212         Int_t glabel=TMath::Abs(mctrack->GetMother());
2213         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2214           mesonPt=mctrackmother->Pt(); //meson pt
2215           bgcategory = -1.;
2216           datamc[1] = bgcategory;
2217           datamc[2] = mesonPt;
2218           mctrackmother->XvYvZv(xr);
2219           datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
2220           datamc[12] = xr[2];
2221
2222           mcpart = mctrackmother->Particle();
2223           if(mcpart){
2224             datamc[15] = mcpart->GetUniqueID();
2225           }
2226           if(glabel>fMCEvent->GetNumberOfPrimaries()) {
2227             bgcategory = -2.;
2228             datamc[1] = bgcategory;
2229             glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2230             if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2231               datamc[3]=mctrackmother->PdgCode();
2232               datamc[4]=mctrackmother->Pt();
2233               if(TMath::Abs(mctrackmother->PdgCode())==310){
2234                bgcategory = -3.;
2235                glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
2236                if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2237                  datamc[5]=mctrackmother->PdgCode();
2238                  glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
2239                  if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
2240                    datamc[6]=mctrackmother->PdgCode();
2241                  }
2242                }
2243               }
2244             }
2245           }
2246         }
2247      }
2248
2249
2250      if(fIsPbPb){
2251        Int_t centBin = GetWeightCentralityBin(fPerCentrality);
2252        if(centBin < 0)return 0.;
2253        weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
2254        for(int ii=0; ii<kBgPtBins; ii++){
2255          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2256            weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
2257            break;
2258          }
2259        }
2260      }
2261      else{
2262        weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
2263        for(int ii=0; ii<kBgPtBins; ii++){
2264          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
2265            weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
2266            break;
2267          }
2268        }
2269      }
2270   }
2271
2272   datamc[13] = weightElecBg;
2273   datamc[16] = Double_t(fContainerStep);
2274
2275   datamc[7] = fHfeImpactR;
2276   datamc[8] = fHfeImpactnsigmaR;
2277
2278   datamc[19] = fRecPt;
2279   datamc[20] = fRecEta;
2280   datamc[21] = fRecPhi;
2281   datamc[22] = fLyrhit;
2282   datamc[23] = fLyrstat;
2283   datamc[24] = fCentrality;
2284
2285   if(fIsDebugStreamerON && fTreeStream){
2286    if(!iBgLevel){
2287     (*fTreeStream)<<"nonhfeQA"<<
2288         "mesonID="<<datamc[0]<<
2289         "bgcategory="<<datamc[1]<<
2290         "mesonPt="<<datamc[2]<<
2291         "mesonMomPdg="<<datamc[3]<<
2292         "mesonMomPt="<<datamc[4]<<
2293         "mesonGMomPdg="<<datamc[5]<<
2294         "mesonGGMomPdg="<<datamc[6]<<
2295         "eIPAbs="<<datamc[7]<<
2296         "eIPSig="<<datamc[8]<<
2297         "eR="<<datamc[9]<<
2298         "eZ="<<datamc[10]<<
2299         "mesonR="<<datamc[11]<<
2300         "mesonZ="<<datamc[12]<<
2301         "weightElecBg="<<datamc[13]<<
2302         "eUniqID="<<datamc[14]<<
2303         "mesonUniqID="<<datamc[15]<<
2304         "containerStep="<<datamc[16]<<
2305         "emcpt="<<datamc[17]<<
2306         "emceta="<<datamc[18]<<
2307         "erecpt="<<datamc[19]<<
2308         "ereceta="<<datamc[20]<<
2309         "erecphi="<<datamc[21]<<
2310         "itshit="<<datamc[22]<<
2311         "itsstat="<<datamc[23]<<
2312         "centrality="<<datamc[24]
2313         << "\n";
2314    }
2315   }
2316
2317   Double_t returnval = bgcategory*weightElecBg;
2318   if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
2319
2320   return returnval;
2321 }
2322
2323 //__________________________________________
2324 Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart){
2325   //
2326   // Wrapper to get the mother label
2327   //
2328   Int_t label = -1;
2329   TClass *type = mcpart->IsA();
2330   if(type == AliMCParticle::Class()){
2331     // ESD analysis
2332     const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
2333     label = emcpart->GetMother();
2334   } else if(type == AliAODMCParticle::Class()){
2335     // AOD analysis
2336     const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
2337     label = amcpart->GetMother();
2338   }
2339   return label;
2340 }
2341 //__________________________________________
2342 Int_t AliHFEmcQA::GetWeightCentralityBin(const Float_t percentile) const {
2343   //
2344   //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
2345   //
2346
2347   Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
2348   Int_t bin = -1;
2349   for(Int_t ibin = 0; ibin < 11; ibin++){
2350     if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
2351       bin = ibin;
2352       break;
2353     }
2354   }
2355   return bin;
2356 }
2357 //__________________________________________
2358 void AliHFEmcQA::AliHists::FillList(TList *l) const {
2359   //
2360   // Fill Histos into a list for output
2361   //
2362   if(fPdgCode) l->Add(fPdgCode);
2363   if(fPt) l->Add(fPt);
2364   if(fY) l->Add(fY);
2365   if(fEta) l->Add(fEta);
2366 }
2367
2368 //__________________________________________
2369 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
2370   //
2371   // Fill Histos into a list for output
2372   //
2373   if(fNq) l->Add(fNq);
2374   if(fProcessID) l->Add(fProcessID);
2375   if(fePtRatio) l->Add(fePtRatio);
2376   if(fPtCorr) l->Add(fPtCorr);
2377   if(fPtCorrDp) l->Add(fPtCorrDp);
2378   if(fPtCorrD0) l->Add(fPtCorrD0);
2379   if(fPtCorrDrest) l->Add(fPtCorrDrest);
2380
2381   if(fPtCorrDinein) l->Add(fPtCorrDinein);
2382   if(fPtCorrDineout) l->Add(fPtCorrDineout);
2383   if(fPtCorrDoutein) l->Add(fPtCorrDoutein);
2384   if(fPtCorrDouteout) l->Add(fPtCorrDouteout);
2385   if(fPtCorrDpDinein) l->Add(fPtCorrDpDinein);
2386   if(fPtCorrDpDineout) l->Add(fPtCorrDpDineout);
2387   if(fPtCorrDpDoutein) l->Add(fPtCorrDpDoutein);
2388   if(fPtCorrDpDouteout) l->Add(fPtCorrDpDouteout);
2389   if(fPtCorrD0Dinein) l->Add(fPtCorrD0Dinein);
2390   if(fPtCorrD0Dineout) l->Add(fPtCorrD0Dineout);
2391   if(fPtCorrD0Doutein) l->Add(fPtCorrD0Doutein);
2392   if(fPtCorrD0Douteout) l->Add(fPtCorrD0Douteout);
2393   if(fPtCorrDrestDinein) l->Add(fPtCorrDrestDinein);
2394   if(fPtCorrDrestDineout) l->Add(fPtCorrDrestDineout);
2395   if(fPtCorrDrestDoutein) l->Add(fPtCorrDrestDoutein);
2396   if(fPtCorrDrestDouteout) l->Add(fPtCorrDrestDouteout);
2397
2398   if(fEtaCorrD) l->Add(fEtaCorrD);
2399   if(fEtaCorrDp) l->Add(fEtaCorrDp);
2400   if(fEtaCorrD0) l->Add(fEtaCorrD0);
2401   if(fEtaCorrDrest) l->Add(fEtaCorrDrest);
2402
2403   if(fEtaCorrGD) l->Add(fEtaCorrGD);
2404   if(fEtaCorrGDp) l->Add(fEtaCorrGDp);
2405   if(fEtaCorrGD0) l->Add(fEtaCorrGD0);
2406   if(fEtaCorrGDrest) l->Add(fEtaCorrGDrest);
2407
2408   if(fEtaCorrB) l->Add(fEtaCorrB);
2409   if(fEtaCorrGB) l->Add(fEtaCorrGB);
2410   if(fPtCorrBinein) l->Add(fPtCorrBinein);
2411   if(fPtCorrBineout) l->Add(fPtCorrBineout);
2412   if(fPtCorrBoutein) l->Add(fPtCorrBoutein);
2413   if(fPtCorrBouteout) l->Add(fPtCorrBouteout);
2414
2415   if(fDePtRatio) l->Add(fDePtRatio);
2416   if(feDistance) l->Add(feDistance);
2417   if(fDeDistance) l->Add(fDeDistance);
2418 }
2419