Package update (Markus)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Jan 2010 01:28:55 +0000 (01:28 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Jan 2010 01:28:55 +0000 (01:28 +0000)
18 files changed:
PWG3/hfe/AliAnalysisTaskHFE.cxx
PWG3/hfe/AliAnalysisTaskHFE.h
PWG3/hfe/AliHFEcollection.cxx
PWG3/hfe/AliHFEcollection.h
PWG3/hfe/AliHFEcuts.cxx
PWG3/hfe/AliHFEcuts.h
PWG3/hfe/AliHFEmcQA.cxx
PWG3/hfe/AliHFEmcQA.h
PWG3/hfe/AliHFEpid.cxx
PWG3/hfe/AliHFEpid.h
PWG3/hfe/AliHFEpidTOF.cxx
PWG3/hfe/AliHFEpidTOF.h
PWG3/hfe/AliHFEpidTPC.cxx
PWG3/hfe/AliHFEpidTPC.h
PWG3/hfe/AliHFEpriVtx.cxx
PWG3/hfe/AliHFEpriVtx.h
PWG3/hfe/AliHFEsecVtx.cxx
PWG3/hfe/AliHFEsecVtx.h

index acc7a15..2ae22fe 100644 (file)
@@ -43,6 +43,9 @@
 #include <TString.h>
 #include <TTree.h>
 
+#include "AliAODInputHandler.h"
+#include "AliAODMCParticle.h"
+#include "AliAODTrack.h"
 #include "AliCFContainer.h"
 #include "AliCFManager.h"
 #include "AliCFEffGrid.h"
 #include "AliStack.h"
 
 #include "AliHFEpid.h"
+#include "AliHFEcollection.h"
 #include "AliHFEcuts.h"
 #include "AliHFEmcQA.h"
+#include "AliHFEpairs.h"
+#include "AliHFEsecVtxs.h"
 #include "AliHFEsecVtx.h"
+#include "AliHFEelecbackground.h"
 #include "AliAnalysisTaskHFE.h"
 
 //____________________________________________________________
@@ -70,21 +77,26 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fQAlevel(0)
   , fPIDdetectors("")
   , fPIDstrategy(0)
-  , fESD(0x0)
-  , fMC(0x0)
-  , fCFM(0x0)
-  , fCorrelation(0x0)
-  , fPIDperformance(0x0)
-  , fPID(0x0)
-  , fCuts(0x0)
-  , fSecVtx(0x0)
-  , fMCQA(0x0)
-  , fNEvents(0x0)
-  , fNElectronTracksEvent(0x0)
-  , fQA(0x0)
-  , fOutput(0x0)
-  , fHistMCQA(0x0)
-  , fHistSECVTX(0x0)
+  , fPlugins(0)
+  , fRecEvent(NULL)
+  , fMC(NULL)
+  , fCFM(NULL)
+  , fCorrelation(NULL)
+  , fPIDperformance(NULL)
+  , fSignalToBackgroundMC(NULL)
+  , fPID(NULL)
+  , fCuts(NULL)
+  , fSecVtx(NULL)
+  , fElecBackGround(NULL)
+  , fMCQA(NULL)
+  , fNEvents(NULL)
+  , fNElectronTracksEvent(NULL)
+  , fQA(NULL)
+  , fOutput(NULL)
+  , fHistMCQA(NULL)
+  , fHistSECVTX(NULL)
+  , fHistELECBACKGROUND(NULL)
+//  , fQAcoll(0x0)
 {
   //
   // Dummy constructor
@@ -93,11 +105,10 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   DefineOutput(0, TH1I::Class());
   DefineOutput(1, TList::Class());
   DefineOutput(2, TList::Class());
+//  DefineOutput(3, TList::Class());
 
   // Initialize cuts
   fPID = new AliHFEpid;
-
-  SetHasMCData();
 }
 
 //____________________________________________________________
@@ -106,34 +117,38 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fQAlevel(0)
   , fPIDdetectors("")
   , fPIDstrategy(0)
-  , fESD(0x0)
-  , fMC(0x0)
-  , fCFM(0x0)
-  , fCorrelation(0x0)
-  , fPIDperformance(0x0)
-  , fPID(0x0)
-  , fCuts(0x0)
-  , fSecVtx(0x0)
-  , fMCQA(0x0)
-  , fNEvents(0x0)
-  , fNElectronTracksEvent(0x0)
-  , fQA(0x0)
-  , fOutput(0x0)
-  , fHistMCQA(0x0)
-  , fHistSECVTX(0x0)
+  , fPlugins(0)
+  , fRecEvent(NULL)
+  , fMC(NULL)
+  , fCFM(NULL)
+  , fCorrelation(NULL)
+  , fPIDperformance(NULL)
+  , fSignalToBackgroundMC(NULL)
+  , fPID(NULL)
+  , fCuts(NULL)
+  , fSecVtx(NULL)
+  , fElecBackGround(NULL)
+  , fMCQA(NULL)
+  , fNEvents(NULL)
+  , fNElectronTracksEvent(NULL)
+  , fQA(NULL)
+  , fOutput(NULL)
+  , fHistMCQA(NULL)
+  , fHistSECVTX(NULL)
+  , fHistELECBACKGROUND(NULL)
+//  , fQAcoll(0x0)
 {
   //
   // Default constructor
-  //
+  // 
   DefineInput(0, TChain::Class());
   DefineOutput(0, TH1I::Class());
   DefineOutput(1, TList::Class());
   DefineOutput(2, TList::Class());
+//  DefineOutput(3, TList::Class());
 
   // Initialize cuts
   fPID = new AliHFEpid;
-
-  SetHasMCData();
 }
 
 //____________________________________________________________
@@ -142,14 +157,17 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fQAlevel(ref.fQAlevel)
   , fPIDdetectors(ref.fPIDdetectors)
   , fPIDstrategy(ref.fPIDstrategy)
-  , fESD(ref.fESD)
+  , fPlugins(ref.fPlugins)
+  , fRecEvent(ref.fRecEvent)
   , fMC(ref.fMC)
   , fCFM(ref.fCFM)
   , fCorrelation(ref.fCorrelation)
   , fPIDperformance(ref.fPIDperformance)
+  , fSignalToBackgroundMC(ref.fSignalToBackgroundMC)
   , fPID(ref.fPID)
   , fCuts(ref.fCuts)
   , fSecVtx(ref.fSecVtx)
+  , fElecBackGround(ref.fElecBackGround)
   , fMCQA(ref.fMCQA)
   , fNEvents(ref.fNEvents)
   , fNElectronTracksEvent(ref.fNElectronTracksEvent)
@@ -157,6 +175,8 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fOutput(ref.fOutput)
   , fHistMCQA(ref.fHistMCQA)
   , fHistSECVTX(ref.fHistSECVTX)
+  , fHistELECBACKGROUND(ref.fHistELECBACKGROUND)
+//  , fQAcoll(ref.fQAcoll)
 {
   //
   // Copy Constructor
@@ -173,14 +193,17 @@ AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref)
   fQAlevel = ref.fQAlevel;
   fPIDdetectors = ref.fPIDdetectors;
   fPIDstrategy = ref.fPIDstrategy;
-  fESD = ref.fESD;
+  fPlugins = ref.fPlugins;
+  fRecEvent = ref.fRecEvent;
   fMC = ref.fMC;
   fCFM = ref.fCFM;
   fCorrelation = ref.fCorrelation;
   fPIDperformance = ref.fPIDperformance;
+  fSignalToBackgroundMC = ref.fSignalToBackgroundMC;
   fPID = ref.fPID;
   fCuts = ref.fCuts;
   fSecVtx = ref.fSecVtx;
+  fElecBackGround = ref.fElecBackGround;
   fMCQA = ref.fMCQA;
   fNEvents = ref.fNEvents;
   fNElectronTracksEvent = ref.fNElectronTracksEvent;
@@ -188,6 +211,9 @@ AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref)
   fOutput = ref.fOutput;
   fHistMCQA = ref.fHistMCQA;
   fHistSECVTX = ref.fHistSECVTX;
+  fHistELECBACKGROUND = ref.fHistELECBACKGROUND;
+  
+//  fQAcoll = ref.fQAcoll;
   return *this;
 }
 
@@ -196,7 +222,7 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
   //
   // Destructor
   //
-  if(fESD) delete fESD;
+  if(fRecEvent) delete fRecEvent;
   if(fMC) delete fMC;
   if(fPID) delete fPID;
   if(fQA){
@@ -215,7 +241,12 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
     fHistSECVTX->Clear();
     delete fHistSECVTX;
   }
+  if(fHistELECBACKGROUND){
+    fHistELECBACKGROUND->Clear();
+    delete fHistELECBACKGROUND;
+  }
   if(fSecVtx) delete fSecVtx;
+  if(fElecBackGround) delete fElecBackGround;
   if(fMCQA) delete fMCQA;
   if(fNEvents) delete fNEvents;
   if(fCorrelation){
@@ -223,6 +254,8 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
     delete fCorrelation;
   }
   if(fPIDperformance) delete fPIDperformance;
+  if(fSignalToBackgroundMC) delete fSignalToBackgroundMC;
+//  if(fQAcoll) delete fQAcoll;
 }
 
 //____________________________________________________________
@@ -239,19 +272,31 @@ void AliAnalysisTaskHFE::ConnectInputData(Option_t *){
     esdchain->SetBranchStatus("Tracks", 1);
   }
 */
-  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
-  if(!esdH){      
-    AliError("No ESD input handler");
+  AliDebug(3, "Connecting input data");
+  AliVEventHandler *inputH = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+  if(!inputH){
+    AliError("Analysis Manager does not have any input event handler");
     return;
-  } else {
-    fESD = esdH->GetEvent();
   }
-  AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-  if(!mcH){       
-    AliError("No MC truth handler");
-    return;
+  if(IsAODanalysis()){
+    AliInfo("Analysis Mode: AOD Analysis");
+    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler *>(inputH);
+    fRecEvent = aodH->GetEvent();
+    if(fRecEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName()))
+      fMC = aodH->MCEvent();
   } else {
-    fMC = mcH->MCEvent();
+    AliInfo("Analysis Mode: ESD Analysis");
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(inputH);
+    fRecEvent = esdH->GetEvent();
+    if(HasMCData()){
+      AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+      if(!mcH){       
+        AliError("No MC truth handler");
+        return;
+      } else {
+        fMC = mcH->MCEvent();
+      }
+    }
   }
 }
 
@@ -267,6 +312,15 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){
   // QA histograms are created if requested
   // Called once per worker
   //
+  AliDebug(3, "Creating Output Objects");
+  printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
+  printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
+
+  // example how to use the AliHFEcollection
+  //fQAcoll = new AliHFEcollection("fQAcoll", "QA");
+  //fQAcoll->CreateTH1F("fNevents", "Number of Events in the Analysis", 2, 0, 2);
+  //fQAcoll->CreateProfile("fNtrdclusters", "Number of TRD clusters as function of momentum; p[GeV/c]", 20, 0, 20);
+
   fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
   fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
   // First Step: TRD alone
@@ -283,14 +337,15 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){
   // Initialize correction Framework and Cuts
   fCFM = new AliCFManager;
   MakeParticleContainer();
-  // Temporary fix: Initialize particle cuts with 0x0
+  // Temporary fix: Initialize particle cuts with NULL
   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
-    fCFM->SetParticleCutsList(istep, 0x0);
+    fCFM->SetParticleCutsList(istep, NULL);
   if(!fCuts){
     AliWarning("Cuts not available. Default cuts will be used");
     fCuts = new AliHFEcuts;
     fCuts->CreateStandardCuts();
   }
+  if(IsAODanalysis()) fCuts->SetAOD();
   fCuts->Initialize(fCFM);
   if(fCuts->IsInDebugMode()) fQA->Add(fCuts->GetQAhistograms());
  
@@ -298,7 +353,8 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){
   fOutput->AddAt(fCFM->GetParticleContainer(), 0);
   fOutput->AddAt(fCorrelation, 1);
   fOutput->AddAt(fPIDperformance, 2);
-  fOutput->AddAt(fNElectronTracksEvent, 3);
+  fOutput->AddAt(fSignalToBackgroundMC, 3);
+  fOutput->AddAt(fNElectronTracksEvent, 4);
 
   // Initialize PID
   if(IsQAOn(kPIDqa)){
@@ -315,7 +371,7 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){
     fPID->InitializePID(fPIDdetectors.Data());     // Only restrictions to TPC allowed 
 
   // mcQA----------------------------------
-  if (IsQAOn(kMCqa)) {
+  if (HasMCData() && IsQAOn(kMCqa)) {
     AliInfo("MC QA on");
     if(!fMCQA) fMCQA = new AliHFEmcQA;
     if(!fHistMCQA) fHistMCQA = new TList();
@@ -346,27 +402,26 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){
   } 
 
   // secvtx----------------------------------
-  if (IsSecVtxOn()) {
+  if (GetPlugin(kSecVtx)) {
     AliInfo("Secondary Vertex Analysis on");
     fSecVtx = new AliHFEsecVtx;
+    fSecVtx->SetHasMCData(HasMCData());
 
     if(!fHistSECVTX) fHistSECVTX = new TList();
-    fHistSECVTX->SetName("SecVtx");
-    fSecVtx->CreateHistograms("secvtx_");
-    TIter next(gDirectory->GetList());
-    TObject *obj;
-    int counter = 0;
-    TString objname;
-    while ((obj = next.Next())) {
-      objname = obj->GetName();
-      TObjArray *toks = objname.Tokenize("_");
-      if (toks->GetEntriesFast()){
-        TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
-        if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj, counter++);
-      }
-    }
+    fSecVtx->CreateHistograms(fHistSECVTX);
     fOutput->Add(fHistSECVTX);
-  } 
+  }
+  
+  // background----------------------------------
+  if (GetPlugin(kIsElecBackGround)) {
+    AliInfo("Electron BackGround Analysis on");
+    fElecBackGround = new AliHFEelecbackground;
+    fElecBackGround->SetHasMCData(HasMCData());
+
+    if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
+    fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
+    fOutput->Add(fHistELECBACKGROUND);
+  }  
 }
 
 //____________________________________________________________
@@ -374,96 +429,177 @@ void AliAnalysisTaskHFE::Exec(Option_t *){
   //
   // Run the analysis
   // 
-  if(!fESD){
-    AliError("No ESD Event");
+  AliDebug(3, "Starting Single Event Analysis");
+  if(!fRecEvent){
+    AliError("Reconstructed Event not available");
     return;
   }
-  if(!fMC){
-    AliError("No MC Event");
-    return;
+  if(HasMCData()){
+    AliDebug(4, Form("MC Event: %p", fMC));
+    if(!fMC){
+      AliError("No MC Event, but MC Data required");
+      return;
+    }
   }
   if(!fCuts){
     AliError("HFE cuts not available");
     return;
   }
 
-  //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
-  Double_t container[6];
-  // container for the output THnSparse
-  Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
-
-  // Loop over the Monte Carlo tracks to see whether we have overlooked any track
-  AliMCParticle *mctrack = 0x0;
-  Int_t nElectrons = 0;
+  if(HasMCData()) ProcessMC();  // Run the MC loop + MC QA in case MC Data are available
 
-  if (IsSecVtxOn()) { 
-    fSecVtx->SetEvent(fESD);
-    fSecVtx->SetStack(fMC->Stack());
-  }
+  if(IsAODanalysis()) ProcessAOD();
+  else ProcessESD();
+  // Done!!!
+  PostData(0, fNEvents);
+  PostData(1, fOutput);
+  PostData(2, fQA);
+//  PostData(3, fQAcoll->GetList());
+}
 
-  // run mc QA ------------------------------------------------
-  if (IsQAOn(kMCqa)) {
-    AliDebug(2, "Running MC QA");
-
-    fMCQA->SetStack(fMC->Stack());
-    fMCQA->Init();
-
-    Int_t nMCTracks = fMC->Stack()->GetNtrack();
-
-    // loop over all tracks for decayed electrons
-    for (Int_t igen = 0; igen < nMCTracks; igen++){
-      TParticle* mcpart = fMC->Stack()->Particle(igen);
-      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, 0); // no accept cut
-      fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
-      if (TMath::Abs(mcpart->Eta()) < 0.9) {
-        fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-        fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-      }
-      if (TMath::Abs(GetRapidity(mcpart)) < 0.5) {
-        fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-        fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-      }
+//____________________________________________________________
+void AliAnalysisTaskHFE::Terminate(Option_t *){
+  //
+  // Terminate not implemented at the moment
+  //
+  if(GetPlugin(kPostProcess)){
+    fOutput = dynamic_cast<TList *>(GetOutputData(1));
+    if(!fOutput){
+      AliError("Results not available");
+      return;
     }
-    fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
-    fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
-
-  } // end of MC QA loop
-  // -----------------------------------------------------------------
+    PostProcess();
+  }
+}
 
+//____________________________________________________________
+void AliAnalysisTaskHFE::Load(TString filename){
   //
-  // Loop MC
+  // Load Results into the task
   //
-  fCFM->SetMCEventInfo(fMC);
-  for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
-    mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(imc));
+  fQA = NULL; fOutput = NULL; fNEvents = NULL;
+  TFile *input = TFile::Open(filename.Data());
+  if(!input || input->IsZombie()){
+    AliError("Cannot read file");
+    return;
+  }
+  TH1 *htmp = dynamic_cast<TH1I *>(input->Get("nEvents"));
+  if(htmp)
+    fNEvents = dynamic_cast<TH1I *>(htmp->Clone());
+  else
+    AliError("Event Counter histogram not found"); 
+  TList *ltmp = dynamic_cast<TList *>(input->Get("Results"));
+  if(ltmp)
+    fOutput = dynamic_cast<TList *>(ltmp->Clone());
+  else
+    AliError("Output Histograms not found");
+  ltmp = dynamic_cast<TList *>(input->Get("QA"));
+  if(ltmp)
+    fQA = dynamic_cast<TList *>(ltmp->Clone());
+  else
+    AliError("QA histograms not found");
+  input->Close();
+  delete input;
+  Int_t nObjects = 0;
+  if(fNEvents) nObjects++;  
+  if(fOutput) nObjects++;
+  if(fQA) nObjects++;
+  AliInfo(Form("Loaded %d Objects into task", nObjects));
+}
 
-    container[0] = mctrack->Pt();
-    container[1] = mctrack->Eta();
-    container[2] = mctrack->Phi();
+//____________________________________________________________
+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");
+  Int_t nElectrons = 0;
+  if(IsESDanalysis()){
+    if (HasMCData() && IsQAOn(kMCqa)) {
+      AliDebug(2, "Running MC QA");
+
+      fMCQA->SetStack(fMC->Stack());
+      fMCQA->SetGenEventHeader(fMC->GenEventHeader());
+      fMCQA->Init();
+
+      Int_t nMCTracks = fMC->Stack()->GetNtrack();
+
+      // loop over all tracks for decayed electrons
+      for (Int_t igen = 0; igen < nMCTracks; igen++){
+        TParticle* mcpart = fMC->Stack()->Particle(igen);
+        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, 0); // no accept cut
+        fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
+        if (TMath::Abs(mcpart->Eta()) < 0.9) {
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
+          }
+        if (TMath::Abs(GetRapidity(mcpart)) < 0.5) {
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
+        }
+      }
+      fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
+      fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 
-    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
-    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
-    if(IsSignalElectron(mctrack)) fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCsignal);
-    (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
-    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
-    // find the label in the vector
-    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
-    nElectrons++;
+    } // end of MC QA loop
+    // -----------------------------------------------------------------
+    fCFM->SetMCEventInfo(fMC);
+    // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
+  } else {
+    fCFM->SetMCEventInfo(fRecEvent);
+  }
+  // Run MC loop
+  for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
+    if(ProcessMCtrack(fMC->GetTrack(imc))) nElectrons++;
   }
-  (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
 
   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
-        
+  (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
+}
 
+//____________________________________________________________
+void AliAnalysisTaskHFE::ProcessESD(){
   //
-  // Loop ESD
+  // 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 *>(fRecEvent);
+  if(!fESD){
+    AliError("ESD Event required for ESD Analysis")
+    return;
+  }
+  if (GetPlugin(kIsElecBackGround)) { 
+    fElecBackGround->SetEvent(fESD);
+  }
+  if (GetPlugin(kSecVtx)) {
+    fSecVtx->SetEvent(fESD);
+    fSecVtx->GetPrimaryCondition();
+  }
+
+  if(HasMCData()){
+    if (GetPlugin(kSecVtx)) { 
+      fSecVtx->SetStack(fMC->Stack());
+    }
+    if (GetPlugin(kIsElecBackGround)) { 
+      fElecBackGround->SetMCEvent(fMC);
+    }
+  }
+
+
+  Double_t container[6];
+  memset(container, 0, sizeof(Double_t) * 6);
+  // container for the output THnSparse
+  Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
   Int_t nElectronCandidates = 0;
-  AliESDtrack *track = 0x0, *htrack = 0x0;
+  AliESDtrack *track = NULL, *htrack = NULL;
+  AliMCParticle *mctrack = NULL;
+  TParticle* mctrack4QA = NULL;
   Int_t pid = 0;
   // For double counted tracks
   LabelContainer cont(fESD->GetNumberOfTracks());
@@ -472,6 +608,9 @@ void AliAnalysisTaskHFE::Exec(Option_t *){
   Bool_t signal = kTRUE;
 
   fCFM->SetRecEventInfo(fESD);
+  //
+  // Loop ESD
+  //
   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
     
     track = fESD->GetTrack(itrack);
@@ -491,15 +630,17 @@ void AliAnalysisTaskHFE::Exec(Option_t *){
     // RecKine: ITSTPC cuts  
     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
     
-    // Check if it is electrons near the vertex
-    if(!(mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
-    TParticle* mctrack4QA = fMC->Stack()->Particle(TMath::Abs(track->GetLabel()));
-
-    container[3] = mctrack->Pt();
-    container[4] = mctrack->Eta();
-    container[5] = mctrack->Phi();
+    if(HasMCData()){
+      // Check if it is electrons near the vertex
+      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
+      mctrack4QA = fMC->Stack()->Particle(TMath::Abs(track->GetLabel()));
+
+      container[3] = mctrack->Pt();
+      container[4] = mctrack->Eta();
+      container[5] = mctrack->Phi();
     
-    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
+      if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
+    }
     
     if(signal) {
       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));
@@ -518,47 +659,28 @@ void AliAnalysisTaskHFE::Exec(Option_t *){
       (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha());    // Check the acceptance without tight cuts
       (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
       (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
+      //fQAcoll->Fill("fNtrdclusters", container[0], track->GetTRDncls());
     }
 
     
     // RecPrim
-    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
-    if(signal) {
-      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1));
-      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim);
-      if(alreadyseen) {
-        fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1);
-      }
-    }
-
+    if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track, container, signal, alreadyseen)) continue;
 
     // HFEcuts: ITS layers cuts
-    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
-    if(signal) {
-      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1));
-      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS);
-      if(alreadyseen) {
-        fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1);
-      }
-    }
+    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track, container, signal, alreadyseen)) continue;
 
     // HFEcuts: Nb of tracklets TRD0
