]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliAnalysisTaskHFE.cxx
Updates for TRD HFE analysis
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFE.cxx
index e3fc8f181a06bed182b0243f7a84c98499479654..3ee5f60bc2cb311030593ee5a2e2636f905918d9 100644 (file)
@@ -43,7 +43,9 @@
 #include <TF1.h>
 #include <TTree.h>
 
+#include "AliESDtrackCuts.h"
 #include "AliAnalysisManager.h"
+#include "AliAnalysisUtils.h"
 #include "AliAODInputHandler.h"
 #include "AliAODMCParticle.h"
 #include "AliAODTrack.h"
 #include "AliOADBContainer.h"
 #include "AliStack.h"
 #include "AliTriggerAnalysis.h"
+#include "AliTRDTriggerAnalysis.h" 
 #include "AliVVertex.h"
 
 #include "AliHFEcollection.h"
 #include "AliHFEcontainer.h"
 #include "AliHFEcuts.h"
 #include "AliHFEelecbackground.h"
+#include "AliHFENonPhotonicElectron.h"
 #include "AliHFEmcQA.h"
 #include "AliHFEpairs.h"
 #include "AliHFEpid.h"
 #include "AliHFEpidQAmanager.h"
-#include "AliHFEpostAnalysis.h"
 #include "AliHFEsecVtxs.h"
 #include "AliHFEsecVtx.h"
 #include "AliHFEsignalCuts.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)
+  , fCollisionSystem(3)
   , fFillSignalOnly(kTRUE)
   , fFillNoCuts(kFALSE)
-  , fUseFlagAOD(kTRUE)
   , fApplyCutAOD(kFALSE)
-  , fFlags(1<<4)
   , fBackGroundFactorApply(kFALSE)
   , fRemovePileUp(kFALSE)
   , fIdentifiedAsPileUp(kFALSE)
@@ -103,14 +109,15 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fRejectKinkMother(kTRUE)
   , fisppMultiBin(kFALSE)
   , fPbPbUserCentralityBinning(kFALSE)
+  , fRemoveFirstEvent(kFALSE)
   , fisNonHFEsystematics(kFALSE)
   , fSpecialTrigger(NULL)
   , fCentralityF(-1)
   , fCentralityPercent(-1)
+  , fCentralityEstimator("V0M")
   , fContributors(0.5)
   , fWeightBackGround(0.)
   , fVz(0.0)
-  , fHadronBackgroundOADB(NULL)
   , fContainer(NULL)
   , fVarManager(NULL)
   , fSignalCuts(NULL)
@@ -118,17 +125,22 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fTriggerAnalysis(NULL)
   , fPID(NULL)
   , fPIDqa(NULL)
+  , fTRDTriggerAnalysis(NULL)
   , fPIDpreselect(NULL)
   , fCuts(NULL)
   , fTaggedTrackCuts(NULL)
   , fCleanTaggedTrack(kFALSE)
   , fVariablesTRDTaggedTrack(kFALSE)
+  , fAnalysisUtils(NULL)
   , fCutspreselect(NULL)
   , fSecVtx(NULL)
   , fElecBackGround(NULL)
   , fMCQA(NULL)
   , fTaggedTrackAnalysis(NULL)
   , fExtraCuts(NULL)
+  , fBackgroundSubtraction(NULL)
+  , fTRDTrigger(kFALSE)
+  , fWhichTRDTrigger(0)
   , fQA(NULL)
   , fOutput(NULL)
   , fHistMCQA(NULL)
@@ -143,21 +155,22 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
-  memset(&fPbPbUserCentralityBinning, 0, sizeof(Float_t) * 12);
-
+  memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
 
+  SetppAnalysis();
 }
 
 //____________________________________________________________
 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   AliAnalysisTaskSE(name)
+  , fAODMCHeader(NULL)
+  , fAODArrayMCInfo(NULL)
   , fQAlevel(0)
   , fPlugins(0)
+  , fCollisionSystem(3)
   , fFillSignalOnly(kTRUE)
   , fFillNoCuts(kFALSE)
-  , fUseFlagAOD(kTRUE)
   , fApplyCutAOD(kFALSE)
-  , fFlags(1<<4)
   , fBackGroundFactorApply(kFALSE)
   , fRemovePileUp(kFALSE)
   , fIdentifiedAsPileUp(kFALSE)
@@ -166,14 +179,15 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fRejectKinkMother(kTRUE)
   , fisppMultiBin(kFALSE)
   , fPbPbUserCentralityBinning(kFALSE)
+  , fRemoveFirstEvent(kFALSE)
   , fisNonHFEsystematics(kFALSE)
   , fSpecialTrigger(NULL)
   , fCentralityF(-1)
   , fCentralityPercent(-1)
+  , fCentralityEstimator("V0M")
   , fContributors(0.5)
   , fWeightBackGround(0.)
   , fVz(0.0)
-  , fHadronBackgroundOADB(NULL)
   , fContainer(NULL)
   , fVarManager(NULL)
   , fSignalCuts(NULL)
@@ -181,17 +195,22 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fTriggerAnalysis(NULL)
   , fPID(NULL)
   , fPIDqa(NULL)
+  , fTRDTriggerAnalysis(NULL)
   , fPIDpreselect(NULL)
   , fCuts(NULL)
   , fTaggedTrackCuts(NULL)
   , fCleanTaggedTrack(kFALSE)
   , fVariablesTRDTaggedTrack(kFALSE)
+  , fAnalysisUtils(NULL)
   , fCutspreselect(NULL)
   , fSecVtx(NULL)
   , fElecBackGround(NULL)
   , fMCQA(NULL)
   , fTaggedTrackAnalysis(NULL)
   , fExtraCuts(NULL)
+  , fBackgroundSubtraction(NULL)
+  , fTRDTrigger(kFALSE)
+  , fWhichTRDTrigger(0)
   , fQA(NULL)
   , fOutput(NULL)
   , fHistMCQA(NULL)
@@ -208,25 +227,29 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   fPID = new AliHFEpid("hfePid");
   fPIDqa = new AliHFEpidQAmanager;
   fVarManager = new AliHFEvarManager("hfeVarManager");
