]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFEmcQA.cxx
Fix of sigmaZ for crossing tracklets from Alex
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEmcQA.cxx
index 418eccc121eb5c1f87f9d565668ab8c3aed344ed..449624ba48c3c82f356f16c8d94eb624218636cd 100644 (file)
@@ -32,6 +32,7 @@
 #include <TH2F.h>
 #include <TList.h>
 #include <TParticle.h>
+#include "TTreeStream.h"
 
 #include <AliLog.h>
 #include <AliMCEvent.h>
 ClassImp(AliHFEmcQA)
 
 //_______________________________________________________________________________________________
-    AliHFEmcQA::AliHFEmcQA() :
-         fMCEvent(NULL)
-        ,fMCHeader(NULL)
-        ,fMCArray(NULL)
-        ,fQAhistos(NULL)
-        ,fMCQACollection(NULL)
-        ,fNparents(0)
-        ,fCentrality(0)
-        ,fIsPbPb(kFALSE)
-        ,fIsppMultiBin(kFALSE)
+AliHFEmcQA::AliHFEmcQA() :
+fMCEvent(NULL)
+  ,fMCHeader(NULL)
+  ,fMCArray(NULL)
+  ,fQAhistos(NULL)
+  ,fMCQACollection(NULL)
+  ,fNparents(0)
+  ,fCentrality(0)
+  ,fPerCentrality(-1)
+  ,fIsPbPb(kFALSE)
+  ,fIsppMultiBin(kFALSE)
+  ,fContainerStep(0)
+  ,fIsDebugStreamerON(kFALSE)
+  ,fRecPt(-999)
+  ,fRecEta(-999)
+  ,fRecPhi(-999)
+  ,fLyrhit(0)
+  ,fLyrstat(0)
+  ,fHfeImpactR(-999)
+  ,fHfeImpactnsigmaR(-999)
+  ,fTreeStream(NULL)
 {
-        // Default constructor
+  // Default constructor
   for(Int_t mom = 0; mom < 9; mom++){
     fhD[mom] = NULL;
   }
@@ -73,18 +85,29 @@ ClassImp(AliHFEmcQA)
 
 //_______________________________________________________________________________________________
 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
-TObject(p)
-        ,fMCEvent(NULL)
-        ,fMCHeader(NULL)
-        ,fMCArray(NULL)
-        ,fQAhistos(p.fQAhistos)
-        ,fMCQACollection(p.fMCQACollection)
-        ,fNparents(p.fNparents)
-        ,fCentrality(0)
-        ,fIsPbPb(kFALSE)
-        ,fIsppMultiBin(kFALSE)
+  TObject(p)
+  ,fMCEvent(NULL)
+  ,fMCHeader(NULL)
+  ,fMCArray(NULL)
+  ,fQAhistos(p.fQAhistos)
+  ,fMCQACollection(p.fMCQACollection)
+  ,fNparents(p.fNparents)
+  ,fCentrality(0)
+  ,fPerCentrality(-1)
+  ,fIsPbPb(kFALSE)
+  ,fIsppMultiBin(kFALSE)
+  ,fContainerStep(0)
+  ,fIsDebugStreamerON(kFALSE)
+  ,fRecPt(-999)
+  ,fRecEta(-999)
+  ,fRecPhi(-999)
+  ,fLyrhit(0)
+  ,fLyrstat(0)
+  ,fHfeImpactR(0)
+  ,fHfeImpactnsigmaR(0)
+  ,fTreeStream(NULL)
 {
-        // Copy constructor
+  // Copy constructor
   for(Int_t mom = 0; mom < 9; mom++){
     fhD[mom] = NULL;
   }
@@ -97,13 +120,12 @@ TObject(p)
   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 }
-
 //_______________________________________________________________________________________________
 AliHFEmcQA&
 AliHFEmcQA::operator=(const AliHFEmcQA &)
 {
   // Assignment operator
-
+  
   AliInfo("Not yet implemented.");
   return *this;
 }
@@ -111,30 +133,28 @@ AliHFEmcQA::operator=(const AliHFEmcQA &)
 //_______________________________________________________________________________________________
 AliHFEmcQA::~AliHFEmcQA()
 {
-        // Destructor
-
-        AliInfo("Analysis Done.");
+  // Destructor
+  
+  if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
+  AliInfo("Analysis Done.");
 }
-
 //_______________________________________________________________________________________________
 void AliHFEmcQA::PostAnalyze() const
 {
-        //
-        // Post analysis
-        //
+  //
+  // Post analysis
+  //
 }
-
 //_______________________________________________________________________________________________
 void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
 {
-   //
-   // copy background weighting factors into data member
-   //
+  //
+  // copy background weighting factors into data member
+  //
   
   memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
   memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
 }
-
 //__________________________________________
 void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
 {      
@@ -143,21 +163,15 @@ void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
   //
   
   if(!qaList) return;
-
+  
   fQAhistos = qaList;
   fQAhistos->SetName("MCqa");
-
-  CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
-  CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
-// prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
+  
+  CreateHistograms(AliHFEmcQA::kCharm);               // create histograms for charm
+  CreateHistograms(AliHFEmcQA::kBeauty);              // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers);              // create histograms for beauty
+  
+  // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
   const Int_t nbspecies = 9;
   TString kDspecies[nbspecies];
   kDspecies[0]="411";   //D+
@@ -210,6 +224,15 @@ void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
       fMCQACollection->CreateTH1Farray(Form("etapspectra_centrbin%i",centbin), "etap yields: MC p_{t} ", iBin[1],kPtRange);
       fMCQACollection->CreateTH1Farray(Form("rhospectra_centrbin%i",centbin), "rho yields: MC p_{t} ", iBin[1],kPtRange);
 
+      fMCQACollection->CreateTH1Farray(Form("pionspectrapri_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
+      fMCQACollection->CreateTH1Farray(Form("etaspectrapri_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
+      fMCQACollection->CreateTH1Farray(Form("pionspectrasec_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
+      fMCQACollection->CreateTH1Farray(Form("etaspectrasec_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
+      fMCQACollection->CreateTH1Farray(Form("pionspectraLogpri_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
+      fMCQACollection->CreateTH1Farray(Form("etaspectraLogpri_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
+      fMCQACollection->CreateTH1Farray(Form("pionspectraLogsec_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[1],kPtRange);
+      fMCQACollection->CreateTH1Farray(Form("etaspectraLogsec_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[1],kPtRange);
+
       fMCQACollection->CreateTH1F(Form("pionspectraLog_centrbin%i",centbin), "pion yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
       fMCQACollection->CreateTH1F(Form("etaspectraLog_centrbin%i",centbin), "eta yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
       fMCQACollection->CreateTH1F(Form("omegaspectraLog_centrbin%i",centbin), "omega yields: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
@@ -229,10 +252,14 @@ void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
 
   fQAhistos->Add(fMCQACollection->GetList());
 
+  if(!fTreeStream && fIsDebugStreamerON){
+   fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
+  }
+
 }
   
 //__________________________________________
-void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt
+void AliHFEmcQA::CreateHistograms(const Int_t kquark) 
 {
   // create histograms
 
@@ -266,6 +293,11 @@ void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
     kqTypeLabel[kMisID-4]="miside";
   }
 
+  TString kqEtaRangeLabel[fgkEtaRanges];
+  kqEtaRangeLabel[0] = "mcqa_";
+  kqEtaRangeLabel[1] = "mcqa_barrel_";
+  kqEtaRangeLabel[2] = "mcqa_unitY_";
+
   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
   const Int_t nptbinning1 = 35;
   Int_t iBin[2];
@@ -283,66 +315,144 @@ void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
     xcorrbin[icorrbin]=icorrbin*0.1;
   }
 
+  Int_t fCentrmax = 0;
+  if(!fIsPbPb&&!fIsppMultiBin) fCentrmax=0;
+  else fCentrmax = kCentBins;
   TString hname; 
   if(kquark == kOthers){
-    for (Int_t iqType = 0; iqType < 4; iqType++ ){
-       hname = hnopt+"Pt_"+kqTypeLabel[iqType];
-       //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL
-       fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25);
-       hname = hnopt+"Y_"+kqTypeLabel[iqType];
-       fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
-       hname = hnopt+"Eta_"+kqTypeLabel[iqType];
-       fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
-       // Fill List
-       if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
-    }
-    return;
+   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
+       for (Int_t iqType = 0; iqType < 4; iqType++ ){
+          for(Int_t icentr = 0; icentr < (fCentrmax+1); icentr++)
+          {
+              hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
+              fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),60,0.25,30.25);
+              hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
+              fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
+              hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
+              fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
+              // Fill List
+              if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
+          }
+       }
+   }
+   return;
   }
-  for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
-     if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
-     hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
-     fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
-     hname = hnopt+"Pt_"+kqTypeLabel[iqType];
-     fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[1],kPtbinning1); // new binning
-     //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); // old log binning
-     //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.,30.); // mj to compare with FONLL
-     //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL
-     hname = hnopt+"Y_"+kqTypeLabel[iqType];
-     fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
-     hname = hnopt+"Eta_"+kqTypeLabel[iqType];
-     fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
-     // Fill List
-     if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
+  for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
+   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
+       if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
+        for(Int_t icentr = 0; icentr<(fCentrmax+1); icentr++)
+          {
+              hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
+              fHist[iq][iqType][icut][icentr].fPdgCode = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;PdgCode",hname.Data(),icentr),20001,-10000.5,10000.5);
+              hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
+              fHist[iq][iqType][icut][icentr].fPt = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;p_{T} (GeV/c)",hname.Data(),icentr),iBin[1],kPtbinning1); // new binning
+              hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
+              fHist[iq][iqType][icut][icentr].fY = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;y",hname.Data(),icentr),150,-7.5,7.5);
+              hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
+              fHist[iq][iqType][icut][icentr].fEta = new TH1F(Form("%sCentr_%i",hname.Data(),icentr),Form("%sCentr_%i;eta",hname.Data(),icentr),150,-7.5,7.5);
+              // Fill List
+              if(fQAhistos) fHist[iq][iqType][icut][icentr].FillList(fQAhistos);
+          }
+   }
   }
 
-  if (icut == 0){ 
-    hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
-    fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5);
-    hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
-    fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
+  for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
+    hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); 
+    hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
+    hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
+    hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
+
+    hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
+    hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
+    hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
+    hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
+
+    if(icut <1){
+      hname = kqEtaRangeLabel[icut]+"PtCorrDinein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDineout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDoutein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDouteout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+
+      hname = kqEtaRangeLabel[icut]+"PtCorrDpDinein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDpDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDpDineout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDpDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDpDoutein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDpDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDpDouteout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDpDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+
+      hname = kqEtaRangeLabel[icut]+"PtCorrD0Dinein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrD0Dinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrD0Dineout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrD0Dineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrD0Doutein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrD0Doutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrD0Douteout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrD0Douteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+
+      hname = kqEtaRangeLabel[icut]+"PtCorrDrestDinein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDrestDinein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDrestDineout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDrestDineout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDrestDoutein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDrestDoutein = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrDrestDouteout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrDrestDouteout = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrD_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrDp_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrD0_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrDrest_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrGD_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrGD = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrGDp_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrGDp = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrGD0_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrGD0 = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrGDrest_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrGDrest = new TH2F(hname,hname+";D Y;e eta",200,-10,10,200,-10,10); 
+
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrB_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10); 
+      hname = kqEtaRangeLabel[icut]+"fEtaCorrGB_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fEtaCorrGB = new TH2F(hname,hname+";B Y;e eta",200,-10,10,200,-10,10); 
+
+      hname = kqEtaRangeLabel[icut]+"PtCorrBinein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrBinein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrBineout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrBineout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrBoutein_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrBoutein = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+      hname = kqEtaRangeLabel[icut]+"PtCorrBouteout_"+kqTypeLabel[kQuark];
+      fHistComm[iq][icut].fPtCorrBouteout = new TH2F(hname,hname+";B p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
+    }
+    if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
   }
-  hname = hnopt+"PtCorr_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
-  hname = hnopt+"PtCorrDp_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
-  hname = hnopt+"PtCorrD0_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
-  hname = hnopt+"PtCorrDrest_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
+
+
+  hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
+  fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
+  hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
+  fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
   
-  hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
-  hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
-  hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
-  hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
-  if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
 }
 
 //__________________________________________
@@ -372,6 +482,7 @@ void AliHFEmcQA::Init()
   fParentSelect[1][5] = 5232; //Ksib0
   fParentSelect[1][6] = 5332; //Omegab-
 
+
 }
 
 //__________________________________________
@@ -388,8 +499,11 @@ void AliHFEmcQA::GetMesonKine()
   Int_t id1=0, id2=0;
 
 
-  if(fCentrality>11) printf("warning centrality out of histogram array limits \n");
-
+  if(fCentrality>=11) {
+      AliWarning(Form("Centrality out of histogram array limits: %d", fCentrality));
+      return;
+  }
+  if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 
   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
      if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
@@ -398,13 +512,21 @@ void AliHFEmcQA::GetMesonKine()
      mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
      if(!mctrack0) continue;
 
-     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
+//     if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
 
-     if(abs(mctrack0->PdgCode()) == 111) // pi0 
+     if(TMath::Abs(mctrack0->PdgCode()) == 111) // pi0 
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("pionspectra_centrbin%i",fCentrality),mctrack0->Pt());
             fMCQACollection->Fill(Form("pionspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
+            if(imc>fMCEvent->GetNumberOfPrimaries()) {
+              fMCQACollection->Fill(Form("pionspectrasec_centrbin%i",fCentrality),mctrack0->Pt());
+              fMCQACollection->Fill(Form("pionspectraLogsec_centrbin%i",fCentrality),mctrack0->Pt());
+            }
+            else {
+              fMCQACollection->Fill(Form("pionspectrapri_centrbin%i",fCentrality),mctrack0->Pt());
+              fMCQACollection->Fill(Form("pionspectraLogpri_centrbin%i",fCentrality),mctrack0->Pt());
+            }
           }
           id1=mctrack0->GetFirstDaughter();
           id2=mctrack0->GetLastDaughter();
@@ -412,15 +534,23 @@ void AliHFEmcQA::GetMesonKine()
           for(int idx=id1; idx<=id2; idx++){
             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
-            if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
+            if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
              fMCQACollection->Fill(Form("piondaughters_centrbin%i",fCentrality),mctrackd->Pt());
           }
        }
-     else if(abs(mctrack0->PdgCode()) == 221) // eta 
+     else if(TMath::Abs(mctrack0->PdgCode()) == 221) // eta 
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("etaspectra_centrbin%i",fCentrality),mctrack0->Pt());
             fMCQACollection->Fill(Form("etaspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
+            if(imc>fMCEvent->GetNumberOfPrimaries()) {
+              fMCQACollection->Fill(Form("etaspectrasec_centrbin%i",fCentrality),mctrack0->Pt());
+              fMCQACollection->Fill(Form("etaspectraLogsec_centrbin%i",fCentrality),mctrack0->Pt());
+            }
+            else {
+              fMCQACollection->Fill(Form("etaspectrapri_centrbin%i",fCentrality),mctrack0->Pt());
+              fMCQACollection->Fill(Form("etaspectraLogpri_centrbin%i",fCentrality),mctrack0->Pt());
+            }
           } 
           id1=mctrack0->GetFirstDaughter();
           id2=mctrack0->GetLastDaughter();
@@ -428,11 +558,11 @@ void AliHFEmcQA::GetMesonKine()
           for(int idx=id1; idx<=id2; idx++){
             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
-            if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
+            if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
              fMCQACollection->Fill(Form("etadaughters_centrbin%i",fCentrality),mctrackd->Pt());
           }
        }
-     else if(abs(mctrack0->PdgCode()) == 223) // omega
+     else if(TMath::Abs(mctrack0->PdgCode()) == 223) // omega
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("omegaspectra_centrbin%i",fCentrality),mctrack0->Pt());
@@ -444,11 +574,11 @@ void AliHFEmcQA::GetMesonKine()
           for(int idx=id1; idx<=id2; idx++){
             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
-            if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
+            if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
              fMCQACollection->Fill(Form("omegadaughters_centrbin%i",fCentrality),mctrackd->Pt());
           }
        }
-     else if(abs(mctrack0->PdgCode()) == 333) // phi 
+     else if(TMath::Abs(mctrack0->PdgCode()) == 333) // phi 
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("phispectra_centrbin%i",fCentrality),mctrack0->Pt());
@@ -460,11 +590,11 @@ void AliHFEmcQA::GetMesonKine()
           for(int idx=id1; idx<=id2; idx++){
             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
-            if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
+            if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
              fMCQACollection->Fill(Form("phidaughters_centrbin%i",fCentrality),mctrackd->Pt());
           }
        }
-     else if(abs(mctrack0->PdgCode()) == 331) // eta prime
+     else if(TMath::Abs(mctrack0->PdgCode()) == 331) // eta prime
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("etapspectra_centrbin%i",fCentrality),mctrack0->Pt());
@@ -476,11 +606,11 @@ void AliHFEmcQA::GetMesonKine()
           for(int idx=id1; idx<=id2; idx++){
             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
-            if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
+            if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
              fMCQACollection->Fill(Form("etapdaughters_centrbin%i",fCentrality),mctrackd->Pt());
           }
        }