-    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
+    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track, container, signal, alreadyseen)) continue;
     if(signal) {
-      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1));
-      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD);
-      if(alreadyseen) {
-        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + (AliHFEcuts::kNcutESDSteps + 1)));
-      }
       // dimensions 3&4&5 : pt,eta,phi (MC)
       ((THnSparseF *)fCorrelation->At(0))->Fill(container);
     }
 
-   if (IsQAOn(kMCqa)) {
-     // mc qa for after the reconstruction cuts  
-     AliDebug(2, "Running MC QA");
-     fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 3);  // charm
-     fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 3); // beauty 
+    if(HasMCData() && IsQAOn(kMCqa)) {
+      // mc qa for after the reconstruction cuts  
+      AliDebug(2, "Running MC QA");
+      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 3);  // charm
+      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 3); // beauty 
     } 
 
     // track accepted, do PID
@@ -569,11 +691,11 @@ void AliAnalysisTaskHFE::Exec(Option_t *){
     if(!fPID->IsSelected(&hfetrack)) continue;
     nElectronCandidates++;
 
-   if (IsQAOn(kMCqa)) {
-     // mc qa for after the reconstruction and pid cuts  
-     AliDebug(2, "Running MC QA");
-     fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 4);  // charm
-     fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 4); // beauty 
+    if (HasMCData() && IsQAOn(kMCqa)) {
+      // mc qa for after the reconstruction and pid cuts  
+      AliDebug(2, "Running MC QA");
+      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 4);  // charm
+      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 4); // beauty 
     } 
 
     // Fill Containers
@@ -587,104 +709,219 @@ void AliAnalysisTaskHFE::Exec(Option_t *){
       ((THnSparseF *)fCorrelation->At(1))->Fill(container);
     }
 
-    // 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 = IsSignalElectron(track);
-      AliDebug(1, Form("Type: %d\n", type));
-      if(type){
-             dataE[4] = type; // beauty[1] or charm[2]
-             dataE[3] = 2;  // signal electron
-      }
-      else{
-             dataE[3] = 1; // not a signal electron
-             dataE[4] = 0;
-      }
-      // pair analysis [mj]
-      if (IsSecVtxOn()) {
-        AliDebug(2, "Running Secondary Vertex Analysis");
-        fSecVtx->InitAnaPair();
+    if(GetPlugin(kSecVtx)) {
+      AliDebug(2, "Running Secondary Vertex Analysis");
+      if(track->Pt()>0.5){
+        fSecVtx->InitHFEpairs();
+        fSecVtx->InitHFEsecvtxs();
+        AliESDtrack *htrack = 0x0; 
         for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
           htrack = fESD->GetTrack(jtrack);
-          if ( itrack == jtrack ) continue;
-          //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue;
-          if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue;
-          fSecVtx->AnaPair(track, htrack, jtrack);
+          if ( itrack == jtrack ) continue; // since it is for tagging single electron, don't need additional condition 
+          if (htrack->Pt()<0.5) continue;
+          if (!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, htrack)) continue;
+          if (!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, htrack)) continue;
+          fSecVtx->PairAnalysis(track, htrack, jtrack); // e-h pairing
         }
-        // based on the partner of e info, you run secandary vertexing function
-        fSecVtx->RunSECVTX(track);
-      } // end of pair analysis
-    } 
-    else {
-      // Fill THnSparse with the information for Fake Electrons
-      dataE[3] = 0;
-      dataE[4] = 0;
+        /*for(int ip=0; ip<fSecVtx->HFEpairs()->GetEntriesFast(); ip++){
+          AliHFEpairs *pair = (AliHFEpairs*) (fSecVtx->HFEpairs()->UncheckedAt(ip));
+          if(HasMCData()){
+            if(!(pair->GetPairCode()>1. && pair->GetPairCode()<4.))  // apply various cuts
+              fSecVtx->HFEpairs()->RemoveAt(ip);
+          }
+        }*/
+        fSecVtx->HFEpairs()->Compress();
+        fSecVtx->RunSECVTX(track); // secondary vertexing with e,h1,h2,.. tracks
+        for(int ip=0; ip<fSecVtx->HFEsecvtxs()->GetEntriesFast(); ip++){
+          AliHFEsecVtxs *secvtx=0x0;
+          secvtx = (AliHFEsecVtxs*) (fSecVtx->HFEsecvtxs()->UncheckedAt(ip));
+          // here you apply cuts, then if it doesn't pass the cut, remove it from the fSecVtx->HFEsecvtxs() 
+        }
+        fSecVtx->DeleteHFEpairs();
+        fSecVtx->DeleteHFEsecvtxs();
+      }
     }
-    // fill the performance THnSparse, if the mc origin could be defined
-    if(dataE[3] > -1){
-      AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4]));
-      fPIDperformance->Fill(dataE);
+
+    if(HasMCData()){
+      // 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 = IsSignalElectron(track);
+        AliDebug(1, Form("Type: %d\n", type));
+        if(type){
+               dataE[4] = type; // beauty[1] or charm[2]
+               dataE[3] = 2;  // signal electron
+        }
+        else{
+               dataE[3] = 1; // not a signal electron
+               dataE[4] = 0;
+        }
+      } 
+      else {
+        // Fill THnSparse with the information for Fake Electrons
+        dataE[3] = 0;
+        dataE[4] = 0;
+      }
+      // fill the performance THnSparse, if the mc origin could be defined
+      if(dataE[3] > -1){
+        AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4]));
+        fPIDperformance->Fill(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
+
   }
-  
   fNEvents->Fill(1);
+  //fQAcoll->Fill("fNevents", 1);
   fNElectronTracksEvent->Fill(nElectronCandidates);
-  
-  // Done!!!
-  PostData(0, fNEvents);
-  PostData(1, fOutput);
-  PostData(2, fQA);
 }
 
 //____________________________________________________________
-void AliAnalysisTaskHFE::Terminate(Option_t *){
+void AliAnalysisTaskHFE::ProcessAOD(){
   //
-  // Terminate not implemented at the moment
+  // Run Analysis in AOD Mode
+  // Function is still in development
   //
-  if(IsRunningPostProcess()){
-    fOutput = dynamic_cast<TList *>(GetOutputData(1));
-    if(!fOutput){
-      AliError("Results not available");
-      return;
+  AliDebug(3, "Processing AOD Event");
+  AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fRecEvent);
+  if(!fAOD){
+    AliError("AOD Event required for AOD Analysis")
+    return;
+  }
+  AliAODTrack *track = NULL;
+  AliAODMCParticle *mctrack = NULL;
+  Double_t container[6]; memset(container, 0, sizeof(Double_t) * 6);
+  Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
+  Int_t nElectronCandidates = 0;
+  Int_t pid;
+  for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
+    track = fAOD->GetTrack(itrack);
+    if(!track) continue;
+    if(track->GetFlags() != 1<<4) continue;  // Only process AOD tracks where the HFE is set
+
+    container[0] = track->Pt();
+    container[1] = track->Eta();
+    container[2] = track->Phi();
+
+    dataE[0] = track->Pt();
+    dataE[1] = track->Eta();
+    dataE[2] = track->Phi();
+    dataE[3] = -1;
+    dataE[4] = -1;
+    
+    if(HasMCData()){
+      Int_t label = TMath::Abs(track->GetLabel());
+      if(label){
+        mctrack = dynamic_cast<AliAODMCParticle *>(fMC->GetTrack(label));
+        container[3] = mctrack->Pt();
+        container[4] = mctrack->Eta();
+        container[5] = mctrack->Phi();
+      }
+    }
+    // track accepted, do PID
+    AliHFEpidObject hfetrack;
+    hfetrack.fAnalysisType = AliHFEpidObject::kAODanalysis;
+    hfetrack.fRecTrack = track;
+    if(HasMCData()) hfetrack.fMCtrack = mctrack;
+    //if(!fPID->IsSelected(&hfetrack)) continue;    // we will do PID here as soon as possible
+    // Particle identified - Fill CF Container
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD +1 + 2*(AliHFEcuts::kNcutESDSteps + 1));
+    nElectronCandidates++;    
+    if(HasMCData()){
+      // Track selected: distinguish between true and fake
+      AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->GetPdgCode()));
+      if((pid = TMath::Abs(mctrack->GetPdgCode())) == 11){
+        Int_t type = IsSignalElectron(track);
+        AliDebug(1, Form("Type: %d\n", type));
+        if(type){
+               dataE[4] = type; // beauty[1] or charm[2]
+               dataE[3] = 2;  // signal electron
+        }
+        else{
+               dataE[3] = 1; // not a signal electron
+               dataE[4] = 0;
+        }
+      } 
+      else {
+        // Fill THnSparse with the information for Fake Electrons
+        dataE[3] = 0;
+        dataE[4] = 0;
+      }
+      // fill the performance THnSparse, if the mc origin could be defined
+      if(dataE[3] > -1){
+        AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4]));
+        fPIDperformance->Fill(dataE);
+      }
     }
-    PostProcess();
   }
+  fNEvents->Fill(1);
+  fNElectronTracksEvent->Fill(nElectronCandidates);
 }
 
 //____________________________________________________________
-void AliAnalysisTaskHFE::Load(TString filename){
+Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
   //
-  // Load Results into the task
+  // Filter the Monte Carlo Track
+  // Additionally Fill a THnSparse for Signal To Background Studies
+  // Works for AOD and MC analysis Type
   //
-  fQA = NULL; fOutput = NULL; fNEvents = NULL;
-  TFile *input = TFile::Open(filename.Data());
-  if(!input || input->IsZombie()){
-    AliError("Cannot read file");
-    return;
+  Double_t container[3], signalContainer[5];
+  Double_t vertex[3]; // 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);
+    container[0] = mctrack->Pt();
+    container[1] = mctrack->Eta();
+    container[2] = mctrack->Phi();
+
+    signalContainer[0] = mctrack->Pt();
+    signalContainer[1] = mctrack->Eta();
+    signalContainer[2] = mctrack->Phi();
+
+    vertex[0] = mctrack->Particle()->Vx();
+    vertex[1] = mctrack->Particle()->Vy();
+  } else {
+    AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
+    container[0] = aodmctrack->Pt();
+    container[1] = aodmctrack->Eta();
+    container[2] = aodmctrack->Phi();
+
+    signalContainer[0] = aodmctrack->Pt();
+    signalContainer[1] = aodmctrack->Eta();
+    signalContainer[2] = aodmctrack->Phi();
+
+    aodmctrack->XvYvZv(vertex);
   }
-  TH1 *htmp = dynamic_cast<TH1I *>(input->Get("nEvents"));
-  if(htmp)
-    fNEvents = dynamic_cast<TH1I *>(htmp->Clone());
-  else
-    AliError("Event Counter histogram not found"); 
-  TList *ltmp = dynamic_cast<TList *>(input->Get("Results"));
-  if(ltmp)
-    fOutput = dynamic_cast<TList *>(ltmp->Clone());
-  else
-    AliError("Output Histograms not found");
-  ltmp = dynamic_cast<TList *>(input->Get("QA"));
-  if(ltmp)
-    fQA = dynamic_cast<TList *>(ltmp->Clone());
-  else
-    AliError("QA histograms not found");
-  input->Close();
-  delete input;
-  Int_t nObjects = 0;
-  if(fNEvents) nObjects++;  
-  if(fOutput) nObjects++;
-  if(fQA) nObjects++;
-  AliInfo(Form("Loaded %d Objects into task", nObjects));
+
+  if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
+  fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
+  if((signalContainer[3] = static_cast<Double_t >(IsSignalElectron(track))) > 1e-3) fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCsignal);
+  signalContainer[4] = 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[4] = 1;
+  } else if (radVertex < 7.5){
+    signalContainer[4] = 2;
+  }
+  fSignalToBackgroundMC->Fill(signalContainer);
+  (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(container[2] - TMath::Pi());
+  //if(IsESDanalysis()){
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
+  //}
+  return kTRUE;
 }
 
 //____________________________________________________________
@@ -879,6 +1116,63 @@ void AliAnalysisTaskHFE::PostProcess(){
 }
 
 //____________________________________________________________
+void AliAnalysisTaskHFE::DrawMCSignal2Background(){
+  //
+  // Draw the MC signal/background plots
+  //
+  fSignalToBackgroundMC = dynamic_cast<THnSparseF *>(fOutput->FindObject("SignalToBackgroundMC"));
+  if(!fSignalToBackgroundMC) return;
+
+  // First Select everything within the first ITS Layer
+  fSignalToBackgroundMC->GetAxis(4)->SetRange(2,2);
+  // For Signal electrons we project axis 3 to everything > 0
+  fSignalToBackgroundMC->GetAxis(3)->SetRange(2,3);
+  TH1 *hSignal = fSignalToBackgroundMC->Projection(0);
+  hSignal->SetName("hSignal");
+  hSignal->SetTitle("Signal Electrons");
+  // For Background studies project axis 3 to bin 0
+  fSignalToBackgroundMC->GetAxis(3)->SetRange(1,1);
+  TH1 *hBackground = fSignalToBackgroundMC->Projection(0); 
+  hBackground->SetName("hBackground");
+  hBackground->SetTitle("Background Electrons");
+  // For All electrons we undo the projection of axis 3
+  fSignalToBackgroundMC->GetAxis(3)->SetRange(0, fSignalToBackgroundMC->GetAxis(3)->GetNbins());
+  TH1 *hAll = fSignalToBackgroundMC->Projection(0);
+  hAll->SetName("hAll");
+  hAll->SetTitle("All Electrons");
+  // Prepare Canvas
+  TCanvas *cEff = new TCanvas("cEff", "MC Sig/Backg studies", 800, 400);
+  cEff->Divide(2);
+  // Plot Efficiency
+  TH1 *hEff = (TH1 *)hSignal->Clone();
+  hEff->Divide(hAll);
+  hEff->SetName("sigEff");
+  hEff->SetTitle("Signal/(Signal+Background)");
+  hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
+  hEff->GetYaxis()->SetTitle("Efficiency");
+  hEff->SetStats(kFALSE);
+  hEff->SetMarkerStyle(22);
+  hEff->SetMarkerColor(kBlue);
+  cEff->cd(1);
+  hEff->Draw("p");
+  // Plot Signal/Background
+  TH1 *hSB = (TH1 *)hSignal->Clone();
+  hSB->Divide(hBackground);
+  hSB->SetName("sigEff");
+  hSB->SetTitle("Signal/Background");
+  hSB->GetXaxis()->SetTitle("p_{T} / GeV/c");
+  hSB->GetYaxis()->SetTitle("Signal/Background");
+  hSB->SetStats(kFALSE);
+  hSB->SetMarkerStyle(22);
+  hSB->SetMarkerColor(kBlue);
+  cEff->cd(2);
+  hSB->Draw("p");
+
+  // Undo projections
+  fSignalToBackgroundMC->GetAxis(4)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetNbins());
+}
+
+//____________________________________________________________
 void AliAnalysisTaskHFE::MakeParticleContainer(){
   //
   // Create the particle container for the correction framework manager and 
@@ -953,9 +1247,12 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
   for(Int_t i=0; i<=nBin[3]; i++) binEdges2[3][i] = i;
   for(Int_t i=0; i<=nBin[4]; i++) binEdges2[4][i] = i;
 
-  fPIDperformance = new THnSparseF("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad] type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin);
-  for(Int_t idim = 0; idim < nDim; idim++)
+  fPIDperformance = new THnSparseF("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin);
+  fSignalToBackgroundMC = new THnSparseF("SignalToBackgroundMC", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; flavor (0 - no, 1 - charm, 2 - bottom); ITS Cluster (0 - no, 1 - first (and maybe second), 2 - second)", nDim, nBin);
+  for(Int_t idim = 0; idim < nDim; idim++){
     fPIDperformance->SetBinEdges(idim, binEdges2[idim]);
+    fSignalToBackgroundMC->SetBinEdges(idim, binEdges2[idim]); 
+  }
 }
 
 //____________________________________________________________
@@ -975,8 +1272,8 @@ void AliAnalysisTaskHFE::PrintStatus() const {
   // Print Analysis status
   //
   printf("\n\tAnalysis Settings\n\t========================================\n\n");
-  printf("\tSecondary Vertex finding: %s\n", IsSecVtxOn() ? "YES" : "NO");
-  printf("\tPrimary Vertex resolution: %s\n", IsPriVtxOn() ? "YES" : "NO");
+  printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
+  printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
   printf("\n");
   printf("\tParticle Identification Detectors:\n");
   TObjArray *detectors = fPIDdetectors.Tokenize(":");
@@ -1047,35 +1344,56 @@ Int_t AliAnalysisTaskHFE::IsSignalElectron(AliVParticle *fTrack) const{
     kCharm = 1,
     kBeauty = 2
   };
-  AliMCParticle *mctrack = NULL;
   TString objname = fTrack->IsA()->GetName();
-  if(!objname.CompareTo("AliESDtrack")){
-    AliDebug(2, "Checking signal for ESD track");
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(fTrack);
-    mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(esdtrack->GetLabel())));
-  }
-  else if(!objname.CompareTo("AliAODtrack")){
-    AliError("AOD Analysis not implemented yet");
-    return kNoSignal;
-  }
-  else if(!objname.CompareTo("AliMCParticle")){
-    AliDebug(2, "Checking signal for MC track");
-    mctrack = dynamic_cast<AliMCParticle *>(fTrack);
-  }
-  else{
-    AliError("Input object not supported");
-    return kNoSignal;
+  Int_t pid = 0;
+  if(IsESDanalysis()){
+    // ESD Analysis
+    AliMCParticle *mctrack = NULL;
+    if(!objname.CompareTo("AliESDtrack")){
+      AliDebug(2, "Checking signal for ESD track");
+      AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(fTrack);
+      mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(esdtrack->GetLabel())));
+    }
+    else if(!objname.CompareTo("AliMCParticle")){
+      AliDebug(2, "Checking signal for MC track");
+      mctrack = dynamic_cast<AliMCParticle *>(fTrack);
+    }
+    else{
+      AliError("Input object not supported");
+      return kNoSignal;
+    }
+    if(!mctrack) return kNoSignal;
+    TParticle *ecand = mctrack->Particle(); 
+    if(TMath::Abs(ecand->GetPdgCode()) != 11) return kNoSignal; // electron candidate not true electron
+    Int_t motherLabel = TMath::Abs(ecand->GetFirstMother());
+    AliDebug(3, Form("mother label: %d\n", motherLabel));
+    if(!motherLabel) return kNoSignal; // mother track unknown
+    AliMCParticle *motherTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(motherLabel));
+    if(!motherTrack) return kNoSignal;
+    TParticle *mparticle = motherTrack->Particle();
+    pid = TMath::Abs(mparticle->GetPdgCode());
+  } else {
+    // AOD Analysis - Different Data handling
+    AliAODMCParticle *aodmc = NULL;
+    if(!objname.CompareTo("AliAODTrack")){
+      AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(fTrack);
+      Int_t aodlabel = TMath::Abs(aodtrack->GetLabel());
+      if(aodlabel >= fMC->GetNumberOfTracks()) return kNoSignal;
+      aodmc = dynamic_cast<AliAODMCParticle *>(fMC->GetTrack(aodlabel));
+    } else if(!objname.CompareTo("AliAODMCParticle")){
+      aodmc = dynamic_cast<AliAODMCParticle *>(fTrack);
+    } else{
+      AliError("Input object not supported");
+      return kNoSignal;
+    }
+    if(!aodmc) return kNoSignal;
+    Int_t motherLabel = TMath::Abs(aodmc->GetMother());
+    AliDebug(3, Form("mother label: %d\n", motherLabel));
+    if(!motherLabel || motherLabel >= fMC->GetNumberOfTracks()) return kNoSignal;
+    AliAODMCParticle *aodmother = dynamic_cast<AliAODMCParticle *>(fMC->GetTrack(motherLabel));
+    pid = aodmother->GetPdgCode();
   }
-  if(!mctrack) return kNoSignal;
-  TParticle *ecand = mctrack->Particle(); 
-  if(TMath::Abs(ecand->GetPdgCode()) != 11) return kNoSignal; // electron candidate not true electron
-  Int_t motherLabel = TMath::Abs(ecand->GetFirstMother());
-  AliDebug(3, Form("mother label: %d\n", motherLabel));
-  if(!motherLabel) return kNoSignal; // mother track unknown
-  AliMCParticle *motherTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(motherLabel));
-  if(!motherTrack) return kNoSignal;
-  TParticle *mparticle = motherTrack->Particle();
-  Int_t pid = TMath::Abs(mparticle->GetPdgCode());
+  // From here the two analysis modes go together
   AliDebug(3, Form("PDG code: %d\n", pid));
 
   // identify signal according to Pdg Code 
@@ -1098,3 +1416,37 @@ Float_t AliAnalysisTaskHFE::GetRapidity(TParticle *part) const
   return rapidity;
 }
 
+//__________________________________________
+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;
+    default: AliError("Unknown Plugin");
+  };
+}
+
+//__________________________________________
+Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track, Double_t *container, Bool_t signal, Bool_t alreadyseen){
+  //
+  // Check single track cuts for a given cut step
+  // Fill the particle container
+  //
+  if(!fCFM->CheckParticleCuts(cutStep, track)) return kFALSE;
+  if(signal) {
+    fCFM->GetParticleContainer()->Fill(container, cutStep + 2*(AliHFEcuts::kNcutESDSteps + 1));
+    fCFM->GetParticleContainer()->Fill(&container[3], cutStep);
+    if(alreadyseen) {
+      fCFM->GetParticleContainer()->Fill(&container[3], cutStep + AliHFEcuts::kNcutESDSteps + 1);
+    }
+  }
+  return kTRUE;
+}
index af5fe54..e07ef9b 100644 (file)
@@ -32,9 +32,10 @@ class AliHFEpid;
 class AliHFEcuts;
 class AliHFEmcQA;
 class AliHFEsecVtx;
+class AliHFEelecbackground;
+class AliHFEcollection;
 class AliCFManager;
-class AliESDEvent;
-class AliESDtrackCuts;
+class AliVEvent;
 class AliMCEvent;
 class AliVParticle;
 class TH1I; 