+  fAnalysisUtils = new AliAnalysisUtils;
+  fTRDTriggerAnalysis = new AliTRDTriggerAnalysis();
 
   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
-  memset(&fPbPbUserCentralityBinning, 0, sizeof(Float_t) * 12);
+  memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
 
+  SetppAnalysis();
 }
 
 //____________________________________________________________
 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   AliAnalysisTaskSE(ref)
+  , fAODMCHeader(NULL)
+  , fAODArrayMCInfo(NULL)
   , fQAlevel(0)
   , fPlugins(0)
+  , fCollisionSystem(ref.fCollisionSystem)
   , fFillSignalOnly(ref.fFillSignalOnly)
   , fFillNoCuts(ref.fFillNoCuts)
-  , fUseFlagAOD(ref.fUseFlagAOD)
   , fApplyCutAOD(ref.fApplyCutAOD)
-  , fFlags(ref.fFlags) 
   , fBackGroundFactorApply(ref.fBackGroundFactorApply)
   , fRemovePileUp(ref.fRemovePileUp)
   , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
@@ -235,14 +258,15 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fRejectKinkMother(ref.fRejectKinkMother)
   , fisppMultiBin(ref.fisppMultiBin)
   , fPbPbUserCentralityBinning(ref.fPbPbUserCentralityBinning)
+  , fRemoveFirstEvent(ref.fRemoveFirstEvent)
   , fisNonHFEsystematics(ref.fisNonHFEsystematics)
   , fSpecialTrigger(ref.fSpecialTrigger)
   , fCentralityF(ref.fCentralityF)
   , fCentralityPercent(ref.fCentralityPercent)
+  , fCentralityEstimator(ref.fCentralityEstimator)
   , fContributors(ref.fContributors)
   , fWeightBackGround(ref.fWeightBackGround)
   , fVz(ref.fVz)
-  , fHadronBackgroundOADB(ref.fHadronBackgroundOADB)
   , fContainer(NULL)
   , fVarManager(NULL)
   , fSignalCuts(NULL)
@@ -250,17 +274,22 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fTriggerAnalysis(NULL)
   , fPID(NULL)
   , fPIDqa(NULL)
+  , fTRDTriggerAnalysis(NULL)
   , fPIDpreselect(NULL)
   , fCuts(NULL)
   , fTaggedTrackCuts(NULL)
   , fCleanTaggedTrack(ref.fCleanTaggedTrack)
   , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
+  , fAnalysisUtils(NULL)
   , fCutspreselect(NULL)
   , fSecVtx(NULL)
   , fElecBackGround(NULL)
   , fMCQA(NULL)
   , fTaggedTrackAnalysis(NULL)
   , fExtraCuts(NULL)
+  , fBackgroundSubtraction(NULL)
+  , fTRDTrigger(ref.fTRDTrigger)
+  , fWhichTRDTrigger(ref.fWhichTRDTrigger)
   , fQA(NULL)
   , fOutput(NULL)
   , fHistMCQA(NULL)
@@ -290,13 +319,14 @@ 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.fCollisionSystem = fCollisionSystem;
   target.fFillSignalOnly = fFillSignalOnly;
   target.fFillNoCuts = fFillNoCuts;
-  target.fUseFlagAOD = fUseFlagAOD;
   target.fApplyCutAOD = fApplyCutAOD;
-  target.fFlags = fFlags;
   target.fBackGroundFactorApply = fBackGroundFactorApply;
   target.fRemovePileUp = fRemovePileUp;
   target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
@@ -305,14 +335,15 @@ void AliAnalysisTaskHFE::Copy(TObject &o) const {
   target.fRejectKinkMother = fRejectKinkMother;
   target.fisppMultiBin =   fisppMultiBin;
   target.fPbPbUserCentralityBinning = fPbPbUserCentralityBinning;
+  target.fRemoveFirstEvent = fRemoveFirstEvent;
   target.fisNonHFEsystematics = fisNonHFEsystematics;
   target.fSpecialTrigger = fSpecialTrigger;
   target.fCentralityF = fCentralityF;
   target.fCentralityPercent = fCentralityPercent;
+  target.fCentralityEstimator = fCentralityEstimator;
   target.fContributors = fContributors;
   target.fWeightBackGround = fWeightBackGround;
   target.fVz = fVz;
-  target.fHadronBackgroundOADB = fHadronBackgroundOADB;
   target.fContainer = fContainer;
   target.fVarManager = fVarManager;
   target.fSignalCuts = fSignalCuts;
@@ -320,17 +351,22 @@ void AliAnalysisTaskHFE::Copy(TObject &o) const {
   target.fTriggerAnalysis = fTriggerAnalysis;
   target.fPID = fPID;
   target.fPIDqa = fPIDqa;
+  target.fTRDTriggerAnalysis = fTRDTriggerAnalysis;
   target.fPIDpreselect = fPIDpreselect;
   target.fCuts = fCuts;
   target.fTaggedTrackCuts = fTaggedTrackCuts;
   target.fCleanTaggedTrack = fCleanTaggedTrack;
   target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
+  target.fAnalysisUtils = fAnalysisUtils;
   target.fCutspreselect = fCutspreselect;
   target.fSecVtx = fSecVtx;
   target.fElecBackGround = fElecBackGround;
   target.fMCQA = fMCQA;
   target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
   target.fExtraCuts = fExtraCuts;
+  target.fBackgroundSubtraction = fBackgroundSubtraction;
+  target.fTRDTrigger = fTRDTrigger;
+  target.fWhichTRDTrigger = fWhichTRDTrigger;
   target.fQA = fQA;
   target.fOutput = fOutput;
   target.fHistMCQA = fHistMCQA;
@@ -347,13 +383,16 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
   if(fPID) delete fPID;
   if(fPIDpreselect) delete fPIDpreselect;
   if(fVarManager) delete fVarManager;
+  if(fTRDTriggerAnalysis) delete fTRDTriggerAnalysis;
   if(fCFM) delete fCFM;
   if(fTriggerAnalysis) delete fTriggerAnalysis;
   if(fSignalCuts) delete fSignalCuts;
   if(fSecVtx) delete fSecVtx;
   if(fMCQA) delete fMCQA;
   if(fElecBackGround) delete fElecBackGround;
+  if(fBackgroundSubtraction) delete fBackgroundSubtraction;
   if(fSpecialTrigger) delete fSpecialTrigger;
+  if(fAnalysisUtils) delete fAnalysisUtils;
   // Delete output objects only if we are not running in PROOF mode because otherwise this produces a crash during merging
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
   if(mgr && mgr->GetAnalysisType() != AliAnalysisManager::kProofAnalysis){
@@ -376,6 +415,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")){
@@ -393,24 +439,15 @@ 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);
   fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
   fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
-  fQACollection->CreateTH2F("TPCdEdxBeforePID", "TPC dE/dx; p (GeV/c); TPC dE/dx (a.u.)", 1000, 0., 10., 200, 0., 200.); 
-  fQACollection->CreateTH2F("TPCnSigmaBeforePID", "TPC dE/dx; p (GeV/c); TPC dE/dx - <TPC dE/dx>|_{el} (#sigma)", 1000, 0., 10., 1000, -10., 10.);
+  fQACollection->CreateTH1F("nTriggerBit", "Histo Trigger Bit", 22, 0, 22);
  
-  InitPIDperformanceQA();
-  InitContaminationQA();
   InitHistoITScluster();
+  InitContaminationQA();
   fQA->Add(fQACollection);
 
   // Initialize PID
@@ -423,6 +460,16 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   }
   fPID->SortDetectors();
 
+  // Background subtraction-------------------------------------------------------------------
+  if (GetPlugin(kNonPhotonicElectron)) {
+    if(!fBackgroundSubtraction) fBackgroundSubtraction = new AliHFENonPhotonicElectron();
+    if(IsAODanalysis()) fBackgroundSubtraction->SetAOD(kTRUE);
+    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;
@@ -475,7 +522,7 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
     fSecVtx->CreateHistograms(fHistSECVTX);
     fOutput->Add(fHistSECVTX);
   }
-  
+
   // background----------------------------------
   if (GetPlugin(kIsElecBackGround)) {
     AliInfo("Electron BackGround Analysis on");
@@ -526,6 +573,8 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
     fQA->Add(fTaggedTrackAnalysis->GetQAcollection());
   }
 
+  //fQA->Print();
+
   PrintStatus();
   // Done!!!
   PostData(1, fOutput);
@@ -537,20 +586,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()){
@@ -558,15 +613,17 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
     fPID->InitializePID(fInputEvent->GetRunNumber());
   }
 
