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