]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliAnalysisTaskHFE.cxx
Update of hfe code
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFE.cxx
index 03adfb44993304d267b153b18a94744b68ae087b..41a261d539aff7d82b83ba498bd4be5eba548fb6 100644 (file)
@@ -43,6 +43,7 @@
 #include <TF1.h>
 #include <TTree.h>
 
+#include "AliESDtrackCuts.h"
 #include "AliAnalysisManager.h"
 #include "AliAODInputHandler.h"
 #include "AliAODMCParticle.h"
@@ -70,6 +71,7 @@
 #include "AliHFEcontainer.h"
 #include "AliHFEcuts.h"
 #include "AliHFEelecbackground.h"
+#include "AliHFENonPhotonicElectron.h"
 #include "AliHFEmcQA.h"
 #include "AliHFEpairs.h"
 #include "AliHFEpid.h"
 #include "AliHFEtools.h"
 #include "AliHFEvarManager.h"
 #include "AliAnalysisTaskHFE.h"
+#include "AliAODMCHeader.h"
+#include "TClonesArray.h"
 
 ClassImp(AliAnalysisTaskHFE)
 
 //____________________________________________________________
 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
-  AliAnalysisTaskSE("PID efficiency Analysis")
+AliAnalysisTaskSE("PID efficiency Analysis")
+  , fAODMCHeader(NULL)
+  , fAODArrayMCInfo(NULL)
   , fQAlevel(0)
   , fPlugins(0)
   , fFillSignalOnly(kTRUE)
   , fFillNoCuts(kFALSE)
+  , fUseFilterAOD(kTRUE)
+  , fApplyCutAOD(kFALSE)
+  , fFilter(1<<4)
   , fBackGroundFactorApply(kFALSE)
   , fRemovePileUp(kFALSE)
   , fIdentifiedAsPileUp(kFALSE)
@@ -99,9 +108,11 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fPassTheEventCut(kFALSE)
   , fRejectKinkMother(kTRUE)
   , fisppMultiBin(kFALSE)
+  , fPbPbUserCentralityBinning(kFALSE)
   , fisNonHFEsystematics(kFALSE)
   , fSpecialTrigger(NULL)
   , fCentralityF(-1)
+  , fCentralityPercent(-1)
   , fContributors(0.5)
   , fWeightBackGround(0.)
   , fVz(0.0)
@@ -124,12 +135,14 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fMCQA(NULL)
   , fTaggedTrackAnalysis(NULL)
   , fExtraCuts(NULL)
+  , fBackgroundSubtraction(NULL)
   , fQA(NULL)
   , fOutput(NULL)
   , fHistMCQA(NULL)
   , fHistSECVTX(NULL)
   , fHistELECBACKGROUND(NULL)
   , fQACollection(NULL)
+  , fQAAODCollection(NULL)
 {
   //
   // Dummy constructor
@@ -138,15 +151,23 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
+  memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
+
+
 }
 
 //____________________________________________________________
 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   AliAnalysisTaskSE(name)
+  , fAODMCHeader(NULL)
+  , fAODArrayMCInfo(NULL)
   , fQAlevel(0)
   , fPlugins(0)
   , fFillSignalOnly(kTRUE)
   , fFillNoCuts(kFALSE)
+  , fUseFilterAOD(kTRUE)
+  , fApplyCutAOD(kFALSE)
+  , fFilter(1<<4)
   , fBackGroundFactorApply(kFALSE)
   , fRemovePileUp(kFALSE)
   , fIdentifiedAsPileUp(kFALSE)
@@ -154,9 +175,11 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fPassTheEventCut(kFALSE)  
   , fRejectKinkMother(kTRUE)
   , fisppMultiBin(kFALSE)
+  , fPbPbUserCentralityBinning(kFALSE)
   , fisNonHFEsystematics(kFALSE)
   , fSpecialTrigger(NULL)
   , fCentralityF(-1)
+  , fCentralityPercent(-1)
   , fContributors(0.5)
   , fWeightBackGround(0.)
   , fVz(0.0)
@@ -179,12 +202,14 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fMCQA(NULL)
   , fTaggedTrackAnalysis(NULL)
   , fExtraCuts(NULL)
+  , fBackgroundSubtraction(NULL)
   , fQA(NULL)
   , fOutput(NULL)
   , fHistMCQA(NULL)
   , fHistSECVTX(NULL)
   , fHistELECBACKGROUND(NULL)
   , fQACollection(0x0)
+  , fQAAODCollection(NULL)
 {
   //
   // Default constructor
@@ -200,15 +225,22 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
+  memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
+
 }
 
 //____________________________________________________________
 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   AliAnalysisTaskSE(ref)
+  , fAODMCHeader(NULL)
+  , fAODArrayMCInfo(NULL)
   , fQAlevel(0)
   , fPlugins(0)
   , fFillSignalOnly(ref.fFillSignalOnly)
   , fFillNoCuts(ref.fFillNoCuts)
+  , fUseFilterAOD(ref.fUseFilterAOD)
+  , fApplyCutAOD(ref.fApplyCutAOD)
+  , fFilter(ref.fFilter) 
   , fBackGroundFactorApply(ref.fBackGroundFactorApply)
   , fRemovePileUp(ref.fRemovePileUp)
   , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
@@ -216,9 +248,11 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fPassTheEventCut(ref.fPassTheEventCut)
   , fRejectKinkMother(ref.fRejectKinkMother)
   , fisppMultiBin(ref.fisppMultiBin)
+  , fPbPbUserCentralityBinning(ref.fPbPbUserCentralityBinning)
   , fisNonHFEsystematics(ref.fisNonHFEsystematics)
   , fSpecialTrigger(ref.fSpecialTrigger)
   , fCentralityF(ref.fCentralityF)
+  , fCentralityPercent(ref.fCentralityPercent)
   , fContributors(ref.fContributors)
   , fWeightBackGround(ref.fWeightBackGround)
   , fVz(ref.fVz)
@@ -241,12 +275,14 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fMCQA(NULL)
   , fTaggedTrackAnalysis(NULL)
   , fExtraCuts(NULL)
+  , fBackgroundSubtraction(NULL)
   , fQA(NULL)
   , fOutput(NULL)
   , fHistMCQA(NULL)
   , fHistSECVTX(NULL)
   , fHistELECBACKGROUND(NULL)
   , fQACollection(NULL)