-     else if(abs(mctrack0->PdgCode()) == 113) // rho
+     else if(TMath::Abs(mctrack0->PdgCode()) == 113) // rho
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("rhospectra_centrbin%i",fCentrality),mctrack0->Pt());
@@ -492,17 +622,17 @@ void AliHFEmcQA::GetMesonKine()
           for(int idx=id1; idx<=id2; idx++){
             if(!(mctrackdaugt = fMCEvent->GetTrack(idx))) continue;
             if(!(mctrackd = dynamic_cast<AliMCParticle *>(mctrackdaugt))) continue;
-            if(abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
+            if(TMath::Abs(mctrackd->PdgCode()) == 11 && TMath::Abs(mctrackd->Eta())<0.8)
              fMCQACollection->Fill(Form("rhodaughters_centrbin%i",fCentrality),mctrackd->Pt());
           }
        }
-     else if(abs(mctrack0->PdgCode()) == 321) // kaon+-
+     else if(TMath::Abs(mctrack0->PdgCode()) == 321) // kaon+-
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("kaonspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
           }
        }
-     else if(abs(mctrack0->PdgCode()) == 130) // k0L
+     else if(TMath::Abs(mctrack0->PdgCode()) == 130) // k0L
        {
           if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8) {
             fMCQACollection->Fill(Form("k0LspectraLog_centrbin%i",fCentrality),mctrack0->Pt());
@@ -527,6 +657,8 @@ void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
       return; 
     }
 
+    if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
+
     AliMCParticle *mctrack = NULL;
     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 
