Fix coding rules violations
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Jun 2009 12:40:35 +0000 (12:40 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Jun 2009 12:40:35 +0000 (12:40 +0000)
13 files changed:
PWG3/hfe/AliAnalysisTaskHFE.cxx
PWG3/hfe/AliAnalysisTaskHFE.h
PWG3/hfe/AliHFEcuts.cxx
PWG3/hfe/AliHFEcuts.h
PWG3/hfe/AliHFEextraCuts.cxx
PWG3/hfe/AliHFEextraCuts.h
PWG3/hfe/AliHFEmcQA.cxx
PWG3/hfe/AliHFEmcQA.h
PWG3/hfe/AliHFEpidTRD.cxx
PWG3/hfe/AliHFEpidTRD.h
PWG3/hfe/AliHFEpriVtx.h
PWG3/hfe/AliHFEsecVtx.cxx
PWG3/hfe/AliHFEsecVtx.h

index d2d0227..c8cd573 100644 (file)
@@ -30,7 +30,6 @@
 #include <TH1F.h>
 #include <TH1I.h>
 #include <TH2F.h>
-#include <THnSparse.h>
 #include <TIterator.h>
 #include <TList.h>
 #include <TLegend.h>
@@ -67,19 +66,23 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fMC(0x0)
   , fCFM(0x0)
   , fCorrelation(0x0)
+  , fFakeElectrons(0x0)
   , fPID(0x0)
   , fCuts(0x0)
   , fSecVtx(0x0)
-  , fAnalysisMCQA(0x0)
+  , fMCQA(0x0)
   , fNEvents(0x0)
   , fQA(0x0)
   , fOutput(0x0)
   , fHistMCQA(0x0)
   , fHistSECVTX(0x0)
 {
-       DefineInput(0, TChain::Class());
-       DefineOutput(0, TH1I::Class());
-       DefineOutput(1, TList::Class());
+  //
+  // Default constructor
+  //
+  DefineInput(0, TChain::Class());
+  DefineOutput(0, TH1I::Class());
+  DefineOutput(1, TList::Class());
   DefineOutput(2, TList::Class());
 
   // Initialize cuts
@@ -88,11 +91,54 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
 }
 
 //____________________________________________________________
+AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
+  AliAnalysisTask("PID efficiency Analysis COPY", "")
+  , fESD(0x0)
+  , fMC(0x0)
+  , fCFM(0x0)
+  , fCorrelation(0x0)
+  , fFakeElectrons(0x0)
+  , fPID(0x0)
+  , fCuts(0x0)
+  , fSecVtx(0x0)
+  , fMCQA(0x0)
+  , fNEvents(0x0)
+  , fQA(0x0)
+  , fOutput(0x0)
+  , fHistMCQA(0x0)
+  , fHistSECVTX(0x0)
+{
+  //
+  // Copy Constructor
+  //
+  DefineInput(0, TChain::Class());
+  DefineOutput(0, TH1I::Class());
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
+
+  ref.Copy(*this);
+}
+
+//____________________________________________________________
+AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
+  //
+  // Assignment operator
+  //
+  if(this != &ref)
+    ref.Copy(*this);
+
+  return *this;
+}
+
+//____________________________________________________________
 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
-       if(fESD) delete fESD;
-       if(fMC) delete fMC;
-       if(fPID) delete fPID;
-       if(fQA){
+  //
+  // Destructor
+  //
+  if(fESD) delete fESD;
+  if(fMC) delete fMC;
+  if(fPID) delete fPID;
+  if(fQA){
     fQA->Clear();
     delete fQA;
   }
@@ -110,49 +156,65 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
   }
   if(fCuts) delete fCuts;
   if(fSecVtx) delete fSecVtx;
-  if(fAnalysisMCQA) delete fAnalysisMCQA;
-       if(fNEvents) delete fNEvents;
+  if(fMCQA) delete fMCQA;
+  if(fNEvents) delete fNEvents;
   if(fCorrelation) delete fCorrelation;
   if(fFakeElectrons) delete fFakeElectrons;
 }
 
 //____________________________________________________________
+void AliAnalysisTaskHFE::Copy(TObject &o) const {
+  //
+  // Make a copy of the task objecy
+  //
+  
+  AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
+
+  // copy objects
+  if(fPID) target.fPID = new AliHFEpid(*fPID);
+  if(fCuts) target.fCuts = new AliHFEcuts(*fCuts);
+
+  if(fESD || fMC) target.ConnectInputData(0x0);
+  if(fOutput || fQA) target.CreateOutputObjects();
+}
+
+//____________________________________________________________
 void AliAnalysisTaskHFE::ConnectInputData(Option_t *){
-       TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
-       if(!esdchain){
-               AliError("ESD chain empty");
-               return;
-       } else {
-               esdchain->SetBranchStatus("Tracks", 1);
-       }
-       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
-       if(!esdH){      
-               AliError("No ESD input handler");
-               return;
-       } else {
-               fESD = esdH->GetEvent();
-       }
-       AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-       if(!mcH){       
-               AliError("No MC truth handler");
-               return;
-       } else {
-               fMC = mcH->MCEvent();
-       }
+  TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
+  if(!esdchain){
+    AliError("ESD chain empty");
+    return;
+  } else {
+    esdchain->SetBranchStatus("Tracks", 1);
+  }
+  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+  if(!esdH){      
+    AliError("No ESD input handler");
+    return;
+  } else {
+    fESD = esdH->GetEvent();
+  }
+  AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+  if(!mcH){       
+    AliError("No MC truth handler");
+    return;
+  } else {
+    fMC = mcH->MCEvent();
+  }
 }
 
 //____________________________________________________________
 void AliAnalysisTaskHFE::CreateOutputObjects(){
-       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
-       // First Step: TRD alone
-       if(!fQA) fQA = new TList;
-       fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
-       fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
-       fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
-       fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
-       fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
-       fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
-       fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
+  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
+  // First Step: TRD alone
+  if(!fQA) fQA = new TList;
+  fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
+  fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
+  fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
+  fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
+  fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
+  fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
+  fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
 
   if(!fOutput) fOutput = new TList;
   // Initialize correction Framework and Cuts
@@ -174,183 +236,393 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){
 
   // Initialize PID
   if(IsQAOn()){
+    AliInfo("QA on for Cuts and PID");
     fPID->SetQAOn();
     fQA->Add(fPID->GetQAhistograms());
   }
   fPID->SetHasMCData(kTRUE);
   fPID->InitializePID("TRD");
 
-  // [mj] mcQA----------------------------------
+  // mcQA----------------------------------
   if (IsMCQAOn()) {
+    AliInfo("MC QA on");
+    if(!fMCQA) fMCQA = new AliHFEmcQA;
     if(!fHistMCQA) fHistMCQA = new TList();
+    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
+    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
+    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
+    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
+    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("mcqa") == 0) fHistMCQA->AddAt(obj_, counter_++);
+      }
+    }
     fQA->Add(fHistMCQA);
-    fAnalysisMCQA->CreateHistograms(AliHFEmcQA::fkCharm,"mcqa_");               // create histograms for charm
-    fAnalysisMCQA->CreateHistograms(AliHFEmcQA::fkBeauty,"mcqa_");              // create histograms for beauty
-    //Rossella
-//     fAnalysisMCQA->CreateHistosHadrons();                   // create histograms for hadrons
+  } 
 
-  }
-  // [mj] secvtx----------------------------
+  // secvtx----------------------------------
   if (IsSecVtxOn()) {
-    fOutput->Add(fHistSECVTX);
+    AliInfo("Secondary Vertex Analysis on");
+    fSecVtx = new AliHFEsecVtx;
+
     if(!fHistSECVTX) fHistSECVTX = new TList();
     fSecVtx->CreateHistograms("secvtx_");
-  }
-
-  TIter next_(gDirectory->GetList());
-  TObject *obj_;
-  int counter_ = 0;
-  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_++);
-      }
-      else if ((fpart->String()).CompareTo("mcqa") == 0){
-        fHistMCQA->AddAt(obj_, counter__++);
+    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__++);
       }
     }
-  }
-  //--------------------------------------- 
+  fOutput->Add(fHistSECVTX);
+  } 
 }
 
 //____________________________________________________________
 void AliAnalysisTaskHFE::Exec(Option_t *){
-       //
-       // Run the analysis
-       //
-       if(!fESD){
-               AliError("No ESD Event");
-               return;
-       }
-       if(!fMC){
-               AliError("No MC Event");
-               return;
-       }
-       fCFM->SetEventInfo(fMC);
+  //
+  // Run the analysis
+  // 
+  if(!fESD){
+    AliError("No ESD Event");
+    return;
+  }
+  if(!fMC){
+    AliError("No MC Event");
+    return;
+  }
+  fCFM->SetEventInfo(fMC);
   fPID->SetMCEvent(fMC);
 
-       //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
+  //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
+  Double_t container[6];
 
-       Double_t container[6];
+  // Loop over the Monte Carlo tracks to see whether we have overlooked any track
+  AliMCParticle *mctrack = 0x0;
+  Int_t nElectrons = 0;
 
-       // Loop over the Monte Carlo tracks to see whether we have overlooked any track
-       AliMCParticle *mctrack = 0x0;
-       Int_t nElectrons = 0;
+  if (IsSecVtxOn()) { 
+    fSecVtx->SetEvent(fESD);
+    fSecVtx->SetStack(fMC->Stack());
+  }
 
-  // [mj] run MC QA ------------------------------------------------
+  // run mc QA ------------------------------------------------
   if (IsMCQAOn()) {
+    AliDebug(2, "Running MC QA");
 
-    fAnalysisMCQA->SetStack(fMC->Stack());
-    fAnalysisMCQA->Init();
+    fMCQA->SetStack(fMC->Stack());
+    fMCQA->Init();
 
     Int_t nPrims = fMC->Stack()->GetNprimary();
     Int_t nMCTracks = fMC->Stack()->GetNtrack();
 
     // loop over primary particles for quark and heavy hadrons
     for (Int_t igen = 0; igen < nPrims; igen++){
-      fAnalysisMCQA->GetQuarkKine(igen, AliHFEmcQA::fkCharm);
-      fAnalysisMCQA->GetQuarkKine(igen, AliHFEmcQA::fkBeauty);
+      fMCQA->GetQuarkKine(igen, AliHFEmcQA::kCharm);
+      fMCQA->GetQuarkKine(igen, AliHFEmcQA::kBeauty);
+      fMCQA->GetHadronKine(igen, AliHFEmcQA::kCharm);
+      fMCQA->GetHadronKine(igen, AliHFEmcQA::kBeauty);
     }
-    fAnalysisMCQA->EndOfEventAna(AliHFEmcQA::fkCharm);
-    fAnalysisMCQA->EndOfEventAna(AliHFEmcQA::fkBeauty);
+    fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
+    fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 
     // loop over all tracks for decayed electrons
     for (Int_t igen = 0; igen < nMCTracks; igen++){
-      fAnalysisMCQA->GetDecayedKine(igen, AliHFEmcQA::fkCharm, AliHFEmcQA::fkElectron);
-      fAnalysisMCQA->GetDecayedKine(igen, AliHFEmcQA::fkBeauty, AliHFEmcQA::fkElectron);
+      fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0);
+      fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0);
+      fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
+      fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
     }
 
   } // end of MC QA loop
   // -----------------------------------------------------------------
 
-       for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
-               mctrack = fMC->GetTrack(imc);
-               container[0] = mctrack->Pt();
-               container[1] = mctrack->Eta();
+  //
+  // Loop MC
+  //
+  for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
+    mctrack = fMC->GetTrack(imc);
+
+    container[0] = mctrack->Pt();
+    container[1] = mctrack->Eta();
     container[2] = mctrack->Phi();
 