+  , fQAAODCollection(NULL)
 {
   //
   // Copy Constructor
@@ -270,10 +306,15 @@ void AliAnalysisTaskHFE::Copy(TObject &o) const {
   // Copy into object o
   //
   AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
+  target.fAODMCHeader = fAODMCHeader;
+  target.fAODArrayMCInfo = fAODArrayMCInfo;
   target.fQAlevel = fQAlevel;
   target.fPlugins = fPlugins;
   target.fFillSignalOnly = fFillSignalOnly;
   target.fFillNoCuts = fFillNoCuts;
+  target.fUseFilterAOD = fUseFilterAOD;
+  target.fApplyCutAOD = fApplyCutAOD;
+  target.fFilter = fFilter;
   target.fBackGroundFactorApply = fBackGroundFactorApply;
   target.fRemovePileUp = fRemovePileUp;
   target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
@@ -281,9 +322,11 @@ void AliAnalysisTaskHFE::Copy(TObject &o) const {
   target.fPassTheEventCut = fPassTheEventCut;
   target.fRejectKinkMother = fRejectKinkMother;
   target.fisppMultiBin =   fisppMultiBin;
+  target.fPbPbUserCentralityBinning = fPbPbUserCentralityBinning;
   target.fisNonHFEsystematics = fisNonHFEsystematics;
   target.fSpecialTrigger = fSpecialTrigger;
   target.fCentralityF = fCentralityF;
+  target.fCentralityPercent = fCentralityPercent;
   target.fContributors = fContributors;
   target.fWeightBackGround = fWeightBackGround;
   target.fVz = fVz;
@@ -306,12 +349,14 @@ void AliAnalysisTaskHFE::Copy(TObject &o) const {
   target.fMCQA = fMCQA;
   target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
   target.fExtraCuts = fExtraCuts;
+  target.fBackgroundSubtraction = fBackgroundSubtraction;
   target.fQA = fQA;
   target.fOutput = fOutput;
   target.fHistMCQA = fHistMCQA;
   target.fHistSECVTX = fHistSECVTX;
   target.fHistELECBACKGROUND = fHistELECBACKGROUND;
   target.fQACollection = fQACollection;
+  target.fQAAODCollection = fQAAODCollection;
 }
 
 //____________________________________________________________
@@ -328,6 +373,7 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
   if(fSecVtx) delete fSecVtx;
   if(fMCQA) delete fMCQA;
   if(fElecBackGround) delete fElecBackGround;
+  if(fBackgroundSubtraction) delete fBackgroundSubtraction;
   if(fSpecialTrigger) delete fSpecialTrigger;
   // Delete output objects only if we are not running in PROOF mode because otherwise this produces a crash during merging
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
@@ -351,6 +397,13 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   // Called once per worker
   //
   AliDebug(3, "Creating Output Objects");
+  
+  // Make lists for Output
+  if(!fQA) fQA = new TList;
+  fQA->SetOwner();
+  if(!fOutput) fOutput = new TList;
+  fOutput->SetOwner();
+  
   // Automatic determination of the analysis mode
   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
@@ -361,6 +414,15 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
       SetHasMCData();
   }
   printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
+  if(IsAODanalysis()) {
+    printf("AOD filter: %s \n", fUseFilterAOD ? "Yes" : "No");
+    if(fUseFilterAOD) printf("AOD filter used: %d \n", fFilter);
+    // First Part: Make QA histograms
+    fQAAODCollection = new AliHFEcollection("TaskQAAOD", "QA histos from the AOD Electron Task");
+    fQAAODCollection->CreateTH1F("Filterorigin", "AOD filter of tracks at the origin", 21, -1, 20);
+    fQAAODCollection->CreateTH1F("Filterend", "AOD filter of tracks after all cuts", 21, -1, 20);
+    fQA->Add(fQAAODCollection);
+  }
   printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
 
   // Enable Trigger Analysis
@@ -368,13 +430,6 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   fTriggerAnalysis->EnableHistograms();
   fTriggerAnalysis->SetAnalyzeMC(HasMCData());
 
-
-  // Make lists for Output
-  if(!fQA) fQA = new TList;
-  fQA->SetOwner();
-  if(!fOutput) fOutput = new TList;
-  fOutput->SetOwner();
-
   // First Part: Make QA histograms
   fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
   fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
@@ -398,6 +453,15 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   }
   fPID->SortDetectors();
 
+  // Background subtraction-------------------------------------------------------------------
+  if (GetPlugin(kNonPhotonicElectron)) {
+    if(!fBackgroundSubtraction) fBackgroundSubtraction = new AliHFENonPhotonicElectron();
+    fBackgroundSubtraction->Init();
+    fOutput->Add(fBackgroundSubtraction->GetListOutput());
+  }
+  //------------------------------------------------------------------------------------------
+
+
   // Initialize correction Framework and Cuts
   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack + AliHFEcuts::kNcutStepsSecvtxTrack;
   fCFM = new AliCFManager;
@@ -431,6 +495,9 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
     fHistMCQA->SetOwner();
     if(IsPbPb()) fMCQA->SetPbPb();
     if(fisppMultiBin) fMCQA->SetPPMultiBin();
+    if(TestBit(kTreeStream)){
+      fMCQA->EnableDebugStreamer();
+    }
     fMCQA->CreatDefaultHistograms(fHistMCQA);
     fMCQA->SetBackgroundWeightFactor(fElecBackgroundFactor[0][0][0],fBinLimit);
     fQA->Add(fHistMCQA);
@@ -447,7 +514,7 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
     fSecVtx->CreateHistograms(fHistSECVTX);
     fOutput->Add(fHistSECVTX);
   }
-  
+
   // background----------------------------------
   if (GetPlugin(kIsElecBackGround)) {
     AliInfo("Electron BackGround Analysis on");
@@ -499,6 +566,9 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   }
 
   PrintStatus();
+  // Done!!!
+  PostData(1, fOutput);
+  PostData(2, fQA);
 }
 
 //____________________________________________________________
@@ -506,20 +576,26 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
   //
   // Run the analysis
   // 
+
+  //printf("test00\n");
+
   AliDebug(3, "Starting Single Event Analysis");
   if(!fInputEvent){
     AliError("Reconstructed Event not available");
+    //printf("Reconstructed Event not available");
     return;
   }
-  if(HasMCData()){
+  if(HasMCData() && IsESDanalysis()){
     AliDebug(4, Form("MC Event: %p", fMCEvent));
     if(!fMCEvent){
       AliError("No MC Event, but MC Data required");
+      //printf("No MC Event, but MC Data required");
       return;
     }
   }
   if(!fCuts){
     AliError("HFE cuts not available");
+    //printf("HFE cuts not available");
     return;
   }
   if(!fPID->IsInitialized()){
@@ -527,6 +603,8 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
     fPID->InitializePID(fInputEvent->GetRunNumber());
   }
 
+  //printf("test0\n");
+
   // Initialize hadronic background from OADB Container
   AliDebug(2, Form("Apply background factors: %s, OADB Container %p", fBackGroundFactorApply ? "Yes" : "No", fHadronBackgroundOADB));
   if(fBackGroundFactorApply && !TestBit(kBackgroundInitialized)){ 
@@ -536,6 +614,8 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
     SetBit(kBackgroundInitialized); 
   }
 
+  //printf("test1\n");
+
   if(IsESDanalysis() && HasMCData()){
     // Protect against missing MC trees
     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
@@ -546,34 +626,75 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
     if(!mcH->InitOk()) return;
     if(!mcH->TreeK()) return;
     if(!mcH->TreeTR()) return;
+    // Background subtraction-------------------------------------------------------------------
+    if(GetPlugin(kNonPhotonicElectron) && fBackgroundSubtraction){
+      fBackgroundSubtraction->SetMCEvent(fMCEvent);
+    }
+    //------------------------------------------------------------------------------------------
+  }
+
+  if(IsAODanalysis() && HasMCData()){
+    // take MC info
+    AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
+    if(!aodE){ 
+      AliError("No AOD Event");
+      return;
+    }
+    fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
+    if(!fAODMCHeader){ 
+      AliError("No AliAODMCHeader");
+      //printf("No AliAODMCHeader");
+      return;
+    }
+    fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+    if(!fAODArrayMCInfo){ 
+      AliError("No AOD MC particles");
+      //printf("No AOD MC particles");
+      return;
+    }
+    fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
+    // Background subtraction-------------------------------------------------------------------
+    if (GetPlugin(kNonPhotonicElectron)) {
+      fBackgroundSubtraction->SetAODArrayMCInfo(fAODArrayMCInfo);
+    }
+    //------------------------------------------------------------------------------------------
   }
 
+  //printf("test2\n");
+
   // need the centrality for everything (MC also)
   fCentralityF = -1;
   if(!ReadCentrality()) fCentralityF = -1;
   //printf("pass centrality\n");
-  //printf("Reading fCentralityF %f\n",fCentralityF);
+  //printf("Reading fCentralityF %d\n",fCentralityF);
   
   // See if pile up and z in the range
   RejectionPileUpVertexRangeEventCut();
 
+  //printf("test3\n");
+
   // Protect agains missing 
   if(HasMCData()){
     //printf("Has MC data\n");
     fSignalCuts->SetMCEvent(fMCEvent);
     ProcessMC();  // Run the MC loop + MC QA in case MC Data are available
   }
-
+  
   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
   if(!pidResponse){
     AliDebug(1, "Using default PID Response");
-    pidResponse = AliHFEtools::GetDefaultPID(HasMCData(), fInputEvent->IsA() == AliAODEvent::Class());
+    pidResponse = AliHFEtools::GetDefaultPID(HasMCData(), fInputEvent->IsA() == AliESDEvent::Class());
   }
   fPID->SetPIDResponse(pidResponse);
   if(fPIDpreselect) fPIDpreselect->SetPIDResponse(pidResponse);
 
+  // Background subtraction-------------------------------------------------------------------
+  if(GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->InitRun(fInputEvent,pidResponse);
+  //------------------------------------------------------------------------------------------
+
   // Event loop
   if(IsAODanalysis()){
+    //printf("test4\n");
     ProcessAOD();
   } else {
     const char *specialTrigger = GetSpecialTrigger(fInputEvent->GetRunNumber());
@@ -586,7 +707,7 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
         return;
       } else AliDebug(2, "Event Selected");
     } else AliDebug(2, "No Special Trigger requested");
-
+    
     ProcessESD();
   }
   // Done!!!
