Removed hardcoded AliESDs.root in AliESDInputHandler(RP)::Notify()

[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFE.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
// The analysis task:
// Filling an AliCFContainer with the quantities pt, eta and phi
// for tracks which survivied the particle cuts (MC resp. ESD tracks)
// Track selection is done using the AliHFE package
// 
// Author:
//  Raphaelle Bailhache <R.Bailhache@gsi.de>
//  Markus Fasel <M.Fasel@gsi.de>
//  Matus Kalisky <matus.kalisky@cern.ch>
//  MinJung Kweon <minjung@physi.uni-heidelberg.de>
//
#include <TAxis.h>
#include <TBits.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TH3D.h>
#include <TIterator.h>
#include <TList.h>
#include <TLegend.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TParticle.h>
#include <TProfile.h>
#include <TString.h>
#include <TF1.h>
#include <TTree.h>

#include "AliAnalysisManager.h"
#include "AliAODInputHandler.h"
#include "AliAODMCParticle.h"
#include "AliAODTrack.h"
#include "AliAODVertex.h"
#include "AliCentrality.h"
#include "AliCFContainer.h"
#include "AliCFManager.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMCParticle.h"
#include "AliMultiplicity.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliOADBContainer.h"
#include "AliStack.h"
#include "AliTriggerAnalysis.h"
#include "AliVVertex.h"

#include "AliHFEcollection.h"
#include "AliHFEcontainer.h"
#include "AliHFEcuts.h"
#include "AliHFEelecbackground.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 "AliHFEtaggedTrackAnalysis.h"
#include "AliHFEtools.h"
#include "AliHFEvarManager.h"
#include "AliAnalysisTaskHFE.h"

ClassImp(AliAnalysisTaskHFE)

//____________________________________________________________
AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  AliAnalysisTaskSE("PID efficiency Analysis")
  , fQAlevel(0)
  , fPlugins(0)
  , fFillSignalOnly(kTRUE)
  , fFillNoCuts(kFALSE)
  , fUseFlagAOD(kTRUE)
  , fApplyCutAOD(kFALSE)
  , fFlags(1<<4)
  , fBackGroundFactorApply(kFALSE)
  , fRemovePileUp(kFALSE)
  , fIdentifiedAsPileUp(kFALSE)
  , fIdentifiedAsOutInz(kFALSE)
  , fPassTheEventCut(kFALSE)
  , fRejectKinkMother(kTRUE)
  , fisppMultiBin(kFALSE)
  , fPbPbUserCentralityBinning(kFALSE)
  , fisNonHFEsystematics(kFALSE)
  , fSpecialTrigger(NULL)
  , fCentralityF(-1)
  , fCentralityPercent(-1)
  , fContributors(0.5)
  , fWeightBackGround(0.)
  , fVz(0.0)
  , fHadronBackgroundOADB(NULL)
  , fContainer(NULL)
  , fVarManager(NULL)
  , fSignalCuts(NULL)
  , fCFM(NULL)
  , fTriggerAnalysis(NULL)
  , fPID(NULL)
  , fPIDqa(NULL)
  , fPIDpreselect(NULL)
  , fCuts(NULL)
  , fTaggedTrackCuts(NULL)
  , fCleanTaggedTrack(kFALSE)
  , fVariablesTRDTaggedTrack(kFALSE)
  , fCutspreselect(NULL)
  , fSecVtx(NULL)
  , fElecBackGround(NULL)
  , fMCQA(NULL)
  , fTaggedTrackAnalysis(NULL)
  , fExtraCuts(NULL)
  , fQA(NULL)
  , fOutput(NULL)
  , fHistMCQA(NULL)
  , fHistSECVTX(NULL)
  , fHistELECBACKGROUND(NULL)
  , fQACollection(NULL)
{
  //
  // Dummy constructor
  //
  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(fCentralityLimits, 0, sizeof(Float_t) * 12);


}

//____________________________________________________________
AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
  AliAnalysisTaskSE(name)
  , fQAlevel(0)
  , fPlugins(0)
  , fFillSignalOnly(kTRUE)
  , fFillNoCuts(kFALSE)
  , fUseFlagAOD(kTRUE)
  , fApplyCutAOD(kFALSE)
  , fFlags(1<<4)
  , fBackGroundFactorApply(kFALSE)
  , fRemovePileUp(kFALSE)
  , fIdentifiedAsPileUp(kFALSE)
  , fIdentifiedAsOutInz(kFALSE)
  , fPassTheEventCut(kFALSE)  
  , fRejectKinkMother(kTRUE)
  , fisppMultiBin(kFALSE)
  , fPbPbUserCentralityBinning(kFALSE)
  , fisNonHFEsystematics(kFALSE)
  , fSpecialTrigger(NULL)
  , fCentralityF(-1)
  , fCentralityPercent(-1)
  , fContributors(0.5)
  , fWeightBackGround(0.)
  , fVz(0.0)
  , fHadronBackgroundOADB(NULL)
  , fContainer(NULL)
  , fVarManager(NULL)
  , fSignalCuts(NULL)
  , fCFM(NULL)
  , fTriggerAnalysis(NULL)
  , fPID(NULL)
  , fPIDqa(NULL)
  , fPIDpreselect(NULL)
  , fCuts(NULL)
  , fTaggedTrackCuts(NULL)
  , fCleanTaggedTrack(kFALSE)
  , fVariablesTRDTaggedTrack(kFALSE)
  , fCutspreselect(NULL)
  , fSecVtx(NULL)
  , fElecBackGround(NULL)
  , fMCQA(NULL)
  , fTaggedTrackAnalysis(NULL)
  , fExtraCuts(NULL)
  , fQA(NULL)
  , fOutput(NULL)
  , fHistMCQA(NULL)
  , fHistSECVTX(NULL)
  , fHistELECBACKGROUND(NULL)
  , fQACollection(0x0)
{
  //
  // Default constructor
  // 
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());

  fPID = new AliHFEpid("hfePid");
  fPIDqa = new AliHFEpidQAmanager;
  fVarManager = new AliHFEvarManager("hfeVarManager");

  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(fCentralityLimits, 0, sizeof(Float_t) * 12);

}