@@ -42,10 +43,16 @@ class TList;
 
 class AliAnalysisTaskHFE : public AliAnalysisTask{
   public:
-  enum{
-    kPIDqa = 0,
-    kMCqa =1 
-  };
+    enum{
+      kPIDqa = 0,
+      kMCqa =1 
+    };
+    enum{
+      kPriVtx = 0,
+      kSecVtx = 1,
+      kIsElecBackGround = 2,
+      kPostProcess = 3
+    };
     AliAnalysisTaskHFE();
     AliAnalysisTaskHFE(const char * name);
     AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref);
@@ -58,31 +65,29 @@ class AliAnalysisTaskHFE : public AliAnalysisTask{
     virtual void Terminate(Option_t *);
 
     Bool_t IsQAOn(Int_t qaLevel) const { return TESTBIT(fQAlevel, qaLevel); };
-    Bool_t IsSecVtxOn() const { return TestBit(kIsSecVtxOn); };
-    Bool_t IsPriVtxOn() const { return TestBit(kIsPriVtxOn); };
-    Bool_t IsRunningPostProcess() const { return TestBit(kIsRunningPostProcess); };
+    Bool_t IsAODanalysis() const { return TestBit(kAODanalysis); };
+    Bool_t IsESDanalysis() const { return !TestBit(kAODanalysis); };
     Bool_t HasMCData() const { return TestBit(kHasMCdata); }
+    Bool_t GetPlugin(Int_t plug) const { return TESTBIT(fPlugins, plug); };
     Int_t IsSignalElectron(AliVParticle *fTrack) const;
     void Load(TString filename = "HFEtask.root");
     void PostProcess();
     void SetHFECuts(AliHFEcuts * const cuts) { fCuts = cuts; };
     void SetQAOn(Int_t qaLevel) { SETBIT(fQAlevel, qaLevel); };
+    void SwitchOnPlugin(Int_t plug);
     void SetHasMCData(Bool_t hasMC = kTRUE) { SetBit(kHasMCdata, hasMC); };
-    void SetPriVtxOn(Bool_t option = kTRUE)        { SetBit(kIsPriVtxOn, option); };
-    void SetSecVtxOn(Bool_t option = kTRUE)        { SetBit(kIsSecVtxOn, option); };
-    void SetRunPostProcess(Bool_t option = kTRUE)  { SetBit(kIsRunningPostProcess, option); };
     void SetPIDdetectors(Char_t * const detectors){ fPIDdetectors = detectors; }
     void SetPIDStrategy(UInt_t strategy) { fPIDstrategy = strategy; }
     void AddPIDdetector(Char_t *detector);
+    void SetAODAnalysis() { SetBit(kAODanalysis, kTRUE); };
+    void SetESDAnalysis() { SetBit(kAODanalysis, kFALSE); };
     void PrintStatus() const;
     Float_t GetRapidity(TParticle *part) const;
  
   private:
     enum{
-      kIsSecVtxOn = BIT(19),
-      kIsPriVtxOn = BIT(20),
-      kIsRunningPostProcess = BIT(21),
-      kHasMCdata = BIT(22)
+      kHasMCdata = BIT(19),
+      kAODanalysis = BIT(20)
     };
     class LabelContainer{
       public:
@@ -104,18 +109,27 @@ class AliAnalysisTaskHFE : public AliAnalysisTask{
         Int_t *fCurrent;      // Current entry to mimic an iterator
     };
     void MakeParticleContainer();
+    void ProcessMC();
+    void ProcessESD();
+    void ProcessAOD();
+    Bool_t ProcessMCtrack(AliVParticle *track);
+    Bool_t ProcessCutStep(Int_t cutStep, AliVParticle *track, Double_t *container, Bool_t signal, Bool_t alreadyseen);
+    void DrawMCSignal2Background();
     
     ULong_t fQAlevel;                     // QA level
     TString fPIDdetectors;                // Detectors for Particle Identification
     UInt_t fPIDstrategy;                  // PID Strategy
-    AliESDEvent *fESD;                    //! The ESD Event
+    UShort_t fPlugins;                    // Enabled Plugins
+    AliVEvent *fRecEvent;                 //! The ESD Event
     AliMCEvent *fMC;                      //! The MC Event
     AliCFManager *fCFM;                   //! Correction Framework Manager
     TList *fCorrelation;                  //! response matrix for unfolding  
     THnSparseF *fPIDperformance;          //! info on contamination and yield of electron spectra
+    THnSparseF *fSignalToBackgroundMC;    //! Signal To Background Studies on pure MC information
     AliHFEpid *fPID;                      //! PID
     AliHFEcuts *fCuts;                    // Cut Collection
     AliHFEsecVtx *fSecVtx;                //! Secondary Vertex Analysis
+    AliHFEelecbackground *fElecBackGround;//! Background analysis
     AliHFEmcQA *fMCQA;                    //! MC QA
     TH1I *fNEvents;                       //! counter for the number of Events
     TH1I *fNElectronTracksEvent;          //! Number of Electron candidates after PID decision per Event
@@ -123,6 +137,8 @@ class AliAnalysisTaskHFE : public AliAnalysisTask{
     TList *fOutput;                       //! Container for Task Output
     TList *fHistMCQA;                     //! Output container for MC QA histograms 
     TList *fHistSECVTX;                   //! Output container for sec. vertexing results
+    TList *fHistELECBACKGROUND;           //! Output container for electron background analysis
+//    AliHFEcollection *fQAcoll;            //! collection class for basic QA histograms
 
     ClassDef(AliAnalysisTaskHFE, 1)       // The electron Analysis Task
 };
index 65c5d41..03d23f3 100644 (file)
 
 #include <TH1F.h>
 #include <TH2F.h>
+#include <THnSparse.h>
+#include <TProfile.h>
 #include <TList.h>
 #include <TString.h>
 #include <TBrowser.h>
-#include <TIterator.h>
+#include <TMath.h>
 
 #include "AliLog.h"
 #include "AliHFEcollection.h"
@@ -40,20 +42,21 @@ ClassImp(AliHFEcollection)
 //___________________________________________________________________
 AliHFEcollection::AliHFEcollection():
   TNamed()
-  , fListE(0x0)
+  , fList(0x0)
 {
+
   //
   // default constructor
   //
 
-  fListE = new TList();
-  if(!fListE){
+  fList = new TList();
+  if(!fList){
     AliError("Initialization of the list failed");
   }
   else{
     // list is owner of the objects. Once list is deleted, the objects
     // it contains will be deleted too
-    fListE->SetOwner(kTRUE);
+    fList->SetOwner(kTRUE);
   }
   //Printf("%s:%d,%p",(char*)__FILE__,__LINE__,fInstance);
   
@@ -61,27 +64,27 @@ AliHFEcollection::AliHFEcollection():
 //___________________________________________________________________
 AliHFEcollection::AliHFEcollection(char* name, char* title):
   TNamed(name, title)
-  , fListE(0x0)
+  , fList(0x0)
 {
  
   //
   // constructor
   //
  
-  fListE = new TList();
-  if(!fListE){
+  fList = new TList();
+  if(!fList){
     AliError("Initialization of the list failed");
   }
   else{
     // list is owner of the objects. Once list is deleted, the objects
     // it contains will be deleted too
-    fListE->SetOwner(kTRUE);
+    fList->SetOwner(kTRUE);
   }
 }
 //___________________________________________________________________
 AliHFEcollection::AliHFEcollection(const AliHFEcollection &c) :
   TNamed(c)
-  , fListE(0x0)
+  , fList(0x0)
 {
 
   //
@@ -104,52 +107,62 @@ AliHFEcollection &AliHFEcollection::operator=(const AliHFEcollection &ref)
 }
 //___________________________________________________________________
 void AliHFEcollection::Copy(TObject &ref) const {
+
   //
   // Performs the copying of the object
   //
+
   AliHFEcollection &target = dynamic_cast<AliHFEcollection &>(ref);
 
-  target.fListE = fListE;          
+  target.fList = fList;          
 }
 //___________________________________________________________________
 AliHFEcollection::~AliHFEcollection(){
+
   //
   // Destructor
   //
+
   AliInfo("DESTRUCTOR");
 }
 //___________________________________________________________________
 Bool_t AliHFEcollection::CreateTH1F(const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax){
+
   //
   // Creates a TH1F histogram for the collection
   //
-  if(!fListE){
+
+  if(!fList){
     AliError("No TList pointer ! ");
     return kFALSE;
   }
   else{
-    fListE->Add(new TH1F(name, title, nBin, nMin, nMax));
+    fList->Add(new TH1F(name, title, nBin, nMin, nMax));
     return CheckObject(name);
   }
 }
 //___________________________________________________________________
 Bool_t AliHFEcollection::CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY){
+
   //
   // Creates a TH2F histogram for the collection
   //
-  if(!fListE){
+
+  if(!fList){
     AliError("No TList pointer ! ");
     return kFALSE;
   }
-  fListE->Add(new TH2F(name, title, nBinX, nMinX, nMaxX, nBinY, nMinY, nMaxY));
+  fList->Add(new TH2F(name, title, nBinX, nMinX, nMaxX, nBinY, nMinY, nMaxY));
   return CheckObject(name); 
 }
 //___________________________________________________________________
 Bool_t AliHFEcollection::CreateTH1Fvector1(Int_t X, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax){
+
   //
   // create a 1 dimensional array of size [X]
   //
-  if(!fListE){
+
+  if(!fList){
     AliError("No TList pointer ! ");
     return kFALSE;
   }
@@ -173,10 +186,12 @@ Bool_t AliHFEcollection::CreateTH1Fvector1(Int_t X, const char* name, const char
 }
 //___________________________________________________________________
 Bool_t AliHFEcollection::CreateTH2Fvector1(Int_t X, const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY){
+
   //
   // create a 1 dimensinal array of TH2F histograms with size [X]
   //
-  if(!fListE){
+
+  if(!fList){
     AliError("No TList pointer !");
     return kFALSE;
   }
@@ -200,10 +215,12 @@ Bool_t AliHFEcollection::CreateTH2Fvector1(Int_t X, const char* name, const char
 }
 //___________________________________________________________________
 Bool_t AliHFEcollection::CreateTH1Fvector2(Int_t X, Int_t Y, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax){
+
   //
   // create a 2 dimensional array of histograms of size [X, Y]
   //
-  if(!fListE){
+
+  if(!fList){
     AliError("No TList pointer ! ");
     return kFALSE;
   }
@@ -225,94 +242,242 @@ Bool_t AliHFEcollection::CreateTH1Fvector2(Int_t X, Int_t Y, const char* name, c
       }
     }
   }
-  return kTRUE;
+  return kTRUE;  
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::CreateProfile(const char* name, const char* title, Int_t nbins, Double_t xmin, Double_t xmax){
+  
+  //
+  // create a simple TProfile
+  //
+
+  if(!fList){
+    AliError("No TList pointer ! ");
+    return kFALSE;
+  }
+  fList->Add(new TProfile(name, title, nbins, xmin, xmax));
+  return CheckObject(name);
+
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::CreateTHnSparse(const char* name, const char* title, Int_t dim, Int_t* nbins, Double_t* xmin, Double_t* xmax){
+
+  //
+  // create 'dim' dimensional THnSparse
+  //
+
+  if(!fList){
+    AliError("No TList pointer ! ");
+    return kFALSE;
+  }
+  fList->Add(new THnSparseF(name, title, dim, nbins, xmin, xmax));
+  return CheckObject(name);
+
+}
+//___________________________________________________________________
+TObject* AliHFEcollection::Get(const char* name){ 
+
+  //
+  // Get histogram with the required name
+  // 
   
+
+  if(!CheckObject(name)){
+    AliError(Form("Not possible to return pointer to the object '%s'\n", name));
+    return 0;
+  }
+
+  return fList->FindObject(name);
   
 }
 //___________________________________________________________________
-TObject* AliHFEcollection::Get(const char* name, Int_t X){
+Bool_t AliHFEcollection::Fill(const char* name, Double_t v){
+
   //
-  // Get the histogram from the vector
-  // Paramter: 
-  //  name - vector name
-  //  X - Number of the desired histogram
+  // fill function for one TH1 histograms
   //
-  TString hname = name;
-  hname.Append(Form("_[%d]", X));
-  if(!CheckObject(hname.Data())){
-    AliError("No such object found in the list");
-    return 0x0;
+
+   if(!CheckObject(name)){
+    AliError(Form("Not possible to return pointer to the object '%s'\n", name));
+    return kFALSE;
   }
-  else{
-    return Get(hname.Data());
+
+  // chack the possible object types
+  if(fList->FindObject(name)->InheritsFrom("TH1")){
+    (dynamic_cast<TH1F*>(fList->FindObject(name)))->Fill(v);
+    return kTRUE;
+  }
+  
+  return kFALSE;
+
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::Fill(const char* name, Int_t X, Double_t v){
+
+  //
+  // fill function for one dimension arrays of TH1
+  //
+
+  const char* n = Form("%s_[%d]", name, X);
+  TObject *o = Get(n);
+  if(!o){
+    return kFALSE;
   }
+  Fill(o->GetName(), v);
+  return kTRUE;
 }
 //___________________________________________________________________
-TObject* AliHFEcollection::Get(const char* name, Int_t X, Int_t Y){
+Bool_t AliHFEcollection::Fill(const char* name, Int_t X, Int_t Y, Double_t v){
+  
+  const char* n = Form("%s_[%d][%d]", name, X, Y);
+  TObject *o = Get(n);
+  if(!o){
+    return kFALSE;
+  }
+  Fill(o->GetName(), v);
+  return kTRUE;
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::Fill(const char* name, Int_t X, Double_t v1, Double_t v2){
+
   //
-  // Get histogram from the 2D vector
-  // Parameters:
-  //  name - Name of the vector
-  //  X,Y - Indices of the histogram
+  // fill function for one dimension array of TH2
   //
-  TString hname = name;
-  hname.Append(Form("_[%d][%d]", X, Y));
-  if(!CheckObject(hname.Data())){
-    AliError("No such object found in the list");
-    AliError(Form("name: %s", hname.Data()));
-    return 0x0;
+
+  const char* n = Form("%s_[%d]", name, X);
+  TObject *o = Get(n);
+  if(!o){
+    return kFALSE;
   }
-  else{
-    return Get(hname.Data());
+  Fill(o->GetName(), v1, v2);
+  
+  return kTRUE;
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::Fill(const char* name, Double_t v1, Double_t v2){
+
+  //
+  // fill function for TH2 objects
+  //
+
+   if(!CheckObject(name)){
+    AliError(Form("Not possible to return pointer to the object '%s'\n", name));
+    return kFALSE;
   }
+
+  // chack the possible object types
+  if(fList->FindObject(name)->InheritsFrom("TH2")){
+    (dynamic_cast<TH2F*>(fList->FindObject(name)))->Fill(v1, v2);
+    return kTRUE;
+  }  
+  if(fList->FindObject(name)->InheritsFrom("TProfile")){
+    (dynamic_cast<TProfile*>(fList->FindObject(name)))->Fill(v1, v2);
+    return kTRUE;
+  }  
+  
+  return kFALSE;
+  
 }
+
 //___________________________________________________________________
 Bool_t AliHFEcollection::CheckObject(const char* name){
+
   //
   // check wheter the creation of the histogram was succesfull
   //
   
-  if(!fListE){
+  if(!fList){
     AliError("No TList pointer ! ");
     return kFALSE;
   }
   
-  if(!fListE->FindObject(name)){
-    AliError("Creating or Finding the object failed");
+  if(!fList->FindObject(name)){
+    AliError(Form("Creating or Finding the object '%s' failed\n", name));
     return kFALSE;
   }
   return kTRUE;
 }
 //___________________________________________________________________
-TObject* AliHFEcollection::Get(const char* name){ 
-  //
-  // Get histogram with the required name
+Bool_t AliHFEcollection::BinLogAxis(const char* name, Int_t dim){
+
   // 
-  return fListE->FindObject(name); 
+  // converts the axis (defined by the dimension) of THx or THnSparse
+  // object to Log scale. Number of bins and bin min and bin max are preserved
+  //
+
+
+  if(!CheckObject(name)){
+    return kFALSE;
+  }
+
+  TObject *o = Get(name);
+  TAxis *axis = 0x0;
+  if(o->InheritsFrom("TH1")){
+    axis = (dynamic_cast<TH1F*>(o))->GetXaxis();
+  }
+  if(o->InheritsFrom("TH2")){
+    if(0 == dim){
+      axis = (dynamic_cast<TH2F*>(o))->GetXaxis();
+    }
+    else if(1 == dim){
+      axis = (dynamic_cast<TH2F*>(o))->GetYaxis();
+    }
+     else{
+       AliError("Only dim = 0 or 1 possible for TH2F");
+     }
+  }
+  if(o->InheritsFrom("THnSparse")){
+    axis = (dynamic_cast<THnSparse*>(o))->GetAxis(dim);
+  }
+  
+  if(!axis){
+    AliError(Form("Axis '%d' could not be identified in the object '%s'\n", dim, name));
+    return kFALSE;
+  }
+  
+  Int_t bins = axis->GetNbins();
+
+  Double_t from = axis->GetXmin();
+  Double_t to = axis->GetXmax();
+  Double_t *newBins = new Double_t[bins+1];
+  newBins[0] = from;
+  Double_t factor = TMath::Power(to/from, 1./bins);
+  for(Int_t i=1; i<=bins; ++i){
+    newBins[i] = factor * newBins[i-1];
+  }
+  axis->Set(bins, newBins);
+  delete newBins;
+
+  return kTRUE;
+
+
 }
 //___________________________________________________________________
 Long64_t AliHFEcollection::Merge(TCollection *list){
+
   //
   // Merge the collections
   //
-  if(!fListE){
+
+  if(!fList){
     AliError("AliHFEcollection::Merge : No TList pointer ! ");
     return 0;
   }
 
-  return fListE->Merge(list);
+  return fList->Merge(list);
   
 }
 //____________________________________________________________________
 void AliHFEcollection::Browse(TBrowser *b)
 {
+
   //
   // Browse the content of the directory.
   //
 
    if (b) {
       TObject *obj = 0;
-      TIter nextin(fListE);
+      TIter nextin(fList);
 
       //Add objects that are only in memory
       while ((obj = nextin())) {
index 7666686..bd9424e 100644 (file)
@@ -52,25 +52,32 @@ class AliHFEcollection : public TNamed{
   Bool_t CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY);
 
   Bool_t CreateTH1Fvector1(Int_t X, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax);
+  Bool_t CreateTH1Fvector2(Int_t X, Int_t Y, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax);
   Bool_t CreateTH2Fvector1(Int_t X, const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY);
+  Bool_t CreateProfile(const char* name, const char* title, Int_t nbins, Double_t xmin, Double_t xmax);
+  Bool_t CreateTHnSparse(const char* name, const char* title, Int_t dim, Int_t* nbins, Double_t* xmin, Double_t* xmax);
 
-  Bool_t CreateTH1Fvector2(Int_t X, Int_t Y, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax);
-  
+  Bool_t BinLogAxis(const char* name, Int_t dim);
+    
 
   Long64_t Merge(TCollection *list);
 
   // Get functions
-  TList* GetList()  const  { return fListE; }
+  TList* GetList()  const  { return fList; }
   TObject* Get(const char* name); 
-  TObject* Get(const char* name, Int_t X);
-  TObject* Get(const char* name, Int_t X, Int_t Y);
 
+  // Fill functions
+  Bool_t Fill(const char* name, Double_t v);
+  Bool_t Fill(const char* name, Int_t X, Double_t v);
+  Bool_t Fill(const char* name, Int_t X, Int_t Y, Double_t v);
+  Bool_t Fill(const char* name, Double_t v1, Double_t v2);
+  Bool_t Fill(const char* name, Int_t X, Double_t v1, Double_t v2);
  private:
   Bool_t CheckObject(const char* name);
    void Copy(TObject &ref) const;
 
  private:
-  TList*                           fListE;      //! Object container
+  TList*                           fList;      //! Object container
 
   ClassDef(AliHFEcollection, 1)
 
index c420738..9f411b1 100644 (file)
@@ -127,11 +127,21 @@ void AliHFEcuts::Initialize(AliCFManager *cfm){
   SetHFElectronITSCuts();
   SetHFElectronTRDCuts();
 
+  // Publish to the cuts which analysis type they are (ESD Analysis by default)
+  if(IsAOD()){
+    AliInfo("Setting AOD Analysis");
+    TObjArray *genCuts = dynamic_cast<TObjArray *>(fCutList->FindObject("fPartGenCuts"));
+    if(genCuts){
+      AliCFParticleGenCuts *myGenCut = dynamic_cast<AliCFParticleGenCuts *>(genCuts->FindObject("fCutsGenMC"));
+      if(myGenCut) myGenCut->SetAODMC(kTRUE);
+    }
+  }
+
   // Connect the event cuts
-  SetEventCutList(kEventStepGenerated);
+  /*SetEventCutList(kEventStepGenerated);
   SetEventCutList(kEventStepReconstructed);
   cfm->SetEventCutsList(kEventStepGenerated, dynamic_cast<TObjArray *>(fCutList->FindObject("fEvGenCuts")));
-  cfm->SetEventCutsList(kEventStepReconstructed, dynamic_cast<TObjArray *>(fCutList->FindObject("fEvRecCuts")));
+  cfm->SetEventCutsList(kEventStepReconstructed, dynamic_cast<TObjArray *>(fCutList->FindObject("fEvRecCuts")));*/
   
   // Connect the particle cuts
   cfm->SetParticleCutsList(kStepMCGenerated, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartGenCuts")));
index 1b558bb..caf8ad7 100644 (file)
@@ -68,8 +68,12 @@ class AliHFEcuts : public TObject{
     TList *GetQAhistograms() const { return fHistQA; }
     
     void SetDebugMode();
+    void SetAOD() { SetBit(kAOD, kTRUE); }
+    void SetESD() { SetBit(kAOD, kFALSE); }
     void UnsetDebugMode() { SetBit(kDebugMode, kFALSE); }
     Bool_t IsInDebugMode() const { return TestBit(kDebugMode); };
+    Bool_t IsAOD() const { return TestBit(kAOD); }
+    Bool_t IsESD() const { return !TestBit(kAOD); }
     
     // Getters
     Bool_t IsRequireITSpixel() const { return TESTBIT(fRequirements, kITSPixel); };
@@ -105,7 +109,8 @@ class AliHFEcuts : public TObject{
 
   private:
     enum{
-      kDebugMode = BIT(14)
+      kDebugMode = BIT(14),
+      kAOD = BIT(15)
     };
     typedef enum{
       kPrimary = 0,
index 79e0bbc..ea6bcaa 100644 (file)
@@ -34,6 +34,7 @@
 
 #include <AliLog.h>
 #include <AliStack.h>
+#include <AliGenEventHeader.h>
 #include <AliAODMCParticle.h>
 
 #include "AliHFEmcQA.h"
@@ -48,6 +49,7 @@ ClassImp(AliHFEmcQA)
 //_______________________________________________________________________________________________
 AliHFEmcQA::AliHFEmcQA() : 
         fStack(0x0) 
+        ,fMCHeader(0x0)
         ,fMCArray(0x0)
         ,fNparents(0) 
 {
@@ -59,6 +61,7 @@ AliHFEmcQA::AliHFEmcQA() :
 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
         TObject(p)
         ,fStack(0x0) 
+        ,fMCHeader(0x0)
         ,fMCArray(0x0)
         ,fNparents(p.fNparents) 
 {
@@ -470,10 +473,17 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
     Int_t maPdgcode = partMother->GetPdgCode();
     Int_t maPdgcodeCopy = maPdgcode;
 
+    // get mc primary vertex
+    TArrayF mcPrimVtx;
+    if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
+
     // get electron production vertex   
     TLorentzVector ePoint;
     mcpart->ProductionVertex(ePoint);
 
+    // calculated production vertex to primary vertex (in xy plane)
+    Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
+
     // if the mother is charmed hadron  
     Bool_t isMotherDirectCharm = kFALSE;
     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { 
@@ -521,11 +531,8 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
                   // ratio between pT of electron and pT of mother B hadron 
                   if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
 
-                  // distance between electron production point and mother hadron production point
-                  TLorentzVector debPoint;
-                  grandMa->ProductionVertex(debPoint);
-                  TLorentzVector dedistance = ePoint - debPoint;
-                  fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),dedistance.P());
+                  // distance between electron production point and primary vertex
+                  fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
                   return;
                 }
              } 
@@ -552,11 +559,8 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
               // ratio between pT of electron and pT of mother B hadron 
               if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
 
-              // distance between electron production point and mother hadron production point
-              TLorentzVector ebPoint;
-              partMotherCopy->ProductionVertex(ebPoint);
-              TLorentzVector edistance = ePoint - ebPoint;
-              fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),edistance.P());
+              // distance between electron production point and primary vertex
+              fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
             }
          }
     } // end of if
@@ -777,11 +781,11 @@ Float_t AliHFEmcQA::GetRapidity(AliAODMCParticle *part) const
 }
 
 //__________________________________________
-Int_t AliHFEmcQA::GetElectronSource(AliAODMCParticle *mcpart)
+Int_t AliHFEmcQA::GetSource(AliAODMCParticle * const mcpart)
 {        
-  // decay electron's origin 
+  // decay particle's origin 
 
-  if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
+  //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
        
   Int_t origin = -1;
   Bool_t isFinalOpenCharm = kFALSE;
@@ -854,11 +858,11 @@ Int_t AliHFEmcQA::GetElectronSource(AliAODMCParticle *mcpart)
 }
 
 //__________________________________________
-Int_t AliHFEmcQA::GetElectronSource(TParticle* mcpart)
+Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
 {
-  // decay electron's origin 
+  // decay particle's origin 
 
-  if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
+  //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
 
   Int_t origin = -1;
   Bool_t isFinalOpenCharm = kFALSE;
index 19de700..767301a 100644 (file)
@@ -35,6 +35,7 @@ class TH2F;
 class TParticle;
 class TString;
 class AliStack;
+class AliGenEventHeader;
 class AliAODMCParticle;
 
 //________________________________________________________________
@@ -58,6 +59,7 @@ class AliHFEmcQA: public TObject {
                 void PostAnalyze() const;
                 void CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt=""); // create histograms for mc qa analysis
                 void SetStack(AliStack* const stack){fStack=stack;} // set stack pointer
+                void SetGenEventHeader(AliGenEventHeader* const mcHeader){fMCHeader=mcHeader;} // set stack pointer
                 void SetMCArray(TClonesArray* const mcarry){fMCArray=mcarry;} // set mcarray pointer
                 void Init();
 
@@ -66,8 +68,8 @@ class AliHFEmcQA: public TObject {
                 void GetDecayedKine(TParticle *part, const Int_t kquark, const Int_t kdecayed, Int_t icut); // get decay electron kinematics distribution
                void GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut); // get decay electron kinematics for AOD 
                 void EndOfEventAna(const Int_t kquark); // run analysis which should be done at the end of the event loop
-               Int_t GetElectronSource(TParticle* mcpart); // return electron source id 
-               Int_t GetElectronSource(AliAODMCParticle *mcpart); // return electron source id for AOD
+               Int_t GetSource(TParticle * const mcpart); // return electron source id 
+               Int_t GetSource(AliAODMCParticle * const mcpart); // return electron source id for AOD
 
         protected:
                 void IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel); // 
@@ -79,7 +81,8 @@ class AliHFEmcQA: public TObject {
                Float_t GetRapidity(AliAODMCParticle *part) const; // return rapidity
 
                 AliStack* fStack; // stack pointer           
-               TClonesArray *fMCArray; // mc array pointer
+                AliGenEventHeader* fMCHeader; // mcheader pointer
+                TClonesArray *fMCArray; // mc array pointer
 
                 static const Int_t fgkGluon; // gluon pdg code
                 static const Int_t fgkMaxGener; // ancester level wanted to be checked 
index 957877b..f6281c2 100644 (file)
@@ -15,6 +15,9 @@
 //
 // PID Steering Class 
 // Interface to the user task
+// Several strategies for Electron identification implemented.
+// In addition users can combine different detectors or use
+// single detector PID
 // 
 // Authors:
 //   Markus Fasel <M.Fasel@gsi.de>
@@ -149,6 +152,7 @@ Bool_t AliHFEpid::InitializePID(TString arg){
       case 3: InitStrategy3(); break;
       case 4: InitStrategy4(); break;
       case 5: InitStrategy5(); break;
+      case 6: InitStrategy6(); break;
       default: strategyStatus = kFALSE;
     }
     return strategyStatus;
@@ -170,6 +174,7 @@ Bool_t AliHFEpid::InitializePID(TString arg){
       fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID");
       SETBIT(fEnabledDetectors, kTRDpid);
     } else if(det->String().CompareTo("TOF") == 0){
+      AliInfo("Doing TOF PID");
       fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOF PID");
       SETBIT(fEnabledDetectors, kTOFpid);
     }
@@ -197,7 +202,7 @@ Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track){
     // MC Event
     return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11);
   }
-  if(fPIDstrategy > 0 && fPIDstrategy < 6){
+  if(fPIDstrategy > 0 && fPIDstrategy < 7){
     Int_t pid = 0;
     switch(fPIDstrategy){
       case 1: pid = IdentifyStrategy1(track); break;
@@ -205,6 +210,7 @@ Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track){
       case 3: pid = IdentifyStrategy3(track); break;
       case 4: pid = IdentifyStrategy4(track); break;
       case 5: pid = IdentifyStrategy5(track); break;
+      case 6: pid = IdentifyStrategy6(track); break;
       default: break;
     }
     return pid;
@@ -236,8 +242,72 @@ Bool_t AliHFEpid::MakePidTpcTof(AliHFEpidObject *track){
   //
   // Combines TPC and TOF PID decision
   //
-  if(fDetectorPID[kTOFpid]->IsSelected(track)) return fDetectorPID[kTPCpid]->IsSelected(track);
-  return kFALSE;
+  if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return kFALSE;
+  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
+  if(!esdTrack) return kFALSE;
+
+  AliHFEpidTOF *tofPID = dynamic_cast<AliHFEpidTOF*>(fDetectorPID[kTOFpid]);
+  if(!tofPID){
+    AliWarning("TOF pid object is NULL");
+    return kFALSE;
+  }
+
+  
+  AliHFEpidTPC *tpcPID = dynamic_cast<AliHFEpidTPC*>(fDetectorPID[kTPCpid]);
+  if(!tpcPID){
+    AliWarning("TPC pid object is NULL");
+    return kFALSE;
+  }
+  
+  // charge of the particle
+  Int_t charge = 0;
+  if(esdTrack->GetOuterParam()) charge = esdTrack->GetOuterParam()->Charge();
+  else charge = esdTrack->GetInnerParam()->Charge();
+  
+  // setup the TPC PID
+  const Int_t cTPCsigma = 2;
+  // set the number of sigmas in the TPC 
+  tpcPID->SetTPCnSigma(cTPCsigma);
+  // turn on the the dEdx line crossing for kaons and protons
+  tpcPID->AddTPCdEdxLineCrossing(3, cTPCsigma);  // for kaons
+  tpcPID->AddTPCdEdxLineCrossing(4, cTPCsigma);  // for protons
+  // request the tpc PID information
+  Int_t pidTPC = tpcPID->IsSelected(track);
+
+  // without TPC pid it makes no sense to look at the TOF pid with this detector combination
+  //if(0 == pidTPC){
+  //  return kFALSE;
+  //}
+
+  // is the TPC in the crossing line region? 0 - no crossing, otherwise the AliPID of the crossing
+  // particle returned - 3 for Kaon and 4 for Proton
+  Int_t crossing = tpcPID->GetCrossingType();
+  // setup the TOF pid
+  const Int_t cTOFsigma = 3;
+  tofPID->SetTOFnSigma(cTOFsigma);
+  // request TOF PID information
+  Int_t pidTOF = tofPID->IsSelected(track);
+  // case that the TOF for some reason does not deliver a PID information 
+  // or the TPC is not in the crossing point region, only the TPC will be used
+  if(0 != pidTOF && 0 != crossing){
+    if(3 == crossing){
+      return (321 == pidTOF) ? kFALSE : kTRUE;
+    }
+    else if(4 == crossing){
+      return (2212 == pidTOF) ? kFALSE : kTRUE;
+    }    
+    else{
+      // something went wrong
+      AliError("Wrong crossing type returned from AliHFEpidTPC - check your code!!");
+      return kFALSE;
+    }
+  }
+  else{
+    // tpc ONLY
+    return (11 == pidTPC) ? kTRUE : kFALSE;   
+  }
 }
 
 //____________________________________________________________
@@ -397,7 +467,7 @@ void AliHFEpid::InitStrategy3(){
   //   
   AliHFEpidTPC *pid = new AliHFEpidTPC("strat3TPCpid");
   pid->SetTPCnSigma(3);
-  pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 3.);
+  pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
   Bool_t status = pid->InitializePID();
   if(IsQAOn() && status) pid->SetQAOn(fQAlist);
   if(HasMCData() && status) pid->SetHasMCData();
@@ -431,6 +501,32 @@ void AliHFEpid::InitStrategy5(){
 }
 
 //____________________________________________________________
+void AliHFEpid::InitStrategy6(){
+  //
+  // Combined TPC-TOF PID, combination is discribed in the funtion MakePidTpcTof
+  //
+  AliHFEpidTPC *tpcpid = new AliHFEpidTPC("strat6TPCpid");
+  AliHFEpidTOF *tofpid = new AliHFEpidTOF("strat6TOFpid");
+  Bool_t status = tpcpid->InitializePID();
+  if(!status)
+    AliError("Initialization of TPC PID failed");
+  Bool_t status1 = tofpid->InitializePID();
+  if(!status1)
+    AliError("Initialization of TOF PID failed");
+  status &= status1;
+  if(IsQAOn() && status){
+    tpcpid->SetQAOn(fQAlist);
+    tofpid->SetQAOn(fQAlist);
+  }
+  if(HasMCData() && status){
+    tpcpid->SetHasMCData();
+    tofpid->SetHasMCData();
+  }
+  fDetectorPID[kTPCpid] = tpcpid;
+  fDetectorPID[kTOFpid] = tofpid;
+}
+
+//____________________________________________________________
 Bool_t AliHFEpid::IdentifyStrategy1(AliHFEpidObject *track){
   return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
 }
@@ -455,3 +551,9 @@ Bool_t AliHFEpid::IdentifyStrategy4(AliHFEpidObject *track){
 Bool_t AliHFEpid::IdentifyStrategy5(AliHFEpidObject *track){
   return MakePidTpcTrd(track);
 }
+
+//____________________________________________________________
+Bool_t AliHFEpid::IdentifyStrategy6(AliHFEpidObject *track){
+  return MakePidTpcTof(track);
+}
+
index 8f904b6..3f40c12 100644 (file)
@@ -62,11 +62,13 @@ class AliHFEpid : public TObject{
     void InitStrategy3();
     void InitStrategy4();
     void InitStrategy5();
+    void InitStrategy6();
     Bool_t IdentifyStrategy1(AliHFEpidObject *track);
     Bool_t IdentifyStrategy2(AliHFEpidObject *track);
     Bool_t IdentifyStrategy3(AliHFEpidObject *track);
     Bool_t IdentifyStrategy4(AliHFEpidObject *track);
     Bool_t IdentifyStrategy5(AliHFEpidObject *track);
+    Bool_t IdentifyStrategy6(AliHFEpidObject *track);
   private:
     enum{
       kIsQAOn = BIT(14),
index 2323840..b3619d6 100644 (file)
@@ -32,6 +32,7 @@
 #include "AliLog.h"
 #include "AliMCParticle.h"
 #include "AliPID.h"
+#include "AliTOFpidESD.h"
 
 #include "AliHFEpidTOF.h"
 #include "AliHFEpidBase.h"
@@ -44,16 +45,23 @@ AliHFEpidTOF::AliHFEpidTOF(const Char_t *name):
   AliHFEpidBase(name)
   , fPID(0x0)
   , fQAList(0x0)
+  , fPIDtofESD(0x0)
+  , fNsigmaTOF(3)
 {
   //
   // Constructor
   //
+
+  fPIDtofESD = new AliTOFpidESD();
+
 }
 //___________________________________________________________________
 AliHFEpidTOF::AliHFEpidTOF(const AliHFEpidTOF &c):
   AliHFEpidBase("")
   , fPID(0x0)
   , fQAList(0x0)
+  , fPIDtofESD(0x0)
+  , fNsigmaTOF(3)
 {  
   // 
   // Copy operator
@@ -79,6 +87,7 @@ AliHFEpidTOF::~AliHFEpidTOF(){
   // Destructor
   //
   if(fPID) delete fPID;
+  if(fPIDtofESD) delete fPIDtofESD;
   if(fQAList){
     fQAList->Delete();
     delete fQAList;
@@ -93,6 +102,7 @@ void AliHFEpidTOF::Copy(TObject &ref) const {
 
   target.fPID = fPID;          
   target.fQAList = fQAList;
+  target.fPIDtofESD = new AliTOFpidESD(*fPIDtofESD);
 
   AliHFEpidBase::Copy(ref);
 }
@@ -114,8 +124,9 @@ Int_t AliHFEpidTOF::IsSelected(AliHFEpidObject *vtrack)
   // the ESD PID will be checked and if necessary improved 
   // in the mean time. Best of luck
   //
-  // returns 10 (== kUnknown)if PID can not be assigned
+  // returns 10 (== kUnknown)if PID can not be assigned for whatever reason
   //
+
   if(vtrack->fAnalysisType == AliHFEpidObject::kESDanalysis){
     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(vtrack->fRecTrack);
     if(!esdTrack) return 0;
@@ -137,25 +148,29 @@ Int_t AliHFEpidTOF::MakePIDesd(AliESDtrack *track, AliMCParticle * /*mcTrack*/){
   Long_t status = 0;
   status = track->GetStatus(); 
 
-  if(!(status & AliESDtrack::kTOFout)) return AliPID::kUnknown;
+  if(!(status & AliESDtrack::kTOFout)) return 0;
   
-  (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(0.);
+  if(IsQAon())(dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(0.);
 
   Double_t tItrackL = track->GetIntegratedLength();
   Double_t tTOFsignal = track->GetTOFsignal();
   
-  if(tItrackL > 0)
-    (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(1.);
+  if(IsQAon()){
+    if(tItrackL > 0)
+      (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(1.);
 
-  if(tTOFsignal > 0)
-    (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(2.);
+    if(tTOFsignal > 0)
+      (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(2.);
+  }
   
 
-  if(tItrackL <=0 || tTOFsignal <=0) return AliPID::kUnknown;
+  if(tItrackL <=0 || tTOFsignal <=0) return 0;
 
-  (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(3.);
-  (dynamic_cast<TH1F *>(fQAList->At(kHistTOFsignal)))->Fill(tTOFsignal/1000.);
-  (dynamic_cast<TH1F *>(fQAList->At(kHistTOFlength)))->Fill(tItrackL);
+  if(IsQAon()){
+    (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(3.);
+    (dynamic_cast<TH1F *>(fQAList->At(kHistTOFsignal)))->Fill(tTOFsignal/1000.);
+    (dynamic_cast<TH1F *>(fQAList->At(kHistTOFlength)))->Fill(tItrackL);
+  }
   // get the TOF pid probabilities
   Double_t tESDpid[5] = {0., 0., 0., 0., 0.};
   Float_t tTOFpidSum = 0.;
@@ -170,20 +185,68 @@ Int_t AliHFEpidTOF::MakePIDesd(AliESDtrack *track, AliMCParticle * /*mcTrack*/){
       tMAXindex = i;
     }
   }
+
+  Int_t pdg = 0;
+
+  switch(tMAXindex){
+    case 0:    pdg = 11; break;
+    case 1:    pdg = 13; break;
+    case 2:    pdg = 211; break;
+    case 3:    pdg = 321; break;
+    case 4:    pdg = 2212; break;
+    default:   pdg = 0;
+  };
+
   
   Double_t p = track->GetOuterParam()->P();
   Double_t beta = (tItrackL/100.)/(TMath::C()*(tTOFsignal/1e12));
   
-  if(TMath::Abs(tTOFpidSum - 1) > 0.01) return AliPID::kUnknown;
-  else{
-    // should be the same as AliPID flags
-    
+  // sanity check, should not be necessary
+  if(TMath::Abs(tTOFpidSum - 1) > 0.01) return 0;
+
+  Double_t nSigmas = fPIDtofESD->GetNumberOfSigmas(track, (AliPID::EParticleType)tMAXindex);
+  if(TMath::Abs(nSigmas) > fNsigmaTOF) return 0;
+
+  
+  // should be the same as AliPID flags
+  
+  if(IsQAon()){
     (dynamic_cast<TH2F *>(fQAList->At(kHistTOFpid0+tMAXindex)))->Fill(beta, p);
     (dynamic_cast<TH2F *>(fQAList->At(kHistTOFpidBetavP)))->Fill(beta, p);
-    return tMAXindex;
   }
+  //return tMAXindex;
+  return pdg;
+  
+}
+//___________________________________________________________________
+Double_t AliHFEpidTOF::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig){
+  
+  //gives probability for track to come from a certain particle species;
+  //no priors considered -> probabilities for equal abundances of all species!
+  //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
+  
+  //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
+  
+  if(!track) return -1.;
+  Bool_t outlier = kTRUE;
+  // Check whether distance from the respective particle line is smaller than r sigma
+  for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
+    if(TMath::Abs(fPIDtofESD->GetNumberOfSigmas(track, (AliPID::EParticleType)hypo)) > rsig)
+      outlier = kTRUE;
+    else {
+      outlier = kFALSE;
+      break;
+    }
+  }
+  if(outlier)
+    return -2.;
+  
+  Double_t tofProb[5];
+  
+  track->GetTOFpid(tofProb);
+  
+  return tofProb[species];
 }
-
 //___________________________________________________________________
 Int_t AliHFEpidTOF::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mcTrack*/){
   AliError("AOD PID not yet implemented");
index 641fb4c..9ae34e4 100644 (file)
@@ -20,6 +20,7 @@ class AliAODTrack;
 class AliAODMCParticle;
 class AliESDtrack;
 class AliMCParticle;
+class AliTOFpidESD;
 
 class AliHFEpidTOF : public AliHFEpidBase{
   public:
@@ -32,7 +33,10 @@ class AliHFEpidTOF : public AliHFEpidBase{
     virtual Int_t     IsSelected(AliHFEpidObject *track);
     virtual Bool_t    HasQAhistos() const { return kTRUE; };
   
-  
+    void SetTOFnSigma(Short_t nSigma) { fNsigmaTOF = nSigma; };
+
+    Double_t Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig = 2.);  
   protected:
     void Copy(TObject &ref) const;
     void AddQAhistograms(TList *qaHist);
@@ -50,12 +54,15 @@ class AliHFEpidTOF : public AliHFEpidBase{
       kHistTOFpid2 = 6,
       kHistTOFpid3 = 7,
       kHistTOFpid4 = 8
-      
     } QAHist_t;
+  
+    AliPID        *fPID;           //! PID Object
+    TList         *fQAList;        //! QA histograms
+    AliTOFpidESD  *fPIDtofESD;     //! TOF pid object
+
+    Short_t fNsigmaTOF;            // TOF sigma band
 
-    AliPID *fPID;           //! PID Object
-    TList *fQAList;         //! QA histograms
-    ClassDef(AliHFEpidTOF, 1) // TOF PID class
+    ClassDef(AliHFEpidTOF, 1)
 };
 
 #endif
index 99b2e57..893e090 100644 (file)
@@ -47,7 +47,8 @@
 //___________________________________________________________________
 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
   // add a list here
-    AliHFEpidBase(name)
+  AliHFEpidBase(name)
+  , fLineCrossingType(0)
   , fLineCrossingsEnabled(0)
   , fNsigmaTPC(3)
   , fRejectionEnabled(0)
@@ -67,7 +68,8 @@ AliHFEpidTPC::AliHFEpidTPC(const char* name) :
 
 //___________________________________________________________________
 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
-    AliHFEpidBase("")
+  AliHFEpidBase("")
+  , fLineCrossingType(0)
   , fLineCrossingsEnabled(0)
   , fNsigmaTPC(2)
   , fRejectionEnabled(0)
@@ -156,22 +158,24 @@ Int_t AliHFEpidTPC::IsSelected(AliHFEpidObject *track)
     return MakePIDaod(aodtrack, aodmctrack);
   }
 }
-
 //___________________________________________________________________
 Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
   //
   //  Doing TPC PID as explained in IsSelected for ESD tracks
   //
   if(IsQAon()) FillTPChistograms(esdTrack, mctrack);
+  Float_t nsigma = fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron);
   // exclude crossing points:
   // Determine the bethe values for each particle species
   Bool_t isLineCrossing = kFALSE;
+  fLineCrossingType = 0;  // default value
   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
     if(ispecies == AliPID::kElectron) continue;
     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
-    if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies) < fLineCrossingSigma[ispecies]){
-      // Point in a line crossing region, no PID possible
-      isLineCrossing = kTRUE;
+    if(TMath::Abs(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
+      // Point in a line crossing region, no PID possible, but !PID still possible ;-)
+      isLineCrossing = kTRUE;      
+      fLineCrossingType = ispecies;
       break;
     }
   }
@@ -180,13 +184,11 @@ Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
   // Check particle rejection
   if(HasParticleRejection()){
     Int_t reject = Reject(esdTrack);
-    AliDebug(1, Form("PID code from Rejection: %d", reject));
     if(reject != 0) return reject;
   }
   // Check whether distance from the electron line is smaller than n-sigma
 
   // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
-  Float_t nsigma = fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron);
   Float_t p = 0.;
   Int_t pdg = 0;
   if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){ 
@@ -195,6 +197,7 @@ Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
     if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
   }
   if(IsQAon() && pdg != 0) (dynamic_cast<TH2I *>(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());
+
   return pdg;
 }
 
@@ -213,11 +216,9 @@ Int_t AliHFEpidTPC::Reject(AliESDtrack *track){
   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
-    AliDebug(1, Form("Particle Rejection enabled for species %d", ispec));
     // Particle rejection enabled
     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
     Double_t sigma = fPIDtpcESD->GetNumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec));
-    AliDebug(1, Form("Sigma %f, min %f, max %f", sigma, fRejection[4*ispec + 1], fRejection[4*ispec+3]));
     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
   }
   return 0;
@@ -351,14 +352,14 @@ void AliHFEpidTPC::AddQAhistograms(TList *qaList){
   fQAList = new TList;
   fQAList->SetName("fTPCqaHistos");
 
-  fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 200, 0, 200), kHistTPCelectron); 
-  fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 200, 0, 200), kHistTPCmuon);
-  fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 200, 0, 200), kHistTPCpion);
-  fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 200, 0, 200), kHistTPCkaon);
-  fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 200, 0, 200), kHistTPCproton);
-  fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 200, 0, 200), kHistTPCothers);
-  fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 200, 0, 200), kHistTPCall);
-  fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 200, 0, 200), kHistTPCselected);
+  fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600), kHistTPCelectron); 
+  fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600), kHistTPCmuon);
+  fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600), kHistTPCpion);
+  fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600), kHistTPCkaon);
+  fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton);
+  fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers);
+  fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall);
+  fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 60, 0, 600), kHistTPCselected);
 
   fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
   fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
index c8c958f..e97d281 100644 (file)
@@ -47,6 +47,8 @@ class AliHFEpidTPC : public AliHFEpidBase{
     virtual Int_t IsSelected(AliHFEpidObject *track);
     virtual Bool_t HasQAhistos() const { return kTRUE; };
 
+    Int_t GetCrossingType() const {return fLineCrossingType; }
+
     void AddTPCdEdxLineCrossing(Int_t species, Double_t sigma);
     Bool_t HasAsymmetricSigmaCut() const { return TestBit(kAsymmetricSigmaCut);}
     Bool_t HasParticleRejection() const { return TestBit(kRejection); }
@@ -99,6 +101,7 @@ class AliHFEpidTPC : public AliHFEpidBase{
       kRejection = BIT(21)
     };
     Double_t fLineCrossingSigma[AliPID::kSPECIES];          // with of the exclusion point
+    Int_t    fLineCrossingType;                             // 0 for no line crossing, otherwise AliPID of the particle crossing the electron dEdx band
     UChar_t fLineCrossingsEnabled;                          // Bitmap showing which line crossing is set
     Float_t fPAsigCut[2];                                   // Momentum region where to perform asymmetric sigma cut
     Float_t fNAsigmaTPC[2];                                 // Asymmetric TPC Sigma band        
index 13bc675..976ef8d 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 //
-//  QA class of primary vertex study for Heavy Flavor electrons
-//
+// QA class of primary vertex study for Heavy Flavor electrons
+// this has functionality to reject electrons from primary vertex
+// and check primary vertex characteristics
+//  
 // Authors:
 //   MinJung Kweon <minjung@physi.uni-heidelberg.de>
 //
index 5959643..2a01d23 100644 (file)
@@ -14,6 +14,8 @@
  **************************************************************************/
 //
 // QA class of primary vertex study for Heavy Flavor electrons
+// this has functionality to reject electrons from primary vertex
+// and check primary vertex characteristics
 //
 
 #ifndef ALIHFEPRIVTX_H
index 774cd7d..928e41f 100644 (file)
@@ -1,4 +1,4 @@
-/**************************************************************************
+ /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
 #include <TParticle.h>
 
 #include <AliESDEvent.h>
+#include <AliAODEvent.h>
+#include <AliVTrack.h>
 #include <AliESDtrack.h>
+#include <AliAODTrack.h>
 #include "AliHFEsecVtx.h"
 #include <AliKFParticle.h>
 #include <AliKFVertex.h>
 #include <AliLog.h>
 #include <AliStack.h>
+#include <AliAODMCParticle.h>
+#include "AliHFEpairs.h"
+#include "AliHFEsecVtxs.h"
 
 ClassImp(AliHFEsecVtx)
 
 //_______________________________________________________________________________________________
 AliHFEsecVtx::AliHFEsecVtx():
-         fESD1(0x0)
-        ,fStack(0x0)
-        ,fNparents(0)
-        ,fHistTagged()
-        ,fPairTagged(0)
-        ,fdistTwoSecVtx(-1)
-        ,fcosPhi(-1)
-        ,fsignedLxy(-1)
-        ,finvmass(-1)
-        ,finvmassSigma(-1)
-        ,fBTagged(0)
-        ,fBElec(0)
+  fESD1(0x0)
+  ,fAOD1(0x0)
+  ,fStack(0x0)
+  ,fUseMCPID(kFALSE)
+  ,fkSourceLabel()
+  ,fNparents(0)
+  ,fParentSelect()
+  ,fNoOfHFEpairs(0)
+  ,fNoOfHFEsecvtxs(0)
+  ,fHFEpairs(0x0)
+  ,fHFEsecvtxs(0x0)
+  ,fMCArray(0x0)
+  ,fPVx(0)
+  ,fPVy(0)
+  ,fCosPhi(-1)
+  ,fSignedLxy(-1)
+  ,fKFchi2(-1)
+  ,fInvmass(-1)
+  ,fInvmassSigma(-1)
+  ,fKFip(0)
+  ,fPairQA(0x0)
+  ,fSecvtxQA(0x0)
+  ,fSecVtxList(0x0)
 { 
-        // Default constructor
-
-        Init();
+  //
+  // Default constructor
+  //
 
+  Init();
 }
 
 //_______________________________________________________________________________________________
 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
-         TObject(p)
-        ,fESD1(0x0)
-        ,fStack(0x0)
-        ,fNparents(p.fNparents)
-        ,fHistTagged()
-        ,fPairTagged(p.fPairTagged)
-        ,fdistTwoSecVtx(p.fdistTwoSecVtx)
-        ,fcosPhi(p.fcosPhi)
-        ,fsignedLxy(p.fsignedLxy)
-        ,finvmass(p.finvmass)
-        ,finvmassSigma(p.finvmassSigma)
-        ,fBTagged(p.fBTagged)
-        ,fBElec(p.fBElec)
+  TObject(p)
+  ,fESD1(0x0)
+  ,fAOD1(0x0)
+  ,fStack(0x0)
+  ,fUseMCPID(p.fUseMCPID)
+  ,fkSourceLabel()
+  ,fNparents(p.fNparents)
+  ,fParentSelect()
+  ,fNoOfHFEpairs(p.fNoOfHFEpairs)
+  ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
+  ,fHFEpairs(0x0)
+  ,fHFEsecvtxs(0x0)
+  ,fMCArray(0x0)
+  ,fPVx(p.fPVx)
+  ,fPVy(p.fPVy)
+  ,fCosPhi(p.fCosPhi)
+  ,fSignedLxy(p.fSignedLxy)
+  ,fKFchi2(p.fKFchi2)
+  ,fInvmass(p.fInvmass)
+  ,fInvmassSigma(p.fInvmassSigma)
+  ,fKFip(p.fKFip)
+  ,fPairQA(0x0)
+  ,fSecvtxQA(0x0)
+  ,fSecVtxList(0x0)
 { 
-        // Copy constructor
+  //
+  // Copy constructor
+  //
 }
 
 //_______________________________________________________________________________________________
 AliHFEsecVtx&
 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
 {
+  //
   // Assignment operator
+  //
 
   AliInfo("Not yet implemented.");
   return *this;
@@ -87,1216 +120,987 @@ AliHFEsecVtx::operator=(const AliHFEsecVtx &)
 //_______________________________________________________________________________________________
 AliHFEsecVtx::~AliHFEsecVtx()
 {
-        // Destructor
+  //
+  // Destructor
+  //
 
-        //cout << "Analysis Done." << endl;
+  //cout << "Analysis Done." << endl;
 }
 
 //__________________________________________
 void AliHFEsecVtx::Init()
 {
-
-        // set pdg code and index
-
-        fNparents = 7;
-
-        fParentSelect[0][0] =  411;
-        fParentSelect[0][1] =  421;
-        fParentSelect[0][2] =  431;
-        fParentSelect[0][3] = 4122;
-        fParentSelect[0][4] = 4132;
-        fParentSelect[0][5] = 4232;
-        fParentSelect[0][6] = 4332;
-
-        fParentSelect[1][0] =  511;
-        fParentSelect[1][1] =  521;
-        fParentSelect[1][2] =  531;
-        fParentSelect[1][3] = 5122;
-        fParentSelect[1][4] = 5132;
-        fParentSelect[1][5] = 5232;
-        fParentSelect[1][6] = 5332;
-
-/*
-        fid[0][0] = 0;
-        fid[0][1] = 1;
-        fid[0][2] = 2;
-
-        fid[1][0] = 0;
-        fid[1][1] = 1;
-        fid[1][2] = 3;
-
-        fid[2][0] = 0;
-        fid[2][1] = 2;
-        fid[2][2] = 3;
-
-        fid[3][0] = 1;
-        fid[3][1] = 2;
-        fid[3][2] = 3;
-
-        fia[0][0][0] = 0;
-        fia[0][0][1] = 1;
-        fia[0][1][0] = 0;
-        fia[0][1][1] = 2;
-        fia[0][2][0] = 1;
-        fia[0][2][1] = 2;
-
-        fia[1][0][0] = 0;
-        fia[1][0][1] = 1;
-        fia[1][1][0] = 0;
-        fia[1][1][1] = 3;
-        fia[1][2][0] = 1;
-        fia[1][2][1] = 3;
-
-        fia[2][0][0] = 0;
-        fia[2][0][1] = 2;
-        fia[2][1][0] = 0;
-        fia[2][1][1] = 3;
-        fia[2][2][0] = 2;
-        fia[2][2][1] = 3;
-
-        fia[3][0][0] = 1;
-        fia[3][0][1] = 2;
-        fia[3][1][0] = 1;
-        fia[3][1][1] = 3;
-        fia[3][2][0] = 2;
-        fia[3][2][1] = 3;
-*/
-
+  //
+  // set pdg code and index
+  //
+
+  fNparents = 7;
+
+  fParentSelect[0][0] =  411;
+  fParentSelect[0][1] =  421;
+  fParentSelect[0][2] =  431;
+  fParentSelect[0][3] = 4122;
+  fParentSelect[0][4] = 4132;
+  fParentSelect[0][5] = 4232;
+  fParentSelect[0][6] = 4332;
+
+  fParentSelect[1][0] =  511;
+  fParentSelect[1][1] =  521;
+  fParentSelect[1][2] =  531;
+  fParentSelect[1][3] = 5122;
+  fParentSelect[1][4] = 5132;
+  fParentSelect[1][5] = 5232;
+  fParentSelect[1][6] = 5332;
 } 
 
-//__________________________________________
-void AliHFEsecVtx::ResetTagVar()
-{
-        // reset tag variables
-
-        fdistTwoSecVtx = -1;
-        fcosPhi = -1;
-        fsignedLxy = -1;
-        finvmass = -1;
-        finvmassSigma = -1;
-        fBTagged = kFALSE;
-        fBElec = kFALSE;
+//_______________________________________________________________________________________________
+void AliHFEsecVtx::CreateHistograms(TList *qaList)
+{ 
+  //
+  // create histograms
+  //
+
+  fSecVtxList = new TList;
+  fSecVtxList->SetName("SecVtx");
+
+  MakeContainer();
+  /*
+  fkSourceLabel[kAll]="all";
+  fkSourceLabel[kDirectCharm]="directCharm";
+  fkSourceLabel[kDirectBeauty]="directBeauty";
+  fkSourceLabel[kBeautyCharm]="beauty2charm";
+  fkSourceLabel[kGamma]="gamma";
+  fkSourceLabel[kPi0]="pi0";
+  fkSourceLabel[kElse]="others";
+  fkSourceLabel[kBeautyGamma]="beauty22gamma";
+  fkSourceLabel[kBeautyPi0]="beauty22pi0";
+  fkSourceLabel[kBeautyElse]="beauty22others";
+
+  TString hname;
+  TString hnopt="secvtx_";
+  for (Int_t isource = 0; isource < 10; isource++ ){
+  }
+  */
+
+  qaList->AddLast(fSecVtxList);
 }
 
-//__________________________________________
-void AliHFEsecVtx::InitAnaPair()
+//_______________________________________________________________________________________________
+void AliHFEsecVtx::GetPrimaryCondition()
 {
-        // initialize pair tagging variables
-
-        fPairTagged = 0;
-        for (Int_t i=0; i<20; i++){
-                fpairedTrackID[i] = -1;
-                fpairedChi2[i] = -1;
-                fpairedInvMass[i] = -1;
-                fpairedSignedLxy[i] = -1;
-        }
-
-        fid[0][0] = 0;
-        fid[0][1] = 1;
-        fid[0][2] = 2;
-
-        fid[1][0] = 0;
-        fid[1][1] = 1;
-        fid[1][2] = 3;
-
-        fid[2][0] = 0;
-        fid[2][1] = 2;
-        fid[2][2] = 3;
-
-        fid[3][0] = 1;
-        fid[3][1] = 2;
-        fid[3][2] = 3;
-
-        fia[0][0][0] = 0;
-        fia[0][0][1] = 1;
-        fia[0][1][0] = 0;
-        fia[0][1][1] = 2;
-        fia[0][2][0] = 1;
-        fia[0][2][1] = 2;
-
-        fia[1][0][0] = 0;
-        fia[1][0][1] = 1;
-        fia[1][1][0] = 0;
-        fia[1][1][1] = 3;
-        fia[1][2][0] = 1;
-        fia[1][2][1] = 3;
-
-        fia[2][0][0] = 0;
-        fia[2][0][1] = 2;
-        fia[2][1][0] = 0;
-        fia[2][1][1] = 3;
-        fia[2][2][0] = 2;
-        fia[2][2][1] = 3;
-
-        fia[3][0][0] = 1;
-        fia[3][0][1] = 2;
-        fia[3][1][0] = 1;
-        fia[3][1][1] = 3;
-        fia[3][2][0] = 2;
-        fia[3][2][1] = 3;
-
+  //
+  // get primary characteristics and set
+  //
+
+  if (fESD1) {
+    AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
+    if( primVtxCopy.GetNDF() <1 ) return;
+    fPVx = primVtxCopy.GetX();
+    fPVx = primVtxCopy.GetY();
+  }
+  else if(fAOD1) {
+    AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
+    if( primVtxCopy.GetNDF() <1 ) return;
+    fPVx = primVtxCopy.GetX();
+    fPVx = primVtxCopy.GetY();
+  }
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::CreateHistograms(TString hnopt)
-{ 
-        // create histograms
-
-        fkSourceLabel[kAll]="all";
-        fkSourceLabel[kDirectCharm]="directCharm";
-        fkSourceLabel[kDirectBeauty]="directBeauty";
-        fkSourceLabel[kBeautyCharm]="beauty2charm";
-        fkSourceLabel[kGamma]="gamma";
-        fkSourceLabel[kPi0]="pi0";
-        fkSourceLabel[kElse]="others";
-        fkSourceLabel[kBeautyGamma]="beauty22gamma";
-        fkSourceLabel[kBeautyPi0]="beauty22pi0";
-        fkSourceLabel[kBeautyElse]="beauty22others";
-
-
-        TString hname;
-        for (Int_t isource = 0; isource < 10; isource++ ){
-
-           hname=hnopt+"InvMass_"+fkSourceLabel[isource];
-           fHistPair[isource].fInvMass = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
-           hname=hnopt+"InvMassCut1_"+fkSourceLabel[isource];
-           fHistPair[isource].fInvMassCut1 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
-           hname=hnopt+"InvMassCut2_"+fkSourceLabel[isource];
-           fHistPair[isource].fInvMassCut2 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
-           hname=hnopt+"KFChi2_"+fkSourceLabel[isource];
-           fHistPair[isource].fKFChi2 = new TH1F(hname,hname,200,0,20);
-           hname=hnopt+"OpenAngle_"+fkSourceLabel[isource];
-           fHistPair[isource].fOpenAngle = new TH1F(hname,hname,100,0,3.14);
-           hname=hnopt+"CosOpenAngle_"+fkSourceLabel[isource];
-           fHistPair[isource].fCosOpenAngle = new TH1F(hname,hname,100,-1.1,1.1);
-           hname=hnopt+"SignedLxy_"+fkSourceLabel[isource];
-           fHistPair[isource].fSignedLxy = new TH2F(hname,hname,1000,-5,5,120,-2,10);
-           hname=hnopt+"KFIP_"+fkSourceLabel[isource];
-           fHistPair[isource].fKFIP = new TH1F(hname,hname,1000,-5,5);
-           hname=hnopt+"IPMax_"+fkSourceLabel[isource];
-           fHistPair[isource].fIPMax= new TH1F(hname,hname,500,0,5);
-
-        }
+void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
+{
+  //
+  // calculate e-h pair characteristics and tag pair 
+  //
 
-        hname=hnopt+"pt_beautyelec";
-        fHistTagged.fPtBeautyElec= new TH1F(hname,hname,250,0,50);
-        hname=hnopt+"pt_taggedelec";
-        fHistTagged.fPtTaggedElec= new TH1F(hname,hname,250,0,50);
-        hname=hnopt+"pt_righttaggedelec";
-        fHistTagged.fPtRightTaggedElec = new TH1F(hname,hname,250,0,50);
-        hname=hnopt+"pt_wrongtaggedelec";
-        fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,250,0,50);
-
-        hname=hnopt+"InvmassBeautyElecSecVtx";
-        fHistTagged.fInvmassBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
-        hname=hnopt+"InvmassNotBeautyElecSecVtx";
-        fHistTagged.fInvmassNotBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
-
-        hname=hnopt+"SignedLxyBeautyElecSecVtx";
-        fHistTagged.fSignedLxyBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
-        hname=hnopt+"SignedLxyNotBeautyElecSecVtx";
-        fHistTagged.fSignedLxyNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
-
-        hname=hnopt+"DistTwoVtxBeautyElecSecVtx";
-        fHistTagged.fDistTwoVtxBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
-        hname=hnopt+"DistTwoVtxNotBeautyElecSecVtx";
-        fHistTagged.fDistTwoVtxNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
-
-        hname=hnopt+"CosPhiBeautyElecSecVtx";
-        fHistTagged.fCosPhiBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
-        hname=hnopt+"CosPhiNotBeautyElecSecVtx";
-        fHistTagged.fCosPhiNotBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
-
-        hname=hnopt+"Chi2BeautyElecSecVtx";
-        fHistTagged.fChi2BeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
-        hname=hnopt+"Chi2NotBeautyElecSecVtx";
-        fHistTagged.fChi2NotBeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
-
-        hname=hnopt+"InvmassBeautyElec2trkVtx";
-        fHistTagged.fInvmassBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
-        hname=hnopt+"InvmassNotBeautyElec2trkVtx";
-        fHistTagged.fInvmassNotBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
-
-        hname=hnopt+"SignedLxyBeautyElec2trkVtx";
-        fHistTagged.fSignedLxyBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
-        hname=hnopt+"SignedLxyNotBeautyElec2trkVtx";
-        fHistTagged.fSignedLxyNotBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
-
-        hname=hnopt+"Chi2BeautyElec2trkVtx";
-        fHistTagged.fChi2BeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
-        hname=hnopt+"Chi2NotBeautyElec2trkVtx";
-        fHistTagged.fChi2NotBeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
+  // get KF particle input pid
+  Int_t pdg1 = GetPDG(track1);
+  Int_t pdg2 = GetPDG(track2);
+        
+  if(pdg1==-1 || pdg2==-1) {
+    //printf("out if considered pid range \n");
+    return;
+  }
+
+  // create KF particle of pair
+  if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
+  else AliKFParticle::SetField(fESD1->GetMagneticField()); 
+  AliKFParticle kfTrack1(*track1, pdg1);
+  AliKFParticle kfTrack2(*track2, pdg2);
+
+  AliKFParticle kfSecondary(kfTrack1,kfTrack2);
+
+  //secondary vertex point from kf particle
+  Double_t kfx = kfSecondary.GetX();
+  Double_t kfy = kfSecondary.GetY();
+  //Double_t kfz = kfSecondary.GetZ();
+        
+  //momentum at the decay point from kf particle
+  Double_t kfpx = kfSecondary.GetPx();
+  Double_t kfpy = kfSecondary.GetPy();
+  //Double_t kfpz = kfSecondary.GetPz();
+
+  Double_t dx = kfx-fPVx;
+  Double_t dy = kfy-fPVy;
+
+  // discriminating variables 
+  // invariant mass of the KF particle
+  Double_t invmass = -1;
+  Double_t invmassSigma = -1;
+  kfSecondary.GetMass(invmass,invmassSigma);
+
+  // chi2 of the KF particle
+  Double_t kfchi2 = -1;
+  if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+
+  // opening angle between two particles in XY plane
+  Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
+  Double_t cosphi = TMath::Cos(phi);
+
+  // projection of kf vertex vector to the kf momentum direction 
+  Double_t signedLxy=-999.;
+  Double_t psqr = kfpx*kfpx+kfpy*kfpy;
+  /*
+  //[the other way to think about]
+  if(psqr>0) {
+    if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);  
+    if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);  
+  }*/
+  if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);  
+
+  // DCA from primary to e-h KF particle (impact parameter of KF particle)
+  Double_t vtx[2]={fPVx, fPVy}; 
+  Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
+
+  /* //later consider the below
+  Float_t dcaR1=-999., dcaR2=-999.;
+  Float_t dcaZ1=-999., dcaZ2=-999.;
+
+  if(IsAODanalysis()){
+    Double_t trkIpPar1[2];
+    Double_t trkIpCov1[3];
+    Double_t trkIpPar2[2];
+    Double_t trkIpCov2[3];
+
+    AliESDtrack esdTrk1(track1);
+    AliESDtrack esdTrk2(track2);
+    esdTrk1.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar1,trkIpCov1);
+    esdTrk2.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar2,trkIpCov2);
+
+    dcaR1=trkIpPar1[0];
+    dcaZ1=trkIpPar1[1];
+    dcaR2=trkIpPar2[0];
+    dcaZ2=trkIpPar2[1];
+  }
+  else {
+    ((AliESDtrack*)track1)->GetImpactParameters(dcaR1,dcaZ1);
+    ((AliESDtrack*)track2)->GetImpactParameters(dcaR2,dcaZ2);
+  }*/
+
+  Int_t paircode = -1;
+  if (HasMCData()) paircode = GetPairCode(track1,track2); 
+
+  AliHFEpairs hfepair; 
+  hfepair.SetTrkLabel(index2);
+  hfepair.SetInvmass(invmass);
+  hfepair.SetKFChi2(kfchi2);
+  hfepair.SetOpenangle(phi);
+  hfepair.SetCosOpenangle(cosphi);
+  hfepair.SetSignedLxy(signedLxy);
+  hfepair.SetKFIP(kfip);
+  hfepair.SetPairCode(paircode);
+  AddHFEpairToArray(&hfepair);
+  fNoOfHFEpairs++; 
+
+  // fill into container for later QA
+  Double_t dataE[6];
+  dataE[0]=invmass;
+  dataE[1]=kfchi2;
+  dataE[2]=phi;
+  dataE[3]=signedLxy;
+  dataE[4]=kfip;
+  dataE[5]=paircode;
+  fPairQA->Fill(dataE);
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t index2)
+void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
 {
-        // calculate e-h pair characteristics and tag pair 
-
-        Int_t sourcePart = PairCode(track1,track2);
-
-        // get KF particle input pid
-        Int_t pdg1 = GetMCPID(track1);
-        Int_t pdg2 = GetMCPID(track2);
-        
-
-        // create KF particle of pair
-        AliKFParticle::SetField(fESD1->GetMagneticField());
-        AliKFParticle kfTrack1(*track1, pdg1);
-        AliKFParticle kfTrack2(*track2, pdg2);
-
-        AliKFParticle kfSecondary(kfTrack1,kfTrack2);
+  //
+  // run secondary vertexing algorithm and do tagging
+  //
 
-        // copy primary vertex from ESD
-        AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
-        if( primVtxCopy.GetNDF() <1 ) return;
-
-        //primary vertex point
-        Double_t pvx = primVtxCopy.GetX();
-        Double_t pvy = primVtxCopy.GetY();
-        //Double_t pvz = primVtxCopy.GetZ();
+  //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
+  FindSECVTXCandid(track);         
+}
 
-        //secondary vertex point from kf particle
-        Double_t kfx = kfSecondary.GetX();
-        Double_t kfy = kfSecondary.GetY();
-        //Double_t kfz = kfSecondary.GetZ();
-        
-        //momentum at the decay point from kf particle
-        Double_t kfpx = kfSecondary.GetPx();
-        Double_t kfpy = kfSecondary.GetPy();
-        //Double_t kfpz = kfSecondary.GetPz();
+//_______________________________________________________________________________________________
+void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track) 
+{
+  //
+  // find secondary vertex candidate and store them into container
+  //
+
+  AliVTrack *htrack[20];
+  Int_t htracklabel[20];
+  Double_t vtxchi2cut=3.; // testing cut 
+  Double_t dataE[6]={-999.,-999.,-999.,-999.,-1.,0};  
+  if (HFEpairs()->GetEntriesFast()>20){
+    AliDebug(3, "number of paired hadron is over maximum(20)");
+    return; 
+  }
         
-
-        Double_t dx = kfx-pvx;
-        Double_t dy = kfy-pvy;
-
-
-
-        // discriminating variables ----------------------------------------------------------
-
-        // invariant mass of the KF particle
-        Double_t invmass = -1;
-        Double_t invmassSigma = -1;
-        kfSecondary.GetMass(invmass,invmassSigma);
-
-        // chi2 of the KF particle
-        Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
-
-        // opening angle between two particles in XY plane
-        Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
-        Double_t cosphi = TMath::Cos(phi);
-
-        // projection of kf vertex vector to the kf momentum direction 
-        Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy);
-        Double_t signedLxy = TMath::Sqrt(dx*dx+dy*dy)*costheta;
-
-        // DCA from primary to e-h KF particle (impact parameter of KF particle)
-        Double_t vtx[2]={pvx, pvy}; 
-        Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
-
-
-             Float_t dcaR=-1; 
-        Float_t dcaR1=-1, dcaR2=-1;
-        Float_t dcaZ1=-1, dcaZ2=-1;
-        track1->GetImpactParameters(dcaR1,dcaZ1);
-        track2->GetImpactParameters(dcaR2,dcaZ2);
-
-             if (TMath::Abs(dcaR1) >= TMath::Abs(dcaR2)) dcaR=dcaR1;
-             else dcaR=dcaR2;
-
-        // fill histograms 
-        fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
-        fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
-        fHistPair[sourcePart].fOpenAngle->Fill(phi);
-        fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
-        fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
-        fHistPair[sourcePart].fKFIP->Fill(kfip);
-        fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR));
-
-        // pair cuts 
-        if( kfchi2 >2. ) return;
-
-        if ( signedLxy > 0.05 && cosphi > 0.5) fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
-
-        // pair tagging if it passed the above cuts
-        if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
-
-
-        // pair tagging condition
-        if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut
-        //if ( signedLxy > 0.06 && cosphi > 0) {
-          fpairedTrackID[fPairTagged] = index2;
-          fpairedChi2[fPairTagged] = kfchi2;
-          fpairedInvMass[fPairTagged] = invmass;
-          fpairedSignedLxy[fPairTagged] = signedLxy;
-          fPairTagged++;
+  // get paired track objects  
+  AliHFEpairs *pair=0x0;
+  if (IsAODanalysis()){
+    for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
+       pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
+       htracklabel[ip] = pair->GetTrkLabel();
+       htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
+    }
+  }
+  else{
+    for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
+       pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
+       htracklabel[ip] = pair->GetTrkLabel();
+       htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
+    }
+  }
+  // in case there is only one paired track with the electron, put pair characteristics into secvtx container
+  // for the moment, I only apply pair vertex chi2 cut
+  if (HFEpairs()->GetEntriesFast() == 1){
+    if (pair->GetKFChi2()<vtxchi2cut) { // you can also put single track cut
+      AliHFEsecVtxs hfesecvtx;
+      hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
+      hfesecvtx.SetTrkLabel2(-999);
+      hfesecvtx.SetInvmass(pair->GetInvmass());
+      hfesecvtx.SetKFChi2(pair->GetKFChi2());
+      hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
+      hfesecvtx.SetKFIP(pair->GetKFIP());
+      AddHFEsecvtxToArray(&hfesecvtx);
+      fNoOfHFEsecvtxs++; 
+
+      dataE[0]=pair->GetInvmass();
+      dataE[1]=pair->GetKFChi2();
+      dataE[2]=pair->GetSignedLxy();
+      dataE[3]=pair->GetKFIP();
+      if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
+      dataE[5]=2;
+      fSecvtxQA->Fill(dataE);
+    }
+    return;
+  }
+
+  // in case there are multiple paired track with the electron, calculate secvtx characteristics
+  // put the secvtx characteristics into container if it passes cuts
+  for (int i=0; i<HFEpairs()->GetEntriesFast()-1; i++){
+     for (int j=i+1; j<HFEpairs()->GetEntriesFast(); j++){
+        CalcSECVTXProperty(track, htrack[i], htrack[j]);
+        if (fKFchi2<vtxchi2cut) {
+          AliHFEsecVtxs hfesecvtx;
+          hfesecvtx.SetTrkLabel1(htracklabel[i]);
+          hfesecvtx.SetTrkLabel2(htracklabel[j]);
+          hfesecvtx.SetKFChi2(fKFchi2);
+          hfesecvtx.SetInvmass(fInvmass);
+          hfesecvtx.SetSignedLxy(fSignedLxy);
+          hfesecvtx.SetKFIP(fKFip);
+          AddHFEsecvtxToArray(&hfesecvtx);
+          fNoOfHFEsecvtxs++; 
+
+          dataE[0]=pair->GetInvmass();
+          dataE[1]=pair->GetKFChi2();
+          dataE[2]=pair->GetSignedLxy();
+          dataE[3]=pair->GetKFIP();
+          if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
+          dataE[5]=3;
+          fSecvtxQA->Fill(dataE);
         }
-
+     }
+  }
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
+void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
 {
-        // run secondary vertexing algorithm and do tagging
+  //
+  // calculate secondary vertex properties
+  //
+
+  // get KF particle input pid
+  Int_t pdg1 = GetPDG(track1);
+  Int_t pdg2 = GetPDG(track2);
+  Int_t pdg3 = GetPDG(track3);
+
+  if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
+    //printf("out if considered pid range \n");
+    return;
+  }
+
+  // create KF particle of pair
+  if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
+  else AliKFParticle::SetField(fESD1->GetMagneticField());
+  AliKFParticle kfTrack1(*track1, pdg1);
+  AliKFParticle kfTrack2(*track2, pdg2);
+  AliKFParticle kfTrack3(*track3, pdg3);
+
+  AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
+        
+  //secondary vertex point from kf particle
+  Double_t kfx = kfSecondary.GetX();
+  Double_t kfy = kfSecondary.GetY();
+  //Double_t kfz = kfSecondary.GetZ();
 
