]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFENonPhotonicElectron.cxx
Fix of sigmaZ for crossing tracklets from Alex
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFENonPhotonicElectron.cxx
index e256dc4046296b7c8c33c80680aff0716d50f4c7..877f743e7ff3246df025c6e43f1ce028d97ea88e 100644 (file)
@@ -6,6 +6,7 @@
  *                                                                                     *
  *     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              *
@@ -45,6 +46,7 @@
 #include "AliESDtrack.h"
 #include "AliESDtrackCuts.h"
 #include "AliPIDResponse.h"
+#include "AliPID.h"
 
 #include "AliKFParticle.h"
 #include "AliKFVertex.h"
 #include "AliHFEpid.h"
 #include "AliHFEpidQAmanager.h"
 #include "AliHFEtools.h"
+#include "AliHFEmcQA.h"
 
 #include "AliHFENonPhotonicElectron.h"
 
 ClassImp(AliHFENonPhotonicElectron)
 //________________________________________________________________________
 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const char *name, const Char_t *title)
-  :TNamed(name, title)
-  ,fMCEvent(NULL)
-  ,fAODArrayMCInfo(NULL)
-  ,fHFEBackgroundCuts(0)
-  ,fPIDBackground(0x0)
-  ,fPIDBackgroundQA(0)
-  ,fAlgorithmMA(kTRUE)
-  ,fUseFilterAOD(kTRUE)
-  ,fFilter(-1)
-  ,fChi2OverNDFCut(3.0)
-  ,fMaxDCA(3.0)
-  ,fMaxOpeningTheta(0.02)
-  ,fMaxOpeningPhi(0.1)
-  ,fMaxOpening3D(TMath::Pi())
-  ,fMaxInvMass(0.6)
-  ,fSetMassConstraint(kFALSE)
-  ,fArraytrack(NULL)
-  ,fCounterPoolBackground(0)
-  ,fListOutput(NULL)
-  ,fMCSourceee(NULL)
-  ,fUSignee(NULL)
-  ,fLSignee(NULL)
-  ,fUSignAngle(NULL)
-  ,fLSignAngle(NULL)
+  :TNamed              (name, title)
+  ,fIsAOD              (kFALSE)
+  ,fMCEvent            (NULL)
+  ,fAODArrayMCInfo     (NULL)
+  ,fLevelBack(-1)
+  ,fHFEBackgroundCuts  (NULL)
+  ,fPIDBackground      (0x0)
+  ,fPIDBackgroundQA    (0)
+  ,fkPIDRespons                (NULL)
+  ,fPtBinning()
+  ,fEtaBinning()
+  ,fAlgorithmMA                (kTRUE)
+  ,fChi2OverNDFCut     (3.0)
+  ,fMaxDCA             (3.0)
+//  ,fMaxOpeningTheta  (0.02)
+//  ,fMaxOpeningPhi    (0.1)
+  ,fMaxOpening3D       (TMath::Pi())
+  ,fMaxInvMass         (1000)
+  ,fSetMassConstraint  (kFALSE)
+  ,fSelectCategory1tracks(kTRUE)
+  ,fSelectCategory2tracks(kFALSE)
+  ,fITSmeanShift(0.)
+  ,fArraytrack         (NULL)
+  ,fCounterPoolBackground      (0)
+  ,fnumberfound                        (0)
+  ,fListOutput         (NULL)
+  ,fAssElectron                (NULL)
+  ,fIncElectron                (NULL)
+  ,fUSign              (NULL)
+  ,fLSign              (NULL)
+  ,fUSmatches(NULL)
+  ,fLSmatches(NULL)
+  ,fHnsigmaITS(NULL)
+  ,fWeightsSource(NULL)
+  ,fIncElectronRadius(NULL)
+  ,fRecElectronRadius(NULL)
+//  ,fUSignAngle       (NULL)
+//  ,fLSignAngle       (NULL)
 {
   //
   // Constructor
   //
-  fPIDBackground = new AliHFEpid("hfePidBackground");
+  fPIDBackground   = new AliHFEpid("hfePidBackground");
   fPIDBackgroundQA = new AliHFEpidQAmanager;
 }
 
 //________________________________________________________________________
 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron()
-  :TNamed()
-  ,fMCEvent(0)
-  ,fAODArrayMCInfo(NULL)
-  ,fHFEBackgroundCuts(NULL)
-  ,fPIDBackground(0x0)
-  ,fPIDBackgroundQA(0)
-  ,fAlgorithmMA(kTRUE)
-  ,fUseFilterAOD(kTRUE)
-  ,fFilter(-1)
-  ,fChi2OverNDFCut(3.0)
-  ,fMaxDCA(3.0)
-  ,fMaxOpeningTheta(0.02)
-  ,fMaxOpeningPhi(0.1)
-  ,fMaxOpening3D(TMath::TwoPi())
-  ,fMaxInvMass(0.6)
-  ,fSetMassConstraint(kFALSE)
-  ,fArraytrack(NULL)
-  ,fCounterPoolBackground(0)
-  ,fListOutput(NULL)
-  ,fMCSourceee(NULL)
-  ,fUSignee(NULL)
-  ,fLSignee(NULL)
-  ,fUSignAngle(NULL)
-  ,fLSignAngle(NULL)
+  :TNamed              ()
+  ,fIsAOD              (kFALSE)
+  ,fMCEvent            (NULL)
+  ,fAODArrayMCInfo     (NULL)
+  ,fLevelBack(-1)
+  ,fHFEBackgroundCuts  (NULL)
+  ,fPIDBackground      (0x0)
+  ,fPIDBackgroundQA    (0)
+  ,fkPIDRespons                (NULL)
+  ,fPtBinning()
+  ,fEtaBinning()
+  ,fAlgorithmMA                (kTRUE)
+  ,fChi2OverNDFCut     (3.0)
+  ,fMaxDCA             (3.0)
+//  ,fMaxOpeningTheta  (0.02)
+//  ,fMaxOpeningPhi    (0.1)
+  ,fMaxOpening3D       (TMath::TwoPi())
+  ,fMaxInvMass         (1000)
+  ,fSetMassConstraint  (kFALSE)
+  ,fSelectCategory1tracks(kTRUE)
+  ,fSelectCategory2tracks(kFALSE)
+  ,fITSmeanShift(0.)
+  ,fArraytrack         (NULL)
+  ,fCounterPoolBackground      (0)
+  ,fnumberfound                        (0)
+  ,fListOutput         (NULL)
+  ,fAssElectron                (NULL)
+  ,fIncElectron                (NULL)
+  ,fUSign              (NULL)
+  ,fLSign              (NULL)
+  ,fUSmatches(NULL)
+  ,fLSmatches(NULL)
+  ,fHnsigmaITS(NULL)
+  ,fWeightsSource(NULL)
+  ,fIncElectronRadius(NULL)
+  ,fRecElectronRadius(NULL)
+//  ,fUSignAngle       (NULL)
+//  ,fLSignAngle       (NULL)
 {
   //
   // Constructor
@@ -128,29 +159,43 @@ AliHFENonPhotonicElectron::AliHFENonPhotonicElectron()
 //________________________________________________________________________
 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const AliHFENonPhotonicElectron &ref)
   :TNamed(ref)
-  ,fMCEvent(NULL)
-  ,fAODArrayMCInfo(NULL)
-  ,fHFEBackgroundCuts(ref.fHFEBackgroundCuts)
-  ,fPIDBackground(ref.fPIDBackground)
-  ,fPIDBackgroundQA(ref.fPIDBackgroundQA)
-  ,fAlgorithmMA(ref.fAlgorithmMA)
-  ,fUseFilterAOD(ref.fUseFilterAOD)
-  ,fFilter(ref.fFilter)
-  ,fChi2OverNDFCut(ref.fChi2OverNDFCut)
-  ,fMaxDCA(ref.fMaxDCA)
-  ,fMaxOpeningTheta(ref.fMaxOpeningTheta)
-  ,fMaxOpeningPhi(ref.fMaxOpeningPhi)
-  ,fMaxOpening3D(ref.fMaxOpening3D)
-  ,fMaxInvMass(ref.fMaxInvMass)
-  ,fSetMassConstraint(ref.fSetMassConstraint)
-  ,fArraytrack(NULL)
-  ,fCounterPoolBackground(0)
-  ,fListOutput(ref.fListOutput)
-  ,fMCSourceee(ref.fMCSourceee)
-  ,fUSignee(ref.fUSignee)
-  ,fLSignee(ref.fLSignee)
-  ,fUSignAngle(ref.fUSignAngle)
-  ,fLSignAngle(ref.fLSignAngle)
+  ,fIsAOD              (ref.fIsAOD)
+  ,fMCEvent            (NULL)
+  ,fAODArrayMCInfo     (NULL)
+  ,fLevelBack           (ref.fLevelBack)
+  ,fHFEBackgroundCuts  (ref.fHFEBackgroundCuts)
+  ,fPIDBackground      (ref.fPIDBackground)
+  ,fPIDBackgroundQA    (ref.fPIDBackgroundQA)
+  ,fkPIDRespons                (ref.fkPIDRespons)
+  ,fPtBinning(ref.fPtBinning)
+  ,fEtaBinning(ref.fEtaBinning)
+  ,fAlgorithmMA                (ref.fAlgorithmMA)
+  ,fChi2OverNDFCut     (ref.fChi2OverNDFCut)
+  ,fMaxDCA             (ref.fMaxDCA)
+//  ,fMaxOpeningTheta  (ref.fMaxOpeningTheta)
+//  ,fMaxOpeningPhi    (ref.fMaxOpeningPhi)
+  ,fMaxOpening3D       (ref.fMaxOpening3D)
+  ,fMaxInvMass         (ref.fMaxInvMass)
+  ,fSetMassConstraint  (ref.fSetMassConstraint)
+  ,fSelectCategory1tracks(ref.fSelectCategory1tracks)
+  ,fSelectCategory2tracks(ref.fSelectCategory2tracks)
+  ,fITSmeanShift(ref.fITSmeanShift)
+  ,fArraytrack         (NULL)
+  ,fCounterPoolBackground      (0)
+  ,fnumberfound                        (0)
+  ,fListOutput         (ref.fListOutput)
+  ,fAssElectron                (ref.fAssElectron)
+  ,fIncElectron                (ref.fIncElectron)
+  ,fUSign              (ref.fUSign)
+  ,fLSign              (ref.fLSign)
+  ,fUSmatches(ref.fUSmatches)
+  ,fLSmatches(ref.fLSmatches)
+  ,fHnsigmaITS(ref.fHnsigmaITS)
+  ,fWeightsSource(ref.fWeightsSource)
+  ,fIncElectronRadius(ref.fIncElectronRadius)
+  ,fRecElectronRadius(ref.fRecElectronRadius)
+//  ,fUSignAngle       (ref.fUSignAngle)
+//  ,fLSignAngle       (ref.fLSignAngle)
 {
   //
   // Copy Constructor
@@ -173,9 +218,10 @@ AliHFENonPhotonicElectron::~AliHFENonPhotonicElectron()
   //
   // Destructor
   //
-  if(fArraytrack) delete fArraytrack;
-  if(fHFEBackgroundCuts) delete fHFEBackgroundCuts;
-  if(fPIDBackground) delete fPIDBackground;
+  if(fArraytrack)              delete fArraytrack;
+  //if(fHFEBackgroundCuts)     delete fHFEBackgroundCuts;
+  if(fPIDBackground)           delete fPIDBackground;
+  if(fPIDBackgroundQA)         delete fPIDBackgroundQA;
 }
 
 //_____________________________________________________________________________________________
