PWGHF/hfe/AliHFEmcQA.cxx

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