-  // Initialize hadronic background from OADB Container
-  AliDebug(2, Form("Apply background factors: %s, OADB Container %p", fBackGroundFactorApply ? "Yes" : "No", fHadronBackgroundOADB));
-  if(fBackGroundFactorApply && !TestBit(kBackgroundInitialized)){ 
-    AliDebug(2, "Initializing Background from OADB");
-    if(!InitializeHadronBackground(fInputEvent->GetRunNumber())) AliError("Failed initializing hadronic background parameterization from OADB");
-    else AliDebug(2, "Successfully loaded Background from OADB");
-    SetBit(kBackgroundInitialized); 
+  if(fRemoveFirstEvent){
+    if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
   }
 
+  AliESDEvent *ev = dynamic_cast<AliESDEvent *>(fInputEvent);
+  if(ev && fTRDTrigger)
+  {
+      if(!CheckTRDTrigger(ev)) return;
+  }
+
+
   if(IsESDanalysis() && HasMCData()){
     // Protect against missing MC trees
     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
@@ -577,48 +634,84 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
     if(!mcH->InitOk()) return;
     if(!mcH->TreeK()) return;
     if(!mcH->TreeTR()) return;
+
+    // Background subtraction-------------------------------------------------------------------
+    if(GetPlugin(kNonPhotonicElectron)) 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("test\n");
+    //printf("test4\n");
     ProcessAOD();
   } else {
     const char *specialTrigger = GetSpecialTrigger(fInputEvent->GetRunNumber());
     // Check Trigger selection
     if(specialTrigger){
       AliDebug(2, Form("Special Trigger requested: %s", specialTrigger));
-      AliESDEvent *ev = dynamic_cast<AliESDEvent *>(fInputEvent);
       if(!(ev && ev->IsTriggerClassFired(specialTrigger))){
         AliDebug(2, "Event not selected"); 
         return;
       } else AliDebug(2, "Event Selected");
     } else AliDebug(2, "No Special Trigger requested");
-
+    
     ProcessESD();
   }
   // Done!!!
@@ -631,49 +724,8 @@ void AliAnalysisTaskHFE::Terminate(Option_t *){
   //
   // Terminate not implemented at the moment
   //
-  if(GetPlugin(kPostProcess)){
-    fOutput = dynamic_cast<TList *>(GetOutputData(1));
-    fQA = dynamic_cast<TList *>(GetOutputData(2));
-    if(!fOutput){
-      AliError("Results not available");
-      return;
-    }
-    if(!fQA){
-      AliError("QA output not available");
-      return;
-    }
-    fContainer = dynamic_cast<AliHFEcontainer *>(fOutput->FindObject("trackContainer")); 
-    if(!fContainer){
-      AliError("Track container not found");
-      return;
-    }
-    AliHFEpostAnalysis postanalysis;
-    postanalysis.SetTaskResults(fContainer);
-    TList *qalist = dynamic_cast<TList *>(fQA->FindObject("list_TaskQA"));
-    if(!qalist){
-      AliError("QA List not found");
-      return;
-    }
-    postanalysis.SetTaskQA(qalist);
-    printf("Running post analysis\n");
-    //if(HasMCData())
-    postanalysis.DrawMCSignal2Background();
-    postanalysis.DrawEfficiency();
-    postanalysis.DrawPIDperformance();
-    postanalysis.DrawCutEfficiency();
-
-    if (GetPlugin(kIsElecBackGround)) {
-      AliHFEelecbackground elecBackGround;
-      TList *oe = 0x0;
-      if(!(oe = (TList*)dynamic_cast<TList *>(fOutput->FindObject("HFEelecbackground")))){
-        return;
-      }
-      elecBackGround.Load(oe);
-      elecBackGround.Plot();
-      elecBackGround.PostProcess();      
-    }
-  }
 }
+
 //_______________________________________________________________
 Bool_t AliAnalysisTaskHFE::IsEventInBinZero() {
   //
@@ -706,8 +758,9 @@ void AliAnalysisTaskHFE::ProcessMC(){
   // In case MC QA is on also MC QA loop is done
   //
   AliDebug(3, "Processing MC Information");
-  Double_t eventContainer [4];
-  eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
+  Double_t eventContainer [4] = {0., 0., 0., 0.};
+  if(IsESDanalysis()) eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
+  else eventContainer[0] = fAODMCHeader->GetVtxZ();
   eventContainer[2] = fCentralityF;
   eventContainer[3] = fContributors;
   fVz = eventContainer[0];
@@ -724,7 +777,7 @@ void AliAnalysisTaskHFE::ProcessMC(){
         fMCQA->SetMCEvent(fMCEvent);
         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
         fMCQA->SetCentrality(fCentralityF);
-        fMCQA->SetPercentrality(fCentralityPercent);
+        fMCQA->SetPercentrality(static_cast<Int_t>(fCentralityPercent));
         if(IsPbPb()) { fMCQA->SetPbPb();}
         else
         {
@@ -761,9 +814,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++;
   }
@@ -778,6 +845,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){
@@ -789,7 +857,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
   if(fTaggedTrackAnalysis) {
     fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
     fTaggedTrackAnalysis->SetCentrality(fCentralityF);
-    if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
+    if(IsHeavyIon()) fTaggedTrackAnalysis->SetPbPb();
     else fTaggedTrackAnalysis->SetPP();
   }
 
@@ -798,7 +866,8 @@ void AliAnalysisTaskHFE::ProcessESD(){
   eventContainer[0] = 0.; 
   if(HasMCData()) eventContainer[0] = fVz;
   else {
-    if(fESD->GetPrimaryVertexSPD()) eventContainer[0] = fESD->GetPrimaryVertexSPD()->GetZ();
+    const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
+    if(vtxESD) eventContainer[0] = vtxESD->GetZ();
   }
   eventContainer[1] = 0.;
   eventContainer[2] = fCentralityF;
@@ -826,7 +895,6 @@ void AliAnalysisTaskHFE::ProcessESD(){
   if(!fPassTheEventCut) return;
   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
 
 
   fContainer->NewEvent();
 
@@ -851,33 +919,43 @@ void AliAnalysisTaskHFE::ProcessESD(){
   Double_t container[10];
   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[8]; // [source, pT, dca, dcaSig, centrality]
-  Double_t eDca[6]; // [source, pT, dca, dcaSig, centrality]
+  Double_t dataDca[6]; // [source, pT, dca, centrality]
   Int_t nElectronCandidates = 0;
   AliESDtrack *track = NULL, *htrack = NULL;
   AliMCParticle *mctrack = NULL;
   AliMCParticle *mctrackmother = NULL;
-  Int_t pid = 0;
 
   Bool_t signal = kTRUE;
 
   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");
-  }   
+  }
   fExtraCuts->SetRecEventInfo(fESD);
 
   // Electron background analysis 
   if (GetPlugin(kIsElecBackGround)) {
-    
+
     AliDebug(2, "Running BackGround Analysis");
-    
+
     fElecBackGround->Reset();
-    
+
   } // end of electron background analysis
+
+
+  // Background subtraction-------------------------------------------------------------------
+  if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->FillPoolAssociatedTracks(fInputEvent, fCentralityF);
+  //------------------------------------------------------------------------------------------
+
   //
   // Loop ESD
   //
@@ -886,7 +964,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
     AliDebug(4, "New ESD track");
     track = fESD->GetTrack(itrack);
     track->SetESDEvent(fESD);
-
+    
     // fill counts of v0-identified particles
     Int_t v0pid = -1;
     if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
@@ -896,14 +974,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;
     }
 
@@ -931,7 +1013,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
         fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
       }
     }
-
+  
     // RecKine: ITSTPC cuts  
     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
     
@@ -980,7 +1062,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
          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){
@@ -1000,6 +1082,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
                  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;
              }
@@ -1020,28 +1103,23 @@ void AliAnalysisTaskHFE::ProcessESD(){
         }
       }
 