-        ResetTagVar();
+  //momentum at the decay point from kf particle
+  Double_t kfpx = kfSecondary.GetPx();
+  Double_t kfpy = kfSecondary.GetPy();
+  //Double_t kfpz = kfSecondary.GetPz();
 
-        Int_t imclabel = TMath::Abs(track->GetLabel());
-        if(imclabel<0) return;
-        TParticle* mcpart = fStack->Particle(imclabel);
-        Int_t esource = GetElectronSource(imclabel);
-        if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
-                fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
-               fBElec = kTRUE;
-        }
+  Double_t dx = kfx-fPVx;
+  Double_t dy = kfy-fPVy;
 
+  // discriminating variables ----------------------------------------------------------
 
-        if (fPairTagged >= 4) {
-          FindSECVTXCandid4Tracks(track);         
-        }
-        else if (fPairTagged == 3) {
-          FindSECVTXCandid3Tracks(track);         
-        }
-        else if (fPairTagged == 2) {
-          FindSECVTXCandid2Tracks(track);         
-        }
-        else if (fPairTagged == 1) {
-          ApplyPairTagCut();      
-        }
+  if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF())); 
 
+  // invariant mass of the KF particle
+  kfSecondary.GetMass(fInvmass,fInvmassSigma);
 
-        if (fBTagged) {
-                fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
-                if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
-                        fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
-                }
-                else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
-        }
+  // DCA from primary to e-h KF particle (impact parameter of KF particle)
+  Double_t vtx[2]={fPVx, fPVy};
+  fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
 
