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 <TH3D.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 <TF1.h>
  42 #include <TTree.h>
  43
  44 #include "AliAODInputHandler.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODVertex.h"
  48 #include "AliCFContainer.h"
  49 #include "AliCFManager.h"
  50 #include "AliESDEvent.h"
  51 #include "AliESDInputHandler.h"
  52 #include "AliESDpid.h"
  53 #include "AliESDtrack.h"
  54 #include "AliCentrality.h"
  55 #include "AliLog.h"
  56 #include "AliAnalysisManager.h"
  57 #include "AliMCEvent.h"
  58 #include "AliMCEventHandler.h"
  59 #include "AliMCParticle.h"
  60 #include "AliPID.h"
  61 #include "AliStack.h"
  62 #include "AliTriggerAnalysis.h"
  63 #include "AliVVertex.h"
  64
  65 #include "AliHFEcollection.h"
  66 #include "AliHFEcontainer.h"
  67 #include "AliHFEcuts.h"
  68 #include "AliHFEelecbackground.h"
  69 #include "AliHFEmcQA.h"
  70 #include "AliHFEpairs.h"
  71 #include "AliHFEpid.h"
  72 #include "AliHFEpidQAmanager.h"
  73 #include "AliHFEpostAnalysis.h"
  74 #include "AliHFEsecVtxs.h"
  75 #include "AliHFEsecVtx.h"
  76 #include "AliHFEsignalCuts.h"
  77 #include "AliHFEtaggedTrackAnalysis.h"
  78 #include "AliHFEtools.h"
  79 #include "AliHFEvarManager.h"
  80 #include "AliAnalysisTaskHFE.h"
  81
  82 ClassImp(AliAnalysisTaskHFE)
  83
  84 //____________________________________________________________
  85 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  86   AliAnalysisTaskSE("PID efficiency Analysis")
  87   , fQAlevel(0)
  88   , fPlugins(0)
  89   , fFillSignalOnly(kTRUE)
  90   , fBackGroundFactorApply(kFALSE)
  91   , fRemovePileUp(kFALSE)
  92   , fIdentifiedAsPileUp(kFALSE)
  93   , fIdentifiedAsOutInz(kFALSE)
  94   , fPassTheEventCut(kFALSE)
  95   , fHasSpecialTriggerSelection(kFALSE)
  96   , fSpecialTrigger("NONE")
  97   , fCentralityF(99.0)
  98   , fContributors(0.5)
  99   , fWeightBackGround(0.)
 100   , fVz(0.0)
 101   , fContainer(NULL)
 102   , fVarManager(NULL)
 103   , fSignalCuts(NULL)
 104   , fCFM(NULL)
 105   , fTriggerAnalysis(NULL)
 106   , fPID(NULL)
 107   , fPIDqa(NULL)
 108   , fPIDpreselect(NULL)
 109   , fCuts(NULL)
 110   , fTaggedTrackCuts(NULL)
 111   , fCleanTaggedTrack(kFALSE)
 112   , fVariablesTRDTaggedTrack(kFALSE)
 113   , fCutspreselect(NULL)
 114   , fSecVtx(NULL)
 115   , fElecBackGround(NULL)
 116   , fMCQA(NULL)
 117   , fTaggedTrackAnalysis(NULL)
 118   , fQA(NULL)
 119   , fOutput(NULL)
 120   , fHistMCQA(NULL)
 121   , fHistSECVTX(NULL)
 122   , fHistELECBACKGROUND(NULL)
 123   , fQACollection(NULL)
 124 {
 125   //
 126   // Dummy constructor
 127   //
 128 }
 129
 130 //____________________________________________________________
 131 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
 132   AliAnalysisTaskSE(name)
 133   , fQAlevel(0)
 134   , fPlugins(0)
 135   , fFillSignalOnly(kTRUE)
 136   , fBackGroundFactorApply(kFALSE)
 137   , fRemovePileUp(kFALSE)
 138   , fIdentifiedAsPileUp(kFALSE)
 139   , fIdentifiedAsOutInz(kFALSE)
 140   , fPassTheEventCut(kFALSE)
 141   , fHasSpecialTriggerSelection(kFALSE)
 142   , fSpecialTrigger("NONE")
 143   , fCentralityF(99.0)
 144   , fContributors(0.5)
 145   , fWeightBackGround(0.)
 146   , fVz(0.0)
 147   , fContainer(NULL)
 148   , fVarManager(NULL)
 149   , fSignalCuts(NULL)
 150   , fCFM(NULL)
 151   , fTriggerAnalysis(NULL)
 152   , fPID(NULL)
 153   , fPIDqa(NULL)
 154   , fPIDpreselect(NULL)
 155   , fCuts(NULL)
 156   , fTaggedTrackCuts(NULL)
 157   , fCleanTaggedTrack(kFALSE)
 158   , fVariablesTRDTaggedTrack(kFALSE)
 159   , fCutspreselect(NULL)
 160   , fSecVtx(NULL)
 161   , fElecBackGround(NULL)
 162   , fMCQA(NULL)
 163   , fTaggedTrackAnalysis(NULL)
 164   , fQA(NULL)
 165   , fOutput(NULL)
 166   , fHistMCQA(NULL)
 167   , fHistSECVTX(NULL)
 168   , fHistELECBACKGROUND(NULL)
 169   , fQACollection(0x0)
 170 {
 171   //
 172   // Default constructor
 173   //
 174   DefineOutput(1, TList::Class());
 175   DefineOutput(2, TList::Class());
 176
 177   fPID = new AliHFEpid("hfePid");
 178   fPIDqa = new AliHFEpidQAmanager;
 179   fVarManager = new AliHFEvarManager("hfeVarManager");
 180 }
 181
 182 //____________________________________________________________
 183 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
 184   AliAnalysisTaskSE(ref)
 185   , fQAlevel(0)
 186   , fPlugins(0)
 187   , fFillSignalOnly(ref.fFillSignalOnly)
 188   , fBackGroundFactorApply(ref.fBackGroundFactorApply)
 189   , fRemovePileUp(ref.fRemovePileUp)
 190   , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
 191   , fIdentifiedAsOutInz(ref.fIdentifiedAsOutInz)
 192   , fPassTheEventCut(ref.fPassTheEventCut)
 193   , fHasSpecialTriggerSelection(ref.fHasSpecialTriggerSelection)
 194   , fSpecialTrigger(ref.fSpecialTrigger)
 195   , fCentralityF(ref.fCentralityF)
 196   , fContributors(ref.fContributors)
 197   , fWeightBackGround(ref.fWeightBackGround)
 198   , fVz(ref.fVz)
 199   , fContainer(NULL)
 200   , fVarManager(NULL)
 201   , fSignalCuts(NULL)
 202   , fCFM(NULL)
 203   , fTriggerAnalysis(NULL)
 204   , fPID(NULL)
 205   , fPIDqa(NULL)
 206   , fPIDpreselect(NULL)
 207   , fCuts(NULL)
 208   , fTaggedTrackCuts(NULL)
 209   , fCleanTaggedTrack(ref.fCleanTaggedTrack)
 210   , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
 211   , fCutspreselect(NULL)
 212   , fSecVtx(NULL)
 213   , fElecBackGround(NULL)
 214   , fMCQA(NULL)
 215   , fTaggedTrackAnalysis(NULL)
 216   , fQA(NULL)
 217   , fOutput(NULL)
 218   , fHistMCQA(NULL)
 219   , fHistSECVTX(NULL)
 220   , fHistELECBACKGROUND(NULL)
 221   , fQACollection(NULL)
 222 {
 223   //
 224   // Copy Constructor
 225   //
 226   ref.Copy(*this);
 227 }
 228
 229 //____________________________________________________________
 230 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
 231   //
 232   // Assignment operator
 233   //
 234   if(this == &ref)
 235     ref.Copy(*this);
 236   return *this;
 237 }
 238
 239 //____________________________________________________________
 240 void AliAnalysisTaskHFE::Copy(TObject &o) const {
 241   //
 242   // Copy into object o
 243   //
 244   AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
 245   target.fQAlevel = fQAlevel;
 246   target.fPlugins = fPlugins;
 247   target.fFillSignalOnly = fFillSignalOnly;
 248   target.fBackGroundFactorApply = fBackGroundFactorApply;
 249   target.fRemovePileUp = fRemovePileUp;
 250   target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
 251   target.fIdentifiedAsOutInz = fIdentifiedAsOutInz;
 252   target.fPassTheEventCut = fPassTheEventCut;
 253   target.fHasSpecialTriggerSelection = fHasSpecialTriggerSelection;
 254   target.fSpecialTrigger = fSpecialTrigger;
 255   target.fCentralityF = fCentralityF;
 256   target.fContributors = fContributors;
 257   target.fWeightBackGround = fWeightBackGround;
 258   target.fVz = fVz;
 259   target.fContainer = fContainer;
 260   target.fVarManager = fVarManager;
 261   target.fSignalCuts = fSignalCuts;
 262   target.fCFM = fCFM;
 263   target.fTriggerAnalysis = fTriggerAnalysis;
 264   target.fPID = fPID;
 265   target.fPIDqa = fPIDqa;
 266   target.fPIDpreselect = fPIDpreselect;
 267   target.fCuts = fCuts;
 268   target.fTaggedTrackCuts = fTaggedTrackCuts;
 269   target.fCleanTaggedTrack = fCleanTaggedTrack;
 270   target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
 271   target.fCutspreselect = fCutspreselect;
 272   target.fSecVtx = fSecVtx;
 273   target.fElecBackGround = fElecBackGround;
 274   target.fMCQA = fMCQA;
 275   target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
 276   target.fQA = fQA;
 277   target.fOutput = fOutput;
 278   target.fHistMCQA = fHistMCQA;
 279   target.fHistSECVTX = fHistSECVTX;
 280   target.fHistELECBACKGROUND = fHistELECBACKGROUND;
 281   target.fQACollection = fQACollection;
 282 }
 283
 284 //____________________________________________________________
 285 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
 286   //
 287   // Destructor
 288   //
 289   if(fPID) delete fPID;
 290   if(fVarManager) delete fVarManager;
 291   if(fPIDqa) delete fPIDqa;
 292   if(fSignalCuts) delete fSignalCuts;
 293   if(fCFM) delete fCFM;
 294   if(fSecVtx) delete fSecVtx;
 295   if(fMCQA) delete fMCQA;
 296   if(fElecBackGround) delete fElecBackGround;
 297   if(fTriggerAnalysis) delete fTriggerAnalysis;
 298   if(fPIDpreselect) delete fPIDpreselect;
 299   if(fQA) delete fQA;
 300   if(fOutput) delete fOutput;
 301 }
 302
 303 //____________________________________________________________
 304 void AliAnalysisTaskHFE::UserCreateOutputObjects(){
 305   //
 306   // Creating output container and output objects
 307   // Here we also Initialize the correction framework container and
 308   // the objects for
 309   // - PID
 310   // - MC QA
 311   // - SecVtx
 312   // QA histograms are created if requested
 313   // Called once per worker
 314   //
 315   AliDebug(3, "Creating Output Objects");
 316   // Automatic determination of the analysis mode
 317   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 318   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
 319     SetAODAnalysis();
 320   } else {
 321     SetESDAnalysis();
 322     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 323       SetHasMCData();
 324   }
 325   printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
 326   printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
 327
 328   // Enable Trigger Analysis
 329   fTriggerAnalysis = new AliTriggerAnalysis;
 330   fTriggerAnalysis->EnableHistograms();
 331   fTriggerAnalysis->SetAnalyzeMC(HasMCData());
 332
 333
 334   // Make lists for Output
 335   if(!fQA) fQA = new TList;
 336   fQA->SetOwner();
 337   if(!fOutput) fOutput = new TList;
 338   fOutput->SetOwner();
 339
 340   // First Part: Make QA histograms
 341   fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
 342   fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
 343   fQACollection->CreateProfile("conr", "Electron PID contamination", 20, 0, 20);
 344   fQACollection->CreateTH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi());
 345   fQACollection->CreateTH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi());
 346   fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
 347   fQACollection->CreateProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20);
 348   fQACollection->CreateProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20);
 349   fQACollection->CreateTH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20);
 350   fQACollection->CreateTH1F("mccharge", "MC Charge", 200, -100, 100);
 351   fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
 352   // Temporary histograms for TPC number of clusters for all signal tracks (MC true electrons) and for selected tracks (Markus Fasel)
 353   fQACollection->CreateTH2F("TPCclusters2_1_Signal", "TPCclusterInfo for findable clusters for 2 neighbors for signal tracks", 30, 0.1, 10., 162, 0., 161.);
 354   fQACollection->CreateTH2F("TPCclusters2_0_Signal", "TPCclusterInfo for the ratio for 2 neighbors for signal tracks", 30, 0.1, 10., 100, 0., 1.);
 355   fQACollection->CreateTH2F("TPCclusters2_1_Selected", "TPCclusterInfo for findable clusters for 2 neighbors for selected tracks", 30, 0.1, 10., 162, 0., 161.);
 356   fQACollection->CreateTH2F("TPCclusters2_0_Selected", "TPCclusterInfo for the ratio for 2 neighbors for selected tracks", 30, 0.1, 10., 110, 0., 1.1);
 357   fQACollection->CreateTH2F("TPCncls_Signal", "TPC Number of clusters for signal tracks", 30, 0.1, 10., 162, 0., 161.);
 358   fQACollection->CreateTH2F("TPCclr_Signal", "TPC cluster ratio for signal tracks", 30, 0.1, 10., 110, 0., 1.1);
 359   fQACollection->BinLogAxis("TPCclusters2_1_Signal", 0);
 360   fQACollection->BinLogAxis("TPCclusters2_0_Signal", 0);
 361   fQACollection->BinLogAxis("TPCclusters2_1_Selected", 0);
 362   fQACollection->BinLogAxis("TPCclusters2_0_Selected", 0);
 363   fQACollection->BinLogAxis("TPCncls_Signal", 0);
 364   fQACollection->BinLogAxis("TPCclr_Signal", 0);
 365
 366   InitPIDperformanceQA();
 367   InitContaminationQA();
 368   fQA->Add(fQACollection->GetList());
 369
 370   // Initialize PID
 371   fPID->SetHasMCData(HasMCData());
 372   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
 373   fPID->InitializePID();
 374   if(IsQAOn(kPIDqa)){
 375     AliInfo("PID QA switched on");
 376     fPIDqa->Initialize(fPID);
 377     fQA->Add(fPIDqa->MakeList("HFEpidQA"));
 378   }
 379
 380   // Initialize correction Framework and Cuts
 381   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack + AliHFEcuts::kNcutStepsSecvtxTrack;
 382   fCFM = new AliCFManager;
 383   fCFM->SetNStepParticle(kNcutSteps);
 384   MakeParticleContainer();
 385   MakeEventContainer();
 386   // Temporary fix: Initialize particle cuts with NULL
 387   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 388     fCFM->SetParticleCutsList(istep, NULL);
 389   if(!fCuts){
 390     AliWarning("Cuts not available. Default cuts will be used");
 391     fCuts = new AliHFEcuts;
 392     fCuts->CreateStandardCuts();
 393   }
 394   if(IsAODanalysis()) fCuts->SetAOD();
 395   // Make clone for V0 tagging step
 396   fCuts->Initialize(fCFM);
 397   if(fCuts->IsQAOn()) fQA->Add(fCuts->GetQAhistograms());
 398   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
 399   fVarManager->SetSignalCuts(fSignalCuts);
 400
 401   // add output objects to the List
 402   fOutput->AddAt(fContainer, 0);
 403   fOutput->AddAt(fCFM->GetEventContainer(), 1);
 404
 405   // mcQA----------------------------------
 406   if (HasMCData() && IsQAOn(kMCqa)) {
 407     AliInfo("MC QA on");
 408     if(!fMCQA) fMCQA = new AliHFEmcQA;
 409     if(!fHistMCQA) fHistMCQA = new TList();
 410     fHistMCQA->SetOwner();
 411     fMCQA->CreatDefaultHistograms(fHistMCQA);
 412     fQA->Add(fHistMCQA);
 413   }
 414
 415   // secvtx----------------------------------
 416   if (GetPlugin(kSecVtx)) {
 417     AliInfo("Secondary Vertex Analysis on");
 418     if(!fSecVtx) fSecVtx = new AliHFEsecVtx;
 419     fSecVtx->SetHasMCData(HasMCData());
 420
 421     if(!fHistSECVTX) fHistSECVTX = new TList();
 422     fHistSECVTX->SetOwner();
 423     fSecVtx->CreateHistograms(fHistSECVTX);
 424     fOutput->Add(fHistSECVTX);
 425   }
 426
 427   // background----------------------------------
 428   if (GetPlugin(kIsElecBackGround)) {
 429     AliInfo("Electron BackGround Analysis on");
 430     if(!fElecBackGround){
 431       AliWarning("ElecBackGround not available. Default elecbackground will be used");
 432       fElecBackGround = new AliHFEelecbackground;
 433     }
 434     fElecBackGround->SetHasMCData(HasMCData());
 435
 436     if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
 437     fHistELECBACKGROUND->SetOwner();
 438     fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
 439     fOutput->Add(fHistELECBACKGROUND);
 440   }
 441
 442   // tagged tracks
 443   if(GetPlugin(kTaggedTrackAnalysis)){
 444     AliInfo("Analysis on V0-tagged tracks enabled");
 445     fTaggedTrackAnalysis = new AliHFEtaggedTrackAnalysis(Form("taggedTrackAnalysis%s", GetName()));
 446     fTaggedTrackAnalysis->SetCuts(fTaggedTrackCuts);
 447     fTaggedTrackAnalysis->SetClean(fCleanTaggedTrack);
 448     AliHFEvarManager *varManager = fTaggedTrackAnalysis->GetVarManager();
 449     TObjArray *array = fVarManager->GetVariables();
 450     Int_t nvars = array->GetEntriesFast();
 451     TString namee;
 452     for(Int_t v = 0; v < nvars; v++) {
 453       AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
 454       if(!variable) continue;
 455       TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
 456       if(!name.CompareTo("source")) namee = TString("species");
 457       else namee = TString(name);
 458       varManager->AddVariable(namee);
 459       //printf("For AliTaggedTrackAnalysis, had the variable %s and the one used %s\n",(const char*)variable->GetName(),(const char*) namee);
 460     }
 461     if(fPIDqa->HasHighResolutionHistos())
 462       fTaggedTrackAnalysis->GetPIDqa()->SetHighResolutionHistos();
 463     fTaggedTrackAnalysis->SetPID(fPID);
 464     fTaggedTrackAnalysis->SetVariablesTRD(fVariablesTRDTaggedTrack);
 465     fTaggedTrackAnalysis->InitContainer();
 466     fOutput->Add(fTaggedTrackAnalysis->GetContainer());
 467     fQA->Add(fTaggedTrackAnalysis->GetPIDQA());
 468     fQA->Add(fTaggedTrackAnalysis->GetCutQA());
 469     fQA->Add(fTaggedTrackAnalysis->GetQAcollection());
 470   }
 471   PrintStatus();
 472 }
 473
 474 //____________________________________________________________
 475 void AliAnalysisTaskHFE::UserExec(Option_t *){
 476   //
 477   // Run the analysis
 478   //
 479   AliDebug(3, "Starting Single Event Analysis");
 480   if(!fInputEvent){
 481     AliError("Reconstructed Event not available");
 482     return;
 483   }
 484   if(HasMCData()){
 485     AliDebug(4, Form("MC Event: %p", fMCEvent));
 486     if(!fMCEvent){
 487       AliError("No MC Event, but MC Data required");
 488       return;
 489     }
 490   }
 491   if(!fCuts){
 492     AliError("HFE cuts not available");
 493     return;
 494   }
 495
 496
 497   if(IsESDanalysis() && HasMCData()){
 498     // Protect against missing MC trees
 499     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 500     if(!mcH){
 501       AliError("No MC Event Handler available");
 502       return;
 503     }
 504     if(!mcH->InitOk()) return;
 505     if(!mcH->TreeK()) return;
 506     if(!mcH->TreeTR()) return;
 507   }
 508
 509   // need the centrality for everything (MC also)
 510   fCentralityF = -100.0;
 511   if(!ReadCentrality()) fCentralityF = -100.0;
 512   //printf("pass centrality\n");
 513   //printf("Reading fCentralityF %f\n",fCentralityF);
 514
 515   // See if pile up and z in the range
 516   RejectionPileUpVertexRangeEventCut();
 517
 518   // Protect agains missing
 519   if(HasMCData()){
 520     //printf("Has MC data\n");
 521     fSignalCuts->SetMCEvent(fMCEvent);
 522     ProcessMC();  // Run the MC loop + MC QA in case MC Data are available
 523   }
 524
 525   if(IsAODanalysis()){
 526     AliAODpidUtil *aodworkingpid = AliHFEtools::GetDefaultAODPID(HasMCData());
 527     fPID->SetAODpid(aodworkingpid);
 528     ProcessAOD();
 529   } else {
 530     // Check Trigger selection
 531     if(fHasSpecialTriggerSelection){
 532       AliESDEvent *ev = dynamic_cast<AliESDEvent *>(fInputEvent);
 533       if(!(ev && ev->IsTriggerClassFired(fSpecialTrigger.Data()))) return;
 534     }
 535     AliESDInputHandler *inH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
 536     if(!inH){
 537       AliError("No ESD Input handler available");
 538       return;
 539     }
 540     AliESDpid *workingPID = inH->GetESDpid();
 541     if(!workingPID){
 542       AliDebug(1, "Using default ESD PID");
 543       workingPID = AliHFEtools::GetDefaultPID(HasMCData());
 544     } else {
 545       AliDebug(1, "Using ESD PID from the input handler");
 546     }
 547     fPID->SetESDpid(workingPID);
 548     if(fPIDpreselect) fPIDpreselect->SetESDpid(workingPID);
 549
 550     ProcessESD();
 551   }
 552   // Done!!!
 553   PostData(1, fOutput);
 554   PostData(2, fQA);
 555 }
 556
 557 //____________________________________________________________
 558 void AliAnalysisTaskHFE::Terminate(Option_t *){
 559   //
 560   // Terminate not implemented at the moment
 561   //
 562   if(GetPlugin(kPostProcess)){
 563     fOutput = dynamic_cast<TList *>(GetOutputData(1));
 564     fQA = dynamic_cast<TList *>(GetOutputData(2));
 565     if(!fOutput){
 566       AliError("Results not available");
 567       return;
 568     }
 569     if(!fQA){
 570       AliError("QA output not available");
 571       return;
 572     }
 573     fContainer = dynamic_cast<AliHFEcontainer *>(fOutput->FindObject("trackContainer"));
 574     if(!fContainer){
 575       AliError("Track container not found");
 576       return;
 577     }
 578     AliHFEpostAnalysis postanalysis;
 579     postanalysis.SetTaskResults(fContainer);
 580     TList *qalist = dynamic_cast<TList *>(fQA->FindObject("list_TaskQA"));
 581     if(!qalist){
 582       AliError("QA List not found");
 583       return;
 584     }
 585     postanalysis.SetTaskQA(qalist);
 586     printf("Running post analysis\n");
 587     //if(HasMCData())
 588     postanalysis.DrawMCSignal2Background();
 589     postanalysis.DrawEfficiency();
 590     postanalysis.DrawPIDperformance();
 591     postanalysis.DrawCutEfficiency();
 592
 593     if (GetPlugin(kIsElecBackGround)) {
 594       AliHFEelecbackground elecBackGround;
 595       TList *oe = 0x0;
 596       if(!(oe = (TList*)dynamic_cast<TList *>(fOutput->FindObject("HFEelecbackground")))){
 597         return;
 598       }
 599       elecBackGround.Load(oe);
 600       elecBackGround.Plot();
 601       elecBackGround.PostProcess();
 602     }
 603   }
 604 }
 605 //_______________________________________________________________
 606 Bool_t AliAnalysisTaskHFE::IsEventInBinZero() {
 607   //
 608   //
 609   //
 610
 611   //printf("test in IsEventInBinZero\n");
 612   if(!fInputEvent){
 613     AliError("Reconstructed Event not available");
 614     return kFALSE;
 615   }
 616
 617   // check vertex
 618   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
 619   if(!vertex) return kTRUE;
 620   //if(vertex) return kTRUE;
 621
 622   // check tracks
 623   if(fInputEvent->GetNumberOfTracks()<=0) return kTRUE;
 624   //if(fInputEvent->GetNumberOfTracks()>0) return kTRUE;
 625
 626
 627   return kFALSE;
 628
 629 }
 630 //____________________________________________________________
 631 void AliAnalysisTaskHFE::ProcessMC(){
 632   //
 633   // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
 634   // In case MC QA is on also MC QA loop is done
 635   //
 636   AliDebug(3, "Processing MC Information");
 637   Double_t eventContainer [4];
 638   eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
 639   eventContainer[2] = fCentralityF;
 640   eventContainer[3] = fContributors;
 641   fVz = eventContainer[0];
 642   //printf("z position is %f\n",eventContainer[0]);
 643   //if(fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent))
 644   fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
 645   Int_t nElectrons = 0;
 646   if(IsESDanalysis()){
 647    if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
 648     if (HasMCData() && IsQAOn(kMCqa)) {
 649       AliDebug(2, "Running MC QA");
 650
 651       if(fMCEvent->Stack()){
 652               fMCQA->SetMCEvent(fMCEvent);
 653         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
 654         fMCQA->Init();
 655
 656         fMCQA->GetMesonKine();
 657
 658         // loop over all tracks for decayed electrons
 659         for (Int_t igen = 0; igen < fMCEvent->GetNumberOfTracks(); igen++){
 660           TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
 661           if(!mcpart) continue;
 662           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
 663           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
 664           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
 665           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
 666           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0); // no accept cut
 667           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 668           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 669           if (TMath::Abs(mcpart->Eta()) < 0.9) {
 670             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 671             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 672             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 673           }
 674           if (TMath::Abs(AliHFEtools::GetRapidity(mcpart)) < 0.5) {
 675             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 676             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 677             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 678           }
 679         }
 680         fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
 681         fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 682       }
 683
 684     } // end of MC QA loop
 685    }
 686    // -----------------------------------------------------------------
 687    fCFM->SetMCEventInfo(fMCEvent);
 688    // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 689   } else {
 690     fCFM->SetMCEventInfo(fInputEvent);
 691   }
 692   // Run MC loop
 693   AliVParticle *mctrack = NULL;
 694   AliDebug(3, Form("Number of Tracks: %d", fMCEvent->GetNumberOfTracks()));
 695   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfTracks(); imc++){
 696     if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
 697     AliDebug(4, "Next MC Track");
 698     if(ProcessMCtrack(mctrack)) nElectrons++;
 699   }
 700
 701   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 702   fQACollection->Fill("nElectron", nElectrons);
 703 }
 704
 705 //____________________________________________________________
 706 void AliAnalysisTaskHFE::ProcessESD(){
 707   //
 708   // Run Analysis of reconstructed event in ESD Mode
 709   // Loop over Tracks, filter according cut steps defined in AliHFEcuts
 710   //
 711   AliDebug(3, "Processing ESD Event");
 712   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 713   if(!fESD){
 714     AliError("ESD Event required for ESD Analysis")
 715     return;
 716   }
 717
 718   // Set magnetic field if V0 task on
 719   if(fTaggedTrackAnalysis) {
 720     fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
 721     fTaggedTrackAnalysis->SetCentrality(fCentralityF);
 722   }
 723
 724   // Do event Normalization
 725   Double_t eventContainer[4];
 726   eventContainer[0] = 0.0;
 727   if(HasMCData()) eventContainer[0] = fVz;
 728   else {
 729     if(fESD->GetPrimaryVertexTracks()) eventContainer[0] = fESD->GetPrimaryVertexTracks()->GetZ();
 730   }
 731   eventContainer[1] = 0.;
 732   eventContainer[2] = fCentralityF;
 733   eventContainer[3] = fContributors;
 734   if(fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0AND))
 735     eventContainer[1] = 1.;
 736
 737   //
 738   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
 739
 740   //
 741   if(fIdentifiedAsPileUp) return;
 742   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
 743
 744   //
 745   if(TMath::Abs(fCentralityF+100.0) < 0.01) return;
 746   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
 747   //printf("In ProcessESD %f\n",fCentralityF);
 748
 749   //
 750   if(fIdentifiedAsOutInz) return;
 751   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
 752
 753   //
 754   if(!fPassTheEventCut) return;
 755   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
 756
 757
 758
 759   fContainer->NewEvent();
 760
 761   if (GetPlugin(kIsElecBackGround)) {
 762     fElecBackGround->SetEvent(fESD);
 763   }
 764   if (GetPlugin(kSecVtx)) {
 765     fSecVtx->SetEvent(fESD);
 766     fSecVtx->GetPrimaryCondition();
 767   }
 768
 769   if(HasMCData()){
 770     if (GetPlugin(kSecVtx)) {
 771       fSecVtx->SetMCEvent(fMCEvent);
 772       fSecVtx->SetMCQA(fMCQA);
 773     }
 774     if (GetPlugin(kIsElecBackGround)) {
 775       fElecBackGround->SetMCEvent(fMCEvent);
 776     }
 777   }
 778
 779   Double_t container[10];
 780   memset(container, 0, sizeof(Double_t) * 10);
 781   // container for the output THnSparse
 782   Double_t dataE[6]; // [pT, eta, Phi, type, 'C' or 'B']
 783   Int_t nElectronCandidates = 0;
 784   AliESDtrack *track = NULL, *htrack = NULL;
 785   AliMCParticle *mctrack = NULL;
 786   TParticle* mctrack4QA = NULL;
 787   Int_t pid = 0;
 788
 789   Bool_t signal = kTRUE;
 790
 791   fCFM->SetRecEventInfo(fESD);
 792   // Electron background analysis
 793   if (GetPlugin(kIsElecBackGround)) {
 794
 795     AliDebug(2, "Running BackGround Analysis");
 796
 797     fElecBackGround->Reset();
 798
 799   } // end of electron background analysis
 800   //
 801   // Loop ESD
 802   //
 803   AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
 804   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 805     AliDebug(4, "New ESD track");
 806     track = fESD->GetTrack(itrack);
 807     track->SetESDEvent(fESD);
 808
 809     // fill counts of v0-identified particles
 810     Int_t v0pid = -1;
 811     if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
 812     else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
 813     else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
 814     // here the tagged track analysis will run
 815     if(fTaggedTrackAnalysis && v0pid > -1){
 816       AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
 817       fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
 818     }
 819
 820     AliDebug(3, Form("Doing track %d, %p", itrack, track));
 821
 822     //////////////////////////////////////
 823     // preselect
 824     /////////////////////////////////////
 825     if(fPIDpreselect && fCutspreselect) {
 826       if(!PreSelectTrack(track)) continue;
 827     }
 828
 829     signal = kTRUE;
 830
 831     // Fill step without any cut
 832
 833     if(HasMCData()){
 834       // Check if it is electrons near the vertex
 835       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
 836       mctrack4QA = mctrack->Particle();
 837
 838       if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
 839       else AliDebug(3, "Signal Electron");
 840     }
 841     // Cache new Track information inside the var manager
 842     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
 843
 844     if(signal) {
 845       fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
 846       fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
 847       if((track->GetStatus() & AliESDtrack::kTPCout)
 848           && (TMath::Abs(track->Eta()) < 0.8)
 849           && (track->GetKinkIndex(0) == 0)){
 850         fQACollection->Fill("TPCclusters2_1_Signal", track->Pt(), track->GetTPCClusterInfo(2,1));
 851         fQACollection->Fill("TPCclusters2_0_Signal", track->Pt(), track->GetTPCNclsF() > 0 ?  track->GetTPCClusterInfo(2,1)/track->GetTPCNclsF() : 0.);
 852         fQACollection->Fill("TPCncls_Signal", track->Pt(), track->GetTPCNcls());
 853         fQACollection->Fill("TPCclr_Signal", track->Pt(), track->GetTPCNclsF() > 0 ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 0.);
 854       }
 855     }
 856
 857     // RecKine: ITSTPC cuts
 858     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 859
 860     // Check TRD criterions (outside the correction framework)
 861     if(track->GetTRDncls()){
 862       fQACollection->Fill("chi2TRD", track->GetTRDchi2()/track->GetTRDncls());
 863       fQACollection->Fill("alpha_rec", track->GetAlpha());
 864       fQACollection->Fill("pidquality", container[0], track->GetTRDpidQuality());
 865       fQACollection->Fill("ntrdclusters", container[0], track->GetTRDncls());
 866     }
 867
 868
 869     // RecPrim
 870     if(track->GetKinkIndex(0) != 0) continue; // Quick and dirty fix to reject both kink mothers and daughters
 871     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 872
 873     // HFEcuts: ITS layers cuts
 874     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 875
 876     if(HasMCData() && IsQAOn(kMCqa)) {
 877       // mc qa for after the reconstruction cuts
 878       AliDebug(2, "Running MC QA");
 879       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 3);  // charm
 880       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 3); // beauty
 881       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kOthers,  AliHFEmcQA::kElectronPDG, 3); // beauty
 882     }
 883
 884     // HFEcuts: Nb of tracklets TRD0
 885     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
 886
 887     // Fill correlation maps before PID
 888     if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
 889       //printf("Fill correlation maps before PID\n");
 890       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
 891     }
 892
 893     if (HasMCData() && IsQAOn(kMCqa)) {
 894       // mc qa for after the reconstruction and pid cuts
 895       AliDebug(2, "Running MC QA");
 896       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 4);  // charm
 897       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 4); // beauty
 898       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kOthers,  AliHFEmcQA::kElectronPDG, 4); // beauty
 899     }
 900
 901     if(HasMCData()){
 902       FillProductionVertex(track);
 903     }
 904
 905     // track accepted, do PID
 906     AliHFEpidObject hfetrack;
 907     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 908     hfetrack.SetRecTrack(track);
 909     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
 910     hfetrack.SetCentrality(fCentralityF);
 911     fPID->SetVarManager(fVarManager);
 912     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
 913     nElectronCandidates++;
 914     fQACollection->Fill("TPCclusters2_1_Selected", track->Pt(), track->GetTPCClusterInfo(2,1));
 915     fQACollection->Fill("TPCclusters2_0_Selected", track->Pt(), track->GetTPCClusterInfo(2,0));
 916
 917     // Fill Histogram for Hadronic Background
 918     if(HasMCData()){
 919       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11))
 920         fVarManager->FillContainer(fContainer, "hadronicBackground", UInt_t(0), kFALSE);
 921     }
 922
 923     // Fill Containers
 924     if(signal) {
 925       // Apply weight for background contamination
 926             if(fBackGroundFactorApply==kTRUE) {
 927               if(IsPbPb()) fWeightBackGround =  fBackGroundFactorArray[(Int_t)fCentralityF]->Eval(TMath::Abs(track->P()));
 928               else    fWeightBackGround =  fBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
 929
 930               if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
 931               else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
 932         // weightBackGround as special weight
 933         fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
 934       }
 935       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
 936     }
 937
 938     if(GetPlugin(kSecVtx)) {
 939       AliDebug(2, "Running Secondary Vertex Analysis");
 940       if(fSecVtx->Process(track) && signal) {
 941         fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
 942         fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
 943       }
 944     }
 945
 946     if(HasMCData()){
 947       dataE[0] = track->Pt();
 948       dataE[1] = track->Eta();
 949       dataE[2] = track->Phi();
 950       dataE[3] = track->Charge();
 951       dataE[4] = -1.;
 952       dataE[5] = -1.;
 953
 954       // Track selected: distinguish between true and fake
 955       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
 956       if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
 957         Int_t type = 0;
 958         if(fSignalCuts->IsCharmElectron(track))
 959           type = 1;
 960         else if(fSignalCuts->IsBeautyElectron(track))
 961           type = 2;
 962         AliDebug(1, Form("Type: %d\n", type));
 963         if(type){
 964                 dataE[5] = type; // beauty[1] or charm[2]
 965                 dataE[4] = 2;  // signal electron
 966         }
 967         else{
 968                 dataE[4] = 1; // not a signal electron
 969                 dataE[5] = 0;
 970         }
 971       }
 972       else {
 973         // Fill THnSparse with the information for Fake Electrons
 974         dataE[4] = 0;
 975         dataE[5] = 0;
 976       }
 977       // fill the performance THnSparse, if the mc origin could be defined
 978       if(dataE[4] > -1){
 979         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
 980         fQACollection->Fill("PIDperformance", dataE);
 981       }
 982     }
 983     // Electron background analysis
 984     if (GetPlugin(kIsElecBackGround)) {
 985
 986       AliDebug(2, "Running BackGround Analysis");
 987
 988       for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
 989               htrack = fESD->GetTrack(jtrack);
 990         if ( itrack == jtrack ) continue;
 991         fElecBackGround->PairAnalysis(track, htrack);
 992       }
 993     } // end of electron background analysis
 994
 995     if (GetPlugin(kDEstep)) {
 996       if(HasMCData()){
 997         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
 998           fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
 999           fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
1000         }
1001       }
1002       // Fill Containers for impact parameter analysis
1003       if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
1004
1005       if(signal) {
1006         fVarManager->FillContainer(fContainer, "recTrackContDEReco", AliHFEcuts::kStepHFEcutsDca, kFALSE);
1007         fVarManager->FillContainer(fContainer, "recTrackContDEMC", AliHFEcuts::kStepHFEcutsDca, kTRUE);
1008         fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterDE"));
1009       }
1010       if(HasMCData()){
1011         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1012           fQACollection->Fill("hadronsAfterIPcut",track->Pt());
1013           fQACollection->Fill("hadronsAfterIPcutMC",mctrack->Pt());
1014         }
1015       }
1016     }
1017
1018   }
1019   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1020 }
1021
1022 //____________________________________________________________
1023 void AliAnalysisTaskHFE::ProcessAOD(){
1024   //
1025   // Run Analysis in AOD Mode
1026   // Function is still in development
1027   //
1028   AliDebug(3, "Processing AOD Event");
1029   Double_t eventContainer[4];
1030   if(HasMCData()) eventContainer[0] = fVz;
1031   else {
1032     eventContainer[0] = fInputEvent->GetPrimaryVertex()->GetZ();
1033   }
1034   eventContainer[1] = 1.; // No Information available in AOD analysis, assume all events have V0AND
1035   eventContainer[2] = fCentralityF;
1036   eventContainer[3] = fContributors;
1037
1038   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1039   if(!fAOD){
1040     AliError("AOD Event required for AOD Analysis")
1041       return;
1042   }
1043
1044   //
1045   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
1046
1047   //
1048   if(fIdentifiedAsPileUp) return;
1049   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
1050
1051   //
1052   if(fIdentifiedAsOutInz) return;
1053   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
1054
1055   //
1056   if(!fPassTheEventCut) return;
1057   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
1058
1059   fContainer->NewEvent();
1060
1061   AliAODTrack *track = NULL;
1062   AliAODMCParticle *mctrack = NULL;
1063   Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
1064   Int_t nElectronCandidates = 0;
1065   Int_t pid;
1066   Bool_t signal;
1067   for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
1068     track = fAOD->GetTrack(itrack);
1069     if(!track) continue;
1070     if(track->GetFlags() != 1<<4) continue;  // Only process AOD tracks where the HFE is set
1071
1072     signal = kTRUE;
1073     if(HasMCData()){
1074
1075       Int_t label = TMath::Abs(track->GetLabel());
1076       if(label)
1077         mctrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(label));
1078         if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
1079     }
1080     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, kTRUE);
1081     // track accepted, do PID
1082     AliHFEpidObject hfetrack;
1083     hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1084     hfetrack.SetRecTrack(track);
1085     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1086     hfetrack.SetCentrality(fCentralityF);
1087     fPID->SetVarManager(fVarManager);
1088     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;    // we will do PID here as soon as possible
1089     // Apply weight for background contamination
1090     Double_t weightBackGround = 1.0;
1091
1092     // not correct treatment for pp
1093     if(fBackGroundFactorApply==kTRUE) {
1094             if(IsPbPb()) weightBackGround =  fBackGroundFactorArray[(Int_t)fCentralityF]->Eval(TMath::Abs(track->P()));
1095       else weightBackGround =  fBackGroundFactorArray[0]->Eval(TMath::Abs(track->P()));
1096
1097       if(weightBackGround < 0.0) weightBackGround = 0.0;
1098             else if(weightBackGround > 1.0) weightBackGround = 1.0;
1099             fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, weightBackGround);
1100     }
1101
1102
1103
1104     nElectronCandidates++;
1105     if(HasMCData()){
1106       dataE[0] = track->Pt();
1107       dataE[1] = track->Eta();
1108       dataE[2] = track->Phi();
1109       dataE[3] = track->Charge();
1110       dataE[4] = -1;
1111       dataE[5] = -1;
1112       // Track selected: distinguish between true and fake
1113       // COVERITY: missing test if mctrack != 0
1114       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->GetPdgCode()));
1115       if((pid = TMath::Abs(mctrack->GetPdgCode())) == 11){
1116
1117         Int_t type = 0;
1118         if(fSignalCuts->IsCharmElectron(track))
1119           type = 1;
1120         else if(fSignalCuts->IsBeautyElectron(track))
1121           type = 2;
1122         AliDebug(1, Form("Type: %d\n", type));
1123         if(type){
1124                 dataE[5] = type; // beauty[1] or charm[2]
1125                 dataE[4] = 2;  // signal electron
1126         }
1127         else{
1128                 dataE[4] = 1; // not a signal electron
1129                 dataE[5] = 0;
1130         }
1131       }
1132       else {
1133         // Fill THnSparse with the information for Fake Electrons
1134         dataE[4] = 0;
1135         dataE[5] = 0;
1136       }
1137       // fill the performance THnSparse, if the mc origin could be defined
1138       if(dataE[4] > -1){
1139         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
1140         fQACollection->Fill("PIDperformance", dataE);
1141       }
1142     }
1143   }
1144   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1145 }
1146
1147 //____________________________________________________________
1148 Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
1149   //
1150   // Filter the Monte Carlo Track
1151   // Additionally Fill a THnSparse for Signal To Background Studies
1152   // Works for AOD and MC analysis Type
1153   //
1154   fVarManager->NewTrack(track, NULL, fCentralityF, -1, kTRUE);
1155   Double_t signalContainer[6];
1156
1157   signalContainer[0] = track->Pt();
1158   signalContainer[1] = track->Eta();
1159   signalContainer[2] = track->Phi();
1160   signalContainer[3] = track->Charge()/3;
1161
1162  Double_t vertex[3]; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
1163   if(IsESDanalysis()){
1164     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
1165     if(mctrack){
1166       vertex[0] = mctrack->Particle()->Vx();
1167       vertex[1] = mctrack->Particle()->Vy();
1168     }
1169   } else {
1170     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
1171     if(aodmctrack) aodmctrack->XvYvZv(vertex);
1172   }
1173
1174   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
1175   fQACollection->Fill("mccharge", signalContainer[3]);
1176   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
1177   signalContainer[4] = 0;
1178   if(fSignalCuts->IsSelected(track)){
1179     //fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCsignal, kFALSE);
1180     // Filling of the Signal/Background histogram using the
1181     // definition of the codes for charm and beauty as below in
1182     // th crearion of the histogram
1183     if(fSignalCuts->IsCharmElectron(track))
1184       signalContainer[4] = 1;
1185     else
1186       signalContainer[4] = 2;
1187   } else {
1188     signalContainer[4] = 0; // (and other background)
1189   }
1190   signalContainer[5] = 0;
1191   // apply cut on the sqrt of the production vertex
1192   Double_t radVertex = TMath::Sqrt(vertex[0]*vertex[0] + vertex[1] * vertex[1]);
1193   if(radVertex < 3.5){
1194     // Within first ITS layer(2) -> Background we cannot reject by ITS cut, let it pass
1195     signalContainer[5] = 1;
1196   } else if (radVertex < 7.5){
1197     signalContainer[5] = 2;
1198   }
1199   fQACollection->Fill("SignalToBackgroundMC", signalContainer);
1200   fQACollection->Fill("alpha_sim", track->Phi() - TMath::Pi());
1201
1202   // Step GeneratedZOutNoPileUp
1203   if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01)) return kFALSE;
1204   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
1205   //printf("In ProcessMCtrack %f\n",fCentralityF);
1206
1207   // Step Generated Event Cut
1208   if(!fPassTheEventCut) return kFALSE;
1209   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedEventCut, kFALSE);
1210
1211   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
1212   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCInAcceptance, kFALSE);
1213   return kTRUE;
1214 }
1215
1216 //____________________________________________________________
1217 Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
1218   //
1219   // Preselect tracks
1220   //
1221
1222
1223   Bool_t survived = kTRUE;
1224
1225   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) {
1226     survived = kFALSE;
1227     //printf("Did not pass AliHFEcuts::kStepRecKineITSTPC\n");
1228   }
1229   //else printf("Pass AliHFEcuts::kStepRecKineITSTPC\n");
1230   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) {
1231     survived = kFALSE;
1232     //printf("Did not pass AliHFEcuts::kStepRecPrim\n");
1233   }
1234   //else printf("Pass AliHFEcuts::kStepRecPrim\n");
1235   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) {
1236     survived = kFALSE;
1237     //printf("Did not pass AliHFEcuts::kStepHFEcutsITS\n");
1238   }
1239   //else printf("Pass AliHFEcuts::kStepHFEcutsITS\n");
1240   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) {
1241     survived = kFALSE;
1242     //printf("Did not pass AliHFEcuts::kStepHFEcutsTRD\n");
1243   }
1244   //else printf("Pass AliHFEcuts::kStepHFEcutsTRD\n");
1245
1246   if(survived){
1247     // Apply PID
1248     AliHFEpidObject hfetrack;
1249     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1250     hfetrack.SetRecTrack(track);
1251     if(!fPIDpreselect->IsSelected(&hfetrack)) {
1252       //printf("Did not pass AliHFEcuts::kPID\n");
1253       survived = kFALSE;
1254     }
1255     //else printf("Pass AliHFEcuts::kPID\n");
1256   }
1257
1258   return survived;
1259
1260 }
1261 //____________________________________________________________
1262 void AliAnalysisTaskHFE::MakeEventContainer(){
1263   //
1264   // Create the event container for the correction framework and link it
1265   // 1st bin: Vertex z-position
1266   // 2nd bin: V0AND decision (normalization to sigma_inel)
1267   // 3rd bin: Centrality class (for pp defined as number of contributors in vertex.)
1268   //
1269
1270   if(IsPbPb()) {
1271
1272     //printf("This is PbPb!!!!!!!!!!!\n");
1273
1274     const Int_t kNvar = 4;  // number of variables on the grid:
1275     Int_t nBins[kNvar] = {120, 2, 11, 2};
1276     Double_t binMin[kNvar] = {-30. , 0., 0.0, 0.};
1277     Double_t binMax[kNvar] = {30., 2., 11.0, 2.};
1278
1279     AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
1280
1281     Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
1282     Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
1283     Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
1284     Double_t *contributorsBins = AliHFEtools::MakeLinearBinning(nBins[3], binMin[3], binMax[3]);
1285     evCont->SetBinLimits(0, vertexBins);
1286     evCont->SetBinLimits(1, v0andBins);
1287     evCont->SetBinLimits(2, centralityBins);
1288     evCont->SetBinLimits(3, contributorsBins);
1289     delete[] vertexBins; delete[] v0andBins; delete[] centralityBins; delete[] contributorsBins;
1290
1291     fCFM->SetEventContainer(evCont);
1292   }
1293   else {
1294
1295     //printf("This is pp!!!!!!!!!!!\n");
1296
1297     const Int_t kNvar = 3;  // number of variables on the grid:
1298     Int_t nBins[kNvar] = {120, 2, 11};
1299     Double_t binMin[kNvar] = {-30. , 0., 0.0};
1300     Double_t binMax[kNvar] = {30., 2., 11.0};
1301
1302     AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
1303
1304     Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
1305     Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
1306     Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
1307     evCont->SetBinLimits(0, vertexBins);
1308     evCont->SetBinLimits(1, v0andBins);
1309     evCont->SetBinLimits(2, centralityBins);
1310     delete[] vertexBins; delete[] v0andBins; delete[] centralityBins;
1311
1312     fCFM->SetEventContainer(evCont);
1313   }
1314
1315 }
1316
1317 //____________________________________________________________
1318 void AliAnalysisTaskHFE::MakeParticleContainer(){
1319   //
1320   // Create the particle container for the correction framework manager and
1321   // link it
1322   //
1323
1324   if(!fContainer) fContainer = new AliHFEcontainer("trackContainer");
1325   fVarManager->DefineVariables(fContainer);
1326
1327   // Create Correction Framework containers
1328   fContainer->CreateContainer("MCTrackCont", "Track Container filled with MC information", AliHFEcuts::kNcutStepsMCTrack);
1329   fContainer->CreateContainer("recTrackContReco", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1330   fContainer->CreateContainer("recTrackContMC", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1331
1332   fContainer->CreateContainer("hadronicBackground", "Container for Hadronic Background", 2);
1333   fContainer->CreateContainer("recTrackContDEReco", "Container for displaced electron analysis with Reco information", 1);
1334   fContainer->CreateContainer("recTrackContDEMC", "Container for displaced electron analysis with MC information", 1);
1335   fContainer->CreateContainer("recTrackContSecvtxReco", "Container for secondary vertexing analysis with Reco information", 1);
1336   fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1);
1337   fContainer->CreateCorrelationMatrix("correlationstepafterPID","THnSparse with correlations");
1338   fContainer->CreateCorrelationMatrix("correlationstepafterDE","THnSparse with correlations");
1339   if(!fVarManager->IsVariableDefined("centrality")) {
1340     //printf("Create the two other correlation maps\n");
1341     fContainer->CreateCorrelationMatrix("correlationstepbeforePID","THnSparse with correlations");
1342     fContainer->CreateCorrelationMatrix("correlationstepafterTOF","THnSparse with correlations");
1343   }
1344
1345   // Define the step names
1346   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsMCTrack; istep++){
1347     fContainer->SetStepTitle("MCTrackCont", AliHFEcuts::MCCutName(istep), istep);
1348   }
1349   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsRecTrack; istep++){
1350     fContainer->SetStepTitle("recTrackContReco", AliHFEcuts::RecoCutName(istep), istep);
1351     fContainer->SetStepTitle("recTrackContMC", AliHFEcuts::RecoCutName(istep), istep);
1352   }
1353   for(UInt_t ipid = 0; ipid < fPID->GetNumberOfPIDdetectors(); ipid++){
1354     fContainer->SetStepTitle("recTrackContReco", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1355     fContainer->SetStepTitle("recTrackContMC", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1356   }
1357 }
1358
1359 //____________________________________________________________
1360 void AliAnalysisTaskHFE::InitPIDperformanceQA(){
1361   // Add a histogram for Fake electrons
1362   const Int_t nDim=6;
1363   Int_t nBin[nDim] = {40, 8, 18, 2, 3, 3};
1364   //number of variables on the grid:pt,eta,phi,charge,
1365   const Double_t kPtbound[2] = {0.1, 20.};
1366   const Double_t kEtabound[2] = {-0.8, 0.8};
1367   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
1368   const Double_t kChargebound[2] = {-1.1, 1.1};
1369   const Double_t kAddInf1bound[2] = {0., 3.};
1370   const Double_t kAddInf2bound[2] = {0., 3.};
1371   Double_t minima[nDim] = {kPtbound[0], kEtabound[0], kPhibound[0], kChargebound[0], kAddInf1bound[0], kAddInf2bound[0]};
1372   Double_t maxima[nDim] = {kPtbound[1], kEtabound[1], kPhibound[1], kChargebound[1], kAddInf1bound[1], kAddInf2bound[1]};
1373
1374   fQACollection->CreateTHnSparse("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin, minima, maxima);
1375   fQACollection->CreateTHnSparse("SignalToBackgroundMC", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; flavor (0 - no, 1 - charm, 2 - bottom); ITS Cluster (0 - no, 1 - first (and maybe second), 2 - second)", nDim, nBin, minima, maxima);
1376
1377   fQACollection->BinLogAxis("PIDperformance", 0);
1378   fQACollection->BinLogAxis("SignalToBackgroundMC", 0);
1379   fQACollection->Sumw2("PIDperformance");
1380   fQACollection->Sumw2("SignalToBackgroundMC");
1381 }
1382
1383 //____________________________________________________________
1384 void AliAnalysisTaskHFE::InitContaminationQA(){
1385   //
1386   // Add QA for Impact Parameter cut
1387   //
1388   const Double_t kPtbound[2] = {0.1, 20.};
1389   Int_t iBin[1];
1390   iBin[0] = 44; // bins in pt
1391   fQACollection->CreateTH1F("hadronsBeforeIPcut", "Hadrons before IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
1392   fQACollection->CreateTH1F("hadronsAfterIPcut", "Hadrons after IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
1393   fQACollection->CreateTH1F("hadronsBeforeIPcutMC", "Hadrons before IP cut: MC p_{t}", iBin[0], kPtbound[0], kPtbound[1], 1);
1394   fQACollection->CreateTH1F("hadronsAfterIPcutMC", "Hadrons after IP cut: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
1395 }
1396
1397 //____________________________________________________________
1398 void AliAnalysisTaskHFE::PrintStatus() const {
1399   //
1400   // Print Analysis status
1401   //
1402   printf("\n\tAnalysis Settings\n\t========================================\n\n");
1403   printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
1404   printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
1405   printf("\tDisplaced electron analysis step: %s\n", GetPlugin(kDEstep) ? "YES" : "NO");
1406   printf("\tTagged Track Analysis: %s\n", GetPlugin(kTaggedTrackAnalysis) ? "YES" : "NO");
1407   printf("\n");
1408   printf("\tParticle Identification Detectors:\n");
1409   fPID->PrintStatus();
1410   printf("\n");
1411   printf("\tQA: \n");
1412   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
1413   printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsQAOn()) ? "YES" : "NO");
1414   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
1415   printf("\n");
1416 }
1417
1418 //____________________________________________________________
1419 Bool_t AliAnalysisTaskHFE::FillProductionVertex(const AliVParticle * const track) const{
1420   //
1421   // Find the production vertex of the associated MC track
1422   //
1423   if(!fMCEvent) return kFALSE;
1424   const AliVParticle *mctrack = NULL;
1425   TString objectType = track->IsA()->GetName();
1426   if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
1427     // Reconstructed track
1428     mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
1429   } else {
1430     // MCParticle
1431     mctrack = track;
1432   }
1433
1434   if(!mctrack) return kFALSE;
1435
1436   Double_t xv = 0.0;
1437   Double_t yv = 0.0;
1438
1439   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
1440     // case MCParticle
1441     const AliMCParticle *mcpart = dynamic_cast<const AliMCParticle *>(mctrack);
1442     if(mcpart){
1443       xv =  mcpart->Xv();
1444       yv =  mcpart->Yv();
1445     }
1446   } else {
1447     // case AODMCParticle
1448     const AliAODMCParticle *mcpart = dynamic_cast<const AliAODMCParticle *>(mctrack);
1449     if(mcpart){
1450       xv =  mcpart->Xv();
1451       yv =  mcpart->Yv();
1452     }
1453   }
1454
1455   //printf("xv %f, yv %f\n",xv,yv);
1456   fQACollection->Fill("radius", TMath::Abs(xv),TMath::Abs(yv));
1457
1458   return kTRUE;
1459
1460 }
1461 //__________________________________________
1462 void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
1463   //
1464   // Switch on Plugin
1465   // Available:
1466   //  - Primary vertex studies
1467   //  - Secondary vertex Studies
1468   //  - Post Processing
1469   //
1470   switch(plug){
1471     case kPriVtx: SETBIT(fPlugins, plug); break;
1472     case kSecVtx: SETBIT(fPlugins, plug); break;
1473     case kIsElecBackGround: SETBIT(fPlugins, plug); break;
1474     case kPostProcess: SETBIT(fPlugins, plug); break;
1475     case kDEstep: SETBIT(fPlugins, plug); break;
1476     case kTaggedTrackAnalysis: SETBIT(fPlugins, plug); break;
1477     default: AliError("Unknown Plugin");
1478   };
1479 }
1480 //__________________________________________
1481 Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track){
1482   //
1483   // Check single track cuts for a given cut step
1484   // Fill the particle container
1485   //
1486   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1487   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1488   if(fVarManager->IsSignalTrack()) {
1489     fVarManager->FillContainer(fContainer, "recTrackContReco", cutStep, kFALSE);
1490     fVarManager->FillContainer(fContainer, "recTrackContMC", cutStep, kTRUE);
1491   }
1492   return kTRUE;
1493 }
1494 //___________________________________________________
1495 Bool_t AliAnalysisTaskHFE::ReadCentrality() {
1496   //
1497   // Recover the centrality of the event from ESD or AOD
1498   //
1499  if(IsAODanalysis()){
1500
1501    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1502    if(!fAOD){
1503      AliError("AOD Event required for AOD Analysis")
1504        return kFALSE;
1505    }
1506
1507    if(IsPbPb()) {
1508      // Centrality
1509      AliCentrality *aodCentrality = fAOD->GetCentrality();
1510      Float_t fCentralityFtemp = aodCentrality->GetCentralityPercentile("V0M");
1511
1512      if( fCentralityFtemp >=  0. && fCentralityFtemp <  10.) fCentralityF =  0;
1513      else if ( fCentralityFtemp >=  10. && fCentralityFtemp <  20.) fCentralityF =  1;
1514      else if ( fCentralityFtemp >= 20. && fCentralityFtemp <  30.) fCentralityF = 2;
1515      else if ( fCentralityFtemp >= 30. && fCentralityFtemp <  40.) fCentralityF = 3;
1516      else if ( fCentralityFtemp >= 40. && fCentralityFtemp <  50.) fCentralityF = 4;
1517      else if ( fCentralityFtemp >= 50. && fCentralityFtemp <  60.) fCentralityF = 5;
1518      else if ( fCentralityFtemp >= 60. && fCentralityFtemp <  90.) fCentralityF = 6;
1519      else if ( fCentralityFtemp >= 90. && fCentralityFtemp <=  100.) fCentralityF = 7;
1520      //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <  95.) fCentralityF = 8;
1521      //else if ( fCentralityF_temp >= 95. && fCentralityF_temp <  90.) fCentralityF = 9;
1522      //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
1523      else return kFALSE;
1524
1525      // contributors
1526      fContributors = 0.5;
1527      Int_t contributorstemp = 0;
1528      const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
1529      if(vtxAOD) contributorstemp = vtxAOD->GetNContributors();
1530
1531      //printf("PbPb contributors_temp %d\n",contributors_temp);
1532
1533      if( contributorstemp <=  0) fContributors =  0.5;
1534      else fContributors = 1.5;
1535
1536
1537
1538    }
1539    else {
1540      fCentralityF = 0;
1541      Int_t centralityFtemp = 0;
1542      const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
1543      if(vtxAOD) centralityFtemp = vtxAOD->GetNContributors();
1544
1545      //printf("pp centralityF_temp %d\n",centralityF_temp);
1546
1547      if( centralityFtemp <=  0) fCentralityF =  0;
1548      else if ( centralityFtemp >  0 && centralityFtemp <  2) fCentralityF = 1;
1549      else if ( centralityFtemp >=  2 && centralityFtemp <  3) fCentralityF = 2;
1550      else if ( centralityFtemp >= 3 && centralityFtemp <  4) fCentralityF = 3;
1551      else if ( centralityFtemp >= 4 && centralityFtemp <  5) fCentralityF = 4;
1552      else if ( centralityFtemp >= 5 && centralityFtemp <  10) fCentralityF = 5;
1553      else if ( centralityFtemp >= 10 && centralityFtemp <  20) fCentralityF = 6;
1554      else if ( centralityFtemp >= 20 && centralityFtemp <  30) fCentralityF = 7;
1555      else if ( centralityFtemp >= 30 && centralityFtemp <  40) fCentralityF = 8;
1556      else if ( centralityFtemp >= 40 && centralityFtemp <  50) fCentralityF = 9;
1557      else if ( centralityFtemp >= 50) fCentralityF = 10;
1558
1559    }
1560
1561    return kTRUE;
1562
1563
1564  } else {
1565
1566    AliDebug(3, "Processing ESD Centrality");
1567    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
1568    if(!fESD){
1569      AliError("ESD Event required for ESD Analysis")
1570        return kFALSE;
1571    }
1572    const char* type = fESD->GetBeamType();
1573
1574    if (strstr(type,"Pb-Pb")) {
1575
1576    // Centrality
1577    AliCentrality *esdCentrality = fESD->GetCentrality();
1578    Float_t fCentralityFtemp = esdCentrality->GetCentralityPercentile("V0M");
1579    //printf("PbPb fCentralityF_temp %f\n",fCentralityF_temp);
1580
1581    if( fCentralityFtemp >=  0. && fCentralityFtemp <   10.) fCentralityF =  0;
1582    else if ( fCentralityFtemp >= 10. && fCentralityFtemp <  20.) fCentralityF =  1;
1583    else if ( fCentralityFtemp >= 20. && fCentralityFtemp <  30.) fCentralityF = 2;
1584    else if ( fCentralityFtemp >= 30. && fCentralityFtemp <  40.) fCentralityF = 3;
1585    else if ( fCentralityFtemp >= 40. && fCentralityFtemp <  50.) fCentralityF = 4;
1586    else if ( fCentralityFtemp >= 50. && fCentralityFtemp <  60.) fCentralityF = 5;
1587    else if ( fCentralityFtemp >= 60. && fCentralityFtemp <  90.) fCentralityF = 6;
1588    else if ( fCentralityFtemp >= 90. && fCentralityFtemp <=  100.) fCentralityF = 7;
1589    //else if ( fCentralityF_temp >= 70. && fCentralityF_temp <  80.) fCentralityF = 8;
1590    //else if ( fCentralityF_temp >= 80. && fCentralityF_temp <  90.) fCentralityF = 9;
1591    //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
1592    else return kFALSE;
1593
1594    //   Float_t fCentralityF_temp10 = esdCentrality->GetCentralityClass10("V0M");
1595    //   printf("PbPb fCentralityF_temp %f %f %f \n",fCentralityF_temp, fCentralityF_temp10, fCentralityF);
1596
1597    // contributors
1598    fContributors = 0.5;
1599    Int_t contributorstemp = 0;
1600    const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
1601    if(vtxESD && vtxESD->GetStatus()) contributorstemp = vtxESD->GetNContributors();
1602
1603    //printf("PbPb contributors_temp %d\n",contributors_temp);
1604
1605    if( contributorstemp <=  0) fContributors =  0.5;
1606    else fContributors = 1.5;
1607
1608    return kTRUE;
1609
1610    }
1611
1612
1613    if (strstr(type,"p-p")) {
1614      fCentralityF = 0;
1615      Int_t centralityFtemp = 0;
1616      const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
1617      if(vtxESD && vtxESD->GetStatus()) centralityFtemp = vtxESD->GetNContributors();
1618
1619      //printf("pp centralityF_temp %d\n",centralityF_temp);
1620
1621      if( centralityFtemp <=  0) fCentralityF =  0;
1622      else if ( centralityFtemp >  0 && centralityFtemp <  2) fCentralityF = 1;
1623      else if ( centralityFtemp >=  2 && centralityFtemp <  3) fCentralityF = 2;
1624      else if ( centralityFtemp >= 3 && centralityFtemp <  4) fCentralityF = 3;
1625      else if ( centralityFtemp >= 4 && centralityFtemp <  5) fCentralityF = 4;
1626      else if ( centralityFtemp >= 5 && centralityFtemp <  10) fCentralityF = 5;
1627      else if ( centralityFtemp >= 10 && centralityFtemp <  20) fCentralityF = 6;
1628      else if ( centralityFtemp >= 20 && centralityFtemp <  30) fCentralityF = 7;
1629      else if ( centralityFtemp >= 30 && centralityFtemp <  40) fCentralityF = 8;
1630      else if ( centralityFtemp >= 40 && centralityFtemp <  50) fCentralityF = 9;
1631      else if ( centralityFtemp >= 50) fCentralityF = 10;
1632
1633      return kTRUE;
1634
1635    }
1636
1637    return kFALSE;
1638
1639   //printf("centrality %f\n",fCentralityF);
1640
1641  }
1642
1643  //printf("centrality %f\n",fCentralityF);
1644
1645 }
1646 //___________________________________________________
1647 void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
1648   //
1649   // Recover the centrality of the event from ESD or AOD
1650   //
1651  if(IsAODanalysis()){
1652
1653    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1654    if(!fAOD){
1655      AliError("AOD Event required for AOD Analysis")
1656        return;
1657    }
1658    // PileUp
1659    if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
1660    // Z vertex
1661    if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
1662    // Event Cut
1663    fPassTheEventCut = kTRUE;
1664    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) fPassTheEventCut = kFALSE;
1665
1666
1667  } else {
1668
1669    AliDebug(3, "Processing ESD Centrality");
1670    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
1671    if(!fESD){
1672      AliError("ESD Event required for ESD Analysis")
1673        return;
1674    }
1675    // PileUp
1676    fIdentifiedAsPileUp = kFALSE;
1677    if(fRemovePileUp && fESD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
1678    // Z vertex
1679    fIdentifiedAsOutInz = kFALSE;
1680    if(fESD->GetPrimaryVertexTracks()){
1681        if(TMath::Abs(fESD->GetPrimaryVertexTracks()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
1682    }
1683    //Event Cut
1684    fPassTheEventCut = kTRUE;
1685    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) fPassTheEventCut = kFALSE;
1686
1687
1688  }
1689
1690 }
1691