-      // check if it is the proton from the lambda decay, and yes fill the dca info
-      Double_t hfeimpactR4p=0., hfeimpactnsigmaR4p=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, hfeimpactR4p, hfeimpactnsigmaR4p);
-          fQACollection->Fill("pionDcaSig",track->Pt(),hfeimpactnsigmaR4p);
-        }
-      }
-      else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 2212)){
-        fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4p, hfeimpactnsigmaR4p);
-        fQACollection->Fill("protonDcaSig",track->Pt(),hfeimpactnsigmaR4p);
-        Int_t glabel=TMath::Abs(mctrack->GetMother());
-        if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
-          if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==3122){
-            fQACollection->Fill("secondaryprotonDcaSig",track->Pt(),hfeimpactnsigmaR4p);
-          }
+          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, hfeimpactR4p, hfeimpactnsigmaR4p);
-          Int_t sourceDca =-1;
+          fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
           if(fSignalCuts->IsCharmElectron(track)){
             sourceDca=1;
           }
@@ -1057,56 +1135,56 @@ void AliAnalysisTaskHFE::ProcessESD(){
           else if(fSignalCuts->IsJpsiElectron(track)){
             sourceDca=5;
           }
-          else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
-            sourceDca=6;
-            fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
-            fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
-          }
           else {
-            sourceDca=7;
+            sourceDca=6;
           }
-          eDca[0]=sourceDca;
-          eDca[1]=track->Pt();
-          eDca[2]=hfeimpactR4p;
-          eDca[3]=hfeimpactnsigmaR4p;
-          eDca[4]=fCentralityF;
-          eDca[5] = track->Charge()/3;
-
-          fQACollection->Fill("eDca", eDca);
+          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){
-      if(track->GetInnerParam())
-        fQACollection->Fill("TPCdEdxBeforePID", track->GetInnerParam()->P(), track->GetTPCsignal());
-      fQACollection->Fill("TPCnSigmaBeforePID", track->P(), fInputHandler->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron));
-    }
-
     AliHFEpidObject hfetrack;
     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
     hfetrack.SetRecTrack(track);
     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
     hfetrack.SetCentrality(fCentralityF);