@@ -185,23 +231,47 @@ void AliHFENonPhotonicElectron::Init()
   // Init
   //
 
+  //printf("Analysis Mode for AliHFENonPhotonicElectron: %s Analysis\n", fIsAOD ? "AOD" : "ESD");
+
   if(!fListOutput) fListOutput = new TList;
   fListOutput->SetName("HFENonPhotonicElectron");
   fListOutput->SetOwner();
 
-  if(!fHFEBackgroundCuts) fHFEBackgroundCuts = new AliESDtrackCuts();
+  if(!fHFEBackgroundCuts) fHFEBackgroundCuts = new AliHFEcuts();
+  if(fIsAOD) fHFEBackgroundCuts->SetAOD();
+  fHFEBackgroundCuts->Initialize();
+  if(fHFEBackgroundCuts->IsQAOn()) {
+    fListOutput->Add(fHFEBackgroundCuts->GetQAhistograms());
+  }
 
   // Initialize PID
   if(!fPIDBackground) fPIDBackground = new AliHFEpid("default pid");
-  if(fMCEvent || fAODArrayMCInfo) fPIDBackground->SetHasMCData(kTRUE);
-  if(!fPIDBackground->GetNumberOfPIDdetectors()) fPIDBackground->AddDetector("TPC", 0);
+  if(fMCEvent || fAODArrayMCInfo) fPIDBackground->SetHasMCData(kTRUE); // does nothing since the fMCEvent are set afterwards at the moment
+  if(!fPIDBackground->GetNumberOfPIDdetectors())
+  {
+    //fPIDBackground->AddDetector("TOF", 0);
+    fPIDBackground->AddDetector("TPC", 0);
+  }
   AliInfo("PID Background QA switched on");
   fPIDBackgroundQA->Initialize(fPIDBackground);
   fListOutput->Add(fPIDBackgroundQA->MakeList("HFENP_PID_Background"));
   fPIDBackground->SortDetectors();
 
-  Int_t nBinsPt = 24;
-  Double_t binLimPt[25] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 3., 3.5, 4., 5., 6.};
+  const Int_t kBinsPtDefault = 35;
+  Double_t binLimPtDefault[kBinsPtDefault+1] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
+  const Int_t kBinsEtaInclusiveDefault = 8;
+  Double_t binLimEtaInclusiveDefault[kBinsEtaInclusiveDefault+1] = {-0.8, -0.6, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8};
+  const Int_t kBinsEtaAssociated = 15;
+  Double_t binLimEtaAssociat[kBinsEtaAssociated+1] = {-1.5,-1.3,-1.1,-0.9,-0.7,-0.5,-0.3,-0.1,0.1,0.3,0.5,0.7,0.9,1.1,1.3,1.5};
+
+  if(!fPtBinning.GetSize()) fPtBinning.Set(kBinsPtDefault+1, binLimPtDefault);
+  if(!fEtaBinning.GetSize()) fEtaBinning.Set(kBinsEtaInclusiveDefault+1, binLimEtaInclusiveDefault);
+
+  //Int_t nBinsP = 400;
+  //Double_t minP = 0.0;
+  //Double_t maxP = 20.0;
+  //Double_t binLimP[nBinsP+1];
+  //for(Int_t i=0; i<=nBinsP; i++) binLimP[i]=(Double_t)minP + (maxP-minP)/nBinsP*(Double_t)i ;
 
   Int_t nBinsC = 11;
   Double_t minC = 0.0;
@@ -215,13 +285,13 @@ void AliHFENonPhotonicElectron::Init()
   Double_t binLimSource[nBinsSource+1];
   for(Int_t i=0; i<=nBinsSource; i++) binLimSource[i]=(Double_t)minSource + (maxSource-minSource)/nBinsSource*(Double_t)i ;
 
-  Int_t nBinsInvMass = 60;
+  Int_t nBinsInvMass = 30;
   Double_t minInvMass = 0.;
-  Double_t maxInvMass = 0.6;
+  Double_t maxInvMass = 0.3;
   Double_t binLimInvMass[nBinsInvMass+1];
   for(Int_t i=0; i<=nBinsInvMass; i++) binLimInvMass[i]=(Double_t)minInvMass + (maxInvMass-minInvMass)/nBinsInvMass*(Double_t)i ;
 
-  Int_t nBinsPhi = 180;
+  Int_t nBinsPhi = 8;
   Double_t minPhi = 0.0;
   Double_t maxPhi = TMath::Pi();
   Double_t binLimPhi[nBinsPhi+1];
@@ -231,66 +301,104 @@ void AliHFENonPhotonicElectron::Init()
     AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
   }
 
-  Int_t nBinsAngle = 180;
+  Int_t nBinsAngle = 72;
   Double_t minAngle = 0.0;
-  Double_t maxAngle = TMath::Pi();
+  Double_t maxAngle = 0.4;
   Double_t binLimAngle[nBinsAngle+1];
   for(Int_t i=0; i<=nBinsAngle; i++)
   {
     binLimAngle[i]=(Double_t)minAngle + (maxAngle-minAngle)/nBinsAngle*(Double_t)i ;
-    AliDebug(2,Form("bin phi is %f for %d",binLimPhi[i],i));
+    AliDebug(2,Form("bin phi is %f for %d",binLimAngle[i],i));
   }
 
-  const Int_t nDimMCSource=3;
-  Int_t nBinMCSource[nDimMCSource] = {nBinsC,nBinsPt,nBinsSource};
-  fMCSourceee = new THnSparseF("fMCSourceee","fMCSourceee",nDimMCSource,nBinMCSource);
-  fMCSourceee->SetBinEdges(0,binLimC);
-  fMCSourceee->SetBinEdges(1,binLimPt);
-  fMCSourceee->SetBinEdges(2,binLimSource);
-/*  fMCSourceee->GetAxis(0)->SetTitle("Multiplicity (arb.unit)");
-  fMCSourceee->GetAxis(1)->SetTitle("#it{p}_{mc,T} (GeV/c)");
-  fMCSourceee->GetAxis(2)->SetTitle("MC");*/
-  fMCSourceee->Sumw2();
-  AliDebug(2,"AliHFENonPhotonicElectron: fMCSourceee");
+  // Constrain histograms
+  const Int_t nDimSingle=4;
+  const Int_t nDimPair=9;
+  Int_t nBinPair[nDimPair] = {nBinsPhi,nBinsC,fPtBinning.GetSize()-1,nBinsInvMass,nBinsSource,nBinsAngle,fPtBinning.GetSize()-1,fEtaBinning.GetSize()-1,kBinsEtaAssociated};
+  
+  // Associated Electron
+  Int_t nBinAssElectron[nDimSingle] = {nBinsC,fPtBinning.GetSize()-1,nBinsSource,kBinsEtaAssociated};
+  fAssElectron = new THnSparseF("fAssElectron","fAssElectron",nDimSingle,nBinAssElectron);
+  fAssElectron->SetBinEdges(0,binLimC);
+  fAssElectron->SetBinEdges(1,fPtBinning.GetArray());
+  fAssElectron->SetBinEdges(2,binLimSource);
+  fAssElectron->SetBinEdges(3,binLimEtaAssociat);
+  fAssElectron->Sumw2();
+  AliDebug(2,"AliHFENonPhotonicElectron: fAssElectron");
+
+  // Inclusive Electron
+  Int_t nBinIncElectron[nDimSingle] = {nBinsC,fPtBinning.GetSize()-1,nBinsSource,fEtaBinning.GetSize()-1};
+  fIncElectron = new THnSparseF("fIncElectron","fIncElectron",nDimSingle,nBinIncElectron);
+  fIncElectron->SetBinEdges(0,binLimC);
+  fIncElectron->SetBinEdges(1,fPtBinning.GetArray());
+  fIncElectron->SetBinEdges(2,binLimSource);
+  fIncElectron->SetBinEdges(3,fEtaBinning.GetArray());
+  fIncElectron->Sumw2();
+  AliDebug(2,"AliHFENonPhotonicElectron: fIncElectron");
 
   // ee invariant mass Unlike Sign