+  // projection of kf vertex vector to the kf momentum direction 
+  Double_t psqr = kfpx*kfpx+kfpy*kfpy;
+  if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);  
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::ApplyPairTagCut()
+Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track) 
 {
-        // apply tagging cut for e-h pair
+  //      
+  // return mc pid
+  //      
 
-        if (fBElec){
-                fHistTagged.fInvmassBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
-                fHistTagged.fSignedLxyBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
-                fHistTagged.fChi2BeautyElec2trkVtx->Fill(fpairedChi2[0]);
-        }
-        else {
-                fHistTagged.fInvmassNotBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
-                fHistTagged.fSignedLxyNotBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
-                fHistTagged.fChi2NotBeautyElec2trkVtx->Fill(fpairedChi2[0]);
-        }
-                
-        if (fpairedChi2[0] > 2.0) return;
-        if (fpairedInvMass[0] < 1.5) return;
-        if (fpairedSignedLxy[0] < 0.05) return;
+  Int_t label = TMath::Abs(track->GetLabel());
+  TParticle* mcpart = fStack->Particle(label);
+  if ( !mcpart ) return 0;
+  Int_t pdgCode = mcpart->GetPdgCode();
 
-        fBTagged = kTRUE;
+  return pdgCode;
 }
 
 //_______________________________________________________________________________________________
-Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id) 
+Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
 {
-        // apply tagging cut for e-h pair of given indexed hadron
-
-        if (fpairedChi2[id] > 2.0) return kFALSE;
-        if (fpairedInvMass[id] < 1.5) return kFALSE;
-        if (fpairedSignedLxy[id] < 0.05) return kFALSE;
-
-        fBTagged = kTRUE;
-        return kTRUE;
-                
+  //
+  // return pdg code of the origin(source) of the pair 
+  // 
+  //
+  // ---*---*---*-----ancester A----- track1
+  //                        |____*______ 
+  //                             |______ track2
+  // => if they originated from same ancester, 
+  //    then return "the absolute value of pdg code of ancester A"
+  //
+  // ---*---*---B hadron-----ancester A----- track1
+  //                               |____*______ 
+  //                                    |______ track2
+  // => if they originated from same ancester, and this ancester originally comes from B hadrons
+  //    then return -1*"the absolute value of pdg code of ancester A"
+  //
+  // caution : it can also return parton pdg code if it originated from same string or gluon spliting. 
+  //           
+
+  if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
+  TParticle* part1 = fStack->Particle(trk1->GetLabel());
+  TParticle* part2 = fStack->Particle(trk2->GetLabel());
+  TParticle* part2cp = part2;
+  if (!(part1) || !(part2)) return 0;
+
+  Int_t srcpdg = 0;
+
+  //if the two tracks' mother's label is same, get the mother info
+  //in case of charm, check if it originated from beauty
+  for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
+     Int_t label1 = part1->GetFirstMother(); 
+     if (label1 < 0) return 0;
+
+     for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
+        Int_t label2 = part2->GetFirstMother(); 
+        if (label2 < 0) break; 
+
+        if (label1 == label2){ //check if two tracks are originated from same mother
+          TParticle* commonmom = fStack->Particle(label2); 
+          srcpdg = abs(commonmom->GetPdgCode()); 
+
+          //check ancester to see if it is originally from beauty 
+          for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
+             Int_t ancesterlabel = commonmom->GetFirstMother();
+             if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg  
+
+             TParticle* commonancester = fStack->Particle(ancesterlabel);
+             Int_t ancesterpdg = abs(commonancester->GetPdgCode());
+
+             for (Int_t l=0; l<fNparents; l++){
+                if (abs(ancesterpdg)==fParentSelect[1][l]){
+                  srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
+                  return srcpdg;
+                }
+             }
+             commonmom = commonancester;
+          }
+        }
+        part2 = fStack->Particle(label2); //if their mother is different, go to earlier generation of 2nd particle
+        if (!(part2)) break;
+     }
+     part1 = fStack->Particle(label1); //if their mother is different, go to earlier generation of 1st particle
+     part2 = part2cp;
+     if (!(part1)) return 0;
+  }
+
+  return srcpdg; 
 }
 
 //_______________________________________________________________________________________________
-Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2) 
+Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
 {
-        // apply tagging cut for secondary vertex
 
-        if (chi2 > 2.0) return kFALSE;
-        if (finvmass < 1.5) return kFALSE;
-        if (fsignedLxy < 0.05) return kFALSE;
-        if (fcosPhi < 0.90) return kFALSE;
-        if (fdistTwoSecVtx > 0.1) return kFALSE;
-        
-        fBTagged = kTRUE;
-        return kTRUE;
+  //
+  // return pdg code of the origin(source) of the pair 
+  // 
+  //
+  // ---*---*---*-----ancester A----- track1
+  //                        |____*______ 
+  //                             |______ track2
+  // => if they originated from same ancester, 
+  //    then return "the absolute value of pdg code of ancester A"
+  //
+  // ---*---*---B hadron-----ancester A----- track1
+  //                               |____*______ 
+  //                                    |______ track2
+  // => if they originated from same ancester, and this ancester originally comes from B hadrons
+  //    then return -1*"the absolute value of pdg code of ancester A"
+  //
+  // caution : it can also return parton pdg code if it originated from same string or gluon spliting. 
+  //           
+
+  if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
+  AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
+  AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
+  AliAODMCParticle *part2cp = part2;
+  if (!(part1) || !(part2)) return 0;
+
+  Int_t srcpdg = 0;
+
+  //if the two tracks' mother's label is same, get the mother info
+  //in case of charm, check if it originated from beauty
+  for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
+     Int_t label1 = part1->GetMother(); 
+     if (label1 < 0) return 0;
+
+     for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
+        Int_t label2 = part2->GetMother(); 
+        if (label2 < 0) return 0; 
+
+        if (label1 == label2){ //check if two tracks are originated from same mother
+          AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
+          srcpdg = abs(commonmom->GetPdgCode()); 
+
+          //check ancester to see if it is originally from beauty 
+          for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
+             Int_t ancesterlabel = commonmom->GetMother();
+             if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg  
+
+             AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
+             Int_t ancesterpdg = abs(commonancester->GetPdgCode());
+
+             for (Int_t l=0; l<fNparents; l++){
+                if (abs(ancesterpdg)==fParentSelect[1][l]){
+                  srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
+                  return srcpdg;
+                }
+             }
+             commonmom = commonancester;
+          }
+        }
+        part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
+        if (!(part2)) break;
+     }
+     part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
+     part2 = part2cp;
+     if (!(part1)) return 0;
+  }
+
+  return srcpdg; 
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2)
+Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
 {
-        // apply tagging cut for e-h pair of given indexed hadron
+  //           
+  // return pair code which is predefinded as:
+  //  kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, 
+  //  kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
+  //           
+
+  Int_t srcpdg = -1;
+  Int_t srccode = kElse;
+
+  if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
+  else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
+
+  if (srcpdg < 0) srccode = kBeautyElse;
+  for (Int_t i=0; i<fNparents; i++){
+    if (abs(srcpdg)==fParentSelect[0][i]) {
+      if (srcpdg>0) srccode = kDirectCharm;
+      if (srcpdg<0) srccode = kBeautyCharm;
+    }
+    if (abs(srcpdg)==fParentSelect[1][i]) {
+      if (srcpdg>0) srccode = kDirectBeauty;
+      if (srcpdg<0)  return kElse;
+    }
+  }
+  if (srcpdg == 22) srccode = kGamma;
+  if (srcpdg == -22) srccode = kBeautyGamma;
+  if (srcpdg == 111) srccode = kPi0;
+  if (srcpdg == -111) srccode = kBeautyPi0;
 
-        if(!ApplyTagCut(chi2)){
-                if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2);
-        }
-        
+  return srccode;
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track) 
+Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack) 
 {
-        // find secondary vertex for >= 4 e-h pairs
+  //
+  // return decay electron's origin 
+  //
 
-        // sort pair in increasing order (kFALSE-increasing order)
-        Int_t index[20];
-        Int_t indexA[4];
-        Int_t indexB[3];
-        Double_t sevchi2[4];
-        AliESDtrack *htrack[4];
+  if (iTrack < 0) {
+    AliDebug(1, "Stack label is negative, return\n");
+    return -1;
+  }
 
-        if(fPairTagged>20) return; // protection
+  TParticle* mcpart = fStack->Particle(iTrack);
 
-        TMath::Sort(fPairTagged,fpairedChi2,index,kFALSE);
+//  if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
 
-        // select 4 partner tracks retruning smallest pair chi2 
+  Int_t iLabel = mcpart->GetFirstMother();
+  if (iLabel<0){
+    AliDebug(1, "Stack label is negative, return\n");
+    return -1;
+  }
 
-        for (Int_t i=0; i<4; i++){
-                htrack[i] = fESD1->GetTrack(fpairedTrackID[index[i]]);
-        }
-        
-        // calculated chi2 with 1 electron and 3 partner tracks 
-        for (Int_t i=0; i<4; i++){
-                sevchi2[i] = GetSecVtxChi2(track, htrack[fid[i][0]], htrack[fid[i][1]], htrack[fid[i][2]]);
-        }
+  Int_t origin = -1;
+  Bool_t isFinalOpenCharm = kFALSE;
 
-        // select 3 partner tracks retruning smallest pair chi2
-        // [think] if two smallest chi2 are similar, have to think about better handling of selection
-        TMath::Sort(4,sevchi2,indexA,kFALSE);
+  TParticle *partMother = fStack->Particle(iLabel);
+  Int_t maPdgcode = partMother->GetPdgCode();
 
-        // calculated chi2 with 1 electron and 2 partner tracks 
-        for (Int_t i=0; i<3; i++){
-                sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
-        }
+  // if the mother is charmed hadron  
+  if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 
-        // select 2 partner tracks retruning smallest pair chi2
-        TMath::Sort(3,sevchi2,indexB,kFALSE);
+    for (Int_t i=0; i<fNparents; i++){
+      if (abs(maPdgcode)==fParentSelect[0][i]){
+        isFinalOpenCharm = kTRUE;
+      }
+    }
+    if (!isFinalOpenCharm) return -1;
 
-        // calculate secondary vertex quality variables with 1 electron and 2 hadrons
-        CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
+    // iterate until you find B hadron as a mother or become top ancester 
+    for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 
-        if (fBElec){
-                fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
-        }
-        else {
-                fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+      Int_t jLabel = partMother->GetFirstMother();
+      if (jLabel == -1){
+        origin = kDirectCharm;
+        return origin;
+      }
+      if (jLabel < 0){ // safety protection
+        AliDebug(1, "Stack label is negative, return\n");
+        return -1;
+      }
+
+      // if there is an ancester
+      TParticle* grandMa = fStack->Particle(jLabel);
+      Int_t grandMaPDG = grandMa->GetPdgCode();
+
+      for (Int_t j=0; j<fNparents; j++){
+        if (abs(grandMaPDG)==fParentSelect[1][j]){
+          origin = kBeautyCharm;
+          return origin;
         }
+      }
+
+      partMother = grandMa;
+    } // end of iteration 
+  } // end of if
+  else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
+    for (Int_t i=0; i<fNparents; i++){
+      if (abs(maPdgcode)==fParentSelect[1][i]){
+        origin = kDirectBeauty;
+        return origin;
+      }
+    }
+  } // end of if
 
-        ApplyVtxTagCut(sevchi2[indexB[0]],index[fia[indexA[0]][indexB[0]][0]],index[fia[indexA[0]][indexB[0]][1]]);
+  //============ gamma ================
+  else if ( abs(maPdgcode) == 22 ) {
+    origin = kGamma;
 
-}
-
-//_______________________________________________________________________________________________
-void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track) 
-{
-        // find secondary vertex for 3 e-h pairs
+    // iterate until you find B hadron as a mother or become top ancester 
+    for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 
-        // sort pair in increasing order (kFALSE-increasing order)
-        Int_t indexA[1] = { 0 };
-        Int_t indexB[3];
-        Double_t sevchi2[3];
-        AliESDtrack *htrack[3];
+      Int_t jLabel = partMother->GetFirstMother();
+      if (jLabel == -1){
+        origin = kGamma;
+        return origin;
+      }
+      if (jLabel < 0){ // safety protection
+        AliDebug(1, "Stack label is negative, return\n");
+        return -1;
+      }
+
+      // if there is an ancester
+      TParticle* grandMa = fStack->Particle(jLabel);
+      Int_t grandMaPDG = grandMa->GetPdgCode();
+
+      for (Int_t j=0; j<fNparents; j++){
+        if (abs(grandMaPDG)==fParentSelect[1][j]){
+          origin = kBeautyGamma;
+          return origin;
+        }
+      }
 
-        // select 4 partner tracks retruning smallest pair chi2 
+      partMother = grandMa;
+    } // end of iteration 
 