@@ -553,10 +685,10 @@ void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
       // if the heavy particle comes from string which is denoted as particle status 12|12|12...12|11,[PYTHIA p.60]
       // in this case, the mother of heavy particle can be one of the fragmented parton of the string
       // should I make a condition that partMother should be quark or diquark? -> not necessary
-      if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
-      //if ( abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
+      if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11) ){
+      //if ( TMath::Abs(partMother->GetPdgCode()) == kquark || (partMother->GetStatusCode() == 11 || partMother->GetStatusCode() == 12) ){
 
-        if ( abs(partMother->GetPdgCode()) != kquark ){
+        if ( TMath::Abs(partMother->GetPdgCode()) != kquark ){
           // search fragmented partons in the same string
           Bool_t isSameString = kTRUE; 
           for (Int_t i=1; i<fgkMaxIter; i++){
@@ -567,13 +699,13 @@ void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
                AliDebug(1, "Stack label is negative, return\n");
                return; 
              }
-             if ( abs(partMother->GetPdgCode()) == kquark ) break;
+             if ( TMath::Abs(partMother->GetPdgCode()) == kquark ) break;
              if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
              if (!isSameString) return; 
           }
         }
         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
-        if (abs(partMother->GetPdgCode()) != kquark) return; 
+        if (TMath::Abs(partMother->GetPdgCode()) != kquark) return; 
 
         if (fIsHeavy[iq] >= 50) return;  
         fHeavyQuark[fIsHeavy[iq]] = partMother;
@@ -581,15 +713,15 @@ void AliHFEmcQA::GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark)
 
         // fill kinematics for heavy parton
         if (partMother->GetPdgCode() > 0) { // quark
-          fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
-          fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
-          fHist[iq][kQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
-          fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
+          fHist[iq][kQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
+          fHist[iq][kQuark][0][fCentrality].fPt->Fill(partMother->Pt());
+          fHist[iq][kQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
+          fHist[iq][kQuark][0][fCentrality].fEta->Fill(partMother->Eta());
         } else{ // antiquark
-          fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
-          fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
-          fHist[iq][kantiQuark][0].fY->Fill(AliHFEtools::GetRapidity(partMother));
-          fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
+          fHist[iq][kantiQuark][0][fCentrality].fPdgCode->Fill(partMother->GetPdgCode());
+          fHist[iq][kantiQuark][0][fCentrality].fPt->Fill(partMother->Pt());
+          fHist[iq][kantiQuark][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMother));
+          fHist[iq][kantiQuark][0][fCentrality].fEta->Fill(partMother->Eta());
         }
 
       } // end of heavy parton slection loop 
@@ -731,6 +863,8 @@ void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
       return;
     }
 
+    if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
+
     TParticle *partCopy = mcpart;
     Int_t pdgcode = mcpart->GetPdgCode();
     Int_t pdgcodeCopy = pdgcode;
@@ -739,7 +873,7 @@ void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
 
     // if the mother is charmed hadron  
     Bool_t isDirectCharm = kFALSE;
-    if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
+    if ( int(TMath::Abs(pdgcode)/100.) == kCharm || int(TMath::Abs(pdgcode)/1000.) == kCharm ) {
 
           // iterate until you find B hadron as a mother or become top ancester 
           for (Int_t i=1; i<fgkMaxIter; i++){
@@ -759,7 +893,7 @@ void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
              Int_t motherPDG = mother->GetPdgCode();
     
              for (Int_t j=0; j<fNparents; j++){
-                if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
+                if (TMath::Abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
              }
 
              mcpart = mother;
@@ -767,13 +901,13 @@ void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
     } // end of if
     if((isDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
          for (Int_t i=0; i<fNparents; i++){
-            if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
+            if (TMath::Abs(pdgcodeCopy)==fParentSelect[iq][i]){
 
               // fill hadron kinematics
-              fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
-              fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
-              fHist[iq][kHadron][0].fY->Fill(AliHFEtools::GetRapidity(partCopy));
-              fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
+              fHist[iq][kHadron][0][fCentrality].fPdgCode->Fill(pdgcodeCopy);
+              fHist[iq][kHadron][0][fCentrality].fPt->Fill(partCopy->Pt());
+              fHist[iq][kHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partCopy));
+              fHist[iq][kHadron][0][fCentrality].fEta->Fill(partCopy->Eta());
 
               if(iq==0) {
                fhD[i]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
@@ -781,17 +915,17 @@ void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
             }
          }
         // I also want to store D* info to compare with D* measurement 
-        if (abs(pdgcodeCopy)==413 && iq==0) { //D*+
+        if (TMath::Abs(pdgcodeCopy)==413 && iq==0) { //D*+
                fhD[7]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
         }
-        if (abs(pdgcodeCopy)==423 && iq==0) { //D*0
+        if (TMath::Abs(pdgcodeCopy)==423 && iq==0) { //D*0
                fhD[8]->Fill(partCopy->Pt(),AliHFEtools::GetRapidity(partCopy));
         }
     } // end of if
 }
 
 //__________________________________________
-void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut
+void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed) 
 {
     // decay electron kinematics
     
@@ -807,14 +941,29 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
       return;
     }
 
+    if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
+
+    Double_t eabsEta = TMath::Abs(mcpart->Eta());
+    Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
+
     if(kquark==kOthers){
       Int_t esource = -1;
-      if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
+      if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
       if(esource==0|| esource==1 || esource==2 || esource==3){
-        fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
-        fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-        fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
+        fHist[iq][esource][0][fCentrality].fPt->Fill(mcpart->Pt());
+        fHist[iq][esource][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+        fHist[iq][esource][0][fCentrality].fEta->Fill(mcpart->Eta());
+        if(eabsEta<0.9){
+          fHist[iq][esource][1][fCentrality].fPt->Fill(mcpart->Pt());
+          fHist[iq][esource][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+          fHist[iq][esource][1][fCentrality].fEta->Fill(mcpart->Eta());
+        }
+        if(eabsY<0.5){
+          fHist[iq][esource][2][fCentrality].fPt->Fill(mcpart->Pt());
+          fHist[iq][esource][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+          fHist[iq][esource][2][fCentrality].fEta->Fill(mcpart->Eta());
+        }
         return; 
       }
       else {
@@ -823,7 +972,7 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
       }
     }
 
-    if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
+    if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
 
     Int_t iLabel = mcpart->GetFirstMother(); 
     if (iLabel<0){
@@ -854,10 +1003,10 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
 
     // if the mother is charmed hadron  
     Bool_t isMotherDirectCharm = kFALSE;
-    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { 
+    if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) { 
 
          for (Int_t i=0; i<fNparents; i++){
-            if (abs(maPdgcode)==fParentSelect[0][i]){
+            if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
               isFinalOpenCharm = kTRUE;
             } 
          }  
@@ -882,26 +1031,95 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
              Int_t grandMaPDG = grandMa->GetPdgCode();
 
              for (Int_t j=0; j<fNparents; j++){
-                if (abs(grandMaPDG)==fParentSelect[1][j]){
+                if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
 
                   if (kquark == kCharm) return;
                   // fill electron kinematics
-                  fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-                  fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
-                  fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-                  fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
+                  fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+                  fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
+                  fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                  fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
 
                   // fill mother hadron kinematics
-                  fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); 
-                  fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
-                  fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
-                  fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
-                 
+                  fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG); 
+                  fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
+                  fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+                  fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
+
+                  if(eabsEta<0.9){
+                    fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+                    fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
+                    fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                    fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
+
+                    // fill mother hadron kinematics
+                    fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG); 
+                    fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
+                    fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+                    fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
+                  }
+
+                  if(eabsY<0.5){
+                    fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+                    fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
+                    fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                    fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
+
+                    // fill mother hadron kinematics
+                    fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG); 
+                    fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
+                    fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+                    fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
+                  }
+
+                  //mj: to calculate B to e eta correlation to calculate total heavy quark cross section
+                  Int_t kLabel0 = grandMa->GetFirstMother();
+                  Bool_t isGGrandmaYes = kFALSE;
+                  Double_t ggmrapidwstmp=0;
+                  if (!(kLabel0 < 0)){ // safety protection
+                    if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel0))))){
+                      TParticle* ggrandMatmp = mctrack->Particle();
+                      Int_t ggrandMaPDGtmp = ggrandMatmp->GetPdgCode();
+                      if ( int(TMath::Abs(ggrandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(ggrandMaPDGtmp)/1000.) == kBeauty) isGGrandmaYes = kTRUE;
+                      ggmrapidwstmp = AliHFEtools::GetRapidity(ggrandMatmp);
+                    }
+                  }
+
+                  Double_t gmrapidwstmp0 = AliHFEtools::GetRapidity(grandMa);
+                  Double_t eetawstmp0 = mcpart->Eta();
+  
+                  Double_t gmrapidtmp0 = TMath::Abs(gmrapidwstmp0);
+                  Double_t eetatmp0 = TMath::Abs(eetawstmp0);
+
+                  fHistComm[iq][0].fEtaCorrB->Fill(gmrapidwstmp0,eetawstmp0);
+                  if(isGGrandmaYes) fHistComm[iq][0].fEtaCorrGB->Fill(ggmrapidwstmp,eetawstmp0);
+                  else fHistComm[iq][0].fEtaCorrGB->Fill(gmrapidwstmp0,eetawstmp0);
+
+                  if(gmrapidtmp0<0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBinein->Fill(grandMa->Pt(),mcpart->Pt());
+                  else if(gmrapidtmp0<0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBineout->Fill(grandMa->Pt(),mcpart->Pt());
+                  else if(gmrapidtmp0>0.5 && eetatmp0<0.5 ) fHistComm[iq][0].fPtCorrBoutein->Fill(grandMa->Pt(),mcpart->Pt());
+                  else if(gmrapidtmp0>0.5 && eetatmp0>0.5 ) fHistComm[iq][0].fPtCorrBouteout->Fill(grandMa->Pt(),mcpart->Pt());
+                  //======================================================================================
+
                   // ratio between pT of electron and pT of mother B hadron 
-                  if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                  if(grandMa->Pt()) {
+                    fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                    if(eabsEta<0.9){
+                      fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                    }
+                    if(eabsY<0.5){
+                      fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                    }
+                  }
 
                   // distance between electron production point and primary vertex
-                  fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  if(eabsEta<0.9){
+                    fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  }
+                  if(eabsY<0.5){
+                    fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  }
                   return;
                 }
              } 
@@ -911,53 +1129,161 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
     } // end of if
     if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
          for (Int_t i=0; i<fNparents; i++){
-            if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
+            if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
+
+              fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+              fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
+              fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+              fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());  
 
-//mj weighting to consider measured spectra!!!
-              Double_t mpt=partMotherCopy->Pt();
-              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
-              //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));
-              // fill electron kinematics
-              if(iq==0){
-                fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode(),wfactor);
-                fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt(),wfactor);
-                fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart),wfactor);
-                fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta(),wfactor);  
-
-                fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-                fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
-                fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-                fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());  
-              } 
-              else{
-                fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-                fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
-                fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-                fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());  
-              }
-//--------------
               // fill mother hadron kinematics