@@ -623,7 +744,7 @@ void AliAnalysisTaskHFE::Terminate(Option_t *){
       return;
     }
     postanalysis.SetTaskQA(qalist);
-    printf("Running post analysis\n");
+    //printf("Running post analysis\n");
     //if(HasMCData())
     postanalysis.DrawMCSignal2Background();
     postanalysis.DrawEfficiency();
@@ -675,7 +796,8 @@ void AliAnalysisTaskHFE::ProcessMC(){
   //
   AliDebug(3, "Processing MC Information");
   Double_t eventContainer [4];
-  eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
+  if(IsESDanalysis()) eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
+  else eventContainer[0] = fAODMCHeader->GetVtxZ();
   eventContainer[2] = fCentralityF;
   eventContainer[3] = fContributors;
   fVz = eventContainer[0];
@@ -692,6 +814,7 @@ void AliAnalysisTaskHFE::ProcessMC(){
         fMCQA->SetMCEvent(fMCEvent);
         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
         fMCQA->SetCentrality(fCentralityF);
+        fMCQA->SetPercentrality(fCentralityPercent);
         if(IsPbPb()) { fMCQA->SetPbPb();}
         else
         {
@@ -710,19 +833,9 @@ void AliAnalysisTaskHFE::ProcessMC(){
           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
-          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0); // no accept cut
-          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
-          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 0); // no accept cut
-          if (TMath::Abs(mcpart->Eta()) < 0.9) {
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-          }
-          if (TMath::Abs(AliHFEtools::GetRapidity(mcpart)) < 0.5) {
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-          }
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG); // no accept cut
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG); // no accept cut
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG); // no accept cut
         }
         //fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
         //fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
@@ -738,9 +851,23 @@ void AliAnalysisTaskHFE::ProcessMC(){
   }
   // Run MC loop
   AliVParticle *mctrack = NULL;
-  AliDebug(3, Form("Number of Tracks: %d", fMCEvent->GetNumberOfTracks()));
-  for(Int_t imc = 0; imc <fMCEvent->GetNumberOfTracks(); imc++){
-    if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
+  Int_t numberofmctracks = 0;
+  if(IsESDanalysis()){
+    numberofmctracks = fMCEvent->GetNumberOfTracks();
+  }
+  else {
+    numberofmctracks = fAODArrayMCInfo->GetEntriesFast();
+  }
+  AliDebug(3, Form("Number of Tracks: %d",numberofmctracks));
+  //printf("Number of MC track %d\n",numberofmctracks);
+  for(Int_t imc = 0; imc <numberofmctracks; imc++){
+    if(IsESDanalysis()) {
+      if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
+    }
+    else {
+      if(!(mctrack = (AliVParticle *) fAODArrayMCInfo->At(imc))) continue;
+    }
+    //printf("Test in ProcessMC\n");
     AliDebug(4, "Next MC Track");
     if(ProcessMCtrack(mctrack)) nElectrons++;
   }
@@ -755,6 +882,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
   // Run Analysis of reconstructed event in ESD Mode
   // Loop over Tracks, filter according cut steps defined in AliHFEcuts
   //
+  AliDebug(1, Form("Task %s", GetName()));
   AliDebug(3, "Processing ESD Event");
   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
   if(!fESD){
@@ -803,7 +931,6 @@ void AliAnalysisTaskHFE::ProcessESD(){
   if(!fPassTheEventCut) return;
   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
 
 
   fContainer->NewEvent();
 
@@ -829,6 +956,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
   memset(container, 0, sizeof(Double_t) * 10);
   // container for the output THnSparse
   Double_t dataE[6]; // [pT, eta, Phi, type, 'C' or 'B']
+  Double_t dataDca[6]; // [source, pT, dca, centrality]
   Int_t nElectronCandidates = 0;
   AliESDtrack *track = NULL, *htrack = NULL;
   AliMCParticle *mctrack = NULL;
@@ -839,6 +967,13 @@ void AliAnalysisTaskHFE::ProcessESD(){
 
   fCFM->SetRecEventInfo(fESD);
 
+  // Get Number of contributors to the primary vertex for multiplicity-dependent correction
+  Int_t ncontribVtx = 0;
+  const AliESDVertex *priVtx = fESD->GetPrimaryVertexTracks();
+  if(priVtx){
+    ncontribVtx = priVtx->GetNContributors();
+  }
+
   // minjung for IP QA(temporary ~ 2weeks)
   if(!fExtraCuts){
     fExtraCuts = new AliHFEextraCuts("hfetmpCuts","HFE tmp Cuts");
@@ -853,6 +988,14 @@ void AliAnalysisTaskHFE::ProcessESD(){
     fElecBackGround->Reset();
     
   } // end of electron background analysis
+  
+
+  // Background subtraction-------------------------------------------------------------------
+  if (GetPlugin(kNonPhotonicElectron)) {
+    fBackgroundSubtraction->FillPoolAssociatedTracks(fInputEvent,fCentralityF);
+  }
+  //------------------------------------------------------------------------------------------
+
   //
   // Loop ESD
   //
@@ -871,14 +1014,18 @@ void AliAnalysisTaskHFE::ProcessESD(){
     if(fTaggedTrackAnalysis && v0pid > -1){ 
       AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
       fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
+      AliDebug(1, "V0 PID done");
     }
  
+
+    //Fill non-HFE source containers at reconstructed events cut step
     AliDebug(3, Form("Doing track %d, %p", itrack, track));
-     
+
+
     //////////////////////////////////////
     // preselect
-    /////////////////////////////////////
-    if(fPIDpreselect && fCutspreselect) {
+    //////////////////////////////////////
+    if(fPIDpreselect || fCutspreselect) {
       if(!PreSelectTrack(track)) continue;
     }
 
@@ -906,7 +1053,44 @@ void AliAnalysisTaskHFE::ProcessESD(){
         fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
       }
     }
-
+        
+    if(fisNonHFEsystematics && IsPbPb())  {
+      //FillProductionVertex(track);     
+      if(fMCQA && signal){ 
+       fMCQA->SetCentrality(fCentralityF);
+       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){
+         Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.};
+          
+         fMCQA->SetTrkKine(track->Pt(),track->Eta(), track->Phi());
+         fMCQA->SetContainerStep(4);
+          
+         weightElecBgV0[0] = fMCQA->GetWeightFactor(mctrack, 0); // positive:conversion e, negative: nonHFE 
+          
+         //Fill additional containers for electron source distinction
+         Int_t elecSource = 0;
+         elecSource = fMCQA->GetElecSource(mctrack->Particle());
+         const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
+         const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
+         Int_t iName = 0;
+         for(Int_t iSource = AliHFEmcQA::kPi0; iSource <=  AliHFEmcQA::kGammaRho0; iSource++){
+           if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
+           if(elecSource == iSource){
+             
+             if(weightElecBgV0[0]>0){ 
+               fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[0]), 4, kTRUE, weightElecBgV0[0]);
+             } 
+             else if(weightElecBgV0[0]<0){ 
+               fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[0]), 4, kTRUE, -1*weightElecBgV0[0]);
+             }           
+             break;
+           }
+           iName++;
+           if(iName == kElecBgSpecies)iName = 0;
+         }
+       }
+      }
+    }
+  
     // RecKine: ITSTPC cuts  
     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
     