-        for (Int_t i=0; i<3; i++){
-                htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
-        }
-        
-        // calculated chi2 with 1 electron and 2 partner tracks 
-        for (Int_t i=0; i<3; i++){
-                sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
-        }
+    return origin;
+  } // end of if
 
-        // select 2 partner tracks retruning smallest pair chi2
-        TMath::Sort(3,sevchi2,indexB,kFALSE);
+  //============ pi0 ================
+  else if ( abs(maPdgcode) == 111 ) {
+    origin = kPi0;
 
-        // calculate secondary vertex quality variables with 1 electron and 2 hadrons
-        CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
+    // iterate until you find B hadron as a mother or become top ancester 
+    for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 
-        if (fBElec){
-                fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
-        }
-        else {
-                fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+      Int_t jLabel = partMother->GetFirstMother();
+      if (jLabel == -1){
+        origin = kPi0;
+        return origin;
+      }
+      if (jLabel < 0){ // safety protection
+        AliDebug(1, "Stack label is negative, return\n");
+        return -1;
+      }
+
+      // if there is an ancester
+      TParticle* grandMa = fStack->Particle(jLabel);
+      Int_t grandMaPDG = grandMa->GetPdgCode();
+
+      for (Int_t j=0; j<fNparents; j++){
+        if (abs(grandMaPDG)==fParentSelect[1][j]){
+          origin = kBeautyPi0;
+          return origin;
         }
+      }
 
-        ApplyVtxTagCut(sevchi2[indexB[0]],fia[indexA[0]][indexB[0]][0],fia[indexA[0]][indexB[0]][1]);
-}
+      partMother = grandMa;
+    } // end of iteration 
 
-//_______________________________________________________________________________________________
-void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track) 
-{
-        // find secondary vertex for 2 e-h pairs
+    return origin;
+  } // end of if
 
-        Double_t sevchi2[1];
-        AliESDtrack *htrack[2];
+  else {
+    origin = kElse;
+    // iterate until you find B hadron as a mother or become top ancester 
+    for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 
-        for (Int_t i=0; i<2; i++){
-                htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
+      Int_t jLabel = partMother->GetFirstMother();
+      if (jLabel == -1){
+        origin = kElse;
+        return origin;
+      }
+      if (jLabel < 0){ // safety protection
+        AliDebug(1, "Stack label is negative, return\n");
+        return -1;
+      }
+
+      // if there is an ancester
+      TParticle* grandMa = fStack->Particle(jLabel);
+      Int_t grandMaPDG = grandMa->GetPdgCode();
+
+      for (Int_t j=0; j<fNparents; j++){
+        if (abs(grandMaPDG)==fParentSelect[1][j]){
+          origin = kBeautyElse;
+          return origin;
         }
-        
-        sevchi2[0] = GetSecVtxChi2(track, htrack[0], htrack[1]);
-
-        // calculate secondary vertex quality variables with 1 electron and 2 hadrons
-        CalcSECVTXProperty(track,htrack[0],htrack[1]);
+      }
 
-        if (fBElec){
-                fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[0]);
-        }
-        else {
-                fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[0]);
-        }
+      partMother = grandMa;
+    } // end of iteration 
+  }
 
-        ApplyVtxTagCut(sevchi2[0],0,1);
+  return origin;
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
+Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
 {
-        // calculate secondary vertex properties
-
-        Int_t pdg1 = GetMCPID(track1);
-        Int_t pdg2 = GetMCPID(track2);
-        Int_t pdg3 = GetMCPID(track3);
-
-        AliKFParticle::SetField(fESD1->GetMagneticField());
-        AliKFParticle kfTrack1(*track1, pdg1);
-        AliKFParticle kfTrack2(*track2, pdg2);
-        AliKFParticle kfTrack3(*track3, pdg3);
-
-        AliKFParticle kfSecondary12(kfTrack1,kfTrack2);
-        AliKFParticle kfSecondary13(kfTrack1,kfTrack3);
-        AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
-
-        // copy primary vertex from ESD
-        AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));        
-        //printf("primary ndf= %d\n",primVtxCopy.GetNDF());
-        if( primVtxCopy.GetNDF() <1 ) return;
-
-        Double_t kdx12 = kfSecondary12.GetX()-primVtxCopy.GetX();
-        Double_t kdy12 = kfSecondary12.GetY()-primVtxCopy.GetY();
-        //Double_t kdz12 = kfSecondary12.GetZ()-primVtxCopy.GetZ();
-
-        Double_t kdx13 = kfSecondary13.GetX()-primVtxCopy.GetX();
-        Double_t kdy13 = kfSecondary13.GetY()-primVtxCopy.GetY();
-        //Double_t kdz13 = kfSecondary13.GetZ()-primVtxCopy.GetZ();
-
-        Double_t kdx = kfSecondary.GetX()-primVtxCopy.GetX();
-        Double_t kdy = kfSecondary.GetY()-primVtxCopy.GetY();
-        //Double_t kdz = kfSecondary.GetZ()-primVtxCopy.GetZ();
-
-        // calculate distance and angle between two secvtxes 
-        fdistTwoSecVtx = TMath::Sqrt((kdx12-kdx13)*(kdx12-kdx13) + (kdy12-kdy13)*(kdy12-kdy13));
-        fcosPhi = ( kdx12*kdx13 + kdy12*kdy13 ) / ( TMath::Sqrt(kdx12*kdx12+kdy12*kdy12)*TMath::Sqrt(kdx13*kdx13+kdy13*kdy13) );
-        //Double_t lengthdiff = TMath::Abs(TMath::Sqrt(kdx12*kdx12+kdy12*kdy12) - TMath::Sqrt(kdx13*kdx13+kdy13*kdy13));
-
-        // calculate angle between secondary vertex vector and secondary particle momentum vector in transverse plane
-        Double_t cosTheta = ( kdx*kfSecondary.GetPx() + kdy*kfSecondary.GetPy()) / TMath::Sqrt(kdx*kdx+kdy*kdy)*TMath::Sqrt(kfSecondary.GetPx()*kfSecondary.GetPx()+kfSecondary.GetPy()*kfSecondary.GetPy());
-        // calculate signed Lxy
-        fsignedLxy = TMath::Sqrt(kdx*kdx+kdy*kdy)*cosTheta;
-
-        // calculate invariant mass of the kf particle
-        kfSecondary.GetMass(finvmass,finvmassSigma);
-
-       if (fBElec){
-               fHistTagged.fInvmassBeautyElecSecVtx->Fill(finvmass);
-               fHistTagged.fSignedLxyBeautyElecSecVtx->Fill(fsignedLxy);
-               fHistTagged.fDistTwoVtxBeautyElecSecVtx->Fill(fdistTwoSecVtx);
-               fHistTagged.fCosPhiBeautyElecSecVtx->Fill(fcosPhi);
-       }
-       else {
-               fHistTagged.fInvmassNotBeautyElecSecVtx->Fill(finvmass);
-               fHistTagged.fSignedLxyNotBeautyElecSecVtx->Fill(fsignedLxy);
-               fHistTagged.fDistTwoVtxNotBeautyElecSecVtx->Fill(fdistTwoSecVtx);
-               fHistTagged.fCosPhiNotBeautyElecSecVtx->Fill(fcosPhi);
-       }
+  //
+  // get KF particle input pdg for mass hypothesis
+  //
+
+  Int_t pdgCode=-1;
+
+  if (fUseMCPID && HasMCData()){
+    pdgCode = GetMCPDG(track);
+  }
+  else if(fESD1){
+    Int_t pid=0;
+    Double_t prob;
+    GetESDPID((AliESDtrack*)track, pid, prob);
+    switch(pid){
+      case 0:  pdgCode = 11; break;
+      case 1:  pdgCode = 13; break;
+      case 2:  pdgCode = 211; break;
+      case 3:  pdgCode = 321; break;
+      case 4:  pdgCode = 2212; break;
+      default: pdgCode = -1;
+    }
+  }
+  else if(fAOD1){
+    Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
+    switch(pid){
+      case 0:  pdgCode = 11; break;
+      case 1:  pdgCode = 13; break;
+      case 2:  pdgCode = 211; break;
+      case 3:  pdgCode = 321; break;
+      case 4:  pdgCode = 2212; break;
+      default: pdgCode = -1;
+    }
+  }
 
+  return pdgCode;
 }
 
 //_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track) 
+Int_t AliHFEsecVtx::GetMCPDG(AliVTrack *track)
 {
-        // return mc pid
-
-        Int_t label = TMath::Abs(track->GetLabel());
-        TParticle* mcpart = fStack->Particle(label);
-        if ( !mcpart ) return 0;
-        Int_t pdgCode = mcpart->GetPdgCode();
-
-        return pdgCode;
+  //
+  // return mc pdg code
+  //
+
+  Int_t label = TMath::Abs(track->GetLabel());
+  Int_t pdgCode; 
+
+  if (IsAODanalysis()) {
+    AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
+    if ( !mcpart ) return 0;
+      pdgCode = mcpart->GetPdgCode();
+  }
+  else {
+    TParticle* mcpart = fStack->Particle(label);
+    if ( !mcpart ) return 0;
+      pdgCode = mcpart->GetPdgCode();
+  }
+
+    return pdgCode;
 }
 
-//_______________________________________________________________________________________________
-Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
+//______________________________________________________________________________
+void AliHFEsecVtx::GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob) 
 {
-        // return 4 track secondary vertex chi2
-
-        Int_t pdg1 = GetMCPID(track1);
-        Int_t pdg2 = GetMCPID(track2);
-        Int_t pdg3 = GetMCPID(track3);
-        Int_t pdg4 = GetMCPID(track4);
+  //
+  // calculate likehood for esd pid
+  // 
 
-        AliKFParticle::SetField(fESD1->GetMagneticField());
-        AliKFParticle kfTrack1(*track1, pdg1);
-        AliKFParticle kfTrack2(*track2, pdg2);
-        AliKFParticle kfTrack3(*track3, pdg3);
-        AliKFParticle kfTrack4(*track4, pdg4);
-        AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3,kfTrack4);
+  recpid = -1;
+  recprob = -1;
 
-        return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+  Int_t ipid=-1;
+  Double_t max=0.;
 
-}
-
-//_______________________________________________________________________________________________
-Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
-{
-        // return 3 track secondary vertex chi2
+  Double_t probability[5];
 
-        Int_t pdg1 = GetMCPID(track1);
-        Int_t pdg2 = GetMCPID(track2);
-        Int_t pdg3 = GetMCPID(track3);
+  // get probability of the diffenrent particle types
+  track->GetESDpid(probability);
 
-        AliKFParticle::SetField(fESD1->GetMagneticField());
-        AliKFParticle kfTrack1(*track1, pdg1);
-        AliKFParticle kfTrack2(*track2, pdg2);
-        AliKFParticle kfTrack3(*track3, pdg3);
-        AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
-
-        return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+  // find most probable particle in ESD pid
+  // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
+  ipid = TMath::LocMax(5,probability);
+  max = TMath::MaxElement(5,probability);
 
+  recpid = ipid;
+  recprob = max;
 }
 
-//_______________________________________________________________________________________________
-Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
+//_____________________________________________________________________________
+void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
 {
-        // return 2 track secondary vertex chi2
-
-        Int_t pdg1 = GetMCPID(track1);
-        Int_t pdg2 = GetMCPID(track2);
-
-        AliKFParticle::SetField(fESD1->GetMagneticField());
-        AliKFParticle kfTrack1(*track1, pdg1);
-        AliKFParticle kfTrack2(*track2, pdg2);
-        AliKFParticle kfSecondary(kfTrack1,kfTrack2);
-
-        return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+  //
+  // Add a HFE pair to the array
+  //
 
+  Int_t n = HFEpairs()->GetEntriesFast();
+  if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
+  new((*HFEpairs())[n]) AliHFEpairs(*pair);
 }
 
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2)
+//_____________________________________________________________________________
+TClonesArray *AliHFEsecVtx::HFEpairs()
 {
-
-        //
-        // return pdg code of the origin(source) of the pair 
-        // 
-        //
-        // ---*---*---*-----ancester A----- track1
-        //                        |____*______ 
-        //                             |______ track2
-        // => if they originated from same ancester, 
-        //    then return "the absolute value of pdg code of ancester A"
-        //
-        // ---*---*---B hadron-----ancester A----- track1
-        //                               |____*______ 
-        //                                    |______ track2
-        // => if they originated from same ancester, and this ancester originally comes from B hadrons
-        //    then return -1*"the absolute value of pdg code of ancester A"
-        //
-        // caution : it can also return parton pdg code if it originated from same string or gluon spliting. 
-        //           
-
-        if (track1->GetLabel()<0 || track2->GetLabel()<0) return 0;
-        TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
-        TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
-        if (!(part1)) return 0;
-        if (!(part2)) return 0;
-
-        TParticle* part1Crtgen = part1; // copy track into current generation particle
-        TParticle* part2Crtgen = part2; // copy track into current generation particle
-
-
-        Int_t sourcePDG = 0;
-
-        // if the two tracks' mother's label is same, get the mother info
-        // in case of charm, check if it originated from beauty
-        for (Int_t i=0; i<100; i++){ // iterate 100
-                Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
-                if (iLabel < 0) return 0;
-
-                for (Int_t j=0; j<100; j++){ // iterate 100
-                        Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
-                        if (jLabel < 0) return 0; // if jLabel == -1
-
-                        if (iLabel == jLabel){ // check if two tracks are originated from same mother
-                                TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info 
-                                sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
-
-                                // check ancester to see if it is originally from beauty 
-                                for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
-                                        Int_t ancesterLabel = thismother->GetFirstMother();
-                                        if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg  
-
-                                        TParticle* thisancester = fStack->Particle(ancesterLabel);
-                                        Int_t ancesterPDG = abs(thisancester->GetPdgCode());
-
-                                        for (Int_t l=0; l<fNparents; l++){
-                                                if (abs(ancesterPDG)==fParentSelect[1][l]){
-                                                        sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
-                                                        return sourcePDG;
-                                                }
-                                        }
-                                        thismother = thisancester;
-                                }
-
-                        }
-                        part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
-                }
-                part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
-
-                // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
-                Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
-                for (Int_t l=0; l<fNparents; l++){
-                        if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
-                }
-
-        }
-
-        return sourcePDG;
-
+  //
+  // Returns the list of HFE pairs
+  //
+
+  if (!fHFEpairs) {
+      fHFEpairs = new TClonesArray("AliHFEpairs", 200);
+  }
+  return fHFEpairs;
 }
 
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
+//_____________________________________________________________________________
+void AliHFEsecVtx::DeleteHFEpairs()
 {
-
-        //           
-        // return pair code which is predefinded as:
-        //  kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
-        //           
-
-        Int_t pairOriginsPDG = PairOrigin(track1,track2);
-
-        Int_t sourcePart = kElse;
-
-        if (pairOriginsPDG < 0) {
-                sourcePart = kBeautyElse;
-        }
-        for (Int_t i=0; i<fNparents; i++){
-                if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
-                        if (pairOriginsPDG>0) sourcePart = kDirectCharm;
-                        if (pairOriginsPDG<0) {
-                                sourcePart = kBeautyCharm;
-                        }
-                }
-                if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
-                        if (pairOriginsPDG>0) {
-                                sourcePart = kDirectBeauty;
-                        }
-                        if (pairOriginsPDG<0)  return kElse;
-                }
-        }
-        if (pairOriginsPDG == 22) sourcePart = kGamma;
-        if (pairOriginsPDG == -22) {
-                sourcePart = kBeautyGamma;
-        }
-        if (pairOriginsPDG == 111) sourcePart = kPi0;
-        if (pairOriginsPDG == -111) {
-                sourcePart = kBeautyPi0;
-        }
-
-        return sourcePart;
-
+  //
+  // Delete the list of HFE pairs
+  //
+
+  if (fHFEpairs) {
+    fHFEpairs->Delete();
+    //delete fHFEpairs;
+  }
 }
 
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack) 
+//_____________________________________________________________________________
+void AliHFEsecVtx::InitHFEpairs()
 {
+  //
+  // Initialization should be done before make all possible pairs for a given electron candidate
+  //
 
-    // return decay electron's origin 
-
-    if (iTrack < 0) {
-      AliDebug(1, "Stack label is negative, return\n");
-      return -1;
-    }
-
-    TParticle* mcpart = fStack->Particle(iTrack);
-
-    if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !!!!!!!!!!!!!!!!!
-
-    Int_t iLabel = mcpart->GetFirstMother();
-    if (iLabel<0){
-      AliDebug(1, "Stack label is negative, return\n");
-      return -1;
-    }
-
-    Int_t origin = -1;
-    Bool_t isFinalOpenCharm = kFALSE;
-
-    TParticle *partMother = fStack->Particle(iLabel);
-    Int_t maPdgcode = partMother->GetPdgCode();
-
-    // if the mother is charmed hadron  
-    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
-
-         for (Int_t i=0; i<fNparents; i++){
-            if (abs(maPdgcode)==fParentSelect[0][i]){
-              isFinalOpenCharm = kTRUE;
-            }
-         }
-         if (!isFinalOpenCharm) return -1;
-
-
-          // iterate until you find B hadron as a mother or become top ancester 
-          for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
-             Int_t jLabel = partMother->GetFirstMother();
-             if (jLabel == -1){
-               origin = kDirectCharm;
-               return origin;
-             }
-             if (jLabel < 0){ // safety protection
-               AliDebug(1, "Stack label is negative, return\n");
-               return -1;
-             }
-
-             // if there is an ancester
-             TParticle* grandMa = fStack->Particle(jLabel);
-             Int_t grandMaPDG = grandMa->GetPdgCode();
-
-             for (Int_t j=0; j<fNparents; j++){
-                if (abs(grandMaPDG)==fParentSelect[1][j]){
-
-                  origin = kBeautyCharm;
-                  return origin;
-                }
-             }
-
-             partMother = grandMa;
-          } // end of iteration 
-    } // end of if
-    else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
-         for (Int_t i=0; i<fNparents; i++){
-            if (abs(maPdgcode)==fParentSelect[1][i]){
-              origin = kDirectBeauty;
-              return origin;
-            }
-         }
-    } // end of if
-
-
-    //============ gamma ================
-    else if ( abs(maPdgcode) == 22 ) {
-         origin = kGamma;
-
-          // iterate until you find B hadron as a mother or become top ancester 
-          for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
-             Int_t jLabel = partMother->GetFirstMother();
-             if (jLabel == -1){
-               origin = kGamma;
-               return origin;
-             }
-             if (jLabel < 0){ // safety protection
-               AliDebug(1, "Stack label is negative, return\n");
-               return -1;
-             }
-
-             // if there is an ancester
-             TParticle* grandMa = fStack->Particle(jLabel);
-             Int_t grandMaPDG = grandMa->GetPdgCode();
-
-             for (Int_t j=0; j<fNparents; j++){
-                if (abs(grandMaPDG)==fParentSelect[1][j]){
-                  origin = kBeautyGamma;
-                  return origin;
-                }
-             }
-
-             partMother = grandMa;
-          } // end of iteration 
-
-         return origin;
-    } // end of if
-
-    //============ pi0 ================
-    else if ( abs(maPdgcode) == 111 ) {
-         origin = kPi0;
-
-          // iterate until you find B hadron as a mother or become top ancester 
-          for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
-             Int_t jLabel = partMother->GetFirstMother();
-             if (jLabel == -1){
-               origin = kPi0;
-               return origin;
-             }
-             if (jLabel < 0){ // safety protection
-               AliDebug(1, "Stack label is negative, return\n");
-               return -1;
-             }
-
-             // if there is an ancester
-             TParticle* grandMa = fStack->Particle(jLabel);
-             Int_t grandMaPDG = grandMa->GetPdgCode();
-
-             for (Int_t j=0; j<fNparents; j++){
-                if (abs(grandMaPDG)==fParentSelect[1][j]){
-                  origin = kBeautyPi0;
-                  return origin;
-                }
-             }
-
-             partMother = grandMa;
-          } // end of iteration 
-
-         return origin;
-    } // end of if
-
-
-    else {
-        origin = kElse;
-
-          // iterate until you find B hadron as a mother or become top ancester 
-          for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
-             Int_t jLabel = partMother->GetFirstMother();
-             if (jLabel == -1){
-               origin = kElse;
-               return origin;
-             }
-             if (jLabel < 0){ // safety protection
-               AliDebug(1, "Stack label is negative, return\n");
-               return -1;
-             }
-
-             // if there is an ancester
-             TParticle* grandMa = fStack->Particle(jLabel);
-             Int_t grandMaPDG = grandMa->GetPdgCode();
-
-             for (Int_t j=0; j<fNparents; j++){
-                if (abs(grandMaPDG)==fParentSelect[1][j]){
-                  origin = kBeautyElse;
-                  return origin;
-                }
-             }
-
-             partMother = grandMa;
-          } // end of iteration 
-
-    }
-
-    return origin;
-
+  fNoOfHFEpairs = 0;
 }
 