-              fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); 
-              fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
-              fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
-              fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
+              fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
+              fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
+              fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+              fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
+
+              if(eabsEta<0.9){
+                fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+                fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
+                fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());  
+
+                // fill mother hadron kinematics
+                fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
+                fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
+                fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+                fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
+              }
+
+              if(eabsY<0.5){
+                fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+                fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
+                fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());  
+
+                // fill mother hadron kinematics
+                fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy); 
+                fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
+                fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+                fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
+              }
 
               // ratio between pT of electron and pT of mother B or direct D hadron 
-              if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
-              fHistComm[iq][icut].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
-              if(TMath::Abs(partMotherCopy->GetPdgCode())==411) fHistComm[iq][icut].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
-              else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) fHistComm[iq][icut].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
-              else fHistComm[iq][icut].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              if(partMotherCopy->Pt()) {
+                 fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
+                 if(eabsEta<0.9){
+                   fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
+                 }
+                 if(eabsY<0.5){
+                   fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
+                 }
+              }
+              fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              if(eabsEta<0.9){
+                fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              }
+              if(eabsY<0.5){
+                fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              }
+              if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
+                fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                if(eabsEta<0.9){
+                  fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+                if(eabsY<0.5){
+                  fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+              }
+              else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
+                fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                if(eabsEta<0.9){
+                  fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+                if(eabsY<0.5){
+                  fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+              }
+              else {
+                fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                if(eabsEta<0.9){
+                  fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+                if(eabsY<0.5){
+                  fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+              }
+
+              //mj: to calculate D to e eta correlation to calculate total heavy quark cross section
+              Int_t kLabel = partMotherCopy->GetFirstMother();
+              Bool_t isGrandmaYes = kFALSE;
+              Double_t gmrapidwstmp =0;
+              if (!(kLabel < 0)){ // safety protection
+                if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(kLabel))))){ 
+                  TParticle* grandMatmp = mctrack->Particle();
+                  Int_t grandMaPDGtmp = grandMatmp->GetPdgCode();
+                  if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kCharm || int(TMath::Abs(grandMaPDGtmp)/1000.) == kCharm ) isGrandmaYes = kTRUE;
+                  if ( int(TMath::Abs(grandMaPDGtmp)/100.) == kBeauty || int(TMath::Abs(grandMaPDGtmp)/1000.) == kBeauty) isGrandmaYes = kTRUE;
+                  gmrapidwstmp = AliHFEtools::GetRapidity(grandMatmp);
+                }
+              }
+
+              Double_t mrapidwstmp = AliHFEtools::GetRapidity(partMotherCopy);
+              Double_t eetawstmp = mcpart->Eta();
+
+              Double_t mrapidtmp = TMath::Abs(mrapidwstmp);
+              Double_t eetatmp = TMath::Abs(eetawstmp);
+
+              fHistComm[iq][0].fEtaCorrD->Fill(mrapidwstmp,eetawstmp);
+              if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD->Fill(gmrapidwstmp,eetawstmp);
+              else fHistComm[iq][0].fEtaCorrGD->Fill(mrapidwstmp,eetawstmp);
+
+              if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
+                fHistComm[iq][0].fEtaCorrDp->Fill(mrapidwstmp,eetawstmp);
+                if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDp->Fill(gmrapidwstmp,eetawstmp);
+                else fHistComm[iq][0].fEtaCorrGDp->Fill(mrapidwstmp,eetawstmp);
+                if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDpDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDpDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              }
+              else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
+                fHistComm[iq][0].fEtaCorrD0->Fill(mrapidwstmp,eetawstmp);
+                if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGD0->Fill(gmrapidwstmp,eetawstmp);
+                else fHistComm[iq][0].fEtaCorrGD0->Fill(mrapidwstmp,eetawstmp);
+                if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Dinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Dineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrD0Doutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrD0Douteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              }
+              else {
+                fHistComm[iq][0].fEtaCorrDrest->Fill(mrapidwstmp,eetawstmp);
+                if(isGrandmaYes) fHistComm[iq][0].fEtaCorrGDrest->Fill(gmrapidwstmp,eetawstmp);
+                else fHistComm[iq][0].fEtaCorrGDrest->Fill(mrapidwstmp,eetawstmp);
+                if(mrapidtmp<0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDinein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp<0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDineout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp>0.5 && eetatmp<0.5 ) fHistComm[iq][0].fPtCorrDrestDoutein->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                else if(mrapidtmp>0.5 && eetatmp>0.5 ) fHistComm[iq][0].fPtCorrDrestDouteout->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              }
 
               // distance between electron production point and primary vertex
-              fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              if(eabsEta<0.9){
+                fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              }
+              if(eabsY<0.5){
+                fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              }
             }
          }
     } // end of if
 }
 
 //____________________________________________________________________
-void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
+void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
 {
   // decay electron kinematics
 
@@ -974,7 +1300,10 @@ void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I
     return;
   }
 
-  if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
+  if ( TMath::Abs(mcpart->GetPdgCode()) != kdecayed ) return;
+
+  Double_t eabsEta = TMath::Abs(mcpart->Eta());
+  Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
 
   // mother
   Int_t iLabel = mcpart->GetMother();
@@ -983,16 +1312,18 @@ void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I
     return;
   }
 
+  if(!fIsPbPb&&!fIsppMultiBin) fCentrality=0;
+
   AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
   AliAODMCParticle *partMotherCopy = partMother;
   Int_t maPdgcode = partMother->GetPdgCode();
   Int_t maPdgcodeCopy = maPdgcode;
 
   Bool_t isMotherDirectCharm = kFALSE;
-  if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
+  if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
 
     for (Int_t i=0; i<fNparents; i++){
-       if (abs(maPdgcode)==fParentSelect[0][i]){
+       if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
          isFinalOpenCharm = kTRUE;
        }
     } 
