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