-  const Int_t nDimUSign=6;
-  Int_t nBinUSign[nDimUSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource,nBinsAngle};
-  fUSignee = new THnSparseF("fUSignee","fUSignee",nDimUSign,nBinUSign);
-  fUSignee->SetBinEdges(0,binLimPhi);
-  fUSignee->SetBinEdges(1,binLimC);
-  fUSignee->SetBinEdges(2,binLimPt);
-  fUSignee->SetBinEdges(3,binLimInvMass);
-  fUSignee->SetBinEdges(4,binLimSource);
-  fUSignee->SetBinEdges(5,binLimAngle);
-/*  fUSignee->GetAxis(0)->SetTitle("#Delta #phi");
-  fUSignee->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fUSignee->GetAxis(2)->SetTitle("#it{p}_{T} (GeV/c)");
-  fUSignee->GetAxis(3)->SetTitle("m_{e^{#pm}e^{#mp}}} (GeV/#it{c}^{2}");
-  fUSignee->GetAxis(4)->SetTitle("MC");
-  fUSignee->GetAxis(5)->SetTitle("#alpha");*/
-  fUSignee->Sumw2();
-  AliDebug(2,"AliHFENonPhotonicElectron: fUSignee");
+  fUSign = new THnSparseF("fUSign","fUSign",nDimPair,nBinPair);
+  fUSign->SetBinEdges(0,binLimPhi);
+  fUSign->SetBinEdges(1,binLimC);
+  fUSign->SetBinEdges(2,fPtBinning.GetArray());
+  fUSign->SetBinEdges(3,binLimInvMass);
+  fUSign->SetBinEdges(4,binLimSource);
+  fUSign->SetBinEdges(5,binLimAngle);
+  fUSign->SetBinEdges(6,fPtBinning.GetArray());
+  fUSign->SetBinEdges(7,fEtaBinning.GetArray());
+  fUSign->SetBinEdges(8,binLimEtaAssociat);
+  fUSign->Sumw2();
+  AliDebug(2,"AliHFENonPhotonicElectron: fUSign");
 
   // ee invariant mass Like Sign
-  const Int_t nDimLSign=6;
-  Int_t nBinLSign[nDimLSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource,nBinsAngle};
-  fLSignee = new THnSparseF("fLSignee","fLSignee",nDimLSign,nBinLSign);
-  fLSignee->SetBinEdges(0,binLimPhi);
-  fLSignee->SetBinEdges(1,binLimC);
-  fLSignee->SetBinEdges(2,binLimPt);
-  fLSignee->SetBinEdges(3,binLimInvMass);
-  fLSignee->SetBinEdges(4,binLimSource);
-  fLSignee->SetBinEdges(5,binLimAngle);
-/*  fLSignee->GetAxis(0)->SetTitle("#Delta #phi");
-  fLSignee->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fLSignee->GetAxis(2)->SetTitle("#it{p}_{T} (GeV/c)");
-  fLSignee->GetAxis(3)->SetTitle("m_{e^{#pm}e^{#pm}}} (GeV/#it{c}^{2}");
-  fLSignee->GetAxis(4)->SetTitle("MC");
-  fLSignee->GetAxis(5)->SetTitle("#alpha");*/
-  fLSignee->Sumw2();
-  AliDebug(2,"AliHFENonPhotonicElectron: fLSignee");
-
+  fLSign = new THnSparseF("fLSign","fLSign",nDimPair,nBinPair);
+  fLSign->SetBinEdges(0,binLimPhi);
+  fLSign->SetBinEdges(1,binLimC);
+  fLSign->SetBinEdges(2,fPtBinning.GetArray());
+  fLSign->SetBinEdges(3,binLimInvMass);
+  fLSign->SetBinEdges(4,binLimSource);
+  fLSign->SetBinEdges(5,binLimAngle);
+  fLSign->SetBinEdges(6,fPtBinning.GetArray());
+  fLSign->SetBinEdges(7,fEtaBinning.GetArray());
+  fLSign->SetBinEdges(8,binLimEtaAssociat);
+  fLSign->Sumw2();
+  AliDebug(2,"AliHFENonPhotonicElectron: fLSign");
+
+  // Histograms counting the number of like sign / unlike sign matches per inclusive track
+  const Int_t nBinsMatches = 50;
+  Double_t binLimMatches[nBinsMatches+1];
+  for(int ib = 0; ib <= nBinsMatches; ib++) binLimMatches[ib] = ib;
+  const Int_t nDimMatches = 3;  // centrality, pt_inc, number of matches 
+  const Int_t nBinsMatchHist[nDimMatches] = {nBinsC, fPtBinning.GetSize()-1, nBinsMatches};
+  fUSmatches = new THnSparseF("fUSmatches", "fUSmatches", nDimMatches, nBinsMatchHist);
+  fUSmatches->SetBinEdges(0,binLimC);
+  fUSmatches->SetBinEdges(1,fPtBinning.GetArray());
+  fUSmatches->SetBinEdges(2,binLimMatches);
+
+  fLSmatches = new THnSparseF("fLSmatches", "fLSmatches", nDimMatches, nBinsMatchHist);
+  fLSmatches->SetBinEdges(0,binLimC);
+  fLSmatches->SetBinEdges(1,fPtBinning.GetArray());
+  fLSmatches->SetBinEdges(2,binLimMatches);
+
+  // Histograms for radius studies
+  Int_t nBinsradius = 50;
+  Double_t minradius = 0.0;
+  Double_t maxradius = 100.0;
+  Double_t binLimradius[nBinsradius+1];
+  for(Int_t i=0; i<=nBinsradius; i++) binLimradius[i]=(Double_t)minradius + (maxradius-minradius)/nBinsradius*(Double_t)i ;
+  const Int_t nDimIncElectronRadius = 3;  // centrality, pt_inc, radius 
+  const Int_t nBinsIncElectronRadius[nDimIncElectronRadius] = {nBinsC, fPtBinning.GetSize()-1, nBinsradius};
+  fIncElectronRadius = new THnSparseF("fIncElectronRadius", "fIncElectronRadius", nDimIncElectronRadius, nBinsIncElectronRadius);
+  fIncElectronRadius->SetBinEdges(0,binLimC);
+  fIncElectronRadius->SetBinEdges(1,fPtBinning.GetArray());
+  fIncElectronRadius->SetBinEdges(2,binLimradius);
+
+  fRecElectronRadius = new THnSparseF("fRecElectronRadius", "fRecElectronRadius", nDimIncElectronRadius, nBinsIncElectronRadius);
+  fRecElectronRadius->SetBinEdges(0,binLimC);
+  fRecElectronRadius->SetBinEdges(1,fPtBinning.GetArray());
+  fRecElectronRadius->SetBinEdges(2,binLimradius);
+
+/*
   // ee angle Unlike Sign
   const Int_t nDimUSignAngle=3;
   Int_t nBinUSignAngle[nDimUSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
@@ -298,9 +406,6 @@ void AliHFENonPhotonicElectron::Init()
   fUSignAngle->SetBinEdges(0,binLimAngle);
   fUSignAngle->SetBinEdges(1,binLimC);
   fUSignAngle->SetBinEdges(2,binLimSource);
-/*  fUSignAngle->GetAxis(0)->SetTitle("#alpha");
-  fUSignAngle->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fUSignAngle->GetAxis(2)->SetTitle("MC");*/
   fUSignAngle->Sumw2();
   AliDebug(2,"AliHFENonPhotonicElectron: fUSignAngle");
 
@@ -311,23 +416,65 @@ void AliHFENonPhotonicElectron::Init()
   fLSignAngle->SetBinEdges(0,binLimAngle);
   fLSignAngle->SetBinEdges(1,binLimC);
   fLSignAngle->SetBinEdges(2,binLimSource);
-/*  fLSignAngle->GetAxis(0)->SetTitle("#alpha");
-  fLSignAngle->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fLSignAngle->GetAxis(2)->SetTitle("MC");*/
   fLSignAngle->Sumw2();
   AliDebug(2,"AliHFENonPhotonicElectron: fLSignAngle");