@@ -1015,20 +1346,47 @@ void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I
        Int_t grandMaPDG = grandMa->GetPdgCode();
 
        for (Int_t j=0; j<fNparents; j++){
-          if (abs(grandMaPDG)==fParentSelect[1][j]){
+          if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
 
             if (kquark == kCharm) return;
             // fill electron kinematics
-            fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-            fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
-            fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-            fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
+            fHist[iq][kElectron2nd][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+            fHist[iq][kElectron2nd][0][fCentrality].fPt->Fill(mcpart->Pt());
+            fHist[iq][kElectron2nd][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+            fHist[iq][kElectron2nd][0][fCentrality].fEta->Fill(mcpart->Eta());
 
             // fill mother hadron kinematics
-            fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
-            fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
-            fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
-            fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
+            fHist[iq][kDeHadron][0][fCentrality].fPdgCode->Fill(grandMaPDG);
+            fHist[iq][kDeHadron][0][fCentrality].fPt->Fill(grandMa->Pt());
+            fHist[iq][kDeHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+            fHist[iq][kDeHadron][0][fCentrality].fEta->Fill(grandMa->Eta());
+
+            if(eabsEta<0.9){
+              // fill electron kinematics
+              fHist[iq][kElectron2nd][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+              fHist[iq][kElectron2nd][1][fCentrality].fPt->Fill(mcpart->Pt());
+              fHist[iq][kElectron2nd][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+              fHist[iq][kElectron2nd][1][fCentrality].fEta->Fill(mcpart->Eta());
+
+              // fill mother hadron kinematics
+              fHist[iq][kDeHadron][1][fCentrality].fPdgCode->Fill(grandMaPDG);
+              fHist[iq][kDeHadron][1][fCentrality].fPt->Fill(grandMa->Pt());
+              fHist[iq][kDeHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+              fHist[iq][kDeHadron][1][fCentrality].fEta->Fill(grandMa->Eta());
+            }
+            if(eabsY<0.5){
+              // fill electron kinematics
+              fHist[iq][kElectron2nd][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+              fHist[iq][kElectron2nd][2][fCentrality].fPt->Fill(mcpart->Pt());
+              fHist[iq][kElectron2nd][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+              fHist[iq][kElectron2nd][2][fCentrality].fEta->Fill(mcpart->Eta());
+
+              // fill mother hadron kinematics
+              fHist[iq][kDeHadron][2][fCentrality].fPdgCode->Fill(grandMaPDG);
+              fHist[iq][kDeHadron][2][fCentrality].fPt->Fill(grandMa->Pt());
+              fHist[iq][kDeHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+              fHist[iq][kDeHadron][2][fCentrality].fEta->Fill(grandMa->Eta());
+            }
 
             return;
           }
@@ -1039,19 +1397,46 @@ void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I
   } // end of if
   if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
     for (Int_t i=0; i<fNparents; i++){
-       if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
+       if (TMath::Abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 
          // fill electron kinematics
-         fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-         fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
-         fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-         fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
+         fHist[iq][kElectron][0][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+         fHist[iq][kElectron][0][fCentrality].fPt->Fill(mcpart->Pt());
+         fHist[iq][kElectron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+         fHist[iq][kElectron][0][fCentrality].fEta->Fill(mcpart->Eta());
 
          // fill mother hadron kinematics
-         fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
-         fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
-         fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
-         fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
+         fHist[iq][keHadron][0][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
+         fHist[iq][keHadron][0][fCentrality].fPt->Fill(partMotherCopy->Pt());
+         fHist[iq][keHadron][0][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+         fHist[iq][keHadron][0][fCentrality].fEta->Fill(partMotherCopy->Eta());
+
+         if(eabsEta<0.9){
+           // fill electron kinematics
+           fHist[iq][kElectron][1][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+           fHist[iq][kElectron][1][fCentrality].fPt->Fill(mcpart->Pt());
+           fHist[iq][kElectron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+           fHist[iq][kElectron][1][fCentrality].fEta->Fill(mcpart->Eta());
+
+           // fill mother hadron kinematics
+           fHist[iq][keHadron][1][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
+           fHist[iq][keHadron][1][fCentrality].fPt->Fill(partMotherCopy->Pt());
+           fHist[iq][keHadron][1][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+           fHist[iq][keHadron][1][fCentrality].fEta->Fill(partMotherCopy->Eta());
+         }
+         if(eabsY<0.5){
+           // fill electron kinematics
+           fHist[iq][kElectron][2][fCentrality].fPdgCode->Fill(mcpart->GetPdgCode());
+           fHist[iq][kElectron][2][fCentrality].fPt->Fill(mcpart->Pt());
+           fHist[iq][kElectron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+           fHist[iq][kElectron][2][fCentrality].fEta->Fill(mcpart->Eta());
+
+           // fill mother hadron kinematics
+           fHist[iq][keHadron][2][fCentrality].fPdgCode->Fill(maPdgcodeCopy);
+           fHist[iq][keHadron][2][fCentrality].fPt->Fill(partMotherCopy->Pt());
+           fHist[iq][keHadron][2][fCentrality].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+           fHist[iq][keHadron][2][fCentrality].fEta->Fill(partMotherCopy->Eta());
+         }
 
        }
     }
@@ -1090,17 +1475,17 @@ void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &moth
            
        motherlabel = -1;
 
-       if (abs(afterinitialrad1->GetPdgCode()) == fgkGluon && abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
+       if (TMath::Abs(afterinitialrad1->GetPdgCode()) == fgkGluon && TMath::Abs(afterinitialrad2->GetPdgCode()) == fgkGluon){
          AliDebug(1,"heavy from gluon gluon pair creation!\n");
          mothertype = -1;
          motherID = fgkGluon;
        }
-       else if (abs(afterinitialrad1->GetPdgCode()) == kquark || abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
+       else if (TMath::Abs(afterinitialrad1->GetPdgCode()) == kquark || TMath::Abs(afterinitialrad2->GetPdgCode()) == kquark){ // one from Q and the other from g
          AliDebug(1,"heavy from flavor exitation!\n");
          mothertype = -1;
          motherID = kquark;
        }
-       else if  (abs(afterinitialrad1->GetPdgCode()) == abs(afterinitialrad2->GetPdgCode())){
+       else if  (TMath::Abs(afterinitialrad1->GetPdgCode()) == TMath::Abs(afterinitialrad2->GetPdgCode())){
          AliDebug(1,"heavy from q-qbar pair creation!\n");
          mothertype = -1;
          motherID = 1;
@@ -1164,11 +1549,11 @@ void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &mo
 }
 
 //__________________________________________
-Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
+Int_t AliHFEmcQA::GetSource(const AliVParticle* const mcpart) const
 {        
   // decay particle's origin 
 
-  //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
+  //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
        
   Int_t origin = -1;
   Bool_t isFinalOpenCharm = kFALSE;
@@ -1179,20 +1564,21 @@ Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
   }
 
   // mother
-  Int_t iLabel = mcpart->GetMother();
+  // Information not in the base class, cast necessary
+  Int_t iLabel = GetMother(mcpart);
   if (iLabel<0){
     AliDebug(1, "Stack label is negative, return\n");
     return -1;
   } 
        
-  AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel);
-  Int_t maPdgcode = partMother->GetPdgCode();
+  const AliVParticle *partMother = fMCEvent->GetTrack(iLabel);
+  Int_t maPdgcode = partMother->PdgCode();
   
   // if the mother is charmed hadron  
-  if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
+  if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
     
     for (Int_t i=0; i<fNparents; i++){
-       if (abs(maPdgcode)==fParentSelect[0][i]){
+       if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
          isFinalOpenCharm = kTRUE;
        }
     }
@@ -1200,8 +1586,8 @@ Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
 
     // iterate until you find B hadron as a mother or become top ancester 
     for (Int_t i=1; i<fgkMaxIter; i++){
-
-       Int_t jLabel = partMother->GetMother();
+      
+       Int_t jLabel = GetMother(partMother);
        if (jLabel == -1){
          origin = kDirectCharm;
          return origin;
@@ -1212,11 +1598,11 @@ Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
        }
 
        // if there is an ancester
-       AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel);
-       Int_t grandMaPDG = grandMa->GetPdgCode();
+       const AliVParticle* grandMa = fMCEvent->GetTrack(jLabel);
+       Int_t grandMaPDG = grandMa->PdgCode();
 
        for (Int_t j=0; j<fNparents; j++){
-          if (abs(grandMaPDG)==fParentSelect[1][j]){
+          if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
             origin = kBeautyCharm;
             return origin;
           }
@@ -1225,19 +1611,19 @@ Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
        partMother = grandMa;
     } // end of iteration 
   } // end of if
-  else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
+  else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
     for (Int_t i=0; i<fNparents; i++){
-       if (abs(maPdgcode)==fParentSelect[1][i]){
+       if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
          origin = kDirectBeauty;
          return origin;
        }
     }
   } // end of if
-  else if ( abs(maPdgcode) == 22 ) {
+  else if ( TMath::Abs(maPdgcode) == 22 ) {
     origin = kGamma;
     return origin;
   } // end of if
-  else if ( abs(maPdgcode) == 111 ) {
+  else if ( TMath::Abs(maPdgcode) == 111 ) {
     origin = kPi0;
     return origin;
   } // end of if
@@ -1246,11 +1632,11 @@ Int_t AliHFEmcQA::GetSource(const AliAODMCParticle * const mcpart)
 }
 
 //__________________________________________
-Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
+Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart) const
 {
   // decay particle's origin 
 
-  //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
+  //if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
 
   Int_t origin = -1;
   Bool_t isFinalOpenCharm = kFALSE;
@@ -1272,10 +1658,10 @@ Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
   Int_t maPdgcode = partMother->GetPdgCode();
 
    // if the mother is charmed hadron  
-   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
+   if ( int(TMath::Abs(maPdgcode)/100.) == kCharm || int(TMath::Abs(maPdgcode)/1000.) == kCharm ) {
 
      for (Int_t i=0; i<fNparents; i++){
-        if (abs(maPdgcode)==fParentSelect[0][i]){
+        if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
           isFinalOpenCharm = kTRUE;
         }
      }
@@ -1300,7 +1686,7 @@ Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
         Int_t grandMaPDG = grandMa->GetPdgCode();
 
         for (Int_t j=0; j<fNparents; j++){
-           if (abs(grandMaPDG)==fParentSelect[1][j]){
+           if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
              origin = kBeautyCharm;
              return origin;
            }
@@ -1309,19 +1695,19 @@ Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
         partMother = grandMa;
      } // end of iteration 
    } // end of if
-   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
+   else if ( int(TMath::Abs(maPdgcode)/100.) == kBeauty || int(TMath::Abs(maPdgcode)/1000.) == kBeauty ) {
      for (Int_t i=0; i<fNparents; i++){
-        if (abs(maPdgcode)==fParentSelect[1][i]){
+        if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
           origin = kDirectBeauty;
           return origin;
         }
      }
    } // end of if
-   else if ( abs(maPdgcode) == 22 ) {
+   else if ( TMath::Abs(maPdgcode) == 22 ) {
      origin = kGamma;
      return origin;
    } // end of if
-   else if ( abs(maPdgcode) == 111 ) {
+   else if ( TMath::Abs(maPdgcode) == 111 ) {
      origin = kPi0;
      return origin;
    } // end of if
@@ -1331,7 +1717,7 @@ Int_t AliHFEmcQA::GetSource(const TParticle * const mcpart)
 }
 
 //__________________________________________
-Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
+Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart) const
 {
   // decay particle's origin 
 
@@ -1340,7 +1726,7 @@ Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
     return -1;
   }
 
-  if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
+  if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
 
   Int_t origin = -1;
   Bool_t isFinalOpenCharm = kFALSE;
@@ -1357,12 +1743,24 @@ Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
   TParticle *partMother = mctrack->Particle();
   TParticle *partMotherCopy = mctrack->Particle();
   Int_t maPdgcode = partMother->GetPdgCode();
+  Int_t grmaPdgcode;
+  Int_t ggrmaPdgcode;
 
    // if the mother is charmed hadron  
-   if ( (int(abs(maPdgcode)/100.)%10) == kCharm || (int(abs(maPdgcode)/1000.)%10) == kCharm ) {
+
+   if(TMath::Abs(maPdgcode)==443){ // J/spi
+      Int_t jLabel = partMother->GetFirstMother();
+      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))){
+        TParticle* grandMa = mctrack->Particle();
+        Int_t grandMaPDG = grandMa->GetPdgCode();
+        if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) return kB2Jpsi;
+      }
+      return kJpsi;   
+   } 
+   else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
 
      for (Int_t i=0; i<fNparents; i++){
-        if (abs(maPdgcode)==fParentSelect[0][i]){
+        if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
           isFinalOpenCharm = kTRUE;
         }
      }
@@ -1373,8 +1771,7 @@ Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
 
         Int_t jLabel = partMother->GetFirstMother();
         if (jLabel == -1){
-          origin = kDirectCharm;
-          return origin;
+          return kDirectCharm;
         }
         if (jLabel < 0){ // safety protection
           AliDebug(1, "Stack label is negative, return\n");
@@ -1387,87 +1784,385 @@ Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
         Int_t grandMaPDG = grandMa->GetPdgCode();
 
         for (Int_t j=0; j<fNparents; j++){
-           if (abs(grandMaPDG)==fParentSelect[1][j]){
-             origin = kBeautyCharm;
-             return origin;
+           if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
+             return kBeautyCharm;
            }
         }
 
         partMother = grandMa;
      } // end of iteration 
    } // end of if
-   else if ( (int(abs(maPdgcode)/100.)%10) == kBeauty || (int(abs(maPdgcode)/1000.)%10) == kBeauty ) {
+   else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
      for (Int_t i=0; i<fNparents; i++){
-        if (abs(maPdgcode)==fParentSelect[1][i]){
-          origin = kDirectBeauty;
-          return origin;
+        if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
+          return kDirectBeauty;
         }
      }
    } // end of if
-   else if ( abs(maPdgcode) == 22 ) { //conversion
+   else if ( TMath::Abs(maPdgcode) == 22 ) { //conversion
 
      tmpMomLabel = partMotherCopy->GetFirstMother();
      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) return -1;
      partMother = mctrack->Particle();
      maPdgcode = partMother->GetPdgCode();
-     if ( abs(maPdgcode) == 111 ) {
-       origin = kGammaPi0;
-       return origin;
+
+     // check if the ligth meson is the decay product of heavy mesons
+     tmpMomLabel = partMother->GetFirstMother();
+     if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
+      partMother = mctrack->Particle();
+      grmaPdgcode = partMother->GetPdgCode();
+
+      if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
+      if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
+      if ( TMath::Abs(grmaPdgcode) == 221 ) return kGammaEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
+      //if ( TMath::Abs(grmaPdgcode) == 221 ) return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
+
+      tmpMomLabel = partMother->GetFirstMother();
+      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
+       partMother = mctrack->Particle();
+       ggrmaPdgcode = partMother->GetPdgCode();
+
+       if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kGammaB2M;
+       if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kGammaD2M;
+      }
+     }
+
+     if ( TMath::Abs(maPdgcode) == 111 ) {
+       return kGammaPi0;
      } 
-     else if ( abs(maPdgcode) == 221 ) {
-       origin = kGammaEta;
-       return origin;
+     else if ( TMath::Abs(maPdgcode) == 221 ) {
+       return kGammaEta;
      } 
-     else if ( abs(maPdgcode) == 223 ) {
-       origin = kGammaOmega;
-       return origin;
+     else if ( TMath::Abs(maPdgcode) == 223 ) {
+       return kGammaOmega;
      } 
-     else if ( abs(maPdgcode) == 333 ) {
-       origin = kGammaPhi;
-       return origin;
+     else if ( TMath::Abs(maPdgcode) == 333 ) {
+       return kGammaPhi;
      }
-     else if ( abs(maPdgcode) == 331 ) {
-       origin = kGammaEtaPrime;
-       return origin; 
+     else if ( TMath::Abs(maPdgcode) == 331 ) {
+       return kGammaEtaPrime; 
      }
-     else if ( abs(maPdgcode) == 113 ) {
-       origin = kGammaRho0;
-       return origin;
+     else if ( TMath::Abs(maPdgcode) == 113 ) {
+       return kGammaRho0;
      }
      else origin = kElse;
-     //origin = kGamma; // finer category above
      return origin;
 
-   } // end of if
-   else if ( abs(maPdgcode) == 111 ) {
-     origin = kPi0;
-     return origin;
-   } // end of if
-   else if ( abs(maPdgcode) == 221 ) {
-     origin = kEta;
-     return origin;
-   } // end of if
-   else if ( abs(maPdgcode) == 223 ) {
-     origin = kOmega;
-     return origin;
-   } // end of if
-   else if ( abs(maPdgcode) == 333 ) {
-     origin = kPhi;
-     return origin;
-   } // end of if
-   else if ( abs(maPdgcode) == 331 ) {
-     origin = kEtaPrime;
-     return origin;
-   } // end of if
-   else if ( abs(maPdgcode) == 113 ) {
-     origin = kRho0;
-     return origin;
-   } // end of if
-   else{ 
-    origin = kElse;
+   } 
+   else {
+
+     // check if the ligth meson is the decay product of heavy mesons
+     tmpMomLabel = partMotherCopy->GetFirstMother();
+     if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
+      partMother = mctrack->Particle();
+      grmaPdgcode = partMother->GetPdgCode();
+
+      if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) return kB2M;
+      if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) return kD2M;
+      if ( TMath::Abs(grmaPdgcode) == 221 ) return kEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
+
+      tmpMomLabel = partMother->GetFirstMother();
+      if((mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tmpMomLabel))))) {
+       partMother = mctrack->Particle();
+       ggrmaPdgcode = partMother->GetPdgCode();
+
+       if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) return kB2M;
+       if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) return kD2M;
+      }
+     }
+
+     if ( TMath::Abs(maPdgcode) == 111 ) {
+       return kPi0;
+     } 
+     else if ( TMath::Abs(maPdgcode) == 221 ) {
+       return kEta;
+     } 
+     else if ( TMath::Abs(maPdgcode) == 223 ) {
+       return kOmega;
+     } 
+     else if ( TMath::Abs(maPdgcode) == 333 ) {
+       return kPhi;
+     } 
+     else if ( TMath::Abs(maPdgcode) == 331 ) {
+       return kEtaPrime;
+     } 
+     else if ( TMath::Abs(maPdgcode) == 113 ) {
+       return kRho0;
+     } 
+     else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
+       return kKe3;
+     }
+     else{ 
+      origin = kElse;
+     }
+
    }
    return origin;
 }
+
+//__________________________________________
+Int_t AliHFEmcQA::GetElecSource(const AliVParticle * const mctrack) const
+{
+  //
+  // decay particle's origin 
+  //
+
+  if(!mctrack){
+    AliDebug(1, "no mcparticle, return\n");
+    return -1;
+  }
+
+  TClass *type = mctrack->IsA();  
+
+  if(type == AliMCParticle::Class()) {
+    const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
+    if(esdmc){
+      TParticle *mcpart =  esdmc->Particle();
+      return GetElecSource(mcpart);
+    }
+    else return -1;
+  }
+  if(type == AliAODMCParticle::Class()) {
+    const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(mctrack);
+    if(aodmc){
+      return GetElecSource(aodmc);
+    }
+    else return -1;
+  }
+  return -1;
+}
+//__________________________________________
+Int_t AliHFEmcQA::GetElecSource(const AliAODMCParticle * const mcpart) const
+{
+  //
+  // Function for AliAODMCParticle
+  //
+
+  if (!mcpart) return -1;
+  if (!fMCArray) return -1;
+  
+  if ( TMath::Abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return kMisID;
+
+  Int_t origin = -1;
+  Bool_t isFinalOpenCharm = kFALSE;
+
+  Int_t iLabel = mcpart->GetMother();
+  if ((iLabel<0) || (iLabel>=fMCArray->GetEntriesFast())){
+    AliDebug(1, "label is out of range, return\n");
+    return -1;
+  }
+
+  AliAODMCParticle *mctrack = NULL; // will change all the time
+  Int_t tmpMomLabel=0;
+  if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return -1; 
+  AliAODMCParticle *partMother = mctrack; 
+  AliAODMCParticle *partMotherCopy = mctrack;
+  Int_t maPdgcode = mctrack->GetPdgCode();
+  Int_t grmaPdgcode;
+  Int_t ggrmaPdgcode;
+
+  // if the mother is charmed hadron  
+
+  if(TMath::Abs(maPdgcode)==443){ 
+    //
+    // J/spi
+    //
+    Int_t jLabel = partMother->GetMother();
+    if ((jLabel>=0) && (jLabel<fMCArray->GetEntriesFast())){
+      if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))){
+       Int_t grandMaPDG = mctrack->GetPdgCode();
+       if((int(TMath::Abs(grandMaPDG)/100.)%10) == kBeauty || (int(TMath::Abs(grandMaPDG)/1000.)%10) == kBeauty) {
+         return kB2Jpsi;
+       }
+      }
+    }
+    return kJpsi;   
+  } 
+  else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(maPdgcode)/1000.)%10) == kCharm ) {
+    //
+    // charm
+    //
+    for (Int_t i=0; i<fNparents; i++){
+      if (TMath::Abs(maPdgcode)==fParentSelect[0][i]){
+       isFinalOpenCharm = kTRUE;
+      }
+    }
+    if (!isFinalOpenCharm) {
+      return -1;
+    }
+    
+    // iterate until you find B hadron as a mother or become top ancester 
+    for (Int_t i=1; i<fgkMaxIter; i++){
+      
+      Int_t jLabel = partMother->GetMother();
+      if (jLabel == -1){
+       return kDirectCharm;
+      }
+      if ((jLabel<0) || (jLabel>=fMCArray->GetEntriesFast())){
+       AliDebug(1, "Stack label is negative, return\n");
+       return -1;
+      }
+      
+      // if there is an ancester
+      if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))) {
+       return -1; 
+      }
+      Int_t grandMaPDG = mctrack->GetPdgCode();
+      for (Int_t j=0; j<fNparents; j++){
+       if (TMath::Abs(grandMaPDG)==fParentSelect[1][j]){
+         return kBeautyCharm;
+       }
+      }
+      partMother = mctrack;
+    } // end of iteration 
+
+  } // end of if
+  else if ( (int(TMath::Abs(maPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(maPdgcode)/1000.)%10) == kBeauty ) {
+    //
+    // beauty
+    //
+    for (Int_t i=0; i<fNparents; i++){
+      if (TMath::Abs(maPdgcode)==fParentSelect[1][i]){
+       return kDirectBeauty;
+      }
+    }
+  } // end of if
+  else if ( TMath::Abs(maPdgcode) == 22 ) { 
+    //
+    //conversion
+    //
+    tmpMomLabel = partMotherCopy->GetMother();
+    if((tmpMomLabel<0) || (tmpMomLabel>=fMCArray->GetEntriesFast())) {
+      return -1;
+    }
+    if(!(mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
+      return -1;
+    }
+    partMother = mctrack;
+    maPdgcode = partMother->GetPdgCode();
+    
+    // check if the ligth meson is the decay product of heavy mesons
+    tmpMomLabel = partMother->GetMother();
+    if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
+      if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
+       partMother = mctrack;
+       grmaPdgcode = partMother->GetPdgCode();
+       
+       if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {         
+         return kGammaB2M;
+       }
+       if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {   
+         return kGammaD2M;
+       }
+       if ( TMath::Abs(grmaPdgcode) == 221 ) {  
+         return kGammaEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
+         //return kElse; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
+       }
+       
+       tmpMomLabel = partMother->GetMother();
+       if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
+         if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
+           partMother = mctrack;
+           ggrmaPdgcode = partMother->GetPdgCode();
+           
+           if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
+             return kGammaB2M;
+           }
+           if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
+             return kGammaD2M;
+           }
+          }
+       }
+      }
+    }
+
+    if ( TMath::Abs(maPdgcode) == 111 ) {
+      return kGammaPi0;
+    } 
+    else if ( TMath::Abs(maPdgcode) == 221 ) {
+      return kGammaEta;
+    } 
+    else if ( TMath::Abs(maPdgcode) == 223 ) {
+      return kGammaOmega;
+    } 
+    else if ( TMath::Abs(maPdgcode) == 333 ) {
+      return kGammaPhi;
+    }
+    else if ( TMath::Abs(maPdgcode) == 331 ) {
+      return kGammaEtaPrime; 
+    }
+    else if ( TMath::Abs(maPdgcode) == 113 ) {
+      return kGammaRho0;
+    }
+    else origin = kElse;
+   
+    return origin;
+    
+  } 
+   else {
+     //
+     // check if the ligth meson is the decay product of heavy mesons
+     //
+     tmpMomLabel = partMotherCopy->GetMother();
+     if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
+       if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
+        partMother = mctrack;
+        grmaPdgcode = partMother->GetPdgCode();
+        
+        if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kBeauty ) {
+          return kB2M;
+        }
+        if ( (int(TMath::Abs(grmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(grmaPdgcode)/1000.)%10) == kCharm ) {
+          return kD2M;
+        }
+        if ( TMath::Abs(grmaPdgcode) == 221 ) {
+          return kEta2Pi0; //mj to remove secondary pi0 decay electrons from eta decays, mainly to eta signal enhance samples
+        }
+        
+        tmpMomLabel = partMother->GetMother();
+        if((tmpMomLabel>=0) && (tmpMomLabel<fMCArray->GetEntriesFast())) {
+          if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(tmpMomLabel))))) {
+            partMother = mctrack;
+            ggrmaPdgcode = partMother->GetPdgCode();
+            
+            if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kBeauty || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kBeauty ) {
+              return kB2M;
+            }
+            if ( (int(TMath::Abs(ggrmaPdgcode)/100.)%10) == kCharm || (int(TMath::Abs(ggrmaPdgcode)/1000.)%10) == kCharm ) {
+              return kD2M;
+            }
+          }
+        }
+       }
+       
+       if ( TMath::Abs(maPdgcode) == 111 ) {
+        return kPi0;
+       } 
+       else if ( TMath::Abs(maPdgcode) == 221 ) {
+        return kEta;
+       } 
+       else if ( TMath::Abs(maPdgcode) == 223 ) {
+        return kOmega;
+       } 
+       else if ( TMath::Abs(maPdgcode) == 333 ) {
+        return kPhi;
+       } 
+       else if ( TMath::Abs(maPdgcode) == 331 ) {
+        return kEtaPrime;
+       } 
+       else if ( TMath::Abs(maPdgcode) == 113 ) {
+        return kRho0;
+       } 
+       else if ( TMath::Abs(maPdgcode) == 321 || TMath::Abs(maPdgcode) == 130 ) {
+        return kKe3;
+       }
+       else{ 
+        origin = kElse;
+       }
+     }
+   }
+  return origin;
+}
 //__________________________________________
 Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel){
   //
@@ -1486,6 +2181,21 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
   else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
   else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
 
+  Double_t datamc[25]={-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999, -999, -999, -999, -999, -999, -999, -999, -999, -999};
+  Double_t xr[3]={-999,-999,-999};
+  datamc[0] = mesonID;
+  datamc[17] = mctrack->Pt(); //electron pt
+  datamc[18] = mctrack->Eta(); //electron eta
+
+  mctrack->XvYvZv(xr);
+  datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+  datamc[10] = xr[2];
+
+  TParticle *mcpart = mctrack->Particle();
+  if(mcpart){
+    datamc[14] = mcpart->GetUniqueID();
+  }
+
   if(!(mArr<0)){
      if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering 
         Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
@@ -1494,6 +2204,37 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
           if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
             mesonPt = mctrackmother->Pt(); //meson pt
             bgcategory = 1.;
+            datamc[1] = bgcategory;
+            datamc[2] = mesonPt;
+            mctrackmother->XvYvZv(xr);
+            datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+            datamc[12] = xr[2];
+
+            mcpart = mctrackmother->Particle();
+            if(mcpart){
+              datamc[15] = mcpart->GetUniqueID();
+            }
+            if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+              bgcategory = 2.;
+              datamc[1] = bgcategory;
+              //printf("I should be gamma meson = %d  mesonlabel= %d  NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries()); 
+              glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+              if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                datamc[3]=mctrackmother->PdgCode();
+                datamc[4]=mctrackmother->Pt();
+                if(TMath::Abs(mctrackmother->PdgCode())==310){
+                  bgcategory = 3.;
+                  glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+                  if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                    datamc[5]=mctrackmother->PdgCode();
+                    glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+                    if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                       datamc[6]=mctrackmother->PdgCode();
+                    }
+                  }
+                }
+              }
+            } 
           } 
         }
      }
