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