PWG3/hfe/AliAnalysisTaskHFE.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // The analysis task:
  17 // Filling an AliCFContainer with the quantities pt, eta and phi
  18 // for tracks which survivied the particle cuts (MC resp. ESD tracks)
  19 // Track selection is done using the AliHFE package
  20 //
  21 // Author:
  22 //  Raphaelle Bailhache <R.Bailhache@gsi.de>
  23 //  Markus Fasel <M.Fasel@gsi.de>
  24 //  Matus Kalisky <matus.kalisky@cern.ch>
  25 //  MinJung Kweon <minjung@physi.uni-heidelberg.de>
  26 //
  27 #include <TAxis.h>
  28 #include <TCanvas.h>
  29 #include <TChain.h>
  30 #include <TDirectory.h>
  31 #include <TFile.h>
  32 #include <TH1D.h>
  33 #include <TH1F.h>
  34 #include <TH1I.h>
  35 #include <TH2F.h>
  36 #include <TIterator.h>
  37 #include <TList.h>
  38 #include <TLegend.h>
  39 #include <TMath.h>
  40 #include <TObjArray.h>
  41 #include <TParticle.h>
  42 #include <TProfile.h>
  43 #include <TString.h>
  44 #include <TTree.h>
  45
  46 #include "AliCFContainer.h"
  47 #include "AliCFManager.h"
  48 #include "AliCFEffGrid.h"
  49 #include "AliESDEvent.h"
  50 #include "AliESDInputHandler.h"
  51 #include "AliESDtrack.h"
  52 #include "AliESDtrackCuts.h"
  53 #include "AliLog.h"
  54 #include "AliAnalysisManager.h"
  55 #include "AliMCEvent.h"
  56 #include "AliMCEventHandler.h"
  57 #include "AliMCParticle.h"
  58 #include "AliPID.h"
  59 #include "AliStack.h"
  60
  61 #include "AliHFEpid.h"
  62 #include "AliHFEcuts.h"
  63 #include "AliHFEmcQA.h"
  64 #include "AliHFEsecVtx.h"
  65 #include "AliAnalysisTaskHFE.h"
  66
  67 //____________________________________________________________
  68 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  69   AliAnalysisTask("PID efficiency Analysis", "")
  70   , fQAlevel(0)
  71   , fPIDdetectors("")
  72   , fPIDstrategy(0)
  73   , fESD(0x0)
  74   , fMC(0x0)
  75   , fCFM(0x0)
  76   , fCorrelation(0x0)
  77   , fPIDperformance(0x0)
  78   , fPID(0x0)
  79   , fCuts(0x0)
  80   , fSecVtx(0x0)
  81   , fMCQA(0x0)
  82   , fNEvents(0x0)
  83   , fNElectronTracksEvent(0x0)
  84   , fQA(0x0)
  85   , fOutput(0x0)
  86   , fHistMCQA(0x0)
  87   , fHistSECVTX(0x0)
  88 {
  89   //
  90   // Dummy constructor
  91   //
  92   DefineInput(0, TChain::Class());
  93   DefineOutput(0, TH1I::Class());
  94   DefineOutput(1, TList::Class());
  95   DefineOutput(2, TList::Class());
  96
  97   // Initialize cuts
  98   fPID = new AliHFEpid;
  99
 100   SetHasMCData();
 101 }
 102
 103 //____________________________________________________________
 104 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
 105   AliAnalysisTask(name, "")
 106   , fQAlevel(0)
 107   , fPIDdetectors("")
 108   , fPIDstrategy(0)
 109   , fESD(0x0)
 110   , fMC(0x0)
 111   , fCFM(0x0)
 112   , fCorrelation(0x0)
 113   , fPIDperformance(0x0)
 114   , fPID(0x0)
 115   , fCuts(0x0)
 116   , fSecVtx(0x0)
 117   , fMCQA(0x0)
 118   , fNEvents(0x0)
 119   , fNElectronTracksEvent(0x0)
 120   , fQA(0x0)
 121   , fOutput(0x0)
 122   , fHistMCQA(0x0)
 123   , fHistSECVTX(0x0)
 124 {
 125   //
 126   // Default constructor
 127   //
 128   DefineInput(0, TChain::Class());
 129   DefineOutput(0, TH1I::Class());
 130   DefineOutput(1, TList::Class());
 131   DefineOutput(2, TList::Class());
 132
 133   // Initialize cuts
 134   fPID = new AliHFEpid;
 135
 136   SetHasMCData();
 137 }
 138
 139 //____________________________________________________________
 140 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
 141   AliAnalysisTask(ref)
 142   , fQAlevel(ref.fQAlevel)
 143   , fPIDdetectors(ref.fPIDdetectors)
 144   , fPIDstrategy(ref.fPIDstrategy)
 145   , fESD(ref.fESD)
 146   , fMC(ref.fMC)
 147   , fCFM(ref.fCFM)
 148   , fCorrelation(ref.fCorrelation)
 149   , fPIDperformance(ref.fPIDperformance)
 150   , fPID(ref.fPID)
 151   , fCuts(ref.fCuts)
 152   , fSecVtx(ref.fSecVtx)
 153   , fMCQA(ref.fMCQA)
 154   , fNEvents(ref.fNEvents)
 155   , fNElectronTracksEvent(ref.fNElectronTracksEvent)
 156   , fQA(ref.fQA)
 157   , fOutput(ref.fOutput)
 158   , fHistMCQA(ref.fHistMCQA)
 159   , fHistSECVTX(ref.fHistSECVTX)
 160 {
 161   //
 162   // Copy Constructor
 163   //
 164 }
 165
 166 //____________________________________________________________
 167 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
 168   //
 169   // Assignment operator
 170   //
 171   if(this == &ref) return *this;
 172   AliAnalysisTask::operator=(ref);
 173   fQAlevel = ref.fQAlevel;
 174   fPIDdetectors = ref.fPIDdetectors;
 175   fPIDstrategy = ref.fPIDstrategy;
 176   fESD = ref.fESD;
 177   fMC = ref.fMC;
 178   fCFM = ref.fCFM;
 179   fCorrelation = ref.fCorrelation;
 180   fPIDperformance = ref.fPIDperformance;
 181   fPID = ref.fPID;
 182   fCuts = ref.fCuts;
 183   fSecVtx = ref.fSecVtx;
 184   fMCQA = ref.fMCQA;
 185   fNEvents = ref.fNEvents;
 186   fNElectronTracksEvent = ref.fNElectronTracksEvent;
 187   fQA = ref.fQA;
 188   fOutput = ref.fOutput;
 189   fHistMCQA = ref.fHistMCQA;
 190   fHistSECVTX = ref.fHistSECVTX;
 191   return *this;
 192 }
 193
 194 //____________________________________________________________
 195 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
 196   //
 197   // Destructor
 198   //
 199   if(fESD) delete fESD;
 200   if(fMC) delete fMC;
 201   if(fPID) delete fPID;
 202   if(fQA){
 203     fQA->Clear();
 204     delete fQA;
 205   }
 206   if(fOutput){
 207     fOutput->Clear();
 208     delete fOutput;
 209   }
 210   if(fHistMCQA){
 211     fHistMCQA->Clear();
 212     delete fHistMCQA;
 213   }
 214   if(fHistSECVTX){
 215     fHistSECVTX->Clear();
 216     delete fHistSECVTX;
 217   }
 218   if(fSecVtx) delete fSecVtx;
 219   if(fMCQA) delete fMCQA;
 220   if(fNEvents) delete fNEvents;
 221   if(fCorrelation){
 222     fCorrelation->Clear();
 223     delete fCorrelation;
 224   }
 225   if(fPIDperformance) delete fPIDperformance;
 226 }
 227
 228 //____________________________________________________________
 229 void AliAnalysisTaskHFE::ConnectInputData(Option_t *){
 230   //
 231   // Connecting the input
 232   // Called once per worker
 233   //
 234 /*  TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
 235   if(!esdchain){
 236     AliError("ESD chain empty");
 237     return;
 238   } else {
 239     esdchain->SetBranchStatus("Tracks", 1);
 240   }
 241 */
 242   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 243   if(!esdH){
 244     AliError("No ESD input handler");
 245     return;
 246   } else {
 247     fESD = esdH->GetEvent();
 248   }
 249   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 250   if(!mcH){
 251     AliError("No MC truth handler");
 252     return;
 253   } else {
 254     fMC = mcH->MCEvent();
 255   }
 256 }
 257
 258 //____________________________________________________________
 259 void AliAnalysisTaskHFE::CreateOutputObjects(){
 260   //
 261   // Creating output container and output objects
 262   // Here we also Initialize the correction framework container and
 263   // the objects for
 264   // - PID
 265   // - MC QA
 266   // - SecVtx
 267   // QA histograms are created if requested
 268   // Called once per worker
 269   //
 270   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
 271   fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
 272   // First Step: TRD alone
 273   if(!fQA) fQA = new TList;
 274   fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
 275   fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
 276   fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
 277   fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
 278   fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
 279   fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
 280   fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
 281
 282   if(!fOutput) fOutput = new TList;
 283   // Initialize correction Framework and Cuts
 284   fCFM = new AliCFManager;
 285   MakeParticleContainer();
 286   // Temporary fix: Initialize particle cuts with 0x0
 287   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
 288     fCFM->SetParticleCutsList(istep, 0x0);
 289   if(!fCuts){
 290     AliWarning("Cuts not available. Default cuts will be used");
 291     fCuts = new AliHFEcuts;
 292     fCuts->CreateStandardCuts();
 293   }
 294   fCuts->Initialize(fCFM);
 295   if(fCuts->IsInDebugMode()) fQA->Add(fCuts->GetQAhistograms());
 296
 297   // add output objects to the List
 298   fOutput->AddAt(fCFM->GetParticleContainer(), 0);
 299   fOutput->AddAt(fCorrelation, 1);
 300   fOutput->AddAt(fPIDperformance, 2);
 301   fOutput->AddAt(fNElectronTracksEvent, 3);
 302
 303   // Initialize PID
 304   if(IsQAOn(kPIDqa)){
 305     AliInfo("PID QA switched on");
 306     //fPID->SetDebugLevel(2);
 307     fPID->SetQAOn();
 308     fQA->Add(fPID->GetQAhistograms());
 309   }
 310   fPID->SetHasMCData(HasMCData());
 311   if(!fPIDdetectors.Length() && ! fPIDstrategy) AddPIDdetector("TPC");
 312   if(fPIDstrategy)
 313     fPID->InitializePID(Form("Strategy%d", fPIDstrategy));
 314   else
 315     fPID->InitializePID(fPIDdetectors.Data());     // Only restrictions to TPC allowed
 316
 317   // mcQA----------------------------------
 318   if (IsQAOn(kMCqa)) {
 319     AliInfo("MC QA on");
 320     if(!fMCQA) fMCQA = new AliHFEmcQA;
 321     if(!fHistMCQA) fHistMCQA = new TList();
 322     fHistMCQA->SetName("MCqa");
 323     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
 324     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
 325     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm
 326     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
 327     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm
 328     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
 329     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_");        // create histograms for charm
 330     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_");       // create histograms for beauty
 331     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_");     // create histograms for charm
 332     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_");    // create histograms for beauty
 333     TIter next(gDirectory->GetList());
 334     TObject *obj;
 335     int counter = 0;
 336     TString objname;
 337     while ((obj = next.Next())) {
 338       objname = obj->GetName();
 339       TObjArray *toks = objname.Tokenize("_");
 340       if (toks->GetEntriesFast()){
 341         TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
 342         if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj, counter++);
 343       }
 344     }
 345     fQA->Add(fHistMCQA);
 346   }
 347
 348   // secvtx----------------------------------
 349   if (IsSecVtxOn()) {
 350     AliInfo("Secondary Vertex Analysis on");
 351     fSecVtx = new AliHFEsecVtx;
 352
 353     if(!fHistSECVTX) fHistSECVTX = new TList();
 354     fHistSECVTX->SetName("SecVtx");
 355     fSecVtx->CreateHistograms("secvtx_");
 356     TIter next(gDirectory->GetList());
 357     TObject *obj;
 358     int counter = 0;
 359     TString objname;
 360     while ((obj = next.Next())) {
 361       objname = obj->GetName();
 362       TObjArray *toks = objname.Tokenize("_");
 363       if (toks->GetEntriesFast()){
 364         TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
 365         if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj, counter++);
 366       }
 367     }
 368     fOutput->Add(fHistSECVTX);
 369   }
 370 }
 371
 372 //____________________________________________________________
 373 void AliAnalysisTaskHFE::Exec(Option_t *){
 374   //
 375   // Run the analysis
 376   //
 377   if(!fESD){
 378     AliError("No ESD Event");
 379     return;
 380   }
 381   if(!fMC){
 382     AliError("No MC Event");
 383     return;
 384   }
 385   if(!fCuts){
 386     AliError("HFE cuts not available");
 387     return;
 388   }
 389
 390   //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
 391   Double_t container[6];
 392   // container for the output THnSparse
 393   Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
 394
 395   // Loop over the Monte Carlo tracks to see whether we have overlooked any track
 396   AliMCParticle *mctrack = 0x0;
 397   Int_t nElectrons = 0;
 398
 399   if (IsSecVtxOn()) {
 400     fSecVtx->SetEvent(fESD);
 401     fSecVtx->SetStack(fMC->Stack());
 402   }
 403
 404   // run mc QA ------------------------------------------------
 405   if (IsQAOn(kMCqa)) {
 406     AliDebug(2, "Running MC QA");
 407
 408     fMCQA->SetStack(fMC->Stack());
 409     fMCQA->Init();
 410
 411     Int_t nMCTracks = fMC->Stack()->GetNtrack();
 412
 413     // loop over all tracks for decayed electrons
 414     for (Int_t igen = 0; igen < nMCTracks; igen++){
 415       TParticle* mcpart = fMC->Stack()->Particle(igen);
 416       fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
 417       fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
 418       fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
 419       fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
 420       fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0); // no accept cut
 421       fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 422       if (TMath::Abs(mcpart->Eta()) < 0.9) {
 423         fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 424         fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 425       }
 426       if (TMath::Abs(GetRapidity(mcpart)) < 0.5) {
 427         fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 428         fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 429       }
 430     }
 431     fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
 432     fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 433
 434   } // end of MC QA loop
 435   // -----------------------------------------------------------------
 436
 437   //
 438   // Loop MC
 439   //
 440   fCuts->SetEventInfo(fMC);
 441   for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
 442     mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(imc));
 443
 444     container[0] = mctrack->Pt();
 445     container[1] = mctrack->Eta();
 446     container[2] = mctrack->Phi();
 447
 448     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
 449     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
 450     if(IsSignalElectron(mctrack)) fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCsignal);
 451     (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
 452     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
 453     // find the label in the vector
 454     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
 455     nElectrons++;
 456   }
 457   (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
 458
 459   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 460
 461
 462   //
 463   // Loop ESD
 464   //
 465   Int_t nElectronCandidates = 0;
 466   AliESDtrack *track = 0x0, *htrack = 0x0;
 467   Int_t pid = 0;
 468   // For double counted tracks
 469   LabelContainer cont(fESD->GetNumberOfTracks());
 470   Bool_t alreadyseen = kFALSE;
 471
 472   Bool_t signal = kTRUE;
 473
 474   fCuts->SetEventInfo(fESD);
 475   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 476
 477     track = fESD->GetTrack(itrack);
 478
 479     container[0] = track->Pt();
 480     container[1] = track->Eta();
 481     container[2] = track->Phi();
 482
 483     dataE[0] = track->Pt();
 484     dataE[1] = track->Eta();
 485     dataE[2] = track->Phi();
 486     dataE[3] = -1;
 487     dataE[4] = -1;
 488
 489     signal = kTRUE;
 490
 491     // RecKine: ITSTPC cuts
 492     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 493
 494     // Check if it is electrons near the vertex
 495     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
 496     TParticle* mctrack4QA = fMC->Stack()->Particle(TMath::Abs(track->GetLabel()));
 497
 498     container[3] = mctrack->Pt();
 499     container[4] = mctrack->Eta();
 500     container[5] = mctrack->Phi();
 501
 502     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
 503
 504     if(signal) {
 505       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));
 506       cont.Append(TMath::Abs(track->GetLabel()));
 507
 508       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecKineITSTPC);
 509       fCFM->GetParticleContainer()->Fill(&container[0], (AliHFEcuts::kStepRecKineITSTPC + 2*(AliHFEcuts::kNcutESDSteps + 1)));
 510       if(alreadyseen) {
 511         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITSTPC + (AliHFEcuts::kNcutESDSteps + 1)));
 512       }
 513     }
 514
 515     // Check TRD criterions (outside the correction framework)
 516     if(track->GetTRDncls()){
 517       (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
 518       (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha());    // Check the acceptance without tight cuts
 519       (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
 520       (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
 521     }
 522
 523
 524     // RecPrim
 525     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
 526     if(signal) {
 527       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1));
 528       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim);
 529       if(alreadyseen) {
 530         fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1);
 531       }
 532     }
 533
 534
 535     // HFEcuts: ITS layers cuts
 536     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 537     if(signal) {
 538       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1));
 539       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS);
 540       if(alreadyseen) {
 541         fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1);
 542       }
 543     }
 544
 545     // HFEcuts: Nb of tracklets TRD0
 546     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
 547     if(signal) {
 548       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1));
 549       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD);
 550       if(alreadyseen) {
 551         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + (AliHFEcuts::kNcutESDSteps + 1)));
 552       }
 553       // dimensions 3&4&5 : pt,eta,phi (MC)
 554       ((THnSparseF *)fCorrelation->At(0))->Fill(container);
 555     }
 556
 557    if (IsQAOn(kMCqa)) {
 558      // mc qa for after the reconstruction cuts
 559      AliDebug(2, "Running MC QA");
 560      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 3);  // charm
 561      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 3); // beauty
 562     }
 563
 564     // track accepted, do PID
 565     AliHFEpidObject hfetrack;
 566     hfetrack.fAnalysisType = AliHFEpidObject::kESDanalysis;
 567     hfetrack.fRecTrack = track;
 568     if(HasMCData()) hfetrack.fMCtrack = mctrack;
 569     if(!fPID->IsSelected(&hfetrack)) continue;
 570     nElectronCandidates++;
 571
 572    if (IsQAOn(kMCqa)) {
 573      // mc qa for after the reconstruction and pid cuts
 574      AliDebug(2, "Running MC QA");
 575      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 4);  // charm
 576      fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 4); // beauty
 577     }
 578
 579     // Fill Containers
 580     if(signal) {
 581       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD +1 + 2*(AliHFEcuts::kNcutESDSteps + 1));
 582       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 1);
 583       if(alreadyseen) {
 584         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + (AliHFEcuts::kNcutESDSteps + 1)));
 585       }
 586       // dimensions 3&4&5 : pt,eta,phi (MC)
 587       ((THnSparseF *)fCorrelation->At(1))->Fill(container);
 588     }
 589
 590     // Track selected: distinguish between true and fake
 591     AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
 592     if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
 593       Int_t type = IsSignalElectron(track);
 594       AliDebug(1, Form("Type: %d\n", type));
 595       if(type){
 596               dataE[4] = type; // beauty[1] or charm[2]
 597               dataE[3] = 2;  // signal electron
 598       }
 599       else{
 600               dataE[3] = 1; // not a signal electron
 601               dataE[4] = 0;
 602       }
 603       // pair analysis [mj]
 604       if (IsSecVtxOn()) {
 605         AliDebug(2, "Running Secondary Vertex Analysis");
 606         fSecVtx->InitAnaPair();
 607         for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
 608           htrack = fESD->GetTrack(jtrack);
 609           if ( itrack == jtrack ) continue;
 610           //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue;
 611           if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue;
 612           fSecVtx->AnaPair(track, htrack, jtrack);
 613         }
 614         // based on the partner of e info, you run secandary vertexing function
 615         fSecVtx->RunSECVTX(track);
 616       } // end of pair analysis
 617     }
 618     else {
 619       // Fill THnSparse with the information for Fake Electrons
 620       dataE[3] = 0;
 621       dataE[4] = 0;
 622     }
 623     // fill the performance THnSparse, if the mc origin could be defined
 624     if(dataE[3] > -1){
 625       AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4]));
 626       fPIDperformance->Fill(dataE);
 627     }
 628   }
 629
 630
 631   fNEvents->Fill(1);
 632   fNElectronTracksEvent->Fill(nElectronCandidates);
 633
 634   // Done!!!
 635   PostData(0, fNEvents);
 636   PostData(1, fOutput);
 637   PostData(2, fQA);
 638 }
 639
 640 //____________________________________________________________
 641 void AliAnalysisTaskHFE::Terminate(Option_t *){
 642   //
 643   // Terminate not implemented at the moment
 644   //
 645   if(IsRunningPostProcess()){
 646     fOutput = dynamic_cast<TList *>(GetOutputData(1));
 647     if(!fOutput){
 648       AliError("Results not available");
 649       return;
 650     }
 651     PostProcess();
 652   }
 653 }
 654
 655 //____________________________________________________________
 656 void AliAnalysisTaskHFE::Load(TString filename){
 657   //
 658   // Load Results into the task
 659   //
 660   fQA = NULL; fOutput = NULL; fNEvents = NULL;
 661   TFile *input = TFile::Open(filename.Data());
 662   if(!input || input->IsZombie()){
 663     AliError("Cannot read file");
 664     return;
 665   }
 666   TH1 *htmp = dynamic_cast<TH1I *>(input->Get("nEvents"));
 667   if(htmp)
 668     fNEvents = dynamic_cast<TH1I *>(htmp->Clone());
 669   else
 670     AliError("Event Counter histogram not found");
 671   TList *ltmp = dynamic_cast<TList *>(input->Get("Results"));
 672   if(ltmp)
 673     fOutput = dynamic_cast<TList *>(ltmp->Clone());
 674   else
 675     AliError("Output Histograms not found");
 676   ltmp = dynamic_cast<TList *>(input->Get("QA"));
 677   if(ltmp)
 678     fQA = dynamic_cast<TList *>(ltmp->Clone());
 679   else
 680     AliError("QA histograms not found");
 681   input->Close();
 682   delete input;
 683   Int_t nObjects = 0;
 684   if(fNEvents) nObjects++;
 685   if(fOutput) nObjects++;
 686   if(fQA) nObjects++;
 687   AliInfo(Form("Loaded %d Objects into task", nObjects));
 688 }
 689
 690 //____________________________________________________________
 691 void AliAnalysisTaskHFE::PostProcess(){
 692   //
 693   // Plotting Ratio histograms
 694   // + All electrons / all candidates (Purity for Electrons)
 695   // + All signal electrons / all electrons (Purity for signals)
 696   // For this the following pt-histograms have to be projected from the THnSparse
 697   // + All Electron candidates
 698   // + All Real electrons
 699   // + All Signal Electrons
 700   // + All misidentified electrons
 701   // Additionally we draw Efficiency histograms:
 702   // + Reconstructible Electrons
 703   // + Signal Electrons
 704   // + Selected Electrons
 705   // + Selected Electrons in Acceptance
 706   //
 707
 708   fPIDperformance = dynamic_cast<THnSparseF *>(fOutput->FindObject("PIDperformance"));
 709   if(!fPIDperformance){
 710     AliError("PID Performance histogram not found in the output container");
 711     return;
 712   }
 713   // Make projection
 714   // always project to pt dimension
 715   // get the histograms under our control
 716   TH1 *allCandidates = NULL, *allElectrons = NULL, *allSignals = NULL, *allFakes = NULL;
 717   allCandidates = fPIDperformance->Projection(0);
 718   allCandidates->SetName("hAllCandidates");
 719   allCandidates->SetTitle("All Candidates");
 720   allCandidates->Sumw2();
 721   Int_t firstDim3 = fPIDperformance->GetAxis(3)->GetFirst(), lastDim3 = fPIDperformance->GetAxis(3)->GetLast();
 722   fPIDperformance->GetAxis(3)->SetRange(firstDim3 + 1, lastDim3);
 723   allElectrons = fPIDperformance->Projection(0);
 724   allElectrons->Sumw2();
 725   allElectrons->SetName("hAllElectrons");
 726   allElectrons->SetTitle("All Electrons");
 727   Int_t firstDim4 = fPIDperformance->GetAxis(4)->GetFirst(), lastDim4 = fPIDperformance->GetAxis(4)->GetLast();
 728   fPIDperformance->GetAxis(4)->SetRange(firstDim4 + 1, lastDim4);
 729   allSignals = fPIDperformance->Projection(0);
 730   allSignals->Sumw2();
 731   allSignals->SetName("hAllSignals");
 732   allSignals->SetTitle("All Signal Electrons");
 733   fPIDperformance->GetAxis(4)->SetRange(firstDim4, lastDim4);       // Reset 4th axis
 734   fPIDperformance->GetAxis(3)->SetRange(firstDim3, firstDim3);  // Select fakes
 735   allFakes = fPIDperformance->Projection(0);
 736   allFakes->Sumw2();
 737   allFakes->SetName("hAllFakes");
 738   allFakes->SetTitle("All Fakes");
 739   fPIDperformance->GetAxis(3)->SetRange(firstDim3, lastDim3);       // Reset also 3rd axis
 740
 741   // Make Ratios
 742   TH1D *electronPurity = dynamic_cast<TH1D *>(allElectrons->Clone());
 743   electronPurity->Divide(allCandidates);
 744   electronPurity->SetName("electronPurity");
 745   electronPurity->SetTitle("Electron Purity");
 746   electronPurity->GetXaxis()->SetTitle("p_{T} / GeV/c");
 747   electronPurity->GetYaxis()->SetTitle("Purity / %");
 748   TH1D *signalPurity = dynamic_cast<TH1D *>(allSignals->Clone());
 749   signalPurity->Divide(allElectrons);
 750   signalPurity->SetName("signalPurity");
 751   signalPurity->SetTitle("Purity of Electrons coming from Heavy flavours");
 752   signalPurity->GetXaxis()->SetTitle("p_{T} / GeV/c");
 753   signalPurity->GetYaxis()->SetTitle("Purity / %");
 754   TH1D *fakeContamination = dynamic_cast<TH1D *>(allFakes->Clone());
 755   fakeContamination->Divide(allCandidates);
 756   fakeContamination->SetName("fakeContamination");
 757   fakeContamination->SetTitle("Contamination of misidentified hadrons");
 758   fakeContamination->GetXaxis()->SetTitle("p_{T} / GeV/c");
 759   fakeContamination->GetYaxis()->SetTitle("Purity / %");
 760
 761   // Graphics settings
 762   const Int_t nHistos = 7;
 763   TH1 *hptr[7] = {allCandidates, allElectrons, allSignals, allFakes, electronPurity, signalPurity, fakeContamination}, *histo;
 764   for(Int_t ihist = 0; ihist < nHistos; ihist++){
 765     (histo = hptr[ihist])->SetStats(kFALSE);
 766     histo->SetLineColor(kBlack);
 767     histo->SetMarkerColor(kBlue);
 768     histo->SetMarkerStyle(22);
 769   }
 770
 771   // Make percent output
 772   electronPurity->Scale(100);
 773   signalPurity->Scale(100);
 774   fakeContamination->Scale(100);
 775
 776   // Draw output
 777   TCanvas *cSpectra = new TCanvas("cSpectra","pt Spectra", 800, 600);
 778   cSpectra->Divide(2,2);
 779   TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600);
 780   cRatios->Divide(2,2);
 781   cSpectra->cd(1);
 782   gPad->SetLogy();
 783   allCandidates->GetXaxis()->SetTitle("p_{T} / GeV/c");
 784   allCandidates->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
 785   allCandidates->Draw("e");
 786   cSpectra->cd(2);
 787   gPad->SetLogy();
 788   allElectrons->GetXaxis()->SetTitle("p_{T} / GeV/c");
 789   allElectrons->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
 790   allElectrons->Draw("e");
 791   cSpectra->cd(3);
 792   gPad->SetLogy();
 793   allSignals->GetXaxis()->SetTitle("p_{T} / GeV/c");
 794   allSignals->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
 795   allSignals->Draw("e");
 796   cSpectra->cd(4);
 797   gPad->SetLogy();
 798   allFakes->GetXaxis()->SetTitle("p_{T} / GeV/c");
 799   allFakes->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
 800   allFakes->Draw("e");
 801   cRatios->cd(1);
 802   electronPurity->Draw("e");
 803   cRatios->cd(2);
 804   signalPurity->Draw("e");
 805   cRatios->cd(3);
 806   fakeContamination->Draw("e");
 807
 808   // Now draw the Efficiency
 809   // We show:
 810   // + InAcceptance / Generated
 811   // + Signal / Generated
 812   // + Selected / Generated
 813   // + Selected / InAcceptance (Reconstructible)
 814   TCanvas *cEff = new TCanvas("cEff", "Efficiency", 800, 600);
 815   cEff->Divide(2,2);
 816   AliCFContainer *effc = dynamic_cast<AliCFContainer *>(fOutput->At(0));
 817   if(!effc){
 818     AliError("Efficiency Container not found in Output Container");
 819     return;
 820   }
 821   AliCFEffGrid *effCalc = new AliCFEffGrid("effCalc", "Efficiency Calculation Grid", *effc);
 822   effCalc->CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
 823   TH1 *effReconstructibleP = effCalc->Project(0);
 824   effReconstructibleP->SetName("effReconstructibleP");
 825   effReconstructibleP->SetTitle("Efficiency of reconstructible tracks");
 826   effReconstructibleP->GetXaxis()->SetTitle("p_{T} / GeV/c");
 827   effReconstructibleP->GetYaxis()->SetTitle("Efficiency");
 828   effReconstructibleP->GetYaxis()->SetRangeUser(0.,1.);
 829   effReconstructibleP->SetMarkerStyle(22);
 830   effReconstructibleP->SetMarkerColor(kBlue);
 831   effReconstructibleP->SetLineColor(kBlack);
 832   effReconstructibleP->SetStats(kFALSE);
 833   cEff->cd(1);
 834   effReconstructibleP->Draw("e");
 835   effCalc->CalculateEfficiency(AliHFEcuts::kStepMCsignal, AliHFEcuts::kStepMCGenerated);
 836   TH1 *effSignal = effCalc->Project(0);
 837   effSignal->SetName("effSignal");
 838   effSignal->SetTitle("Efficiency of Signal Electrons");
 839   effSignal->GetXaxis()->SetTitle("p_{T} / GeV/c");
 840   effSignal->GetYaxis()->SetTitle("Efficiency");
 841   effSignal->GetYaxis()->SetRangeUser(0., 1.);
 842   effSignal->SetMarkerStyle(22);
 843   effSignal->SetMarkerColor(kBlue);
 844   effSignal->SetLineColor(kBlack);
 845   effSignal->SetStats(kFALSE);
 846   cEff->cd(2);
 847   effSignal->Draw("e");
 848   effCalc->CalculateEfficiency(AliHFEcuts::kStepHFEcutsTRD + 1, AliHFEcuts::kStepMCGenerated);
 849   TH1 *effPIDP = effCalc->Project(0);
 850   effPIDP->SetName("effPIDP");
 851   effPIDP->SetTitle("Efficiency of selected tracks");
 852   effPIDP->GetXaxis()->SetTitle("p_{T} / GeV/c");
 853   effPIDP->GetYaxis()->SetTitle("Efficiency");
 854   effPIDP->GetYaxis()->SetRangeUser(0.,1.);
 855   effPIDP->SetMarkerStyle(22);
 856   effPIDP->SetMarkerColor(kBlue);
 857   effPIDP->SetLineColor(kBlack);
 858   effPIDP->SetStats(kFALSE);
 859   cEff->cd(3);
 860   effPIDP->Draw("e");
 861   effCalc->CalculateEfficiency(AliHFEcuts::kStepHFEcutsTRD + 1, AliHFEcuts::kStepMCInAcceptance);
 862   TH1 *effPIDAcc = effCalc->Project(0);
 863   effPIDAcc->SetName("effPIDAcc");
 864   effPIDAcc->SetTitle("Efficiency of selected tracks in acceptance");
 865   effPIDAcc->GetXaxis()->SetTitle("p_{T} / GeV/c");
 866   effPIDAcc->GetYaxis()->SetTitle("Efficiency");
 867   effPIDAcc->GetYaxis()->SetRangeUser(0.,1.);
 868   effPIDAcc->SetMarkerStyle(22);
 869   effPIDAcc->SetMarkerColor(kBlue);
 870   effPIDAcc->SetLineColor(kBlack);
 871   effPIDAcc->SetStats(kFALSE);
 872   cEff->cd(4);
 873   effPIDAcc->Draw("e");
 874   delete effCalc;
 875
 876   // cleanup
 877   //delete allCandidates; delete allElectrons; delete allSignals; delete allFakes;
 878   //delete electronPurity; delete signalPurity; delete fakeContamination;
 879 }
 880
 881 //____________________________________________________________
 882 void AliAnalysisTaskHFE::MakeParticleContainer(){
 883   //
 884   // Create the particle container for the correction framework manager and
 885   // link it
 886   //
 887   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 888   const Double_t kPtmin = 0.1, kPtmax = 10.;
 889   const Double_t kEtamin = -0.9, kEtamax = 0.9;
 890   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 891
 892   //arrays for the number of bins in each dimension
 893   Int_t iBin[kNvar];
 894   iBin[0] = 40; //bins in pt
 895   iBin[1] =  8; //bins in eta
 896   iBin[2] = 18; // bins in phi
 897
 898   //arrays for lower bounds :
 899   Double_t* binEdges[kNvar];
 900   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 901     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 902
 903   //values for bin lower bounds
 904   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);
 905   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 906   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 907
 908   //one "container" for MC
 909   AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutSteps + 1 + 2*(AliHFEcuts::kNcutESDSteps + 1)), kNvar, iBin);
 910
 911   //setting the bin limits
 912   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 913     container -> SetBinLimits(ivar, binEdges[ivar]);
 914   fCFM->SetParticleContainer(container);
 915
 916   //create correlation matrix for unfolding
 917   Int_t thnDim[2*kNvar];
 918   for (int k=0; k<kNvar; k++) {
 919     //first half  : reconstructed
 920     //second half : MC
 921     thnDim[k]      = iBin[k];
 922     thnDim[k+kNvar] = iBin[k];
 923   }
 924
 925   if(!fCorrelation) fCorrelation = new TList();
 926   fCorrelation->SetName("correlation");
 927
 928   THnSparseF *correlation0 = new THnSparseF("correlationstepbeforePID","THnSparse with correlations",2*kNvar,thnDim);
 929   THnSparseF *correlation1 = new THnSparseF("correlationstepafterPID","THnSparse with correlations",2*kNvar,thnDim);
 930   for (int k=0; k<kNvar; k++) {
 931     correlation0->SetBinEdges(k,binEdges[k]);
 932     correlation0->SetBinEdges(k+kNvar,binEdges[k]);
 933     correlation1->SetBinEdges(k,binEdges[k]);
 934     correlation1->SetBinEdges(k+kNvar,binEdges[k]);
 935   }
 936   correlation0->Sumw2();
 937   correlation1->Sumw2();
 938
 939   fCorrelation->AddAt(correlation0,0);
 940   fCorrelation->AddAt(correlation1,1);
 941
 942   // Add a histogram for Fake electrons
 943   const Int_t nDim=5;
 944   Int_t nBin[nDim] = {40, 8, 18, 3, 3};
 945   Double_t* binEdges2[nDim];
 946   for(Int_t ivar = 0; ivar < nDim; ivar++)
 947     binEdges2[ivar] = new Double_t[nBin[ivar] + 1];
 948
 949   //values for bin lower bounds
 950   for(Int_t i=0; i<=nBin[0]; i++) binEdges2[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/nBin[0]*(Double_t)i);
 951   for(Int_t i=0; i<=nBin[1]; i++) binEdges2[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/nBin[1]*(Double_t)i;
 952   for(Int_t i=0; i<=nBin[2]; i++) binEdges2[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/nBin[2]*(Double_t)i;
 953   for(Int_t i=0; i<=nBin[3]; i++) binEdges2[3][i] = i;
 954   for(Int_t i=0; i<=nBin[4]; i++) binEdges2[4][i] = i;
 955
 956   fPIDperformance = new THnSparseF("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad] type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin);
 957   for(Int_t idim = 0; idim < nDim; idim++)
 958     fPIDperformance->SetBinEdges(idim, binEdges2[idim]);
 959 }
 960
 961 //____________________________________________________________
 962 void AliAnalysisTaskHFE::AddPIDdetector(Char_t *detector){
 963   //
 964   // Adding PID detector to the task
 965   //
 966   if(!fPIDdetectors.Length())
 967     fPIDdetectors = detector;
 968   else
 969     fPIDdetectors += Form(":%s", detector);
 970 }
 971
 972 //____________________________________________________________
 973 void AliAnalysisTaskHFE::PrintStatus() const {
 974   //
 975   // Print Analysis status
 976   //
 977   printf("\n\tAnalysis Settings\n\t========================================\n\n");
 978   printf("\tSecondary Vertex finding: %s\n", IsSecVtxOn() ? "YES" : "NO");
 979   printf("\tPrimary Vertex resolution: %s\n", IsPriVtxOn() ? "YES" : "NO");
 980   printf("\n");
 981   printf("\tParticle Identification Detectors:\n");
 982   TObjArray *detectors = fPIDdetectors.Tokenize(":");
 983   for(Int_t idet = 0; idet < detectors->GetEntries(); idet++)
 984     printf("\t\t%s\n", (dynamic_cast<TObjString *>(detectors->At(idet)))->String().Data());
 985   printf("\n");
 986   printf("\tQA: \n");
 987   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
 988   printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsInDebugMode()) ? "YES" : "NO");
 989   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
 990   printf("\n");
 991 }
 992
 993 //____________________________________________________________
 994 AliAnalysisTaskHFE::LabelContainer::LabelContainer(Int_t capacity):
 995   fContainer(NULL),
 996   fBegin(NULL),
 997   fEnd(NULL),
 998   fLast(NULL),
 999   fCurrent(NULL)