@@ -1501,16 +2242,50 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
         Int_t glabel=TMath::Abs(mctrack->GetMother()); 
         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
           mesonPt=mctrackmother->Pt(); //meson pt
-          bgcategory = -1.;
+          if(mesonID==kEta) bgcategory = -1.41; // to consider new branching ratio for the eta Dalitz decay
+          else bgcategory = -1.;
+          datamc[1] = bgcategory;
+          datamc[2] = mesonPt;
+          mctrackmother->XvYvZv(xr);
+          datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+          datamc[12] = xr[2];
+
+          mcpart = mctrackmother->Particle();
+          if(mcpart){
+            datamc[15] = mcpart->GetUniqueID();
+          }
+          if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+            if(mesonID==kEta) bgcategory = -2.82; // to consider new branching ratio for the eta Dalitz decay (1.41)
+            else bgcategory = -2.;
+            datamc[1] = bgcategory;
+            glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+            if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+              datamc[3]=mctrackmother->PdgCode();
+              datamc[4]=mctrackmother->Pt();
+              if(TMath::Abs(mctrackmother->PdgCode())==310){
+               bgcategory = -3.;
+               glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+               if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                 datamc[5]=mctrackmother->PdgCode();
+                 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+                 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                   datamc[6]=mctrackmother->PdgCode();
+                 }
+               }
+              }  
+            }
+          }
         }
      }
 