+*/
+
+  // control histogram for ITS PID
+  fHnsigmaITS = new TH2F("fHnsigmaITS", "Number of sigmas in the ITS", 30, 0., 0.3, 1200, -10., 10.);
+
+  // control histogram for weights sources
+  fWeightsSource = new TH2F("fWeightsSource", "Source code for weights", 11, -1.5, 9.5, 29, -1.5, 27.5);
+
+  fListOutput->Add(fAssElectron);
+  fListOutput->Add(fIncElectron);
+  fListOutput->Add(fUSign);
+  fListOutput->Add(fLSign);
+  fListOutput->Add(fUSmatches);
+  fListOutput->Add(fLSmatches);
+  fListOutput->Add(fHnsigmaITS);
+  fListOutput->Add(fWeightsSource);
+  fListOutput->Add(fIncElectronRadius);
+  fListOutput->Add(fRecElectronRadius);
+//  fListOutput->Add(fUSignAngle);
+//  fListOutput->Add(fLSignAngle);
+
+}
+
+//_____________________________________________________________________________________________
+void AliHFENonPhotonicElectron::SetWithWeights(Int_t levelBack)
+{
+  //
+  // Init the HFE level
+  //
+  if(levelBack >= 0) fLevelBack = levelBack;
 
+}
 
-  fListOutput->Add(fMCSourceee);
-  fListOutput->Add(fUSignee);
-  fListOutput->Add(fLSignee);
-  fListOutput->Add(fUSignAngle);
-  fListOutput->Add(fLSignAngle);
+//_____________________________________________________________________________________________
+void AliHFENonPhotonicElectron::SetMCEvent(AliMCEvent *mcEvent)
+{
+  //
+  // Pass the mcEvent
+  //
+  
+  fMCEvent = mcEvent;
+  
+}
 
+//_____________________________________________________________________________________________
+void AliHFENonPhotonicElectron::SetAODArrayMCInfo(TClonesArray *aodArrayMCInfo)
+{
+  //
+  // Pass the mcEvent info
+  //
+
+  fAODArrayMCInfo = aodArrayMCInfo;
+  
 }
 
 //_____________________________________________________________________________________________