@@ -937,13 +1121,25 @@ void AliAnalysisTaskHFE::ProcessESD(){
     if(HasMCData()){
       FillProductionVertex(track);
 
-      if(fMCQA){
+      if(fMCQA && signal){
         fMCQA->SetCentrality(fCentralityF);
         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){
          Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.};
+         Double_t hfeimpactRtmp=0., hfeimpactnsigmaRtmp=0.;
+         fExtraCuts->GetHFEImpactParameters(track, hfeimpactRtmp, hfeimpactnsigmaRtmp);
+         UChar_t itsPixel = track->GetITSClusterMap();
+         Double_t ilyrhit=0, ilyrstat=0;
+         for(Int_t ilyr=0; ilyr<6; ilyr++){
+           if(TESTBIT(itsPixel, ilyr)) ilyrhit += TMath::Power(2,ilyr);
+           if(fExtraCuts->CheckITSstatus(fExtraCuts->GetITSstatus(track,ilyr))) ilyrstat += TMath::Power(2,ilyr);
+         }
+         fMCQA->SetITSInfo(ilyrhit,ilyrstat);
+         fMCQA->SetHFEImpactParameters(hfeimpactRtmp, hfeimpactnsigmaRtmp);
+         fMCQA->SetTrkKine(track->Pt(),track->Eta(), track->Phi());
+         fMCQA->SetContainerStep(3);
          for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
            weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE 
-           if(!fisNonHFEsystematics)break;   
+           if(!fisNonHFEsystematics || IsPbPb())break;   
          }
          
          if(fisNonHFEsystematics){
@@ -957,8 +1153,13 @@ void AliAnalysisTaskHFE::ProcessESD(){
              if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
              if(elecSource == iSource){
                for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
-                 if(weightElecBgV0[iLevel]>0){ fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);}
-                 else if(weightElecBgV0[iLevel]<0){ fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);}
+                 if(weightElecBgV0[iLevel]>0){ 
+                   fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);
+                 } 
+                 else if(weightElecBgV0[iLevel]<0){ 
+                   fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);
+                 }
+                 if(IsPbPb())break;
                }
                break;
              }
@@ -967,15 +1168,67 @@ void AliAnalysisTaskHFE::ProcessESD(){
            }
          }
          //else{
-           if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]);
-           else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]);
+           if(weightElecBgV0[0]>0) {
+            fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]);
+            fVarManager->FillContainer(fContainer, "conversionElecs", 4, kTRUE, weightElecBgV0[0]);
+          }
+           else if(weightElecBgV0[0]<0) {
+            fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]);
+            fVarManager->FillContainer(fContainer, "mesonElecs", 4, kTRUE, -1*weightElecBgV0[0]);
+          }
            //}
         }
       }
+
+      Double_t hfeimpactR4all=0., hfeimpactnsigmaR4all=0.;
+      Int_t sourceDca =-1;
+      if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 211)){
+        if(track->Pt()>4.){
+          fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
+          dataDca[0]=0; //pion
+          dataDca[1]=track->Pt();
+          dataDca[2]=hfeimpactR4all;
+          dataDca[3]=fCentralityF;
+          dataDca[4] = v0pid;
+          dataDca[5] = double(track->Charge());
+          fQACollection->Fill("Dca", dataDca);
+        }
+      }
+      else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){ // to increas statistics for Martin
+        if(signal){
+          fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
+          if(fSignalCuts->IsCharmElectron(track)){
+            sourceDca=1;
+          }
+          else if(fSignalCuts->IsBeautyElectron(track)){
+            sourceDca=2;
+          }
+          else if(fSignalCuts->IsGammaElectron(track)){
+            sourceDca=3;
+          }
+          else if(fSignalCuts->IsNonHFElectron(track)){
+            sourceDca=4;
+          }
+          else if(fSignalCuts->IsJpsiElectron(track)){
+            sourceDca=5;
+          }
+          else {
+            sourceDca=6;
+          }
+          dataDca[0]=sourceDca;
+          dataDca[1]=track->Pt();
+          dataDca[2]=hfeimpactR4all;
+          dataDca[3]=fCentralityF;
+          dataDca[4] = v0pid;
+          dataDca[5] = double(track->Charge());
+          if(signal) fQACollection->Fill("Dca", dataDca);
+        }
+      }
     }
 
     if(TMath::Abs(track->Eta()) < 0.5){
-      fQACollection->Fill("TPCdEdxBeforePID", track->P(), track->GetTPCsignal());
+      if(track->GetInnerParam())
+        fQACollection->Fill("TPCdEdxBeforePID", track->GetInnerParam()->P(), track->GetTPCsignal());
       fQACollection->Fill("TPCnSigmaBeforePID", track->P(), fInputHandler->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron));
     }
 