-    if(IsPbPb()) hfetrack.SetPbPb();
+    hfetrack.SetMulitplicity(ncontribVtx);
+    if(IsHeavyIon()) 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()) {
-            TBits shared = track->GetTPCSharedMap();
+      TBits shared = track->GetTPCSharedMap();
            Int_t sharebit=0;
-            if(shared.CountBits() >= 2) sharebit=1;
+      if(shared.CountBits() >= 2) sharebit=1;
 
            Double_t itschi2percluster = 0.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};
-            fQACollection->Fill("fChi2perITScluster", itsChi2);
+      Double_t itsChi2[7] = {track->Pt(),track->Eta(), track->Phi(),
+                       static_cast<Double_t>(fCentralityF),static_cast<Double_t>(track->GetTPCsignalN()), static_cast<Double_t>(sharebit),itschi2percluster};
+      fQACollection->Fill("fChi2perITScluster", itsChi2);
     }
     else{
 
@@ -1114,7 +1192,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
       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};
+      Double_t itsChi2[3] = {track->Pt(), static_cast<Double_t>(fCentralityF), itschi2percluster};
       fQACollection->Fill("fChi2perITScluster", itsChi2);
     }
 
@@ -1137,12 +1215,12 @@ void AliAnalysisTaskHFE::ProcessESD(){
     // Fill Containers
     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(fBackGroundFactorApply) {
+        if(IsHeavyIon() && 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;
+        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);
       }
@@ -1159,44 +1237,6 @@ void AliAnalysisTaskHFE::ProcessESD(){
       }
     }
 
-    if(HasMCData()){
-      dataE[0] = track->Pt();
-      dataE[1] = track->Eta();
-      dataE[2] = track->Phi();
-      dataE[3] = track->Charge();
-      dataE[4] = -1.;
-      dataE[5] = -1.;
-
-      // Track selected: distinguish between true and fake
-      AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
-      if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
-        Int_t type = 0;
-        if(fSignalCuts->IsCharmElectron(track))
-          type = 1;
-        else if(fSignalCuts->IsBeautyElectron(track))
-          type = 2;
-        AliDebug(1, Form("Type: %d\n", type));
-        if(type){
-                dataE[5] = type; // beauty[1] or charm[2]
-                dataE[4] = 2;  // signal electron
-        }
-        else{
-                dataE[4] = 1; // not a signal electron
-                dataE[5] = 0;
-        }
-      } 
-      else {
-        // Fill THnSparse with the information for Fake Electrons
-        dataE[4] = 0;
-        dataE[5] = 0;
-      }
-      // fill the performance THnSparse, if the mc origin could be defined
-      if(dataE[4] > -1){
-        AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
-        fQACollection->Fill("PIDperformance", dataE);
-      }
-    }
-
     // Electron background analysis 
     if (GetPlugin(kIsElecBackGround)) {
       
@@ -1209,47 +1249,22 @@ void AliAnalysisTaskHFE::ProcessESD(){
       }
     } // end of electron background analysis
 
-
-    Int_t sourceDca =-1;
     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);
-      sourceDca=0;
       if(HasMCData())
       {
-       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 if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
-         sourceDca=6;
-          fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
-          fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
-        }
-       else {
-         sourceDca=7;
-       }
-
+        if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
+            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){
@@ -1264,6 +1279,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;
               }
@@ -1288,33 +1304,33 @@ void AliAnalysisTaskHFE::ProcessESD(){
         }
       } // end of MC
 
-      dataDca[0]=sourceDca;
+      dataDca[0]=-1; //for data, don't know the origin
       dataDca[1]=track->Pt();
       dataDca[2]=hfeimpactR;
-      dataDca[3]=hfeimpactnsigmaR;
-      dataDca[4]=fCentralityF;
-      dataDca[5] = 49;
-      Double_t xr[3]={49,49,49};
-      if(HasMCData()) {
-        mctrack->XvYvZv(xr);
-        dataDca[5] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
-      } 
-      dataDca[6] = v0pid;
-      dataDca[7] = track->Charge()/3;
-
- //     printf("Entries dca: [%.3f|%.3f|%.3f|%f|%f]\n", dataDca[0],dataDca[1],dataDca[2],dataDca[3],dataDca[4]);
+      dataDca[3]=fCentralityF;
+      dataDca[4] = v0pid;
+      dataDca[5] = double(track->Charge());
       if (!HasMCData()) fQACollection->Fill("Dca", dataDca);
-      else if(signal) 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 && 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             
@@ -1328,6 +1344,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;
               }
@@ -1355,12 +1372,16 @@ void AliAnalysisTaskHFE::ProcessESD(){
       if(HasMCData()){
         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
           fQACollection->Fill("hadronsAfterIPcut",track->Pt());
-          fQACollection->Fill("hadronsAfterIPcutMC",mctrack->Pt());
         }
       }
     }
 
   }
+
+  // Background subtraction-------------------------------------------------------------------
+  if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->CountPoolAssociated(fInputEvent, fCentralityF);
+  //------------------------------------------------------------------------------------------
+
   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
 }
 
@@ -1373,9 +1394,10 @@ void AliAnalysisTaskHFE::ProcessAOD(){
   //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; 
@@ -1407,31 +1429,65 @@ void AliAnalysisTaskHFE::ProcessAOD(){
   //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(fUseFlagAOD){
-      if(track->GetFlags() != fFlags) 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){
@@ -1441,12 +1497,25 @@ void AliAnalysisTaskHFE::ProcessAOD(){
     }
 
     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;
       
@@ -1466,6 +1535,53 @@ void AliAnalysisTaskHFE::ProcessAOD(){
       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
@@ -1474,22 +1590,41 @@ void AliAnalysisTaskHFE::ProcessAOD(){
     hfetrack.SetRecTrack(track);
     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
     hfetrack.SetCentrality(fCentralityF);
-    if(IsPbPb()) hfetrack.SetPbPb();
+    hfetrack.SetMulitplicity(ncontribVtx); // for correction
+    if(IsHeavyIon()) 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() && mctrack)  mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
+      fBackgroundSubtraction->LookAtNonHFE(itrack, track, fInputEvent, 1, fCentralityF, -1,mcsource, indexmother);
+    }
+    //---------------------------------------------------------------------------------------------------------------------
+
+    // end AOD QA
+    fQACollection->Fill("Filterend", -1);  
+    for(Int_t k=0; k<20; k++) {
+      Int_t u = 1<<k;
+      if((track->TestFilterBit(u))) {
+             fQACollection->Fill("Filterend", k);
+      }
+    }
+       
     // Apply weight for background contamination
     //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(IsHeavyIon() && 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;
+             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);
       }