//____________________________________________________________
AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
  AliAnalysisTaskSE(ref)
  , fQAlevel(0)
  , fPlugins(0)
  , fFillSignalOnly(ref.fFillSignalOnly)
  , fFillNoCuts(ref.fFillNoCuts)
  , fUseFlagAOD(ref.fUseFlagAOD)
  , fApplyCutAOD(ref.fApplyCutAOD)
  , fFlags(ref.fFlags) 
  , fBackGroundFactorApply(ref.fBackGroundFactorApply)
  , fRemovePileUp(ref.fRemovePileUp)
  , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
  , fIdentifiedAsOutInz(ref.fIdentifiedAsOutInz)
  , 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)
  , fHadronBackgroundOADB(ref.fHadronBackgroundOADB)
  , fContainer(NULL)
  , fVarManager(NULL)
  , fSignalCuts(NULL)
  , fCFM(NULL)
  , fTriggerAnalysis(NULL)
  , fPID(NULL)
  , fPIDqa(NULL)
  , fPIDpreselect(NULL)
  , fCuts(NULL)
  , fTaggedTrackCuts(NULL)
  , fCleanTaggedTrack(ref.fCleanTaggedTrack)
  , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
  , fCutspreselect(NULL)
  , fSecVtx(NULL)
  , fElecBackGround(NULL)
  , fMCQA(NULL)
  , fTaggedTrackAnalysis(NULL)
  , fExtraCuts(NULL)
  , fQA(NULL)
  , fOutput(NULL)
  , fHistMCQA(NULL)
  , fHistSECVTX(NULL)
  , fHistELECBACKGROUND(NULL)
  , fQACollection(NULL)
{
  //
  // Copy Constructor
  //
  ref.Copy(*this);
}

//____________________________________________________________
AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
  //
  // Assignment operator
  //
  if(this == &ref) 
    ref.Copy(*this);
  return *this;
}

//____________________________________________________________
void AliAnalysisTaskHFE::Copy(TObject &o) const {
  // 
  // Copy into object o
  //
  AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
  target.fQAlevel = fQAlevel;
  target.fPlugins = fPlugins;
  target.fFillSignalOnly = fFillSignalOnly;
  target.fFillNoCuts = fFillNoCuts;
  target.fUseFlagAOD = fUseFlagAOD;
  target.fApplyCutAOD = fApplyCutAOD;
  target.fFlags = fFlags;
  target.fBackGroundFactorApply = fBackGroundFactorApply;
  target.fRemovePileUp = fRemovePileUp;
  target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
  target.fIdentifiedAsOutInz = fIdentifiedAsOutInz;
  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;
  target.fHadronBackgroundOADB = fHadronBackgroundOADB;
  target.fContainer = fContainer;
  target.fVarManager = fVarManager;
  target.fSignalCuts = fSignalCuts;
  target.fCFM = fCFM;
  target.fTriggerAnalysis = fTriggerAnalysis;
  target.fPID = fPID;
  target.fPIDqa = fPIDqa;
  target.fPIDpreselect = fPIDpreselect;
  target.fCuts = fCuts;
  target.fTaggedTrackCuts = fTaggedTrackCuts;
  target.fCleanTaggedTrack = fCleanTaggedTrack;
  target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
  target.fCutspreselect = fCutspreselect;
  target.fSecVtx = fSecVtx;
  target.fElecBackGround = fElecBackGround;
  target.fMCQA = fMCQA;
  target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
  target.fExtraCuts = fExtraCuts;
  target.fQA = fQA;
  target.fOutput = fOutput;
  target.fHistMCQA = fHistMCQA;
  target.fHistSECVTX = fHistSECVTX;
  target.fHistELECBACKGROUND = fHistELECBACKGROUND;
  target.fQACollection = fQACollection;
}