+
      if(fIsPbPb){
-       if(fCentrality < 0)return 0.;
-       weightElecBg=fElecBackgroundFactor[iBgLevel][fCentrality][mArr][kBgPtBins-1];                        
+       Int_t centBin = GetWeightCentralityBin(fPerCentrality);
+       if(centBin < 0)return 0.;
+       weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];                        
        for(int ii=0; ii<kBgPtBins; ii++){              
         if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
-          weightElecBg = fElecBackgroundFactor[iBgLevel][fCentrality][mArr][ii];
+          weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
           break;
         }
        }
@@ -1525,9 +2300,165 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
        }
      }    
   }
-  return bgcategory*weightElecBg;
+
+  datamc[13] = weightElecBg;
+  datamc[16] = Double_t(fContainerStep);
+
+  datamc[7] = fHfeImpactR;
+  datamc[8] = fHfeImpactnsigmaR;
+
+  datamc[19] = fRecPt;
+  datamc[20] = fRecEta;
+  datamc[21] = fRecPhi;
+  datamc[22] = fLyrhit;
+  datamc[23] = fLyrstat;
+  datamc[24] = fCentrality;
+
+  if(fIsDebugStreamerON && fTreeStream){
+   if(!iBgLevel){
+    (*fTreeStream)<<"nonhfeQA"<<
+        "mesonID="<<datamc[0]<<
+        "bgcategory="<<datamc[1]<<
+        "mesonPt="<<datamc[2]<<
+        "mesonMomPdg="<<datamc[3]<<
+        "mesonMomPt="<<datamc[4]<<
+        "mesonGMomPdg="<<datamc[5]<<
+        "mesonGGMomPdg="<<datamc[6]<<
+        "eIPAbs="<<datamc[7]<<
+        "eIPSig="<<datamc[8]<<
+        "eR="<<datamc[9]<<
+        "eZ="<<datamc[10]<<
+        "mesonR="<<datamc[11]<<
+        "mesonZ="<<datamc[12]<<
+        "weightElecBg="<<datamc[13]<< 
+        "eUniqID="<<datamc[14]<<
+        "mesonUniqID="<<datamc[15]<<
+        "containerStep="<<datamc[16]<<
+        "emcpt="<<datamc[17]<<
+        "emceta="<<datamc[18]<<
+        "erecpt="<<datamc[19]<<
+        "ereceta="<<datamc[20]<<
+        "erecphi="<<datamc[21]<< 
+        "itshit="<<datamc[22]<<
+       "itsstat="<<datamc[23]<<
+        "centrality="<<datamc[24]
+        << "\n";
+   }
+  }
+
+  Double_t returnval = bgcategory*weightElecBg;
+  if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
+
+  return returnval;
 }
 