-void AliHFENonPhotonicElectron::InitRun(const AliVEvent *inputEvent,AliPIDResponse *pidResponse)
+void AliHFENonPhotonicElectron::InitRun(const AliVEvent *inputEvent,const AliPIDResponse *pidResponse)
 {
   //
   // Init run
@@ -353,115 +500,94 @@ void AliHFENonPhotonicElectron::InitRun(const AliVEvent *inputEvent,AliPIDRespon
 }
 
 //_____________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent,Int_t binct)
+Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent, Int_t binct)
 {
   //
   // Fill the pool of associated tracks
   // Return the number of associated tracks
   //
 
+  fnumberfound = 0;
+
+  fHFEBackgroundCuts->SetRecEvent(inputEvent);
   Int_t nbtracks = inputEvent->GetNumberOfTracks();
 
-  if( fArraytrack )
-  {
+  if( fArraytrack ){
     fArraytrack->~TArrayI();
     new(fArraytrack) TArrayI(nbtracks);
-  }
-  else
-  {
+  } else {
     fArraytrack = new TArrayI(nbtracks);
   }
 
   fCounterPoolBackground = 0;
 
-  for(Int_t k = 0; k < nbtracks; k++)
-  {
+  Bool_t isSelected(kFALSE);
+  Bool_t isAOD = (dynamic_cast<AliAODEvent *>(inputEvent) != NULL);
+  AliDebug(2, Form("isAOD: %s", isAOD ? "yes" : "no"));
+  for(Int_t k = 0; k < nbtracks; k++) {
     AliVTrack *track = (AliVTrack *) inputEvent->GetTrack(k);
     if(!track) continue;
-
-    // Track cuts
-    Bool_t survivedbackground = kTRUE;
-    AliAODEvent *aodeventu = dynamic_cast<AliAODEvent *>(inputEvent);
-
-    if(aodeventu)
-    {
-      // AOD
-      AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
-      if(aodtrack)
-      {
-       // filter
-       if(fUseFilterAOD)
-       {
-         if(!(aodtrack->TestFilterBit(fFilter))) survivedbackground = kFALSE;
-       }
-
-       // additional cuts
-       if(survivedbackground)
-       {
-         AliESDtrack esdTrack(aodtrack);
-
-         // set the TPC cluster info
-         esdTrack.SetTPCClusterMap(aodtrack->GetTPCClusterMap());
-         esdTrack.SetTPCSharedMap(aodtrack->GetTPCSharedMap());
-         esdTrack.SetTPCPointsF(aodtrack->GetTPCNclsF());
-         AliAODVertex *vAOD = aodeventu->GetPrimaryVertex();
-         Double_t bfield = aodeventu->GetMagneticField();
-         Double_t pos[3],cov[6];
-         vAOD->GetXYZ(pos);
-         vAOD->GetCovarianceMatrix(cov);
-         const AliESDVertex vESD(pos,cov,100.,100);
-         esdTrack.RelateToVertex(&vESD,bfield,3.);
-
-         if(!fHFEBackgroundCuts->IsSelected(&esdTrack))
-         {
-           survivedbackground = kFALSE;
-         }
-       }
-      }
-    }
-    else
-    {
-      // ESD
-      AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-      if(esdtrack)
-      {
-       if(!fHFEBackgroundCuts->IsSelected(esdtrack)) survivedbackground = kFALSE;
-      }
+    
+    //
+    isSelected = kFALSE;
+    if(fSelectCategory1tracks && FilterCategory1Track(track, isAOD, binct)) isSelected = kTRUE;
+    else if(fSelectCategory2tracks && FilterCategory2Track(track, isAOD)) isSelected = kTRUE;
+
+    if(isSelected){
+           AliDebug(2,Form("fCounterPoolBackground %d, track %d",fCounterPoolBackground,k));
+      fArraytrack->AddAt(k,fCounterPoolBackground);
+      fCounterPoolBackground++;
     }
+  } // loop tracks
 
-    // PID
-    if(survivedbackground)
-    {
-      // PID track cuts
-      AliHFEpidObject hfetrack2;
-
-      if(!aodeventu) hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
-      else hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
-
-      hfetrack2.SetRecTrack(track);
-      if(binct>-1)
-      {
-       hfetrack2.SetCentrality((Int_t)binct);
-       AliDebug(2,Form("centrality %d and %d",binct,hfetrack2.GetCentrality()));
-       hfetrack2.SetPbPb();
-      }
+  //printf(Form("Associated Pool: Tracks %d, fCounterPoolBackground %d \n", nbtracks, fCounterPoolBackground));
+
+  return fCounterPoolBackground;
 
-      if(fPIDBackground->IsSelected(&hfetrack2,0x0,"recTrackCont",fPIDBackgroundQA))
-      {
-       fArraytrack->AddAt(k,fCounterPoolBackground);
-       fCounterPoolBackground++;
-       AliDebug(2,Form("fCounterPoolBackground %d, track %d",fCounterPoolBackground,k));
+}
+
+//_____________________________________________________________________________________________
+Int_t AliHFENonPhotonicElectron::CountPoolAssociated(AliVEvent *inputEvent, Int_t binct)
+{
+  //
+  // Count the pool of assiocated tracks
+  //
+
+
+  if(fnumberfound > 0) //!count only events with an inclusive electron
+  {
+    Double_t valueAssElectron[4] = { binct, -1, -1};           //Centrality    Pt      Source
+    Int_t iTrack2 = 0;
+    Int_t indexmother2 = -1;
+    AliVTrack *track2 = 0x0;
+
+    for(Int_t ii = 0; ii < fCounterPoolBackground; ii++){
+      iTrack2 = fArraytrack->At(ii);
+      AliDebug(2,Form("track %d",iTrack2));
+      track2 = (AliVTrack *)inputEvent->GetTrack(iTrack2);
+
+      if(!track2){
+             //printf("ERROR: Could not receive track %d", iTrack2);
+             continue;
       }
-    }
-  } // loop tracks
 
-  //printf(Form("Test pool: nbtrack %d, binct %d \n fCounterPoolBackground %d \n",nbtracks,binct,fCounterPoolBackground));
+      // if MC look
+      if(fMCEvent || fAODArrayMCInfo) valueAssElectron[2] = FindMother(TMath::Abs(track2->GetLabel()), indexmother2) ;
 
-  return fCounterPoolBackground;
+      fkPIDRespons = fPIDBackground->GetPIDResponse();
+
+      valueAssElectron[1] = track2->Pt() ;
+      valueAssElectron[3] = track2->Eta() ;
+
+      fAssElectron->Fill( valueAssElectron) ;
+    }
+  //printf(Form("Associated Pool: fCounterPoolBackground %d \n", fCounterPoolBackground));
+  }
+  return fnumberfound;
 }
 
 //_____________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *vEvent, Double_t weight, Int_t binct, Double_t deltaphi, Int_t source, Int_t indexmother)
+Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *vEvent, Double_t weight, Int_t binct, Double_t deltaphi, Int_t source, Int_t indexmother,Int_t mcQAsource)
 {
   //
   // Look At Non HFE
@@ -486,237 +612,172 @@ Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1,
    *                                                                                   *
    ***********************************************************************************/
 
+  //printf("weight %f and source %d\n",weight,source);
+
+  
   AliAODEvent *aodeventu = dynamic_cast<AliAODEvent*>(vEvent);
   Int_t taggedphotonic = -1;
 
   AliDebug(2,Form("fCounterPoolBackground %d in LookAtNonHFE!!!",fCounterPoolBackground));
   if(!fArraytrack) return taggedphotonic;
   AliDebug(2,Form("process track %d",iTrack1));
+  AliDebug(1,Form("Inclusive source is %d\n", source));
 
-  //Set Fill-Arrays for THnSparse
-  Double_t value_MCSource[3] = { binct, track1->Pt(), source};         //Centrality    Pt              Source
-  Double_t valueee[6] = { deltaphi, binct, track1->Pt(), -1, source,-1};       //DeltaPhi      Centrality      Pt      InvariantMass   Source
-  Double_t valueAngle[3] = {   -1,     binct,          source};                //Angle         Centrality      Source
+  fkPIDRespons = fPIDBackground->GetPIDResponse();
 
-  Bool_t USignPhotonic = kFALSE;
-  Bool_t LSignPhotonic = kFALSE;
-  Bool_t hasdcaT1 = kFALSE;
-  Bool_t hasdcaT2 = kFALSE;
+  //Set Fill-Arrays for THnSparse
+  Double_t valueIncElectron[4] = { binct, track1->Pt(), source, track1->Eta()};        //Centrality    Pt      Source  P       
+  Double_t valueSign[9]                = { deltaphi, binct, track1->Pt(), -1, source, -1, -1, track1->Eta(), -1};                      //DeltaPhi      Centrality      Pt      InvariantMass   Source  Angle   Pt
+  //Double_t valueAngle[3]     = { -1, binct, source}; 
+  Double_t valueradius[3]      = { binct, track1->Pt(), 0.};                                                   //Angle         Centrality      Source
 
   Int_t pdg1 = CheckPdg(TMath::Abs(track1->GetLabel()));
-  Int_t pdg2 = -100;
-  Int_t numberfound = 0;
+  Double_t radius = Radius(TMath::Abs(track1->GetLabel()));
+  AliKFParticle::SetField(vEvent->GetMagneticField());
+  AliKFVertex primV(*(vEvent->GetPrimaryVertex()));
+  valueradius[2] = radius;
+
+  AliVTrack *track2(NULL);
   Int_t iTrack2 = 0;
-  Int_t source2 = 0;
   Int_t indexmother2 = -1;
-  Int_t fPDGtrack1 = 0;
-  Int_t fPDGtrack2 = 0;
-
-  Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
-  Double_t bfield = vEvent->GetMagneticField();
-  Double_t xt1 = 0; //radial position track 1 at the DCA point
-  Double_t xt2 = 0; //radial position track 2 at the DCA point
-  Double_t p1[3] = {0,0,0};
-  Double_t p2[3] = {0,0,0};
-  Double_t dca12 = 0;
-
-  Double_t angleESD = 0;
-  Double_t invmassESD = 0;
-
-  Double_t chi2OverNDF = 0;
-  Double_t width = 0;
-  Double_t angleAOD = 0;
-  Double_t invmassAOD = 0;
-
-  Float_t fCharge1 = 0;
+  Int_t pdg2 = -100;
+  Int_t source2 = -1;
   Float_t fCharge2 = 0;
 
-  TLorentzVector electron1;
-  TLorentzVector electron2;
-  TLorentzVector mother;
+  // count number of matches with opposite/same sign track in the given mass range
+  Int_t countsMatchLikesign(0),
+        countsMatchUnlikesign(0);
 
-  AliVTrack *track2;
-  AliESDtrack *esdtrack1;
-  AliESDtrack *esdtrack2;
-  AliKFParticle *ktrack1;
-  AliKFParticle *ktrack2;
-  AliKFVertex primV(*(vEvent->GetPrimaryVertex()));
+  // Values to fill
+  Double_t angle(-1.);
+  Double_t invmass(-1);
 
-  // FILL MCsource TODO: correct?
-  fMCSourceee->Fill(&value_MCSource[0],weight);
+  Float_t fCharge1 = track1->Charge();                                                 //Charge from track1
 
+  Bool_t kUSignPhotonic = kFALSE;
+  Bool_t kLSignPhotonic = kFALSE;
 
-  for(Int_t idex = 0; idex < fCounterPoolBackground; idex++)
-  {
+  //! FILL Inclusive Electron
+  fWeightsSource->Fill(source,mcQAsource);
+  fIncElectron->Fill(valueIncElectron,weight);
+  fnumberfound++;
+  if(source == kElectronfromconversion) {
+    fIncElectronRadius->Fill(valueradius,weight);
+    //printf("radius %f\n",radius);
+  }
+  //printf(Form("Inclusive Pool: TrackNr. %d, fnumberfound %d \n", iTrack1, fnumberfound));
+
+  for(Int_t idex = 0; idex < fCounterPoolBackground; idex++){
     iTrack2 = fArraytrack->At(idex);
     AliDebug(2,Form("track %d",iTrack2));
     track2 = (AliVTrack *)vEvent->GetTrack(iTrack2);
 
-    if(!track2)
-    {
-      printf("ERROR: Could not receive track %d", iTrack2);
+    if(!track2){
+      //printf("ERROR: Could not receive track %d", iTrack2);
       continue;
     }
 
-    if(iTrack2==iTrack1) continue;
-    AliDebug(2,"Different");
-
-    fCharge1 = track1->Charge();               //Charge from track1
     fCharge2 = track2->Charge();               //Charge from track2
 
     // Reset the MC info
-    valueAngle[2] = source;
-    valueee[4] = source;
+    //valueAngle[2] = source;
+    valueSign[4] = source;
+    valueSign[6] = track2->Pt();
+    valueSign[8] = track2->Eta();
 
     // track cuts and PID already done
 
+    // Checking if it is the same Track!
+    if(iTrack2==iTrack1) continue;
+    AliDebug(2,"Different");
+
     // if MC look
-    if(fMCEvent || fAODArrayMCInfo)
-    {
+    if(fMCEvent || fAODArrayMCInfo){
+      AliDebug(2, "Checking for source");
       source2   = FindMother(TMath::Abs(track2->GetLabel()), indexmother2);
+      AliDebug(2, Form("source is %d", source2));
       pdg2      = CheckPdg(TMath::Abs(track2->GetLabel()));
 
-      if(source2 >=0 )
-      {
-       if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0))
-       {
-         if(source == kElectronfromconversion)
-         {
-           valueAngle[2] = kElectronfromconversionboth;
-           valueee[4] = kElectronfromconversionboth;
-           numberfound++;
-         }
-
-         if(source == kElectronfrompi0)
-         {
-           valueAngle[2] = kElectronfrompi0both;
-           valueee[4] = kElectronfrompi0both;
-         }
-
-         if(source == kElectronfrometa)
-         {
-           valueAngle[2] = kElectronfrometaboth;
-           valueee[4] = kElectronfrometaboth;
-         }
-       }
+      if(source == kElectronfromconversion){
+        AliDebug(2, Form("Electron from conversion (source %d), paired with source %d", source, source2));
+        AliDebug(2, Form("Index of the mothers: incl %d, associated %d", indexmother, indexmother2));
+        AliDebug(2, Form("PDGs: incl %d, associated %d", pdg1, pdg2));
       }
-    }
-
-    if(fAlgorithmMA && (!aodeventu)) 
-    {
-      /**                              *
-       *       ESD-Analysis            *
-       **                              */
-
-      esdtrack1 = dynamic_cast<AliESDtrack *>(track1);                 //ESD-track1
-      esdtrack2 = dynamic_cast<AliESDtrack *>(track2);                 //ESD-track2
-      if((!esdtrack1) || (!esdtrack2)) continue;
-
-      dca12 = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);             //DCA track1-track2
 
-      if(dca12 > fMaxDCA) continue;                                    //! Cut on DCA
-
-      //Momento of the track extrapolated to DCA track-track
-      hasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);                //Track1
-      hasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);                //Track2
-      if(!hasdcaT1 || !hasdcaT2) AliWarning("It could be a problem in the extrapolation");
-
-      electron1.SetXYZM(esdtrack1->Px(), esdtrack1->Py(), esdtrack1->Pz(), eMass);
-      electron2.SetXYZM(esdtrack2->Px(), esdtrack2->Py(), esdtrack2->Pz(), eMass);
-
-      mother      = electron1 + electron2;
-      invmassESD  = mother.M();
-      angleESD    = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
-
-      valueAngle[0] = angleESD;
-      valueee[5] = angleESD;
-
-      if((fCharge1*fCharge2)>0.0)      fLSignAngle->Fill(&valueAngle[0],weight);
-      else                             fUSignAngle->Fill(&valueAngle[0],weight);
-
-      if(angleESD > fMaxOpening3D) continue;                            //! Cut on Opening Angle
-      valueee[3] = invmassESD;
-
-      if((fCharge1*fCharge2)>0.0)      fLSignee->Fill(&valueee[0],weight);
-      else                             fUSignee->Fill(&valueee[0],weight);
-
-      if(invmassESD > fMaxInvMass) continue;                           //! Cut on Invariant Mass
-
-      if((fCharge1*fCharge2)>0.0)      LSignPhotonic=kTRUE;
-      else                             USignPhotonic=kTRUE;
-    }
-    else
-    {
-      /**                              *
-       *       AliKF-Analysis          *
-       **                              */
-
-      fPDGtrack1 = 11;
-      fPDGtrack2 = 11;
-
-      if(fCharge1>0) fPDGtrack1 = -11;
-      if(fCharge2>0) fPDGtrack2 = -11;
-
-      ktrack1 = new AliKFParticle(*track1, fPDGtrack1);
-      ktrack2 = new AliKFParticle(*track2, fPDGtrack2);
-      AliKFParticle recoGamma(*ktrack1,*ktrack2);
-
-      if(recoGamma.GetNDF()<1) continue;                               //! Cut on Reconstruction
-
-      chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
-      if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) continue;
-
-      // DCA
-      //Double_t dca12 = ktrack1.GetDistanceFromParticle(ktrack2);
-      //if(dca12 > fMaxDCA) continue;
-
-      // if set mass constraint
-      if(fSetMassConstraint) //&& pVtx)
-      {
-       primV += recoGamma;
-       primV -= *ktrack1;
-       primV -= *ktrack2;
-       recoGamma.SetProductionVertex(primV);
-       recoGamma.SetMassConstraint(0,0.0001);
+      if(source2 >=0 ){
+             if((indexmother2 == indexmother) && (source == source2) && ((pdg1*pdg2)<0.0)){
+          AliDebug(1, "Real pair");
+          switch(source){
+            case kElectronfromconversion: 
+                 valueSign[4] = kElectronfromconversionboth; 
+                 break;
+            case kElectronfrompi0: 
+                 valueSign[4] = kElectronfrompi0both; 
+                 break;
+            case kElectronfrometa:
+                 valueSign[4] = kElectronfrometaboth;
+                 break;
+          };
+        }
       }
+    }
 
-      recoGamma.GetMass(invmassAOD,width);
-      angleAOD = ktrack1->GetAngle(*ktrack2);
-
-      valueAngle[0] = angleAOD;
-      valueee[5] = angleAOD;
-
-      if((fCharge1*fCharge2)>0.0)      fLSignAngle->Fill(&valueAngle[0],weight);
-      else                             fUSignAngle->Fill(&valueAngle[0],weight);
-
-      if(angleAOD > fMaxOpening3D) continue;                           //! Cut on Opening Angle
-
-      valueee[3] = invmassAOD;
-
-      if((fCharge1*fCharge2)>0.0)      fLSignee->Fill(&valueee[0],weight);
-      else                             fUSignee->Fill(&valueee[0],weight);
-
-      if(invmassAOD > fMaxInvMass) continue;                           //! Cut on Invariant Mass
+    if(fAlgorithmMA){
+      // Use TLorentzVector
+      if(!MakePairDCA(track1, track2, vEvent, (aodeventu != NULL), invmass, angle)) continue;
+    } else {
+      // Use AliKF package
+      if(!MakePairKF(track1, track2, primV, invmass, angle)) continue;
+    }
 
-      if((fCharge1*fCharge2)>0.0)      LSignPhotonic=kTRUE;
-      else                             USignPhotonic=kTRUE;
+    valueSign[3] = invmass;
+    valueSign[5] = angle;
+
+    //if((fCharge1*fCharge2)>0.0)      fLSignAngle->Fill(&valueAngle[0],weight);
+    //else                             fUSignAngle->Fill(&valueAngle[0],weight);
+
+    if(angle > fMaxOpening3D) continue;                                 //! Cut on Opening Angle
+    if(invmass > fMaxInvMass) continue;                                //! Cut on Invariant Mass
+
+    if((fCharge1*fCharge2)>0.0){       
+      if(invmass < 1.0) fLSign->Fill( valueSign, weight);
+      // count like-sign background matched pairs per inclusive based on mass cut
+      if(invmass < 0.14) countsMatchLikesign++;
+      AliDebug(1, "Selected Like sign");
+    } else {
+      if(invmass < 1.0)fUSign->Fill( valueSign, weight);
+      // count unlike-sign matched pairs per inclusive based on mass cut
+      if(invmass < 0.14) {
+       countsMatchUnlikesign++;
+       if(source == kElectronfromconversionboth) fRecElectronRadius->Fill(valueradius,weight);
+      }
+      AliDebug(1, "Selected Unike sign");
     }
+
+    if((fCharge1*fCharge2)>0.0)        kLSignPhotonic=kTRUE;
+    else                               kUSignPhotonic=kTRUE;
   }
 
-  if( USignPhotonic &&  LSignPhotonic) taggedphotonic = 6;
-  if(!USignPhotonic &&  LSignPhotonic) taggedphotonic = 4;
-  if( USignPhotonic && !LSignPhotonic) taggedphotonic = 2;
+  // Fill counted
+  Double_t valCountsLS[3] = {binct, track1->Pt(), countsMatchLikesign},
+           valCountsUS[3] = {binct, track1->Pt(), countsMatchUnlikesign}; 
+  fUSmatches->Fill(valCountsUS);
+  fLSmatches->Fill(valCountsLS);
+
+  if( kUSignPhotonic &&  kLSignPhotonic) taggedphotonic = 6;
+  if(!kUSignPhotonic &&  kLSignPhotonic) taggedphotonic = 4;
+  if( kUSignPhotonic && !kLSignPhotonic) taggedphotonic = 2;
 
   return taggedphotonic;
 }
 
 //_________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::FindMother(Int_t tr, Int_t &indexmother){
+Int_t AliHFENonPhotonicElectron::FindMother(Int_t tr, Int_t &indexmother) const {
   //
   // Find the mother if MC
   //
 
-  if(!fMCEvent && !fAODArrayMCInfo) return 0;
+  if(!fMCEvent && !fAODArrayMCInfo) return -1;
 
   Int_t pdg = CheckPdg(tr);
   if(TMath::Abs(pdg)!= 11)
@@ -740,7 +801,7 @@ Int_t AliHFENonPhotonicElectron::FindMother(Int_t tr, Int_t &indexmother){
 }
 
 //________________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::CheckPdg(Int_t tr) {
+Int_t AliHFENonPhotonicElectron::CheckPdg(Int_t tr) const {
 
   //
   // Return the pdg of the particle
@@ -749,336 +810,402 @@ Int_t AliHFENonPhotonicElectron::CheckPdg(Int_t tr) {
   Int_t pdgcode = -1;
   if(tr < 0) return pdgcode;
 
-  if(fMCEvent)
-  {
+  AliMCParticle *mctrackesd = NULL; AliAODMCParticle *mctrackaod = NULL;
+  if(fMCEvent){
     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
-    if(!mctrack) return -1;
-    AliMCParticle *mctrackesd = NULL;
-    if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return pdgcode;
-    pdgcode = mctrackesd->PdgCode();
-  }
-
-  if(fAODArrayMCInfo)
-  {
-    if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
-    AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
-    if(!mctrackaod) return pdgcode;
-    pdgcode = mctrackaod->GetPdgCode();
+    if(mctrack){
+      if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) pdgcode = mctrackesd->PdgCode();
+      else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) pdgcode = mctrackaod->GetPdgCode(); 
+    }
+  } else if(fAODArrayMCInfo) {
+    if(tr < fAODArrayMCInfo->GetEntriesFast()){
+      mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
+      if(mctrackaod) return pdgcode = mctrackaod->GetPdgCode();
+    }
   }
 
   return pdgcode;
 }
 
-//_______________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::IsMotherGamma(Int_t tr) {
+//________________________________________________________________________________________________
+Double_t AliHFENonPhotonicElectron::Radius(Int_t tr) const {
 
   //
-  // Return the lab of gamma mother or -1 if not gamma
+  // Return the production vertex radius
   //
 
-  if(tr < 0) return -1;
+  Double_t radius = 0.;
+  if(tr < 0) return radius;
 
-  if(fMCEvent)
-  {
+  AliMCParticle *mctrackesd = NULL; AliAODMCParticle *mctrackaod = NULL;
+  if(fMCEvent){
     AliVParticle *mctrack = fMCEvent->GetTrack(tr);
-    if(!mctrack) return -1;
-    AliMCParticle *mctrackesd = NULL;
-    if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
-    TParticle *particle = 0x0;
-    particle = mctrackesd->Particle();
-
-    // Take mother
-    if(!particle) return -1;
-    Int_t imother   = particle->GetFirstMother();
-    if(imother < 0) return -1;
-    AliMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
-    TParticle * mother = mothertrack->Particle();
-    if(!mother) return -1;
-
-    // Check gamma
-    Int_t pdg = mother->GetPdgCode();
-    if(TMath::Abs(pdg) == 22) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
-      return IsMotherGamma(imother);
+    if(mctrack){
+      if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) radius = TMath::Sqrt(mctrackesd->Xv()*mctrackesd->Xv()+mctrackesd->Yv()*mctrackesd->Yv());
+      else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) radius = TMath::Sqrt(mctrackaod->Xv()*mctrackaod->Xv()+mctrackaod->Yv()*mctrackaod->Yv());
     }
+  } else if(fAODArrayMCInfo) {
+    if(tr < fAODArrayMCInfo->GetEntriesFast()){
+      mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
+      if(mctrackaod) return radius = TMath::Sqrt(mctrackaod->Xv()*mctrackaod->Xv()+mctrackaod->Yv()*mctrackaod->Yv());
+    }
+  }
 
-    return -1;
+  return radius;
+}
+
+//_______________________________________________________________________________________________
+Int_t AliHFENonPhotonicElectron::GetMotherPDG(Int_t tr, Int_t &motherIndex) const {
+  //
+  // Returns the mother PDG of the track (return value) and the index of the
+  // mother track 
+  //
+  if(tr < 0) return -1;
+
+  Int_t pdg(-1);
+  AliMCParticle *mctrackesd(NULL); AliAODMCParticle *mctrackaod(NULL);
+
+  motherIndex = -1;
+  if(fMCEvent) {
+    AliDebug(2, "Using MC Event");
+    AliVParticle *mctrack = fMCEvent->GetTrack(tr);
+    if(mctrack){
+      if((mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))){
+        // Case ESD
+        TParticle *particle = mctrackesd->Particle();
+
+        // Take mother
+        if(particle){
+          motherIndex   = particle->GetFirstMother();
+          if(motherIndex >= 0){
+            AliMCParticle *mothertrack = NULL;
+            if((mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(motherIndex))))){
+              TParticle * mother = mothertrack->Particle();
+              pdg = mother->GetPdgCode();
+            }
+          }
+        }
+      } else if((mctrackaod = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))){
+        // Case AOD
+        // Take mother
+        motherIndex = mctrackaod->GetMother();
+        if(motherIndex >= 0){
+          AliAODMCParticle *mothertrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(motherIndex)); 
+          if(mothertrack) pdg = mothertrack->GetPdgCode();
+        }
+      }
+    }
+  } else if(fAODArrayMCInfo) {
+    AliDebug(2, "Using AOD list");
+    if(tr < fAODArrayMCInfo->GetEntriesFast()){ 
+      mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
+
+      // Take mother
+      if(mctrackaod){
+        motherIndex = mctrackaod->GetMother();
+        if(motherIndex >= 0 && motherIndex < fAODArrayMCInfo->GetEntriesFast()){
+          AliAODMCParticle *mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(motherIndex)));
+          if(mothertrack) pdg = mothertrack->GetPdgCode();
+        }
+      }
+    }
   }