-               if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
-               fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
-               (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++;
-       }
-       (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
-
-       // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
-       AliESDtrack *track = 0x0, *htrack = 0x0;
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
+    (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++;
+  }
+  (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
+
+  // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
+        
+
+  //
+  // Loop ESD
+  //
+
+  Int_t nbtrackrec = fESD->GetNumberOfTracks();
+  AliESDtrack *track = 0x0, *htrack = 0x0;
   Int_t pid = 0;
-       for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
-               track = fESD->GetTrack(itrack);
-               container[0] = track->Pt();
-               container[1] = track->Eta();
+  // For double counted tracks
+  Int_t *indexRecKineTPC = new Int_t[nbtrackrec];
+  Int_t *indexRecKineITS = new Int_t[nbtrackrec];
+  Int_t *indexRecPrim = new Int_t[nbtrackrec];
+  Int_t *indexHFEcutsITS = new Int_t[nbtrackrec];
+  Int_t *indexHFEcutsTPC = new Int_t[nbtrackrec];
+  Int_t *indexHFEcutsTRD = new Int_t[nbtrackrec];
+  Int_t *indexPid = new Int_t[nbtrackrec];
+
+  for(Int_t nb = 0; nb < nbtrackrec; nb++){
+    indexRecKineTPC[nb] = -1;
+    indexRecKineITS[nb] = -1;
+    indexRecPrim[nb]    = -1;
+    indexHFEcutsITS[nb] = -1;
+    indexHFEcutsTPC[nb] = -1;
+    indexHFEcutsTRD[nb] = -1;
+    indexPid[nb]        = -1;
+  } 
+
+  Bool_t alreadyseen = kFALSE;
+
+  for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
+    
+    track = fESD->GetTrack(itrack);
+          
+    container[0] = track->Pt();
+    container[1] = track->Eta();
     container[2] = track->Phi();
-               if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKine, track)) continue;
-               fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKine);
-
-               // Check TRD criterions (outside the correction framework)
-               if(track->GetTRDncls()){
-                       (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
-                       (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());
-               }
-               if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
+
+    Bool_t signal = kTRUE;
+          
+    // RecKine: TPC cuts  
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineTPC, track)) continue;
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineTPC);
+    
+
+    // Check if it is signal electrons
+    if(!(mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())))) continue;
+    
+    container[3] = mctrack->Pt();
+    container[4] = mctrack->Eta();
+    container[5] = mctrack->Phi();
+    
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
+    
+    if(signal) {
+      alreadyseen = kFALSE;
+      for(Int_t nb = 0; nb < itrack; nb++){
+        if(indexRecKineTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
+      }
+      indexRecKineTPC[itrack] = TMath::Abs(track->GetLabel());
+      
+      fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineTPC + AliHFEcuts::kNcutESDSteps + 1));
+      fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 2*(AliHFEcuts::kNcutESDSteps + 1)));
+      if(alreadyseen) {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+      else {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+    }
+
+    
+    // RecKine: ITS cuts
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITS, track)) continue;
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineITS);
+    if(signal) {
+      alreadyseen = kFALSE;
+      for(Int_t nb = 0; nb < itrack; nb++){
+        if(indexRecKineITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
+      }
+      indexRecKineITS[itrack] = TMath::Abs(track->GetLabel());
+      fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineITS + AliHFEcuts::kNcutESDSteps + 1));
+      fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 2*(AliHFEcuts::kNcutESDSteps + 1)));
+      if(alreadyseen) {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+      else {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+    }
+
+
+    // Check TRD criterions (outside the correction framework)
+    if(track->GetTRDncls()){
+      (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
+      (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());
+    }
+
+    
+    // RecPrim
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
-    if(fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcuts, track)) continue;
-    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts);
+    if(signal) {
+      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1);
+      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1));
+      alreadyseen = kFALSE;
+      for(Int_t nb = 0; nb < itrack; nb++){
+        if(indexRecPrim[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
+      }
+      indexRecPrim[itrack] = TMath::Abs(track->GetLabel());
+      if(alreadyseen) {
+        fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 4*(AliHFEcuts::kNcutESDSteps + 1));
+      }
+      else {
+        fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 3*(AliHFEcuts::kNcutESDSteps + 1));
+      }
+    }
+
+
+    // HFEcuts: ITS layers cuts
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS);
+    if(signal) {
+      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1);
+      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1));
+      alreadyseen = kFALSE;
+      for(Int_t nb = 0; nb < itrack; nb++){
+        if(indexHFEcutsITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
+      }
+      indexHFEcutsITS[itrack] = TMath::Abs(track->GetLabel());
+      if(alreadyseen) {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+      else {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+    }
+
+    // HFEcuts: TPC ratio clusters
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC);
+    if(signal) {
+      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC + AliHFEcuts::kNcutESDSteps + 1);
+      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTPC + 2*(AliHFEcuts::kNcutESDSteps + 1));
+      alreadyseen = kFALSE;
+      for(Int_t nb = 0; nb < itrack; nb++){
+        if(indexHFEcutsTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
+      }
+      indexHFEcutsTPC[itrack] = TMath::Abs(track->GetLabel());    
+      if(alreadyseen) {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+      else {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+    }
+    
+    // HFEcuts: Nb of tracklets TRD
+    if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD);
+    if(signal) {
+      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 1);
+      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1));
+      alreadyseen = kFALSE;
+      for(Int_t nb = 0; nb < itrack; nb++){
+        if(indexHFEcutsTRD[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
+      }
+      indexHFEcutsTRD[itrack] = TMath::Abs(track->GetLabel());
+      if(alreadyseen) {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 4*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+      else {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 3*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+    }
+
+
     // track accepted, do PID
-               if(!fPID->IsSelected(track)) continue;
-    // Track selected: distinguish between true and fake
-    if(!(mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())))) continue;
-    if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
-         fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts + 1);
+    if(!fPID->IsSelected(track)) continue;
+
+
+    // Fill Containers
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 1);
+    if(signal) {
+      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 2);
+      fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1)+1);
+      alreadyseen = kFALSE;
+      for(Int_t nb = 0; nb < itrack; nb++){
+        if(indexPid[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
+      }
+      indexPid[itrack] = TMath::Abs(track->GetLabel());
+      if(alreadyseen) {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
+      else {
+        fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 3*(AliHFEcuts::kNcutESDSteps + 1)));
+      }
       // dimensions 3&4&5 : pt,eta,phi (MC)
-      container[3] = mctrack->Pt();
-      container[4] = mctrack->Eta();
-      container[5] = mctrack->Phi();
       fCorrelation->Fill(container);
+    }
+
+
+    // Track selected: distinguish between true and fake
+    if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
       // pair analysis [mj]
       if (IsSecVtxOn()) {
+        AliDebug(2, "Running Secondary Vertex Analysis");
+        fSecVtx->InitAnaPair();
         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, fSecVtx->PairCode(track,htrack), jtrack);
+          fSecVtx->AnaPair(track, htrack, jtrack);
         }
-      }
+        // 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
       switch(pid){
-        case 13:    container[3] = AliPID::kMuon;
-        case 211:   container[3] = AliPID::kPion;
-        case 321:   container[3] = AliPID::kKaon;
-        case 2212:  container[3] = AliPID::kProton;
+      case 13:    container[3] = AliPID::kMuon;
+      case 211:   container[3] = AliPID::kPion;
+      case 321:   container[3] = AliPID::kKaon;
+      case 2212:  container[3] = AliPID::kProton;
       }
       fFakeElectrons->Fill(container);
     }
-       }
-       fNEvents->Fill(1);
-
-       // Done!!!
-       PostData(0, fNEvents);
-       PostData(1, fOutput);
-       PostData(2, fQA);
+  }
+  fNEvents->Fill(1);
+  
+  // release memory
+  delete[] indexRecKineTPC;
+  delete[] indexRecKineITS;
+  delete[] indexRecPrim;
+  delete[] indexHFEcutsITS;
+  delete[] indexHFEcutsTPC;
+  delete[] indexHFEcutsTRD;
+  delete[] indexPid;
+  
+  // Done!!!
+  PostData(0, fNEvents);
+  PostData(1, fOutput);
+  PostData(2, fQA);
 }
 
 //____________________________________________________________
 void AliAnalysisTaskHFE::Terminate(Option_t *){
-       //
-       // Terminate not implemented at the moment
-       //
+  //
+  // Terminate not implemented at the moment
+  //
 }
 
 //____________________________________________________________
@@ -360,13 +632,13 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
   // link it
   //
   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
-  const Double_t kPtmin = 0., kPtmax = 10.;
+  const Double_t kPtmin = 0.1, kPtmax = 10.;
   const Double_t kEtamin = -0.9, kEtamax = 0.9;
   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 
   //arrays for the number of bins in each dimension
   Int_t iBin[kNvar];
-  iBin[0] = 20; //bins in pt
+  iBin[0] = 40; //bins in pt
   iBin[1] =  8; //bins in eta 
   iBin[2] = 18; // bins in phi
 
@@ -376,14 +648,16 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 
   //values for bin lower bounds
-  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/iBin[0]*(Double_t)i; 
+  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);  
   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 
   //one "container" for MC
-  AliCFContainer* container = new AliCFContainer("container","container for tracks", AliHFEcuts::kNcutSteps + 1, kNvar, iBin);
+  AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutSteps + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)), kNvar, iBin);
+
   //setting the bin limits
-  for(Int_t ivar = 0; ivar < kNvar; ivar++)  container -> SetBinLimits(0, binEdges[ivar]);
+  for(Int_t ivar = 0; ivar < kNvar; ivar++)
+    container -> SetBinLimits(ivar, binEdges[ivar]);
   fCFM->SetParticleContainer(container);
 
   //create correlation matrix for unfolding
@@ -408,3 +682,4 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
   for(Int_t idim = 0; idim < kNvar; idim++)
     fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
 }
+
index f158dbb..56af78a 100644 (file)
 #ifndef ALIANALYSISTASKHFE_H
 #define ALIANALYSISTASKHFE_H
 
+#ifndef ALIANALYSISTASK_H
 #include "AliAnalysisTask.h"
+#endif
+
+#ifndef ROOT_THnSparse
+#include <THnSparse.h>
+#endif
 
 class AliHFEpid;
 class AliHFEcuts;
@@ -26,49 +32,57 @@ class AliESDEvent;
 class AliESDtrackCuts;
 class AliMCEvent;
 class TH1I; 
-class THnSparse;
 class TList;
 
 class AliAnalysisTaskHFE : public AliAnalysisTask{
   enum{
-    kIsQAOn = BIT(14),
-    kIsMCQAOn = BIT(15),
-    kIsSecVtxOn = BIT(16)
+    kIsQAOn = BIT(18),
+    kIsMCQAOn = BIT(19),
+    kIsSecVtxOn = BIT(20),
+    kIsPriVtxOn = BIT(21)
   };
-       public:
-               AliAnalysisTaskHFE();
-               ~AliAnalysisTaskHFE();
+  public:
+    AliAnalysisTaskHFE();
+    AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref);
+    AliAnalysisTaskHFE& operator=(const AliAnalysisTaskHFE &ref);
+    ~AliAnalysisTaskHFE();
 
-               virtual void ConnectInputData(Option_t *);
-               virtual void CreateOutputObjects();
-               virtual void Exec(Option_t *);
-               virtual void Terminate(Option_t *);
+    virtual void ConnectInputData(Option_t *);
+    virtual void CreateOutputObjects();
+    virtual void Exec(Option_t *);
+    virtual void Terminate(Option_t *);
 
     Bool_t IsQAOn() const     { return TestBit(kIsQAOn); };
-    Bool_t IsMCQAOn() const   { return TestBit(kIsMCQAOn); }
+    Bool_t IsMCQAOn() const   { return TestBit(kIsMCQAOn); };
     Bool_t IsSecVtxOn() const { return TestBit(kIsSecVtxOn); };
+    Bool_t IsPriVtxOn() const { return TestBit(kIsPriVtxOn); };
     void SetQAOn()            { SetBit(kIsQAOn, kTRUE); };
     void SetMCQAOn()          { SetBit(kIsMCQAOn, kTRUE); };
+    void SetPriVtxOn()        { SetBit(kIsPriVtxOn, kTRUE); };
     void SetSecVtxOn()        { SetBit(kIsSecVtxOn, kTRUE); };
+  protected:
+    void Copy(TObject &o) const;
 
-       private:
+  private:
     void MakeParticleContainer();
 
-               AliESDEvent *fESD;                              //! The ESD Event
-               AliMCEvent *fMC;                                  //! The MC Event
-               AliCFManager *fCFM;                             //! Correction Framework Manager
+    AliESDEvent *fESD;                    //! The ESD Event
+    AliMCEvent *fMC;                      //! The MC Event
+    AliCFManager *fCFM;                   //! Correction Framework Manager
     THnSparseF *fCorrelation;             //! response matrix for unfolding  
     THnSparseF *fFakeElectrons;           //! Contamination from Fake Electrons
-               AliHFEpid *fPID;                      //! PID
+    AliHFEpid *fPID;                      //! PID
     AliHFEcuts *fCuts;                    //! Cut Collection
     AliHFEsecVtx *fSecVtx;                //! Secondary Vertex Analysis
-    AliHFEmcQA *fAnalysisMCQA;            //! MC QA
-               TH1I *fNEvents;                                   //! counter for the number of Events
-               TList *fQA;                                           //! QA histos for the cuts
+    AliHFEmcQA *fMCQA;                    //! MC QA
+    TH1I *fNEvents;                       //! counter for the number of Events
+    TList *fQA;                           //! QA histos for the cuts
     TList *fOutput;                       //! Container for Task Output
     TList *fHistMCQA;                     //! Output container for MC QA histograms 
     TList *fHistSECVTX;                   //! Output container for sec. vertexing results
 
-       ClassDef(AliAnalysisTaskHFE, 1)    // The electron Analysis Task
+    ClassDef(AliAnalysisTaskHFE, 1)       // The electron Analysis Task
 };
 #endif
+
index 5dfa40a..608c05b 100644 (file)
@@ -18,6 +18,7 @@
  * ALICE Heavy Flavour Electron Group                                   *
  *                                                                      *
  * Authors:                                                             *