1000 {
1001   //
1002   // Default constructor
1003   //
1004   fContainer = new Int_t[capacity];
1005   fBegin = &fContainer[0];
1006   fEnd = &fContainer[capacity - 1];
1007   fLast = fCurrent = fBegin;
1008 }
1009
1010 //____________________________________________________________
1011 Bool_t AliAnalysisTaskHFE::LabelContainer::Append(Int_t label){
1012   //
1013   // Add Label to the container
1014   //
1015   if(fLast > fEnd) return kFALSE;
1016   *fLast++ = label;
1017   return kTRUE;
1018 }
1019
1020 //____________________________________________________________
1021 Bool_t AliAnalysisTaskHFE::LabelContainer::Find(Int_t label) const {
1022   //
1023   // Find track in the list of labels
1024   //
1025   for(Int_t *entry = fBegin; entry <= fLast; entry++)
1026     if(*entry == label) return kTRUE;
1027   return kFALSE;
1028 }
1029
1030 //____________________________________________________________
1031 Int_t AliAnalysisTaskHFE::LabelContainer::Next(){
1032   //
1033   // Mimic iterator
1034   //
1035   if(fCurrent > fLast) return -1;
1036   return *fCurrent++;
1037 }
1038
1039 //____________________________________________________________
1040 Int_t AliAnalysisTaskHFE::IsSignalElectron(AliVParticle *fTrack) const{
1041   //
1042   // Checks whether the identified electron track is coming from heavy flavour
1043   // returns 0 in case of no signal, 1 in case of charm and 2 in case of Bottom
1044   //
1045   enum{
1046     kNoSignal = 0,
1047     kCharm = 1,
1048     kBeauty = 2
1049   };
1050   AliMCParticle *mctrack = NULL;
1051   TString objname = fTrack->IsA()->GetName();
1052   if(!objname.CompareTo("AliESDtrack")){
1053     AliDebug(2, "Checking signal for ESD track");
1054     AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(fTrack);
1055     mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(esdtrack->GetLabel())));
1056   }
1057   else if(!objname.CompareTo("AliAODtrack")){
1058     AliError("AOD Analysis not implemented yet");
1059     return kNoSignal;
1060   }
1061   else if(!objname.CompareTo("AliMCParticle")){
1062     AliDebug(2, "Checking signal for MC track");
1063     mctrack = dynamic_cast<AliMCParticle *>(fTrack);
1064   }
1065   else{
1066     AliError("Input object not supported");
1067     return kNoSignal;
1068   }
1069   if(!mctrack) return kNoSignal;
1070   TParticle *ecand = mctrack->Particle();
1071   if(TMath::Abs(ecand->GetPdgCode()) != 11) return kNoSignal; // electron candidate not true electron
1072   Int_t motherLabel = TMath::Abs(ecand->GetFirstMother());
1073   AliDebug(3, Form("mother label: %d\n", motherLabel));
1074   if(!motherLabel) return kNoSignal; // mother track unknown
1075   AliMCParticle *motherTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(motherLabel));
1076   if(!motherTrack) return kNoSignal;
1077   TParticle *mparticle = motherTrack->Particle();
1078   Int_t pid = TMath::Abs(mparticle->GetPdgCode());
1079   AliDebug(3, Form("PDG code: %d\n", pid));
1080
1081   // identify signal according to Pdg Code
1082   if((pid % 1000) / 100 == 4) return kCharm;    // charmed meson, 3rd position in pdg code == 4
1083   if(pid / 1000 == 4) return kCharm;            // charmed baryon, 4th position in pdg code == 4
1084   if((pid % 1000) / 100 == 5) return kBeauty;   // beauty meson, 3rd position in pdg code == 5
1085   if(pid / 1000 == 5) return kBeauty;           // beauty baryon, 4th position in pdg code == 5
1086   return kNoSignal;
1087 }
1088
1089 //__________________________________________
1090 Float_t AliAnalysisTaskHFE::GetRapidity(TParticle *part) const
1091 {
1092   //
1093   // return rapidity
1094   //
1095   Float_t rapidity;
1096   if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
1097   else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
1098   return rapidity;
1099 }
1100