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
  36 #include <AliLog.h>
  37 #include <AliMCEvent.h>
  38 #include <AliGenEventHeader.h>
  39 #include <AliAODMCParticle.h>
  40 #include <AliStack.h>
  41
  42 #include "AliHFEmcQA.h"
  43 #include "AliHFEtools.h"
  44 #include "AliHFEcollection.h"
  45
  46 ClassImp(AliHFEmcQA)
  47
  48 //_______________________________________________________________________________________________
  49     AliHFEmcQA::AliHFEmcQA() :
  50          fMCEvent(NULL)
  51         ,fMCHeader(NULL)
  52         ,fMCArray(NULL)
  53         ,fQAhistos(NULL)
  54         ,fMCQACollection(NULL)
  55         ,fNparents(0)
  56         ,fCentrality(0)
  57         ,fIsPbPb(kFALSE)
  58         ,fIsppMultiBin(kFALSE)
  59 {
  60         // Default constructor
  61   for(Int_t mom = 0; mom < 9; mom++){
  62     fhD[mom] = NULL;
  63     fhDLogbin[mom] = NULL;
  64   }
  65   for(Int_t mom = 0; mom < 50; mom++){
  66     fHeavyQuark[mom] = NULL;
  67   }
  68   for(Int_t mom = 0; mom < 2; mom++){
  69     fIsHeavy[mom] = 0;
  70   }
  71   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
  72   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
  73 }
  74
  75 //_______________________________________________________________________________________________
  76 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
  77 TObject(p)
  78         ,fMCEvent(NULL)
  79         ,fMCHeader(NULL)
  80         ,fMCArray(NULL)
  81         ,fQAhistos(p.fQAhistos)
  82         ,fMCQACollection(p.fMCQACollection)
  83         ,fNparents(p.fNparents)
  84         ,fCentrality(0)
  85         ,fIsPbPb(kFALSE)
  86         ,fIsppMultiBin(kFALSE)
  87 {
  88         // Copy constructor
  89   for(Int_t mom = 0; mom < 9; mom++){
  90     fhD[mom] = NULL;
  91     fhDLogbin[mom] = NULL;
  92   }
  93   for(Int_t mom = 0; mom < 50; mom++){
  94     fHeavyQuark[mom] = NULL;
  95   }
  96   for(Int_t mom = 0; mom < 2; mom++){
  97     fIsHeavy[mom] = 0;
  98   }
  99   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 100   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 101 }
 102
 103 //_______________________________________________________________________________________________
 104 AliHFEmcQA&
 105 AliHFEmcQA::operator=(const AliHFEmcQA &)
 106 {
 107   // Assignment operator
 108
 109   AliInfo("Not yet implemented.");
 110   return *this;
 111 }
 112
 113 //_______________________________________________________________________________________________
 114 AliHFEmcQA::~AliHFEmcQA()
 115 {
 116         // Destructor
 117
 118         AliInfo("Analysis Done.");
 119 }
 120
 121 //_______________________________________________________________________________________________
 122 void AliHFEmcQA::PostAnalyze() const
 123 {
 124         //
 125         // Post analysis
 126         //
 127 }
 128
 129 //_______________________________________________________________________________________________
 130 void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
 131 {
 132    //
 133    // copy background weighting factors into data member
 134    //
 135
 136   memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 137   memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
 138 }
 139
 140 //__________________________________________
 141 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
 142 {
 143   //
 144   // make default histograms
 145   //
 146
 147   if(!qaList) return;
 148
 149   fQAhistos = qaList;
 150   fQAhistos->SetName("MCqa");
 151
 152   CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
 153   CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
 154   CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
 155   CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm
 156   CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
 157   CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
 158   CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm
 159   CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
 160   CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
 161
 162 // check D spectra
 163   TString kDspecies[9];
 164   kDspecies[0]="411";
 165   kDspecies[1]="421";
 166   kDspecies[2]="431";
 167   kDspecies[3]="4122";
 168   kDspecies[4]="4132";
 169   kDspecies[5]="4232";
 170   kDspecies[6]="4332";
 171   kDspecies[7]="413";
 172   kDspecies[8]="423";
 173
 174   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 175   Int_t iBin[2];
 176   iBin[0] = 44; // bins in pt
 177   iBin[1] = 23; // bins in pt
 178   Double_t* binEdges[1];
 179   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 180
 181   // bin size is chosen to consider ALICE D measurement
 182   Double_t xbins[15]={0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,8.0,12,16,20,30,40,50};
 183   Double_t ybins[10]={-7.5,-1.0,-0.9,-0.8,-0.5,0.5,0.8,0.9,1.0,7.5};
 184   TString hname;
 185   for (Int_t iDmeson=0; iDmeson<9; iDmeson++){
 186      hname = "Dmeson"+kDspecies[iDmeson];
 187      fhD[iDmeson] = new TH2F(hname,hname+";p_{T} (GeV/c)",14,xbins,9,ybins);
 188      hname = "DmesonLogbin"+kDspecies[iDmeson];
 189      fhDLogbin[iDmeson] = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
 190      if(fQAhistos) fQAhistos->Add(fhD[iDmeson]);
 191      if(fQAhistos) fQAhistos->Add(fhDLogbin[iDmeson]);
 192   }
 193
 194   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
 195
 196   Int_t kNcent;
 197   if(fIsPbPb) kNcent=11;
 198   else
 199   {
 200       if(fIsppMultiBin) kNcent=8;
 201       else kNcent = 1;
 202   }
 203
 204   fMCQACollection = new AliHFEcollection("TaskMCQA", "MC QA histos for meason pt spectra");
 205
 206   for(Int_t centbin=0; centbin<kNcent; centbin++)
 207   {
 208       fMCQACollection->CreateTH1Farray(Form("pionspectra_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
 209       fMCQACollection->CreateTH1Farray(Form("etaspectra_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
 210       fMCQACollection->CreateTH1Farray(Form("omegaspectra_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[1],kPtRange);
 211       fMCQACollection->CreateTH1Farray(Form("phispectra_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[1],kPtRange);
 212       fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
 213       fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
 214
 215       fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 216       fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 217       fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 218       fMCQACollection->CreateTH1F(Form("phispectraLog_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 219       fMCQACollection->CreateTH1F(Form("etapspectraLog_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 220       fMCQACollection->CreateTH1F(Form("rhospectraLog_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 221       fMCQACollection->CreateTH1F(Form("kaonspectraLog_centrbin%i",centbin), "kaon yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 222       fMCQACollection->CreateTH1F(Form("k0LspectraLog_centrbin%i",centbin), "k0L yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 223
 224       fMCQACollection->CreateTH1F(Form("piondaughters_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 225       fMCQACollection->CreateTH1F(Form("etadaughters_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 226       fMCQACollection->CreateTH1F(Form("omegadaughters_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 227       fMCQACollection->CreateTH1F(Form("phidaughters_centrbin%i",centbin), "phi yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 228       fMCQACollection->CreateTH1F(Form("etapdaughters_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 229       fMCQACollection->CreateTH1F(Form("rhodaughters_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
 230   }
 231
 232   fQAhistos->Add(fMCQACollection->GetList());
 233
 234 }
 235
 236 //__________________________________________
 237 void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
 238 {
 239   // create histograms
 240
 241   if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)) {
 242     AliDebug(1, "This task is only for heavy quark QA, return\n");
 243     return;
 244   }
 245   Int_t iq = kquark - kCharm;
 246
 247   TString kqTypeLabel[fgkqType];
 248   if (kquark == kCharm){
 249     kqTypeLabel[kQuark]="c";
 250     kqTypeLabel[kantiQuark]="cbar";
 251     kqTypeLabel[kHadron]="cHadron";
 252     kqTypeLabel[keHadron]="ceHadron";
 253     kqTypeLabel[kDeHadron]="nullHadron";
 254     kqTypeLabel[kElectron]="ce";
 255     kqTypeLabel[kElectron2nd]="nulle";
 256   } else if (kquark == kBeauty){
 257     kqTypeLabel[kQuark]="b";
 258     kqTypeLabel[kantiQuark]="bbar";
 259     kqTypeLabel[kHadron]="bHadron";
 260     kqTypeLabel[keHadron]="beHadron";
 261     kqTypeLabel[kDeHadron]="bDeHadron";
 262     kqTypeLabel[kElectron]="be";
 263     kqTypeLabel[kElectron2nd]="bce";
 264   } else if (kquark == kOthers){
 265     kqTypeLabel[kGamma-4]="gammae";
 266     kqTypeLabel[kPi0-4]="pi0e";
 267     kqTypeLabel[kElse-4]="elsee";
 268     kqTypeLabel[kMisID-4]="miside";
 269   }
 270
 271   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
 272   Int_t iBin[1];
 273   iBin[0] = 44; // bins in pt
 274   Double_t* binEdges[1];
 275   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 276
 277
 278   Double_t xcorrbin[501];
 279   for (int icorrbin = 0; icorrbin< 501; icorrbin++){
 280     xcorrbin[icorrbin]=icorrbin*0.1;
 281   }
 282
 283   TString hname;
 284   if(kquark == kOthers){
 285     for (Int_t iqType = 0; iqType < 4; iqType++ ){
 286        hname = hnopt+"Pt_"+kqTypeLabel[iqType];
 287        //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL
 288        fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25);
 289        hname = hnopt+"Y_"+kqTypeLabel[iqType];
 290        fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
 291        hname = hnopt+"Eta_"+kqTypeLabel[iqType];
 292        fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
 293        // Fill List
 294        if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
 295     }
 296     return;
 297   }
 298   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
 299      if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
 300      hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
 301      fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
 302      hname = hnopt+"Pt_"+kqTypeLabel[iqType];
 303      fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
 304      //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.,30.); // mj to compare with FONLL
 305      //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL
 306      hname = hnopt+"Y_"+kqTypeLabel[iqType];
 307      fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
 308      hname = hnopt+"Eta_"+kqTypeLabel[iqType];
 309      fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
 310      // Fill List
 311      if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
 312   }
 313
 314   if (icut == 0){
 315     hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
 316     fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5);
 317     hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
 318     fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
 319   }
 320   hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
 321   fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
 322   hname = hnopt+"PtCorr_"+kqTypeLabel[kQuark];
 323   fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
 324   //fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
 325   hname = hnopt+"PtCorrDp_"+kqTypeLabel[kQuark];
 326   fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
 327   //fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
 328   hname = hnopt+"PtCorrD0_"+kqTypeLabel[kQuark];
 329   fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
 330   //fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
 331   hname = hnopt+"PtCorrDrest_"+kqTypeLabel[kQuark];
 332   fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",500,xcorrbin,44,binEdges[0]);
 333   //fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",200,0,20,200,0,20);
 334   hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
 335   fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
 336   hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
 337   fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
 338   hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
 339   fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
 340   if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
 341 }
 342
 343 //__________________________________________
 344 void AliHFEmcQA::Init()
 345 {
 346   // called at begining every event
 347
 348   for (Int_t i=0; i<2; i++){
 349      fIsHeavy[i] = 0;
 350   }
 351
 352   fNparents = 7;
 353
 354   fParentSelect[0][0] =  411;
 355   fParentSelect[0][1] =  421;
 356   fParentSelect[0][2] =  431;
 357   fParentSelect[0][3] = 4122;
 358   fParentSelect[0][4] = 4132;
 359   fParentSelect[0][5] = 4232;
 360   fParentSelect[0][6] = 4332;
 361
 362   fParentSelect[1][0] =  511;
 363   fParentSelect[1][1] =  521;
 364   fParentSelect[1][2] =  531;
 365   fParentSelect[1][3] = 5122;
 366   fParentSelect[1][4] = 5132;
 367   fParentSelect[1][5] = 5232;
 368   fParentSelect[1][6] = 5332;
 369
 370 }
 371
 372 //__________________________________________
 373 void AliHFEmcQA::GetMesonKine()
 374 {
 375   //
 376   // get meson pt spectra
 377   //
 378
 379   AliVParticle *mctrack2 = NULL;
 380   AliMCParticle *mctrack0 = NULL;
 381   AliVParticle *mctrackdaugt= NULL;
 382   AliMCParticle *mctrackd= NULL;
 383   Int_t id1=0, id2=0;
 384
 385
 386   if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
 387
 388
 389   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
 390      if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
 391      TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
 392      if(!mcpart0) continue;
 393      mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
 394      if(!mctrack0) continue;
 395
 396      if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 397
 398      if(abs(mctrack0->PdgCode()) == 111) // pi0
 399        {
 400           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 401             fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
 402             fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 403           }
 404           id1=mctrack0->GetFirstDaughter();
 405           id2=mctrack0->GetLastDaughter();
 406           if(!((id2-id1)==2)) continue;
 407           for(int idx=id1; idx<=id2; idx++){
 408             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 409             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 410             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 411              fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
 412           }
 413        }
 414      else if(abs(mctrack0->PdgCode()) == 221) // eta
 415        {
 416           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 417             fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
 418             fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 419           }
 420           id1=mctrack0->GetFirstDaughter();
 421           id2=mctrack0->GetLastDaughter();
 422           if(!((id2-id1)==2||(id2-id1)==3)) continue;
 423           for(int idx=id1; idx<=id2; idx++){
 424             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 425             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 426             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 427              fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
 428           }
 429        }
 430      else if(abs(mctrack0->PdgCode()) == 223) // omega
 431        {
 432           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 433             fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
 434             fMCQACollection->Fill(Form("omegaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 435           }
 436           id1=mctrack0->GetFirstDaughter();
 437           id2=mctrack0->GetLastDaughter();
 438           if(!((id2-id1)==1||(id2-id1)==2)) continue;
 439           for(int idx=id1; idx<=id2; idx++){
 440             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 441             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 442             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 443              fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
 444           }
 445        }
 446      else if(abs(mctrack0->PdgCode()) == 333) // phi
 447        {
 448           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 449             fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
 450             fMCQACollection->Fill(Form("phispectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 451           }
 452           id1=mctrack0->GetFirstDaughter();
 453           id2=mctrack0->GetLastDaughter();
 454           if(!((id2-id1)==1)) continue;
 455           for(int idx=id1; idx<=id2; idx++){
 456             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 457             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 458             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 459              fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
 460           }
 461        }
 462      else if(abs(mctrack0->PdgCode()) == 331) // eta prime
 463        {
 464           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 465             fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
 466             fMCQACollection->Fill(Form("etapspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 467           }
 468           id1=mctrack0->GetFirstDaughter();
 469           id2=mctrack0->GetLastDaughter();
 470           if(!((id2-id1)==2||(id2-id1)==3)) continue;
 471           for(int idx=id1; idx<=id2; idx++){
 472             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 473             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 474             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 475              fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
 476           }
 477        }
 478      else if(abs(mctrack0->PdgCode()) == 113) // rho
 479        {
 480           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 481             fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
 482             fMCQACollection->Fill(Form("rhospectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 483           }
 484           id1=mctrack0->GetFirstDaughter();
 485           id2=mctrack0->GetLastDaughter();
 486           if(!((id2-id1)==1)) continue;
 487           for(int idx=id1; idx<=id2; idx++){
 488             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
 489             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
 490             if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
 491              fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
 492           }
 493        }
 494      else if(abs(mctrack0->PdgCode()) == 321) // kaon+-
 495        {
 496           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 497             fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 498           }
 499        }
 500      else if(abs(mctrack0->PdgCode()) == 130) // k0L
 501        {
 502           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
 503             fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
 504           }
 505        }
 506      }
 507
 508 }
 509 //__________________________________________
 510 void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
 511 {
 512   // get heavy quark kinematics
 513
 514     if (kquark != kCharm && kquark != kBeauty) {
 515       AliDebug(1, "This task is only for heavy quark QA, return\n");
 516       return;
 517     }
 518     Int_t iq = kquark - kCharm;
 519
 520     if (iTrack < 0 || !part) {
 521       AliDebug(1, "Stack label is negative or no mcparticle, return\n");
 522       return;
 523     }
 524
 525     AliMCParticle *mctrack = NULL;
 526     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 527
 528     // select heavy hadron or not fragmented heavy quark
 529     if ( int(partPdgcode/100.)==kquark || int(partPdgcode/1000.)==kquark || (partPdgcode==kquark && (part->GetNDaughters()==0 && iTrack>5)) ){
 530
 531       TParticle *partMother;
 532       Int_t iLabel;
 533
 534       if (partPdgcode == kquark){ // in case of not fragmented heavy quark
 535         partMother = part;
 536         iLabel = iTrack;
 537       } else{ // in case of heavy hadron, start to search for mother heavy parton
 538         iLabel = part->GetFirstMother();
 539         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 540         if (iLabel>-1) { partMother = mctrack->Particle(); }
 541         else {
 542           AliDebug(1, "Stack label is negative, return\n");
 543           return;
 544         }
 545       }
 546
 547       // heavy parton selection as a mother of heavy hadron
 548       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
 549       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
 550       // should I make a condition that partMother should be quark or diquark? -> not necessary
 551       if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
 552       //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
 553
 554         if ( abs(partMother->GetPdgCode()) != kquark ){
 555           // search fragmented partons in the same string
 556           Bool_t isSameString = kTRUE;
 557           for (Int_t i=1; i<fgkMaxIter; i++){
 558              iLabel = iLabel - 1;
 559              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 560              if (iLabel>-1) { partMother = mctrack->Particle(); }
 561              else {
 562                AliDebug(1, "Stack label is negative, return\n");
 563                return;
 564              }
 565              if ( abs(partMother->GetPdgCode()) == kquark ) break;
 566              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
 567              if (!isSameString) return;
 568           }
 569         }
 570         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
 571         if (abs(partMother->GetPdgCode()) != kquark) return;
 572
 573         if (fIsHeavy[iq] >= 50) return;
 574         fHeavyQuark[fIsHeavy[iq]] = partMother;
 575         fIsHeavy[iq]++;
 576
 577         // fill kinematics for heavy parton
 578         if (partMother->GetPdgCode() > 0) { // quark
 579           fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 580           fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
 581           fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
 582           fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
 583         } else{ // antiquark
 584           fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
 585           fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
 586           fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
 587           fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
 588         }
 589
 590       } // end of heavy parton slection loop
 591
 592     } // end of heavy hadron or quark selection
 593
 594 }
 595
 596 //__________________________________________
 597 void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
 598 {
 599   // end of event analysis
 600
 601   if (kquark != kCharm && kquark != kBeauty) {
 602     AliDebug(1, "This task is only for heavy quark QA, return\n");
 603     return;
 604   }
 605   Int_t iq = kquark - kCharm;
 606
 607
 608   // # of heavy quark per event
 609   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
 610   fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
 611
 612   Int_t motherID[fgkMaxGener];
 613   Int_t motherType[fgkMaxGener];
 614   Int_t motherLabel[fgkMaxGener];
 615   Int_t ancestorPdg[fgkMaxGener];
 616   Int_t ancestorLabel[fgkMaxGener];
 617
 618   for (Int_t i = 0; i < fgkMaxGener; i++){ // initialization
 619      motherID[i] = 0;
 620      motherType[i] = 0;
 621      motherLabel[i] = 0;
 622      ancestorPdg[i] = 0;
 623      ancestorLabel[i] = 0;
 624   }
 625
 626
 627   // check history of found heavy quarks
 628   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 629
 630      if(!fHeavyQuark[i]) return;
 631
 632      ancestorLabel[0] = i;
 633      ancestorPdg[0] = fHeavyQuark[i]->GetPdgCode();
 634      ancestorLabel[1] = fHeavyQuark[i]->GetFirstMother();
 635
 636      AliDebug(1,Form("pdg code= %d\n",ancestorPdg[0]));
 637      AliDebug(1,Form("ancestor label= %d\n",ancestorLabel[1]));
 638
 639      Int_t ig = 1;
 640      while (ancestorLabel[ig] != -1){
 641           // in case there is mother, get mother's pdg code and grandmother's label
 642           IdentifyMother(ancestorLabel[ig], ancestorPdg[ig], ancestorLabel[ig+1]);
 643           // if mother is still heavy, find again mother's ancestor
 644           if (ancestorPdg[ig-1] == ancestorPdg[ig]) {
 645             ig++;
 646             continue; // if it is from same heavy
 647           }
 648           // if the heavy's mother is not heavy, check the mother's label to know if it comes from inital or final parton shower
 649           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 650           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
 651           // if it is not the above case, something is strange
 652           ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
 653           break;
 654      }
 655      if (ancestorLabel[ig] == -1){ // from hard scattering
 656        HardScattering(kquark, motherID[i],motherType[i], motherLabel[i]);
 657      }
 658
 659   } // end of found heavy quark loop
 660
 661
 662   // check process type
 663   Int_t processID = 0;
 664   for (Int_t i = 0; i < fIsHeavy[iq]; i++){
 665      AliDebug(1,Form("Mother ID= %d type= %d label= %d\n",motherID[i],motherType[i],motherLabel[i]));
 666   }
 667
 668
 669   Int_t nheavypair = Int_t(fIsHeavy[iq]/2.);
 670   for (Int_t ipair = 0; ipair < nheavypair; ipair++){
 671
 672      Int_t id1 = ipair*2;
 673      Int_t id2 = ipair*2 + 1;
 674
 675      if (motherType[id1] == 2 && motherType[id2] == 2){
 676        if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
 677        else processID = -9;
 678      }
 679      else if (motherType[id1] == -1 && motherType[id2] == -1) {
 680        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
 681          if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
 682          else processID = kPairCreationFromq; // q-qbar pair creation
 683        }
 684        else processID = -8;
 685      }
 686      else if (motherType[id1] == -1 || motherType[id2] == -1) {
 687        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
 688          if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation
 689          else processID = kLightQuarkShower;
 690        }
 691        else processID = -7;
 692      }
 693      else if (motherType[id1] == -2 || motherType[id2] == -2) {
 694        if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
 695        else processID = -6;
 696
 697      }
 698      else processID = -5;
 699
 700      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
 701      else fHistComm[iq][0].fProcessID->Fill(processID);
 702      AliDebug(1,Form("Process ID = %d\n",processID));
 703   } // end of # heavy quark pair loop
 704
 705 }
 706
 707 //__________________________________________
 708 void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
 709 {
 710     // decay electron kinematics
 711
 712     if (kquark != kCharm && kquark != kBeauty) {
 713       AliDebug(1, "This task is only for heavy quark QA, return\n");
 714       return;
 715     }
 716     Int_t iq = kquark - kCharm;
 717
 718     if(!mcpart){
 719       AliDebug(1, "no mc particle, return\n");
 720       return;
 721     }
 722
 723     Int_t iLabel = mcpart->GetFirstMother();
 724     if (iLabel<0){
 725       AliDebug(1, "Stack label is negative, return\n");
 726       return;
 727     }
 728
 729     TParticle *partCopy = mcpart;
 730     Int_t pdgcode = mcpart->GetPdgCode();
 731     Int_t pdgcodeCopy = pdgcode;
 732
 733     AliMCParticle *mctrack = NULL;
 734
 735     // if the mother is charmed hadron
 736     Bool_t isDirectCharm = kFALSE;
 737     if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
 738
 739           // iterate until you find B hadron as a mother or become top ancester
 740           for (Int_t i=1; i<fgkMaxIter; i++){
 741
 742              Int_t jLabel = mcpart->GetFirstMother();
 743              if (jLabel == -1){
 744                isDirectCharm = kTRUE;
 745                break; // if there is no ancester
 746              }
 747              if (jLabel < 0){ // safety protection
 748                AliDebug(1, "Stack label is negative, return\n");
 749                return;
 750              }
 751              // if there is an ancester
 752              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
 753              TParticle* mother = mctrack->Particle();
 754              Int_t motherPDG = mother->GetPdgCode();
 755
 756              for (Int_t j=0; j<fNparents; j++){
 757                 if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
 758              }
 759
 760              mcpart = mother;
 761           } // end of iteration
 762     } // end of if
 763     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 764          for (Int_t i=0; i<fNparents; i++){
 765             if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
 766
 767               // fill hadron kinematics
 768               fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
 769               fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
 770               fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
 771               fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
 772
 773               if(iq==0) {
 774                fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 775                if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[i]->Fill(partCopy->Pt());
 776               }
 777             }
 778          }
 779          // I also want to store D* info to compare with D* measurement
 780          if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
 781                fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 782                if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[7]->Fill(partCopy->Pt());
 783          }
 784          if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
 785                fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
 786                if(TMath::Abs(AliHFEtools::GetRapidity(partCopy))<0.5) fhDLogbin[8]->Fill(partCopy->Pt());
 787          }
 788     } // end of if
 789 }
 790
 791 //__________________________________________
 792 void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
 793 {
 794     // decay electron kinematics
 795
 796     if (!(kquark == kCharm || kquark == kBeauty || kquark == kOthers)){
 797       AliDebug(1, "This task is only for heavy quark QA, return\n");
 798       return;
 799     }
 800     Int_t iq = kquark - kCharm;
 801     Bool_t isFinalOpenCharm = kFALSE;
 802
 803     if(!mcpart){
 804       AliDebug(1, "no mcparticle, return\n");
 805       return;
 806     }
 807
 808     if(kquark==kOthers){
 809       Int_t esource = -1;
 810       if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
 811       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
 812       if(esource==0|| esource==1 || esource==2 || esource==3){
 813         fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
 814         fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 815         fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
 816         return;
 817       }
 818       else {
 819         AliDebug(1, "e source is out of defined ranges, return\n");
 820         return;
 821       }
 822     }
 823
 824     if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 825
 826     Int_t iLabel = mcpart->GetFirstMother();
 827     if (iLabel<0){
 828       AliDebug(1, "Stack label is negative, return\n");
 829       return;
 830     }
 831
 832     AliMCParticle *mctrack = NULL;
 833     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return;
 834     TParticle *partMother = mctrack->Particle();
 835     TParticle *partMotherCopy = partMother;
 836     Int_t maPdgcode = partMother->GetPdgCode();
 837     Int_t maPdgcodeCopy = maPdgcode;
 838
 839     // get mc primary vertex
 840     /*
 841     TArrayF mcPrimVtx;
 842     if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
 843
 844     // get electron production vertex
 845     TLorentzVector ePoint;
 846     mcpart->ProductionVertex(ePoint);
 847
 848     // calculated production vertex to primary vertex (in xy plane)
 849     Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
 850     */
 851     Float_t decayLxy = 0;
 852
 853     // if the mother is charmed hadron
 854     Bool_t isMotherDirectCharm = kFALSE;
 855     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 856
 857          for (Int_t i=0; i<fNparents; i++){
 858             if (abs(maPdgcode)==fParentSelect[0][i]){
 859               isFinalOpenCharm = kTRUE;
 860             }
 861          }
 862          if (!isFinalOpenCharm) return ;
 863
 864           // iterate until you find B hadron as a mother or become top ancester
 865           for (Int_t i=1; i<fgkMaxIter; i++){
 866
 867              Int_t jLabel = partMother->GetFirstMother();
 868              if (jLabel == -1){
 869                isMotherDirectCharm = kTRUE;
 870                break; // if there is no ancester
 871              }
 872              if (jLabel < 0){ // safety protection
 873                AliDebug(1, "Stack label is negative, return\n");
 874                return;
 875              }
 876
 877              // if there is an ancester
 878              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return;
 879              TParticle* grandMa = mctrack->Particle();
 880              Int_t grandMaPDG = grandMa->GetPdgCode();
 881
 882              for (Int_t j=0; j<fNparents; j++){
 883                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 884
 885                   if (kquark == kCharm) return;
 886                   // fill electron kinematics
 887                   fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 888                   fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
 889                   fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 890                   fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
 891
 892                   // fill mother hadron kinematics
 893                   fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
 894                   fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
 895                   fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
 896                   fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
 897
 898                   // ratio between pT of electron and pT of mother B hadron
 899                   if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 900
 901                   // distance between electron production point and primary vertex
 902                   fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
 903                   return;
 904                 }
 905              }
 906
 907              partMother = grandMa;
 908           } // end of iteration
 909     } // end of if
 910     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
 911          for (Int_t i=0; i<fNparents; i++){
 912             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 913
 914 //mj weighting to consider measured spectra!!!
 915               Double_t mpt=partMotherCopy->Pt();
 916               Double_t wfactor=2*(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225)); // 2* considering in pythia I used particle + antiparticle differently from the measurement
 917               //Double_t wfactor=(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225));
 918               // fill electron kinematics
 919               if(iq==0){
 920                 fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode(),wfactor);
 921                 fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt(),wfactor);
 922                 fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart),wfactor);
 923                 fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta(),wfactor);
 924
 925                 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 926                 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 927                 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 928                 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 929               }
 930               else{
 931                 fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
 932                 fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
 933                 fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
 934                 fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
 935               }
 936 //--------------
 937               // fill mother hadron kinematics
 938               fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
 939               fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
 940               fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
 941               fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
 942
 943               // ratio between pT of electron and pT of mother B or direct D hadron
 944               if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
 945               fHistComm[iq][icut].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
 946               if(TMath::Abs(partMotherCopy->GetPdgCode())==411) fHistComm[iq][icut].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
 947               else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) fHistComm[iq][icut].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
 948               else fHistComm[iq][icut].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
 949
 950               // distance between electron production point and primary vertex
 951               fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
 952             }
 953          }
 954     } // end of if
 955 }
 956
 957 //____________________________________________________________________
 958 void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
 959 {
 960   // decay electron kinematics
 961
 962   if (kquark != kCharm && kquark != kBeauty) {
 963     AliDebug(1, "This task is only for heavy quark QA, return\n");
 964     return;
 965   }
 966
 967   Int_t iq = kquark - kCharm;
 968   Bool_t isFinalOpenCharm = kFALSE;
 969
 970   if(!mcpart){
 971     AliDebug(1, "no mcparticle, return\n");
 972     return;
 973   }
 974
 975   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 976
 977   // mother
 978   Int_t iLabel = mcpart->GetMother();
 979   if (iLabel<0){
 980     AliDebug(1, "Stack label is negative, return\n");
 981     return;
 982   }
 983
 984   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
 985   AliAODMCParticle *partMotherCopy = partMother;
 986   Int_t maPdgcode = partMother->GetPdgCode();
 987   Int_t maPdgcodeCopy = maPdgcode;
 988
 989   Bool_t isMotherDirectCharm = kFALSE;
 990   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 991
 992     for (Int_t i=0; i<fNparents; i++){
 993        if (abs(maPdgcode)==fParentSelect[0][i]){
 994          isFinalOpenCharm = kTRUE;
 995        }
 996     }
 997     if (!isFinalOpenCharm) return;
 998
 999     for (Int_t i=1; i<fgkMaxIter; i++){
1000
1001        Int_t jLabel = partMother->GetMother();
1002        if (jLabel == -1){
1003          isMotherDirectCharm = kTRUE;
1004          break; // if there is no ancester
1005        }
1006        if (jLabel < 0){ // safety protection
1007          AliDebug(1, "Stack label is negative, return\n");
1008          return;
1009        }
1010
1011        // if there is an ancester
1012        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1013        Int_t grandMaPDG = grandMa->GetPdgCode();
1014
1015        for (Int_t j=0; j<fNparents; j++){
1016           if (abs(grandMaPDG)==fParentSelect[1][j]){
1017
1018             if (kquark == kCharm) return;
1019             // fill electron kinematics
1020             fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
1021             fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
1022             fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1023             fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
1024
1025             // fill mother hadron kinematics
1026             fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
1027             fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
1028             fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
1029             fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
1030
1031             return;
1032           }
1033        }
1034
1035        partMother = grandMa;
1036     } // end of iteration
1037   } // end of if
1038   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
1039     for (Int_t i=0; i<fNparents; i++){
1040        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
1041
1042          // fill electron kinematics
1043          fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
1044          fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
1045          fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
1046          fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
1047
1048          // fill mother hadron kinematics
1049          fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
1050          fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
1051          fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
1052          fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
1053
1054        }
1055     }
1056   } // end of if
1057
1058 }
1059
1060 //__________________________________________
1061 void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel)
1062 {
1063        // find mother pdg code and label
1064
1065        if (motherlabel < 0) {
1066          AliDebug(1, "Stack label is negative, return\n");
1067          return;
1068        }
1069        AliMCParticle *mctrack = NULL;
1070        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherlabel))))) return;
1071        TParticle *heavysMother = mctrack->Particle();
1072        motherpdg = heavysMother->GetPdgCode();
1073        grandmotherlabel = heavysMother->GetFirstMother();
1074        AliDebug(1,Form("ancestor pdg code= %d\n",motherpdg));
1075 }
1076
1077 //__________________________________________
1078 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1079 {
1080        // mothertype -1 means this heavy quark coming from hard vertex
1081
1082        AliMCParticle *mctrack1 = NULL;
1083        AliMCParticle *mctrack2 = NULL;
1084        if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(4))))) return;
1085        if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(5))))) return;
1086        TParticle *afterinitialrad1  = mctrack1->Particle();
1087        TParticle *afterinitialrad2  = mctrack2->Particle();
1088
1089        motherlabel = -1;
1090
1091        if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
1092          AliDebug(1,"heavy from gluon gluon pair creation!\n");
1093          mothertype = -1;
1094          motherID = fgkGluon;
1095        }
1096        else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
1097          AliDebug(1,"heavy from flavor exitation!\n");
1098          mothertype = -1;
1099          motherID = kquark;
1100        }
1101        else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
1102          AliDebug(1,"heavy from q-qbar pair creation!\n");
1103          mothertype = -1;
1104          motherID = 1;
1105        }
1106        else {
1107          AliDebug(1,"something strange!\n");
1108          mothertype = -999;
1109          motherlabel = -999;
1110          motherID = -999;
1111        }
1112 }
1113
1114 //__________________________________________
1115 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1116 {
1117        // mothertype -2 means this heavy quark coming from initial state
1118
1119        AliMCParticle *mctrack = NULL;
1120        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
1121          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1122          TParticle *heavysMother = mctrack->Particle();
1123          motherID = heavysMother->GetPdgCode();
1124          mothertype = -2; // there is mother before initial state radiation
1125          motherlabel = inputmotherlabel;
1126          AliDebug(1,"initial parton shower! \n");
1127
1128          return kTRUE;
1129        }
1130
1131        return kFALSE;
1132 }
1133
1134 //__________________________________________
1135 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1136 {
1137        // mothertype 2 means this heavy quark coming from final state
1138
1139        AliMCParticle *mctrack = NULL;
1140        if (inputmotherlabel > 5){ // mother exist after hard scattering
1141          if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(inputmotherlabel))))) return kFALSE;
1142          TParticle *heavysMother = mctrack->Particle();
1143          motherID = heavysMother->GetPdgCode();
1144          mothertype = 2; //
1145          motherlabel = inputmotherlabel;
1146          AliDebug(1,Form("heavy quark from %d after hard scattering! \n",motherID));
1147
1148          return kTRUE;
1149        }
1150        return kFALSE;
1151 }
1152
1153 //__________________________________________
1154 void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
1155 {
1156       // mark strange behavior
1157
1158        mothertype = -888;
1159        motherlabel = -888;
1160        motherID = -888;
1161        AliDebug(1,"something strange!\n");
1162 }
1163
1164 //__________________________________________
1165 Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
1166 {
1167   // decay particle's origin
1168
1169   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1170
1171   Int_t origin = -1;
1172   Bool_t isFinalOpenCharm = kFALSE;
1173
1174   if(!mcpart){
1175     AliDebug(1, "Stack label is negative or no mcparticle, return\n");
1176     return -1;
1177   }
1178
1179   // mother
1180   Int_t iLabel = mcpart->GetMother();
1181   if (iLabel<0){
1182     AliDebug(1, "Stack label is negative, return\n");
1183     return -1;
1184   }
1185
1186   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
1187   Int_t maPdgcode = partMother->GetPdgCode();
1188
1189   // if the mother is charmed hadron
1190   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1191
1192     for (Int_t i=0; i<fNparents; i++){
1193        if (abs(maPdgcode)==fParentSelect[0][i]){
1194          isFinalOpenCharm = kTRUE;
1195        }
1196     }
1197     if (!isFinalOpenCharm) return -1;
1198
1199     // iterate until you find B hadron as a mother or become top ancester
1200     for (Int_t i=1; i<fgkMaxIter; i++){
1201
1202        Int_t jLabel = partMother->GetMother();
1203        if (jLabel == -1){
1204          origin = kDirectCharm;
1205          return origin;
1206        }
1207        if (jLabel < 0){ // safety protection
1208          AliDebug(1, "Stack label is negative, return\n");
1209          return -1;
1210        }
1211
1212        // if there is an ancester
1213        AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
1214        Int_t grandMaPDG = grandMa->GetPdgCode();
1215
1216        for (Int_t j=0; j<fNparents; j++){
1217           if (abs(grandMaPDG)==fParentSelect[1][j]){
1218             origin = kBeautyCharm;
1219             return origin;
1220           }
1221        }
1222
1223        partMother = grandMa;
1224     } // end of iteration
1225   } // end of if
1226   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1227     for (Int_t i=0; i<fNparents; i++){
1228        if (abs(maPdgcode)==fParentSelect[1][i]){
1229          origin = kDirectBeauty;
1230          return origin;
1231        }
1232     }
1233   } // end of if
1234   else if ( abs(maPdgcode) == 22 ) {
1235     origin = kGamma;
1236     return origin;
1237   } // end of if
1238   else if ( abs(maPdgcode) == 111 ) {
1239     origin = kPi0;
1240     return origin;
1241   } // end of if
1242
1243   return origin;
1244 }
1245
1246 //__________________________________________
1247 Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
1248 {
1249   // decay particle's origin
1250
1251   //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
1252
1253   Int_t origin = -1;
1254   Bool_t isFinalOpenCharm = kFALSE;
1255
1256   if(!mcpart){
1257     AliDebug(1, "no mcparticle, return\n");
1258     return -1;
1259   }
1260
1261   Int_t iLabel = mcpart->GetFirstMother();
1262   if (iLabel<0){
1263     AliDebug(1, "Stack label is negative, return\n");
1264     return -1;
1265   }
1266
1267   AliMCParticle *mctrack = NULL;
1268   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1269   TParticle *partMother = mctrack->Particle();
1270   Int_t maPdgcode = partMother->GetPdgCode();
1271
1272    // if the mother is charmed hadron
1273    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1274
1275      for (Int_t i=0; i<fNparents; i++){
1276         if (abs(maPdgcode)==fParentSelect[0][i]){
1277           isFinalOpenCharm = kTRUE;
1278         }
1279      }
1280      if (!isFinalOpenCharm) return -1;
1281
1282      // iterate until you find B hadron as a mother or become top ancester
1283      for (Int_t i=1; i<fgkMaxIter; i++){
1284
1285         Int_t jLabel = partMother->GetFirstMother();
1286         if (jLabel == -1){
1287           origin = kDirectCharm;
1288           return origin;
1289         }
1290         if (jLabel < 0){ // safety protection
1291           AliDebug(1, "Stack label is negative, return\n");
1292           return -1;
1293         }
1294
1295         // if there is an ancester
1296         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1297         TParticle* grandMa = mctrack->Particle();
1298         Int_t grandMaPDG = grandMa->GetPdgCode();
1299
1300         for (Int_t j=0; j<fNparents; j++){
1301            if (abs(grandMaPDG)==fParentSelect[1][j]){
1302              origin = kBeautyCharm;
1303              return origin;
1304            }
1305         }
1306
1307         partMother = grandMa;
1308      } // end of iteration
1309    } // end of if
1310    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1311      for (Int_t i=0; i<fNparents; i++){
1312         if (abs(maPdgcode)==fParentSelect[1][i]){
1313           origin = kDirectBeauty;
1314           return origin;
1315         }
1316      }
1317    } // end of if
1318    else if ( abs(maPdgcode) == 22 ) {
1319      origin = kGamma;
1320      return origin;
1321    } // end of if
1322    else if ( abs(maPdgcode) == 111 ) {
1323      origin = kPi0;
1324      return origin;
1325    } // end of if
1326    else origin = kElse;
1327
1328    return origin;
1329 }
1330
1331 //__________________________________________
1332 Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
1333 {
1334   // decay particle's origin
1335
1336   if(!mcpart){
1337     AliDebug(1, "no mcparticle, return\n");
1338     return -1;
1339   }
1340
1341   if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
1342
1343   Int_t origin = -1;
1344   Bool_t isFinalOpenCharm = kFALSE;
1345
1346   Int_t iLabel = mcpart->GetFirstMother();
1347   if (iLabel<0){
1348     AliDebug(1, "Stack label is negative, return\n");
1349     return -1;
1350   }
1351
1352   AliMCParticle *mctrack = NULL;
1353   Int_t tmpMomLabel=0;
1354   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1355   TParticle *partMother = mctrack->Particle();
1356   TParticle *partMotherCopy = mctrack->Particle();
1357   Int_t maPdgcode = partMother->GetPdgCode();
1358
1359    // if the mother is charmed hadron
1360    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1361
1362      for (Int_t i=0; i<fNparents; i++){
1363         if (abs(maPdgcode)==fParentSelect[0][i]){
1364           isFinalOpenCharm = kTRUE;
1365         }
1366      }
1367      if (!isFinalOpenCharm) return -1;
1368
1369      // iterate until you find B hadron as a mother or become top ancester
1370      for (Int_t i=1; i<fgkMaxIter; i++){
1371
1372         Int_t jLabel = partMother->GetFirstMother();
1373         if (jLabel == -1){
1374           origin = kDirectCharm;
1375           return origin;
1376         }
1377         if (jLabel < 0){ // safety protection
1378           AliDebug(1, "Stack label is negative, return\n");
1379           return -1;
1380         }
1381
1382         // if there is an ancester
1383         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1384         TParticle* grandMa = mctrack->Particle();
1385         Int_t grandMaPDG = grandMa->GetPdgCode();
1386
1387         for (Int_t j=0; j<fNparents; j++){
1388            if (abs(grandMaPDG)==fParentSelect[1][j]){
1389              origin = kBeautyCharm;
1390              return origin;
1391            }
1392         }
1393
1394         partMother = grandMa;
1395      } // end of iteration
1396    } // end of if
1397    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1398      for (Int_t i=0; i<fNparents; i++){
1399         if (abs(maPdgcode)==fParentSelect[1][i]){
1400           origin = kDirectBeauty;
1401           return origin;
1402         }
1403      }
1404    } // end of if
1405    else if ( abs(maPdgcode) == 22 ) { //conversion
1406
1407      tmpMomLabel = partMotherCopy->GetFirstMother();
1408      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
1409      partMother = mctrack->Particle();
1410      maPdgcode = partMother->GetPdgCode();
1411      if ( abs(maPdgcode) == 111 ) {
1412        origin = kGammaPi0;
1413        return origin;
1414      }
1415      else if ( abs(maPdgcode) == 221 ) {
1416        origin = kGammaEta;
1417        return origin;
1418      }
1419      else if ( abs(maPdgcode) == 223 ) {
1420        origin = kGammaOmega;
1421        return origin;
1422      }
1423      else if ( abs(maPdgcode) == 333 ) {
1424        origin = kGammaPhi;
1425        return origin;
1426      }
1427      else if ( abs(maPdgcode) == 331 ) {
1428        origin = kGammaEtaPrime;
1429        return origin;
1430      }
1431      else if ( abs(maPdgcode) == 113 ) {
1432        origin = kGammaRho0;
1433        return origin;
1434      }
1435      else origin = kElse;
1436      //origin = kGamma; // finer category above
1437      return origin;
1438
1439    } // end of if
1440    else if ( abs(maPdgcode) == 111 ) {
1441      origin = kPi0;
1442      return origin;
1443    } // end of if
1444    else if ( abs(maPdgcode) == 221 ) {
1445      origin = kEta;
1446      return origin;
1447    } // end of if
1448    else if ( abs(maPdgcode) == 223 ) {
1449      origin = kOmega;
1450      return origin;
1451    } // end of if
1452    else if ( abs(maPdgcode) == 333 ) {
1453      origin = kPhi;
1454      return origin;
1455    } // end of if
1456    else if ( abs(maPdgcode) == 331 ) {
1457      origin = kEtaPrime;
1458      return origin;
1459    } // end of if
1460    else if ( abs(maPdgcode) == 113 ) {
1461      origin = kRho0;
1462      return origin;
1463    } // end of if
1464    else origin = kElse;
1465
1466    return origin;
1467 }
1468 //__________________________________________
1469 Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
1470   //
1471   // 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)
1472   //
1473   AliMCParticle *mctrackmother = NULL;
1474   Double_t weightElecBg = 0.;
1475   Double_t mesonPt = 0.;
1476   Double_t bgcategory = 0.;
1477   Int_t mArr = -1;
1478   Int_t mesonID = GetElecSource(mctrack->Particle());
1479   if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0;                //pion
1480   else if(mesonID==kGammaEta || mesonID==kEta) mArr=1;           //eta
1481   else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2;       //omega
1482   else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3;           //phi
1483   else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
1484   else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
1485
1486   if(!(mArr<0)){
1487      if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering
1488         Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
1489         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1490           glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's label
1491           if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1492             mesonPt = mctrackmother->Pt(); //meson pt
1493             bgcategory = 1.;
1494           }
1495         }
1496      }
1497      else{ // nonHFE except for the conversion electron
1498         Int_t glabel=TMath::Abs(mctrack->GetMother());
1499         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1500           mesonPt=mctrackmother->Pt(); //meson pt
1501           bgcategory = -1.;
1502         }
1503      }
1504
1505      if(fIsPbPb){
1506        if(fCentrality < 0)return 0.;
1507        weightElecBg=fElecBackgroundFactor[iBgLevel][fCentrality][mArr][kBgPtBins-1];
1508        for(int ii=0; ii<kBgPtBins; ii++){
1509          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1510            weightElecBg = fElecBackgroundFactor[iBgLevel][fCentrality][mArr][ii];
1511            break;
1512          }
1513        }
1514      }
1515      else{
1516        weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
1517        for(int ii=0; ii<kBgPtBins; ii++){
1518          if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
1519            weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
1520            break;
1521          }
1522        }
1523      }
1524   }
1525   return bgcategory*weightElecBg;
1526 }
1527
1528 //__________________________________________
1529 void AliHFEmcQA::AliHists::FillList(TList *l) const {
1530   //
1531   // Fill Histos into a list for output
1532   //
1533   if(fPdgCode) l->Add(fPdgCode);
1534   if(fPt) l->Add(fPt);
1535   if(fY) l->Add(fY);
1536   if(fEta) l->Add(fEta);
1537 }
1538
1539 //__________________________________________
1540 void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
1541   //
1542   // Fill Histos into a list for output
1543   //
1544   if(fNq) l->Add(fNq);
1545   if(fProcessID) l->Add(fProcessID);
1546   if(fePtRatio) l->Add(fePtRatio);
1547   if(fPtCorr) l->Add(fPtCorr);
1548   if(fPtCorrDp) l->Add(fPtCorrDp);
1549   if(fPtCorrD0) l->Add(fPtCorrD0);
1550   if(fPtCorrDrest) l->Add(fPtCorrDrest);
1551   if(fDePtRatio) l->Add(fDePtRatio);
1552   if(feDistance) l->Add(feDistance);
1553   if(fDeDistance) l->Add(fDeDistance);
1554 }
1555