@@ -984,11 +1237,23 @@ void AliAnalysisTaskHFE::ProcessESD(){
     hfetrack.SetRecTrack(track);
     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
     hfetrack.SetCentrality(fCentralityF);
+    hfetrack.SetMulitplicity(ncontribVtx);
     if(IsPbPb()) hfetrack.SetPbPb();
     else hfetrack.SetPP();
     fPID->SetVarManager(fVarManager);
     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
     nElectronCandidates++;
+    
+    // Background subtraction------------------------------------------------------------------------------------------
+    if (GetPlugin(kNonPhotonicElectron)) {
+      Int_t indexmother = -1;
+      Int_t mcsource = -1;
+      if(HasMCData()){
+       mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
+      }
+      fBackgroundSubtraction->LookAtNonHFE(itrack, track, fInputEvent, 1, fCentralityF, -1, mcsource, indexmother);
+    }
+    //-----------------------------------------------------------------------------------------------------------------
 
     // Temporary histogram for chi2/ITS cluster
     if(IsPbPb()) {
@@ -997,16 +1262,18 @@ void AliAnalysisTaskHFE::ProcessESD(){
             if(shared.CountBits() >= 2) sharebit=1;
 
            Double_t itschi2percluster = 0.0;
-           if(track->GetNcls(0) > 0) itschi2percluster = track->GetITSchi2()/static_cast<Double_t>(track->GetNcls(0));
+           Double_t itsnbcls = static_cast<Double_t>(track->GetNcls(0));
+           if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
 
             Double_t itsChi2[7] = {track->Pt(),track->Eta(), track->Phi(),
-                                  fCentralityF,track->GetTPCsignalN(), sharebit, itschi2percluster};
+                                  fCentralityF,track->GetTPCsignalN(), sharebit,itschi2percluster};
             fQACollection->Fill("fChi2perITScluster", itsChi2);
     }
     else{
-      
+
       Double_t itschi2percluster = 0.0;
-      if(track->GetNcls(0) > 0) itschi2percluster = track->GetITSchi2()/static_cast<Double_t>(track->GetNcls(0));
+      Double_t itsnbcls = static_cast<Double_t>(track->GetNcls(0));
+      if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
 
       Double_t itsChi2[3] = {track->Pt(), fCentralityF, itschi2percluster};
       fQACollection->Fill("fChi2perITScluster", itsChi2);
@@ -1021,9 +1288,9 @@ void AliAnalysisTaskHFE::ProcessESD(){
         Int_t glabel=TMath::Abs(mctrack->GetMother());
         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
           if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==321)
-            fQACollection->Fill("Ke3Kecorr",mctrackmother->Pt(),mctrack->Pt());
+            fQACollection->Fill("Ke3Kecorr",mctrack->Pt(),mctrackmother->Pt());
           else if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==130)
-            fQACollection->Fill("Ke3K0Lecorr",mctrackmother->Pt(),mctrack->Pt());
+            fQACollection->Fill("Ke3K0Lecorr",mctrack->Pt(),mctrackmother->Pt());
         }
       }
     }
@@ -1103,46 +1370,22 @@ void AliAnalysisTaskHFE::ProcessESD(){
       }
     } // end of electron background analysis
 
-
-
     if (GetPlugin(kDEstep)) { 
       Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.,};
       Int_t elecSource = 0;
-      // minjung for IP QA(temporary ~ 2weeks)
       Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
       fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
-      fQACollection->Fill("dataDca",track->Pt(),hfeimpactR);
-      fQACollection->Fill("dataDcaSig",track->Pt(),hfeimpactnsigmaR);
-      fQACollection->Fill("dataDcaSigDca",hfeimpactR,hfeimpactnsigmaR);
-      if(HasMCData()){
-        // minjung for IP QA(temporary ~ 2weeks)
-        if(fSignalCuts->IsCharmElectron(track)){
-          fQACollection->Fill("charmDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("charmDcaSig",track->Pt(),hfeimpactnsigmaR);
-        }
-        else if(fSignalCuts->IsBeautyElectron(track)){
-          fQACollection->Fill("beautyDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("beautyDcaSig",track->Pt(),hfeimpactnsigmaR);
-        }
-        else if(fSignalCuts->IsGammaElectron(track)){
-          fQACollection->Fill("conversionDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("conversionDcaSig",track->Pt(),hfeimpactnsigmaR);
-          fQACollection->Fill("conversionDcaSigDca",hfeimpactR,hfeimpactnsigmaR);
-        }
-        else if(fSignalCuts->IsNonHFElectron(track)){
-          fQACollection->Fill("nonhfeDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("nonhfeDcaSig",track->Pt(),hfeimpactnsigmaR);
-        }
-
+      if(HasMCData())
+      {
         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
-          fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
-          fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
-        }
-        if(fMCQA) {
+            fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
+        } 
+        if(fMCQA && signal) {
           
+          fMCQA->SetContainerStep(0);
           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE 
-            if(!fisNonHFEsystematics)break;        
+            if(!fisNonHFEsystematics || IsPbPb())break;        
           }
           
           if(fisNonHFEsystematics){
@@ -1157,6 +1400,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
                 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
                   if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, weightElecBgV0[iLevel]);
                   else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, -1*weightElecBgV0[iLevel]);
+                  if(IsPbPb())break;
                 }
                 break;
               }
@@ -1165,22 +1409,49 @@ void AliAnalysisTaskHFE::ProcessESD(){
             }
           }
           //else{
-          if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]);
-          else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]);
+          if(weightElecBgV0[0]>0) {
+           fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]);
+           fVarManager->FillContainer(fContainer, "conversionElecs", 5, kTRUE, weightElecBgV0[0]);
+         }
+          else if(weightElecBgV0[0]<0) {
+           fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]);
+           fVarManager->FillContainer(fContainer, "mesonElecs", 5, kTRUE, -1*weightElecBgV0[0]);
+         }  
           //}
           if(bTagged){ // bg estimation for the secondary vertex tagged signals
             if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 2, kFALSE, weightElecBgV0[0]);
             else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 2, kFALSE, -1*weightElecBgV0[0]);
           }
         }
-      }
+      } // end of MC
+
+      dataDca[0]=-1; //for data, don't know the origin
+      dataDca[1]=track->Pt();
+      dataDca[2]=hfeimpactR;
+      dataDca[3]=fCentralityF;
+      dataDca[4] = v0pid;
+      dataDca[5] = double(track->Charge());
+      if (!HasMCData()) fQACollection->Fill("Dca", dataDca);
+
       // Fill Containers for impact parameter analysis
       if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
+      if(signal) {
+        // Apply weight for background contamination after ip cut
+        if(fBackGroundFactorApply) {
+              fWeightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
+              if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
+              else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
+              // weightBackGround as special weight
+              fVarManager->FillContainer(fContainer, "hadronicBackground", 2, kFALSE, fWeightBackGround);
+        }
+      }
+
       if(HasMCData()){
-        if(fMCQA) {
+        if(fMCQA && signal) {
+          fMCQA->SetContainerStep(1);
           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE 
-            if(!fisNonHFEsystematics)break;        
+            if(!fisNonHFEsystematics || IsPbPb())break;        
           }       
           if(fisNonHFEsystematics){
             //Fill additional containers for electron source distinction             
@@ -1194,6 +1465,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
                 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
                   if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, weightElecBgV0[iLevel]);
                   else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, -1*weightElecBgV0[iLevel]);
+                  if(IsPbPb())break;
                 }
                 break;
               }
@@ -1202,8 +1474,14 @@ void AliAnalysisTaskHFE::ProcessESD(){
             }
           }
           // else{
-            if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]);
-            else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]);
+            if(weightElecBgV0[0]>0) {
+             fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]);
+             fVarManager->FillContainer(fContainer, "conversionElecs", 6, kTRUE, weightElecBgV0[0]);
+           }
+            else if(weightElecBgV0[0]<0) {
+             fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]);
+             fVarManager->FillContainer(fContainer, "mesonElecs", 6, kTRUE, -1*weightElecBgV0[0]);
+            }
             //}
         }
       }
