PWGHF/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 <TBits.h>
  29 #include <TCanvas.h>
  30 #include <TChain.h>
  31 #include <TDirectory.h>
  32 #include <TFile.h>
  33 #include <TH3D.h>
  34 #include <TIterator.h>
  35 #include <TList.h>
  36 #include <TLegend.h>
  37 #include <TMath.h>
  38 #include <TObjArray.h>
  39 #include <TObjString.h>
  40 #include <TParticle.h>
  41 #include <TProfile.h>
  42 #include <TString.h>
  43 #include <TF1.h>
  44 #include <TTree.h>
  45
  46 #include "AliAnalysisManager.h"
  47 #include "AliAODInputHandler.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAODTrack.h"
  50 #include "AliAODVertex.h"
  51 #include "AliCentrality.h"
  52 #include "AliCFContainer.h"
  53 #include "AliCFManager.h"
  54 #include "AliESDEvent.h"
  55 #include "AliESDInputHandler.h"
  56 #include "AliESDtrack.h"
  57 #include "AliLog.h"
  58 #include "AliMCEvent.h"
  59 #include "AliMCEventHandler.h"
  60 #include "AliMCParticle.h"
  61 #include "AliMultiplicity.h"
  62 #include "AliPID.h"
  63 #include "AliPIDResponse.h"
  64 #include "AliOADBContainer.h"
  65 #include "AliStack.h"
  66 #include "AliTriggerAnalysis.h"
  67 #include "AliVVertex.h"
  68 #include "TTreeStream.h"
  69 #include "AliESDtrackCuts.h"
  70
  71 #include "AliHFEcollection.h"
  72 #include "AliHFEcontainer.h"
  73 #include "AliHFEcuts.h"
  74 #include "AliHFEelecbackground.h"
  75 #include "AliHFEmcQA.h"
  76 #include "AliHFEpairs.h"
  77 #include "AliHFEpid.h"
  78 #include "AliHFEpidQAmanager.h"
  79 #include "AliHFEpostAnalysis.h"
  80 #include "AliHFEsecVtxs.h"
  81 #include "AliHFEsecVtx.h"
  82 #include "AliHFEsignalCuts.h"
  83 #include "AliHFEtaggedTrackAnalysis.h"
  84 #include "AliHFEtools.h"
  85 #include "AliHFEvarManager.h"
  86 #include "AliAnalysisTaskHFE.h"
  87
  88 ClassImp(AliAnalysisTaskHFE)
  89
  90 //____________________________________________________________
  91 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  92   AliAnalysisTaskSE("PID efficiency Analysis")
  93   , fQAlevel(0)
  94   , fPlugins(0)
  95   , fFillSignalOnly(kTRUE)
  96   , fFillNoCuts(kFALSE)
  97   , fBackGroundFactorApply(kFALSE)
  98   , fRemovePileUp(kFALSE)
  99   , fIdentifiedAsPileUp(kFALSE)
 100   , fIdentifiedAsOutInz(kFALSE)
 101   , fPassTheEventCut(kFALSE)
 102   , fRejectKinkMother(kTRUE)
 103   , fisppMultiBin(kFALSE)
 104   , fisNonHFEsystematics(kFALSE)
 105   , fSpecialTrigger(NULL)
 106   , fCentralityF(-1)
 107   , fContributors(0.5)
 108   , fWeightBackGround(0.)
 109   , fVz(0.0)
 110   , fHadronBackgroundOADB(NULL)
 111   , fContainer(NULL)
 112   , fVarManager(NULL)
 113   , fSignalCuts(NULL)
 114   , fCFM(NULL)
 115   , fTriggerAnalysis(NULL)
 116   , fPID(NULL)
 117   , fPIDqa(NULL)
 118   , fPIDpreselect(NULL)
 119   , fCuts(NULL)
 120   , fTaggedTrackCuts(NULL)
 121   , fCleanTaggedTrack(kFALSE)
 122   , fVariablesTRDTaggedTrack(kFALSE)
 123   , fCutspreselect(NULL)
 124   , fSecVtx(NULL)
 125   , fElecBackGround(NULL)
 126   , fMCQA(NULL)
 127   , fTaggedTrackAnalysis(NULL)
 128   , fExtraCuts(NULL)
 129   , fQA(NULL)
 130   , fOutput(NULL)
 131   , fHistMCQA(NULL)
 132   , fHistSECVTX(NULL)
 133   , fHistELECBACKGROUND(NULL)
 134   , fQACollection(NULL)
 135   , fDebugLevel(0)
 136   , fTreeStream(NULL)
 137 {
 138   //
 139   // Dummy constructor
 140   //
 141   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 142   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
 143   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 144   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
 145 }
 146
 147 //____________________________________________________________
 148 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
 149   AliAnalysisTaskSE(name)
 150   , fQAlevel(0)
 151   , fPlugins(0)
 152   , fFillSignalOnly(kTRUE)
 153   , fFillNoCuts(kFALSE)
 154   , fBackGroundFactorApply(kFALSE)
 155   , fRemovePileUp(kFALSE)
 156   , fIdentifiedAsPileUp(kFALSE)
 157   , fIdentifiedAsOutInz(kFALSE)
 158   , fPassTheEventCut(kFALSE)
 159   , fRejectKinkMother(kTRUE)
 160   , fisppMultiBin(kFALSE)
 161   , fisNonHFEsystematics(kFALSE)
 162   , fSpecialTrigger(NULL)
 163   , fCentralityF(-1)
 164   , fContributors(0.5)
 165   , fWeightBackGround(0.)
 166   , fVz(0.0)
 167   , fHadronBackgroundOADB(NULL)
 168   , fContainer(NULL)
 169   , fVarManager(NULL)
 170   , fSignalCuts(NULL)
 171   , fCFM(NULL)
 172   , fTriggerAnalysis(NULL)
 173   , fPID(NULL)
 174   , fPIDqa(NULL)
 175   , fPIDpreselect(NULL)
 176   , fCuts(NULL)
 177   , fTaggedTrackCuts(NULL)
 178   , fCleanTaggedTrack(kFALSE)
 179   , fVariablesTRDTaggedTrack(kFALSE)
 180   , fCutspreselect(NULL)
 181   , fSecVtx(NULL)
 182   , fElecBackGround(NULL)
 183   , fMCQA(NULL)
 184   , fTaggedTrackAnalysis(NULL)
 185   , fExtraCuts(NULL)
 186   , fQA(NULL)
 187   , fOutput(NULL)
 188   , fHistMCQA(NULL)
 189   , fHistSECVTX(NULL)
 190   , fHistELECBACKGROUND(NULL)
 191   , fQACollection(0x0)
 192   , fDebugLevel(0)
 193   , fTreeStream(NULL)
 194 {
 195   //
 196   // Default constructor
 197   //
 198   DefineOutput(1, TList::Class());
 199   DefineOutput(2, TList::Class());
 200
 201   fPID = new AliHFEpid("hfePid");
 202   fPIDqa = new AliHFEpidQAmanager;
 203   fVarManager = new AliHFEvarManager("hfeVarManager");
 204
 205   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 206   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
 207   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 208   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
 209 }
 210
 211 //____________________________________________________________
 212 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
 213   AliAnalysisTaskSE(ref)
 214   , fQAlevel(0)
 215   , fPlugins(0)
 216   , fFillSignalOnly(ref.fFillSignalOnly)
 217   , fFillNoCuts(ref.fFillNoCuts)
 218   , fBackGroundFactorApply(ref.fBackGroundFactorApply)
 219   , fRemovePileUp(ref.fRemovePileUp)
 220   , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
 221   , fIdentifiedAsOutInz(ref.fIdentifiedAsOutInz)
 222   , fPassTheEventCut(ref.fPassTheEventCut)
 223   , fRejectKinkMother(ref.fRejectKinkMother)
 224   , fisppMultiBin(ref.fisppMultiBin)
 225   , fisNonHFEsystematics(ref.fisNonHFEsystematics)
 226   , fSpecialTrigger(ref.fSpecialTrigger)
 227   , fCentralityF(ref.fCentralityF)
 228   , fContributors(ref.fContributors)
 229   , fWeightBackGround(ref.fWeightBackGround)
 230   , fVz(ref.fVz)
 231   , fHadronBackgroundOADB(ref.fHadronBackgroundOADB)
 232   , fContainer(NULL)
 233   , fVarManager(NULL)
 234   , fSignalCuts(NULL)
 235   , fCFM(NULL)
 236   , fTriggerAnalysis(NULL)
 237   , fPID(NULL)
 238   , fPIDqa(NULL)
 239   , fPIDpreselect(NULL)
 240   , fCuts(NULL)
 241   , fTaggedTrackCuts(NULL)
 242   , fCleanTaggedTrack(ref.fCleanTaggedTrack)
 243   , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
 244   , fCutspreselect(NULL)
 245   , fSecVtx(NULL)
 246   , fElecBackGround(NULL)
 247   , fMCQA(NULL)
 248   , fTaggedTrackAnalysis(NULL)
 249   , fExtraCuts(NULL)
 250   , fQA(NULL)
 251   , fOutput(NULL)
 252   , fHistMCQA(NULL)
 253   , fHistSECVTX(NULL)
 254   , fHistELECBACKGROUND(NULL)
 255   , fQACollection(NULL)
 256   , fDebugLevel(ref.fDebugLevel)
 257   , fTreeStream(NULL)
 258 {
 259   //
 260   // Copy Constructor
 261   //
 262   ref.Copy(*this);
 263 }
 264
 265 //____________________________________________________________
 266 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
 267   //
 268   // Assignment operator
 269   //
 270   if(this == &ref)
 271     ref.Copy(*this);
 272   return *this;
 273 }
 274
 275 //____________________________________________________________
 276 void AliAnalysisTaskHFE::Copy(TObject &o) const {
 277   //
 278   // Copy into object o
 279   //
 280   AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
 281   target.fQAlevel = fQAlevel;
 282   target.fPlugins = fPlugins;
 283   target.fFillSignalOnly = fFillSignalOnly;
 284   target.fFillNoCuts = fFillNoCuts;
 285   target.fBackGroundFactorApply = fBackGroundFactorApply;
 286   target.fRemovePileUp = fRemovePileUp;
 287   target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
 288   target.fIdentifiedAsOutInz = fIdentifiedAsOutInz;
 289   target.fPassTheEventCut = fPassTheEventCut;
 290   target.fRejectKinkMother = fRejectKinkMother;
 291   target.fisppMultiBin =   fisppMultiBin;
 292   target.fisNonHFEsystematics = fisNonHFEsystematics;
 293   target.fSpecialTrigger = fSpecialTrigger;
 294   target.fCentralityF = fCentralityF;
 295   target.fContributors = fContributors;
 296   target.fWeightBackGround = fWeightBackGround;
 297   target.fVz = fVz;
 298   target.fHadronBackgroundOADB = fHadronBackgroundOADB;
 299   target.fContainer = fContainer;
 300   target.fVarManager = fVarManager;
 301   target.fSignalCuts = fSignalCuts;
 302   target.fCFM = fCFM;
 303   target.fTriggerAnalysis = fTriggerAnalysis;
 304   target.fPID = fPID;
 305   target.fPIDqa = fPIDqa;
 306   target.fPIDpreselect = fPIDpreselect;
 307   target.fCuts = fCuts;
 308   target.fTaggedTrackCuts = fTaggedTrackCuts;
 309   target.fCleanTaggedTrack = fCleanTaggedTrack;
 310   target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
 311   target.fCutspreselect = fCutspreselect;
 312   target.fSecVtx = fSecVtx;
 313   target.fElecBackGround = fElecBackGround;
 314   target.fMCQA = fMCQA;
 315   target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
 316   target.fExtraCuts = fExtraCuts;
 317   target.fQA = fQA;
 318   target.fOutput = fOutput;
 319   target.fHistMCQA = fHistMCQA;
 320   target.fHistSECVTX = fHistSECVTX;
 321   target.fHistELECBACKGROUND = fHistELECBACKGROUND;
 322   target.fQACollection = fQACollection;
 323   target.fDebugLevel = fDebugLevel;
 324   target.fTreeStream = fTreeStream;
 325 }
 326
 327 //____________________________________________________________
 328 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
 329   //
 330   // Destructor
 331   //
 332   if(fPID) delete fPID;
 333   if(fPIDpreselect) delete fPIDpreselect;
 334   if(fVarManager) delete fVarManager;
 335   if(fCFM) delete fCFM;
 336   if(fTriggerAnalysis) delete fTriggerAnalysis;
 337   if(fSignalCuts) delete fSignalCuts;
 338   if(fSecVtx) delete fSecVtx;
 339   if(fMCQA) delete fMCQA;
 340   if(fElecBackGround) delete fElecBackGround;
 341   if(fSpecialTrigger) delete fSpecialTrigger;
 342   // Delete output objects only if we are not running in PROOF mode because otherwise this produces a crash during merging
 343   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 344   if(mgr && mgr->GetAnalysisType() != AliAnalysisManager::kProofAnalysis){
 345     if(fPIDqa) delete fPIDqa;
 346     if(fOutput) delete fOutput;
 347     if(fQA) delete fQA;
 348     if(fTreeStream) delete fTreeStream;
 349   }
 350 }
 351
 352 //____________________________________________________________
 353 void AliAnalysisTaskHFE::UserCreateOutputObjects(){
 354   //
 355   // Creating output container and output objects
 356   // Here we also Initialize the correction framework container and
 357   // the objects for
 358   // - PID
 359   // - MC QA
 360   // - SecVtx
 361   // QA histograms are created if requested
 362   // Called once per worker
 363   //
 364   AliDebug(3, "Creating Output Objects");
 365   // Automatic determination of the analysis mode
 366   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 367   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
 368     SetAODAnalysis();
 369   } else {
 370     SetESDAnalysis();
 371     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 372       SetHasMCData();
 373   }
 374   printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
 375   printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
 376
 377   // Enable Trigger Analysis
 378   fTriggerAnalysis = new AliTriggerAnalysis;
 379   fTriggerAnalysis->EnableHistograms();
 380   fTriggerAnalysis->SetAnalyzeMC(HasMCData());
 381
 382
 383   // Make lists for Output
 384   if(!fQA) fQA = new TList;
 385   fQA->SetOwner();
 386   if(!fOutput) fOutput = new TList;
 387   fOutput->SetOwner();
 388
 389   // First Part: Make QA histograms
 390   fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
 391   fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
 392   fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
 393   fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
 394
 395   InitPIDperformanceQA();
 396   InitContaminationQA();
 397   InitHistoITScluster();
 398   fQA->Add(fQACollection);
 399
 400   // Initialize PID
 401   fPID->SetHasMCData(HasMCData());
 402   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
 403   if(IsQAOn(kPIDqa)){
 404     AliInfo("PID QA switched on");
 405     fPIDqa->Initialize(fPID);
 406     fQA->Add(fPIDqa->MakeList("HFEpidQA"));
 407   }
 408   fPID->SortDetectors();
 409
 410   // Initialize correction Framework and Cuts
 411   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack + AliHFEcuts::kNcutStepsSecvtxTrack;
 412   fCFM = new AliCFManager;
 413   fCFM->SetNStepParticle(kNcutSteps);
 414   MakeParticleContainer();
 415   MakeEventContainer();
 416   // Temporary fix: Initialize particle cuts with NULL
 417   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 418     fCFM->SetParticleCutsList(istep, NULL);
 419   if(!fCuts){
 420     AliWarning("Cuts not available. Default cuts will be used");
 421     fCuts = new AliHFEcuts;
 422     fCuts->CreateStandardCuts();
 423   }
 424   if(IsAODanalysis()) fCuts->SetAOD();
 425   // Make clone for V0 tagging step
 426   fCuts->Initialize(fCFM);
 427   if(fCuts->IsQAOn()) fQA->Add(fCuts->GetQAhistograms());
 428   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
 429   fVarManager->SetSignalCuts(fSignalCuts);
 430
 431   // add output objects to the List
 432   fOutput->AddAt(fContainer, 0);
 433   fOutput->AddAt(fCFM->GetEventContainer(), 1);
 434
 435   // mcQA----------------------------------
 436   if (HasMCData() && IsQAOn(kMCqa)) {
 437     AliInfo("MC QA on");
 438     if(!fMCQA) fMCQA = new AliHFEmcQA;
 439     if(!fHistMCQA) fHistMCQA = new TList();
 440     fHistMCQA->SetOwner();
 441     if(IsPbPb()) fMCQA->SetPbPb();
 442     if(fisppMultiBin) fMCQA->SetPPMultiBin();
 443     fMCQA->CreatDefaultHistograms(fHistMCQA);
 444     if(!fFillSignalOnly) fMCQA->SetBackgroundWeightFactor(fElecBackgroundFactor[0][0][0],fBinLimit);
 445     fQA->Add(fHistMCQA);
 446   }
 447
 448   // secvtx----------------------------------
 449   if (GetPlugin(kSecVtx)) {
 450     AliInfo("Secondary Vertex Analysis on");
 451     if(!fSecVtx) fSecVtx = new AliHFEsecVtx;
 452     fSecVtx->SetHasMCData(HasMCData());
 453
 454     if(!fHistSECVTX) fHistSECVTX = new TList();
 455     fHistSECVTX->SetOwner();
 456     fSecVtx->CreateHistograms(fHistSECVTX);
 457     fOutput->Add(fHistSECVTX);
 458   }
 459
 460   // background----------------------------------
 461   if (GetPlugin(kIsElecBackGround)) {
 462     AliInfo("Electron BackGround Analysis on");
 463     if(!fElecBackGround){
 464       AliWarning("ElecBackGround not available. Default elecbackground will be used");
 465       fElecBackGround = new AliHFEelecbackground;
 466     }
 467     fElecBackGround->SetHasMCData(HasMCData());
 468
 469     if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
 470     fHistELECBACKGROUND->SetOwner();
 471     fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
 472     fOutput->Add(fHistELECBACKGROUND);
 473   }
 474
 475   // tagged tracks
 476   if(GetPlugin(kTaggedTrackAnalysis)){
 477     AliInfo("Analysis on V0-tagged tracks enabled");
 478     fTaggedTrackAnalysis = new AliHFEtaggedTrackAnalysis(Form("taggedTrackAnalysis%s", GetName()));
 479     fTaggedTrackAnalysis->SetCuts(fTaggedTrackCuts);
 480     fTaggedTrackAnalysis->SetClean(fCleanTaggedTrack);
 481     AliHFEvarManager *varManager = fTaggedTrackAnalysis->GetVarManager();
 482     TObjArray *array = fVarManager->GetVariables();
 483     Int_t nvars = array->GetEntriesFast();
 484     TString namee;
 485     for(Int_t v = 0; v < nvars; v++) {
 486       AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
 487       if(!variable) continue;
 488       TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
 489       if(!name.CompareTo("source")) namee = TString("species");
 490       else namee = TString(name);
 491       Int_t nbins = variable->GetNumberOfBins();
 492       if(variable->HasUserDefinedBinning()){
 493         varManager->AddVariable(namee, nbins, variable->GetBinning());
 494       } else {
 495         varManager->AddVariable(namee, nbins, variable->GetMinimum(), variable->GetMaximum());
 496       }
 497       //printf("For AliTaggedTrackAnalysis, had the variable %s and the one used %s\n",(const char*)variable->GetName(),(const char*) namee);
 498     }
 499     if(fPIDqa->HasHighResolutionHistos())
 500       fTaggedTrackAnalysis->GetPIDqa()->SetHighResolutionHistos();
 501     fTaggedTrackAnalysis->SetPID(fPID);
 502     fTaggedTrackAnalysis->SetVariablesTRD(fVariablesTRDTaggedTrack);
 503     fTaggedTrackAnalysis->InitContainer();
 504     fOutput->Add(fTaggedTrackAnalysis->GetContainer());
 505     fQA->Add(fTaggedTrackAnalysis->GetPIDQA());
 506     fQA->Add(fTaggedTrackAnalysis->GetCutQA());
 507     fQA->Add(fTaggedTrackAnalysis->GetQAcollection());
 508   }
 509
 510   Bool_t isProof = AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == AliAnalysisManager::kProofAnalysis;
 511   if(fDebugLevel && !isProof){
 512     AliDebug(1,"Create OutputStream");
 513     fTreeStream = new TTreeSRedirector(Form("HFEdebugTree%s.root", GetName()));
 514   }
 515
 516   PrintStatus();
 517 }
 518
 519 //____________________________________________________________
 520 void AliAnalysisTaskHFE::UserExec(Option_t *){
 521   //
 522   // Run the analysis
 523   //
 524   AliDebug(3, "Starting Single Event Analysis");
 525   if(!fInputEvent){
 526     AliError("Reconstructed Event not available");
 527     return;
 528   }
 529   if(HasMCData()){
 530     AliDebug(4, Form("MC Event: %p", fMCEvent));
 531     if(!fMCEvent){
 532       AliError("No MC Event, but MC Data required");
 533       return;
 534     }
 535   }
 536   if(!fCuts){
 537     AliError("HFE cuts not available");
 538     return;
 539   }
 540   if(!fPID->IsInitialized()){
 541     // Initialize PID with the given run number
 542     fPID->InitializePID(fInputEvent->GetRunNumber());
 543   }
 544
 545   // Initialize hadronic background from OADB Container
 546   AliDebug(2, Form("Apply background factors: %s, OADB Container %p", fBackGroundFactorApply ? "Yes" : "No", fHadronBackgroundOADB));
 547   if(fBackGroundFactorApply && !TestBit(kBackgroundInitialized)){
 548     AliDebug(2, "Initializing Background from OADB");
 549     if(!InitializeHadronBackground(fInputEvent->GetRunNumber())) AliError("Failed initializing hadronic background parameterization from OADB");
 550     else AliDebug(2, "Successfully loaded Background from OADB");
 551     SetBit(kBackgroundInitialized);
 552   }
 553
 554   if(IsESDanalysis() && HasMCData()){
 555     // Protect against missing MC trees
 556     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 557     if(!mcH){
 558       AliError("No MC Event Handler available");
 559       return;
 560     }
 561     if(!mcH->InitOk()) return;
 562     if(!mcH->TreeK()) return;
 563     if(!mcH->TreeTR()) return;
 564   }
 565
 566   // need the centrality for everything (MC also)
 567   fCentralityF = -1;
 568   if(!ReadCentrality()) fCentralityF = -1;
 569   //printf("pass centrality\n");
 570   //printf("Reading fCentralityF %f\n",fCentralityF);
 571
 572   // See if pile up and z in the range
 573   RejectionPileUpVertexRangeEventCut();
 574
 575   // Protect agains missing
 576   if(HasMCData()){
 577     //printf("Has MC data\n");
 578     fSignalCuts->SetMCEvent(fMCEvent);
 579     ProcessMC();  // Run the MC loop + MC QA in case MC Data are available
 580   }
 581
 582   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 583   if(!pidResponse){
 584     AliDebug(1, "Using default PID Response");
 585     pidResponse = AliHFEtools::GetDefaultPID(HasMCData(), fInputEvent->IsA() == AliAODEvent::Class());
 586   }
 587   fPID->SetPIDResponse(pidResponse);
 588   if(fPIDpreselect) fPIDpreselect->SetPIDResponse(pidResponse);
 589
 590   // Event loop
 591   if(IsAODanalysis()){
 592     ProcessAOD();
 593   } else {
 594     const char *specialTrigger = GetSpecialTrigger(fInputEvent->GetRunNumber());
 595     // Check Trigger selection
 596     if(specialTrigger){
 597       AliDebug(2, Form("Special Trigger requested: %s", specialTrigger));
 598       AliESDEvent *ev = dynamic_cast<AliESDEvent *>(fInputEvent);
 599       if(!(ev && ev->IsTriggerClassFired(specialTrigger))){
 600         AliDebug(2, "Event not selected");
 601         return;
 602       } else AliDebug(2, "Event Selected");
 603     } else AliDebug(2, "No Special Trigger requested");
 604
 605     ProcessESD();
 606   }
 607   // Done!!!
 608   PostData(1, fOutput);
 609   PostData(2, fQA);
 610 }
 611
 612 //____________________________________________________________
 613 void AliAnalysisTaskHFE::Terminate(Option_t *){
 614   //
 615   // Terminate not implemented at the moment
 616   //
 617   if(GetPlugin(kPostProcess)){
 618     fOutput = dynamic_cast<TList *>(GetOutputData(1));
 619     fQA = dynamic_cast<TList *>(GetOutputData(2));
 620     if(!fOutput){
 621       AliError("Results not available");
 622       return;
 623     }
 624     if(!fQA){
 625       AliError("QA output not available");
 626       return;
 627     }
 628     fContainer = dynamic_cast<AliHFEcontainer *>(fOutput->FindObject("trackContainer"));
 629     if(!fContainer){
 630       AliError("Track container not found");
 631       return;
 632     }
 633     AliHFEpostAnalysis postanalysis;
 634     postanalysis.SetTaskResults(fContainer);
 635     TList *qalist = dynamic_cast<TList *>(fQA->FindObject("list_TaskQA"));
 636     if(!qalist){
 637       AliError("QA List not found");
 638       return;
 639     }
 640     postanalysis.SetTaskQA(qalist);
 641     printf("Running post analysis\n");
 642     //if(HasMCData())
 643     postanalysis.DrawMCSignal2Background();
 644     postanalysis.DrawEfficiency();
 645     postanalysis.DrawPIDperformance();
 646     postanalysis.DrawCutEfficiency();
 647
 648     if (GetPlugin(kIsElecBackGround)) {
 649       AliHFEelecbackground elecBackGround;
 650       TList *oe = 0x0;
 651       if(!(oe = (TList*)dynamic_cast<TList *>(fOutput->FindObject("HFEelecbackground")))){
 652         return;
 653       }
 654       elecBackGround.Load(oe);
 655       elecBackGround.Plot();
 656       elecBackGround.PostProcess();
 657     }
 658   }
 659 }
 660 //_______________________________________________________________
 661 Bool_t AliAnalysisTaskHFE::IsEventInBinZero() {
 662   //
 663   //
 664   //
 665
 666   //printf("test in IsEventInBinZero\n");
 667   if(!fInputEvent){
 668     AliError("Reconstructed Event not available");
 669     return kFALSE;
 670   }
 671
 672   // check vertex
 673   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
 674   if(!vertex) return kTRUE;
 675   //if(vertex) return kTRUE;
 676
 677   // check tracks
 678   if(fInputEvent->GetNumberOfTracks()<=0) return kTRUE;
 679   //if(fInputEvent->GetNumberOfTracks()>0) return kTRUE;
 680
 681
 682   return kFALSE;
 683
 684 }
 685 //____________________________________________________________
 686 void AliAnalysisTaskHFE::ProcessMC(){
 687   //
 688   // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
 689   // In case MC QA is on also MC QA loop is done
 690   //
 691   AliDebug(3, "Processing MC Information");
 692   Double_t eventContainer [4];
 693   eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
 694   eventContainer[2] = fCentralityF;
 695   eventContainer[3] = fContributors;
 696   fVz = eventContainer[0];
 697   //printf("z position is %f\n",eventContainer[0]);
 698   //if(fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent))
 699   fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
 700   Int_t nElectrons = 0;
 701   if(IsESDanalysis()){
 702    if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
 703     if (HasMCData() && IsQAOn(kMCqa)) {
 704       AliDebug(2, "Running MC QA");
 705
 706       if(fMCEvent->Stack()){
 707         fMCQA->SetMCEvent(fMCEvent);
 708         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
 709         fMCQA->SetCentrality(fCentralityF);
 710         if(IsPbPb()) { fMCQA->SetPbPb();}
 711         else
 712         {
 713             if(fisppMultiBin) fMCQA->SetPPMultiBin();
 714             else fMCQA->SetPP();
 715         }
 716         fMCQA->Init();
 717
 718         fMCQA->GetMesonKine();
 719
 720         // loop over all tracks for decayed electrons
 721         for (Int_t igen = 0; igen < fMCEvent->GetNumberOfTracks(); igen++){
 722           TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
 723           if(!mcpart) continue;
 724           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
 725           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
 726           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
 727           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
 728           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0); // no accept cut
 729           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 730           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 731           if (TMath::Abs(mcpart->Eta()) < 0.9) {
 732             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 733             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 734             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 735           }
 736           if (TMath::Abs(AliHFEtools::GetRapidity(mcpart)) < 0.5) {
 737             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 738             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 739             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 740           }
 741         }
 742         //fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
 743         //fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 744       }
 745
 746     } // end of MC QA loop
 747    }
 748    // -----------------------------------------------------------------
 749    fCFM->SetMCEventInfo(fMCEvent);
 750    // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 751   } else {
 752     fCFM->SetMCEventInfo(fInputEvent);
 753   }
 754   // Run MC loop
 755   AliVParticle *mctrack = NULL;
 756   AliDebug(3, Form("Number of Tracks: %d", fMCEvent->GetNumberOfTracks()));
 757   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfTracks(); imc++){
 758     if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
 759     AliDebug(4, "Next MC Track");
 760     if(ProcessMCtrack(mctrack)) nElectrons++;
 761   }
 762
 763   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 764   fQACollection->Fill("nElectron", nElectrons);
 765 }
 766
 767 //____________________________________________________________
 768 void AliAnalysisTaskHFE::ProcessESD(){
 769   //
 770   // Run Analysis of reconstructed event in ESD Mode
 771   // Loop over Tracks, filter according cut steps defined in AliHFEcuts
 772   //
 773   AliDebug(3, "Processing ESD Event");
 774   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 775   if(!fESD){
 776     AliError("ESD Event required for ESD Analysis");
 777     return;
 778   }
 779
 780   // Set magnetic field if V0 task on
 781   if(fTaggedTrackAnalysis) {
 782     fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
 783     fTaggedTrackAnalysis->SetCentrality(fCentralityF);
 784     if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
 785     else fTaggedTrackAnalysis->SetPP();
 786   }
 787
 788   // Do event Normalization
 789   Double_t eventContainer[4];
 790   eventContainer[0] = 0.;
 791   if(HasMCData()) eventContainer[0] = fVz;
 792   else {
 793     if(fESD->GetPrimaryVertexSPD()) eventContainer[0] = fESD->GetPrimaryVertexSPD()->GetZ();
 794   }
 795   eventContainer[1] = 0.;
 796   eventContainer[2] = fCentralityF;
 797   eventContainer[3] = fContributors;
 798   if(fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0AND))
 799     eventContainer[1] = 1.;
 800
 801   //
 802   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
 803
 804   //
 805   if(fIdentifiedAsPileUp) return;
 806   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
 807
 808   //
 809   if(TMath::Abs(fCentralityF) < 0) return;
 810   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
 811   //printf("In ProcessESD %f\n",fCentralityF);
 812
 813   //
 814   if(fIdentifiedAsOutInz) return;
 815   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
 816
 817   //
 818   if(!fPassTheEventCut) return;
 819   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
 820
 821
 822
 823   fContainer->NewEvent();
 824
 825   if (GetPlugin(kIsElecBackGround)) {
 826     fElecBackGround->SetEvent(fESD);
 827   }
 828   if (GetPlugin(kSecVtx)) {
 829     fSecVtx->SetEvent(fESD);
 830     fSecVtx->GetPrimaryCondition();
 831   }
 832
 833   if(HasMCData()){
 834     if (GetPlugin(kSecVtx)) {
 835       fSecVtx->SetMCEvent(fMCEvent);
 836       fSecVtx->SetMCQA(fMCQA);
 837     }
 838     if (GetPlugin(kIsElecBackGround)) {
 839       fElecBackGround->SetMCEvent(fMCEvent);
 840     }
 841   }
 842
 843   Double_t container[10];
 844   memset(container, 0, sizeof(Double_t) * 10);
 845   // container for the output THnSparse
 846   Double_t dataE[6]; // [pT, eta, Phi, type, 'C' or 'B']
 847   Int_t nElectronCandidates = 0;
 848   AliESDtrack *track = NULL, *htrack = NULL;
 849   AliMCParticle *mctrack = NULL;
 850   AliMCParticle *mctrackmother = NULL;
 851   Int_t pid = 0;
 852
 853   Bool_t signal = kTRUE;
 854
 855   fCFM->SetRecEventInfo(fESD);
 856
 857   // minjung for IP QA(temporary ~ 2weeks)
 858   if(!fExtraCuts){
 859     fExtraCuts = new AliHFEextraCuts("hfetmpCuts","HFE tmp Cuts");
 860   }
 861   fExtraCuts->SetRecEventInfo(fESD);
 862
 863   // Electron background analysis
 864   if (GetPlugin(kIsElecBackGround)) {
 865
 866     AliDebug(2, "Running BackGround Analysis");
 867
 868     fElecBackGround->Reset();
 869
 870   } // end of electron background analysis
 871   //
 872   // Loop ESD
 873   //
 874   AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
 875   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 876     AliDebug(4, "New ESD track");
 877     track = fESD->GetTrack(itrack);
 878     track->SetESDEvent(fESD);
 879
 880     // fill counts of v0-identified particles
 881     Int_t v0pid = -1;
 882     if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
 883     else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
 884     else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
 885     // here the tagged track analysis will run
 886     if(fTaggedTrackAnalysis && v0pid > -1){
 887       AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
 888       fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
 889     }
 890
 891     AliDebug(3, Form("Doing track %d, %p", itrack, track));
 892
 893     //////////////////////////////////////
 894     // preselect
 895     /////////////////////////////////////
 896     if(fPIDpreselect && fCutspreselect) {
 897       if(!PreSelectTrack(track)) continue;
 898     }
 899
 900     signal = kTRUE;
 901
 902     // Fill step without any cut
 903
 904     if(HasMCData()){
 905       // Check if it is electrons near the vertex
 906       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
 907
 908       if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
 909       else AliDebug(3, "Signal Electron");
 910
 911       // Fill K pt for Ke3 contributions
 912       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==321)) fQACollection->Fill("Kptspectra",mctrack->Pt());
 913       else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==130)) fQACollection->Fill("K0Lptspectra",mctrack->Pt());
 914     }
 915     // Cache new Track information inside the var manager
 916     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
 917
 918     if(fFillNoCuts) {
 919       if(signal || !fFillSignalOnly){
 920         fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
 921         fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
 922       }
 923     }
 924
 925     // RecKine: ITSTPC cuts
 926     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 927
 928
 929     // RecPrim
 930     if(fRejectKinkMother) {
 931       if(track->GetKinkIndex(0) != 0) continue; } // Quick and dirty fix to reject both kink mothers and daughters
 932     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 933
 934     if(fTreeStream && fDebugLevel >= 2){
 935       // Debug streaming of PID-related quantities
 936       Double_t nSigmaTOF = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 937       Double_t nSigmaTPC = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 938       if(TMath::Abs(nSigmaTOF) < 5 && TMath::Abs(nSigmaTPC) < 5){
 939         // we are not interested in tracks which are more than 5 sigma away from the electron hypothesis in either TOF or TPC
 940         Double_t charge = track->Charge() > 0 ? 1. : -1.;
 941         Char_t myv0pid = v0pid;
 942         Double_t momentum = track->P() * charge;
 943         Double_t transversemomentum = track->Pt() * charge;
 944         Int_t run = fInputEvent->GetRunNumber();
 945         Double_t eta = track->Eta();
 946         Double_t phi = track->Phi();
 947         UChar_t ntracklets = track->GetTRDntrackletsPID();
 948         UChar_t nclustersTPCPID = track->GetTPCsignalN();
 949         UChar_t nclustersTPCshared = 0;
 950         const TBits &sharedTPC = track->GetTPCSharedMap();
 951         for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
 952         UChar_t nclustersTPC = track->GetTPCncls();
 953         UChar_t nclustersITS = track->GetITSclusters(NULL);
 954         UChar_t nclustersTRD = track->GetTRDncls();
 955         UChar_t hasClusterITS[6], hasTrackletTRD[6];
 956         UChar_t itsPixel = track->GetITSClusterMap();
 957         for(Int_t icl = 0; icl < 6; icl++) hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
 958         for(Int_t itl = 0; itl < 6; itl++){
 959           Int_t nSliceNonZero = 0;
 960           for(Int_t islice = 0; islice < 8; islice++){
 961             if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
 962           }
 963           hasTrackletTRD[itl] = nSliceNonZero ? 1 : 0;
 964         }
 965         Double_t pidprobs[5];
 966         track->GetTRDpid(pidprobs);
 967         Double_t likeEleTRD = pidprobs[0];
 968         Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
 969         (*fTreeStream) << "PIDdebug"
 970               << "signal="              << signal
 971               << "v0pid="               << myv0pid
 972               << "run="                 << run
 973               << "p="                   << momentum
 974               << "pt="                  << transversemomentum
 975               << "eta="                 << eta
 976               << "phi="                 << phi
 977               << "ntracklets="          << ntracklets
 978               << "nclustersTPC="        << nclustersTPC
 979               << "nclustersTPCPID="     << nclustersTPCPID
 980               << "nclustersTPCshared="  << nclustersTPCshared
 981               << "nclustersITS="        << nclustersITS
 982               << "nclusters="           << nclustersTRD
 983               << "its0="                << hasClusterITS[0]
 984               << "its1="                << hasClusterITS[1]
 985               << "its2="                << hasClusterITS[2]
 986               << "its3="                << hasClusterITS[3]
 987               << "its4="                << hasClusterITS[4]
 988               << "its5="                << hasClusterITS[5]
 989               << "trd0="                << hasTrackletTRD[0]
 990               << "trd1="                << hasTrackletTRD[1]
 991               << "trd2="                << hasTrackletTRD[2]
 992               << "trd3="                << hasTrackletTRD[3]
 993               << "trd4="                << hasTrackletTRD[4]
 994               << "trd5="                << hasTrackletTRD[5]
 995               << "TOFsigmaEl="          << nSigmaTOF
 996               << "TPCsigmaEl="          << nSigmaTPC
 997               << "TRDlikeEl="           << likeEleTRD
 998               << "TRDlikeEln="          << likeEleTRDn
 999               << "\n";
1000       }
1001     }
1002
1003     // HFEcuts: ITS layers cuts
1004     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
1005
1006     // HFE cuts: TOF PID and mismatch flag
1007     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) continue;
1008
1009     // HFE cuts: TPC PID cleanup
1010     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
1011
1012     // HFEcuts: Nb of tracklets TRD0
1013     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
1014
1015     // Fill correlation maps before PID
1016     if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
1017       //printf("Fill correlation maps before PID\n");
1018       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
1019     }
1020
1021     if(HasMCData()){
1022       FillProductionVertex(track);
1023
1024       if(fMCQA){
1025         fMCQA->SetCentrality(fCentralityF);
1026         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){
1027          Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.};
1028          for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1029            weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1030            if(!fisNonHFEsystematics)break;
1031          }
1032
1033          if(fisNonHFEsystematics){
1034            //Fill additional containers for electron source distinction
1035            Int_t elecSource = 0;
1036            elecSource = fMCQA->GetElecSource(mctrack->Particle());
1037            const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1038            const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1039            Int_t iName = 0;
1040            for(Int_t iSource = AliHFEmcQA::kPi0; iSource <=  AliHFEmcQA::kGammaRho0; iSource++){
1041              if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1042              if(elecSource == iSource){
1043                for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1044                  if(weightElecBgV0[iLevel]>0){ fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);}
1045                  else if(weightElecBgV0[iLevel]<0){ fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);}
1046                }
1047                break;
1048              }
1049              iName++;
1050              if(iName == kElecBgSpecies)iName = 0;
1051            }
1052          }
1053          //else{
1054            if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]);
1055            else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]);
1056            //}
1057         }
1058       }
1059     }
1060
1061     AliHFEpidObject hfetrack;
1062     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1063     hfetrack.SetRecTrack(track);
1064     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1065     hfetrack.SetCentrality(fCentralityF);
1066     if(IsPbPb()) hfetrack.SetPbPb();
1067     else hfetrack.SetPP();
1068     fPID->SetVarManager(fVarManager);
1069     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
1070     nElectronCandidates++;
1071
1072     // Temporary histogram for chi2/ITS cluster
1073     if(IsPbPb()) {
1074             TBits shared = track->GetTPCSharedMap();
1075           Int_t sharebit=0;
1076             if(shared.CountBits() >= 2) sharebit=1;
1077
1078             Double_t itsChi2[7] = {track->Pt(),track->Eta(), track->Phi(),
1079             fCentralityF,track->GetTPCsignalN(), sharebit,
1080             track->GetITSchi2()/static_cast<Double_t>(track->GetNcls(0))};
1081             fQACollection->Fill("fChi2perITScluster", itsChi2);
1082     }
1083     else{
1084             Double_t itsChi2[3] = {track->Pt(), fCentralityF, track->GetITSchi2()/static_cast<Double_t>(track->GetNcls(0))};
1085             fQACollection->Fill("fChi2perITScluster", itsChi2);
1086     }
1087
1088     // Fill Histogram for Hadronic Background
1089     if(HasMCData()){
1090       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11))
1091         fVarManager->FillContainer(fContainer, "hadronicBackground", UInt_t(0), kFALSE);
1092       else if(mctrack){
1093         // Fill Ke3 contributions
1094         Int_t glabel=TMath::Abs(mctrack->GetMother());
1095         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1096           if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==321)
1097             fQACollection->Fill("Ke3Kecorr",mctrackmother->Pt(),mctrack->Pt());
1098           else if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==130)
1099             fQACollection->Fill("Ke3K0Lecorr",mctrackmother->Pt(),mctrack->Pt());
1100         }
1101       }
1102     }
1103
1104     // Fill Containers
1105     if(signal) {
1106       // Apply weight for background contamination
1107             if(fBackGroundFactorApply) {
1108               if(IsPbPb() && fCentralityF >= 0) fWeightBackGround =  fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
1109               else    fWeightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
1110
1111               if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
1112               else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
1113         // weightBackGround as special weight
1114         fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
1115       }
1116       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
1117     }
1118
1119     Bool_t bTagged=kFALSE;
1120     if(GetPlugin(kSecVtx)) {
1121       AliDebug(2, "Running Secondary Vertex Analysis");
1122       if(fSecVtx->Process(track) && signal) {
1123         fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
1124         fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
1125         bTagged=kTRUE;
1126       }
1127     }
1128
1129     if(HasMCData()){
1130       dataE[0] = track->Pt();
1131       dataE[1] = track->Eta();
1132       dataE[2] = track->Phi();
1133       dataE[3] = track->Charge();
1134       dataE[4] = -1.;
1135       dataE[5] = -1.;
1136
1137       // Track selected: distinguish between true and fake
1138       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
1139       if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
1140         Int_t type = 0;
1141         if(fSignalCuts->IsCharmElectron(track))
1142           type = 1;
1143         else if(fSignalCuts->IsBeautyElectron(track))
1144           type = 2;
1145         AliDebug(1, Form("Type: %d\n", type));
1146         if(type){
1147                 dataE[5] = type; // beauty[1] or charm[2]
1148                 dataE[4] = 2;  // signal electron
1149         }
1150         else{
1151                 dataE[4] = 1; // not a signal electron
1152                 dataE[5] = 0;
1153         }
1154       }
1155       else {
1156         // Fill THnSparse with the information for Fake Electrons
1157         dataE[4] = 0;
1158         dataE[5] = 0;
1159       }
1160       // fill the performance THnSparse, if the mc origin could be defined
1161       if(dataE[4] > -1){
1162         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
1163         fQACollection->Fill("PIDperformance", dataE);
1164       }
1165     }
1166
1167     // Electron background analysis
1168     if (GetPlugin(kIsElecBackGround)) {
1169
1170       AliDebug(2, "Running BackGround Analysis");
1171
1172       for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
1173         htrack = fESD->GetTrack(jtrack);
1174         if ( itrack == jtrack ) continue;
1175         fElecBackGround->PairAnalysis(track, htrack);
1176       }
1177     } // end of electron background analysis
1178
1179
1180     // high dca track study [for temporary] ---------
1181     if (fTreeStream && fDebugLevel >= 1){
1182       Float_t b[2] = {0.,0.};
1183       Float_t bCov[3] = {0.,0.,0.};
1184       track->GetImpactParameters(b,bCov);
1185       Double_t dataD[11];
1186       dataD[0] = b[0]; // impact parameter xy
1187       dataD[1] = b[1]; // impact parameter z
1188       dataD[2] = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); // impact parameter space
1189       dataD[3]=0; dataD[4]=0;
1190       if(bCov[0]>0) dataD[3] = b[0]/TMath::Sqrt(bCov[0]); // normalised impact parameter xy
1191       if(bCov[2]>0) dataD[4] = b[1]/TMath::Sqrt(bCov[2]); // normalised impact parameter z
1192       dataD[5] = AliESDtrackCuts::GetSigmaToVertex(track); // n_sigma
1193       dataD[6] = track->GetTPCclusters(0x0);
1194       dataD[7] = track->GetITSclusters(0x0);
1195       dataD[8] = track->Eta();
1196       dataD[9] = track->Pt();
1197       Double_t p = track->GetP();
1198       Double_t phi = track->Phi();
1199       if(HasMCData()){
1200         if(fSignalCuts->IsCharmElectron(track)) dataD[10] = 0;
1201         else if(fSignalCuts->IsBeautyElectron(track)) dataD[10] = 1;
1202         else if(fSignalCuts->IsGammaElectron(track)) dataD[10] = 2;
1203         else if(fSignalCuts->IsNonHFElectron(track)) dataD[10] = 3;
1204         else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)) dataD[10] = 4;
1205         else dataD[10] = 5;
1206       }
1207       Double_t vtx[3];
1208       fInputEvent->GetPrimaryVertex()->GetXYZ(vtx);
1209       Double_t nt = fInputEvent->GetPrimaryVertex()->GetNContributors();
1210       Double_t nSigmaTPC = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
1211       Double_t runn = (Double_t)fInputEvent->GetRunNumber();
1212
1213       //printf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",dataD[0],dataD[1],dataD[2],dataD[3],dataD[4],dataD[5],dataD[6],dataD[7],dataD[8],dataD[9]);
1214       (*fTreeStream)<<"dcaQA"<<
1215         "dcaR="<<dataD[0]<<
1216         "dcaZ="<<dataD[1]<<
1217         "dca="<<dataD[2]<<
1218         "dcaSR="<<dataD[3]<<
1219         "dcaSZ="<<dataD[4]<<
1220         "dcaS="<<dataD[5]<<
1221         "nTPCclus="<<dataD[6]<<
1222         "nITSclus="<<dataD[7]<<
1223         "eta="<<dataD[8]<<
1224         "pt="<<dataD[9]<<
1225               "p="<< p <<
1226               "phi="<< phi <<
1227         "source="<<dataD[10] <<
1228         "vx=" << vtx[0] <<
1229         "vy=" << vtx[1] <<
1230         "vz=" << vtx[2] <<
1231               "nt=" << nt <<
1232         "TPCnSigma=" << nSigmaTPC <<
1233               "run=" << runn
1234         << "\n";
1235     }
1236     //-------------------------------------------
1237
1238     if (GetPlugin(kDEstep)) {
1239       Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.,};
1240       Int_t elecSource = 0;
1241       // minjung for IP QA(temporary ~ 2weeks)
1242       Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
1243       fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
1244       fQACollection->Fill("dataDca",track->Pt(),hfeimpactR);
1245       fQACollection->Fill("dataDcaSig",track->Pt(),hfeimpactnsigmaR);
1246       fQACollection->Fill("dataDcaSigDca",hfeimpactR,hfeimpactnsigmaR);
1247       if(HasMCData()){
1248         // minjung for IP QA(temporary ~ 2weeks)
1249         if(fSignalCuts->IsCharmElectron(track)){
1250           fQACollection->Fill("charmDca",track->Pt(),hfeimpactR);
1251           fQACollection->Fill("charmDcaSig",track->Pt(),hfeimpactnsigmaR);
1252         }
1253         else if(fSignalCuts->IsBeautyElectron(track)){
1254           fQACollection->Fill("beautyDca",track->Pt(),hfeimpactR);
1255           fQACollection->Fill("beautyDcaSig",track->Pt(),hfeimpactnsigmaR);
1256         }
1257         else if(fSignalCuts->IsGammaElectron(track)){
1258           fQACollection->Fill("conversionDca",track->Pt(),hfeimpactR);
1259           fQACollection->Fill("conversionDcaSig",track->Pt(),hfeimpactnsigmaR);
1260           fQACollection->Fill("conversionDcaSigDca",hfeimpactR,hfeimpactnsigmaR);
1261         }
1262         else if(fSignalCuts->IsNonHFElectron(track)){
1263           fQACollection->Fill("nonhfeDca",track->Pt(),hfeimpactR);
1264           fQACollection->Fill("nonhfeDcaSig",track->Pt(),hfeimpactnsigmaR);
1265         }
1266
1267         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1268           fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
1269           fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
1270         }
1271         if(fMCQA && !fFillSignalOnly) {
1272
1273           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1274             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1275             if(!fisNonHFEsystematics)break;
1276           }
1277
1278           if(fisNonHFEsystematics){
1279             //Fill additional containers for electron source distinction
1280             elecSource = fMCQA->GetElecSource(mctrack->Particle());
1281             const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1282             const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1283             Int_t iName = 0;
1284             for(Int_t iSource = AliHFEmcQA::kPi0; iSource <=  AliHFEmcQA::kGammaRho0; iSource++){
1285               if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1286               if(elecSource == iSource){
1287                 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1288                   if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, weightElecBgV0[iLevel]);
1289                   else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, -1*weightElecBgV0[iLevel]);
1290                 }
1291                 break;
1292               }
1293               iName++;
1294               if(iName == kElecBgSpecies)iName = 0;
1295             }
1296           }
1297           //else{
1298             if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]);
1299             else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]);
1300             //}
1301         }
1302       }
1303       // Fill Containers for impact parameter analysis
1304       if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
1305       if(HasMCData()){
1306         if(fMCQA && !fFillSignalOnly) {
1307           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1308             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1309             if(!fisNonHFEsystematics)break;
1310           }
1311           if(fisNonHFEsystematics){
1312             //Fill additional containers for electron source distinction
1313             elecSource = fMCQA->GetElecSource(mctrack->Particle());
1314             const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1315             const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1316             Int_t iName = 0;
1317             for(Int_t iSource = AliHFEmcQA::kPi0; iSource <=  AliHFEmcQA::kGammaRho0; iSource++){
1318               if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1319               if(elecSource == iSource){
1320                 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1321                   if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, weightElecBgV0[iLevel]);
1322                   else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, -1*weightElecBgV0[iLevel]);
1323                 }
1324                 break;
1325               }
1326               iName++;
1327               if(iName == kElecBgSpecies)iName = 0;
1328             }
1329           }
1330           // else{
1331             if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]);
1332             else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]);
1333             //}
1334         }
1335       }
1336       if(signal) {
1337         fVarManager->FillContainer(fContainer, "recTrackContDEReco", AliHFEcuts::kStepHFEcutsDca, kFALSE);
1338         fVarManager->FillContainer(fContainer, "recTrackContDEMC", AliHFEcuts::kStepHFEcutsDca, kTRUE);
1339         fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterDE"));
1340       }
1341       if(HasMCData()){
1342         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1343           fQACollection->Fill("hadronsAfterIPcut",track->Pt());
1344           fQACollection->Fill("hadronsAfterIPcutMC",mctrack->Pt());
1345         }
1346       }
1347     }
1348
1349   }
1350   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1351 }
1352
1353 //____________________________________________________________
1354 void AliAnalysisTaskHFE::ProcessAOD(){
1355   //
1356   // Run Analysis in AOD Mode
1357   // Function is still in development
1358   //
1359   AliDebug(3, "Processing AOD Event");
1360   Double_t eventContainer[4];
1361   if(HasMCData()) eventContainer[0] = fVz;
1362   else {
1363     eventContainer[0] = fInputEvent->GetPrimaryVertex()->GetZ();
1364   }
1365   eventContainer[1] = 1.; // No Information available in AOD analysis, assume all events have V0AND
1366   eventContainer[2] = fCentralityF;
1367   eventContainer[3] = fContributors;
1368
1369   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1370   if(!fAOD){
1371     AliError("AOD Event required for AOD Analysis");
1372       return;
1373   }
1374
1375   //
1376   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
1377
1378   //
1379   if(fIdentifiedAsPileUp) return;
1380   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
1381
1382   //
1383   if(fIdentifiedAsOutInz) return;
1384   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
1385
1386   //
1387   if(!fPassTheEventCut) return;
1388   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
1389
1390   fContainer->NewEvent();
1391
1392   AliAODTrack *track = NULL;
1393   AliAODMCParticle *mctrack = NULL;
1394   Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
1395   Int_t nElectronCandidates = 0;
1396   Int_t pid;
1397   Bool_t signal;
1398   for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
1399     track = fAOD->GetTrack(itrack); mctrack = NULL;
1400     if(!track) continue;
1401     if(track->GetFlags() != 1<<4) continue;  // Only process AOD tracks where the HFE is set
1402
1403     signal = kTRUE;
1404     if(HasMCData()){
1405
1406       Int_t label = TMath::Abs(track->GetLabel());
1407       if(label)
1408         mctrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(label));
1409         if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
1410     }
1411     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, kTRUE);
1412     // track accepted, do PID
1413     AliHFEpidObject hfetrack;
1414     hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1415     hfetrack.SetRecTrack(track);
1416     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1417     hfetrack.SetCentrality(fCentralityF);
1418     fPID->SetVarManager(fVarManager);
1419     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;    // we will do PID here as soon as possible
1420     // Apply weight for background contamination
1421     Double_t weightBackGround = 1.0;
1422
1423     // not correct treatment for pp
1424     if(fBackGroundFactorApply==kTRUE) {
1425             if(IsPbPb()) weightBackGround =  fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
1426       else weightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P()));
1427
1428       if(weightBackGround < 0.0) weightBackGround = 0.0;
1429             else if(weightBackGround > 1.0) weightBackGround = 1.0;
1430             fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, weightBackGround);
1431     }
1432
1433     nElectronCandidates++;
1434     if(HasMCData()){
1435       dataE[0] = track->Pt();
1436       dataE[1] = track->Eta();
1437       dataE[2] = track->Phi();
1438       dataE[3] = track->Charge();
1439       dataE[4] = -1;
1440       dataE[5] = -1;
1441       // Track selected: distinguish between true and fake
1442       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack ? mctrack->GetPdgCode(): -1));
1443       if(mctrack && ((pid = TMath::Abs(mctrack->GetPdgCode())) == 11)){
1444
1445         Int_t type = 0;
1446         if(fSignalCuts->IsCharmElectron(track))
1447           type = 1;
1448         else if(fSignalCuts->IsBeautyElectron(track))
1449           type = 2;
1450         AliDebug(1, Form("Type: %d\n", type));
1451         if(type){
1452                 dataE[5] = type; // beauty[1] or charm[2]
1453                 dataE[4] = 2;  // signal electron
1454         }
1455         else{
1456                 dataE[4] = 1; // not a signal electron
1457                 dataE[5] = 0;
1458         }
1459       }
1460       else {
1461         // Fill THnSparse with the information for Fake Electrons
1462         dataE[4] = 0;
1463         dataE[5] = 0;
1464       }
1465       // fill the performance THnSparse, if the mc origin could be defined
1466       if(dataE[4] > -1){
1467         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
1468         fQACollection->Fill("PIDperformance", dataE);
1469       }
1470     }
1471   }
1472   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1473 }
1474
1475 //____________________________________________________________
1476 Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
1477   //
1478   // Filter the Monte Carlo Track
1479   // Additionally Fill a THnSparse for Signal To Background Studies
1480   // Works for AOD and MC analysis Type
1481   //
1482   fVarManager->NewTrack(track, NULL, fCentralityF, -1, kTRUE);
1483   Double_t signalContainer[6];
1484
1485   signalContainer[0] = track->Pt();
1486   signalContainer[1] = track->Eta();
1487   signalContainer[2] = track->Phi();
1488   signalContainer[3] = track->Charge()/3;
1489
1490  Double_t vertex[3]; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
1491   if(IsESDanalysis()){
1492     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
1493     if(mctrack){
1494       vertex[0] = mctrack->Particle()->Vx();
1495       vertex[1] = mctrack->Particle()->Vy();
1496     }
1497   } else {
1498     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
1499     if(aodmctrack) aodmctrack->XvYvZv(vertex);
1500   }
1501
1502   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
1503   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
1504   signalContainer[4] = 0;
1505   if(fSignalCuts->IsSelected(track)){
1506     //fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCsignal, kFALSE);
1507     // Filling of the Signal/Background histogram using the
1508     // definition of the codes for charm and beauty as below in
1509     // th crearion of the histogram
1510     if(fSignalCuts->IsCharmElectron(track))
1511       signalContainer[4] = 1;
1512     else
1513       signalContainer[4] = 2;
1514   } else {
1515     signalContainer[4] = 0; // (and other background)
1516   }
1517   signalContainer[5] = 0;
1518   // apply cut on the sqrt of the production vertex
1519   Double_t radVertex = TMath::Sqrt(vertex[0]*vertex[0] + vertex[1] * vertex[1]);
1520   if(radVertex < 3.5){
1521     // Within first ITS layer(2) -> Background we cannot reject by ITS cut, let it pass
1522     signalContainer[5] = 1;
1523   } else if (radVertex < 7.5){
1524     signalContainer[5] = 2;
1525   }
1526   fQACollection->Fill("SignalToBackgroundMC", signalContainer);
1527
1528   // Step GeneratedZOutNoPileUp
1529   if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0)) return kFALSE;
1530   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
1531   //printf("In ProcessMCtrack %f\n",fCentralityF);
1532
1533   // Step Generated Event Cut
1534   if(!fPassTheEventCut) return kFALSE;
1535   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedEventCut, kFALSE);
1536
1537   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
1538   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCInAcceptance, kFALSE);
1539   return kTRUE;
1540 }
1541
1542 //____________________________________________________________
1543 Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
1544   //
1545   // Preselect tracks
1546   //
1547
1548
1549   Bool_t survived = kTRUE;
1550
1551   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) {
1552     survived = kFALSE;
1553     //printf("Did not pass AliHFEcuts::kStepRecKineITSTPC\n");
1554   }
1555   //else printf("Pass AliHFEcuts::kStepRecKineITSTPC\n");
1556   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) {
1557     survived = kFALSE;
1558     //printf("Did not pass AliHFEcuts::kStepRecPrim\n");
1559   }
1560   //else printf("Pass AliHFEcuts::kStepRecPrim\n");
1561   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) {
1562     survived = kFALSE;
1563     //printf("Did not pass AliHFEcuts::kStepHFEcutsITS\n");
1564   }
1565   //else printf("Pass AliHFEcuts::kStepHFEcutsITS\n");
1566   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTOF, track)) {
1567     survived = kFALSE;
1568     //printf("Did not pass AliHFEcuts::kStepHFEcutsTOF\n");
1569   }
1570   //else printf("Pass AliHFEcuts::kStepHFEcutsTOF\n");
1571   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) {
1572     survived = kFALSE;
1573     //printf("Did not pass AliHFEcuts::kStepHFEcutsTRD\n");
1574   }
1575   //else printf("Pass AliHFEcuts::kStepHFEcutsTRD\n");
1576
1577   if(survived){
1578     // Apply PID
1579     AliHFEpidObject hfetrack;
1580     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1581     hfetrack.SetRecTrack(track);
1582     if(!fPIDpreselect->IsSelected(&hfetrack)) {
1583       //printf("Did not pass AliHFEcuts::kPID\n");
1584       survived = kFALSE;
1585     }
1586     //else printf("Pass AliHFEcuts::kPID\n");
1587   }
1588
1589   return survived;
1590
1591 }
1592 //____________________________________________________________
1593 void AliAnalysisTaskHFE::MakeEventContainer(){
1594   //
1595   // Create the event container for the correction framework and link it
1596   // 1st bin: Vertex z-position
1597   // 2nd bin: V0AND decision (normalization to sigma_inel)
1598   // 3rd bin: Centrality class (for pp defined as number of contributors in vertex.)
1599   // 4th bin: Number of contributors > 0
1600   //
1601
1602   const Int_t kNvar = 4;  // number of variables on the grid:
1603   Int_t nBins[kNvar] = {120, 2, 11, 2};
1604   Double_t binMin[kNvar] = {-30. , 0., 0.0, 0.};
1605   Double_t binMax[kNvar] = {30., 2., 11.0, 2.};
1606
1607   AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
1608
1609   Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
1610   Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
1611   Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
1612   Double_t *contributorsBins = AliHFEtools::MakeLinearBinning(nBins[3], binMin[3], binMax[3]);
1613   evCont->SetBinLimits(0, vertexBins);
1614   evCont->SetBinLimits(1, v0andBins);
1615   evCont->SetBinLimits(2, centralityBins);
1616   evCont->SetBinLimits(3, contributorsBins);
1617   delete[] vertexBins; delete[] v0andBins; delete[] centralityBins; delete[] contributorsBins;
1618
1619   fCFM->SetEventContainer(evCont);
1620 }
1621
1622 //____________________________________________________________
1623 void AliAnalysisTaskHFE::MakeParticleContainer(){
1624   //
1625   // Create the particle container for the correction framework manager and
1626   // link it
1627   //
1628   if(!fContainer) fContainer = new AliHFEcontainer("trackContainer");
1629   fVarManager->DefineVariables(fContainer);
1630
1631   // Create Correction Framework containers
1632   fContainer->CreateContainer("MCTrackCont", "Track Container filled with MC information", AliHFEcuts::kNcutStepsMCTrack);
1633   fContainer->CreateContainer("recTrackContReco", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1634   fContainer->CreateContainer("recTrackContMC", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1635
1636   fContainer->CreateContainer("hadronicBackground", "Container for Hadronic Background", 2);
1637   fContainer->CreateContainer("recTrackContDEReco", "Container for displaced electron analysis with Reco information", 1);
1638   fContainer->CreateContainer("recTrackContDEMC", "Container for displaced electron analysis with MC information", 1);
1639   fContainer->CreateContainer("recTrackContSecvtxReco", "Container for secondary vertexing analysis with Reco information", 1);
1640   fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1);
1641
1642   if(HasMCData()){
1643     fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",4);
1644     fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",4);
1645
1646     if(fisNonHFEsystematics){
1647       const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1648       const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1649       for(Int_t iSource = 0; iSource < kElecBgSpecies; iSource++){
1650         for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1651           fContainer->CreateContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted conversion electrons from %s grandm., %s level",sourceName[iSource],levelName[iLevel]),4);
1652           fContainer->CreateContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted electrons from %s decays, %s level",sourceName[iSource],levelName[iLevel]),4);
1653         }
1654       }
1655     }
1656     //fContainer->CreateContainer("charmElecs", "Container for weighted charm electrons",2);
1657   }
1658
1659   fContainer->CreateCorrelationMatrix("correlationstepafterPID","THnSparse with correlations");
1660   fContainer->CreateCorrelationMatrix("correlationstepafterDE","THnSparse with correlations");
1661   if(!fVarManager->IsVariableDefined("centrality")) {
1662     //printf("Create the two other correlation maps\n");
1663     fContainer->CreateCorrelationMatrix("correlationstepbeforePID","THnSparse with correlations");
1664     fContainer->CreateCorrelationMatrix("correlationstepafterTOF","THnSparse with correlations");
1665   }
1666
1667   // Define the step names
1668   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsMCTrack; istep++){
1669     fContainer->SetStepTitle("MCTrackCont", AliHFEcuts::MCCutName(istep), istep);
1670   }
1671   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsRecTrack; istep++){
1672     fContainer->SetStepTitle("recTrackContReco", AliHFEcuts::RecoCutName(istep), istep);
1673     fContainer->SetStepTitle("recTrackContMC", AliHFEcuts::RecoCutName(istep), istep);
1674   }
1675   for(UInt_t ipid = 0; ipid < fPID->GetNumberOfPIDdetectors(); ipid++){
1676     fContainer->SetStepTitle("recTrackContReco", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1677     fContainer->SetStepTitle("recTrackContMC", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1678   }
1679 }
1680
1681 //____________________________________________________________
1682 void AliAnalysisTaskHFE::InitPIDperformanceQA(){
1683   // Add a histogram for Fake electrons
1684   const Int_t nDim=6;
1685   Int_t nBin[nDim] = {40, 8, 18, 2, 3, 3};
1686   //number of variables on the grid:pt,eta,phi,charge,
1687   const Double_t kPtbound[2] = {0.1, 20.};
1688   const Double_t kEtabound[2] = {-0.8, 0.8};
1689   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
1690   const Double_t kChargebound[2] = {-1.1, 1.1};
1691   const Double_t kAddInf1bound[2] = {0., 3.};
1692   const Double_t kAddInf2bound[2] = {0., 3.};
1693   Double_t minima[nDim] = {kPtbound[0], kEtabound[0], kPhibound[0], kChargebound[0], kAddInf1bound[0], kAddInf2bound[0]};
1694   Double_t maxima[nDim] = {kPtbound[1], kEtabound[1], kPhibound[1], kChargebound[1], kAddInf1bound[1], kAddInf2bound[1]};
1695
1696   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);
1697   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);
1698
1699   fQACollection->BinLogAxis("PIDperformance", 0);
1700   fQACollection->BinLogAxis("SignalToBackgroundMC", 0);
1701   fQACollection->Sumw2("PIDperformance");
1702   fQACollection->Sumw2("SignalToBackgroundMC");
1703 }
1704
1705 //____________________________________________________________
1706 void AliAnalysisTaskHFE::InitContaminationQA(){
1707   //
1708   // Add QA for Impact Parameter cut
1709   //
1710   const Double_t kPtbound[2] = {0.1, 20.};
1711   Int_t iBin[1];
1712   iBin[0] = 44; // bins in pt
1713   fQACollection->CreateTH1F("hadronsBeforeIPcut", "Hadrons before IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
1714   fQACollection->CreateTH1F("hadronsAfterIPcut", "Hadrons after IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
1715   fQACollection->CreateTH1F("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", iBin[0], kPtbound[0], kPtbound[1], 1);
1716   fQACollection->CreateTH1F("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
1717
1718   fQACollection->CreateTH2F("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1);
1719   fQACollection->CreateTH2F("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1);
1720   fQACollection->CreateTH1F("Kptspectra", "Charged Kaons: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
1721   fQACollection->CreateTH1F("K0Lptspectra", "K0L: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
1722
1723   const Double_t kDCAbound[2] = {-5., 5.};
1724   const Double_t kDCAsigbound[2] = {-50., 50.};
1725
1726   fQACollection->CreateTH2F("dataDcaSig", "data dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
1727   fQACollection->CreateTH2F("charmDcaSig", "charm dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
1728   fQACollection->CreateTH2F("beautyDcaSig", "beauty dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
1729   fQACollection->CreateTH2F("conversionDcaSig", "conversion dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
1730   fQACollection->CreateTH2F("nonhfeDcaSig", "nonhfe dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
1731
1732   fQACollection->CreateTH2F("dataDca", "data dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
1733   fQACollection->CreateTH2F("charmDca", "charm dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
1734   fQACollection->CreateTH2F("beautyDca", "beauty dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
1735   fQACollection->CreateTH2F("conversionDca", "conversion dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
1736   fQACollection->CreateTH2F("nonhfeDca", "nonhfe dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
1737
1738   fQACollection->CreateTH2F("dataDcaSigDca", "data dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0);
1739   fQACollection->CreateTH2F("conversionDcaSigDca", "conversion dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0);
1740 }
1741
1742 //____________________________________________________________
1743 void AliAnalysisTaskHFE::InitHistoITScluster(){
1744   //
1745     // Initialize a temporary histogram to monitor the chi2/ITS cluster
1746     if(IsPbPb()) {
1747         const Int_t kNDim = 7;
1748         const Int_t kNBins[kNDim] = {88, 20,90,11, 160, 2, 1000};
1749         const Double_t kMin[kNDim] = {0.1, -1,0,  0.,0., 0,  0.};
1750         const Double_t kMax[kNDim] = {20., 1, 2.*TMath::Pi(), 11.,160, 2, 100.};
1751         fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c);eta;phi; centrality class;nclus;sharebit; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
1752         fQACollection->BinLogAxis("fChi2perITScluster", 0);
1753     }
1754     else
1755     {
1756         const Int_t kNDim = 3;
1757         const Int_t kNBins[kNDim] = {44, 11, 1000};
1758         const Double_t kMin[kNDim] = {0.1, 0., 0.};
1759         const Double_t kMax[kNDim] = {20., 11., 100.};
1760         fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c); centrality class; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
1761         fQACollection->BinLogAxis("fChi2perITScluster", 0);
1762     }
1763 }
1764
1765 //____________________________________________________________
1766 void AliAnalysisTaskHFE::SelectSpecialTrigger(const Char_t *trgclust, Int_t runMin, Int_t runMax){
1767   //
1768   // Select only events triggered by a special trigeer cluster
1769   //
1770   if(!fSpecialTrigger) fSpecialTrigger = new AliOADBContainer("SpecialTrigger");
1771   fSpecialTrigger->AppendObject(new TObjString(trgclust), runMin, runMax);
1772 }
1773
1774 //____________________________________________________________
1775 const Char_t * AliAnalysisTaskHFE::GetSpecialTrigger(Int_t run){
1776   //
1777   // Derive selected trigger string for given run
1778   //
1779   if(!fSpecialTrigger) return NULL;
1780   TObjString *trg = dynamic_cast<TObjString *>(fSpecialTrigger->GetObject(run));
1781   if(!trg) return NULL;
1782   return trg->String().Data();
1783 }
1784
1785 //____________________________________________________________
1786 void AliAnalysisTaskHFE::PrintStatus() const {
1787   //
1788   // Print Analysis status
1789   //
1790   printf("\n\tAnalysis Settings\n\t========================================\n\n");
1791   printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
1792   printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
1793   printf("\tDisplaced electron analysis step: %s\n", GetPlugin(kDEstep) ? "YES" : "NO");
1794   printf("\tTagged Track Analysis: %s\n", GetPlugin(kTaggedTrackAnalysis) ? "YES" : "NO");
1795   printf("\n");
1796   printf("\tParticle Identification Detectors:\n");
1797   fPID->PrintStatus();
1798   printf("\n");
1799   printf("\tQA: \n");
1800   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
1801   printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsQAOn()) ? "YES" : "NO");
1802   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
1803   printf("\n");
1804 }
1805
1806 //____________________________________________________________
1807 Bool_t AliAnalysisTaskHFE::FillProductionVertex(const AliVParticle * const track) const{
1808   //
1809   // Find the production vertex of the associated MC track
1810   //
1811   if(!fMCEvent) return kFALSE;
1812   const AliVParticle *mctrack = NULL;
1813   TString objectType = track->IsA()->GetName();
1814   if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
1815     // Reconstructed track
1816     mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
1817   } else {
1818     // MCParticle
1819     mctrack = track;
1820   }
1821
1822   if(!mctrack) return kFALSE;
1823
1824   Double_t xv = 0.0;
1825   Double_t yv = 0.0;
1826
1827   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
1828     // case MCParticle
1829     const AliMCParticle *mcpart = dynamic_cast<const AliMCParticle *>(mctrack);
1830     if(mcpart){
1831       xv =  mcpart->Xv();
1832       yv =  mcpart->Yv();
1833     }
1834   } else {
1835     // case AODMCParticle
1836     const AliAODMCParticle *mcpart = dynamic_cast<const AliAODMCParticle *>(mctrack);
1837     if(mcpart){
1838       xv =  mcpart->Xv();
1839       yv =  mcpart->Yv();
1840     }
1841   }
1842
1843   //printf("xv %f, yv %f\n",xv,yv);
1844   fQACollection->Fill("radius", TMath::Abs(xv),TMath::Abs(yv));
1845
1846   return kTRUE;
1847
1848 }
1849 //__________________________________________
1850 void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
1851   //
1852   // Switch on Plugin
1853   // Available:
1854   //  - Primary vertex studies
1855   //  - Secondary vertex Studies
1856   //  - Post Processing
1857   //
1858   switch(plug){
1859     case kPriVtx: SETBIT(fPlugins, plug); break;
1860     case kSecVtx: SETBIT(fPlugins, plug); break;
1861     case kIsElecBackGround: SETBIT(fPlugins, plug); break;
1862     case kPostProcess: SETBIT(fPlugins, plug); break;
1863     case kDEstep: SETBIT(fPlugins, plug); break;
1864     case kTaggedTrackAnalysis: SETBIT(fPlugins, plug); break;
1865     default: AliError("Unknown Plugin");
1866   };
1867 }
1868 //__________________________________________
1869 Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track){
1870   //
1871   // Check single track cuts for a given cut step
1872   // Fill the particle container
1873   //
1874   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1875   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1876   if(fVarManager->IsSignalTrack()) {
1877     fVarManager->FillContainer(fContainer, "recTrackContReco", cutStep, kFALSE);
1878     fVarManager->FillContainer(fContainer, "recTrackContMC", cutStep, kTRUE);
1879   }
1880   return kTRUE;
1881 }
1882 //___________________________________________________
1883 Bool_t AliAnalysisTaskHFE::ReadCentrality() {
1884   //
1885   // Recover the centrality of the event from ESD or AOD
1886   //
1887   Int_t bin = -1;
1888   if(IsPbPb()) {
1889     // Centrality
1890     AliCentrality *centrality = fInputEvent->GetCentrality();
1891     Float_t fCentralityFtemp = centrality->GetCentralityPercentile("V0M");
1892     Float_t centralityLimits[12] = {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
1893     for(Int_t ibin = 0; ibin < 11; ibin++){
1894       if(fCentralityFtemp >= centralityLimits[ibin] && fCentralityFtemp < centralityLimits[ibin+1]){
1895           bin = ibin;
1896         break;
1897       }
1898     }
1899     if(bin == -1) bin = 11; // Overflow
1900   } else {
1901     // PP: Tracklet multiplicity, use common definition
1902     Int_t itsMultiplicity = GetITSMultiplicity(fInputEvent);
1903     Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
1904     for(Int_t ibin = 0; ibin < 7; ibin++){
1905       if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){
1906         bin = ibin;
1907         break;
1908       }
1909     }
1910     if(bin == -1) bin = 7;  // Overflow
1911   }
1912   fCentralityF = bin;
1913   AliDebug(2, Form("Centrality class %d\n", fCentralityF));
1914
1915   // contributors, to be outsourced
1916   const AliVVertex *vtx;
1917   if(IsAODanalysis()){
1918     AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1919     if(!fAOD){
1920       AliError("AOD Event required for AOD Analysis");
1921       return kFALSE;
1922     }
1923     vtx = fAOD->GetPrimaryVertex();
1924   } else {
1925     AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
1926     if(!fESD){
1927       AliError("ESD Event required for ESD Analysis");
1928       return kFALSE;
1929     }
1930     vtx = fESD->GetPrimaryVertexSPD();
1931   }
1932   if(!vtx){
1933     fContributors = 0.5;
1934     return kFALSE;
1935   }
1936   else {
1937     Int_t contributorstemp = vtx->GetNContributors();
1938     if( contributorstemp <=  0) fContributors =  0.5;
1939     else fContributors = 1.5;
1940   }
1941   return kTRUE;
1942 }
1943
1944 //___________________________________________________
1945 Int_t AliAnalysisTaskHFE::GetITSMultiplicity(AliVEvent *ev){
1946   //
1947   // Definition of the Multiplicity according to the JPSI group (F. Kramer)
1948   //
1949   Int_t nTracklets = 0;
1950   Int_t nAcc = 0;
1951   Double_t etaRange = 1.6;
1952
1953   if (ev->IsA() == AliAODEvent::Class()) {
1954     AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
1955     nTracklets = tracklets->GetNumberOfTracklets();
1956     for (Int_t nn = 0; nn < nTracklets; nn++) {
1957       Double_t theta = tracklets->GetTheta(nn);
1958       Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
1959       if (TMath::Abs(eta) < etaRange) nAcc++;
1960     }
1961   } else if (ev->IsA() == AliESDEvent::Class()) {
1962     nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
1963     for (Int_t nn = 0; nn < nTracklets; nn++) {
1964        Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
1965       if (TMath::Abs(eta) < etaRange) nAcc++;
1966     }
1967   } else return -1;
1968
1969   return nAcc;
1970 }
1971
1972 //___________________________________________________
1973 Bool_t AliAnalysisTaskHFE::InitializeHadronBackground(Int_t run){
1974   //
1975   // Initialize background factors array from the AliOADBContainer
1976   // The container is expected to provide a TObjArray with the name
1977   // "hadronBackground" and the size 12 for the given run number
1978   //
1979   AliDebug(1, "Deriving hadronic background parameterization from OADB container");
1980   TObjArray *functions = dynamic_cast<TObjArray *>(fHadronBackgroundOADB->GetObject(run, "hadronBackground"));
1981   if(!functions){
1982     AliDebug(1, "Content in the OADB Container is not a TObjArray");
1983     fBackGroundFactorApply = kFALSE;
1984     return kFALSE;
1985   }
1986   if(functions->GetSize() < 12){
1987     AliDebug(1, Form("Size not matching: 12 expected, %d provided", functions->GetSize()));
1988     fBackGroundFactorApply = kFALSE;
1989     return kFALSE;
1990   }
1991   for(Int_t icent = 0; icent < 12; icent++) fkBackGroundFactorArray[icent] = dynamic_cast<const TF1 *>(functions->UncheckedAt(icent));
1992   return kTRUE;
1993 }
1994
1995 //___________________________________________________
1996 void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
1997   //
1998   // Recover the centrality of the event from ESD or AOD
1999   //
2000  if(IsAODanalysis()){
2001
2002    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
2003    if(!fAOD){
2004      AliError("AOD Event required for AOD Analysis");
2005        return;
2006    }
2007    // PileUp
2008    if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
2009    // Z vertex
2010    if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2011    // Event Cut
2012    fPassTheEventCut = kTRUE;
2013    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) fPassTheEventCut = kFALSE;
2014
2015
2016  } else {
2017
2018    AliDebug(3, "Processing ESD Centrality");
2019    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
2020    if(!fESD){
2021      AliError("ESD Event required for ESD Analysis");
2022        return;
2023    }
2024    // PileUp
2025    fIdentifiedAsPileUp = kFALSE;
2026    if(fRemovePileUp && fESD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
2027    // Z vertex
2028    fIdentifiedAsOutInz = kFALSE;
2029    if (IsPbPb()) {
2030      //printf("PbPb\n");
2031      if(fESD->GetPrimaryVertex()){
2032        if(TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2033      }
2034    }
2035    else {
2036      //printf("pp\n");
2037      if(fESD->GetPrimaryVertexTracks()){
2038        if(TMath::Abs(fESD->GetPrimaryVertexTracks()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2039      }
2040    }
2041    //Event Cut
2042    fPassTheEventCut = kTRUE;
2043    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) fPassTheEventCut = kFALSE;
2044
2045
2046  }
2047
2048 }
2049