+  return pdg;
+}
 
-  if(fAODArrayMCInfo)
-  {
-    if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
-    AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
-    if(!mctrackaod) return -1;
-
-    // Take mother
-    Int_t imother = mctrackaod->GetMother();
-    if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
-    AliAODMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
-
-    // Check gamma
-    Int_t pdg = mothertrack->GetPdgCode();
-    if(TMath::Abs(pdg) == 22) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
+//_______________________________________________________________________________________________
+Int_t AliHFENonPhotonicElectron::IsMotherGamma(Int_t tr) const {
+
+  //
+  // Return the lab of gamma mother or -1 if not gamma
+  //
+
+  Int_t imother(-1), pdg(-1);
+  pdg = GetMotherPDG(tr, imother);
+  
+  // Check gamma
+  if(imother >= 0){
+    if(TMath::Abs(pdg) == 22){
+      AliDebug(2, "Gamma Mother selected");
+      return imother;
+    }
+    if(TMath::Abs(pdg) == 11){
+      AliDebug(2, "Mother is electron - look further in hierarchy");
       return IsMotherGamma(imother);
     }
-
+    AliDebug(2, "Nothing selected");
     return -1;
   }
-
+  AliDebug(2, "Not mother");
   return -1;
 }
 
 //________________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::IsMotherPi0(Int_t tr) {
+Int_t AliHFENonPhotonicElectron::IsMotherPi0(Int_t tr) const {
 
   //
   // Return the lab of pi0 mother or -1 if not pi0
   //
 
-  if(tr < 0) return -1;
+  Int_t imother(-1), pdg(-1);
+  pdg = GetMotherPDG(tr, imother);
 
-  if(fMCEvent)
-  {
-    AliVParticle *mctrack = fMCEvent->GetTrack(tr);
-    if(!mctrack) return -1;
-    AliMCParticle *mctrackesd = NULL;
-    if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
-    TParticle *particle = 0x0;
-    particle = mctrackesd->Particle();
-    // Take mother
-    if(!particle) return -1;
-    Int_t imother   = particle->GetFirstMother();
-    if(imother < 0) return -1;
-    AliMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
-    TParticle * mother = mothertrack->Particle();
-    if(!mother) return -1;
-    // Check pi0
-    Int_t pdg = mother->GetPdgCode();
-    if(TMath::Abs(pdg) == 111) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
-      return IsMotherPi0(imother);
+  // Check pi0
+  if(imother >= 0){
+    if(TMath::Abs(pdg) == 111){
+      AliDebug(2, "Pi0 Mother selected");
+      return imother;
     }
-
-    return -1;
-  }
-
-  if(fAODArrayMCInfo)  {
-
-    if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
-    AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
-    if(!mctrackaod) return -1;
-
-    // Take mother
-    Int_t imother = mctrackaod->GetMother();
-    if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
-    AliAODMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
-    // Check pi0
-    Int_t pdg = mothertrack->GetPdgCode();
-    if(TMath::Abs(pdg) == 111) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
+    if(TMath::Abs(pdg) == 11){
+      AliDebug(2, "Mother is electron - look further in hierarchy");
       return IsMotherPi0(imother);
     }
-
+    AliDebug(2, "Nothing selected");
     return -1;
   }
-
+  AliDebug(2, "Not mother");
   return -1;
 }
 //________________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::IsMotherC(Int_t tr) {
+Int_t AliHFENonPhotonicElectron::IsMotherC(Int_t tr) const {
 
   //
   // Return the lab of signal mother or -1 if not from C
   //
 
-  if(tr < 0) return -1;
-
-  if(fMCEvent)
-  {
-    AliVParticle *mctrack = fMCEvent->GetTrack(tr);
-    if(!mctrack) return -1;
-    AliMCParticle *mctrackesd = NULL;
-    if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
-    TParticle *particle = 0x0;
-    particle = mctrackesd->Particle();
-
-    // Take mother
-    if(!particle) return -1;
-    Int_t imother   = particle->GetFirstMother();
-    if(imother < 0) return -1;
-    AliMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
-    TParticle * mother = mothertrack->Particle();
-    if(!mother) return -1;
+  Int_t imother(-1), pdg(-1);
+  pdg = GetMotherPDG(tr, imother);
 
     // Check C
-    Int_t pdg = mother->GetPdgCode();
-    if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
+  if(imother >= 0){
+    if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)){ 
+      AliDebug(2, "Charm Mother selected");
+      return imother;
+    }
+    if(TMath::Abs(pdg) == 11){
+      AliDebug(2, "Mother is electron - look further in hierarchy");
       return IsMotherC(imother);
     }
-
+    AliDebug(2, "Nothing selected");
     return -1;
   }
+  AliDebug(2, "Not mother");
+  return -1;
+}
 
-  if(fAODArrayMCInfo)
-  {
-    if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
-    AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
-    if(!mctrackaod) return -1;
+//_______________________________________________________________________________________________
+Int_t AliHFENonPhotonicElectron::IsMotherB(Int_t tr) const {
+
+  //
+  // Return the lab of signal mother or -1 if not B
+  //
 
-    // Take mother
-    Int_t imother = mctrackaod->GetMother();
-    if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
-    AliAODMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
+  Int_t imother(-1), pdg(-1);
+  pdg = GetMotherPDG(tr, imother);
 
-    // Check C
-    Int_t pdg = mothertrack->GetPdgCode();
-    if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
-      return IsMotherC(imother);
+  // Check B
+  if(imother >= 0){  
+    if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)){
+      AliDebug(2, "Bottom Mother selected");
+      return imother;
     }
-
+    if(TMath::Abs(pdg) == 11){
+      return IsMotherB(imother);
+      AliDebug(2, "Mother is electron - look further in hierarchy");
+    }
+    AliDebug(2, "Nothing selected");
     return -1;
   }
-
+  AliDebug(2, "Not mother");
   return -1;
 }