@@ -1215,7 +1493,6 @@ void AliAnalysisTaskHFE::ProcessESD(){
       if(HasMCData()){
         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
           fQACollection->Fill("hadronsAfterIPcut",track->Pt());
-          fQACollection->Fill("hadronsAfterIPcutMC",mctrack->Pt());
         }
       }
     }
@@ -1230,25 +1507,30 @@ void AliAnalysisTaskHFE::ProcessAOD(){
   // Run Analysis in AOD Mode
   // Function is still in development
   //
+  //printf("Process AOD\n");
   AliDebug(3, "Processing AOD Event");
   Double_t eventContainer[4];
+  eventContainer[0] = 0.0;
   if(HasMCData()) eventContainer[0] = fVz;
   else {
-    eventContainer[0] = fInputEvent->GetPrimaryVertex()->GetZ();
+    if(fInputEvent->GetPrimaryVertex()) eventContainer[0] = fInputEvent->GetPrimaryVertex()->GetZ();
   }
   eventContainer[1] = 1.; // No Information available in AOD analysis, assume all events have V0AND
   eventContainer[2] = fCentralityF; 
   eventContainer[3] = fContributors; 
   
+  //printf("value event container %f, %f, %f, %f\n",eventContainer[0],eventContainer[1],eventContainer[2],eventContainer[3]);
+
   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
   if(!fAOD){
     AliError("AOD Event required for AOD Analysis");
       return;
   }
   
+  //printf("Will fill\n");
   //
   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
-
+  //printf("Fill\n");
   //
   if(fIdentifiedAsPileUp) return; 
   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
@@ -1260,50 +1542,231 @@ void AliAnalysisTaskHFE::ProcessAOD(){
   //
   if(!fPassTheEventCut) return;
   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
+  //printf("pass\n");
 
   fContainer->NewEvent();
+
+  fCFM->SetRecEventInfo(fAOD);
+
+  if(!fExtraCuts){
+    fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
+  }
+  fExtraCuts->SetRecEventInfo(fAOD);
+
+  // Get Number of contributors to the primary vertex for multiplicity-dependent correction
+  Int_t ncontribVtx = 0;
+  AliAODVertex *priVtx = fAOD->GetPrimaryVertex();
+  if(priVtx){
+    ncontribVtx = priVtx->GetNContributors();
+  }
+
+  // Look for kink mother
+  Int_t numberofvertices = fAOD->GetNumberOfVertices();
+  Double_t listofmotherkink[numberofvertices];
+  Int_t numberofmotherkink = 0;
+  for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
+    AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
+    if(!aodvertex) continue;
+    if(aodvertex->GetType()==AliAODVertex::kKink) {
+      AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
+      if(!mother) continue;
+      Int_t idmother = mother->GetID();
+      listofmotherkink[numberofmotherkink] = idmother;
+      //printf("ID %d\n",idmother);
+      numberofmotherkink++;
+    }
+  }
+  //printf("Number of kink mother in the events %d\n",numberofmotherkink);
+
+  // Background subtraction-------------------------------------------------------------------
+  if (GetPlugin(kNonPhotonicElectron)) {
+    fBackgroundSubtraction->FillPoolAssociatedTracks(fInputEvent,fCentralityF);
+  }
+  //------------------------------------------------------------------------------------------
+
+  // Loop over tracks
   AliAODTrack *track = NULL;
   AliAODMCParticle *mctrack = NULL;
   Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
+  Double_t dataDca[6]; // [source, pT, dca, centrality]
   Int_t nElectronCandidates = 0;
   Int_t pid;
   Bool_t signal;
+  //printf("Number of track %d\n",(Int_t) fAOD->GetNumberOfTracks());
   for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
     track = fAOD->GetTrack(itrack); mctrack = NULL;
     if(!track) continue;
-    if(track->GetFlags() != 1<<4) continue;  // Only process AOD tracks where the HFE is set
+    // Begining
+    Bool_t passone = kFALSE;  
+    fQAAODCollection->Fill("Filterorigin", -1);
+    for(Int_t k=0; k<20; k++) {
+      Int_t u = 1<<k;     
+      if((track->TestFilterBit(u))) {
+       fQAAODCollection->Fill("Filterorigin", k);
+       passone = kTRUE;
+      }
+    }
+    //if(!passone) printf("what is the filter %d\n",track->GetFilterMap());
 
+    if(fUseFilterAOD){
+      //printf("Filter of the track  %d\n",track->GetFilterMap());
+      if(!(track->TestFilterBit(fFilter))) continue;  // Only process AOD tracks where the HFE is set
+    }
+    //printf("Pass the flag\n");
+    
     signal = kTRUE;
     if(HasMCData()){
 
       Int_t label = TMath::Abs(track->GetLabel());
-      if(label)
-        mctrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(label));
+      if(label && label < fAODArrayMCInfo->GetEntriesFast())
+        mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(label));
         if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
     }
-    fVarManager->NewTrack(track, mctrack, fCentralityF, -1, kTRUE);
+    
+    fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
+    
+    if(fFillNoCuts) {
+      if(signal || !fFillSignalOnly){
+        fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
+        fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
+      }
+    }
+
+    if(fApplyCutAOD) {
+      //printf("Apply cuts\n");
+      // RecKine: ITSTPC cuts  
+      if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
+
+      // Reject kink mother
+      if(fRejectKinkMother) {
+       Bool_t kinkmotherpass = kTRUE;
+       for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
+         if(track->GetID() == listofmotherkink[kinkmother]) {
+           kinkmotherpass = kFALSE;
+           continue;
+         }
+       }
+       if(!kinkmotherpass) continue;
+      }
+      
+      // RecPrim
+      if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
+      
+      // HFEcuts: ITS layers cuts
+      if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
+      
+      // HFE cuts: TOF PID and mismatch flag
+      if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) continue;
+      
+      // HFE cuts: TPC PID cleanup
+      if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
+      
+      // HFEcuts: Nb of tracklets TRD0
+      if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
+    }
+
+    // Fill correlation maps before PID
+    if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
+      //printf("Fill correlation maps before PID\n");
+      fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
+    }
+
+    if(HasMCData()){
+      Double_t hfeimpactR4all=0., hfeimpactnsigmaR4all=0.;
+      Int_t sourceDca =-1;
+      if(mctrack && (TMath::Abs(mctrack->GetPdgCode()) == 211)){
+        if(track->Pt()>4.){
+          fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
+          dataDca[0]=0; //pion
+          dataDca[1]=track->Pt();
+          dataDca[2]=hfeimpactR4all;
+          dataDca[3]=fCentralityF;
+          dataDca[4] = -1; // not store V0 for the moment
+          dataDca[5] = double(track->Charge());
+          fQACollection->Fill("Dca", dataDca);
+        }
+      }
+      else if(mctrack && (TMath::Abs(mctrack->GetPdgCode()) == 11)){ // to increas statistics for Martin
+        if(signal){
+          fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
+          if(fSignalCuts->IsCharmElectron(track)){
+            sourceDca=1;
+          }
+          else if(fSignalCuts->IsBeautyElectron(track)){
+            sourceDca=2;
+          }
+          else if(fSignalCuts->IsGammaElectron(track)){
+            sourceDca=3;
+          }
+          else if(fSignalCuts->IsNonHFElectron(track)){
+            sourceDca=4;
+          }
+          else if(fSignalCuts->IsJpsiElectron(track)){
+            sourceDca=5;
+          }
+          else {
+            sourceDca=6;
+          }
+          dataDca[0]=sourceDca;
+          dataDca[1]=track->Pt();
+          dataDca[2]=hfeimpactR4all;
+          dataDca[3]=fCentralityF;
+          dataDca[4] = -1; // not store V0 for the moment
+          dataDca[5] = double(track->Charge());
+          if(signal) fQACollection->Fill("Dca", dataDca);
+        }
+      }
+    }
+
+    //printf("Will process to PID\n");
+
     // track accepted, do PID
     AliHFEpidObject hfetrack;
     hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
     hfetrack.SetRecTrack(track);
     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
     hfetrack.SetCentrality(fCentralityF);
