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 //_______________________________________________________________________________________________
 125 AliHFEmcQA&
 126 AliHFEmcQA::operator=(const AliHFEmcQA &)
 127 {
 128   // Assignment operator
 129
 130   AliInfo("Not yet implemented.");
 131   return *this;
 132 }
 133
 134 //_______________________________________________________________________________________________
 135 AliHFEmcQA::~AliHFEmcQA()
 136 {
 137         // Destructor
 138
 139         if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
 140         AliInfo("Analysis Done.");
 141 }
 142
 143 //_______________________________________________________________________________________________
 144 void AliHFEmcQA::PostAnalyze() const
 145 {
 146         //
 147         // Post analysis
 148         //
 149 }
 150
 151 //_______________________________________________________________________________________________
 152 void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
 153 {
 154    //
 155    // copy background weighting factors into data member
 156    //
 157
 158   memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 159   memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
 160 }
 161
 162 //__________________________________________
 163 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
 164 {
 165   //
 166   // make default histograms
 167   //
 168
 169   if(!qaList) return;
 170
 171   fQAhistos = qaList;
 172   fQAhistos->SetName("MCqa");
 173
 174   CreateHistograms(AliHFEmcQA::kCharm);               // create histograms for charm
 175   CreateHistograms(AliHFEmcQA::kBeauty);              // create histograms for beauty
 176   CreateHistograms(AliHFEmcQA::kOthers);              // create histograms for beauty
 177
 178 // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
 179   const Int_t nbspecies = 9;
 180   TString kDspecies[nbspecies];
 181   kDspecies[0]="411";   //D+
 182   kDspecies[1]="421";   //D0
 183   kDspecies[2]="431";   //Ds+
 184   kDspecies[3]="4122";  //Lambdac+
 185   kDspecies[4]="4132";  //Ksic0
 186   kDspecies[5]="4232";  //Ksic+
 187   kDspecies[6]="4332";  //OmegaC0
 188   kDspecies[7]="413";   //D*(2010)+
 189   kDspecies[8]="423";   //D*(2007)0
 190
 191   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 192   Int_t iBin[2];
 193   iBin[0] = 44; // bins in pt for log binning
 194   iBin[1] = 23; // bins in pt for pi0 measurement binning
 195   Double_t* binEdges[1];
 196   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 197
 198   // bin size is chosen to consider ALICE D measurement
 199   const Int_t nptbins = 15;
 200   const Int_t nybins = 9;
 201   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
 202   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
 203   TString hname;
 204   for (Int_t iDmeson=0; iDmeson<nbspecies; iDmeson++){
 205      hname = "Dmeson"+kDspecies[iDmeson];
 206      fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",nptbins,xbins,nybins,ybins);
 207      if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
 208   }
 209
 210   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
 211
 212   Int_t kNcent;
 213   if(fIsPbPb) kNcent=11;
 214   else
 215   {
 216       if(fIsppMultiBin) kNcent=8;
 217       else kNcent = 1;
 218   }
 219
 220   fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
 221
 222   for(Int_t centbin=0; centbin<kNcent; centbin++)
 223   {
 224       fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
 225       fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
 226       fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
 227       fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
 228       fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
 229       fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
 230
 231       fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 232       fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 233       fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 234       fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 235       fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 236       fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 237       fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 238       fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 239
 240       fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 241       fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 242       fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 243       fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 244       fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 245       fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 246   }
 247
 248   fQAhistos->Add(fMCQACollection->GetList());
 249
 250   if(!fTreeStream && fIsDebugStreamerON){
 251    fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
 252   }
 253
 254 }
 255
 256 //__________________________________________
 257 void AliHFEmcQA::CreateHistograms(const Int_t kquark)
 258 {
 259   // create histograms
 260
 261   if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
 262     AliDebug(1, "This task is only for heavy quark QA, return\n");
 263     return;
 264   }
 265   Int_t iq = kquark - kCharm;
 266
 267   TString kqTypeLabel[fgkqType];
 268   if (kquark == kCharm){
 269     kqTypeLabel[kQuark]="c";
 270     kqTypeLabel[kantiQuark]="cbar";
 271     kqTypeLabel[kHadron]="cHadron";
 272     kqTypeLabel[keHadron]="ceHadron";
 273     kqTypeLabel[kDeHadron]="nullHadron";
 274     kqTypeLabel[kElectron]="ce";
 275     kqTypeLabel[kElectron2nd]="nulle";
 276   } else if (kquark == kBeauty){
 277     kqTypeLabel[kQuark]="b";
 278     kqTypeLabel[kantiQuark]="bbar";
 279     kqTypeLabel[kHadron]="bHadron";
 280     kqTypeLabel[keHadron]="beHadron";
 281     kqTypeLabel[kDeHadron]="bDeHadron";
 282     kqTypeLabel[kElectron]="be";
 283     kqTypeLabel[kElectron2nd]="bce";
 284   } else if (kquark == kOthers){
 285     kqTypeLabel[kGamma-4]="gammae";
 286     kqTypeLabel[kPi0-4]="pi0e";
 287     kqTypeLabel[kElse-4]="elsee";
 288     kqTypeLabel[kMisID-4]="miside";
 289   }
 290
 291   TString kqEtaRangeLabel[fgkEtaRanges];
 292   kqEtaRangeLabel[0] = "mcqa_";
 293   kqEtaRangeLabel[1] = "mcqa_barrel_";
 294   kqEtaRangeLabel[2] = "mcqa_unitY_";
 295
 296   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 297   const Int_t nptbinning1 = 35;
 298   Int_t iBin[2];
 299   iBin[0] = 44; // bins in pt
 300   iBin[1] = nptbinning1; // bins in pt
 301   Double_t* binEdges[1];
 302   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 303
 304   // new binning for final electron analysis
 305   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.};
 306
 307   const Int_t ndptbins = 500;
 308   Double_t xcorrbin[ndptbins+1];
 309   for (int icorrbin = 0; icorrbin< ndptbins+1; icorrbin++){
 310     xcorrbin[icorrbin]=icorrbin*0.1;
 311   }
 312
 313   Int_t fCentrmax = 0;
 314   if(!fIsPbPb&&!fIsppMultiBin) fCentrmax=0;
 315   else fCentrmax = kCentBins;
 316   TString hname;
 317   if(kquark == kOthers){
 318    for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
 319        for (Int_t iqType = 0; iqType < 4; iqType++ ){
 320            for(Int_t icentr = 0; icentr < (fCentrmax+1); icentr++)
 321            {
 322                hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
 323                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);
 324                hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
 325                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);
 326                hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
 327                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);
 328                // Fill List
 329                if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
 330            }
 331        }
 332    }
 333    return;
 334   }
 335   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
 336    for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
 337        if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
 338         for(Int_t icentr = 0; icentr<(fCentrmax+1); icentr++)
 339            {
 340                hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
 341                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);
 342                hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
 343                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
 344                hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
 345                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);
 346                hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
 347                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);
 348                // Fill List
 349                if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
 350            }
 351    }
 352   }
 353
 354   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
 355     hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
 356     fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 357     hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
 358     fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 359     hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
 360     fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 361     hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
 362     fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
 363
 364     hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
 365     fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
 366     hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
 367     fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
 368     hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
 369     fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
 370     hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
 371     fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
 372     if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
 373   }
 374
 375   hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
 376   fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
 377   hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
 378   fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
 379
 380 }
 381
 382 //__________________________________________
 383 void AliHFEmcQA::Init()
 384 {
 385   // called at begining every event
 386
 387   for (Int_t i=0; i<2; i++){
 388      fIsHeavy[i] = 0;
 389   }
 390
 391   fNparents = 7;
 392
 393   fParentSelect[0][0] =  411; //D+
 394   fParentSelect[0][1] =  421; //D0
 395   fParentSelect[0][2] =  431; //Ds+
 396   fParentSelect[0][3] = 4122; //Lambdac+
 397   fParentSelect[0][4] = 4132; //Ksic0
 398   fParentSelect[0][5] = 4232; //Ksic+
 399   fParentSelect[0][6] = 4332; //OmegaC0
 400
 401   fParentSelect[1][0] =  511; //B0
 402   fParentSelect[1][1] =  521; //B+
 403   fParentSelect[1][2] =  531; //Bs0
 404   fParentSelect[1][3] = 5122; //Lambdab0
 405   fParentSelect[1][4] = 5132; //Ksib-
 406   fParentSelect[1][5] = 5232; //Ksib0
 407   fParentSelect[1][6] = 5332; //Omegab-
 408
 409
 410 }
 411
 412 //__________________________________________
 413 void AliHFEmcQA::GetMesonKine()
 414 {
 415   //
 416   // get meson pt spectra
 417   //
 418
 419   AliVParticle *mctrack2 = NULL;
 420   AliMCParticle *mctrack0 = NULL;
 421   AliVParticle *mctrackdaugt= NULL;
 422   AliMCParticle *mctrackd= NULL;
 423   Int_t id1=0, id2=0;
 424
 425
 426   if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
 427   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 428
 429   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
 430      if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
 431      TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
 432      if(!mcpart0) continue;
 433      mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
 434      if(!mctrack0) continue;
 435
 436 //     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 437
 438      if(abs(mctrack0->PdgCode()) == 111) // pi0
 439        {
 440           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 441             fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
 442             fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 443           }
 444           id1=mctrack0->GetFirstDaughter();
 445           id2=mctrack0->GetLastDaughter();
 446           if(!((id2-id1)==2)) continue;
 447           for(int idx=id1; idx<=id2; idx++){
 448             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 449             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 450             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 451              fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
 452           }
 453        }
 454      else if(abs(mctrack0->PdgCode()) == 221) // eta
 455        {
 456           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 457             fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
 458             fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 459           }
 460           id1=mctrack0->GetFirstDaughter();
 461           id2=mctrack0->GetLastDaughter();
 462           if(!((id2-id1)==2||(id2-id1)==3)) continue;
 463           for(int idx=id1; idx<=id2; idx++){
 464             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 465             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 466             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 467              fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
 468           }
 469        }
 470      else if(abs(mctrack0->PdgCode()) == 223) // omega
 471        {
 472           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 473             fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
 474             fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 475           }
 476           id1=mctrack0->GetFirstDaughter();
 477           id2=mctrack0->GetLastDaughter();
 478           if(!((id2-id1)==1||(id2-id1)==2)) continue;
 479           for(int idx=id1; idx<=id2; idx++){
 480             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 481             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 482             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 483              fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
 484           }
 485        }
 486      else if(abs(mctrack0->PdgCode()) == 333) // phi
 487        {
 488           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 489             fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
 490             fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 491           }
 492           id1=mctrack0->GetFirstDaughter();
 493           id2=mctrack0->GetLastDaughter();
 494           if(!((id2-id1)==1)) continue;
 495           for(int idx=id1; idx<=id2; idx++){
 496             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 497             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 498             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 499              fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
 500           }
 501        }
 502      else if(abs(mctrack0->PdgCode()) == 331) // eta prime
 503        {
 504           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 505             fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
 506             fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 507           }
 508           id1=mctrack0->GetFirstDaughter();
 509           id2=mctrack0->GetLastDaughter();
 510           if(!((id2-id1)==2||(id2-id1)==3)) continue;
 511           for(int idx=id1; idx<=id2; idx++){
 512             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 513             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 514             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 515              fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
 516           }
 517        }
 518      else if(abs(mctrack0->PdgCode()) == 113) // rho
 519        {
 520           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 521             fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
 522             fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 523           }
 524           id1=mctrack0->GetFirstDaughter();
 525           id2=mctrack0->GetLastDaughter();
 526           if(!((id2-id1)==1)) continue;
 527           for(int idx=id1; idx<=id2; idx++){
 528             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 529             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 530             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 531              fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
 532           }
 533        }
 534      else if(abs(mctrack0->PdgCode()) == 321) // kaon+-
 535        {
 536           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 537             fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 538           }
 539        }
 540      else if(abs(mctrack0->PdgCode()) == 130) // k0L
 541        {
 542           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 543             fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 544           }
 545        }
 546      }
 547
 548 }
 549 //__________________________________________
 550 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
 551 {
 552   // get heavy quark kinematics
 553
 554     if (kquark != kCharm && kquark != kBeauty) {
 555       AliDebug(1, "This task is only for heavy quark QA, return\n");
 556       return;
 557     }
 558     Int_t iq = kquark - kCharm;
 559
 560     if (iTrack < 0 || !part) {
 561       AliDebug(1, "Stack label is negative or no mcparticle, return\n");
 562       return;
 563     }
 564
 565     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 566
 567     AliMCParticle *mctrack = NULL;
 568     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 569
 570     // select heavy hadron or not fragmented heavy quark
 571     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
 572
 573       TParticle *partMother;
 574       Int_t iLabel;
 575
 576       if (partPdgcode == kquark){ // in case of not fragmented heavy quark
 577         partMother = part;
 578         iLabel = iTrack;
 579       } else{ // in case of heavy hadron, start to search for mother heavy parton
 580         iLabel = part->GetFirstMother();
 581         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 582         if (iLabel>-1) { partMother = mctrack->Particle(); }
 583         else {
 584           AliDebug(1, "Stack label is negative, return\n");
 585           return;
 586         }
 587       }
 588
 589       // heavy parton selection as a mother of heavy hadron
 590       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
 591       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
 592       // should I make a condition that partMother should be quark or diquark? -> not necessary
 593       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
 594       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
 595
 596         if ( abs(partMother->GetPdgCode()) != kquark ){
 597           // search fragmented partons in the same string
 598           Bool_t isSameString = kTRUE;
 599           for (Int_t i=1; i<fgkMaxIter; i++){
 600              iLabel = iLabel - 1;
 601              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 602              if (iLabel>-1) { partMother = mctrack->Particle(); }
 603              else {
 604                AliDebug(1, "Stack label is negative, return\n");
 605                return;
 606              }
 607              if ( abs(partMother->GetPdgCode()) == kquark ) break;
 608              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
 609              if (!isSameString) return;
 610           }
 611         }
 612         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
 613         if (abs(partMother->GetPdgCode()) != kquark) return;
 614
 615         if (fIsHeavy[iq] >= 50) return;
 616         fHeavyQuark[fIsHeavy[iq]] = partMother;
 617         fIsHeavy[iq]++;
 618
 619         // fill kinematics for heavy parton
 620         if (partMother->GetPdgCode() > 0) { // quark
 621           fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
 622           fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
 623           fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
 624           fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
 625         } else{ // antiquark
 626           fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
 627           fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
 628           fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
 629           fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
 630         }
 631
 632       } // end of heavy parton slection loop
 633
 634     } // end of heavy hadron or quark selection
 635
 636 }
 637
 638 //__________________________________________
 639 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
 640 {
 641   // end of event analysis
 642
 643   if (kquark != kCharm && kquark != kBeauty) {
 644     AliDebug(1, "This task is only for heavy quark QA, return\n");
 645     return;
 646   }
 647   Int_t iq = kquark - kCharm;
 648
 649
 650   // # of heavy quark per event
 651   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
 652   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
 653
 654   Int_t motherID[fgkMaxGener];
 655   Int_t motherType[fgkMaxGener];
 656   Int_t motherLabel[fgkMaxGener];
 657   Int_t ancestorPdg[fgkMaxGener];
 658   Int_t ancestorLabel[fgkMaxGener];
 659
 660   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
 661      motherID[i] = 0;
 662      motherType[i] = 0;
 663      motherLabel[i] = 0;
 664      ancestorPdg[i] = 0;
 665      ancestorLabel[i] = 0;
 666   }
 667
 668
 669   // check history of found heavy quarks
 670   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 671
 672      if(!fHeavyQuark[i]) return;
 673
 674      ancestorLabel[0] = i;
 675      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
 676      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
 677
 678      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
 679      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
 680
 681      Int_t ig = 1;
 682      while (ancestorLabel[ig] != -1){
 683           // in case there is mother, get mother's pdg code and grandmother's label
 684           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
 685           // if mother is still heavy, find again mother's ancestor
 686           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
 687             ig++;
 688             continue; // if it is from same heavy
 689           }
 690           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
 691           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 692           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 693           // if it is not the above case, something is strange
 694           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
 695           break;
 696      }
 697      if (ancestorLabel[ig] == -1){ // from hard scattering
 698        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
 699      }
 700
 701   } // end of found heavy quark loop
 702
 703
 704   // check process type
 705   Int_t processID = 0;
 706   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 707      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
 708   }
 709
 710
 711   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
 712   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
 713
 714      Int_t id1 = ipair*2;
 715      Int_t id2 = ipair*2 + 1;
 716
 717      if (motherType[id1] == 2 && motherType[id2] == 2){
 718        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
 719        else processID = -9;
 720      }
 721      else if (motherType[id1] == -1 && motherType[id2] == -1) {
 722        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
 723          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
 724          else processID = kPairCreationFromq; // q-qbar pair creation
 725        }
 726        else processID = -8;
 727      }
 728      else if (motherType[id1] == -1 || motherType[id2] == -1) {
 729        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
 730          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
 731          else processID = kLightQuarkShower;
 732        }
 733        else processID = -7;
 734      }
 735      else if (motherType[id1] == -2 || motherType[id2] == -2) {
 736        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
 737        else processID = -6;
 738
 739      }
 740      else processID = -5;
 741
 742      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
 743      else fHistComm[iq][0].fProcessID->Fill(processID);
 744      AliDebug(1,Form("Process ID = %d\n",processID));
 745   } // end of # heavy quark pair loop
 746
 747 }
 748
 749 //__________________________________________
 750 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
 751 {
 752     // decay electron kinematics
 753
 754     if (kquark != kCharm && kquark != kBeauty) {
 755       AliDebug(1, "This task is only for heavy quark QA, return\n");
 756       return;
 757     }
 758     Int_t iq = kquark - kCharm;
 759
 760     if(!mcpart){
 761       AliDebug(1, "no mc particle, return\n");
 762       return;
 763     }
 764
 765     Int_t iLabel = mcpart->GetFirstMother();
 766     if (iLabel<0){
 767       AliDebug(1, "Stack label is negative, return\n");
 768       return;
 769     }
 770
 771     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 772
 773     TParticle *partCopy = mcpart;
 774     Int_t pdgcode = mcpart->GetPdgCode();
 775     Int_t pdgcodeCopy = pdgcode;
 776
 777     AliMCParticle *mctrack = NULL;
 778
 779     // if the mother is charmed hadron
 780     Bool_t isDirectCharm = kFALSE;
 781     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
 782
 783           // iterate until you find B hadron as a mother or become top ancester
 784           for (Int_t i=1; i<fgkMaxIter; i++){
 785
 786              Int_t jLabel = mcpart->GetFirstMother();
 787              if (jLabel == -1){
 788                isDirectCharm = kTRUE;
 789                break; // if there is no ancester
 790              }
 791              if (jLabel < 0){ // safety protection
 792                AliDebug(1, "Stack label is negative, return\n");
 793                return;
 794              }
 795              // if there is an ancester
 796              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
 797              TParticle* mother = mctrack->Particle();
 798              Int_t motherPDG = mother->GetPdgCode();
 799
 800              for (Int_t j=0; j<fNparents; j++){
 801                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
 802              }
 803
 804              mcpart = mother;
 805           } // end of iteration
 806     } // end of if
 807     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 808          for (Int_t i=0; i<fNparents; i++){
 809             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
 810
 811               // fill hadron kinematics
 812               fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
 813               fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
 814               fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
 815               fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
 816
 817               if(iq==0) {
 818                fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 819               }
 820             }
 821          }
 822          // I also want to store D* info to compare with D* measurement
 823          if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
 824                fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 825          }
 826          if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
 827                fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 828          }
 829     } // end of if
 830 }
 831
 832 //__________________________________________
 833 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed)
 834 {
 835     // decay electron kinematics
 836
 837     if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
 838       AliDebug(1, "This task is only for heavy quark QA, return\n");
 839       return;
 840     }
 841     Int_t iq = kquark - kCharm;
 842     Bool_t isFinalOpenCharm = kFALSE;
 843
 844     if(!mcpart){
 845       AliDebug(1, "no mcparticle, return\n");
 846       return;
 847     }
 848
 849     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 850
 851     Double_t eabsEta = TMath::Abs(mcpart->Eta());
 852     Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
 853
 854     if(kquark==kOthers){
 855       Int_t esource = -1;
 856       if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
 857       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
 858       if(esource==0|| esource==1 || esource==2 || esource==3){
 859         fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
 860         fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 861         fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
 862         if(eabsEta<0.9){
 863           fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
 864           fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 865           fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
 866         }
 867         if(eabsY<0.5){
 868           fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
 869           fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 870           fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
 871         }
 872         return;
 873       }
 874       else {
 875         AliDebug(1, "e source is out of defined ranges, return\n");
 876         return;
 877       }
 878     }
 879
 880     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 881
 882     Int_t iLabel = mcpart->GetFirstMother();
 883     if (iLabel<0){
 884       AliDebug(1, "Stack label is negative, return\n");
 885       return;
 886     }
 887
 888     AliMCParticle *mctrack = NULL;
 889     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 890     TParticle *partMother = mctrack->Particle();
 891     TParticle *partMotherCopy = partMother;
 892     Int_t maPdgcode = partMother->GetPdgCode();
 893     Int_t maPdgcodeCopy = maPdgcode;
 894
 895     // get mc primary vertex
 896     /*
 897     TArrayF mcPrimVtx;
 898     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
 899
 900     // get electron production vertex
 901     TLorentzVector ePoint;
 902     mcpart->ProductionVertex(ePoint);
 903
 904     // calculated production vertex to primary vertex (in xy plane)
 905     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
 906     */
 907     Float_t decayLxy = 0;
 908
 909     // if the mother is charmed hadron
 910     Bool_t isMotherDirectCharm = kFALSE;
 911     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 912
 913          for (Int_t i=0; i<fNparents; i++){
 914             if (abs(maPdgcode)==fParentSelect[0][i]){
 915               isFinalOpenCharm = kTRUE;
 916             }
 917          }
 918          if (!isFinalOpenCharm) return ;
 919
 920           // iterate until you find B hadron as a mother or become top ancester
 921           for (Int_t i=1; i<fgkMaxIter; i++){
 922
 923              Int_t jLabel = partMother->GetFirstMother();
 924              if (jLabel == -1){
 925                isMotherDirectCharm = kTRUE;
 926                break; // if there is no ancester
 927              }
 928              if (jLabel < 0){ // safety protection
 929                AliDebug(1, "Stack label is negative, return\n");
 930                return;
 931              }
 932
 933              // if there is an ancester
 934              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
 935              TParticle* grandMa = mctrack->Particle();
 936              Int_t grandMaPDG = grandMa->GetPdgCode();
 937
 938              for (Int_t j=0; j<fNparents; j++){
 939                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 940
 941                   if (kquark == kCharm) return;
 942                   // fill electron kinematics
 943                   fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
 944                   fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
 945                   fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 946                   fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
 947
 948                   // fill mother hadron kinematics
 949                   fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
 950                   fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
 951                   fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
 952                   fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
 953
 954                   if(eabsEta<0.9){
 955                     fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
 956                     fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
 957                     fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 958                     fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
 959
 960                     // fill mother hadron kinematics
 961                     fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
 962                     fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
 963                     fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
 964                     fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
 965                   }
 966
 967                   if(eabsY<0.5){
 968                     fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
 969                     fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
 970                     fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 971                     fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
 972
 973                     // fill mother hadron kinematics
 974                     fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
 975                     fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
 976                     fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
 977                     fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
 978                   }
 979
 980                   // ratio between pT of electron and pT of mother B hadron
 981                   if(grandMa->Pt()) {
 982                     fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 983                     if(eabsEta<0.9){
 984                       fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 985                     }
 986                     if(eabsY<0.5){
 987                       fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 988                     }
 989                   }
 990
 991                   // distance between electron production point and primary vertex
 992                   fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
 993                   if(eabsEta<0.9){
 994                     fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
 995                   }
 996                   if(eabsY<0.5){
 997                     fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
 998                   }
 999                   return;
1000                 }
1001              }
1002
1003              partMother = grandMa;
1004           } // end of iteration
1005     } // end of if
1006     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1007          for (Int_t i=0; i<fNparents; i++){
1008             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1009
1010               fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1011               fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1012               fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1013               fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1014
1015               // fill mother hadron kinematics
1016               fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1017               fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1018               fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1019               fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1020
1021               if(eabsEta<0.9){
1022                 fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1023                 fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1024                 fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1025                 fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1026
1027                 // fill mother hadron kinematics
1028                 fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1029                 fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1030                 fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1031                 fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1032               }
1033
1034               if(eabsY<0.5){
1035                 fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1036                 fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1037                 fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1038                 fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1039
1040                 // fill mother hadron kinematics
1041                 fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1042                 fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1043                 fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1044                 fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1045               }
1046
1047               // ratio between pT of electron and pT of mother B or direct D hadron
1048               if(partMotherCopy->Pt()) {
1049                  fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1050                  if(eabsEta<0.9){
1051                    fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1052                  }
1053                  if(eabsY<0.5){
1054                    fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
1055                  }
1056               }
1057               fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1058               if(eabsEta<0.9){
1059                 fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1060               }
1061               if(eabsY<0.5){
1062                 fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
1063               }
1064               if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
1065                 fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1066                 if(eabsEta<0.9){
1067                   fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1068                 }
1069                 if(eabsY<0.5){
1070                   fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
1071                 }
1072               }
1073               else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
1074                 fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1075                 if(eabsEta<0.9){
1076                   fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1077                 }
1078                 if(eabsY<0.5){
1079                   fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
1080                 }
1081               }
1082               else {
1083                 fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1084                 if(eabsEta<0.9){
1085                   fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1086                 }
1087                 if(eabsY<0.5){
1088                   fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
1089                 }
1090               }
1091
1092               // distance between electron production point and primary vertex
1093               fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1094               if(eabsEta<0.9){
1095                 fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1096               }
1097               if(eabsY<0.5){
1098                 fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
1099               }
1100             }
1101          }
1102     } // end of if
1103 }
1104
1105 //____________________________________________________________________
1106 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
1107 {
1108   // decay electron kinematics
1109
1110   if (kquark != kCharm && kquark != kBeauty) {
1111     AliDebug(1, "This task is only for heavy quark QA, return\n");
1112     return;
1113   }
1114
1115   Int_t iq = kquark - kCharm;
1116   Bool_t isFinalOpenCharm = kFALSE;
1117
1118   if(!mcpart){
1119     AliDebug(1, "no mcparticle, return\n");
1120     return;
1121   }
1122
1123   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
1124
1125   Double_t eabsEta = TMath::Abs(mcpart->Eta());
1126   Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
1127
1128   // mother
1129   Int_t iLabel = mcpart->GetMother();
1130   if (iLabel<0){
1131     AliDebug(1, "Stack label is negative, return\n");
1132     return;
1133   }
1134
1135   if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
1136
1137   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1138   AliAODMCParticle *partMotherCopy = partMother;
1139   Int_t maPdgcode = partMother->GetPdgCode();
1140   Int_t maPdgcodeCopy = maPdgcode;
1141
1142   Bool_t isMotherDirectCharm = kFALSE;
1143   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1144
1145     for (Int_t i=0; i<fNparents; i++){
1146        if (abs(maPdgcode)==fParentSelect[0][i]){
1147          isFinalOpenCharm = kTRUE;
1148        }
1149     }
1150     if (!isFinalOpenCharm) return;
1151
1152     for (Int_t i=1; i<fgkMaxIter; i++){
1153
1154        Int_t jLabel = partMother->GetMother();
1155        if (jLabel == -1){
1156          isMotherDirectCharm = kTRUE;
1157          break; // if there is no ancester
1158        }
1159        if (jLabel < 0){ // safety protection
1160          AliDebug(1, "Stack label is negative, return\n");
1161          return;
1162        }
1163
1164        // if there is an ancester
1165        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1166        Int_t grandMaPDG = grandMa->GetPdgCode();
1167
1168        for (Int_t j=0; j<fNparents; j++){
1169           if (abs(grandMaPDG)==fParentSelect[1][j]){
1170
1171             if (kquark == kCharm) return;
1172             // fill electron kinematics
1173             fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1174             fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
1175             fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1176             fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
1177
1178             // fill mother hadron kinematics
1179             fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
1180             fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
1181             fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1182             fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
1183
1184             if(eabsEta<0.9){
1185               // fill electron kinematics
1186               fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1187               fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
1188               fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1189               fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
1190
1191               // fill mother hadron kinematics
1192               fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
1193               fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
1194               fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1195               fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
1196             }
1197             if(eabsY<0.5){
1198               // fill electron kinematics
1199               fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1200               fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
1201               fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1202               fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
1203
1204               // fill mother hadron kinematics
1205               fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
1206               fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
1207               fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1208               fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
1209             }
1210
1211             return;
1212           }
1213        }
1214
1215        partMother = grandMa;
1216     } // end of iteration
1217   } // end of if
1218   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1219     for (Int_t i=0; i<fNparents; i++){
1220        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1221
1222          // fill electron kinematics
1223          fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1224          fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
1225          fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1226          fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
1227
1228          // fill mother hadron kinematics
1229          fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1230          fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
1231          fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1232          fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
1233
1234          if(eabsEta<0.9){
1235            // fill electron kinematics
1236            fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1237            fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
1238            fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1239            fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
1240
1241            // fill mother hadron kinematics
1242            fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1243            fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
1244            fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1245            fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
1246          }
1247          if(eabsY<0.5){
1248            // fill electron kinematics
1249            fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
1250            fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
1251            fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1252            fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
1253
1254            // fill mother hadron kinematics
1255            fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
1256            fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
1257            fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1258            fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
1259          }
1260
1261        }
1262     }
1263   } // end of if
1264
1265 }
1266
1267 //__________________________________________
1268 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
1269 {
1270        // find mother pdg code and label
1271
1272        if (motherlabel < 0) {
1273          AliDebug(1, "Stack label is negative, return\n");
1274          return;
1275        }
1276        AliMCParticle *mctrack = NULL;
1277        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1278        TParticle *heavysMother = mctrack->Particle();
1279        motherpdg = heavysMother->GetPdgCode();
1280        grandmotherlabel = heavysMother->GetFirstMother();
1281        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
1282 }
1283
1284 //__________________________________________
1285 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1286 {
1287        // mothertype -1 means this heavy quark coming from hard vertex
1288
1289        AliMCParticle *mctrack1 = NULL;
1290        AliMCParticle *mctrack2 = NULL;
1291        if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1292        if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1293        TParticle *afterinitialrad1  = mctrack1->Particle();
1294        TParticle *afterinitialrad2  = mctrack2->Particle();
1295
1296        motherlabel = -1;
1297
1298        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1299          AliDebug(1,"heavy from gluon gluon pair creation!\n");
1300          mothertype = -1;
1301          motherID = fgkGluon;
1302        }
1303        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1304          AliDebug(1,"heavy from flavor exitation!\n");
1305          mothertype = -1;
1306          motherID = kquark;
1307        }
1308        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
1309          AliDebug(1,"heavy from q-qbar pair creation!\n");
1310          mothertype = -1;
1311          motherID = 1;
1312        }
1313        else {
1314          AliDebug(1,"something strange!\n");
1315          mothertype = -999;
1316          motherlabel = -999;
1317          motherID = -999;
1318        }
1319 }
1320
1321 //__________________________________________
1322 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1323 {
1324        // mothertype -2 means this heavy quark coming from initial state
1325
1326        AliMCParticle *mctrack = NULL;
1327        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
1328          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1329          TParticle *heavysMother = mctrack->Particle();
1330          motherID = heavysMother->GetPdgCode();
1331          mothertype = -2; // there is mother before initial state radiation
1332          motherlabel = inputmotherlabel;
1333          AliDebug(1,"initial parton shower! \n");
1334
1335          return kTRUE;
1336        }
1337
1338        return kFALSE;
1339 }
1340
1341 //__________________________________________
1342 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1343 {
1344        // mothertype 2 means this heavy quark coming from final state
1345
1346        AliMCParticle *mctrack = NULL;
1347        if (inputmotherlabel > 5){ // mother exist after hard scattering
1348          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1349          TParticle *heavysMother = mctrack->Particle();
1350          motherID = heavysMother->GetPdgCode();
1351          mothertype = 2; //
1352          motherlabel = inputmotherlabel;
1353          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1354
1355          return kTRUE;
1356        }
1357        return kFALSE;
1358 }
1359
1360 //__________________________________________
1361 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1362 {
1363       // mark strange behavior
1364
1365        mothertype = -888;
1366        motherlabel = -888;
1367        motherID = -888;
1368        AliDebug(1,"something strange!\n");
1369 }
1370
1371 //__________________________________________
1372 Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart)
1373 {
1374   // decay particle's origin
1375
1376   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1377
1378   Int_t origin = -1;
1379   Bool_t isFinalOpenCharm = kFALSE;
1380
1381   if(!mcpart){
1382     AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1383     return -1;
1384   }
1385
1386   // mother
1387   // Information not in the base class, cast necessary
1388   Int_t iLabel = GetMother(mcpart);
1389   if (iLabel<0){
1390     AliDebug(1, "Stack label is negative, return\n");
1391     return -1;
1392   }
1393
1394   const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
1395   Int_t maPdgcode = partMother->PdgCode();
1396
1397   // if the mother is charmed hadron
1398   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1399
1400     for (Int_t i=0; i<fNparents; i++){
1401        if (abs(maPdgcode)==fParentSelect[0][i]){
1402          isFinalOpenCharm = kTRUE;
1403        }
1404     }
1405     if (!isFinalOpenCharm) return -1;
1406
1407     // iterate until you find B hadron as a mother or become top ancester
1408     for (Int_t i=1; i<fgkMaxIter; i++){
1409
1410        Int_t jLabel = GetMother(partMother);
1411        if (jLabel == -1){
1412          origin = kDirectCharm;
1413          return origin;
1414        }
1415        if (jLabel < 0){ // safety protection
1416          AliDebug(1, "Stack label is negative, return\n");
1417          return -1;
1418        }
1419
1420        // if there is an ancester
1421        const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
1422        Int_t grandMaPDG = grandMa->PdgCode();
1423
1424        for (Int_t j=0; j<fNparents; j++){
1425           if (abs(grandMaPDG)==fParentSelect[1][j]){
1426             origin = kBeautyCharm;
1427             return origin;
1428           }
1429        }
1430
1431        partMother = grandMa;
1432     } // end of iteration
1433   } // end of if
1434   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1435     for (Int_t i=0; i<fNparents; i++){
1436        if (abs(maPdgcode)==fParentSelect[1][i]){
1437          origin = kDirectBeauty;
1438          return origin;
1439        }
1440     }
1441   } // end of if
1442   else if ( abs(maPdgcode) == 22 ) {
1443     origin = kGamma;
1444     return origin;
1445   } // end of if
1446   else if ( abs(maPdgcode) == 111 ) {
1447     origin = kPi0;
1448     return origin;
1449   } // end of if
1450
1451   return origin;
1452 }
1453
1454 //__________________________________________
1455 Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
1456 {
1457   // decay particle's origin
1458
1459   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1460
1461   Int_t origin = -1;
1462   Bool_t isFinalOpenCharm = kFALSE;
1463
1464   if(!mcpart){
1465     AliDebug(1, "no mcparticle, return\n");
1466     return -1;
1467   }
1468
1469   Int_t iLabel = mcpart->GetFirstMother();
1470   if (iLabel<0){
1471     AliDebug(1, "Stack label is negative, return\n");
1472     return -1;
1473   }
1474
1475   AliMCParticle *mctrack = NULL;
1476   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1477   TParticle *partMother = mctrack->Particle();
1478   Int_t maPdgcode = partMother->GetPdgCode();
1479
1480    // if the mother is charmed hadron
1481    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1482
1483      for (Int_t i=0; i<fNparents; i++){
1484         if (abs(maPdgcode)==fParentSelect[0][i]){
1485           isFinalOpenCharm = kTRUE;
1486         }
1487      }
1488      if (!isFinalOpenCharm) return -1;
1489
1490      // iterate until you find B hadron as a mother or become top ancester
1491      for (Int_t i=1; i<fgkMaxIter; i++){
1492
1493         Int_t jLabel = partMother->GetFirstMother();
1494         if (jLabel == -1){
1495           origin = kDirectCharm;
1496           return origin;
1497         }
1498         if (jLabel < 0){ // safety protection
1499           AliDebug(1, "Stack label is negative, return\n");
1500           return -1;
1501         }
1502
1503         // if there is an ancester
1504         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1505         TParticle* grandMa = mctrack->Particle();
1506         Int_t grandMaPDG = grandMa->GetPdgCode();
1507
1508         for (Int_t j=0; j<fNparents; j++){
1509            if (abs(grandMaPDG)==fParentSelect[1][j]){
1510              origin = kBeautyCharm;
1511              return origin;
1512            }
1513         }
1514
1515         partMother = grandMa;
1516      } // end of iteration
1517    } // end of if
1518    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1519      for (Int_t i=0; i<fNparents; i++){
1520         if (abs(maPdgcode)==fParentSelect[1][i]){
1521           origin = kDirectBeauty;
1522           return origin;
1523         }
1524      }
1525    } // end of if
1526    else if ( abs(maPdgcode) == 22 ) {
1527      origin = kGamma;
1528      return origin;
1529    } // end of if
1530    else if ( abs(maPdgcode) == 111 ) {
1531      origin = kPi0;
1532      return origin;
1533    } // end of if
1534    else origin = kElse;
1535
1536    return origin;
1537 }
1538
1539 //__________________________________________
1540 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1541 {
1542   // decay particle's origin
1543
1544   if(!mcpart){
1545     AliDebug(1, "no mcparticle, return\n");
1546     return -1;
1547   }
1548
1549   if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1550
1551   Int_t origin = -1;
1552   Bool_t isFinalOpenCharm = kFALSE;
1553
1554   Int_t iLabel = mcpart->GetFirstMother();
1555   if (iLabel<0){
1556     AliDebug(1, "Stack label is negative, return\n");
1557     return -1;
1558   }
1559
1560   AliMCParticle *mctrack = NULL;
1561   Int_t tmpMomLabel=0;
1562   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1563   TParticle *partMother = mctrack->Particle();
1564   TParticle *partMotherCopy = mctrack->Particle();
1565   Int_t maPdgcode = partMother->GetPdgCode();
1566
1567    // if the mother is charmed hadron
1568
1569    if(abs(maPdgcode)==443){ // J/spi
1570       Int_t jLabel = partMother->GetFirstMother();
1571       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return kJpsi;
1572       TParticle* grandMa = mctrack->Particle();
1573       Int_t grandMaPDG = grandMa->GetPdgCode();
1574       if((int(abs(grandMaPDG)/100.)%10) == kBeauty || (int(abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
1575       else return -1;
1576    }
1577    else if ( (int(abs(maPdgcode)/100.)%10) == kCharm || (int(abs(maPdgcode)/1000.)%10) == kCharm ) {
1578
1579      for (Int_t i=0; i<fNparents; i++){
1580         if (abs(maPdgcode)==fParentSelect[0][i]){
1581           isFinalOpenCharm = kTRUE;
1582         }
1583      }
1584      if (!isFinalOpenCharm) return -1;
1585
1586      // iterate until you find B hadron as a mother or become top ancester
1587      for (Int_t i=1; i<fgkMaxIter; i++){
1588
1589         Int_t jLabel = partMother->GetFirstMother();
1590         if (jLabel == -1){
1591           return kDirectCharm;
1592         }
1593         if (jLabel < 0){ // safety protection
1594           AliDebug(1, "Stack label is negative, return\n");
1595           return -1;
1596         }
1597
1598         // if there is an ancester
1599         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1600         TParticle* grandMa = mctrack->Particle();
1601         Int_t grandMaPDG = grandMa->GetPdgCode();
1602
1603         for (Int_t j=0; j<fNparents; j++){
1604            if (abs(grandMaPDG)==fParentSelect[1][j]){
1605              return kBeautyCharm;
1606            }
1607         }
1608
1609         partMother = grandMa;
1610      } // end of iteration
1611    } // end of if
1612    else if ( (int(abs(maPdgcode)/100.)%10) == kBeauty || (int(abs(maPdgcode)/1000.)%10) == kBeauty ) {
1613      for (Int_t i=0; i<fNparents; i++){
1614         if (abs(maPdgcode)==fParentSelect[1][i]){
1615           return kDirectBeauty;
1616         }
1617      }
1618    } // end of if
1619    else if ( abs(maPdgcode) == 22 ) { //conversion
1620
1621      tmpMomLabel = partMotherCopy->GetFirstMother();
1622      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1623      partMother = mctrack->Particle();
1624      maPdgcode = partMother->GetPdgCode();
1625      if ( abs(maPdgcode) == 111 ) {
1626        return kGammaPi0;
1627      }
1628      else if ( abs(maPdgcode) == 221 ) {
1629        return kGammaEta;
1630      }
1631      else if ( abs(maPdgcode) == 223 ) {
1632        return kGammaOmega;
1633      }
1634      else if ( abs(maPdgcode) == 333 ) {
1635        return kGammaPhi;
1636      }
1637      else if ( abs(maPdgcode) == 331 ) {
1638        return kGammaEtaPrime;
1639      }
1640      else if ( abs(maPdgcode) == 113 ) {
1641        return kGammaRho0;
1642      }
1643      else origin = kElse;
1644      return origin;
1645
1646    }
1647    else if ( abs(maPdgcode) == 111 ) {
1648      return kPi0;
1649    }
1650    else if ( abs(maPdgcode) == 221 ) {
1651      return kEta;
1652    }
1653    else if ( abs(maPdgcode) == 223 ) {
1654      return kOmega;
1655    }
1656    else if ( abs(maPdgcode) == 333 ) {
1657      return kPhi;
1658    }
1659    else if ( abs(maPdgcode) == 331 ) {
1660      return kEtaPrime;
1661    }
1662    else if ( abs(maPdgcode) == 113 ) {
1663      return kRho0;
1664    }
1665    else if ( abs(maPdgcode) == 321 || abs(maPdgcode) == 130 ) {
1666      return kKe3;
1667    }
1668    else{
1669     origin = kElse;
1670    }
1671    return origin;
1672 }
1673 //__________________________________________
1674 Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
1675   //
1676   // 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)
1677   //
1678   AliMCParticle *mctrackmother = NULL;
1679   Double_t weightElecBg = 0.;
1680   Double_t mesonPt = 0.;
1681   Double_t bgcategory = 0.;
1682   Int_t mArr = -1;
1683   Int_t mesonID = GetElecSource(mctrack->Particle());
1684   if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0;                //pion
1685   else if(mesonID==kGammaEta || mesonID==kEta) mArr=1;           //eta
1686   else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2;       //omega
1687   else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3;           //phi
1688   else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
1689   else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
1690
1691   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};
1692   Double_t xr[3]={-999,-999,-999};
1693   datamc[0] = mesonID;
1694   datamc[17] = mctrack->Pt(); //electron pt
1695   datamc[18] = mctrack->Eta(); //electron eta
1696
1697   mctrack->XvYvZv(xr);
1698   datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1699   datamc[10] = xr[2];
1700
1701   TParticle *mcpart = mctrack->Particle();
1702   if(mcpart){
1703     datamc[14] = mcpart->GetUniqueID();
1704   }
1705
1706   if(!(mArr<0)){
1707      if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering
1708         Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
1709         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1710           glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
1711           if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1712             mesonPt = mctrackmother->Pt(); //meson pt
1713             bgcategory = 1.;
1714             datamc[1] = bgcategory;
1715             datamc[2] = mesonPt;
1716             mctrackmother->XvYvZv(xr);
1717             datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1718             datamc[12] = xr[2];
1719
1720             mcpart = mctrackmother->Particle();
1721             if(mcpart){
1722               datamc[15] = mcpart->GetUniqueID();
1723             }
1724             if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1725               bgcategory = 2.;
1726               datamc[1] = bgcategory;
1727               //printf("I should be gamma meson = %d  mesonlabel= %d  NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries());
1728               glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1729               if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1730                 datamc[3]=mctrackmother->PdgCode();
1731                 datamc[4]=mctrackmother->Pt();
1732                 if(TMath::Abs(mctrackmother->PdgCode())==310){
1733                   bgcategory = 3.;
1734                   glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1735                   if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1736                     datamc[5]=mctrackmother->PdgCode();
1737                     glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1738                     if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1739                        datamc[6]=mctrackmother->PdgCode();
1740                     }
1741                   }
1742                 }
1743               }
1744             }
1745           }
1746         }
1747      }
1748      else{ // nonHFE except for the conversion electron
1749         Int_t glabel=TMath::Abs(mctrack->GetMother());
1750         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1751           mesonPt=mctrackmother->Pt(); //meson pt
1752           bgcategory = -1.;
1753           datamc[1] = bgcategory;
1754           datamc[2] = mesonPt;
1755           mctrackmother->XvYvZv(xr);
1756           datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1757           datamc[12] = xr[2];
1758
1759           mcpart = mctrackmother->Particle();
1760           if(mcpart){
1761             datamc[15] = mcpart->GetUniqueID();
1762           }
1763           if(glabel>fMCEvent->GetNumberOfPrimaries()) {
1764             bgcategory = -2.;
1765             datamc[1] = bgcategory;
1766             glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1767             if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1768               datamc[3]=mctrackmother->PdgCode();
1769               datamc[4]=mctrackmother->Pt();
1770               if(TMath::Abs(mctrackmother->PdgCode())==310){
1771                bgcategory = -3.;
1772                glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
1773                if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1774                  datamc[5]=mctrackmother->PdgCode();
1775                  glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
1776                  if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1777                    datamc[6]=mctrackmother->PdgCode();
1778                  }
1779                }
1780               }
1781             }
1782           }
1783         }
1784      }
1785
1786
1787      if(fIsPbPb){
1788        Int_t centBin = GetWeightCentralityBin(fPerCentrality);
1789        if(centBin < 0)return 0.;
1790        weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
1791        for(int ii=0; ii<kBgPtBins; ii++){
1792          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1793            weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
1794            break;
1795          }
1796        }
1797      }
1798      else{
1799        weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
1800        for(int ii=0; ii<kBgPtBins; ii++){
1801          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1802            weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
1803            break;
1804          }
1805        }
1806      }
1807   }
1808
1809   datamc[13] = weightElecBg;
1810   datamc[16] = Double_t(fContainerStep);
1811
1812   datamc[7] = fHfeImpactR;
1813   datamc[8] = fHfeImpactnsigmaR;
1814
1815   datamc[19] = fRecPt;
1816   datamc[20] = fRecEta;
1817   datamc[21] = fRecPhi;
1818   datamc[22] = fLyrhit;
1819   datamc[23] = fLyrstat;
1820   datamc[24] = fCentrality;
1821
1822   if(fIsDebugStreamerON && fTreeStream){
1823    if(!iBgLevel){
1824     (*fTreeStream)<<"nonhfeQA"<<
1825         "mesonID="<<datamc[0]<<
1826         "bgcategory="<<datamc[1]<<
1827         "mesonPt="<<datamc[2]<<
1828         "mesonMomPdg="<<datamc[3]<<
1829         "mesonMomPt="<<datamc[4]<<
1830         "mesonGMomPdg="<<datamc[5]<<
1831         "mesonGGMomPdg="<<datamc[6]<<
1832         "eIPAbs="<<datamc[7]<<
1833         "eIPSig="<<datamc[8]<<
1834         "eR="<<datamc[9]<<
1835         "eZ="<<datamc[10]<<
1836         "mesonR="<<datamc[11]<<
1837         "mesonZ="<<datamc[12]<<
1838         "weightElecBg="<<datamc[13]<<
1839         "eUniqID="<<datamc[14]<<
1840         "mesonUniqID="<<datamc[15]<<
1841         "containerStep="<<datamc[16]<<
1842         "emcpt="<<datamc[17]<<
1843         "emceta="<<datamc[18]<<
1844         "erecpt="<<datamc[19]<<
1845         "ereceta="<<datamc[20]<<
1846         "erecphi="<<datamc[21]<<
1847         "itshit="<<datamc[22]<<
1848         "itsstat="<<datamc[23]<<
1849         "centrality="<<datamc[24]
1850         << "\n";
1851    }
1852   }
1853
1854   Double_t returnval = bgcategory*weightElecBg;
1855   if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
1856
1857   return returnval;
1858 }
1859
1860 //__________________________________________
1861 Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart){
1862   //
1863   // Wrapper to get the mother label
1864   //
1865   Int_t label = -1;
1866   TClass *type = mcpart->IsA();
1867   if(type == AliMCParticle::Class()){
1868     // ESD analysis
1869     const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
1870     label = emcpart->GetMother();
1871   } else if(type == AliAODMCParticle::Class()){
1872     // AOD analysis
1873     const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
1874     label = amcpart->GetMother();
1875   }
1876   return label;
1877 }
1878 //__________________________________________
1879 Int_t AliHFEmcQA::GetWeightCentralityBin(const Float_t percentile) const {
1880   //
1881   //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
1882   //
1883
1884   Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
1885   Int_t bin = -1;
1886   for(Int_t ibin = 0; ibin < 11; ibin++){
1887     if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
1888       bin = ibin;
1889       break;
1890     }
1891   }
1892   return bin;
1893 }
1894 //__________________________________________
1895 void AliHFEmcQA::AliHists::FillList(TList *l) const {
1896   //
1897   // Fill Histos into a list for output
1898   //
1899   if(fPdgCode) l->Add(fPdgCode);
1900   if(fPt) l->Add(fPt);
1901   if(fY) l->Add(fY);
1902   if(fEta) l->Add(fEta);
1903 }
1904
1905 //__________________________________________
1906 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
1907   //
1908   // Fill Histos into a list for output
1909   //
1910   if(fNq) l->Add(fNq);
1911   if(fProcessID) l->Add(fProcessID);
1912   if(fePtRatio) l->Add(fePtRatio);
1913   if(fPtCorr) l->Add(fPtCorr);
1914   if(fPtCorrDp) l->Add(fPtCorrDp);
1915   if(fPtCorrD0) l->Add(fPtCorrD0);
1916   if(fPtCorrDrest) l->Add(fPtCorrDrest);
1917   if(fDePtRatio) l->Add(fDePtRatio);
1918   if(feDistance) l->Add(feDistance);
1919   if(fDeDistance) l->Add(fDeDistance);
1920 }
1921