+ *   Raphaelle Bailhache <R.Bailhache@gsi.de>                           *
  *   Markus Fasel <M.Fasel@gsi.de>                                      *
  *   Markus Heide <mheide@uni-muenster.de>                              *
  *   Matus Kalisky <m.kalisky@uni-muenster.de>                          *
@@ -42,8 +43,8 @@
 
 ClassImp(AliHFEcuts)
 
-const Int_t AliHFEcuts::kNcutSteps = 5;
-
+const Int_t AliHFEcuts::kNcutSteps = 8;
+const Int_t AliHFEcuts::kNcutESDSteps = 6;
 //__________________________________________________________________
 AliHFEcuts::AliHFEcuts():
   fRequirements(0),
@@ -54,7 +55,8 @@ AliHFEcuts::AliHFEcuts():
   fMinClusterRatioTPC(0.),
   fSigmaToVtx(0.),
   fHistQA(0x0),
-  fCutList(0x0)
+  fCutList(0x0),
+  fDebugLevel(0)
 {
   //
   // Default Constructor
@@ -77,7 +79,8 @@ AliHFEcuts::AliHFEcuts(const AliHFEcuts &c):
   fMinClusterRatioTPC(c.fMinClusterRatioTPC),
   fSigmaToVtx(c.fSigmaToVtx),
   fHistQA(0x0),
-  fCutList(0x0)
+  fCutList(0x0),
+  fDebugLevel(0)
 {
   //
   // Copy Constructor
@@ -105,19 +108,28 @@ AliHFEcuts::~AliHFEcuts(){
 
 //__________________________________________________________________
 void AliHFEcuts::Initialize(AliCFManager *cfm){
+
   // Call all the setters for the cuts
   SetParticleGenCutList();
   SetAcceptanceCutList();
-  SetRecKineCutList();
+  SetRecKineTPCCutList();
+  SetRecKineITSCutList();
   SetRecPrimaryCutList();
-  SetHFElectronCuts();
+  SetHFElectronITSCuts();
+  SetHFElectronTPCCuts();
+  SetHFElectronTRDCuts();
+
   
   // Connect the cuts
   cfm->SetParticleCutsList(kStepMCGenerated, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartGenCuts")));
   cfm->SetParticleCutsList(kStepMCInAcceptance, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartAccCuts")));
-  cfm->SetParticleCutsList(kStepRecKine, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartRecCuts")));
+  cfm->SetParticleCutsList(kStepRecKineTPC, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartRecKineTPCCuts")));
+  cfm->SetParticleCutsList(kStepRecKineITS, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartRecKineITSCuts")));
   cfm->SetParticleCutsList(kStepRecPrim, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartPrimCuts")));
-  cfm->SetParticleCutsList(kStepHFEcuts, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartHFECuts")));
+  cfm->SetParticleCutsList(kStepHFEcutsITS, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartHFECutsITS")));
+  cfm->SetParticleCutsList(kStepHFEcutsTPC, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartHFECutsTPC")));
+  cfm->SetParticleCutsList(kStepHFEcutsTRD, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartHFECutsTRD")));
+
 }
 
 //__________________________________________________________________
@@ -125,9 +137,13 @@ void AliHFEcuts::Initialize(){
   // Call all the setters for the cuts
   SetParticleGenCutList();
   SetAcceptanceCutList();
-  SetRecKineCutList();
+  SetRecKineTPCCutList();
+  SetRecKineITSCutList();
   SetRecPrimaryCutList();
-  SetHFElectronCuts();
+  SetHFElectronITSCuts();
+  SetHFElectronTPCCuts();
+  SetHFElectronTRDCuts();
+
 }
 
 //__________________________________________________________________
@@ -186,7 +202,7 @@ void AliHFEcuts::SetAcceptanceCutList(){
 }
 
 //__________________________________________________________________
-void AliHFEcuts::SetRecKineCutList(){
+void AliHFEcuts::SetRecKineTPCCutList(){
   //
   // Track Kinematics and Quality cuts (Based on the Standard cuts from PWG0)
   //
@@ -199,7 +215,6 @@ void AliHFEcuts::SetRecKineCutList(){
   // Min. Number of Clusters:
   //  TPC: 50
   // RefitRequired:
-  //  ITS
   //  TPC
   // Chi2 per TPC cluster: 3.5
   //
@@ -210,7 +225,7 @@ void AliHFEcuts::SetRecKineCutList(){
   AliCFTrackQualityCuts *trackQuality = new AliCFTrackQualityCuts("fCutsQualityRec","REC Track Quality Cuts");
   trackQuality->SetMinNClusterTPC(fMinClustersTPC);
   trackQuality->SetMaxChi2PerClusterTPC(fMaxChi2clusterTPC);
-  trackQuality->SetStatus(AliESDtrack::kTPCrefit & AliESDtrack::kITSrefit);
+  trackQuality->SetStatus(AliESDtrack::kTPCrefit);
   trackQuality->SetMaxCovDiagonalElements(2., 2., 0.5, 0.5, 2); 
 
   AliCFTrackKineCuts *kineCuts = new AliCFTrackKineCuts("fCutsKineRec", "REC Kine Cuts");
@@ -223,13 +238,34 @@ void AliHFEcuts::SetRecKineCutList(){
   }
   
   TObjArray *recCuts = new TObjArray;
-  recCuts->SetName("fPartRecCuts");
+  recCuts->SetName("fPartRecKineTPCCuts");
   recCuts->AddLast(trackQuality);
   recCuts->AddLast(kineCuts);
   fCutList->AddLast(recCuts);
 }
 
 //__________________________________________________________________
+void AliHFEcuts::SetRecKineITSCutList(){
+  //
+  // Track Kinematics and Quality cuts (Based on the Standard cuts from PWG0)
+  //
+  // RefitRequired:
+  //  ITS
+  //
+  AliCFTrackQualityCuts *trackQuality = new AliCFTrackQualityCuts("fCutsQualityRec","REC Track Quality Cuts");
+  trackQuality->SetStatus(AliESDtrack::kITSrefit);
+  
+  if(IsInDebugMode()){
+    trackQuality->SetQAOn(fHistQA);
+  }
+  
+  TObjArray *recCuts = new TObjArray;
+  recCuts->SetName("fPartRecKineITSCuts");
+  recCuts->AddLast(trackQuality);
+  fCutList->AddLast(recCuts);
+}
+
+//__________________________________________________________________
 void AliHFEcuts::SetRecPrimaryCutList(){
   //
   // Primary cuts (based on standard cuts from PWG0):
@@ -239,11 +275,13 @@ void AliHFEcuts::SetRecPrimaryCutList(){
   //  No Kink daughters
   //
   AliCFTrackIsPrimaryCuts *primaryCut = new AliCFTrackIsPrimaryCuts("fCutsPrimaryCuts", "REC Primary Cuts");
+#ifdef V4_17_00
   if(IsRequireDCAToVertex()){
     primaryCut->SetDCAToVertex2D(kTRUE);
     primaryCut->SetMaxDCAToVertexXY(fDCAtoVtx[0]);
     primaryCut->SetMaxDCAToVertexZ(fDCAtoVtx[1]);
   }
+#endif
   if(IsRequireSigmaToVertex()){
     primaryCut->SetRequireSigmaToVertex(kTRUE);
     primaryCut->SetMaxNSigmaToVertex(fSigmaToVtx);
@@ -258,24 +296,55 @@ void AliHFEcuts::SetRecPrimaryCutList(){
 }
 
 //__________________________________________________________________
-void AliHFEcuts::SetHFElectronCuts(){
+void AliHFEcuts::SetHFElectronITSCuts(){
   //
-  // Special Cuts introduced by the HFElectron Group
+  // Special Cuts introduced by the HFElectron Group: ITS
   //
-  AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroup","Extra cuts from the HFE group");
+  AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupPixels","Extra cuts from the HFE group");
   if(IsRequireITSpixel()){
     hfecuts->SetRequireITSpixel(AliHFEextraCuts::ITSPixel_t(fCutITSPixel));
   }
-/*  if(IsRequireDCAToVertex()){
+#ifndef V4_17_00
+  if(IsRequireDCAToVertex()){
     hfecuts->SetMaxImpactParamR(fDCAtoVtx[0]);
     hfecuts->SetMaxImpactParamZ(fDCAtoVtx[1]);
-  }*/
-  if(fMinTrackletsTRD) hfecuts->SetMinTrackletsTRD(fMinTrackletsTRD);
+  }
+#endif
+  
+  if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA);
+  
+  TObjArray *hfeCuts = new TObjArray;
+  hfeCuts->SetName("fPartHFECutsITS");
+  hfeCuts->AddLast(hfecuts);
+  fCutList->AddLast(hfeCuts);
+}
+
+//__________________________________________________________________
+void AliHFEcuts::SetHFElectronTPCCuts(){
+  //
+  // Special Cuts introduced by the HFElectron Group: TPC
+  //
+  AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupTPC","Extra cuts from the HFE group");
   if(fMinClusterRatioTPC > 0.) hfecuts->SetClusterRatioTPC(fMinClusterRatioTPC);
   if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA);
   
   TObjArray *hfeCuts = new TObjArray;
-  hfeCuts->SetName("fPartHFECuts");
+  hfeCuts->SetName("fPartHFECutsTPC");
+  hfeCuts->AddLast(hfecuts);
+  fCutList->AddLast(hfeCuts);
+}
+
+//__________________________________________________________________
+void AliHFEcuts::SetHFElectronTRDCuts(){
+  //
+  // Special Cuts introduced by the HFElectron Group: TRD
+  //
+  AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupTRD","Extra cuts from the HFE group");
+  if(fMinTrackletsTRD > 0.) hfecuts->SetMinTrackletsTRD(fMinTrackletsTRD);
+  if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA);
+  
+  TObjArray *hfeCuts = new TObjArray;
+  hfeCuts->SetName("fPartHFECutsTRD");
   hfeCuts->AddLast(hfecuts);
   fCutList->AddLast(hfeCuts);
 }
index b0bce92..9bf43f5 100644 (file)
@@ -47,12 +47,16 @@ class AliHFEcuts : public TObject{
   typedef enum{
     kStepMCGenerated = 0,
     kStepMCInAcceptance = 1,
-    kStepRecKine = 2,
-    kStepRecPrim = 3,
-    kStepHFEcuts = 4
-  } CutStep_t;
+    kStepRecKineTPC = 2,
+    kStepRecKineITS = 3,
+    kStepRecPrim = 4,
+    kStepHFEcutsITS = 5,
+    kStepHFEcutsTPC = 6,
+    kStepHFEcutsTRD = 7
+   } CutStep_t;
 
   static const Int_t kNcutSteps;
+  static const Int_t kNcutESDSteps;
 
   AliHFEcuts();
   AliHFEcuts(const AliHFEcuts &c);
@@ -97,14 +101,20 @@ class AliHFEcuts : public TObject{
   void SetRequireITSPixel() { SETBIT(fRequirements, kITSPixel); }
   void SetRequireProdVertex() { SETBIT(fRequirements, kProductionVertex); };
   void SetRequireSigmaToVertex() { SETBIT(fRequirements, kSigmaToVertex); };
+
+  void SetDebugLevel(Int_t level) { fDebugLevel = level; };
+  Int_t GetDebugLevel() const { return fDebugLevel; };
     
  private:
   void SetParticleGenCutList();
   void SetAcceptanceCutList();
-  void SetRecKineCutList();
+  void SetRecKineTPCCutList();
+  void SetRecKineITSCutList();
   void SetRecPrimaryCutList();
-  void SetHFElectronCuts();
-    
+  void SetHFElectronITSCuts();
+  void SetHFElectronTPCCuts();
+  void SetHFElectronTRDCuts();
+  
   ULong64_t fRequirements;     // Bitmap for requirements
   Double_t fDCAtoVtx[2];       // DCA to Vertex
   Double_t fProdVtx[4];        // Production Vertex
@@ -118,6 +128,8 @@ class AliHFEcuts : public TObject{
     
   TList *fHistQA;              //! QA Histograms
   TObjArray *fCutList; //! List of cut objects(Correction Framework Manager)
+
+  Int_t fDebugLevel;    // Debug Level
     
   ClassDef(AliHFEcuts, 1)   // Container for HFE cuts
     };
@@ -166,7 +178,8 @@ void AliHFEcuts::CreateStandardCuts(){
   fDCAtoVtx[1] = 10.;
   fMinClustersTPC = 50;
   fMinTrackletsTRD = 6;
-  fCutITSPixel = AliHFEextraCuts::kAny;
+  SetRequireITSPixel();
+  fCutITSPixel = AliHFEextraCuts::kBoth;
   fMaxChi2clusterTPC = 3.5;
   fMinClusterRatioTPC = 0.6;
   fPtRange[0] = 0.1;
index 0043891..dedbd79 100644 (file)
@@ -47,7 +47,8 @@ AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
   fClusterRatioTPC(0.),
   fMinTrackletsTRD(0),
   fPixelITS(0),
-  fQAlist(0x0)
+  fQAlist(0x0),
+  fDebugLevel(0)
 {
   //
   // Default Constructor
@@ -63,7 +64,8 @@ AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
   fClusterRatioTPC(c.fClusterRatioTPC),
   fMinTrackletsTRD(c.fMinTrackletsTRD),
   fPixelITS(c.fPixelITS),
-  fQAlist(0x0)
+  fQAlist(0x0),
+  fDebugLevel(c.fDebugLevel)
 {
   //
   // Copy constructor
@@ -89,6 +91,7 @@ AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
     fClusterRatioTPC = c.fClusterRatioTPC;
     fMinTrackletsTRD = c.fMinTrackletsTRD;
     fPixelITS = c.fPixelITS;
+    fDebugLevel = c.fDebugLevel;
 
     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
     if(IsQAOn()){
@@ -144,6 +147,10 @@ Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){
   trdTracklets = track->GetTRDpidQuality();
   #endif  
   UChar_t itsPixel = track->GetITSClusterMap();
+  Int_t det, status1, status2;
+  Float_t xloc, zloc;
+  track->GetITSModuleIndexInfo(0, det, status1, xloc, zloc);
+  track->GetITSModuleIndexInfo(1, det, status2, xloc, zloc);
   if(TESTBIT(fRequirements, kMinImpactParamR)){
     // cut on min. Impact Parameter in Radial direction
     if(impact_r >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
@@ -170,18 +177,38 @@ Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){
   }
   if(TESTBIT(fRequirements, kPixelITS)){
     // cut on ITS pixel layers
-    if(fPixelITS == kFirst)
+    if(fDebugLevel > 0){
+      printf("ITS cluster Map: ");
+      PrintBitMap(itsPixel);
+    }
     switch(fPixelITS){
-      case kFirst: if(itsPixel & BIT(0)) SETBIT(survivedCut, kPixelITS);
-                  break;
-      case kSecond: if(itsPixel & BIT(1)) SETBIT(survivedCut, kPixelITS);
-                  break;
-      case kBoth: if((itsPixel & BIT(0)) && (itsPixel & BIT(1))) SETBIT(survivedCut, kPixelITS);
-                 break;
-      case kAny: if((itsPixel & BIT(0)) || (itsPixel & BIT(1))) SETBIT(survivedCut, kPixelITS);
-                break;
+      case kFirst: 
+          if(itsPixel & BIT(0)) 
+            SETBIT(survivedCut, kPixelITS);
+          else if(!CheckITSstatus(status1))
+            SETBIT(survivedCut, kPixelITS);
+                   break;
+      case kSecond: 
+          if(itsPixel & BIT(1)) 
+            SETBIT(survivedCut, kPixelITS);
+          else if(!CheckITSstatus(status2))
+            SETBIT(survivedCut, kPixelITS);
+                   break;
+      case kBoth: 
+          if((itsPixel & BIT(0)) && (itsPixel & BIT(1))) 
+            SETBIT(survivedCut, kPixelITS);
+          else if(!CheckITSstatus(status1) && !CheckITSstatus(status2))
+            SETBIT(survivedCut, kPixelITS);
+                   break;
+      case kAny: 
+          if((itsPixel & BIT(0)) || (itsPixel & BIT(1))) 
+            SETBIT(survivedCut, kPixelITS);
+          else if(!CheckITSstatus(status1) || !CheckITSstatus(status2))
+            SETBIT(survivedCut, kPixelITS);
+                   break;
       default: break;
     }
+    if(fDebugLevel > 0)printf("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO");
   }
   if(fRequirements == survivedCut){
     //
@@ -321,3 +348,25 @@ void AliHFEextraCuts::AddQAHistograms(TList *qaList){
   }
   qaList->AddLast(fQAlist);
 }
+
+//______________________________________________________
+void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
+  for(Int_t ibit = 32; ibit--; )
+    printf("%d", bitmap & BIT(ibit) ? 1 : 0);
+  printf("\n");
+}
+
+//______________________________________________________
+Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus){
+  //
+  // Check whether ITS area is dead
+  //
+  Bool_t status;
+  switch(itsStatus){
+    case 2: status = kFALSE; break;
+    case 3: status = kFALSE; break;
+    case 7: status = kFALSE; break;
+    default: status = kTRUE;
+  }
+  return status;
+}
index 8a22fff..6aadae9 100644 (file)
@@ -48,14 +48,19 @@ class AliHFEextraCuts : public AliCFCutBase{
     inline void SetMinImpactParamZ(Double_t impactParam);
     inline void SetMaxImpactParamZ(Double_t impactParam);
     inline void SetMinTrackletsTRD(Int_t minTracklets);
+
+    void SetDebugLevel(Int_t level) { fDebugLevel = level; };
+    Int_t GetDebugLevel() const { return fDebugLevel; };
     
   protected:
     virtual void AddQAHistograms(TList *qaList);
     Bool_t CheckESDCuts(AliESDtrack *track);
     Bool_t CheckMCCuts(AliMCParticle */*track*/);
+    Bool_t CheckITSstatus(Int_t itsStatus);
     void FillQAhistosESD(AliESDtrack *track, UInt_t when);
 //     void FillQAhistosMC(AliMCParticle *track, UInt_t when);
     void FillCutCorrelation(ULong64_t survivedCut);
+    void PrintBitMap(Int_t bitmap);
     
   private:
     typedef enum{
@@ -84,6 +89,8 @@ class AliHFEextraCuts : public AliCFCutBase{
     
     TList *fQAlist;                    //! Directory for QA histograms
   
+    Int_t fDebugLevel;    // Debug Level
+  
   ClassDef(AliHFEextraCuts, 1)      // Additional cuts implemented by the ALICE HFE group
 };
 
index bde5e4e..e359c93 100644 (file)
 /**************************************************************************
  *                                                                        *
  * QA class of Heavy Flavor quark and fragmeted/decayed particles         *
+ * -Check kinematics of Heavy Quarks/hadrons, and decayed leptons         *
+ *    pT, rapidity                                                        *
+ *    decay lepton kinematics w/wo acceptance                             *
+ *    heavy hadron decay length, electron pT fraction carried from decay  *
+ * -Check yield of Heavy Quarks/hadrons                                   *
+ *    Number of produced heavy quark                                      *
+ *    Number of produced hadron of given pdg code                         *
+ *                                                                        *
  *                                                                        *
  * Authors:                                                               *
  *   MinJung Kweon <minjung@physi.uni-heidelberg.de>                      *
 
 const Int_t AliHFEmcQA::fgkGluon=21;
 const Int_t AliHFEmcQA::fgkMaxGener=10;
-const Int_t AliHFEmcQA::fgkMaxIter=1000;
-const Int_t AliHFEmcQA::fgkqType = 6;    // number of species waiting for QA done
+const Int_t AliHFEmcQA::fgkMaxIter=100;
+const Int_t AliHFEmcQA::fgkqType = 7;    // number of species waiting for QA done
 
 ClassImp(AliHFEmcQA)
 
 //_______________________________________________________________________________________________
 AliHFEmcQA::AliHFEmcQA() : 
-        fVerbos(kFALSE) 
-        ,fStack(0x0) 
+        fStack(0x0) 
         ,fNparents(0) 
 {
         // Default constructor
         
-        if (fVerbos) AliInfo("***** Warning! fVerbos is set to TRUE! ******");
 }
 
 //_______________________________________________________________________________________________
 AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
         TObject(p)
-        ,fVerbos(p.fVerbos)
         ,fStack(0x0) 
         ,fNparents(p.fNparents) 
 {
@@ -90,63 +95,68 @@ AliHFEmcQA::~AliHFEmcQA()
 }
 
 //_______________________________________________________________________________________________
-void AliHFEmcQA::PostAnalyze()
+const void AliHFEmcQA::PostAnalyze()
 {
 }
 
 //__________________________________________
-void AliHFEmcQA::CreateHistograms(const Int_t kquark, TString hnopt) 
+void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) 
 {
   // create histograms
 
-  if (kquark != fkCharm && kquark != fkBeauty) {
+  if (kquark != kCharm && kquark != kBeauty) {
     AliDebug(1, "This task is only for heavy quark QA, return\n");
     return; 
   }
-  Int_t iq = kquark - fkCharm; 
+  Int_t iq = kquark - kCharm; 
 
   TString kqTypeLabel[fgkqType];
-  if (kquark == fkCharm){
-    kqTypeLabel[fkQuark]="c";
-    kqTypeLabel[fkantiQuark]="cbar";
-    kqTypeLabel[fkElectron]="ce";
-    kqTypeLabel[fkElectron2nd]="nulle";
-    kqTypeLabel[fkeHadron]="ceHadron";
-    kqTypeLabel[fkDeHadron]="nullHadron";
-  } else if (kquark == fkBeauty){
-    kqTypeLabel[fkQuark]="b";
-    kqTypeLabel[fkantiQuark]="bbar";
-    kqTypeLabel[fkElectron]="be";
-    kqTypeLabel[fkElectron2nd]="bce";
-    kqTypeLabel[fkeHadron]="beHadron";
-    kqTypeLabel[fkDeHadron]="bDeHadron";
+  if (kquark == kCharm){
+    kqTypeLabel[kQuark]="c";
+    kqTypeLabel[kantiQuark]="cbar";
+    kqTypeLabel[kHadron]="cHadron";
+    kqTypeLabel[keHadron]="ceHadron";
+    kqTypeLabel[kDeHadron]="nullHadron";
+    kqTypeLabel[kElectron]="ce";
+    kqTypeLabel[kElectron2nd]="nulle";
+  } else if (kquark == kBeauty){
+    kqTypeLabel[kQuark]="b";
+    kqTypeLabel[kantiQuark]="bbar";
+    kqTypeLabel[kHadron]="bHadron";
+    kqTypeLabel[keHadron]="beHadron";
+    kqTypeLabel[kDeHadron]="bDeHadron";
+    kqTypeLabel[kElectron]="be";
+    kqTypeLabel[kElectron2nd]="bce";
   }
 
 
   TString hname; 
   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
+     if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
      hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
-     fHist[iq][iqType].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
+     fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
      hname = hnopt+"Pt_"+kqTypeLabel[iqType];
-     fHist[iq][iqType].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",150,0,30);
+     fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",150,0,30);
      hname = hnopt+"Y_"+kqTypeLabel[iqType];
-     fHist[iq][iqType].fY = new TH1F(hname,hname,150,-7.5,7.5);
+     fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
      hname = hnopt+"Eta_"+kqTypeLabel[iqType];
-     fHist[iq][iqType].fEta = new TH1F(hname,hname,150,-7.5,7.5);
+     fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
   }
 
-  hname = hnopt+"Nq_"+kqTypeLabel[fkQuark];
-  fHistComm[iq].fNq = new TH1F(hname,hname,10,-0.5,9.5);
-  hname = hnopt+"ProcessID_"+kqTypeLabel[fkQuark];
-  fHistComm[iq].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
-  hname = hnopt+"ePtRatio_"+kqTypeLabel[fkQuark];
-  fHistComm[iq].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
-  hname = hnopt+"DePtRatio_"+kqTypeLabel[fkQuark];
-  fHistComm[iq].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
-  hname = hnopt+"eDistance_"+kqTypeLabel[fkQuark];
-  fHistComm[iq].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
-  hname = hnopt+"DeDistance_"+kqTypeLabel[fkQuark];
-  fHistComm[iq].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
+  if (icut == 0){ 
+    hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
+    hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
+  }
+  hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
+  fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
+  hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
+  fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1);
+  hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
+  fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
+  hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
+  fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2);
 
 }
 
@@ -184,11 +194,11 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark)
 {
   // get heavy quark kinematics
 
-    if (kquark != fkCharm && kquark != fkBeauty) {
+    if (kquark != kCharm && kquark != kBeauty) {
       AliDebug(1, "This task is only for heavy quark QA, return\n");
       return; 
     }
-    Int_t iq = kquark - fkCharm; 
+    Int_t iq = kquark - kCharm; 
 
 
     if (iTrack < 0) { 
@@ -226,7 +236,7 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark)
 
         if ( abs(partMother->GetPdgCode()) != kquark ){
           // search fragmented partons in the same string
-          Bool_t IsSameString = kTRUE; 
+          Bool_t isSameString = kTRUE; 
           for (Int_t i=1; i<fgkMaxIter; i++){
              iLabel = iLabel - 1;
              if (iLabel>-1) { partMother = fStack->Particle(iLabel); }
@@ -235,8 +245,8 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark)
                return; 
              }
              if ( abs(partMother->GetPdgCode()) == kquark ) break;
-             if ( partMother->GetStatusCode() != 12 ) IsSameString = kFALSE;
-             if (!IsSameString) return; 
+             if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE;
+             if (!isSameString) return; 
           }
         }
         AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n");
@@ -248,15 +258,15 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark)
 
         // fill kinematics for heavy parton
         if (partMother->GetPdgCode() > 0) { // quark
-          fHist[iq][fkQuark].fPdgCode->Fill(partMother->GetPdgCode());
-          fHist[iq][fkQuark].fPt->Fill(partMother->Pt());
-          fHist[iq][fkQuark].fY->Fill(GetRapidity(partMother));
-          fHist[iq][fkQuark].fEta->Fill(partMother->Eta());
+          fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
+          fHist[iq][kQuark][0].fPt->Fill(partMother->Pt());
+          fHist[iq][kQuark][0].fY->Fill(GetRapidity(partMother));
+          fHist[iq][kQuark][0].fEta->Fill(partMother->Eta());
         } else{ // antiquark
-          fHist[iq][fkantiQuark].fPdgCode->Fill(partMother->GetPdgCode());
-          fHist[iq][fkantiQuark].fPt->Fill(partMother->Pt());
-          fHist[iq][fkantiQuark].fY->Fill(GetRapidity(partMother));
-          fHist[iq][fkantiQuark].fEta->Fill(partMother->Eta());
+          fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode());
+          fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt());
+          fHist[iq][kantiQuark][0].fY->Fill(GetRapidity(partMother));
+          fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta());
         }
 
       } // end of heavy parton slection loop 
@@ -270,16 +280,16 @@ void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
 {
   // end of event analysis
 
-  if (kquark != fkCharm && kquark != fkBeauty) {
+  if (kquark != kCharm && kquark != kBeauty) {
     AliDebug(1, "This task is only for heavy quark QA, return\n");
     return; 
   }
-  Int_t iq = kquark - fkCharm; 
+  Int_t iq = kquark - kCharm; 
 
 
   // # of heavy quark per event
   AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq]));
-  fHistComm[iq].fNq->Fill(fIsHeavy[iq]);
+  fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]);
 
   Int_t motherID[fgkMaxGener];
   Int_t motherType[fgkMaxGener];
@@ -318,7 +328,7 @@ void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
           if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
           if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break;
           // if it is not the above case, something is strange
-          reportStrangeness(motherID[i],motherType[i],motherLabel[i]);
+          ReportStrangeness(motherID[i],motherType[i],motherLabel[i]);
           break;
      } 
      if (ancestorLabel[ig] == -1){ // from hard scattering
@@ -342,47 +352,117 @@ void AliHFEmcQA::EndOfEventAna(const Int_t kquark)
      Int_t id2 = ipair*2 + 1;
 
      if (motherType[id1] == 2 && motherType[id2] == 2){
-       if (motherLabel[id1] == motherLabel[id2]) processID = fkGluonSplitting; // gluon spliting
+       if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting
        else processID = -9;
      }
      else if (motherType[id1] == -1 && motherType[id2] == -1) {
        if (motherLabel[id1] == -1 && motherLabel[id2] == -1) {
-         if (motherID[id1] == fgkGluon) processID = fkPairCreationFromg; // gluon fusion
-         else processID = fkPairCreationFromq; // q-qbar pair creation
+         if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion
+         else processID = kPairCreationFromq; // q-qbar pair creation
        }
        else processID = -8;
      }
      else if (motherType[id1] == -1 || motherType[id2] == -1) {
        if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) {
-         if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = fkFlavourExitation; // flavour exitation 
-         else processID = fkLightQuarkShower;
+         if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation 
+         else processID = kLightQuarkShower;
        }
        else processID = -7;
      }
      else if (motherType[id1] == -2 || motherType[id2] == -2) {
-       if (motherLabel[id1] == motherLabel[id2]) processID = fkInitialPartonShower; // initial parton shower
+       if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower
        else processID = -6;
        
      }
      else processID = -5;
 
      if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID));
-     else fHistComm[iq].fProcessID->Fill(processID);
+     else fHistComm[iq][0].fProcessID->Fill(processID);
      AliDebug(1,Form("Process ID = %d\n",processID));
   } // end of # heavy quark pair loop
 
 }
 
 //__________________________________________
-void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed, Bool_t isbarrel) 
+void AliHFEmcQA::GetHadronKine(Int_t iTrack, const Int_t kquark)
+{
+    // decay electron kinematics
+
+    if (kquark != kCharm && kquark != kBeauty) {
+      AliDebug(1, "This task is only for heavy quark QA, return\n");
+      return;
+    }
+    Int_t iq = kquark - kCharm;
+
+    if (iTrack < 0) {
+      AliDebug(1, "Stack label is negative, return\n");
+      return;
+    }
+
+    TParticle* mcpart = fStack->Particle(iTrack);
+
+    Int_t iLabel = mcpart->GetFirstMother();
+    if (iLabel<0){
+      AliDebug(1, "Stack label is negative, return\n");
+      return;
+    }
+
+    TParticle *partCopy = mcpart;
+    Int_t pdgcode = mcpart->GetPdgCode();
+    Int_t pdgcodeCopy = pdgcode;
+
+    // if the mother is charmed hadron  
+    Bool_t IsDirectCharm = kFALSE;
+    if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) {
+
+          // iterate until you find B hadron as a mother or become top ancester 
+          for (Int_t i=1; i<fgkMaxIter; i++){
+
+             Int_t jLabel = mcpart->GetFirstMother();
+             if (jLabel == -1){
+               IsDirectCharm = kTRUE;
+               break; // if there is no ancester
+             }
+             if (jLabel < 0){ // safety protection
+               AliDebug(1, "Stack label is negative, return\n");
+               return;
+             }
+             // if there is an ancester
+             TParticle* mother = fStack->Particle(jLabel);
+             Int_t motherPDG = mother->GetPdgCode();
+    
+             for (Int_t j=0; j<fNparents; j++){
+                if (abs(motherPDG)==fParentSelect[1][j]) return; // return if this hadron is originated from b
+             }
+
+             mcpart = mother;
+          } // end of iteration 
+    } // end of if
+    if((IsDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
+         for (Int_t i=0; i<fNparents; i++){
+            if (abs(pdgcodeCopy)==fParentSelect[iq][i]){
+
+              // fill hadron kinematics
+              fHist[iq][kHadron][0].fPdgCode->Fill(pdgcodeCopy);
+              fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt());
+              fHist[iq][kHadron][0].fY->Fill(GetRapidity(partCopy));
+              fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta());
+
+            }
+         }
+    } // end of if
+}
+
+//__________________________________________
+void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed, Int_t icut, Bool_t isbarrel) 
 {
     // decay electron kinematics
     
-    if (kquark != fkCharm && kquark != fkBeauty) {
+    if (kquark != kCharm && kquark != kBeauty) {
       AliDebug(1, "This task is only for heavy quark QA, return\n");
       return; 
     }
-    Int_t iq = kquark - fkCharm; 
+    Int_t iq = kquark - kCharm; 
 
     if (iTrack < 0) { 
       AliDebug(1, "Stack label is negative, return\n");
@@ -401,25 +481,24 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed
     }
 
     TParticle *partMother = fStack->Particle(iLabel);
+    TParticle *partMotherCopy = partMother;
     Int_t maPdgcode = partMother->GetPdgCode();
-    TParticle *partMotherCopy = fStack->Particle(iLabel);
-    Int_t maPdgcodeCopy = partMotherCopy->GetPdgCode();
+    Int_t maPdgcodeCopy = maPdgcode;
 
     // get electron production vertex   
     TLorentzVector ePoint;
     mcpart->ProductionVertex(ePoint);
 
-
     // if the mother is charmed hadron  
-    Bool_t IsMotherDirectCharm = kFALSE;
-    if ( int(abs(maPdgcode)/100.) == fkCharm || int(abs(maPdgcode)/1000.) == fkCharm ) { 
+    Bool_t isMotherDirectCharm = kFALSE;
+    if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { 
 
           // iterate until you find B hadron as a mother or become top ancester 
           for (Int_t i=1; i<fgkMaxIter; i++){
 
              Int_t jLabel = partMother->GetFirstMother(); 
              if (jLabel == -1){
-               IsMotherDirectCharm = kTRUE;
+               isMotherDirectCharm = kTRUE;
                break; // if there is no ancester
              }
              if (jLabel < 0){ // safety protection
@@ -434,27 +513,27 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed
              for (Int_t j=0; j<fNparents; j++){
                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 
-                  if (kquark == fkCharm) return;
+                  if (kquark == kCharm) return;
                   // fill electron kinematics
-                  fHist[iq][fkElectron2nd].fPdgCode->Fill(mcpart->GetPdgCode());
-                  fHist[iq][fkElectron2nd].fPt->Fill(mcpart->Pt());
-                  fHist[iq][fkElectron2nd].fY->Fill(GetRapidity(mcpart));
-                  fHist[iq][fkElectron2nd].fEta->Fill(mcpart->Eta());
+                  fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
+                  fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
+                  fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart));
+                  fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
 
                   // fill mother hadron kinematics
-                  fHist[iq][fkDeHadron].fPdgCode->Fill(grandMaPDG); 
-                  fHist[iq][fkDeHadron].fPt->Fill(grandMa->Pt());
-                  fHist[iq][fkDeHadron].fY->Fill(GetRapidity(grandMa));
-                  fHist[iq][fkDeHadron].fEta->Fill(grandMa->Eta());
+                  fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); 
+                  fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
+                  fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa));
+                  fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
                  
                   // ratio between pT of electron and pT of mother B hadron 
-                  if(grandMa->Pt()) fHistComm[iq].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                  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].fDeDistance->Fill(grandMa->Pt(),dedistance.P());
+                  fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),dedistance.P());
                   return;
                 }
              }
@@ -462,30 +541,30 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed
              partMother = grandMa;
           } // end of iteration 
     } // end of if
-    if((IsMotherDirectCharm == kTRUE && kquark == fkCharm) || kquark == fkBeauty) {
+    if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) {
          for (Int_t i=0; i<fNparents; i++){
             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 
               // fill electron kinematics
-              fHist[iq][fkElectron].fPdgCode->Fill(mcpart->GetPdgCode());
-              fHist[iq][fkElectron].fPt->Fill(mcpart->Pt());
-              fHist[iq][fkElectron].fY->Fill(GetRapidity(mcpart));
-              fHist[iq][fkElectron].fEta->Fill(mcpart->Eta());  
+              fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
+              fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
+              fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart));
+              fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());  
 
               // fill mother hadron kinematics
-              fHist[iq][fkeHadron].fPdgCode->Fill(maPdgcodeCopy); 
-              fHist[iq][fkeHadron].fPt->Fill(partMotherCopy->Pt());
-              fHist[iq][fkeHadron].fY->Fill(GetRapidity(partMotherCopy));
-              fHist[iq][fkeHadron].fEta->Fill(partMotherCopy->Eta());
+              fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); 
+              fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
+              fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy));
+              fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
 
               // ratio between pT of electron and pT of mother B hadron 
-              if(partMotherCopy->Pt()) fHistComm[iq].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
+              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].feDistance->Fill(partMotherCopy->Pt(),edistance.P());
+              fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),edistance.P());
             }
          }
     } // end of if
@@ -495,6 +574,8 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed
 //__________________________________________
 void AliHFEmcQA::IdentifyMother(Int_t mother_label, Int_t &mother_pdg, Int_t &grandmother_label)
 {
+       // find mother pdg code and label 
+
        if (mother_label < 0) { 
          AliDebug(1, "Stack label is negative, return\n");
          return; 
@@ -506,9 +587,10 @@ void AliHFEmcQA::IdentifyMother(Int_t mother_label, Int_t &mother_pdg, Int_t &gr
 }
 
 //__________________________________________
-// mothertype -1 means this heavy quark coming from hard vertex
 void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 {
+       // mothertype -1 means this heavy quark coming from hard vertex
+
        TParticle *afterinitialrad1  = fStack->Particle(4);
        TParticle *afterinitialrad2  = fStack->Particle(5);
            
@@ -538,9 +620,10 @@ void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &moth
 }
 
 //__________________________________________
-// mothertype -2 means this heavy quark coming from initial state 
 Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 {
+       // mothertype -2 means this heavy quark coming from initial state 
+
        if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation
          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
          motherID = heavysMother->GetPdgCode(); 
@@ -555,9 +638,10 @@ Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID,
 }
 
 //__________________________________________
-// mothertype 2 means this heavy quark coming from final state 
 Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 {
+       // mothertype 2 means this heavy quark coming from final state 
+
        if (inputmotherlabel > 5){ // mother exist after hard scattering
          TParticle *heavysMother = fStack->Particle(inputmotherlabel);
          motherID = heavysMother->GetPdgCode(); 
@@ -571,8 +655,10 @@ Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, In
 }
 
 //__________________________________________
-void AliHFEmcQA::reportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
+void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel)
 {
+      // mark strange behavior  
+
        mothertype = -888;
        motherlabel = -888;
        motherID = -888;
@@ -580,8 +666,10 @@ void AliHFEmcQA::reportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &mo
 }
 
 //__________________________________________
-Float_t AliHFEmcQA::GetRapidity(TParticle *part)
+const Float_t AliHFEmcQA::GetRapidity(TParticle *part)
 {
+      // return rapidity
+
        Float_t rapidity;        
        if(part->Energy() - part->Pz() == 0 || part->Energy() + part->Pz() == 0) rapidity=-999; 
        else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); 
index df12d8b..de8a243 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-
 /**************************************************************************
  *                                                                        *
  * QA class of Heavy Flavor quark and fragmeted/decayed particles         *
+ * -Check kinematics of Heavy Quarks/hadrons, and decayed leptons         *
+ *    pT, rapidity                                                        *
+ *    decay lepton kinematics w/wo acceptance                             *
+ *    heavy hadron decay length, electron pT fraction carried from decay  *
+ * -Check yield of Heavy Quarks/hadrons                                   *
+ *    Number of produced heavy quark                                      *
+ *    Number of produced hadron of given pdg code                         *
  *                                                                        *
  **************************************************************************/
 
-#ifndef _ALIHFEMCQA_H_
-#define _ALIHFEMCQA_H_
+#ifndef ALIHFEMCQA_H
+#define ALIHFEMCQA_H
 
 #ifndef ROOT_TObject
 #include <TObject.h>
@@ -36,49 +42,45 @@ class AliStack;
 class AliHFEmcQA: public TObject {
 
         public: 
-                enum heavyType {fkCharm=4, fkBeauty=5};
-                enum qType {fkQuark, fkantiQuark, fkElectron, fkElectron2nd, fkeHadron, fkDeHadron};
+                enum heavyType {kCharm=4, kBeauty=5, kElectronPDG=11};
+                enum qType {kQuark, kantiQuark, kHadron, keHadron, kDeHadron, kElectron, kElectron2nd};
                 AliHFEmcQA();
                 AliHFEmcQA(const AliHFEmcQA &p); // copy constructor
                 AliHFEmcQA &operator=(const AliHFEmcQA &); // assignment operator
 
                 virtual ~AliHFEmcQA();
 
-                void PostAnalyze();
-                void SetVerbosity(Bool_t verb=kTRUE) {fVerbos=verb;}
-
-        protected:
-                void SetStack(AliStack* stack){fStack=stack;} // set stack pointer
-                void CreateHistograms(const Int_t kquark, TString hnopt=""); // create histograms for mc qa analysis
+                const void PostAnalyze();
+                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 Init();
+
                 void GetQuarkKine(Int_t iTrack, const Int_t kquark); // get heavy quark kinematics distribution
-                void GetDecayedKine(Int_t iTrack, const Int_t kquark, const Int_t kdecayed, Bool_t isbarrel=kFALSE); // get decay electron kinematics distribution
+                void GetHadronKine(Int_t iTrack, const Int_t kquark); // get heavy hadron kinematics distribution
+                void GetDecayedKine(Int_t iTrack, const Int_t kquark, const Int_t kdecayed, Int_t icut, Bool_t isbarrel=kFALSE); // get decay electron kinematics distribution
                 void EndOfEventAna(const Int_t kquark); // run analysis which should be done at the end of the event loop
+
+        protected:
                 void IdentifyMother(Int_t mother_label, Int_t &mother_pdg, Int_t &grandmother_label); // 
                 void HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // check if the quark is produced from hard scattering
-                void reportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // report if the quark production process is unknown
+                void ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // report if the quark production process is unknown
                 Bool_t IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // check if the quark is produced from initial parton shower 
                 Bool_t IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // check if the quark is produced from final parton shower
-                Float_t GetRapidity(TParticle *part); // return rapidity
+                const Float_t GetRapidity(TParticle *part); // return rapidity
 
-                Bool_t fVerbos; //
+                AliStack* fStack; // stack pointer           
 
-                AliStack* fStack; //            
-
-                static const Int_t fgkGluon; //
-                static const Int_t fgkPDGInterested; //
-                static const Int_t fgkMaxGener; //
-                static const Int_t fgkMaxIter; //
-                static const Int_t fgkqType; //
+                static const Int_t fgkGluon; // gluon pdg code
+                static const Int_t fgkMaxGener; // ancester level wanted to be checked 
+                static const Int_t fgkMaxIter; // number of iteration to find out matching particle 
+                static const Int_t fgkqType; // number of particle type to be checked
 
 
                 enum ProcessType_t
                         {
-                        fkPairCreationFromq,  fkPairCreationFromg,  fkFlavourExitation,  fkGluonSplitting, fkInitialPartonShower, fkLightQuarkShower
+                        kPairCreationFromq,  kPairCreationFromg,  kFlavourExitation,  kGluonSplitting, kInitialPartonShower, kLightQuarkShower
                         };
 
-                TString fkqTypeLabel[6]; //
-
                 struct hists{
                         TH1F *fPdgCode; // histogram to store particle pdg code
                         TH1F *fPt; // histogram to store pt
@@ -88,22 +90,22 @@ class AliHFEmcQA: public TObject {
                 struct histsComm {
                         TH1F *fNq; // histogram to store number of quark
                         TH1F *fProcessID; // histogram to store process id 
-                        TH2F *fePtRatio; //
-                        TH2F *fDePtRatio; //
-                        TH2F *feDistance; //
-                        TH2F *fDeDistance; //
+                        TH2F *fePtRatio; // fraction of electron pT from D or B hadron
+                        TH2F *fDePtRatio; // fraction of D electron pT from B hadron 
+                        TH2F *feDistance; // distance between electron production point to mother particle 
+                        TH2F *fDeDistance; // distance between D electron production point to mother particle
                 };
 
-                hists fHist[2][6]; //
-                histsComm fHistComm[2]; //
+                hists fHist[2][7][5]; // struct of histograms to store kinematics of given particles
+                histsComm fHistComm[2][5]; // struct of additional histograms of given particles
 
-                TParticle *fHeavyQuark[50]; //
-                Int_t fIsHeavy[2]; //
-                Int_t fNparents; //
-                Int_t fParentSelect[2][7]; //
+                TParticle *fHeavyQuark[50]; // store pointer of heavy flavour quark 
+                Int_t fIsHeavy[2]; // count of heavy flavour
+                Int_t fNparents; // number of heavy hadrons to be considered
+                Int_t fParentSelect[2][7]; // heavy hadron species
 
 
-               ClassDef(AliHFEmcQA,0); // HFE Monte Carlo quality assurance
+        ClassDef(AliHFEmcQA,0);
 };
 
 #endif
index 89973d5..68e9d29 100644 (file)
@@ -44,28 +44,23 @@ ClassImp(AliHFEpidTRD)
 //___________________________________________________________________
 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
     AliHFEpidBase(name)
-  , fThresholdFile("TRD.PIDthresholds.root")
   , fPIDMethod(kNN)
-  , fTRDthresholds(0x0)
-  , fTRDelectronEfficiencies(0x0)
 {
   //
   // default  constructor
   // 
-  fTRDthresholds = new TMap();
+  memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
 }
 
 //___________________________________________________________________
 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
     AliHFEpidBase("")
-  , fThresholdFile("")
   , fPIDMethod(kLQ)
-  , fTRDthresholds(0x0)
-  , fTRDelectronEfficiencies(0x0)
 {
   //
   // Copy constructor
   //
+  memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
   ref.Copy(*this);
 }
 
@@ -87,14 +82,8 @@ void AliHFEpidTRD::Copy(TObject &ref) const {
   //
   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
 
-  target.fThresholdFile = fThresholdFile;
   target.fPIDMethod = fPIDMethod;
-  if(fTRDthresholds){
-    target.fTRDthresholds = dynamic_cast<TMap *>(fTRDthresholds->Clone());
-  }
-  if(fTRDelectronEfficiencies){
-    target.fTRDelectronEfficiencies = new TAxis(*fTRDelectronEfficiencies);
-  }
+  memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
   AliHFEpidBase::Copy(ref);
 }
 
@@ -103,11 +92,6 @@ AliHFEpidTRD::~AliHFEpidTRD(){
   //
   // Destructor
   //
-  if(fTRDthresholds){
-    fTRDthresholds->Delete();
-    delete fTRDthresholds;
-  }
-  if(fTRDelectronEfficiencies) delete fTRDelectronEfficiencies;
 }
 
 //______________________________________________________
@@ -116,44 +100,7 @@ Bool_t AliHFEpidTRD::InitializePID(){
   // InitializePID: Load TRD thresholds and create the electron efficiency axis
   // to navigate 
   //
-  LoadTRDthresholds();
-  // Fill the electron efficiencies axis for the TRD alone PID
-  Int_t nEffs = fTRDthresholds->GetEntries() + 1;
-  TIterator* thres =fTRDthresholds->MakeIterator();
-  TObjString *key = 0x0;
-  Double_t *tmp = new Double_t[nEffs], *titer = tmp;
-  while((key = dynamic_cast<TObjString *>((*thres)()))){
-    (*titer++) = static_cast<Double_t>(key->String().Atoi())/100.;
-  }
-  delete thres;
-  *titer = 1.;
-  // Sort the electron efficiencies and put them into the TAxis for later navigation
-  Int_t *ind = new Int_t[nEffs], *iiter = ind;
-  TMath::Sort(nEffs, tmp, ind, kFALSE);
-  Double_t *eleffs = new Double_t[nEffs], *eiter = eleffs;
-  while(eiter < &eleffs[nEffs]) *(eiter++) = tmp[*(iiter++)];
-  // print the content
-  Int_t cnt = 0;
-  if(GetDebugLevel() > 1){
-    printf("Printing electron efficiency bins:\n");
-    eiter = eleffs;
-    thres = fTRDthresholds->MakeIterator();
-    TObject *object = 0x0;
-    while(eiter < &eleffs[nEffs - 1]){
-      printf("eleffs[%d] = %f", cnt++, *(eiter++));
-      key = dynamic_cast<TObjString *>(thres->Next());
-      object = fTRDthresholds->GetValue(key->String().Data());
-      printf(", Content: %p\n", (void *)object);
-    }
-  }
-  delete[] tmp; delete[] ind;
-  fTRDelectronEfficiencies = new TAxis(nEffs - 1, eleffs);
-  if(GetDebugLevel() > 1){
-    printf("Printing axis content:\n");
-    for(Int_t ibin = fTRDelectronEfficiencies->GetFirst(); ibin <= fTRDelectronEfficiencies->GetLast(); ibin++)
-      printf("%d.) minimum: %f, maximum? %f\n", ibin, fTRDelectronEfficiencies->GetBinLowEdge(ibin), fTRDelectronEfficiencies->GetBinUpEdge(ibin));
-  }
-  delete[] eleffs;
+  InitParameters();
   return kTRUE;
 }
 
@@ -169,71 +116,60 @@ Int_t AliHFEpidTRD::IsSelected(AliVParticle *track){
   Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
   if(p < 0.6) return 0;
 
-  // Get the Histograms
-  TH1 *threshist = GetTRDthresholds(0.91);
-  Int_t bin = 0;
-  if(p > threshist->GetXaxis()->GetXmax()) 
-    bin = threshist->GetXaxis()->GetLast();
-  else if(p < threshist->GetXaxis()->GetXmin()) 
-    bin = threshist->GetXaxis()->GetFirst();
-  else
-    bin = threshist->GetXaxis()->FindBin(p);
-
   Double_t pidProbs[AliPID::kSPECIES];
   esdTrack->GetTRDpid(pidProbs);
-  if(pidProbs[AliPID::kElectron] > threshist->GetBinContent(bin)) return 11;
+  if(pidProbs[AliPID::kElectron] > GetTRDthresholds(0.91, p)) return 11;
   return 0;
 }
 
 //___________________________________________________________________
-void AliHFEpidTRD::LoadTRDthresholds(){
-  //
-  // Load TRD threshold histograms from File
-  //
-  TFile *mythresholds = TFile::Open(fThresholdFile);
-  TKey *object = 0x0;
-  TString electron_eff;
-  TObjArray *histos = 0x0;
-  Float_t eff;
-  TH1F *refhist = 0x0;
-  TIterator *keyIterator = mythresholds->GetListOfKeys()->MakeIterator();
-  TString histnames[2] = {"fHistThreshLQ", "fHistThreshNN"};
-  gROOT->cd();
-  while((object = dynamic_cast<TKey *>((*keyIterator)()))){
-    // Get the electron efficiency bin this histogram was taken with
-    electron_eff = object->GetName();
-    electron_eff = electron_eff.Remove(0,3);
-    eff = static_cast<Float_t>(electron_eff.Atoi())/100.;
-
-    // Get the threshold according to the selected 
-    histos = dynamic_cast<TObjArray *>(object->ReadObj());
-    refhist = dynamic_cast<TH1F *>(histos->FindObject(histnames[fPIDMethod].Data()));
-    SetTRDthresholds(refhist, eff);
-    histos->Delete();
-    delete histos;
-  }
-  delete keyIterator;
-  mythresholds->Close();
-  delete mythresholds;
+Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){ 
+  //
+  // Return momentum dependent and electron efficiency dependent TRD thresholds
+  // 
+  Double_t params[4];
+  GetParameters(electronEff, params);
+  return 1. - params[0] * p - params[1] * TMath::Exp(-params[2] * p) - params[3];
 }
 
 //___________________________________________________________________
-void AliHFEpidTRD::SetTRDthresholds(TH1F *thresholds, Float_t electronEff){
-  //
-  // Set the threshold histogram for the TRD pid
-  //
-  fTRDthresholds->Add(new TObjString(Form("%d", TMath::Nint(electronEff * 100.))), new TH1F(*thresholds));
+void AliHFEpidTRD::InitParameters(){
+  //
+  // Fill the Parameters into an array
+  //
+
+  // Parameters for 6 Layers
+  fThreshParams[0] = -0.001839; // 0.7 electron eff
+  fThreshParams[1] = 0.000276;
+  fThreshParams[2] = 0.044902; 
+  fThreshParams[3] = 1.726751;
+  fThreshParams[4] = -0.002405; // 0.75 electron eff
+  fThreshParams[5] = 0.000372;
+  fThreshParams[6] = 0.061775;
+  fThreshParams[7] = 1.739371;
+  fThreshParams[8] = -0.003178; // 0.8 electron eff
+  fThreshParams[9] = 0.000521;
+  fThreshParams[10] = 0.087585;
+  fThreshParams[11] = 1.749154;
+  fThreshParams[12] = -0.004058; // 0.85 electron eff
+  fThreshParams[13] = 0.000748;
+  fThreshParams[14] = 0.129583;
+  fThreshParams[15] = 1.782323;
+  fThreshParams[16] = -0.004967; // 0.9 electron eff
+  fThreshParams[17] = 0.001216;
+  fThreshParams[18] = 0.210128;
+  fThreshParams[19] = 1.807665;
+  fThreshParams[20] = -0.000996; // 0.95 electron eff
+  fThreshParams[21] = 0.002627;
+  fThreshParams[22] = 0.409099;
+  fThreshParams[23] = 1.787076;
 }
 
 //___________________________________________________________________
-TH1F *AliHFEpidTRD::GetTRDthresholds(Float_t electronEff){ 
-  Int_t bin = 0;
-  if(electronEff < fTRDelectronEfficiencies->GetXmin()) bin = fTRDelectronEfficiencies->GetFirst();
-  else if(electronEff > fTRDelectronEfficiencies->GetXmax()) bin = fTRDelectronEfficiencies->GetLast();
-  else bin = fTRDelectronEfficiencies->FindBin(electronEff);
- TObjString keyname = Form("%d", TMath::Nint(fTRDelectronEfficiencies->GetBinLowEdge(bin)* 100.));
-/*  printf("Key: %s\n", keyname.String().Data());*/
-  TH1F *thresholds = dynamic_cast<TH1F *>((dynamic_cast<TPair *>(fTRDthresholds->FindObject(&keyname)))->Value());
-/*  printf("thresholds: %p\n", thresholds);*/
-  return thresholds;
+void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
+  //
+  // return parameter set for the given efficiency bin
+  //
+  Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.3 * 6.);
+  memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
 }
index e9757d4..b2e4ee4 100644 (file)
  #include "AliHFEpidBase.h"
  #endif
 
-class TAxis;
-class TH1F;
-class TMap;
-class TString;
 class AliVParticle;
 
 class AliHFEpidTRD : public AliHFEpidBase{
@@ -31,6 +27,9 @@ class AliHFEpidTRD : public AliHFEpidBase{
       kLQ = 0,
       kNN = 1
     } PIDMethodTRD_t;
+    enum{
+      kThreshParams = 24
+    };
     AliHFEpidTRD(const Char_t *name);
     AliHFEpidTRD(const AliHFEpidTRD &ref);
     AliHFEpidTRD& operator=(const AliHFEpidTRD &ref);
@@ -40,22 +39,15 @@ class AliHFEpidTRD : public AliHFEpidBase{
     virtual Int_t IsSelected(AliVParticle *track);
     virtual Bool_t HasQAhistos() const { return kFALSE; };
 
-    void LoadTRDthresholds();
-
     void SetPIDMethod(PIDMethodTRD_t method) { fPIDMethod = method; };
-    void SetThresholdFile(Char_t *thresholdFile) { fThresholdFile = thresholdFile; };
-
   protected:
     void Copy(TObject &ref) const;
-    TH1F *GetTRDthresholds(Float_t electronEff);
-    void SetTRDthresholds(TH1F *thresholds, Float_t electronEff);
-
+    Double_t GetTRDthresholds(Double_t electronEff, Double_t p);
+    void InitParameters();
+    void GetParameters(Double_t electronEff, Double_t *parameters);
   private:
-    TString fThresholdFile;                                 // Threshold file name
     PIDMethodTRD_t fPIDMethod;                              // PID Method: 2D Likelihood or Neural Network
-    TMap *fTRDthresholds;                                   //! TRD Thresholds
-    TAxis *fTRDelectronEfficiencies;                        //! Electron Efficiencies corresponding to reference histos
-
+    Double_t fThreshParams[kThreshParams];                  // Threshold parametrisation
   ClassDef(AliHFEpidTRD, 1)     // TRD electron ID class
 };
 
index c01be6e..ef255c4 100644 (file)
@@ -43,11 +43,7 @@ class AliHFEpriVtx : public TObject {
                 AliHFEpriVtx &operator=(const AliHFEpriVtx &); // assignment operator
                 virtual ~AliHFEpriVtx();
 
-
-
-        protected:
                 void CreateHistograms(TString hnopt=""); // create histograms
-
                 void Init();
                 void SetEvent(AliESDEvent* ESD){fESD1=ESD;}; // set ESD pointer
                 void SetStack(AliStack* stack){fStack=stack;} // set stack pointer
@@ -87,7 +83,7 @@ class AliHFEpriVtx : public TObject {
                 TH2F *fDiffDCAvsNt; // histogram to fill DCA difference as a function of pT
 
 
-               ClassDef(AliHFEpriVtx,0); // primary vertex QA for HFE
+        ClassDef(AliHFEpriVtx,0);
 };
 
 #endif
index e06460c..85d33d8 100644 (file)
@@ -15,6 +15,8 @@
 /**************************************************************************
  *                                                                        *
  * Secondary vertexing construction Class                                 *
+ *  Construct secondary vertex from Beauty hadron with electron and       *
+ *  hadrons, then apply selection criteria                                *
  *                                                                        *
  * Authors:                                                               *
  *   MinJung Kweon <minjung@physi.uni-heidelberg.de>                      *
@@ -52,6 +54,7 @@ AliHFEsecVtx::AliHFEsecVtx():
         ,finvmass(-1)
         ,finvmassSigma(-1)
         ,fBTagged(0)
+        ,fBElec(0)
 { 
         // Default constructor
 
@@ -73,6 +76,7 @@ AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
         ,finvmass(p.finvmass)
         ,finvmassSigma(p.finvmassSigma)
         ,fBTagged(p.fBTagged)
+        ,fBElec(p.fBElec)
 { 
         // Copy constructor
 }
@@ -99,6 +103,8 @@ AliHFEsecVtx::~AliHFEsecVtx()
 void AliHFEsecVtx::Init()
 {
 
+        // set pdg code and index
+
         fNparents = 7;
 
         fParentSelect[0][0] =  411;
@@ -161,26 +167,27 @@ void AliHFEsecVtx::Init()
         fia[3][2][0] = 2;
         fia[3][2][1] = 3;
 
-
-  //fid[4][3] = {0,1,2, 0,1,3, 0,2,3, 1,2,3};
-  //fia[4][3][2] = {0,1, 0,2, 1,2,  0,1, 0,3, 1,3,  0,2, 0,3, 2,3,  1,2, 1,3, 2,3};
-
 } 
 
 //__________________________________________
 void AliHFEsecVtx::ResetTagVar()
 {
+        // reset tag variables
+
         fdistTwoSecVtx = -1;
         fcosPhi = -1;
         fsignedLxy = -1;
         finvmass = -1;
         finvmassSigma = -1;
         fBTagged = kFALSE;
+        fBElec = kFALSE;
 }
 
 //__________________________________________
 void AliHFEsecVtx::InitAnaPair()
 {
+        // initialize pair tagging variables
+
         fPairTagged = 0;
         for (Int_t i=0; i<20; i++){
                 fpairedTrackID[i] = -1;
@@ -196,16 +203,16 @@ void AliHFEsecVtx::CreateHistograms(TString hnopt)
 { 
         // create histograms
 
-        fkSourceLabel[fkAll]="all";
-        fkSourceLabel[fkDirectCharm]="directCharm";
-        fkSourceLabel[fkDirectBeauty]="directBeauty";
-        fkSourceLabel[fkBeautyCharm]="beauty2charm";
-        fkSourceLabel[fkGamma]="gamma";
-        fkSourceLabel[fkPi0]="pi0";
-        fkSourceLabel[fkElse]="others";
-        fkSourceLabel[fkBeautyGamma]="beauty22gamma";
-        fkSourceLabel[fkBeautyPi0]="beauty22pi0";
-        fkSourceLabel[fkBeautyElse]="beauty22others";
+        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;
@@ -219,12 +226,16 @@ void AliHFEsecVtx::CreateHistograms(TString hnopt)
            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);
 
         }
 
@@ -237,11 +248,53 @@ void AliHFEsecVtx::CreateHistograms(TString hnopt)
         hname=hnopt+"pt_wrongtaggedelec";
         fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,150,0,30);
 
+        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);
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourcePart, Int_t index2)
+void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t index2)
 {
+        // calculate e-h pair characteristics and tag pair 
+
+        Int_t sourcePart = PairCode(track1,track2);
 
         // get KF particle input pid
         Int_t pdg1 = GetMCPID(track1);
@@ -291,7 +344,8 @@ void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourc
         Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 
         // opening angle between two particles in XY plane
-        Double_t cosphi = TMath::Cos(kfTrack1.GetAngleXY(kfTrack2));
+        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);
@@ -302,32 +356,35 @@ void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourc
         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);
 
-
-        // pair cuts 
-        if( kfchi2 >2. ) return;
+             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);
-
-        if ( signedLxy > 0.05 && cosphi > 0.5) {
-          fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
-        }
+        fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR));
 
         // pair cuts 
-        if (signedLxy < 0.05) return;
-        if (cosphi < 0.0) return;
+        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
-        fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
+        if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
 
 
         // pair tagging condition
-        if ( signedLxy > 0.0 && cosphi > 0) {
+        if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut
         //if ( signedLxy > 0.06 && cosphi > 0) {
           fpairedTrackID[fPairTagged] = index2;
           fpairedChi2[fPairTagged] = kfchi2;
@@ -341,9 +398,20 @@ void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourc
 //_______________________________________________________________________________________________
 void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
 {
+        // run secondary vertexing algorithm and do tagging
 
         ResetTagVar();
 
+        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;
+        }
+
+
         if (fPairTagged >= 4) {
           FindSECVTXCandid4Tracks(track);         
         }
@@ -357,19 +425,10 @@ void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
           ApplyPairTagCut();      
         }
 
-        
-        Int_t imclabel = TMath::Abs(track->GetLabel());
-        if(imclabel<0) return;
-        TParticle* mcpart = fStack->Particle(imclabel);
-        Int_t esource = GetElectronSource(imclabel);
-        if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){
-                fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
-        }
-
 
         if (fBTagged) {
                 fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
-                if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){
+                if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
                         fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
                 }
                 else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
@@ -380,22 +439,34 @@ void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
 //_______________________________________________________________________________________________
 void AliHFEsecVtx::ApplyPairTagCut()
 {
+        // apply tagging cut for e-h pair
 
+        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.5) return;
+        if (fpairedSignedLxy[0] < 0.05) return;
 
         fBTagged = kTRUE;
-                
 }
 
 //_______________________________________________________________________________________________
 Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id) 
 {
+        // 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.5) return kFALSE;
+        if (fpairedSignedLxy[id] < 0.05) return kFALSE;
 
         fBTagged = kTRUE;
         return kTRUE;
@@ -405,10 +476,11 @@ Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id)
 //_______________________________________________________________________________________________
 Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2) 
 {
+        // apply tagging cut for secondary vertex
 
         if (chi2 > 2.0) return kFALSE;
         if (finvmass < 1.5) return kFALSE;
-        if (fsignedLxy < 0.5) return kFALSE;
+        if (fsignedLxy < 0.05) return kFALSE;
         if (fcosPhi < 0.90) return kFALSE;
         if (fdistTwoSecVtx > 0.1) return kFALSE;
         
@@ -417,12 +489,12 @@ Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2)
 }
 
 //_______________________________________________________________________________________________
-//void AliHFEsecVtx::ApplyTagCut(Double_t chi2, AliESDtrack *track, AliESDtrack *htrack1, AliESDtrack *htrack2)
-void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2)
+void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2)
 {
+        // apply tagging cut for e-h pair of given indexed hadron
 
         if(!ApplyTagCut(chi2)){
-                if(!ApplyPairTagCut(0)) ApplyPairTagCut(1);
+                if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2);
         }
         
 }
@@ -430,6 +502,7 @@ void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2)
 //_______________________________________________________________________________________________
 void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track) 
 {
+        // find secondary vertex for >= 4 e-h pairs
 
         // sort pair in increasing order (kFALSE-increasing order)
         Int_t index[20];
@@ -466,14 +539,21 @@ void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track)
         // 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]]);
 
-        ApplyVtxTagCut(sevchi2[indexB[0]]);
-        //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
+        if (fBElec){
+                fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+        }
+        else {
+                fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+        }
+
+        ApplyVtxTagCut(sevchi2[indexB[0]],index[fia[indexA[0]][indexB[0]][0]],index[fia[indexA[0]][indexB[0]][1]]);
 
 }
 
 //_______________________________________________________________________________________________
 void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track) 
 {
+        // find secondary vertex for 3 e-h pairs
 
         // sort pair in increasing order (kFALSE-increasing order)
         Int_t indexA[1] = { 0 };
@@ -498,13 +578,20 @@ void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track)
         // 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]]);
 
-        ApplyVtxTagCut(sevchi2[indexB[0]]);
-        //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
+        if (fBElec){
+                fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+        }
+        else {
+                fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+        }
+
+        ApplyVtxTagCut(sevchi2[indexB[0]],fia[indexA[0]][indexB[0]][0],fia[indexA[0]][indexB[0]][1]);
 }
 
 //_______________________________________________________________________________________________
 void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track) 
 {
+        // find secondary vertex for 2 e-h pairs
 
         Double_t sevchi2[1];
         AliESDtrack *htrack[2];
@@ -518,13 +605,20 @@ void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track)
         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
         CalcSECVTXProperty(track,htrack[0],htrack[1]);
 
-        ApplyVtxTagCut(sevchi2[0]);
-        //ApplyTagCut(sevchi2[0],track,htrack[0],htrack[1]);
+        if (fBElec){
+                fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[0]);
+        }
+        else {
+                fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[0]);
+        }
+
+        ApplyVtxTagCut(sevchi2[0],0,1);
 }
 
 //_______________________________________________________________________________________________
 void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
 {
+        // calculate secondary vertex properties
 
         Int_t pdg1 = GetMCPID(track1);
         Int_t pdg2 = GetMCPID(track2);
@@ -568,11 +662,26 @@ void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2,
 
         // 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);
+       }
+
 }
 
 //_______________________________________________________________________________________________
 Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track) 
 {
+        // return mc pid
 
         Int_t label = TMath::Abs(track->GetLabel());
         TParticle* mcpart = fStack->Particle(label);
@@ -585,6 +694,7 @@ Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
 //_______________________________________________________________________________________________
 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
 {
+        // return 4 track secondary vertex chi2
 
         Int_t pdg1 = GetMCPID(track1);
         Int_t pdg2 = GetMCPID(track2);
@@ -605,6 +715,7 @@ Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, A
 //_______________________________________________________________________________________________
 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
 {
+        // return 3 track secondary vertex chi2
 
         Int_t pdg1 = GetMCPID(track1);
         Int_t pdg2 = GetMCPID(track2);
@@ -623,6 +734,7 @@ Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, A
 //_______________________________________________________________________________________________
 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
 {
+        // return 2 track secondary vertex chi2
 
         Int_t pdg1 = GetMCPID(track1);
         Int_t pdg2 = GetMCPID(track2);
@@ -723,38 +835,38 @@ Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
 {
 
         //           
-        // return pair code which is predefind as:
-        //  fkDirectCharm, fkDirectBeauty, fkBeautyCharm, fkGamma, fkPi0, fkElse, fkBeautyGamma, fkBeautyPi0, fkBeautyElse
+        // 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 = fkElse;
+        Int_t sourcePart = kElse;
 
         if (pairOriginsPDG < 0) {
-                sourcePart = fkBeautyElse;
+                sourcePart = kBeautyElse;
         }
         for (Int_t i=0; i<fNparents; i++){
                 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
-                        if (pairOriginsPDG>0) sourcePart = fkDirectCharm;
+                        if (pairOriginsPDG>0) sourcePart = kDirectCharm;
                         if (pairOriginsPDG<0) {
-                                sourcePart = fkBeautyCharm;
+                                sourcePart = kBeautyCharm;
                         }
                 }
                 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
                         if (pairOriginsPDG>0) {
-                                sourcePart = fkDirectBeauty;
+                                sourcePart = kDirectBeauty;
                         }
-                        if (pairOriginsPDG<0)  return fkElse;
+                        if (pairOriginsPDG<0)  return kElse;
                 }
         }
-        if (pairOriginsPDG == 22) sourcePart = fkGamma;
+        if (pairOriginsPDG == 22) sourcePart = kGamma;
         if (pairOriginsPDG == -22) {
-                sourcePart = fkBeautyGamma;
+                sourcePart = kBeautyGamma;
         }
-        if (pairOriginsPDG == 111) sourcePart = fkPi0;
+        if (pairOriginsPDG == 111) sourcePart = kPi0;
         if (pairOriginsPDG == -111) {
-                sourcePart = fkBeautyPi0;
+                sourcePart = kBeautyPi0;
         }
 
         return sourcePart;
@@ -762,6 +874,188 @@ Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
 }
 
 //_______________________________________________________________________________________________
+Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack) 
+{
+
+    // 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 i=0; i<fNparents; i++){
+                if (abs(grandMaPDG)==fParentSelect[1][i]){
+
+                  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 i=0; i<fNparents; i++){
+                if (abs(grandMaPDG)==fParentSelect[1][i]){
+                  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 i=0; i<fNparents; i++){
+                if (abs(grandMaPDG)==fParentSelect[1][i]){
+                  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 i=0; i<fNparents; i++){
+                if (abs(grandMaPDG)==fParentSelect[1][i]){
+                  origin = kBeautyElse;
+                  return origin;
+                }
+             }
+
+             partMother = grandMa;
+          } // end of iteration 
+
+    }
+
+    return origin;
+
+}
+
+/*
+//_______________________________________________________________________________________________
 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) 
 {
 
@@ -782,16 +1076,16 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
         }
 
         Int_t origin = -1;
-        Bool_t IsFinalOpenCharm = kFALSE;
+        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.) == fkBeauty || int(abs(maPdgcode)/1000.) == fkBeauty ) {
+        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 = fkDirectBeauty;
+                                origin = kDirectBeauty;
                                 return origin;
                         }
                         else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
@@ -799,20 +1093,20 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
         } // end of if
 
         //charm --------------------------
-        else if ( int(abs(maPdgcode)/100.) == fkCharm || int(abs(maPdgcode)/1000.) == fkCharm ) {
+        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;
+                                isFinalOpenCharm = kTRUE;
                 }
-                if (!IsFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms   
+                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 = fkDirectCharm;
+                                origin = kDirectCharm;
                                 return origin;
                         }
                         if (jLabel < 0){ // safety protection even though not really necessary here
@@ -826,7 +1120,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 
                         for (Int_t j=0; j<fNparents; j++){
                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                origin = fkBeautyCharm;
+                                origin = kBeautyCharm;
                                 return origin;
                                 }
                         } 
@@ -837,14 +1131,14 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 
         //gamma --------------------------
         else if ( abs(maPdgcode) == 22 ) {
-                origin = fkGamma;
+                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 = fkGamma;
+                                origin = kGamma;
                                 return origin;
                         }
                         if (jLabel < 0){ // safety protection
@@ -858,7 +1152,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 
                         for (Int_t j=0; j<fNparents; j++){
                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                        origin = fkBeautyGamma;
+                                        origin = kBeautyGamma;
                                         return origin;
                                 }
                         }
@@ -871,14 +1165,14 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 
         //pi0 --------------------------
         else if ( abs(maPdgcode) == 111 ) {
-                origin = fkPi0;
+                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 = fkPi0;
+                                origin = kPi0;
                                 return origin;
                         }
                         if (jLabel < 0){ // safety protection
@@ -892,7 +1186,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 
                         for (Int_t j=0; j<fNparents; j++){
                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                        origin = fkBeautyPi0;
+                                        origin = kBeautyPi0;
                                         return origin;
                                 }
                         }
@@ -906,14 +1200,14 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 
         //else --------------------------
         else {
-                origin = fkElse;
+                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 = fkElse;
+                                origin = kElse;
                                 return origin;
                         }
                         if (jLabel < 0){ // safety protection
@@ -927,7 +1221,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 
                         for (Int_t j=0; j<fNparents; j++){
                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
-                                        origin = fkBeautyElse;
+                                        origin = kBeautyElse;
                                         return origin;
                                 }
                         }
@@ -940,10 +1234,12 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
         return origin;
 
 }
+*/
 
 //_______________________________________________________________________________________________
 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track)
 {
+        // test cuts 
 
         //if (track->Pt() < 1.0) return kFALSE;
         //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
index 2edeb14..484f311 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-
 /**************************************************************************
  *                                                                        *
  * Secondary vertexing construction Class                                 *
+ *  Construct secondary vertex from Beauty hadron with electron and       *
+ *  hadrons, then apply selection criteria                                *
  *                                                                        *
  **************************************************************************/
 
 
-#ifndef _ALIHFESECVTX_H_
-#define _ALIHFESECVTX_H_
+#ifndef ALIHFESECVTX_H
+#define ALIHFESECVTX_H
 
 #ifndef ROOT_TObject
 #include <TObject.h>
@@ -43,33 +44,35 @@ class AliHFEsecVtx : public TObject {
                 AliHFEsecVtx &operator=(const AliHFEsecVtx &); // assignment operator
                 virtual ~AliHFEsecVtx();
 
-
-        protected:
                 void CreateHistograms(TString hnopt="");
 
-                void Init();
+                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 SetEvent(AliESDEvent* ESD){fESD1=ESD;}; // set ESD pointer
-                void SetStack(AliStack* stack){fStack=stack;} // set stack pointer
-                void AnaPair(AliESDtrack* ESDtrack1, AliESDtrack* ESDtrack2, Int_t sourcePart, Int_t index2); // do e-h analysis
+                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 GetElectronSource(Int_t mclabel); // return origin of the electron
+
+                Bool_t SingleTrackCut(AliESDtrack* track1); // single track cut
+
+        protected:
+
+                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); // apply tagging condition for 3 track secondary vertex
-                //void ApplyTagCut(Double_t chi2, AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // apply tagging condition for 3 track 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
 
                 Bool_t ApplyPairTagCut(Int_t id); // apply strict tagging condition for 1 pair secondary vertex
                 Bool_t ApplyTagCut(Double_t chi2); // apply tagging condition
-                Bool_t SingleTrackCut(AliESDtrack* track1); // single track cut
 
-                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 GetElectronSource(Int_t mclabel); // return origin of the electron
 
                 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
@@ -81,22 +84,24 @@ class AliHFEsecVtx : public TObject {
                 AliESDEvent* fESD1; // ESD pointer             
                 AliStack* fStack; // stack pointer              
 
-                Int_t fParentSelect[2][7]; //
-                Int_t fNparents; // 
+                Int_t fParentSelect[2][7]; // heavy hadron species
+                Int_t fNparents; // number of heavy hadrons to be considered
 
-                TString fkSourceLabel[10]; //
+                TString fkSourceLabel[10]; // electron source label
 
-                enum {fkAll, fkDirectCharm, fkDirectBeauty, fkBeautyCharm, fkGamma, fkPi0, fkElse, fkBeautyGamma, fkBeautyPi0, fkBeautyElse};
-                enum {fkCharm=4, fkBeauty=5};
+                enum {kAll, kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse};
+                enum {kCharm=4, kBeauty=5};
 
                 struct histspair{
                         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 *fCosOpenAngle; // histogram to store pair opening angle
+                       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
                 };
 
                 struct histstag{
@@ -104,29 +109,46 @@ class AliHFEsecVtx : public TObject {
                         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
                 };
 
-                histspair fHistPair[10]; // 
-                histstag fHistTagged; //
+                histspair fHistPair[10]; // struct of above given histogram for different electron sources
+                histstag fHistTagged; // struct of above given histogram
 
-                Int_t fPairTagged; //
-                Int_t fpairedTrackID[20]; //
-                Double_t fpairedChi2[20]; //
-                Double_t fpairedInvMass[20]; //
-                Double_t fpairedSignedLxy[20]; //
+                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 fid[4][3]; //
-                Int_t fia[4][3][2]; //
+                Int_t fid[4][3]; // index to store sorting result
+                Int_t fia[4][3][2]; // index to store sorting result
 
-                Double_t fdistTwoSecVtx; // 
-                Double_t fcosPhi; // 
+                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
 
                 Bool_t fBTagged; // if b tagged, set true
+                Bool_t fBElec; // if set true for b electron, set true
 
-               ClassDef(AliHFEsecVtx,0); // secondary vertex for electrons with AliKFParticle package
+        ClassDef(AliHFEsecVtx,0);
 };
 
 #endif