+    hfetrack.SetMulitplicity(ncontribVtx); // for correction
+    if(IsPbPb()) hfetrack.SetPbPb();
+    else hfetrack.SetPP();
     fPID->SetVarManager(fVarManager);
-    if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;    // we will do PID here as soon as possible
+    if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;   
+    // we will do PID here as soon as possible
+
+    // Background subtraction----------------------------------------------------------------------------------------------
+    if (GetPlugin(kNonPhotonicElectron)) {
+      Int_t indexmother = -1;
+      Int_t mcsource = -1;
+      if(HasMCData())  mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
+      fBackgroundSubtraction->LookAtNonHFE(itrack, track, fInputEvent, 1, fCentralityF, -1,mcsource, indexmother);
+    }
+    //---------------------------------------------------------------------------------------------------------------------
+
+    // end AOD QA
+    fQAAODCollection->Fill("Filterend", -1);  
+    for(Int_t k=0; k<20; k++) {
+      Int_t u = 1<<k;
+      if((track->TestFilterBit(u))) {
+       fQAAODCollection->Fill("Filterend", k);
+      }
+    }
+       
     // Apply weight for background contamination
-    Double_t weightBackGround = 1.0;
-
-    // not correct treatment for pp
-    if(fBackGroundFactorApply==kTRUE) {
-            if(IsPbPb()) weightBackGround =  fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
-      else weightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P()));
-
-      if(weightBackGround < 0.0) weightBackGround = 0.0;
-            else if(weightBackGround > 1.0) weightBackGround = 1.0;
-            fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, weightBackGround);
+    //Double_t weightBackGround = 1.0;
+    if(signal) {
+      // Apply weight for background contamination
+      if(fBackGroundFactorApply) {
+       if(IsPbPb() && fCentralityF >= 0) fWeightBackGround =  fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
+       else    fWeightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
+       
+       if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
+       else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
+        // weightBackGround as special weight
+        fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
+      }
+      fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
     }
-
+    
     nElectronCandidates++;    
     if(HasMCData()){
       dataE[0] = track->Pt();
@@ -1342,6 +1805,33 @@ void AliAnalysisTaskHFE::ProcessAOD(){
         fQACollection->Fill("PIDperformance", dataE);
       }
     }
+
+    if (GetPlugin(kDEstep)) {
+      if (!HasMCData()){
+        Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
+        fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
+        dataDca[0]=-1; //for data, don't know the origin
+        dataDca[1]=track->Pt();
+        dataDca[2]=hfeimpactR;
+        dataDca[3]=fCentralityF;
+        dataDca[4] = -1; // not store V0 for the moment
+        dataDca[5] = double(track->Charge());
+        fQACollection->Fill("Dca", dataDca);
+      }
+
+      // Fill Containers for impact parameter analysis
+      if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
+      if(signal) {
+        // Apply weight for background contamination after ip cut
+        if(fBackGroundFactorApply) {
+              fWeightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
+              if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
+              else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
+              // weightBackGround as special weight
+              fVarManager->FillContainer(fContainer, "hadronicBackground", 2, kFALSE, fWeightBackGround);
+        }
+      }
+    }
   }
   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
 }
@@ -1373,7 +1863,9 @@ Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
     if(aodmctrack) aodmctrack->XvYvZv(vertex);
   }
 
+  //printf("MC Generated\n");
   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