+//__________________________________________
+Double_t AliHFEmcQA::GetWeightFactor(const AliAODMCParticle * const mcpart, const Int_t iBgLevel){
+  //
+  // Get weighting factor for the realistic background estimation, for three possible background yield levels, indicated by the argument "iLevel": the best estimate (0), the lower uncertainty level (1), and the upper uncertainty level (2)
+  //
+
+  if (!mcpart) return 0;
+  if (!fMCArray) return 0;
+
+  Double_t weightElecBg = 0.;
+  Double_t mesonPt = 0.;
+  Double_t bgcategory = 0.;
+  Int_t mArr = -1;
+  Int_t mesonID = GetElecSource(mcpart);
+  if(mesonID==kGammaPi0 || mesonID==kPi0) mArr=0;                //pion
+  else if(mesonID==kGammaEta || mesonID==kEta) mArr=1;           //eta
+  else if(mesonID==kGammaOmega || mesonID==kOmega) mArr=2;       //omega
+  else if(mesonID==kGammaPhi || mesonID==kPhi) mArr=3;           //phi
+  else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
+  else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
+
+  if(!(mArr<0)){
+
+     AliAODMCParticle *mctrackmother = NULL; // will change all the time
+
+     if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering
+        Int_t iLabel = mcpart->GetMother(); //gamma label
+        if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
+        iLabel = mctrackmother->GetMother(); //gamma's mother's label
+        if(!(mctrackmother= dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
+        mesonPt = mctrackmother->Pt(); //meson pt
+        bgcategory = 1.;
+
+        //if(glabel>fMCEvent->GetNumberOfPrimaries()) bgcategory = 2.;
+     }
+     else{ // nonHFE except for the conversion electron
+        Int_t iLabel = mcpart->GetMother(); //meson label
+        if(!(mctrackmother = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(iLabel))))) return 0;
+        mesonPt = mctrackmother->Pt(); //meson pt
+        if(mesonID==kEta) bgcategory = -1.41; // to consider new branching ratio for the eta Dalitz decay
+        else bgcategory = -1.;
+
+        //if(glabel>fMCEvent->GetNumberOfPrimaries()) bgcategory = -2.;
+     }
+
+     if(fIsPbPb){
+       Int_t centBin = GetWeightCentralityBin(fPerCentrality);
+       if(centBin < 0)return 0.;
+       weightElecBg=fElecBackgroundFactor[iBgLevel][centBin][mArr][kBgPtBins-1];
+       for(int ii=0; ii<kBgPtBins; ii++){
+         if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
+           weightElecBg = fElecBackgroundFactor[iBgLevel][centBin][mArr][ii];
+           break;
+         }
+       }
+     }
+     else{
+       weightElecBg=fElecBackgroundFactor[iBgLevel][0][mArr][kBgPtBins-1];
+       for(int ii=0; ii<kBgPtBins; ii++){
+         if((mesonPt > fBinLimit[ii]) && (mesonPt < fBinLimit[ii+1])){
+           weightElecBg = fElecBackgroundFactor[iBgLevel][0][mArr][ii];
+           break;
+         }
+       }
+     }
+  }
+
+  Double_t returnval = bgcategory*weightElecBg;
+  if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
+
+  return returnval;
+}
+
+//__________________________________________
+Int_t AliHFEmcQA::GetMother(const AliVParticle * const mcpart) const {
+  //
+  // Wrapper to get the mother label
+  //
+  Int_t label = -1; 
+  TClass *type = mcpart->IsA();
+  if(type == AliMCParticle::Class()){
+    // ESD analysis
+    const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mcpart);
+    label = emcpart->GetMother();
+  } else if(type == AliAODMCParticle::Class()){
+    // AOD analysis
+    const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mcpart);
+    label = amcpart->GetMother();
+  }
+  return label;
+}
+//__________________________________________
+Int_t AliHFEmcQA::GetWeightCentralityBin(const Float_t percentile) const {
+  //
+  //translate the centrality percentile into the centrality bin of the reference weighting histograms for electron background
+  //
+
+  Float_t centralityLimits[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
+  Int_t bin = -1;
+  for(Int_t ibin = 0; ibin < 11; ibin++){
+    if(percentile >= centralityLimits[ibin] && percentile < centralityLimits[ibin+1]){
+      bin = ibin;
+      break;
+    }
+  } 
+  return bin;
+}
 //__________________________________________
 void AliHFEmcQA::AliHists::FillList(TList *l) const {
   //
@@ -1551,6 +2482,41 @@ void AliHFEmcQA::AliHistsComm::FillList(TList *l) const {
   if(fPtCorrDp) l->Add(fPtCorrDp);
   if(fPtCorrD0) l->Add(fPtCorrD0);
   if(fPtCorrDrest) l->Add(fPtCorrDrest);
+
+  if(fPtCorrDinein) l->Add(fPtCorrDinein);
+  if(fPtCorrDineout) l->Add(fPtCorrDineout);
+  if(fPtCorrDoutein) l->Add(fPtCorrDoutein);
+  if(fPtCorrDouteout) l->Add(fPtCorrDouteout);
+  if(fPtCorrDpDinein) l->Add(fPtCorrDpDinein);
+  if(fPtCorrDpDineout) l->Add(fPtCorrDpDineout);
+  if(fPtCorrDpDoutein) l->Add(fPtCorrDpDoutein);
+  if(fPtCorrDpDouteout) l->Add(fPtCorrDpDouteout);
+  if(fPtCorrD0Dinein) l->Add(fPtCorrD0Dinein);
+  if(fPtCorrD0Dineout) l->Add(fPtCorrD0Dineout);
+  if(fPtCorrD0Doutein) l->Add(fPtCorrD0Doutein);
+  if(fPtCorrD0Douteout) l->Add(fPtCorrD0Douteout);
+  if(fPtCorrDrestDinein) l->Add(fPtCorrDrestDinein);
+  if(fPtCorrDrestDineout) l->Add(fPtCorrDrestDineout);
+  if(fPtCorrDrestDoutein) l->Add(fPtCorrDrestDoutein);
+  if(fPtCorrDrestDouteout) l->Add(fPtCorrDrestDouteout);
+
+  if(fEtaCorrD) l->Add(fEtaCorrD);
+  if(fEtaCorrDp) l->Add(fEtaCorrDp);
+  if(fEtaCorrD0) l->Add(fEtaCorrD0);
+  if(fEtaCorrDrest) l->Add(fEtaCorrDrest);
+
+  if(fEtaCorrGD) l->Add(fEtaCorrGD);
+  if(fEtaCorrGDp) l->Add(fEtaCorrGDp);
+  if(fEtaCorrGD0) l->Add(fEtaCorrGD0);
+  if(fEtaCorrGDrest) l->Add(fEtaCorrGDrest);
+
+  if(fEtaCorrB) l->Add(fEtaCorrB);
+  if(fEtaCorrGB) l->Add(fEtaCorrGB);
+  if(fPtCorrBinein) l->Add(fPtCorrBinein);
+  if(fPtCorrBineout) l->Add(fPtCorrBineout);
+  if(fPtCorrBoutein) l->Add(fPtCorrBoutein);
+  if(fPtCorrBouteout) l->Add(fPtCorrBouteout);
+
   if(fDePtRatio) l->Add(fDePtRatio);
   if(feDistance) l->Add(feDistance);
   if(fDeDistance) l->Add(fDeDistance);