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