//____________________________________________________________
AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
  //
  // Destructor
  //
  if(fPID) delete fPID;
  if(fPIDpreselect) delete fPIDpreselect;
  if(fVarManager) delete fVarManager;
  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(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();
  if(mgr && mgr->GetAnalysisType() != AliAnalysisManager::kProofAnalysis){
    if(fPIDqa) delete fPIDqa;
    if(fOutput) delete fOutput;
    if(fQA) delete fQA;
  }
}

//____________________________________________________________
void AliAnalysisTaskHFE::UserCreateOutputObjects(){
  //
  // Creating output container and output objects
  // Here we also Initialize the correction framework container and 
  // the objects for
  // - PID
  // - MC QA
  // - SecVtx
  // QA histograms are created if requested
  // Called once per worker
  //
  AliDebug(3, "Creating Output Objects");
  // Automatic determination of the analysis mode
  AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
    SetAODAnalysis();
  } else {
    SetESDAnalysis();
    if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
      SetHasMCData();
  }
  printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
  printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");

  // Enable Trigger Analysis
  fTriggerAnalysis = new AliTriggerAnalysis;
  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.);
 
  InitPIDperformanceQA();
  InitContaminationQA();
  InitHistoITScluster();
  fQA->Add(fQACollection);

  // Initialize PID
  fPID->SetHasMCData(HasMCData());
  if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
  if(IsQAOn(kPIDqa)){
    AliInfo("PID QA switched on");
    fPIDqa->Initialize(fPID);
    fQA->Add(fPIDqa->MakeList("HFEpidQA"));
  }
  fPID->SortDetectors();

  // Initialize correction Framework and Cuts
  const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack + AliHFEcuts::kNcutStepsSecvtxTrack;
  fCFM = new AliCFManager;
  fCFM->SetNStepParticle(kNcutSteps);
  MakeParticleContainer();
  MakeEventContainer();
  // Temporary fix: Initialize particle cuts with NULL
  for(Int_t istep = 0; istep < kNcutSteps; istep++)
    fCFM->SetParticleCutsList(istep, NULL);
  if(!fCuts){
    AliWarning("Cuts not available. Default cuts will be used");
    fCuts = new AliHFEcuts;
    fCuts->CreateStandardCuts();
  }
  if(IsAODanalysis()) fCuts->SetAOD();
  // Make clone for V0 tagging step
  fCuts->Initialize(fCFM);
  if(fCuts->IsQAOn()) fQA->Add(fCuts->GetQAhistograms());
  fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
  fVarManager->SetSignalCuts(fSignalCuts);
 
  // add output objects to the List
  fOutput->AddAt(fContainer, 0);
  fOutput->AddAt(fCFM->GetEventContainer(), 1);
  
  // mcQA----------------------------------
  if (HasMCData() && IsQAOn(kMCqa)) {
    AliInfo("MC QA on");
    if(!fMCQA) fMCQA = new AliHFEmcQA;
    if(!fHistMCQA) fHistMCQA = new TList();
    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);
  } 

  // secvtx----------------------------------
  if (GetPlugin(kSecVtx)) {
    AliInfo("Secondary Vertex Analysis on");
    if(!fSecVtx) fSecVtx = new AliHFEsecVtx;
    fSecVtx->SetHasMCData(HasMCData());

    if(!fHistSECVTX) fHistSECVTX = new TList();
    fHistSECVTX->SetOwner();
    fSecVtx->CreateHistograms(fHistSECVTX);
    fOutput->Add(fHistSECVTX);
  }
  
  // background----------------------------------
  if (GetPlugin(kIsElecBackGround)) {
    AliInfo("Electron BackGround Analysis on");
    if(!fElecBackGround){
      AliWarning("ElecBackGround not available. Default elecbackground will be used");
      fElecBackGround = new AliHFEelecbackground;
    }
    fElecBackGround->SetHasMCData(HasMCData());

    if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
    fHistELECBACKGROUND->SetOwner();
    fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
    fOutput->Add(fHistELECBACKGROUND);
  }  

  // tagged tracks
  if(GetPlugin(kTaggedTrackAnalysis)){
    AliInfo("Analysis on V0-tagged tracks enabled");
    fTaggedTrackAnalysis = new AliHFEtaggedTrackAnalysis(Form("taggedTrackAnalysis%s", GetName()));
    fTaggedTrackAnalysis->SetCuts(fTaggedTrackCuts);
    fTaggedTrackAnalysis->SetClean(fCleanTaggedTrack);
    AliHFEvarManager *varManager = fTaggedTrackAnalysis->GetVarManager();
    TObjArray *array = fVarManager->GetVariables();
    Int_t nvars = array->GetEntriesFast();
    TString namee;
    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("source")) namee = TString("species");
      else namee = TString(name);
      Int_t nbins = variable->GetNumberOfBins();
      if(variable->HasUserDefinedBinning()){
        varManager->AddVariable(namee, nbins, variable->GetBinning());
      } else {
        varManager->AddVariable(namee, nbins, variable->GetMinimum(), variable->GetMaximum());
      }
      //printf("For AliTaggedTrackAnalysis, had the variable %s and the one used %s\n",(const char*)variable->GetName(),(const char*) namee);
    }
    if(fPIDqa->HasHighResolutionHistos()) 
      fTaggedTrackAnalysis->GetPIDqa()->SetHighResolutionHistos();
    fTaggedTrackAnalysis->SetPID(fPID);
    fTaggedTrackAnalysis->SetVariablesTRD(fVariablesTRDTaggedTrack);
    fTaggedTrackAnalysis->InitContainer();
    fOutput->Add(fTaggedTrackAnalysis->GetContainer());
    fQA->Add(fTaggedTrackAnalysis->GetPIDQA());
    fQA->Add(fTaggedTrackAnalysis->GetCutQA());
    fQA->Add(fTaggedTrackAnalysis->GetQAcollection());
  }

  PrintStatus();
  // Done!!!
  PostData(1, fOutput);
  PostData(2, fQA);
}