@@ -1497,44 +1632,39 @@ void AliAnalysisTaskHFE::ProcessAOD(){
     }
     
     nElectronCandidates++;    
-    if(HasMCData()){
-      dataE[0] = track->Pt();
-      dataE[1] = track->Eta();
-      dataE[2] = track->Phi();
-      dataE[3] = track->Charge();
-      dataE[4] = -1;
-      dataE[5] = -1;
-      // Track selected: distinguish between true and fake
-      AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack ? mctrack->GetPdgCode(): -1));
-      if(mctrack && ((pid = TMath::Abs(mctrack->GetPdgCode())) == 11)){
-      
-        Int_t type = 0;
-        if(fSignalCuts->IsCharmElectron(track))
-          type = 1;
-        else if(fSignalCuts->IsBeautyElectron(track))
-          type = 2;
-        AliDebug(1, Form("Type: %d\n", type));
-        if(type){
-                dataE[5] = type; // beauty[1] or charm[2]
-                dataE[4] = 2;  // signal electron
-        }
-        else{
-                dataE[4] = 1; // not a signal electron
-                dataE[5] = 0;
-        }
-      } 
-      else {
-        // Fill THnSparse with the information for Fake Electrons
-        dataE[4] = 0;
-        dataE[5] = 0;
+
+    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 the performance THnSparse, if the mc origin could be defined
-      if(dataE[4] > -1){
-        AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
-        fQACollection->Fill("PIDperformance", dataE);
+
+      // 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);
+        }
       }
     }
   }
+
+  // Background subtraction-------------------------------------------------------------------
+  if (GetPlugin(kNonPhotonicElectron)) fBackgroundSubtraction->CountPoolAssociated(fInputEvent, fCentralityF);
+  //------------------------------------------------------------------------------------------
+
   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
 }
 
@@ -1546,12 +1676,7 @@ Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
   // Works for AOD and MC analysis Type
   //
   fVarManager->NewTrack(track, NULL, fCentralityF, -1, kTRUE);
-  Double_t signalContainer[6];
 
-  signalContainer[0] = track->Pt();
-  signalContainer[1] = track->Eta();
-  signalContainer[2] = track->Phi();
-  signalContainer[3] = track->Charge()/3;
 
   Double_t vertex[3] = {0.,0.,0.}; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
   if(IsESDanalysis()){
@@ -1565,31 +1690,10 @@ 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)){
-    //fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCsignal, kFALSE);
-    // Filling of the Signal/Background histogram using the 
-    // definition of the codes for charm and beauty as below in
-    // th crearion of the histogram
-    if(fSignalCuts->IsCharmElectron(track))
-      signalContainer[4] = 1;
-    else 
-      signalContainer[4] = 2;
-  } else {
-    signalContainer[4] = 0; // (and other background)
-  }
-  signalContainer[5] = 0;
-  // apply cut on the sqrt of the production vertex
-  Double_t radVertex = TMath::Sqrt(vertex[0]*vertex[0] + vertex[1] * vertex[1]);
-  if(radVertex < 3.5){
-    // Within first ITS layer(2) -> Background we cannot reject by ITS cut, let it pass
-    signalContainer[5] = 1;
-  } else if (radVertex < 7.5){
-    signalContainer[5] = 2;
-  }
-  fQACollection->Fill("SignalToBackgroundMC", signalContainer);
 
   // Step GeneratedZOutNoPileUp
   if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0)) return kFALSE;
@@ -1613,34 +1717,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);
@@ -1699,7 +1783,7 @@ 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);
@@ -1716,10 +1800,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;
         }
       }
     }
@@ -1747,31 +1832,6 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
     fContainer->SetStepTitle("recTrackContMC", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
   }
 }