+  //printf("MC Generated pass\n");
   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
   signalContainer[4] = 0;
   if(fSignalCuts->IsSelected(track)){
@@ -1421,34 +1913,14 @@ Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
   
 
   Bool_t survived = kTRUE;
-  
-  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) {
-    survived = kFALSE;
-    //printf("Did not pass AliHFEcuts::kStepRecKineITSTPC\n");
-  }
-  //else printf("Pass AliHFEcuts::kStepRecKineITSTPC\n");
-  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) {
-    survived = kFALSE;
-    //printf("Did not pass AliHFEcuts::kStepRecPrim\n");
-  }
-  //else printf("Pass AliHFEcuts::kStepRecPrim\n");
-  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) {
-    survived = kFALSE;
-    //printf("Did not pass AliHFEcuts::kStepHFEcutsITS\n");
-  }
-  //else printf("Pass AliHFEcuts::kStepHFEcutsITS\n");
-  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTOF, track)) {
-    survived = kFALSE;
-    //printf("Did not pass AliHFEcuts::kStepHFEcutsTOF\n");
-  }
-  //else printf("Pass AliHFEcuts::kStepHFEcutsTOF\n");
-  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) {
-    survived = kFALSE;
-    //printf("Did not pass AliHFEcuts::kStepHFEcutsTRD\n");
+
+  if(fCutspreselect) {
+    //printf("test preselect\n");
+    if(!fCutspreselect->IsSelected(track)) survived=kFALSE;
   }
-  //else printf("Pass AliHFEcuts::kStepHFEcutsTRD\n");
+  //printf("survived %d\n",(Int_t)survived);
   
-  if(survived){
+  if(survived && fPIDpreselect){
     // Apply PID
     AliHFEpidObject hfetrack;
     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
@@ -1507,15 +1979,15 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
   fContainer->CreateContainer("recTrackContReco", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
   fContainer->CreateContainer("recTrackContMC", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
   
-  fContainer->CreateContainer("hadronicBackground", "Container for Hadronic Background", 2);
+  fContainer->CreateContainer("hadronicBackground", "Container for Hadronic Background", 3);
   fContainer->CreateContainer("recTrackContDEReco", "Container for displaced electron analysis with Reco information", 1);
   fContainer->CreateContainer("recTrackContDEMC", "Container for displaced electron analysis with MC information", 1);
   fContainer->CreateContainer("recTrackContSecvtxReco", "Container for secondary vertexing analysis with Reco information", 1);
   fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1);
 
   if(HasMCData()){
-    fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",4);
-    fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",4);
+    fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",7);
+    fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",7);
     fContainer->Sumw2("conversionElecs");
     fContainer->Sumw2("mesonElecs");
    
@@ -1524,10 +1996,11 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
       const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
       for(Int_t iSource = 0; iSource < kElecBgSpecies; iSource++){
         for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
-          fContainer->CreateContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted conversion electrons from %s grandm., %s level",sourceName[iSource],levelName[iLevel]),4);
-          fContainer->CreateContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted electrons from %s decays, %s level",sourceName[iSource],levelName[iLevel]),4);
+          fContainer->CreateContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted conversion electrons from %s grandm., %s level",sourceName[iSource],levelName[iLevel]),5);
+          fContainer->CreateContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted electrons from %s decays, %s level",sourceName[iSource],levelName[iLevel]),5);
           fContainer->Sumw2(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
           fContainer->Sumw2(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
+          if(IsPbPb())break;
         }
       }
     }
@@ -1585,36 +2058,50 @@ void AliAnalysisTaskHFE::InitContaminationQA(){
   // 
   // Add QA for Impact Parameter cut
   //
-  const Double_t kPtbound[2] = {0.1, 20.};
-  Int_t iBin[1];
-  iBin[0] = 44; // bins in pt
-  fQACollection->CreateTH1F("hadronsBeforeIPcut", "Hadrons before IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("hadronsAfterIPcut", "Hadrons after IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", iBin[0], kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
-
-  fQACollection->CreateTH2F("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1);
-  fQACollection->CreateTH2F("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1);
-  fQACollection->CreateTH1F("Kptspectra", "Charged Kaons: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("K0Lptspectra", "K0L: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
-
-  const Double_t kDCAbound[2] = {-5., 5.};
-  const Double_t kDCAsigbound[2] = {-50., 50.};
-
-  fQACollection->CreateTH2F("dataDcaSig", "data dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("charmDcaSig", "charm dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("beautyDcaSig", "beauty dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("conversionDcaSig", "conversion dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("nonhfeDcaSig", "nonhfe dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-
-  fQACollection->CreateTH2F("dataDca", "data dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("charmDca", "charm dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("beautyDca", "beauty dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("conversionDca", "conversion dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("nonhfeDca", "nonhfe dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-
-  fQACollection->CreateTH2F("dataDcaSigDca", "data dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("conversionDcaSigDca", "conversion dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0);
+
+  TObjArray *array = fVarManager->GetVariables();
+  Int_t nvars = array->GetEntriesFast();
+  for(Int_t v = 0; v < nvars; v++) {
+    AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
+    if(!variable) continue;
+    TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
+    if(!name.CompareTo("pt")) {
+      const Int_t nBinPt  = variable->GetNumberOfBins();
+      const Double_t *kPtRange = variable->GetBinning();
+
+      fQACollection->CreateTH1Farray("hadronsBeforeIPcut", "Hadrons before IP cut", nBinPt, kPtRange);
+      fQACollection->CreateTH1Farray("hadronsAfterIPcut", "Hadrons after IP cut", nBinPt, kPtRange);
+
+      fQACollection->CreateTH2Farray("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
+      fQACollection->CreateTH2Farray("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
+      fQACollection->CreateTH1Farray("Kptspectra", "Charged Kaons: MC p_{t} ", nBinPt, kPtRange);
+      fQACollection->CreateTH1Farray("K0Lptspectra", "K0L: MC p_{t} ", nBinPt, kPtRange);
+
+      const Double_t kDCAbound[2] = {-0.2, 0.2};
+
+      const Int_t nDimDca=6;
+      const Int_t nBinDca[nDimDca] = { 8, nBinPt, 800, 12,  6, 2};
+      Double_t minimaDca[nDimDca]  = { -1., 0., kDCAbound[0], -1., -1, -1.1};
+      Double_t maximaDca[nDimDca]  = { 7., 20., kDCAbound[1], 11.,  5, 1.1};
+
+      Double_t *sourceBins = AliHFEtools::MakeLinearBinning(nBinDca[0], minimaDca[0], maximaDca[0]);
+      Double_t *dcaBins = AliHFEtools::MakeLinearBinning(nBinDca[2], minimaDca[2], maximaDca[2]);
+      Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBinDca[3], minimaDca[3], maximaDca[3]);
+      Double_t *v0PIDBins = AliHFEtools::MakeLinearBinning(nBinDca[4], minimaDca[4], maximaDca[4]);
+      Double_t *chargeBins = AliHFEtools::MakeLinearBinning(nBinDca[5], minimaDca[5], maximaDca[5]);
+
+      fQACollection->CreateTHnSparseNoLimits("Dca", "Dca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; centrality bin; v0pid; charge", nDimDca, nBinDca);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(0, sourceBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(1, kPtRange);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(2, dcaBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(3, centralityBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(4, v0PIDBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(5, chargeBins);
+
+      break;
+    }  
+  }
+
 }
 
 //____________________________________________________________
@@ -1740,6 +2227,7 @@ void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
     case kPostProcess: SETBIT(fPlugins, plug); break;
     case kDEstep: SETBIT(fPlugins, plug); break;
     case kTaggedTrackAnalysis: SETBIT(fPlugins, plug); break;
+    case kNonPhotonicElectron: SETBIT(fPlugins, plug); break; 
     default: AliError("Unknown Plugin");
   };
 }
@@ -1762,18 +2250,28 @@ Bool_t AliAnalysisTaskHFE::ReadCentrality() {
   //
   // Recover the centrality of the event from ESD or AOD
   //
+  
+  Float_t fCentralityLimitstemp[12];
+  Float_t fCentralityLimitsdefault[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
+  if(!fPbPbUserCentralityBinning) memcpy(fCentralityLimitstemp,fCentralityLimitsdefault,sizeof(fCentralityLimitsdefault));
+  else memcpy(fCentralityLimitstemp,fCentralityLimits,sizeof(fCentralityLimitsdefault));
+  
+
   Int_t bin = -1;
   if(IsPbPb()) {
     // Centrality
     AliCentrality *centrality = fInputEvent->GetCentrality();
-    Float_t fCentralityFtemp = centrality->GetCentralityPercentile("V0M");
-    Float_t centralityLimits[12] = {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
+    fCentralityPercent = centrality->GetCentralityPercentile("V0M");
+    //printf("centrality %f\n",fCentralityPercent);
+
     for(Int_t ibin = 0; ibin < 11; ibin++){
-      if(fCentralityFtemp >= centralityLimits[ibin] && fCentralityFtemp < centralityLimits[ibin+1]){
-          bin = ibin;
+      if(fCentralityPercent >= fCentralityLimitstemp[ibin] && fCentralityPercent < fCentralityLimitstemp[ibin+1]){
+        bin = ibin;
+       //printf("test bin %f, low %f, high %f, %d\n",fCentralityPercent,fCentralityLimitstemp[ibin],fCentralityLimitstemp[ibin+1],ibin);
         break;
       }
-    } 
+    }
+    
     if(bin == -1) bin = 11; // Overflow
   } else {
     // PP: Tracklet multiplicity, use common definition
@@ -1789,6 +2287,7 @@ Bool_t AliAnalysisTaskHFE::ReadCentrality() {
   }
   fCentralityF = bin;
   AliDebug(2, Form("Centrality class %d\n", fCentralityF));
+
  
   // contributors, to be outsourced
   const AliVVertex *vtx;
@@ -1813,8 +2312,12 @@ Bool_t AliAnalysisTaskHFE::ReadCentrality() {
   }
   else {
     Int_t contributorstemp = vtx->GetNContributors();
-    if( contributorstemp <=  0) fContributors =  0.5;
+    if( contributorstemp <=  0) {
+      fContributors =  0.5;
+      //printf("Number of contributors %d and vz %f\n",contributorstemp,vtx->GetZ());
+    }
     else fContributors = 1.5;
+    //printf("Number of contributors %d\n",contributorstemp);
   }
   return kTRUE;
 }
@@ -1883,8 +2386,11 @@ void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
        return;
    }
    // PileUp
+   fIdentifiedAsPileUp = kFALSE;
    if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE; 
    // Z vertex
+   fIdentifiedAsOutInz = kFALSE;
+   //printf("Z vertex %f and out %f\n",fAOD->GetPrimaryVertex()->GetZ(),fCuts->GetVertexRange());
    if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
    // Event Cut
    fPassTheEventCut = kTRUE;