//____________________________________________________________
void AliAnalysisTaskHFE::UserExec(Option_t *){
  //
  // Run the analysis
  // 
  AliDebug(3, "Starting Single Event Analysis");
  if(!fInputEvent){
    AliError("Reconstructed Event not available");
    return;
  }
  if(HasMCData()){
    AliDebug(4, Form("MC Event: %p", fMCEvent));
    if(!fMCEvent){
      AliError("No MC Event, but MC Data required");
      return;
    }
  }
  if(!fCuts){
    AliError("HFE cuts not available");
    return;
  }
  if(!fPID->IsInitialized()){
    // Initialize PID with the given run number
    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(IsESDanalysis() && HasMCData()){
    // Protect against missing MC trees
    AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if(!mcH){ 
      AliError("No MC Event Handler available");
      return;
    }
    if(!mcH->InitOk()) return;
    if(!mcH->TreeK()) return;
    if(!mcH->TreeTR()) return;
  }

  // need the centrality for everything (MC also)
  fCentralityF = -1;
  if(!ReadCentrality()) fCentralityF = -1;
  //printf("pass centrality\n");
  //printf("Reading fCentralityF %f\n",fCentralityF);
  
  // See if pile up and z in the range
  RejectionPileUpVertexRangeEventCut();

  // 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());
  }
  fPID->SetPIDResponse(pidResponse);
  if(fPIDpreselect) fPIDpreselect->SetPIDResponse(pidResponse);

  // Event loop
  if(IsAODanalysis()){
    //printf("test\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!!!
  PostData(1, fOutput);
  PostData(2, fQA);
}

//____________________________________________________________
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() {
  //
  //
  //

  //printf("test in IsEventInBinZero\n");
  if(!fInputEvent){
    AliError("Reconstructed Event not available");
    return kFALSE;
  }

  // check vertex
  const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
  if(!vertex) return kTRUE;
  //if(vertex) return kTRUE;

  // check tracks
  if(fInputEvent->GetNumberOfTracks()<=0) return kTRUE;
  //if(fInputEvent->GetNumberOfTracks()>0) return kTRUE;
  
  
  return kFALSE;
  
}
//____________________________________________________________
void AliAnalysisTaskHFE::ProcessMC(){
  //
  // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
  // 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();
  eventContainer[2] = fCentralityF;
  eventContainer[3] = fContributors;
  fVz = eventContainer[0];
  //printf("z position is %f\n",eventContainer[0]);
  //if(fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent)) 
  fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
  Int_t nElectrons = 0;
  if(IsESDanalysis()){
   if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
    if (HasMCData() && IsQAOn(kMCqa)) {
      AliDebug(2, "Running MC QA");

      if(fMCEvent->Stack()){
        fMCQA->SetMCEvent(fMCEvent);
        fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
        fMCQA->SetCentrality(fCentralityF);
        fMCQA->SetPercentrality(fCentralityPercent);
        if(IsPbPb()) { fMCQA->SetPbPb();}
        else
        {
            if(fisppMultiBin) fMCQA->SetPPMultiBin();
            else fMCQA->SetPP();
        }
        fMCQA->Init();

        fMCQA->GetMesonKine();

        // loop over all tracks for decayed electrons
        for (Int_t igen = 0; igen < fMCEvent->GetNumberOfTracks(); igen++){
          TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
          if(!mcpart) continue;
          fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
          fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
          fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
          fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
          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);
      }

    } // end of MC QA loop
   }
   // -----------------------------------------------------------------
   fCFM->SetMCEventInfo(fMCEvent);
   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
  } else {
    fCFM->SetMCEventInfo(fInputEvent);
  }
  // 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;
    AliDebug(4, "Next MC Track");
    if(ProcessMCtrack(mctrack)) nElectrons++;
  }

  // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
  fQACollection->Fill("nElectron", nElectrons);
}