+
 //_______________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::IsMotherB(Int_t tr) {
+Int_t AliHFENonPhotonicElectron::IsMotherEta(Int_t tr) const {
 
   //
-  // Return the lab of signal mother or -1 if not B
+  // Return the lab of eta mother or -1 if not eta
   //
 
-  if(tr < 0) return -1;
+  Int_t imother(-1), pdg(-1);
+  pdg = GetMotherPDG(tr, imother);
 
-  if(fMCEvent)
-  {
-    AliVParticle *mctrack = fMCEvent->GetTrack(tr);
-    if(!mctrack) return -1;
-    AliMCParticle *mctrackesd = NULL;
-    if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
-    TParticle *particle = 0x0;
-    particle = mctrackesd->Particle();
-
-    // Take mother
-    if(!particle) return -1;
-    Int_t imother   = particle->GetFirstMother();
-    if(imother < 0) return -1;
-    AliMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
-    TParticle * mother = mothertrack->Particle();
-    if(!mother) return -1;
-
-    // Check B
-    Int_t pdg = mother->GetPdgCode();
-    if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
-      return IsMotherB(imother);
+  // Check eta
+  if(imother >= 0){  
+    if(TMath::Abs(pdg) == 221){
+      AliDebug(2, "Eta mother selected");
+      return imother;
     }
-
+    if(TMath::Abs(pdg) == 11){
+      AliDebug(2, "Mother is electron - look further in hierarchy");
+      return IsMotherEta(imother);
+    }
+    AliDebug(2, "Nothing selected");
     return -1;
   }
+  AliDebug(2, "Not mother");
+  return -1;
+}
 
-  if(fAODArrayMCInfo)
-  {
-    if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
-    AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
-    if(!mctrackaod) return -1;
-
-    // Take mother
-    Int_t imother = mctrackaod->GetMother();
-    if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
-    AliAODMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
-    // Check B
-    Int_t pdg = mothertrack->GetPdgCode();
-    if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
-      return IsMotherB(imother);
-    }
+//_______________________________________________________________________________________________
+Bool_t AliHFENonPhotonicElectron::MakePairDCA(const AliVTrack *inclusive, const AliVTrack *associated, AliVEvent *vEvent, Bool_t isAOD, Double_t &invMass, Double_t &angle) const {
+  //
+  // Make Pairs of electrons using TLorentzVector
+  //
+  Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
+  Double_t bfield = vEvent->GetMagneticField();
 
-    return -1;
+  AliESDtrack *esdtrack1, *esdtrack2;
+  if(isAOD){
+    // call copy constructor for AODs
+    esdtrack1 = new AliESDtrack(inclusive);
+    esdtrack2 = new AliESDtrack(associated);
+  } else {
+    // call copy constructor for ESDs
+    esdtrack1 = new AliESDtrack(*(static_cast<const AliESDtrack *>(inclusive)));
+    esdtrack2 = new AliESDtrack(*(static_cast<const AliESDtrack *>(associated)));
+  }
+  if((!esdtrack1) || (!esdtrack2)){ 
+    delete esdtrack1;
+    delete esdtrack2;
+    return kFALSE;
   }
 
-  return -1;
+  Double_t xt1 = 0; //radial position track 1 at the DCA point
+  Double_t xt2 = 0; //radial position track 2 at the DCA point
+  Double_t dca = esdtrack2->GetDCA(esdtrack1,bfield,xt2,xt1);          //DCA track1-track2
+  if(dca > fMaxDCA){
+    // Apply DCA cut already in the function
+    delete esdtrack1;
+    delete esdtrack2;
+    return kFALSE;
+  }
+
+  //Momenta of the track extrapolated to DCA track-track
+  Double_t p1[3] = {0,0,0};
+  Double_t p2[3] = {0,0,0};
+  Bool_t kHasdcaT1 = esdtrack1->GetPxPyPzAt(xt1,bfield,p1);            //Track1
+  Bool_t kHasdcaT2 = esdtrack2->GetPxPyPzAt(xt2,bfield,p2);            //Track2
+  if(!kHasdcaT1 || !kHasdcaT2) AliWarning("It could be a problem in the extrapolation");
+
+  TLorentzVector electron1, electron2, mother;
+  electron1.SetXYZM(p1[0], p1[1], p1[2], eMass);
+  electron2.SetXYZM(p2[0], p2[1], p2[2], eMass);
+
+  mother      = electron1 + electron2;
+  invMass  = mother.M();
+  angle    = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
+
+  delete esdtrack1;
+  delete esdtrack2;
+  return kTRUE;
 }
 
 //_______________________________________________________________________________________________
-Int_t AliHFENonPhotonicElectron::IsMotherEta(Int_t tr) {
-
+Bool_t AliHFENonPhotonicElectron::MakePairKF(const AliVTrack *inclusive, const AliVTrack *associated, AliKFVertex &primV, Double_t &invMass, Double_t &angle) const {
   //
-  // Return the lab of eta mother or -1 if not eta
+  // Make pairs of electrons using the AliKF package
   //
 
-  if(tr < 0) return -1;
+  //printf("AOD HFE non photonic\n");
 
-  if(fMCEvent)
-  {
-    AliVParticle *mctrack = fMCEvent->GetTrack(tr);
-    if(!mctrack) return -1;
-    AliMCParticle *mctrackesd = NULL;
-    if(!(mctrackesd = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(tr))))) return -1;
-    TParticle *particle = 0x0;
-    particle = mctrackesd->Particle();
-
-    // Take mother
-    if(!particle) return -1;
-    Int_t imother   = particle->GetFirstMother();
-    if(imother < 0) return -1;
-    AliMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(imother))))) return -1;
-    TParticle * mother = mothertrack->Particle();
-    if(!mother) return -1;
-
-    // Check eta
-    Int_t pdg = mother->GetPdgCode();
-    if(TMath::Abs(pdg) == 221) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
-      return IsMotherEta(imother);
-    }
+  Int_t fPDGtrack1 = 11;
+  if(inclusive->Charge()>0) fPDGtrack1 = -11;
+  Int_t fPDGtrack2 = 11;
+  if(associated->Charge()>0) fPDGtrack2 = -11;
 
