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  *  MinJung Kweon <minjung@physi.uni-heidelberg.de>
  25  */
  26 #include <TAxis.h>
  27 #include <TCanvas.h>
  28 #include <TChain.h>
  29 #include <TDirectory.h>
  30 #include <TH1F.h>
  31 #include <TH1I.h>
  32 #include <TH2F.h>
  33 #include <TIterator.h>
  34 #include <TList.h>
  35 #include <TLegend.h>
  36 #include <TMath.h>
  37 #include <TObjArray.h>
  38 #include <TParticle.h>
  39 #include <TProfile.h>
  40 #include <TString.h>
  41 #include <TTree.h>
  42
  43 #include "AliCFContainer.h"
  44 #include "AliCFManager.h"
  45 #include "AliESDEvent.h"
  46 #include "AliESDInputHandler.h"
  47 #include "AliESDtrack.h"
  48 #include "AliESDtrackCuts.h"
  49 #include "AliLog.h"
  50 #include "AliAnalysisManager.h"
  51 #include "AliMCEvent.h"
  52 #include "AliMCEventHandler.h"
  53 #include "AliMCParticle.h"
  54 #include "AliPID.h"
  55 #include "AliStack.h"
  56
  57 #include "AliHFEpid.h"
  58 #include "AliHFEcuts.h"
  59 #include "AliHFEmcQA.h"
  60 #include "AliHFEsecVtx.h"
  61 #include "AliAnalysisTaskHFE.h"
  62
  63 //____________________________________________________________
  64 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  65   AliAnalysisTask("PID efficiency Analysis", "")
  66   , fQAlevel(0)
  67   , fPIDdetectors("")
  68   , fESD(0x0)
  69   , fMC(0x0)
  70   , fCFM(0x0)
  71   , fCorrelation(0x0)
  72   , fFakeElectrons(0x0)
  73   , fPID(0x0)
  74   , fCuts(0x0)
  75   , fSecVtx(0x0)
  76   , fMCQA(0x0)
  77   , fNEvents(0x0)
  78   , fNElectronTracksEvent(0x0)
  79   , fQA(0x0)
  80   , fOutput(0x0)
  81   , fHistMCQA(0x0)
  82   , fHistSECVTX(0x0)
  83 {
  84   //
  85   // Default constructor
  86   //
  87   DefineInput(0, TChain::Class());
  88   DefineOutput(0, TH1I::Class());
  89   DefineOutput(1, TList::Class());
  90   DefineOutput(2, TList::Class());
  91
  92   // Initialize cuts
  93   fCuts = new AliHFEcuts;
  94   fPID = new AliHFEpid;
  95 }
  96
  97 //____________________________________________________________
  98 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
  99   AliAnalysisTask(ref)
 100   , fQAlevel(ref.fQAlevel)
 101   , fPIDdetectors(ref.fPIDdetectors)
 102   , fESD(ref.fESD)
 103   , fMC(ref.fMC)
 104   , fCFM(ref.fCFM)
 105   , fCorrelation(ref.fCorrelation)
 106   , fFakeElectrons(ref.fFakeElectrons)
 107   , fPID(ref.fPID)
 108   , fCuts(ref.fCuts)
 109   , fSecVtx(ref.fSecVtx)
 110   , fMCQA(ref.fMCQA)
 111   , fNEvents(ref.fNEvents)
 112   , fNElectronTracksEvent(ref.fNElectronTracksEvent)
 113   , fQA(ref.fQA)
 114   , fOutput(ref.fOutput)
 115   , fHistMCQA(ref.fHistMCQA)
 116   , fHistSECVTX(ref.fHistSECVTX)
 117 {
 118   //
 119   // Copy Constructor
 120   //
 121 }
 122
 123 //____________________________________________________________
 124 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
 125   //
 126   // Assignment operator
 127   //
 128   if(this == &ref) return *this;
 129   AliAnalysisTask::operator=(ref);
 130   fQAlevel = ref.fQAlevel;
 131   fPIDdetectors = ref.fPIDdetectors;
 132   fESD = ref.fESD;
 133   fMC = ref.fMC;
 134   fCFM = ref.fCFM;
 135   fCorrelation = ref.fCorrelation;
 136   fFakeElectrons = ref.fFakeElectrons;
 137   fPID = ref.fPID;
 138   fCuts = ref.fCuts;
 139   fSecVtx = ref.fSecVtx;
 140   fMCQA = ref.fMCQA;
 141   fNEvents = ref.fNEvents;
 142   fNElectronTracksEvent = ref.fNElectronTracksEvent;
 143   fQA = ref.fQA;
 144   fOutput = ref.fOutput;
 145   fHistMCQA = ref.fHistMCQA;
 146   fHistSECVTX = ref.fHistSECVTX;
 147   return *this;
 148 }
 149
 150 //____________________________________________________________
 151 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
 152   //
 153   // Destructor
 154   //
 155   if(fESD) delete fESD;
 156   if(fMC) delete fMC;
 157   if(fPID) delete fPID;
 158   if(fQA){
 159     fQA->Clear();
 160     delete fQA;
 161   }
 162   if(fOutput){
 163     fOutput->Clear();
 164     delete fOutput;
 165   }
 166   if(fHistMCQA){
 167     fHistMCQA->Clear();
 168     delete fHistMCQA;
 169   }
 170   if(fHistSECVTX){
 171     fHistSECVTX->Clear();
 172     delete fHistSECVTX;
 173   }
 174   if(fCuts) delete fCuts;
 175   if(fSecVtx) delete fSecVtx;
 176   if(fMCQA) delete fMCQA;
 177   if(fNEvents) delete fNEvents;
 178   if(fCorrelation) delete fCorrelation;
 179   if(fFakeElectrons) delete fFakeElectrons;
 180 }
 181
 182 //____________________________________________________________
 183 void AliAnalysisTaskHFE::ConnectInputData(Option_t *){
 184 /*  TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
 185   if(!esdchain){
 186     AliError("ESD chain empty");
 187     return;
 188   } else {
 189     esdchain->SetBranchStatus("Tracks", 1);
 190   }
 191 */
 192   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 193   if(!esdH){
 194     AliError("No ESD input handler");
 195     return;
 196   } else {
 197     fESD = esdH->GetEvent();
 198   }
 199   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 200   if(!mcH){
 201     AliError("No MC truth handler");
 202     return;
 203   } else {
 204     fMC = mcH->MCEvent();
 205   }
 206 }
 207
 208 //____________________________________________________________
 209 void AliAnalysisTaskHFE::CreateOutputObjects(){
 210   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
 211   fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
 212   // First Step: TRD alone
 213   if(!fQA) fQA = new TList;
 214   fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
 215   fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
 216   fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
 217   fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
 218   fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
 219   fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
 220   fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
 221
 222   if(!fOutput) fOutput = new TList;
 223   // Initialize correction Framework and Cuts
 224   fCFM = new AliCFManager;
 225   MakeParticleContainer();
 226   // Temporary fix: Initialize particle cuts with 0x0
 227   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
 228     fCFM->SetParticleCutsList(istep, 0x0);
 229   if(IsQAOn(kCUTqa)){
 230     AliInfo("QA on for Cuts");
 231     fCuts->SetDebugMode();
 232     fQA->Add(fCuts->GetQAhistograms());
 233   }
 234   fCuts->CreateStandardCuts();
 235   fCuts->SetMinNTrackletsTRD(0);  // Minimal requirement to get a minimum biased electron sample (only TPC and ITS requirements allowed)
 236   fCuts->Initialize(fCFM);
 237   // add output objects to the List
 238   fOutput->AddAt(fCFM->GetParticleContainer(), 0);
 239   fOutput->AddAt(fCorrelation, 1);
 240   fOutput->AddAt(fFakeElectrons, 2);
 241   fOutput->AddAt(fNElectronTracksEvent, 3);
 242
 243   // Initialize PID
 244   if(IsQAOn(kPIDqa)){
 245     AliInfo("PID QA switched on");
 246     fPID->SetQAOn();
 247     fQA->Add(fPID->GetQAhistograms());
 248   }
 249   fPID->SetHasMCData(kTRUE);
 250   if(!fPIDdetectors.Length()) AddPIDdetector("TPC");
 251   fPID->InitializePID(fPIDdetectors.Data());     // Only restrictions to TPC allowed
 252
 253   // mcQA----------------------------------
 254   if (IsQAOn(kMCqa)) {
 255     AliInfo("MC QA on");
 256     if(!fMCQA) fMCQA = new AliHFEmcQA;
 257     if(!fHistMCQA) fHistMCQA = new TList();
 258     fHistMCQA->SetName("MCqa");
 259     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
 260     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
 261     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm
 262     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
 263     TIter next(gDirectory->GetList());
 264     TObject *obj;
 265     int counter = 0;
 266     TString objname;
 267     while ((obj = next.Next())) {
 268       objname = obj->GetName();
 269       TObjArray *toks = objname.Tokenize("_");
 270       if (toks->GetEntriesFast()){
 271         TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
 272         if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj, counter++);
 273       }
 274     }
 275     fQA->Add(fHistMCQA);
 276   }
 277
 278   // secvtx----------------------------------
 279   if (IsSecVtxOn()) {
 280     AliInfo("Secondary Vertex Analysis on");
 281     fSecVtx = new AliHFEsecVtx;
 282
 283     if(!fHistSECVTX) fHistSECVTX = new TList();
 284     fHistSECVTX->SetName("SecVtx");
 285     fSecVtx->CreateHistograms("secvtx_");
 286     TIter next(gDirectory->GetList());
 287     TObject *obj;
 288     int counter = 0;
 289     TString objname;
 290     while ((obj = next.Next())) {
 291       objname = obj->GetName();
 292       TObjArray *toks = objname.Tokenize("_");
 293       if (toks->GetEntriesFast()){
 294         TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
 295         if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj, counter++);
 296       }
 297     }
 298     fOutput->Add(fHistSECVTX);
 299   }
 300 }
 301
 302 //____________________________________________________________
 303 void AliAnalysisTaskHFE::Exec(Option_t *){
 304   //
 305   // Run the analysis
 306   //
 307   if(!fESD){
 308     AliError("No ESD Event");
 309     return;
 310   }
 311   if(!fMC){
 312     AliError("No MC Event");
 313     return;
 314   }
 315   fCFM->SetEventInfo(fMC);
 316   fPID->SetMCEvent(fMC);
 317
 318   //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
 319   Double_t container[6];
 320
 321   // Loop over the Monte Carlo tracks to see whether we have overlooked any track
 322   AliMCParticle *mctrack = 0x0;
 323   Int_t nElectrons = 0;
 324
 325   if (IsSecVtxOn()) {
 326     fSecVtx->SetEvent(fESD);
 327     fSecVtx->SetStack(fMC->Stack());
 328   }
 329
 330   // run mc QA ------------------------------------------------
 331   if (IsQAOn(kMCqa)) {
 332     AliDebug(2, "Running MC QA");
 333
 334     fMCQA->SetStack(fMC->Stack());
 335     fMCQA->Init();
 336
 337     Int_t nPrims = fMC->Stack()->GetNprimary();
 338     Int_t nMCTracks = fMC->Stack()->GetNtrack();
 339
 340     // loop over primary particles for quark and heavy hadrons
 341     for (Int_t igen = 0; igen < nPrims; igen++){
 342       fMCQA->GetQuarkKine(igen, AliHFEmcQA::kCharm);
 343       fMCQA->GetQuarkKine(igen, AliHFEmcQA::kBeauty);
 344       fMCQA->GetHadronKine(igen, AliHFEmcQA::kCharm);
 345       fMCQA->GetHadronKine(igen, AliHFEmcQA::kBeauty);
 346     }
 347     fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
 348     fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 349
 350     // loop over all tracks for decayed electrons
 351     for (Int_t igen = 0; igen < nMCTracks; igen++){
 352       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0);
 353       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0);
 354       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
 355       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
 356     }
 357
 358   } // end of MC QA loop
 359   // -----------------------------------------------------------------
 360
 361   //
 362   // Loop MC
 363   //
 364   for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
 365     mctrack = fMC->GetTrack(imc);
 366
 367     container[0] = mctrack->Pt();
 368     container[1] = mctrack->Eta();
 369     container[2] = mctrack->Phi();
 370
 371     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
 372     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
 373     (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
 374     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
 375     // find the label in the vector
 376     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
 377     nElectrons++;
 378   }
 379   (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
 380
 381   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 382
 383
 384   //
 385   // Loop ESD
 386   //
 387
 388   Int_t nbtrackrec = fESD->GetNumberOfTracks();
 389   Int_t nElectronCandidates = 0;
 390   AliESDtrack *track = 0x0, *htrack = 0x0;
 391   Int_t pid = 0;
 392   // For double counted tracks
 393   Int_t *indexRecKineTPC = new Int_t[nbtrackrec];
 394   Int_t *indexRecKineITS = new Int_t[nbtrackrec];
 395   Int_t *indexRecPrim = new Int_t[nbtrackrec];
 396   Int_t *indexHFEcutsITS = new Int_t[nbtrackrec];
 397   Int_t *indexHFEcutsTPC = new Int_t[nbtrackrec];
 398   Int_t *indexHFEcutsTRD = new Int_t[nbtrackrec];
 399   Int_t *indexPid = new Int_t[nbtrackrec];
 400
 401   for(Int_t nb = 0; nb < nbtrackrec; nb++){
 402     indexRecKineTPC[nb] = -1;
 403     indexRecKineITS[nb] = -1;
 404     indexRecPrim[nb]    = -1;
 405     indexHFEcutsITS[nb] = -1;
 406     indexHFEcutsTPC[nb] = -1;
 407     indexHFEcutsTRD[nb] = -1;
 408     indexPid[nb]        = -1;
 409   }
 410
 411   Bool_t alreadyseen = kFALSE;
 412
 413   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 414
 415     track = fESD->GetTrack(itrack);
 416
 417     container[0] = track->Pt();
 418     container[1] = track->Eta();
 419     container[2] = track->Phi();
 420
 421     Bool_t signal = kTRUE;
 422
 423     // RecKine: TPC cuts
 424     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineTPC, track)) continue;
 425     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineTPC);
 426
 427
 428     // Check if it is signal electrons
 429     if(!(mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())))) continue;
 430
 431     container[3] = mctrack->Pt();
 432     container[4] = mctrack->Eta();
 433     container[5] = mctrack->Phi();
 434
 435     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
 436
 437     if(signal) {
 438       alreadyseen = kFALSE;
 439       for(Int_t nb = 0; nb < itrack; nb++){
 440         if(indexRecKineTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
 441       }
 442       indexRecKineTPC[itrack] = TMath::Abs(track->GetLabel());
 443
 444       fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineTPC + AliHFEcuts::kNcutESDSteps + 1));
 445       fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 2*(AliHFEcuts::kNcutESDSteps + 1)));
 446       if(alreadyseen) {
 447         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
 448       }
 449       else {
 450         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
 451       }
 452     }
 453
 454
 455     // RecKine: ITS cuts
 456     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITS, track)) continue;
 457     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineITS);
 458     if(signal) {
 459       alreadyseen = kFALSE;
 460       for(Int_t nb = 0; nb < itrack; nb++){
 461         if(indexRecKineITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
 462       }
 463       indexRecKineITS[itrack] = TMath::Abs(track->GetLabel());
 464       fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineITS + AliHFEcuts::kNcutESDSteps + 1));
 465       fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 2*(AliHFEcuts::kNcutESDSteps + 1)));
 466       if(alreadyseen) {
 467         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
 468       }
 469       else {
 470         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
 471       }
 472     }
 473
 474
 475     // Check TRD criterions (outside the correction framework)
 476     if(track->GetTRDncls()){
 477       (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
 478       (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha());    // Check the acceptance without tight cuts
 479       (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
 480       (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
 481     }
 482
 483
 484     // RecPrim
 485     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
 486     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
 487     if(signal) {
 488       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1);
 489       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1));
 490       alreadyseen = kFALSE;
 491       for(Int_t nb = 0; nb < itrack; nb++){
 492         if(indexRecPrim[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
 493       }
 494       indexRecPrim[itrack] = TMath::Abs(track->GetLabel());
 495       if(alreadyseen) {
 496         fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 4*(AliHFEcuts::kNcutESDSteps + 1));
 497       }
 498       else {
 499         fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 3*(AliHFEcuts::kNcutESDSteps + 1));
 500       }
 501     }
 502
 503
 504     // HFEcuts: ITS layers cuts
 505     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 506     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS);
 507     if(signal) {
 508       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1);
 509       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1));
 510       alreadyseen = kFALSE;
 511       for(Int_t nb = 0; nb < itrack; nb++){
 512         if(indexHFEcutsITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
 513       }
 514       indexHFEcutsITS[itrack] = TMath::Abs(track->GetLabel());
 515       if(alreadyseen) {
 516         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
 517       }
 518       else {
 519         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
 520       }
 521     }
 522
 523     // HFEcuts: TPC ratio clusters
 524     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 525     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC);
 526     if(signal) {
 527       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC + AliHFEcuts::kNcutESDSteps + 1);
 528       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTPC + 2*(AliHFEcuts::kNcutESDSteps + 1));
 529       alreadyseen = kFALSE;
 530       for(Int_t nb = 0; nb < itrack; nb++){
 531         if(indexHFEcutsTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
 532       }
 533       indexHFEcutsTPC[itrack] = TMath::Abs(track->GetLabel());
 534       if(alreadyseen) {
 535         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
 536       }
 537       else {
 538         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
 539       }
 540     }
 541
 542     // HFEcuts: Nb of tracklets TRD
 543     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
 544     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD);
 545     if(signal) {
 546       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 1);
 547       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1));
 548       alreadyseen = kFALSE;
 549       for(Int_t nb = 0; nb < itrack; nb++){
 550         if(indexHFEcutsTRD[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
 551       }
 552       indexHFEcutsTRD[itrack] = TMath::Abs(track->GetLabel());
 553       if(alreadyseen) {
 554         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 4*(AliHFEcuts::kNcutESDSteps + 1)));
 555       }
 556       else {
 557         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 3*(AliHFEcuts::kNcutESDSteps + 1)));
 558       }
 559     }
 560
 561
 562     // track accepted, do PID
 563     if(!fPID->IsSelected(track)) continue;
 564     nElectronCandidates++;
 565
 566
 567     // Fill Containers
 568     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 1);
 569     if(signal) {
 570       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 2);
 571       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1)+1);
 572       alreadyseen = kFALSE;
 573       for(Int_t nb = 0; nb < itrack; nb++){
 574         if(indexPid[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
 575       }
 576       indexPid[itrack] = TMath::Abs(track->GetLabel());
 577       if(alreadyseen) {
 578         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)));
 579       }
 580       else {
 581         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 3*(AliHFEcuts::kNcutESDSteps + 1)));
 582       }
 583       // dimensions 3&4&5 : pt,eta,phi (MC)
 584       fCorrelation->Fill(container);
 585     }
 586
 587
 588     // Track selected: distinguish between true and fake
 589     if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
 590       // pair analysis [mj]
 591       if (IsSecVtxOn()) {
 592         AliDebug(2, "Running Secondary Vertex Analysis");
 593         fSecVtx->InitAnaPair();
 594         for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
 595           htrack = fESD->GetTrack(jtrack);
 596           if ( itrack == jtrack ) continue;
 597           //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue;
 598           if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue;
 599           fSecVtx->AnaPair(track, htrack, jtrack);
 600         }
 601         // based on the partner of e info, you run secandary vertexing function
 602         fSecVtx->RunSECVTX(track);
 603       } // end of pair analysis
 604     } else {
 605       // Fill THnSparse with the information for Fake Electrons
 606       switch(pid){
 607       case 13:    container[3] = AliPID::kMuon;
 608       case 211:   container[3] = AliPID::kPion;
 609       case 321:   container[3] = AliPID::kKaon;
 610       case 2212:  container[3] = AliPID::kProton;
 611       }
 612       fFakeElectrons->Fill(container);
 613     }
 614   }
 615   fNEvents->Fill(1);
 616   fNElectronTracksEvent->Fill(nElectronCandidates);
 617
 618   // release memory
 619   delete[] indexRecKineTPC;
 620   delete[] indexRecKineITS;
 621   delete[] indexRecPrim;
 622   delete[] indexHFEcutsITS;
 623   delete[] indexHFEcutsTPC;
 624   delete[] indexHFEcutsTRD;
 625   delete[] indexPid;
 626
 627   // Done!!!
 628   PostData(0, fNEvents);
 629   PostData(1, fOutput);
 630   PostData(2, fQA);
 631 }
 632
 633 //____________________________________________________________
 634 void AliAnalysisTaskHFE::Terminate(Option_t *){
 635   //
 636   // Terminate not implemented at the moment
 637   //
 638 }
 639
 640 //____________________________________________________________
 641 void AliAnalysisTaskHFE::MakeParticleContainer(){
 642   //
 643   // Create the particle container for the correction framework manager and
 644   // link it
 645   //
 646   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 647   const Double_t kPtmin = 0.1, kPtmax = 10.;
 648   const Double_t kEtamin = -0.9, kEtamax = 0.9;
 649   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 650
 651   //arrays for the number of bins in each dimension
 652   Int_t iBin[kNvar];
 653   iBin[0] = 40; //bins in pt
 654   iBin[1] =  8; //bins in eta
 655   iBin[2] = 18; // bins in phi
 656
 657   //arrays for lower bounds :
 658   Double_t* binEdges[kNvar];
 659   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 660     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 661
 662   //values for bin lower bounds
 663   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);
 664   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 665   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 666
 667   //one "container" for MC
 668   AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutSteps + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)), kNvar, iBin);
 669
 670   //setting the bin limits
 671   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 672     container -> SetBinLimits(ivar, binEdges[ivar]);
 673   fCFM->SetParticleContainer(container);
 674
 675   //create correlation matrix for unfolding
 676   Int_t thnDim[2*kNvar];
 677   for (int k=0; k<kNvar; k++) {
 678     //first half  : reconstructed
 679     //second half : MC
 680     thnDim[k]      = iBin[k];
 681     thnDim[k+kNvar] = iBin[k];
 682   }
 683
 684   fCorrelation = new THnSparseF("correlation","THnSparse with correlations",2*kNvar,thnDim);
 685   for (int k=0; k<kNvar; k++) {
 686     fCorrelation->SetBinEdges(k,binEdges[k]);
 687     fCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
 688   }
 689   fCorrelation->Sumw2();
 690
 691   // Add a histogram for Fake electrons
 692   thnDim[3] = AliPID::kSPECIES;
 693   fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
 694   for(Int_t idim = 0; idim < kNvar; idim++)
 695     fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
 696 }
 697
 698 void AliAnalysisTaskHFE::AddPIDdetector(Char_t *detector){
 699   if(!fPIDdetectors.Length())
 700     fPIDdetectors = detector;
 701   else
 702     fPIDdetectors += Form(":%s", detector);
 703 }
 704
 705 //____________________________________________________________
 706 void AliAnalysisTaskHFE::PrintStatus(){
 707   printf("\n\tAnalysis Settings\n\t========================================\n\n");
 708   printf("\tSecondary Vertex finding: %s\n", IsSecVtxOn() ? "YES" : "NO");
 709   printf("\tPrimary Vertex resolution: %s\n", IsPriVtxOn() ? "YES" : "NO");
 710   printf("\n");
 711   printf("\tParticle Identification Detectors:\n");
 712   TObjArray *detectors = fPIDdetectors.Tokenize(":");
 713   for(Int_t idet = 0; idet < detectors->GetEntries(); idet++)
 714     printf("\t\t%s\n", (dynamic_cast<TObjString *>(detectors->At(idet)))->String().Data());
 715   printf("\n");
 716   printf("\tQA: \n");
 717   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
 718   printf("\t\tCUTS: %s\n", IsQAOn(kCUTqa) ? "YES" : "NO");
 719   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
 720   printf("\n");
 721 }