//____________________________________________________________
void AliAnalysisTaskHFE::ProcessESD(){
  //
  // Run Analysis of reconstructed event in ESD Mode
  // Loop over Tracks, filter according cut steps defined in AliHFEcuts
  //
  AliDebug(3, "Processing ESD Event");
  AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
  if(!fESD){
    AliError("ESD Event required for ESD Analysis");
    return;
  }

  // Set magnetic field if V0 task on
  if(fTaggedTrackAnalysis) {
    fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
    fTaggedTrackAnalysis->SetCentrality(fCentralityF);
    if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
    else fTaggedTrackAnalysis->SetPP();
  }

  // Do event Normalization
  Double_t eventContainer[4];
  eventContainer[0] = 0.; 
  if(HasMCData()) eventContainer[0] = fVz;
  else {
    if(fESD->GetPrimaryVertexSPD()) eventContainer[0] = fESD->GetPrimaryVertexSPD()->GetZ();
  }
  eventContainer[1] = 0.;
  eventContainer[2] = fCentralityF;
  eventContainer[3] = fContributors;
  if(fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0AND))
    eventContainer[1] = 1.;

  //
  fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);

  //
  if(fIdentifiedAsPileUp) return; 
  fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);

  //
  if(TMath::Abs(fCentralityF) < 0) return; 
  fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
  //printf("In ProcessESD %f\n",fCentralityF);

  //
  if(fIdentifiedAsOutInz) return;
  fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);  

  //
  if(!fPassTheEventCut) return;
  fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);

 

  fContainer->NewEvent();

  if (GetPlugin(kIsElecBackGround)) { 
    fElecBackGround->SetEvent(fESD);
  }
  if (GetPlugin(kSecVtx)) {
    fSecVtx->SetEvent(fESD);
    fSecVtx->GetPrimaryCondition();
  }

  if(HasMCData()){
    if (GetPlugin(kSecVtx)) { 
      fSecVtx->SetMCEvent(fMCEvent);
      fSecVtx->SetMCQA(fMCQA); 
    }
    if (GetPlugin(kIsElecBackGround)) { 
      fElecBackGround->SetMCEvent(fMCEvent);
    }
  }

  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]
  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);

  // 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
  //
  // Loop ESD
  //
  AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
  for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
    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;
    else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
    else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
    // here the tagged track analysis will run
    if(fTaggedTrackAnalysis && v0pid > -1){ 
      AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
      fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
    }
 
    AliDebug(3, Form("Doing track %d, %p", itrack, track));
     
    //////////////////////////////////////
    // preselect
    /////////////////////////////////////
    if(fPIDpreselect && fCutspreselect) {
      if(!PreSelectTrack(track)) continue;
    }

    signal = kTRUE;
    
    // Fill step without any cut
          
    if(HasMCData()){
      // Check if it is electrons near the vertex
      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;

      if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE; 
      else AliDebug(3, "Signal Electron");

      // Fill K pt for Ke3 contributions
      if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==321)) fQACollection->Fill("Kptspectra",mctrack->Pt());
      else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==130)) fQACollection->Fill("K0Lptspectra",mctrack->Pt());
    } 
    // Cache new Track information inside the var manager
    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);
      }
    }

    // RecKine: ITSTPC cuts  
    if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
    
    
    // RecPrim
    if(fRejectKinkMother) { 
      if(track->GetKinkIndex(0) != 0) continue; } // Quick and dirty fix to reject both kink mothers and daughters
    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()){
      FillProductionVertex(track);

      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){
           //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){
               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]);
                 }
               }
               break;
             }
             iName++;
             if(iName == kElecBgSpecies)iName = 0;
           }
         }
         //else{
           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]);
	   }
           //}
        }
      }

      // check if it is the proton from the lambda decay, and yes fill the dca info
      Double_t hfeimpactR4p=0., hfeimpactnsigmaR4p=0.;
      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);
          }
        }
      }
      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;
          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;
          }
          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);
        }
      }
    }

    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();
    else hfetrack.SetPP();
    fPID->SetVarManager(fVarManager);
    if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
    nElectronCandidates++;

    // Temporary histogram for chi2/ITS cluster
    if(IsPbPb()) {
            TBits shared = track->GetTPCSharedMap();
	    Int_t sharebit=0;
            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);
    }
    else{

      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[3] = {track->Pt(), fCentralityF, itschi2percluster};
      fQACollection->Fill("fChi2perITScluster", itsChi2);
    }

    // Fill Histogram for Hadronic Background
    if(HasMCData()){
      if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11))
        fVarManager->FillContainer(fContainer, "hadronicBackground", UInt_t(0), kFALSE);
      else if(mctrack){
        // Fill Ke3 contributions
        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",mctrack->Pt(),mctrackmother->Pt());
          else if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==130)
            fQACollection->Fill("Ke3K0Lecorr",mctrack->Pt(),mctrackmother->Pt());
        }
      }
    }

    // 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(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"));
    }

    Bool_t bTagged=kFALSE;
    if(GetPlugin(kSecVtx)) {
      AliDebug(2, "Running Secondary Vertex Analysis");
      if(fSecVtx->Process(track) && signal) {
        fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
        fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
        bTagged=kTRUE;
      }
    }

    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)) {
      
      AliDebug(2, "Running BackGround Analysis");
      
      for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
        htrack = fESD->GetTrack(jtrack);
        if ( itrack == jtrack ) continue;  
        fElecBackGround->PairAnalysis(track, htrack); 
      }
    } // 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(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){
            //Fill additional containers for electron source distinction           
            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){
                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]);
                }
                break;
              }
              iName++;
              if(iName == kElecBgSpecies)iName = 0;
            }
          }
          //else{
          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]=sourceDca;
      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]);
      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(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){
            //Fill additional containers for electron source distinction             
            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){
                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]);
                }
                break;
              }
              iName++;
              if(iName == kElecBgSpecies)iName = 0;
            }
          }
          // else{
            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]);
            }
            //}
        }
      }
      if(signal) {
        fVarManager->FillContainer(fContainer, "recTrackContDEReco", AliHFEcuts::kStepHFEcutsDca, kFALSE);
        fVarManager->FillContainer(fContainer, "recTrackContDEMC", AliHFEcuts::kStepHFEcutsDca, kTRUE);
        fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterDE"));
      }
      if(HasMCData()){
        if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
          fQACollection->Fill("hadronsAfterIPcut",track->Pt());
          fQACollection->Fill("hadronsAfterIPcutMC",mctrack->Pt());
        }
      }
    }

  }
  fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
}