-    return -1;
+  AliKFParticle ktrack1(*inclusive, fPDGtrack1);
+  AliKFParticle ktrack2(*associated, fPDGtrack2);
+  AliKFParticle recoGamma(ktrack1,ktrack2);
+
+  if(recoGamma.GetNDF()<1) return kFALSE;                              //! Cut on Reconstruction
+
+  Double_t chi2OverNDF = recoGamma.GetChi2()/recoGamma.GetNDF();
+  if(TMath::Sqrt(TMath::Abs(chi2OverNDF))>fChi2OverNDFCut) return kFALSE;
+
+  if(fSetMassConstraint){
+         primV += recoGamma;
+         primV -= ktrack1;
+         primV -= ktrack2;
+         recoGamma.SetProductionVertex(primV);
+         recoGamma.SetMassConstraint(0,0.0001);
   }
 
-  if(fAODArrayMCInfo)
-  {
-    if((tr+1)>fAODArrayMCInfo->GetEntriesFast()) return -1;
-    AliAODMCParticle *mctrackaod = (AliAODMCParticle *) fAODArrayMCInfo->At(tr);
-    if(!mctrackaod) return -1;
-
-    // Take mother
-    Int_t imother = mctrackaod->GetMother();
-    if(imother < 0 || ((imother+1)>fAODArrayMCInfo->GetEntriesFast())) return -1;
-    AliAODMCParticle *mothertrack = NULL;
-    if(!(mothertrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(TMath::Abs(imother))))) return -1;
-
-    // Check eta
-    Int_t pdg = mothertrack->GetPdgCode();
-    if(TMath::Abs(pdg) == 221) return imother;
-    if(TMath::Abs(pdg) == 11)
-    {
-      return IsMotherEta(imother);
+  Double_t width(0.);
+  recoGamma.GetMass(invMass,width);
+  angle = ktrack1.GetAngle(ktrack2);
+  return kTRUE;
+}
+
+//_______________________________________________________________________________________________
+Bool_t AliHFENonPhotonicElectron::FilterCategory1Track(const AliVTrack * const track, Bool_t isAOD, Int_t binct){
+  //
+  // Selection of good associated tracks for the pool
+  // selection is done using strong cuts
+  // Tracking in the TPC and the ITS is a minimal requirement
+  // 
+  Bool_t survivedbackground = kTRUE;
+
+  if(!fHFEBackgroundCuts->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *) track)) survivedbackground = kFALSE;
+  AliDebug(3, Form("First cut: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
+  if(!fHFEBackgroundCuts->CheckParticleCuts(AliHFEcuts::kStepRecPrim       + AliHFEcuts::kNcutStepsMCTrack, (TObject *) track)) survivedbackground = kFALSE;
+  AliDebug(3, Form("Second cut: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
+
+  if(survivedbackground){
+    // PID track cuts
+    AliHFEpidObject hfetrack2;
+
+    if(!isAOD) hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
+    else               hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
+
+    hfetrack2.SetRecTrack(track);
+    if(binct>-1){
+     hfetrack2.SetCentrality((Int_t)binct);
+     AliDebug(3,Form("centrality %d and %d",binct,hfetrack2.GetCentrality()));
+     hfetrack2.SetPbPb();
     }
 
-    return -1;
+    if(fPIDBackground->IsSelected(&hfetrack2,0x0,"recTrackCont",fPIDBackgroundQA)){
+      survivedbackground = kTRUE;
+    } else survivedbackground = kFALSE;
   }
+  AliDebug(3, Form("PID: %s\n", survivedbackground == kTRUE ? "yes" : "no"));
+  return survivedbackground;
+}
 
-  return -1;
+//_______________________________________________________________________________________________
+Bool_t AliHFENonPhotonicElectron::FilterCategory2Track(const AliVTrack * const track, Bool_t isAOD){
+  //
+  // Selection of category 2 tracks: These tracks are exclusively low pt tracks below
+  // 300 MeV/c which have at least 2 hits in the 4 outer ITS layers and are identified as
+  // electron candidates by the ITS
+  //
+  if(TMath::Abs(track->Pt()) > 0.3) return kFALSE;
+  if(TMath::Abs(track->Pt()) < 0.1) return kFALSE;
+  Int_t nclustersITS(0), nclustersOuter(0);
+  if(isAOD){
+    const AliAODTrack *aodtrack = static_cast<const AliAODTrack *>(track);
+    if(!(aodtrack->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA) || aodtrack->TestFilterBit(AliAODTrack::kTrkITSsa))) return kFALSE;
+    if(!aodtrack->IsOn(AliAODTrack::kITSrefit)) return kFALSE;
+    nclustersITS = aodtrack->GetITSNcls();
+    for(int ily = 2; ily < 5; ily++)
+      if(aodtrack->HasPointOnITSLayer(ily)) nclustersOuter++;
+  } else {
+    const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track);
+    if(esdtrack->GetStatus() & AliESDtrack::kITSpureSA) return kFALSE;
+    if(esdtrack->GetStatus() & AliESDtrack::kTPCin) return kFALSE;
+    if(!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
+    nclustersITS = esdtrack->GetITSclusters(NULL);
+    for(int ily = 2; ily < 5; ily++)
+      if(esdtrack->HasPointOnITSLayer(ily)) nclustersOuter++;
+  }
+  if(nclustersITS < 4) return kFALSE;
+  if(nclustersOuter < 3) return kFALSE;
+
+  // Do ITS PID
+  Double_t nsigmaITS = fPIDBackground->GetPIDResponse()->NumberOfSigmasITS(track, AliPID::kElectron);
+  fHnsigmaITS->Fill(track->Pt(), nsigmaITS);
+  if(TMath::Abs(nsigmaITS - fITSmeanShift) > 3.) return kFALSE;
+  // if global track, we apply also TPC PID
+  return kTRUE;
 }