-
-//____________________________________________________________
-void AliAnalysisTaskHFE::InitPIDperformanceQA(){
-  // Add a histogram for Fake electrons
-  const Int_t nDim=6;
-  Int_t nBin[nDim] = {40, 8, 18, 2, 3, 3};
-  //number of variables on the grid:pt,eta,phi,charge,
-  const Double_t kPtbound[2] = {0.1, 20.};
-  const Double_t kEtabound[2] = {-0.8, 0.8};
-  const Double_t kPhibound[2] = {0., 2. * TMath::Pi()}; 
-  const Double_t kChargebound[2] = {-1.1, 1.1};
-  const Double_t kAddInf1bound[2] = {0., 3.};
-  const Double_t kAddInf2bound[2] = {0., 3.};
-  Double_t minima[nDim] = {kPtbound[0], kEtabound[0], kPhibound[0], kChargebound[0], kAddInf1bound[0], kAddInf2bound[0]}; 
-  Double_t maxima[nDim] = {kPtbound[1], kEtabound[1], kPhibound[1], kChargebound[1], kAddInf1bound[1], kAddInf2bound[1]}; 
-  
-  fQACollection->CreateTHnSparse("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin, minima, maxima);
-  fQACollection->CreateTHnSparse("SignalToBackgroundMC", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; flavor (0 - no, 1 - charm, 2 - bottom); ITS Cluster (0 - no, 1 - first (and maybe second), 2 - second)", nDim, nBin, minima, maxima);
-
-  fQACollection->BinLogAxis("PIDperformance", 0);
-  fQACollection->BinLogAxis("SignalToBackgroundMC", 0);
-  fQACollection->Sumw2("PIDperformance");
-  fQACollection->Sumw2("SignalToBackgroundMC");
-}
-
 //____________________________________________________________
 void AliAnalysisTaskHFE::InitContaminationQA(){
   // 
@@ -1790,53 +1850,32 @@ void AliAnalysisTaskHFE::InitContaminationQA(){
 
       fQACollection->CreateTH1Farray("hadronsBeforeIPcut", "Hadrons before IP cut", nBinPt, kPtRange);
       fQACollection->CreateTH1Farray("hadronsAfterIPcut", "Hadrons after IP cut", nBinPt, kPtRange);
-      fQACollection->CreateTH1Farray("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", nBinPt, kPtRange);
-      fQACollection->CreateTH1Farray("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", 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] = {-5., 5.};
-      const Double_t kDCAsigbound[2] = {-50., 50.};
+      const Double_t kDCAbound[2] = {-0.2, 0.2};
 
-      const Int_t nDimDca=8;
-      const Int_t nBinDca[nDimDca] = { 9, nBinPt, 2000, 2000, 12, 500, 6, 2};
-      const Int_t nDimeDca=6;
-      const Int_t nBineDca[nDimeDca] = { 9, nBinPt, 2000, 2000, 12, 2};
-      Double_t minimaDca[nDimDca]  = { -1., 0., kDCAbound[0], kDCAsigbound[0], -1., 0, -1, -1.1};
-      Double_t maximaDca[nDimDca]  = { 8., 20., kDCAbound[1], kDCAsigbound[1], 11., 50, 5, 1.1};
+      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 *dcaSigBins = AliHFEtools::MakeLinearBinning(nBinDca[3], minimaDca[3], maximaDca[3]);
-      Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBinDca[4], minimaDca[4], maximaDca[4]);
-      Double_t *eProdRBins = AliHFEtools::MakeLinearBinning(nBinDca[5], minimaDca[5], maximaDca[5]);
-      Double_t *v0PIDBins = AliHFEtools::MakeLinearBinning(nBinDca[6], minimaDca[6], maximaDca[6]);
-      Double_t *chargeBins = AliHFEtools::MakeLinearBinning(nBinDca[7], minimaDca[7], maximaDca[7]);
+      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; dcasig; centrality bin; eProdR; v0pid; charge", nDimDca, nBinDca);
+      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, dcaSigBins);
-      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(4, centralityBins);
-      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(5, eProdRBins);
-      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(6, v0PIDBins);
-      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(7, chargeBins);
-
-      fQACollection->CreateTHnSparseNoLimits("eDca", "eDca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; dcasig; centrality bin; charge", nDimeDca, nBineDca);
-      ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(0, sourceBins);
-      ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(1, kPtRange);
-      ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(2, dcaBins);
-      ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(3, dcaSigBins);
-      ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(4, centralityBins);
-      ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(5, chargeBins);
-
-      fQACollection->CreateTH2Farray("secondaryprotonDcaSig", "secondary proton dca significance: dca sig ", nBinPt, kPtRange, 2000, kDCAsigbound[0], kDCAsigbound[1]);
-      fQACollection->CreateTH2Farray("protonDcaSig", "proton dca significance: dca sig ", nBinPt, kPtRange, 2000, kDCAsigbound[0], kDCAsigbound[1]);
-      fQACollection->CreateTH2Farray("pionDcaSig", "pion dca significance: dca sig ", nBinPt, kPtRange, 2000, kDCAsigbound[0], kDCAsigbound[1]);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(3, centralityBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(4, v0PIDBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(5, chargeBins);
 
       break;
     }  
@@ -1967,6 +2006,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");
   };
 }
@@ -1997,10 +2037,10 @@ Bool_t AliAnalysisTaskHFE::ReadCentrality() {
   
 
   Int_t bin = -1;
-  if(IsPbPb()) {
+  if(IsHeavyIon()) {
     // Centrality
     AliCentrality *centrality = fInputEvent->GetCentrality();
-    fCentralityPercent = centrality->GetCentralityPercentile("V0M");
+    fCentralityPercent = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
     //printf("centrality %f\n",fCentralityPercent);
 
     for(Int_t ibin = 0; ibin < 11; ibin++){
@@ -2043,7 +2083,7 @@ Bool_t AliAnalysisTaskHFE::ReadCentrality() {
       AliError("ESD Event required for ESD Analysis");
       return kFALSE;
     }
-    vtx = fESD->GetPrimaryVertexSPD();
+    vtx = fESD->GetPrimaryVertex() ;
   }
   if(!vtx){ 
     fContributors = 0.5;
@@ -2051,7 +2091,10 @@ 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);
   }
@@ -2086,29 +2129,6 @@ Int_t AliAnalysisTaskHFE::GetITSMultiplicity(AliVEvent *ev){
   return nAcc;
 }
 
-//___________________________________________________
-Bool_t AliAnalysisTaskHFE::InitializeHadronBackground(Int_t run){
-  //
-  // Initialize background factors array from the AliOADBContainer
-  // The container is expected to provide a TObjArray with the name
-  // "hadronBackground" and the size 12 for the given run number
-  //
-  AliDebug(1, "Deriving hadronic background parameterization from OADB container");
-  TObjArray *functions = dynamic_cast<TObjArray *>(fHadronBackgroundOADB->GetObject(run, "hadronBackground"));
-  if(!functions){
-    AliDebug(1, "Content in the OADB Container is not a TObjArray");
-    fBackGroundFactorApply = kFALSE;
-    return kFALSE;
-  } 
-  if(functions->GetSize() < 12){
-    AliDebug(1, Form("Size not matching: 12 expected, %d provided", functions->GetSize()));
-    fBackGroundFactorApply = kFALSE;
-    return kFALSE;
-  }
-  for(Int_t icent = 0; icent < 12; icent++) fkBackGroundFactorArray[icent] = dynamic_cast<const TF1 *>(functions->UncheckedAt(icent));
-  return kTRUE;
-}
-
 //___________________________________________________
 void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
   //
@@ -2144,26 +2164,211 @@ void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
    // PileUp
    fIdentifiedAsPileUp = kFALSE;
    if(fRemovePileUp && fESD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE; 
+   
+
+
    // Z vertex
    fIdentifiedAsOutInz = kFALSE;
-   if (IsPbPb()) {
-     //printf("PbPb\n");
-     if(fESD->GetPrimaryVertex()){
-       if(TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
-     }
-   }
-   else {
-     //printf("pp\n");
-     if(fESD->GetPrimaryVertexTracks()){
-       if(TMath::Abs(fESD->GetPrimaryVertexTracks()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
-     }
+   Bool_t findvertex = kTRUE;
+   const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
+   if((!vtxESD) || (vtxESD->GetNContributors() <= 0)) findvertex = kFALSE;
+   if(findvertex) {
+     if(TMath::Abs(vtxESD->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
    }
+   
    //Event Cut
    fPassTheEventCut = kTRUE;
    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) fPassTheEventCut = kFALSE;   
   
-
  }
 
 }
 
+//___________________________________________________
+Bool_t AliAnalysisTaskHFE::CheckTRDTrigger(AliESDEvent *ev) {
+//
+// Check TRD trigger; pPb settings
+//
+    Bool_t cint8=kFALSE;
+    Bool_t cint7=kFALSE;
+    Bool_t cint5=kFALSE;
+    Bool_t trdtrgevent=kFALSE;
+    fTRDTriggerAnalysis->CalcTriggers(ev); 
+
+    // HSE no cleanup
+    if(fWhichTRDTrigger==1)
+    {
+       cint8= ev->IsTriggerClassFired("CINT8WUHSE-B-NOPF-CENT");
+       cint7= ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT");
+       cint5= (ev->IsTriggerClassFired("CINT5WU-B-NOPF-ALL")) &&
+           (ev->GetHeader()->GetL1TriggerInputs() & (1 << 10));
+       if((cint7==kFALSE)&&(cint8==kFALSE)&&(cint5==kFALSE)) return kFALSE;
+       else
+       {
+           DrawTRDTrigger(ev);
+           return kTRUE;
+       }
+    }
+
+    // HSE cleanup
+    if(fWhichTRDTrigger==2)
+    {
+       if(!fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHSE))
+       {
+           return kFALSE;
+       }
+       else
+       {
+           DrawTRDTrigger(ev);
+           return kTRUE;
+       }
+    }
+
+    //HQU no cleanup
+    if(fWhichTRDTrigger==3)
+    {
+       cint8= ev->IsTriggerClassFired("CINT8WUHQU-B-NOPF-CENT");
+       cint7= ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT");
+       cint5= (ev->IsTriggerClassFired("CINT5WU-B-NOPF-ALL")) &&
+           (ev->GetHeader()->GetL1TriggerInputs() & (1 << 12));
+       if((cint7==kFALSE)&&(cint8==kFALSE)&&(cint5==kFALSE)) return kFALSE;
+       else
+       {
+           DrawTRDTrigger(ev);
+           return kTRUE;
+       }
+    }
+
+    // HQU cleanup
+    if(fWhichTRDTrigger==4)
+    {
+
+       if(!fTRDTriggerAnalysis->IsFired(AliTRDTriggerAnalysis::kHQU))
+       {
+           return kFALSE;
+       }
+       else
+       {
+           DrawTRDTrigger(ev);
+           return kTRUE;
+       }
+    }
+   
+    return trdtrgevent;
+
+}
+
+void AliAnalysisTaskHFE::DrawTRDTrigger(AliESDEvent *ev) {
+
+    Int_t ntriggerbit=0;
+    fQACollection->Fill("nTriggerBit",ntriggerbit);
+    if(ev->IsTriggerClassFired("CINT7-B-NOPF-ALLNOTRD"))
+    {
+       ntriggerbit=2;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+    }
+    if(ev->IsTriggerClassFired("CINT7WU-B-NOPF-ALL"))
+    {
+       ntriggerbit=3;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+       if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
+           ntriggerbit=18;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+       }
+       if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")) {
+           ntriggerbit=19;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+       }
+    }
+    if(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-CENT"))
+    {
+       ntriggerbit=4;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+
+       if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
+           ntriggerbit=13;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+       }
+       if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")) {
+           ntriggerbit=14;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+           if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
+               ntriggerbit=17;
+               fQACollection->Fill("nTriggerBit",ntriggerbit);
+           }
+       }
+    }
+    if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT"))
+    {
+       ntriggerbit=5;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+       if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
+           ntriggerbit=11;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+       }
+       if((!(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")))&&(!(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-CENT")))) {
+           ntriggerbit=21;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+
+            /*
+           Int_t nTrdTracks = ev->GetNumberOfTrdTracks();
+           for (Int_t iTrack = 0; iTrack < nTrdTracks; ++iTrack) {
+               AliESDTrdTrack* trdTrack = ev->GetTrdTrack(iTrack);
+               printf("GTU track %3i: pt = %5.1f, PID = %3i\n", iTrack, trdTrack->Pt(), trdTrack->GetPID());
+           }*/
+
+
+       }
+
+    }
+    if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT"))
+    {
+       ntriggerbit=6;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+       if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")) {
+           ntriggerbit=12;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+       }
+       if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-FAST")){
+           ntriggerbit=15;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+
+           if((!(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-CENT")))&&(!(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-CENT")))) {
+               ntriggerbit=20;
+               fQACollection->Fill("nTriggerBit",ntriggerbit);
+               /*
+                Int_t nTrdTracks = ev->GetNumberOfTrdTracks();
+                for (Int_t iTrack = 0; iTrack < nTrdTracks; ++iTrack) {
+                AliESDTrdTrack* trdTrack = ev->GetTrdTrack(iTrack);
+                printf("HSE GTU track %3i: pt = %5.1f, PID = %3i\n", iTrack, trdTrack->Pt(), trdTrack->GetPID());
+                }                          */
+
+           }
+
+       }
+
+    }
+    if(ev->IsTriggerClassFired("CEMC7WUHEE-B-NOPF-CENT")) {
+       ntriggerbit=7;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+    }
+    if(ev->IsTriggerClassFired("CINT7WUHJT-B-NOPF-FAST")){
+       ntriggerbit=8;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+    }
+    if(ev->IsTriggerClassFired("CINT7WUHQU-B-NOPF-FAST")){
+       ntriggerbit=9;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+    }
+    if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-FAST")){
+       ntriggerbit=10;
+       fQACollection->Fill("nTriggerBit",ntriggerbit);
+       if(ev->IsTriggerClassFired("CINT7WUHSE-B-NOPF-CENT")) {
+           ntriggerbit=16;
+           fQACollection->Fill("nTriggerBit",ntriggerbit);
+       }
+    }
+    if(ntriggerbit==0) fQACollection->Fill("nTriggerBit",1);
+
+}
+