//____________________________________________________________
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];
  if(HasMCData()) eventContainer[0] = fVz;
  else {
    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);

  //
  if(fIdentifiedAsOutInz) return;
  fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);  

  //
  if(!fPassTheEventCut) return;
  fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
  //printf("pass\n");

  fContainer->NewEvent();
 
  AliAODTrack *track = NULL;
  AliAODMCParticle *mctrack = NULL;
  Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
  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(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
    }
    fVarManager->NewTrack(track, mctrack, fCentralityF, -1, kTRUE);
    
    if(fFillNoCuts) {
      if(signal || !fFillSignalOnly){
        fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
        fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
      }
    }

    if(fApplyCutAOD) {
      // RecKine: ITSTPC cuts  
      if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) 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"));
    }

    //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);
    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

    
    // 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(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();
      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;
      }
      // 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);
      }
    }
  }
  fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
}

//____________________________________________________________
Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
  //
  // Filter the Monte Carlo Track
  // Additionally Fill a THnSparse for Signal To Background Studies
  // 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()){
    AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
    if(mctrack){
      vertex[0] = mctrack->Particle()->Vx();
      vertex[1] = mctrack->Particle()->Vy();
    }
  } else {
    AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
    if(aodmctrack) aodmctrack->XvYvZv(vertex);
  }

  if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
  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;
  fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
  //printf("In ProcessMCtrack %f\n",fCentralityF);

  // Step Generated Event Cut
  if(!fPassTheEventCut) return kFALSE;
  fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedEventCut, kFALSE);

  if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
  fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCInAcceptance, kFALSE);
  return kTRUE;
}

//____________________________________________________________
Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
  //
  // Preselect tracks
  //
  

  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");
  }
  //else printf("Pass AliHFEcuts::kStepHFEcutsTRD\n");
  
  if(survived){
    // Apply PID
    AliHFEpidObject hfetrack;
    hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
    hfetrack.SetRecTrack(track);
    if(!fPIDpreselect->IsSelected(&hfetrack)) {
      //printf("Did not pass AliHFEcuts::kPID\n");
      survived = kFALSE;
    }
    //else printf("Pass AliHFEcuts::kPID\n");
  }

  return survived; 
      
}
//____________________________________________________________
void AliAnalysisTaskHFE::MakeEventContainer(){
  //
  // Create the event container for the correction framework and link it
  // 1st bin: Vertex z-position
  // 2nd bin: V0AND decision (normalization to sigma_inel)
  // 3rd bin: Centrality class (for pp defined as number of contributors in vertex.)
  // 4th bin: Number of contributors > 0
  //
  
  const Int_t kNvar = 4;  // number of variables on the grid: 
  Int_t nBins[kNvar] = {120, 2, 11, 2};
  Double_t binMin[kNvar] = {-30. , 0., 0.0, 0.};
  Double_t binMax[kNvar] = {30., 2., 11.0, 2.};
  
  AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
  
  Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
  Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
  Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
  Double_t *contributorsBins = AliHFEtools::MakeLinearBinning(nBins[3], binMin[3], binMax[3]);
  evCont->SetBinLimits(0, vertexBins);
  evCont->SetBinLimits(1, v0andBins);
  evCont->SetBinLimits(2, centralityBins);
  evCont->SetBinLimits(3, contributorsBins);
  delete[] vertexBins; delete[] v0andBins; delete[] centralityBins; delete[] contributorsBins;
    
  fCFM->SetEventContainer(evCont);
}