-/*
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) 
+//_____________________________________________________________________________
+void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
 {
+  //
+  // Add a HFE secondary vertex to the array
+  //
 
-        //           
-        // decay electron's origin 
-        //           
-             
-        if (iTrackLabel < 0) {
-                AliDebug(1, "Stack label is negative, return\n");
-                return -1; 
-        }           
-
-        TParticle* mcpart = fStack->Particle(iTrackLabel);
-        Int_t iLabel = mcpart->GetFirstMother();
-        if (iLabel<0){
-                AliDebug(1, "Stack label is negative, return\n");
-                return -1;
-        }
-
-        Int_t origin = -1;
-        Bool_t isFinalOpenCharm = kFALSE;
-
-        TParticle *partMother = fStack->Particle(iLabel);
-        Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
-
-        //beauty --------------------------
-        if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
-                for (Int_t i=0; i<fNparents; i++){
-                        if (abs(maPdgcode)==fParentSelect[1][i]){
-                                origin = kDirectBeauty;
-                                return origin;
-                        }
-                        else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
-                }
-        } // end of if
+  Int_t n = HFEsecvtxs()->GetEntriesFast();
+  if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
+  new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
+}
 
-        //charm --------------------------
-        else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
-        
-                for (Int_t i=0; i<fNparents; i++){
-                        if (abs(maPdgcode)==fParentSelect[0][i])
-                                isFinalOpenCharm = kTRUE;
-                }
-                if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms   
-                                                  // to prevent any possible double counting  
-
-                for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
-
-                        Int_t jLabel = partMother->GetFirstMother();
-                        if (jLabel == -1){
-                                origin = kDirectCharm;
-                                return origin;
-                        }
-                        if (jLabel < 0){ // safety protection even though not really necessary here
-                                AliDebug(1, "Stack label is negative, return\n");
-                                return -1;
-                        }
-
-                        // if there is an ancester, check if it in the final B hadron list 
-                        TParticle* grandMa = fStack->Particle(jLabel);
-                        Int_t grandMaPDG = grandMa->GetPdgCode();
-
-                        for (Int_t j=0; j<fNparents; j++){
-                                if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                origin = kBeautyCharm;
-                                return origin;
-                                }
-                        } 
-
-                        partMother = grandMa;
-                } // end of iteration 
-        } // end of if
-
-        //gamma --------------------------
-        else if ( abs(maPdgcode) == 22 ) {
-                origin = kGamma;
-
-                // iterate until you find B hadron as a mother or become top ancester 
-                for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
-
-                        Int_t jLabel = partMother->GetFirstMother();
-                        if (jLabel == -1){
-                                origin = kGamma;
-                                return origin;
-                        }
-                        if (jLabel < 0){ // safety protection
-                                AliDebug(1, "Stack label is negative, return\n");
-                                return -1;
-                        }
-
-                        // if there is an ancester
-                        TParticle* grandMa = fStack->Particle(jLabel);
-                        Int_t grandMaPDG = grandMa->GetPdgCode();
-
-                        for (Int_t j=0; j<fNparents; j++){
-                                if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                        origin = kBeautyGamma;
-                                        return origin;
-                                }
-                        }
-
-                        partMother = grandMa;
-                } // end of iteration 
-
-                return origin;
-        } // end of if
-
-        //pi0 --------------------------
-        else if ( abs(maPdgcode) == 111 ) {
-                origin = kPi0;
-
-                // iterate until you find B hadron as a mother or become top ancester 
-                for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
-
-                        Int_t jLabel = partMother->GetFirstMother();
-                        if (jLabel == -1){
-                                origin = kPi0;
-                                return origin;
-                        }
-                        if (jLabel < 0){ // safety protection
-                                AliDebug(1, "Stack label is negative, return\n");
-                                return -1;
-                        }
-
-                        // if there is an ancester
-                        TParticle* grandMa = fStack->Particle(jLabel);
-                        Int_t grandMaPDG = grandMa->GetPdgCode();
-
-                        for (Int_t j=0; j<fNparents; j++){
-                                if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                        origin = kBeautyPi0;
-                                        return origin;
-                                }
-                        }
-
-                        partMother = grandMa;
-                } // end of iteration 
-
-                return origin;
-        } // end of if
-
-
-        //else --------------------------
-        else {
-                origin = kElse;
-
-                // iterate until you find B hadron as a mother or become top ancester 
-                for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
-
-                        Int_t jLabel = partMother->GetFirstMother();
-                        if (jLabel == -1){
-                                origin = kElse;
-                                return origin;
-                        }
-                        if (jLabel < 0){ // safety protection
-                                AliDebug(1, "Stack label is negative, return\n");
-                                return -1;
-                        }
-
-                        // if there is an ancester
-                        TParticle* grandMa = fStack->Particle(jLabel);
-                        Int_t grandMaPDG = grandMa->GetPdgCode();
-
-                        for (Int_t j=0; j<fNparents; j++){
-                                if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                        origin = kBeautyElse;
-                                        return origin;
-                                }
-                        }
-
-                        partMother = grandMa;
-                } // end of iteration 
+//_____________________________________________________________________________
+TClonesArray *AliHFEsecVtx::HFEsecvtxs()
+{
+  //
+  // Returns the list of HFE secvtx
+  //
+
+  if (!fHFEsecvtxs) {
+      fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
+  }
+  return fHFEsecvtxs;
+}
 
-        }
+//_____________________________________________________________________________
+void AliHFEsecVtx::DeleteHFEsecvtxs()
+{
+  //
+  // Delete the list of HFE pairs
+  //
+
+  if (fHFEsecvtxs) {
+    fHFEsecvtxs->Delete();
+    //delete fHFEsecvtx;
+  }
+}
 
-        return origin;
+//_____________________________________________________________________________
+void AliHFEsecVtx::InitHFEsecvtxs()
+{
+  //
+  // Initialization should be done 
+  //
 
+  fNoOfHFEsecvtxs = 0;
 }
-*/
 
 //_______________________________________________________________________________________________
 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
 {
-        // test cuts 
-
-        //if (track->Pt() < 1.0) return kFALSE;
-        //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
-        //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
-        //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
-        if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
-        //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
-
+  //if (track->Pt() < 1.0) return kFALSE;
+  //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
+  //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
+  //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
+  if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
+  //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
 
 /*
-        Float_t dcaR=-1;
-        Float_t dcaZ=-1;
-        track->GetImpactParameters(dcaR,dcaZ);
-        if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
-        if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
+  Float_t dcaR=-1;
+  Float_t dcaZ=-1;
+  track->GetImpactParameters(dcaR,dcaZ);
+  if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
+  if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
 */
-        return kTRUE;
+  return kTRUE;
+}
+//____________________________________________________________
+void AliHFEsecVtx::MakeContainer(){
+
+  //
+  // make container
+  //
+
+  const Int_t nDimPair=6;
+  Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11};
+  const Double_t kInvmassmin = 0., kInvmassmax = 20.;
+  const Double_t kKFChi2min = 0, kKFChi2max= 50;
+  const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
+  const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
+  const Double_t kKFIPmin = -10, kKFIPmax= 10;
+  const Double_t kPairCodemin = -1, kPairCodemax= 10;
+
+  Double_t* binEdgesPair[nDimPair];
+  for(Int_t ivar = 0; ivar < nDimPair; ivar++)
+    binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
+
+  for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
+  for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
+  for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
+  for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
+  for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
+  for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
+
+  fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code ", nDimPair, nBinPair);
+  for(Int_t idim = 0; idim < nDimPair; idim++){
+    fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
+  }
+
+  fSecVtxList->AddAt(fPairQA,0);
+
+  const Int_t nDimSecvtx=6;
+  Double_t* binEdgesSecvtx[nDimSecvtx];
+  Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 11, 3};
+  const Double_t kNtrksmin = 0, kNtrksmax= 3;
+  for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
+    binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
+
+  for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
+  for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
+  for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
+  for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
+  for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
+  for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
+
+  fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
+  for(Int_t idim = 0; idim < nDimSecvtx; idim++){
+    fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
+  }
+
+  fSecVtxList->AddAt(fSecvtxQA,1);
 }
index 50a26e9..0772917 100644 (file)
 //#include <TObject.h>
 #endif
 
+#ifndef ROOT_THnSparse
+#include <THnSparse.h>
+#endif
+
 class TH1F;
 class TH2F;
 class TString;
 class AliESDEvent;
+class AliAODEvent;
+class AliVTrack;
 class AliESDtrack;
+class AliAODTrack;
 class AliStack;
+class AliHFEpairs;
+class AliHFEsecVtxs;
 
 //________________________________________________________________
 class AliHFEsecVtx : public TObject {
@@ -41,197 +50,96 @@ class AliHFEsecVtx : public TObject {
                 AliHFEsecVtx &operator=(const AliHFEsecVtx &); // assignment operator
                 virtual ~AliHFEsecVtx();
 
-                void CreateHistograms(TString hnopt="");
+                void CreateHistograms(TList *qaList);
+
+                Bool_t HasMCData() const { return TestBit(kHasMCData); };
+                Bool_t IsAODanalysis() const { return TestBit(kAODanalysis); };
+                Bool_t IsESDanalysis() const { return !TestBit(kAODanalysis); };
+
+                           void SetHasMCData(Bool_t hasMCdata = kTRUE) { SetBit(kHasMCData,hasMCdata); };
+                void SetAODAnalysis() { SetBit(kAODanalysis, kTRUE); };
+                void SetESDAnalysis() { SetBit(kAODanalysis, kFALSE); };
+                void SetEvent(AliESDEvent* const ESD){fESD1=ESD;};    // set ESD pointer
+                void SetEventAOD(AliAODEvent* const AOD){fAOD1=AOD;}; // set ESD pointer
+                void SetStack(AliStack* const stack){fStack=stack;};  // set stack pointer
+                void SetMCArray(TClonesArray* const mcarry){fMCArray=mcarry;} // set mcarray pointer
+                void SetUseMCPID(Bool_t usemcpid){fUseMCPID=usemcpid;};
 
-                void SetEvent(AliESDEvent* const ESD){fESD1=ESD;}; // set ESD pointer
-                void SetStack(AliStack* const stack){fStack=stack;} // set stack pointer
-                void InitAnaPair(); // initialize default parameters 
-                void AnaPair(AliESDtrack* ESDtrack1, AliESDtrack* ESDtrack2, Int_t index2); // do e-h analysis
-                void RunSECVTX(AliESDtrack *track); // run secondary vertexing algorithm
 
                 Int_t GetMCPID(AliESDtrack *track); // return MC pid
-                Int_t PairOrigin(AliESDtrack* track1, AliESDtrack* track2); // return pair origin as a pdg code
-                Int_t PairCode(AliESDtrack* track1,AliESDtrack* track2); // return corresponding pair code to pdg code
+                           Int_t GetMCPDG(AliVTrack *track);   // return MC pid
+                Int_t GetPairOriginESD(AliESDtrack* track1, AliESDtrack* track2); // return pair origin as a pdg code
+                Int_t GetPairOriginAOD(AliAODTrack* track1, AliAODTrack* track2); // return pair origin as a pdg code
+                Int_t GetPairCode(const AliVTrack* const track1, const AliVTrack* const track2); // return corresponding pair code to pdg code
                 Int_t GetElectronSource(Int_t mclabel); // return origin of the electron
+                Int_t GetPDG(AliVTrack *track);     // return pdg 
+                           void GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob); //return esd pid likelihood
+                void GetPrimaryCondition();
 
-                Bool_t SingleTrackCut(AliESDtrack* track1) const; // single track cut
+                TClonesArray *HFEpairs();
+                TClonesArray *HFEsecvtxs();
 
-        protected:
+                void AddHFEpairToArray(const AliHFEpairs* const pair);
+                void AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx);
 
-                void Init();
-                void FindSECVTXCandid4Tracks(AliESDtrack* track1); // find secondary vertex for 4 tracks
-                void FindSECVTXCandid3Tracks(AliESDtrack* track1); // find secondary vertex for 3 tracks        
-                void FindSECVTXCandid2Tracks(AliESDtrack* track1); // find secondary vertex for 2 tracks
-                void CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // calculated distinctive variables
-                void ApplyPairTagCut(); // apply strict tagging condition for 1 pair secondary vertex
-                void ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2); // apply tagging condition for 3 track secondary vertex
-                void ResetTagVar(); // reset tagging variables
+                void InitHFEpairs();
+                void InitHFEsecvtxs();
 
-                Bool_t ApplyPairTagCut(Int_t id); // apply strict tagging condition for 1 pair secondary vertex
-                Bool_t ApplyTagCut(Double_t chi2); // apply tagging condition
+                void DeleteHFEpairs();
+                void DeleteHFEsecvtxs();
+
+                Bool_t SingleTrackCut(AliESDtrack* track1) const; // single track cut
+                void PairAnalysis(AliVTrack* ESDtrack1, AliVTrack* ESDtrack2, Int_t index2); // do e-h analysis
+                void RunSECVTX(AliVTrack *track); // run secondary vertexing algorithm
 
+                void MakeContainer(); // make containers
 
-                Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4); // get secondary vertex chi2 for 4 tracks
-                Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // get secondary vertex chi2 for 3 tracks
-                Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2); // get secondary vertex chi2 for 2 tracks
+        protected:
+                void Init();
+                void FindSECVTXCandid(AliVTrack *track);
+                void CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3); // calculated distinctive variables
 
 
         private:
+                enum{
+                  kHasMCData = BIT(15),     // bitset for mc data usage
+                  kAODanalysis = BIT(16)    // bitset for aod analysis
+                };
+                enum {kAll, kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse}; // electron origin 
+                enum {kCharm=4, kBeauty=5}; // quark flavor
 
                 AliESDEvent* fESD1; // ESD pointer             
-                AliStack* fStack; // stack pointer              
+                AliAODEvent* fAOD1; // AOD pointer             
+                AliStack* fStack;   // stack pointer              
 
-                Int_t fParentSelect[2][7]; // heavy hadron species
-                Int_t fNparents; // number of heavy hadrons to be considered
+                Bool_t fUseMCPID;   // if use MC pid 
 
                 TString fkSourceLabel[10]; // electron source label
 
-                enum {kAll, kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse};
-                enum {kCharm=4, kBeauty=5};
-
-                struct AliHistspair{
-                        TH2F *fInvMass; // histogram to store pair invariant mass
-                        TH2F *fInvMassCut1; // histogram to store pair invariant mass after cut1
-                        TH2F *fInvMassCut2; // histogram to store pair invariant mass after cut2
-                        TH1F *fKFChi2; // histogram to store pair vertex chi2
-                       TH1F *fOpenAngle; // histogram to store pair opening angle
-                        TH1F *fCosOpenAngle; // histogram to store cosine of the pair opening angle
-                        TH2F *fSignedLxy; // histogram to store signed Lxy
-                        TH1F *fKFIP; // histogram to store pair impact parameter
-                        TH1F *fIPMax; // histogram to store larger impact parameter among pair tracks
-
-                       AliHistspair()
-                       : fInvMass()
-                       , fInvMassCut1()
-                       , fInvMassCut2()
-                       , fKFChi2()
-                       , fOpenAngle()
-                       , fCosOpenAngle()
-                       , fSignedLxy()
-                       , fKFIP()
-                       , fIPMax()
-                       {
-                         // define constructor
-                       }
-                       AliHistspair(const AliHistspair & p)
-                       : fInvMass(p.fInvMass)
-                       , fInvMassCut1(p.fInvMassCut1)
-                       , fInvMassCut2(p.fInvMassCut2)
-                       , fKFChi2(p.fKFChi2)
-                       , fOpenAngle(p.fOpenAngle)
-                       , fCosOpenAngle(p.fCosOpenAngle)
-                       , fSignedLxy(p.fSignedLxy)
-                       , fKFIP(p.fKFIP)
-                       , fIPMax(p.fIPMax)
-                       {
-                         // copy constructor
-                       }
-                       AliHistspair &operator=(const AliHistspair &)
-                       {
-                         // assignment operator, not yet implemented
-                         return *this;
-                       }
-                };
-
-                struct AliHiststag{
-                        TH1F *fPtBeautyElec; // histogram for electrons of single track cut passed 
-                        TH1F *fPtTaggedElec; // histogram for total b tagged electrons
-                        TH1F *fPtRightTaggedElec; // histogram for right b tagged electrons
-                        TH1F *fPtWrongTaggedElec; // histogram for wrong b tagged electrons
-                        TH1F *fInvmassBeautyElecSecVtx;  // histogram for right-tagged b invariant mass
-                        TH1F *fInvmassNotBeautyElecSecVtx; // histogram for mis-tagged b invariant mass
-                        TH1F *fSignedLxyBeautyElecSecVtx; // histogram for right-tagged b signed Lxy
-                        TH1F *fSignedLxyNotBeautyElecSecVtx; // histogram for mis-tagged b signed Lxy
-                        TH1F *fDistTwoVtxBeautyElecSecVtx; // histogram for right-tagged b two vertex distance
-                        TH1F *fDistTwoVtxNotBeautyElecSecVtx; // histogram for mis-tagged b two vertex distance
-                        TH1F *fCosPhiBeautyElecSecVtx; // histogram for right-tagged b cos of opening angle
-                        TH1F *fCosPhiNotBeautyElecSecVtx; // histogram for mis-tagged b cos of opening angle
-                        TH1F *fChi2BeautyElecSecVtx; // histogram for right-tagged b chi2 of secondary vertex 
-                        TH1F *fChi2NotBeautyElecSecVtx; // histogram for mis-tagged b chi2 of secondary vertex
-                        TH1F *fInvmassBeautyElec2trkVtx; // histogram for right-tagged b invariant mass of two track vertex
-                        TH1F *fInvmassNotBeautyElec2trkVtx; // histogram for mis-tagged b invariant mass of two track vertex
-                        TH1F *fSignedLxyBeautyElec2trkVtx; // histogram for right-tagged b signed Lxy of two track vertex
-                        TH1F *fSignedLxyNotBeautyElec2trkVtx; // histogram for mis-tagged b signed Lxy of two track vertex
-                        TH1F *fChi2BeautyElec2trkVtx; // histogram for right-tagged b chi2 of two track vertex
-                        TH1F *fChi2NotBeautyElec2trkVtx; // histogram for mis-tagged b chi2 of two track vertex
-
-                       AliHiststag()
-                       : fPtBeautyElec()
-                       , fPtTaggedElec()
-                       , fPtRightTaggedElec()
-                       , fPtWrongTaggedElec()
-                       , fInvmassBeautyElecSecVtx()
-                       , fInvmassNotBeautyElecSecVtx()
-                       , fSignedLxyBeautyElecSecVtx()
-                       , fSignedLxyNotBeautyElecSecVtx()
-                       , fDistTwoVtxBeautyElecSecVtx()
-                       , fDistTwoVtxNotBeautyElecSecVtx()
-                       , fCosPhiBeautyElecSecVtx()
-                       , fCosPhiNotBeautyElecSecVtx()
-                       , fChi2BeautyElecSecVtx()
-                       , fChi2NotBeautyElecSecVtx()
-                       , fInvmassBeautyElec2trkVtx()
-                       , fInvmassNotBeautyElec2trkVtx()
-                       , fSignedLxyBeautyElec2trkVtx()
-                       , fSignedLxyNotBeautyElec2trkVtx()
-                       , fChi2BeautyElec2trkVtx()
-                       , fChi2NotBeautyElec2trkVtx()
-                       {
-                         // define constructor
-                       }
-                       AliHiststag(const AliHiststag & p)
-                       : fPtBeautyElec(p.fPtBeautyElec)
-                       , fPtTaggedElec(p.fPtTaggedElec)
-                       , fPtRightTaggedElec(p.fPtRightTaggedElec)
-                       , fPtWrongTaggedElec(p.fPtWrongTaggedElec)
-                       , fInvmassBeautyElecSecVtx(p.fInvmassBeautyElecSecVtx)
-                       , fInvmassNotBeautyElecSecVtx(p.fInvmassNotBeautyElecSecVtx)
-                       , fSignedLxyBeautyElecSecVtx(p.fSignedLxyBeautyElecSecVtx)
-                       , fSignedLxyNotBeautyElecSecVtx(p.fSignedLxyNotBeautyElecSecVtx)
-                       , fDistTwoVtxBeautyElecSecVtx(p.fDistTwoVtxBeautyElecSecVtx)
-                       , fDistTwoVtxNotBeautyElecSecVtx(p.fDistTwoVtxNotBeautyElecSecVtx)
-                       , fCosPhiBeautyElecSecVtx(p.fCosPhiBeautyElecSecVtx)
-                       , fCosPhiNotBeautyElecSecVtx(p.fCosPhiNotBeautyElecSecVtx)
-                       , fChi2BeautyElecSecVtx(p.fChi2BeautyElecSecVtx)
-                       , fChi2NotBeautyElecSecVtx(p.fChi2NotBeautyElecSecVtx)
-                       , fInvmassBeautyElec2trkVtx(p.fInvmassBeautyElec2trkVtx)
-                       , fInvmassNotBeautyElec2trkVtx(p.fInvmassNotBeautyElec2trkVtx)
-                       , fSignedLxyBeautyElec2trkVtx(p.fSignedLxyBeautyElec2trkVtx)
-                       , fSignedLxyNotBeautyElec2trkVtx(p.fSignedLxyNotBeautyElec2trkVtx)
-                       , fChi2BeautyElec2trkVtx(p.fChi2BeautyElec2trkVtx)
-                       , fChi2NotBeautyElec2trkVtx(p.fChi2NotBeautyElec2trkVtx)
-                       {
-                         // copy constructor
-                       }
-                       AliHiststag &operator=(const AliHiststag &)
-                       {
-                         // assignment operator, not yet implemented
-                         return *this;
-                       }
-                };
-
-                AliHistspair fHistPair[10]; // struct of above given histogram for different electron sources
-                AliHiststag fHistTagged; // struct of above given histogram
+                Int_t fNparents;           // number of heavy hadrons to be considered
+                Int_t fParentSelect[2][7]; // heavy hadron species
 
-                Int_t fPairTagged; // number of tagged e-h pairs
-                Int_t fpairedTrackID[20]; // paird hadron track track 
-                Double_t fpairedChi2[20]; // pair chi2
-                Double_t fpairedInvMass[20]; // pair invariant mass
-                Double_t fpairedSignedLxy[20]; // pair signed Lxy
+                Int_t fNoOfHFEpairs;       // number of e-h pairs  
+                Int_t fNoOfHFEsecvtxs;     // number of secondary vertexes
 
-                Int_t fid[4][3]; // index to store sorting result
-                Int_t fia[4][3][2]; // index to store sorting result
+                TClonesArray *fHFEpairs;   //! Array of pair 
+                TClonesArray *fHFEsecvtxs; //! Array of secondary vertexes 
+                TClonesArray *fMCArray;    //! mc array pointer
 
-                Double_t fdistTwoSecVtx; // distance between two pair vertex
-                Double_t fcosPhi; // cos of opening angle of two pair vertex
-                Double_t fsignedLxy; // signed Lxy of secondary vertex
-                Double_t finvmass; // invariant mass of secondary vertex
-                Double_t finvmassSigma; // invariant mass sigma of secondary vertex
+                           Double_t fPVx;          // primary vertex copy x 
+                           Double_t fPVy;          // primary vertex copy y
+                Double_t fCosPhi;       // cos of opening angle of two pair vertex
+                Double_t fSignedLxy;    // signed Lxy of secondary vertex
+                Double_t fKFchi2;       // chi2 of secondary vertex
+                Double_t fInvmass;      // invariant mass of secondary vertex
+                Double_t fInvmassSigma; // invariant mass sigma of secondary vertex
+                Double_t fKFip;         // impact parameter of secondary vertex track
 
-                Bool_t fBTagged; // if b tagged, set true
-                Bool_t fBElec; // if set true for b electron, set true
+                THnSparseF *fPairQA;    // qa histos for pair analysis 
+                THnSparseF *fSecvtxQA;  // qa histos for secvtx
+                TList *fSecVtxList;     // list for secondary vertexing outputs
 
-        ClassDef(AliHFEsecVtx,0);
+      ClassDef(AliHFEsecVtx,0);
 };
 
 #endif