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