//____________________________________________________________
void AliAnalysisTaskHFE::MakeParticleContainer(){
  //
  // Create the particle container for the correction framework manager and 
  // link it
  //
  if(!fContainer) fContainer = new AliHFEcontainer("trackContainer");
  fVarManager->DefineVariables(fContainer);

  // Create Correction Framework containers
  fContainer->CreateContainer("MCTrackCont", "Track Container filled with MC information", AliHFEcuts::kNcutStepsMCTrack);
  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("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",7);
    fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",7);
    fContainer->Sumw2("conversionElecs");
    fContainer->Sumw2("mesonElecs");
   
    if(fisNonHFEsystematics){
      const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
      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->Sumw2(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
          fContainer->Sumw2(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
        }
      }
    }
    //fContainer->CreateContainer("charmElecs", "Container for weighted charm electrons",2);
  }

  fContainer->CreateCorrelationMatrix("correlationstepafterPID","THnSparse with correlations");
  fContainer->CreateCorrelationMatrix("correlationstepafterDE","THnSparse with correlations");
  if(!fVarManager->IsVariableDefined("centrality")) {
    //printf("Create the two other correlation maps\n");
    fContainer->CreateCorrelationMatrix("correlationstepbeforePID","THnSparse with correlations");
    fContainer->CreateCorrelationMatrix("correlationstepafterTOF","THnSparse with correlations");
  }

  // Define the step names
  for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsMCTrack; istep++){
    fContainer->SetStepTitle("MCTrackCont", AliHFEcuts::MCCutName(istep), istep);
  }
  for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsRecTrack; istep++){
    fContainer->SetStepTitle("recTrackContReco", AliHFEcuts::RecoCutName(istep), istep);
    fContainer->SetStepTitle("recTrackContMC", AliHFEcuts::RecoCutName(istep), istep);
  }
  for(UInt_t ipid = 0; ipid < fPID->GetNumberOfPIDdetectors(); ipid++){
    fContainer->SetStepTitle("recTrackContReco", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
    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(){
  // 
  // Add QA for Impact Parameter cut
  //

  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->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 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};

      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]);

      fQACollection->CreateTHnSparseNoLimits("Dca", "Dca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; dcasig; centrality bin; eProdR; 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]);

      break;
    }  
  }

}

//____________________________________________________________
void AliAnalysisTaskHFE::InitHistoITScluster(){
  //
    // Initialize a temporary histogram to monitor the chi2/ITS cluster
    if(IsPbPb()) {
        const Int_t kNDim = 7;
        const Int_t kNBins[kNDim] = {88, 20,90,11, 160, 2, 1000};
        const Double_t kMin[kNDim] = {0.1, -1,0,  0.,0., 0,  0.};
        const Double_t kMax[kNDim] = {20., 1, 2.*TMath::Pi(), 11.,160, 2, 100.};
        fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c);eta;phi; centrality class;nclus;sharebit; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
        fQACollection->BinLogAxis("fChi2perITScluster", 0);
    }
    else
    {
        const Int_t kNDim = 3;
        const Int_t kNBins[kNDim] = {44, 11, 1000};
        const Double_t kMin[kNDim] = {0.1, 0., 0.};
        const Double_t kMax[kNDim] = {20., 11., 100.};
        fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c); centrality class; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
        fQACollection->BinLogAxis("fChi2perITScluster", 0);
    }
}

//____________________________________________________________
void AliAnalysisTaskHFE::SelectSpecialTrigger(const Char_t *trgclust, Int_t runMin, Int_t runMax){
  //
  // Select only events triggered by a special trigeer cluster
  //
  if(!fSpecialTrigger) fSpecialTrigger = new AliOADBContainer("SpecialTrigger");
  fSpecialTrigger->AppendObject(new TObjString(trgclust), runMin, runMax);
}

//____________________________________________________________
const Char_t * AliAnalysisTaskHFE::GetSpecialTrigger(Int_t run){
  //
  // Derive selected trigger string for given run
  //
  if(!fSpecialTrigger) return NULL;
  TObjString *trg = dynamic_cast<TObjString *>(fSpecialTrigger->GetObject(run));
  if(!trg) return NULL;
  return trg->String().Data();
}

//____________________________________________________________
void AliAnalysisTaskHFE::PrintStatus() const {
  //
  // Print Analysis status
  //
  printf("\n\tAnalysis Settings\n\t========================================\n\n");
  printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
  printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
  printf("\tDisplaced electron analysis step: %s\n", GetPlugin(kDEstep) ? "YES" : "NO");
  printf("\tTagged Track Analysis: %s\n", GetPlugin(kTaggedTrackAnalysis) ? "YES" : "NO");
  printf("\n");
  printf("\tParticle Identification Detectors:\n");
  fPID->PrintStatus();
  printf("\n");
  printf("\tQA: \n");
  printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
  printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsQAOn()) ? "YES" : "NO");
  printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
  printf("\n");
}

//____________________________________________________________
Bool_t AliAnalysisTaskHFE::FillProductionVertex(const AliVParticle * const track) const{
  //
  // Find the production vertex of the associated MC track
  //
  if(!fMCEvent) return kFALSE;
  const AliVParticle *mctrack = NULL;
  TString objectType = track->IsA()->GetName();
  if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
    // Reconstructed track
    mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
  } else {
    // MCParticle
    mctrack = track;
  }

  if(!mctrack) return kFALSE;

  Double_t xv = 0.0;
  Double_t yv = 0.0;
 
  if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
    // case MCParticle
    const AliMCParticle *mcpart = dynamic_cast<const AliMCParticle *>(mctrack);
    if(mcpart){
      xv =  mcpart->Xv();
      yv =  mcpart->Yv();
    }
  } else {
    // case AODMCParticle
    const AliAODMCParticle *mcpart = dynamic_cast<const AliAODMCParticle *>(mctrack);
    if(mcpart){
      xv =  mcpart->Xv();
      yv =  mcpart->Yv();
    }
  }

  //printf("xv %f, yv %f\n",xv,yv);
  fQACollection->Fill("radius", TMath::Abs(xv),TMath::Abs(yv));

  return kTRUE;

}
//__________________________________________
void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
  //
  // Switch on Plugin
  // Available:
  //  - Primary vertex studies
  //  - Secondary vertex Studies
  //  - Post Processing
  //
  switch(plug){
    case kPriVtx: SETBIT(fPlugins, plug); break;
    case kSecVtx: SETBIT(fPlugins, plug); break;
    case kIsElecBackGround: SETBIT(fPlugins, plug); break;
    case kPostProcess: SETBIT(fPlugins, plug); break;
    case kDEstep: SETBIT(fPlugins, plug); break;
    case kTaggedTrackAnalysis: SETBIT(fPlugins, plug); break;
    default: AliError("Unknown Plugin");
  };
}
//__________________________________________
Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track){
  //
  // Check single track cuts for a given cut step
  // Fill the particle container
  //
  const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
  if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
  if(fVarManager->IsSignalTrack()) {
    fVarManager->FillContainer(fContainer, "recTrackContReco", cutStep, kFALSE);
    fVarManager->FillContainer(fContainer, "recTrackContMC", cutStep, kTRUE);
  }
  return kTRUE;
}
//___________________________________________________
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();
    fCentralityPercent = centrality->GetCentralityPercentile("V0M");
    //printf("centrality %f\n",fCentralityPercent);

    for(Int_t ibin = 0; ibin < 11; 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
    Int_t itsMultiplicity = GetITSMultiplicity(fInputEvent);
    Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
    for(Int_t ibin = 0; ibin < 7; ibin++){  
      if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){
        bin = ibin;
        break;
      }
    }
    if(bin == -1) bin = 7;  // Overflow
  }
  fCentralityF = bin;
  AliDebug(2, Form("Centrality class %d\n", fCentralityF));

 
  // contributors, to be outsourced
  const AliVVertex *vtx;
  if(IsAODanalysis()){
    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
    if(!fAOD){
      AliError("AOD Event required for AOD Analysis");
      return kFALSE;
    }
    vtx = fAOD->GetPrimaryVertex();
  } else {
    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
    if(!fESD){
      AliError("ESD Event required for ESD Analysis");
      return kFALSE;
    }
    vtx = fESD->GetPrimaryVertexSPD();
  }
  if(!vtx){ 
    fContributors = 0.5;
    return kFALSE;
  }
  else {
    Int_t contributorstemp = vtx->GetNContributors();
    if( contributorstemp <=  0) fContributors =  0.5;
    else fContributors = 1.5;
    //printf("Number of contributors %d\n",contributorstemp);
  }
  return kTRUE;
}

//___________________________________________________
Int_t AliAnalysisTaskHFE::GetITSMultiplicity(AliVEvent *ev){
  //
  // Definition of the Multiplicity according to the JPSI group (F. Kramer)
  //
  Int_t nTracklets = 0;
  Int_t nAcc = 0;
  Double_t etaRange = 1.6;

  if (ev->IsA() == AliAODEvent::Class()) {
    AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
    nTracklets = tracklets->GetNumberOfTracklets();
    for (Int_t nn = 0; nn < nTracklets; nn++) {
      Double_t theta = tracklets->GetTheta(nn);
      Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
      if (TMath::Abs(eta) < etaRange) nAcc++;
    }
  } else if (ev->IsA() == AliESDEvent::Class()) {
    nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
    for (Int_t nn = 0; nn < nTracklets; nn++) {
       Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
      if (TMath::Abs(eta) < etaRange) nAcc++;
    }
  } else return -1;

  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() {
  //
  // Recover the centrality of the event from ESD or AOD
  //
 if(IsAODanalysis()){

   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
   if(!fAOD){
     AliError("AOD Event required for AOD Analysis");
       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;
   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) fPassTheEventCut = kFALSE; 
   
   
 } else {
   
   AliDebug(3, "Processing ESD Centrality");
   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
   if(!fESD){
     AliError("ESD Event required for ESD Analysis");
       return;
   }
   // 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;
     }
   }
   //Event Cut
   fPassTheEventCut = kTRUE;
   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) fPassTheEventCut = kFALSE;   
  

 }

}