PWGHF/hfe/AliAnalysisTaskEMCalHFEpA.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 ////////////////////////////////////////////////////////////////////////
  17 //                                                                    //
  18 //      Task for Heavy-flavour electron analysis in pPb collisions    //
  19 //      (+ Electron-Hadron Jetlike Azimuthal Correlation)             //
  20 //                                                                                                                                        //
  21 //              v1.0                                                                                                              //
  22 //                                                                    //
  23 //          Authors                                                                               //
  24 //              Elienos Pereira de Oliveira Filho (epereira@cern.ch)          //
  25 //              Cristiane Jahnke                (cristiane.jahnke@cern.ch)                    //
  26 //                                                                    //
  27 ////////////////////////////////////////////////////////////////////////
  28
  29 #include "TChain.h"
  30 #include "TTree.h"
  31 #include "TNtuple.h"
  32 #include "TH1F.h"
  33 #include "TH2F.h"
  34 #include "TCanvas.h"
  35 #include "AliAnalysisTask.h"
  36 #include "AliAnalysisManager.h"
  37 #include "AliESDEvent.h"
  38 #include "AliAODEvent.h"
  39 #include "AliVEvent.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliESDtrackCuts.h"
  42 #include "AliESDCaloCluster.h"
  43 #include "AliESDCaloCells.h"
  44 #include "AliEMCALTrack.h"
  45 #include "AliExternalTrackParam.h"
  46 #include "AliPhysicsSelection.h"
  47 #include "TGeoGlobalMagField.h"
  48 #include "AliMagF.h"
  49 #include "AliLog.h"
  50 #include "AliStack.h"
  51 #include "AliCentrality.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliAODMCHeader.h"
  54 #include "AliPID.h"
  55 #include "AliPIDResponse.h"
  56 #include "AliHFEcontainer.h"
  57 #include "AliHFEcuts.h"
  58 #include "AliHFEpid.h"
  59 #include "AliHFEpidBase.h"
  60 #include "AliHFEpidQAmanager.h"
  61 #include "AliHFEtools.h"
  62 #include "AliCFContainer.h"
  63 #include "AliCFManager.h"
  64 #include "AliSelectNonHFE.h"
  65 #include "AliHFEpidTPC.h"
  66 #include "AliAnalysisTaskEMCalHFEpA.h"
  67 #include "TMath.h"
  68 #include "THnSparse.h"
  69 #include "TLorentzVector.h"
  70 #include "TString.h"
  71 #include "TFile.h"
  72 #include "AliESDHandler.h"
  73 #include "AliMCEventHandler.h"
  74 #include "AliMCEvent.h"
  75 #include "AliStack.h"
  76 #include "TParticle.h"
  77 #include "AliLog.h"
  78 #include "AliAnalysisTaskSE.h"
  79 #include "TRefArray.h"
  80 #include "TVector.h"
  81 #include "stdio.h"
  82 #include "TGeoManager.h"
  83 #include "iostream"
  84 #include "fstream"
  85 #include "AliKFParticle.h"
  86 #include "AliKFVertex.h"
  87 #include "AliVParticle.h"
  88 #include "AliVTrack.h"
  89 #include "AliEventPoolManager.h"
  90 #include "TObjArray.h"
  91 //______________________________________________________________________
  92
  93 //______________________________________________________________________
  94 ClassImp(AliAnalysisTaskEMCalHFEpA)
  95
  96 //______________________________________________________________________
  97 AliAnalysisTaskEMCalHFEpA::AliAnalysisTaskEMCalHFEpA(const char *name)
  98   : AliAnalysisTaskSE(name)
  99         ,fCorrelationFlag(0)
 100         ,fIsMC(0)
 101         ,fUseEMCal(kFALSE)
 102         ,fEMCEG1(kFALSE)
 103         ,fEMCEG2(kFALSE)
 104         ,fIsHFE1(kFALSE)
 105         ,fIsHFE2(kFALSE)
 106         ,fIsNonHFE(kFALSE)
 107         ,fIsFromD(kFALSE)
 108         ,fIsFromB(kFALSE)
 109         ,fIsFromPi0(kFALSE)
 110         ,fIsFromEta(kFALSE)
 111         ,fIsFromGamma(kFALSE)
 112         ,fESD(0)
 113         ,fAOD(0)
 114         ,fVevent(0)
 115         ,fPartnerCuts(new AliESDtrackCuts())
 116         ,fOutputList(0)
 117         ,fPidResponse(0)
 118         ,fNonHFE(new AliSelectNonHFE())
 119         ,fIsAOD(kFALSE)
 120         ,fCentrality(0)
 121         ,fCentralityMin(0)
 122         ,fCentralityMax(100)
 123         ,fHasCentralitySelection(kFALSE)
 124         ,fCentralityHist(0)
 125         ,fCentralityHistPass(0)
 126         ,fZvtx(0)
 127         ,fEstimator(0)
 128         ,fClus(0)
 129         ,fNevent(0)
 130         ,fPtElec_Inc(0)
 131         ,fPtElec_ULS(0)
 132         ,fPtElec_LS(0)
 133         ,fpid(0)
 134         ,fEoverP_pt(0)
 135         ,fEoverP_tpc(0)
 136         ,fTPC_pt(0)
 137         ,fTPC_p(0)
 138         ,fTPCnsigma_pt(0)
 139         ,fTPCnsigma_p(0)
 140     ,fTPCnsigma_pt_2D(0)
 141         ,fTPCnsigma_eta(0)
 142         ,fTPCnsigma_phi(0)
 143         ,fECluster(0)
 144         ,fEtaPhi(0)
 145         ,fVtxZ(0)
 146         ,fNTracks(0)
 147         ,fNClusters(0)
 148         ,fTPCNcls_EoverP(0)
 149         ,fEta(0)
 150         ,fPhi(0)
 151         ,fR(0)
 152         ,fR_EoverP(0)
 153         ,fNcells(0)
 154         ,fNcells_EoverP(0)
 155         ,fNcells_electrons(0)
 156         ,fNcells_hadrons(0)
 157         ,fECluster_ptbins(0)
 158         ,fEoverP_ptbins(0)
 159         ,fEoverP_wSSCut(0)
 160         ,fM02_EoverP(0)
 161         ,fM20_EoverP(0)
 162         ,fTPCnsigma_eta_electrons(0)
 163         ,fTPCnsigma_eta_hadrons(0)
 164         ,fEoverP_pt_pions(0)
 165     ,ftpc_p_EoverPcut(0)
 166     ,fnsigma_p_EoverPcut(0)
 167         ,fEoverP_pt_pions2(0)
 168         ,fNcells_pt(0)
 169         ,fEoverP_pt_hadrons(0)
 170         ,fCEtaPhi_Inc(0)
 171         ,fCEtaPhi_ULS(0)
 172         ,fCEtaPhi_LS(0)
 173         ,fCEtaPhi_ULS_NoP(0)
 174         ,fCEtaPhi_LS_NoP(0)
 175         ,fCEtaPhi_ULS_Weight(0)
 176         ,fCEtaPhi_LS_Weight(0)
 177         ,fCEtaPhi_ULS_NoP_Weight(0)
 178         ,fCEtaPhi_LS_NoP_Weight(0)
 179         ,fInvMass(0)
 180         ,fInvMassBack(0)
 181         ,fDCA(0)
 182         ,fDCABack(0)
 183         ,fOpAngle(0)
 184         ,fOpAngleBack(0)
 185         ,fMassCut(0.1)
 186         ,fEtaCutMin(-0.9)
 187         ,fEtaCutMax(0.9)
 188         ,fEoverPCutMin(0.8)
 189         ,fEoverPCutMax(1.2)
 190         ,fAngleCut(999)
 191         ,fChi2Cut(3.5)
 192         ,fDCAcut(999)
 193         ,fMassCutFlag(kTRUE)
 194         ,fAngleCutFlag(kFALSE)
 195         ,fChi2CutFlag(kFALSE)
 196         ,fDCAcutFlag(kFALSE)
 197         ,fPtBackgroundBeforeReco(0)
 198         ,fPtBackgroundAfterReco(0)
 199         ,fPtMCparticleAll(0)
 200         ,fPtMCparticleReco(0)
 201         ,fPtMCparticleAllHfe1(0)
 202         ,fPtMCparticleRecoHfe1(0)
 203         ,fPtMCparticleAllHfe2(0)
 204         ,fPtMCparticleRecoHfe2(0)
 205         ,fPtMCelectronAfterAll(0)
 206         ,fPtMCpi0(0)
 207         ,fPtMC_EMCal_All(0)
 208         ,fPtMC_EMCal_Selected(0)
 209         ,fPtMC_TPC_All(0)
 210         ,fPtMC_TPC_Selected(0)
 211         ,fPtMCWithLabel(0)
 212         ,fPtMCWithoutLabel(0)
 213         ,fPtIsPhysicaPrimary(0)
 214         ,fCuts(0)
 215         ,fCFM(0)
 216         ,fPID(new AliHFEpid("hfePid"))
 217         ,fPIDqa(0)
 218         ,fMCstack(0)
 219         ,fRejectKinkMother(kFALSE)
 220         ,fMCtrack(0)
 221         ,fMCtrackMother(0)
 222         ,fMCtrackGMother(0)
 223         ,fMCtrackGGMother(0)
 224         ,fMCtrackGGGMother(0)
 225         ,fMCarray(0)
 226         ,fMCheader(0)
 227         ,fMCparticle(0)
 228         ,fMCparticleMother(0)
 229         ,fMCparticleGMother(0)
 230         ,fMCparticleGGMother(0)
 231         ,fMCparticleGGGMother(0)
 232         ,fEventHandler(0)
 233         ,fMCevent(0)
 234         ,fPoolMgr(0)
 235         ,fPool(0)
 236         ,fTracksClone(0)
 237         ,fTracks(0)
 238         ,fCEtaPhi_Inc_EM(0)
 239         ,fCEtaPhi_ULS_EM(0)
 240         ,fCEtaPhi_LS_EM(0)
 241         ,fCEtaPhi_ULS_Weight_EM(0)
 242         ,fCEtaPhi_LS_Weight_EM(0)
 243         ,fPoolNevents(0)
 244         ,fEventMixingFlag(0)
 245         ,fCEtaPhi_Inc_DiHadron(0)
 246         ,fPtTrigger_Inc(0)
 247 {
 248   //Named constructor
 249   // Define input and output slots here
 250   // Input slot #0 works with a TChain
 251   DefineInput(0, TChain::Class());
 252   // Output slot #0 id reserved by the base class for AOD
 253   // Output slot #1 writes into a TH1 container
 254   // DefineOutput(1, TH1I::Class());
 255   DefineOutput(1, TList::Class());
 256   //  DefineOutput(3, TTree::Class());
 257 }
 258
 259 //________________________________________________________________________
 260 AliAnalysisTaskEMCalHFEpA::AliAnalysisTaskEMCalHFEpA()
 261   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCalHFEpA")
 262         ,fCorrelationFlag(0)
 263         ,fIsMC(0)
 264         ,fUseEMCal(kFALSE)
 265         ,fEMCEG1(kFALSE)
 266         ,fEMCEG2(kFALSE)
 267         ,fIsHFE1(kFALSE)
 268         ,fIsHFE2(kFALSE)
 269         ,fIsNonHFE(kFALSE)
 270         ,fIsFromD(kFALSE)
 271         ,fIsFromB(kFALSE)
 272         ,fIsFromPi0(kFALSE)
 273         ,fIsFromEta(kFALSE)
 274         ,fIsFromGamma(kFALSE)
 275         ,fESD(0)
 276         ,fAOD(0)
 277         ,fVevent(0)
 278         ,fPartnerCuts(new AliESDtrackCuts())
 279         ,fOutputList(0)
 280         ,fPidResponse(0)
 281         ,fNonHFE(new AliSelectNonHFE())
 282         ,fIsAOD(kFALSE)
 283         ,fCentrality(0)
 284         ,fCentralityMin(0)
 285         ,fCentralityMax(100)
 286         ,fHasCentralitySelection(kFALSE)
 287         ,fCentralityHist(0)
 288         ,fCentralityHistPass(0)
 289         ,fZvtx(0)
 290         ,fEstimator(0)
 291         ,fClus(0)
 292         ,fNevent(0)
 293         ,fPtElec_Inc(0)
 294         ,fPtElec_ULS(0)
 295         ,fPtElec_LS(0)
 296         ,fpid(0)
 297         ,fEoverP_pt(0)
 298         ,fEoverP_tpc(0)
 299         ,fTPC_pt(0)
 300         ,fTPC_p(0)
 301         ,fTPCnsigma_pt(0)
 302         ,fTPCnsigma_p(0)
 303     ,fTPCnsigma_pt_2D(0)
 304         ,fTPCnsigma_eta(0)
 305         ,fTPCnsigma_phi(0)
 306         ,fECluster(0)
 307         ,fEtaPhi(0)
 308         ,fVtxZ(0)
 309         ,fNTracks(0)
 310         ,fNClusters(0)
 311         ,fTPCNcls_EoverP(0)
 312         ,fEta(0)
 313         ,fPhi(0)
 314         ,fR(0)
 315         ,fR_EoverP(0)
 316         ,fNcells(0)
 317         ,fNcells_EoverP(0)
 318         ,fNcells_electrons(0)
 319         ,fNcells_hadrons(0)
 320         ,fECluster_ptbins(0)
 321         ,fEoverP_ptbins(0)
 322         ,fEoverP_wSSCut(0)
 323         ,fM02_EoverP(0)
 324         ,fM20_EoverP(0)
 325         ,fTPCnsigma_eta_electrons(0)
 326         ,fTPCnsigma_eta_hadrons(0)
 327         ,fEoverP_pt_pions(0)
 328     ,ftpc_p_EoverPcut(0)
 329     ,fnsigma_p_EoverPcut(0)
 330         ,fEoverP_pt_pions2(0)
 331         ,fNcells_pt(0)
 332         ,fEoverP_pt_hadrons(0)
 333         ,fCEtaPhi_Inc(0)
 334         ,fCEtaPhi_ULS(0)
 335         ,fCEtaPhi_LS(0)
 336         ,fCEtaPhi_ULS_NoP(0)
 337         ,fCEtaPhi_LS_NoP(0)
 338         ,fCEtaPhi_ULS_Weight(0)
 339         ,fCEtaPhi_LS_Weight(0)
 340         ,fCEtaPhi_ULS_NoP_Weight(0)
 341         ,fCEtaPhi_LS_NoP_Weight(0)
 342         ,fInvMass(0)
 343         ,fInvMassBack(0)
 344         ,fDCA(0)
 345         ,fDCABack(0)
 346         ,fOpAngle(0)
 347         ,fOpAngleBack(0)
 348         ,fMassCut(0.1)
 349         ,fEtaCutMin(-0.9)
 350         ,fEtaCutMax(0.9)
 351         ,fEoverPCutMin(0.8)
 352         ,fEoverPCutMax(1.2)
 353         ,fAngleCut(999)
 354         ,fChi2Cut(3.5)
 355         ,fDCAcut(999)
 356         ,fMassCutFlag(kTRUE)
 357         ,fAngleCutFlag(kFALSE)
 358         ,fChi2CutFlag(kFALSE)
 359         ,fDCAcutFlag(kFALSE)
 360         ,fPtBackgroundBeforeReco(0)
 361         ,fPtBackgroundAfterReco(0)
 362         ,fPtMCparticleAll(0)
 363         ,fPtMCparticleReco(0)
 364         ,fPtMCparticleAllHfe1(0)
 365         ,fPtMCparticleRecoHfe1(0)
 366         ,fPtMCparticleAllHfe2(0)
 367         ,fPtMCparticleRecoHfe2(0)
 368         ,fPtMCelectronAfterAll(0)
 369         ,fPtMCpi0(0)
 370         ,fPtMC_EMCal_All(0)
 371         ,fPtMC_EMCal_Selected(0)
 372         ,fPtMC_TPC_All(0)
 373         ,fPtMC_TPC_Selected(0)
 374         ,fPtMCWithLabel(0)
 375         ,fPtMCWithoutLabel(0)
 376         ,fPtIsPhysicaPrimary(0)
 377         ,fCuts(0)
 378         ,fCFM(0)
 379         ,fPID(new AliHFEpid("hfePid"))
 380         ,fPIDqa(0)
 381         ,fMCstack(0)
 382         ,fRejectKinkMother(kFALSE)
 383         ,fMCtrack(0)
 384         ,fMCtrackMother(0)
 385         ,fMCtrackGMother(0)
 386         ,fMCtrackGGMother(0)
 387         ,fMCtrackGGGMother(0)
 388         ,fMCarray(0)
 389         ,fMCheader(0)
 390         ,fMCparticle(0)
 391         ,fMCparticleMother(0)
 392         ,fMCparticleGMother(0)
 393         ,fMCparticleGGMother(0)
 394         ,fMCparticleGGGMother(0)
 395         ,fEventHandler(0)
 396         ,fMCevent(0)
 397         ,fPoolMgr(0)
 398         ,fPool(0)
 399         ,fTracksClone(0)
 400         ,fTracks(0)
 401         ,fCEtaPhi_Inc_EM(0)
 402         ,fCEtaPhi_ULS_EM(0)
 403         ,fCEtaPhi_LS_EM(0)
 404         ,fCEtaPhi_ULS_Weight_EM(0)
 405         ,fCEtaPhi_LS_Weight_EM(0)
 406         ,fPoolNevents(0)
 407         ,fEventMixingFlag(0)
 408         ,fCEtaPhi_Inc_DiHadron(0)
 409         ,fPtTrigger_Inc(0)
 410 {
 411         // Constructor
 412         // Define input and output slots here
 413         // Input slot #0 works with a TChain
 414         DefineInput(0, TChain::Class());
 415         // Output slot #0 id reserved by the base class for AOD
 416         // Output slot #1 writes into a TH1 container
 417         // DefineOutput(1, TH1I::Class());
 418         DefineOutput(1, TList::Class());
 419         //DefineOutput(3, TTree::Class());
 420 }
 421
 422 //______________________________________________________________________
 423 AliAnalysisTaskEMCalHFEpA::~AliAnalysisTaskEMCalHFEpA()
 424 {
 425         //Destructor
 426         delete fOutputList;
 427         delete fPID;
 428         delete fCFM;
 429         delete fPIDqa;
 430 }
 431
 432 //______________________________________________________________________
 433 //Create Output Objects
 434 //Here we can define the histograms and others output files
 435 //Called once
 436 void AliAnalysisTaskEMCalHFEpA::UserCreateOutputObjects()
 437 {
 438 //______________________________________________________________________
 439 //Initialize PID
 440         if(!fPID->GetNumberOfPIDdetectors())
 441     {
 442                 fPID->AddDetector("TPC", 0);
 443     }
 444
 445         fPID->SortDetectors();
 446
 447         fPIDqa = new AliHFEpidQAmanager();
 448         fPIDqa->Initialize(fPID);
 449 //______________________________________________________________________
 450
 451 //______________________________________________________________________
 452 //Initialize correction Framework and Cuts
 453         fCFM = new AliCFManager;
 454         const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 455         fCFM->SetNStepParticle(kNcutSteps);
 456         for(Int_t istep = 0; istep < kNcutSteps; istep++) fCFM->SetParticleCutsList(istep, NULL);
 457
 458         if(!fCuts)
 459         {
 460                 AliWarning("Cuts not available. Default cuts will be used");
 461                 fCuts = new AliHFEcuts;
 462                 fCuts->CreateStandardCuts();
 463         }
 464
 465         fCuts->Initialize(fCFM);
 466 //______________________________________________________________________
 467
 468 ///______________________________________________________________________
 469 ///Output Tlist
 470 //Create TList
 471         fOutputList = new TList();
 472         fOutputList->SetOwner();
 473
 474 //PIDqa
 475         fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 476
 477 //Store the number of events
 478         //Define the histo
 479         fNevent = new TH1F("fNevent","Number of Events",5,-0.5,4.5);
 480         //And then, add to the output list
 481         fOutputList->Add(fNevent);
 482
 483         fpid = new TH1F("fpid","PID flag",5,0,5);
 484         fOutputList->Add(fpid);
 485
 486 //pt Distribution
 487         fPtElec_Inc = new TH1F("fPtElec_Inc","Inclusive Electrons; p_{t} (GeV/c); Count",200,0,20);
 488         fPtElec_ULS = new TH1F("fPtElec_ULS","Inclusive Electrons; p_{t} (GeV/c); Count",200,0,20);
 489         fPtElec_LS = new TH1F("fPtElec_LS","Inclusive Electrons; p_{t} (GeV/c); Count",200,0,20);
 490         fPtTrigger_Inc = new TH1F("fPtTrigger_Inc","pT dist for Hadron Contamination; p_{t} (GeV/c); Count",200,0,20);
 491         fTPCnsigma_pt_2D = new TH2F("fTPCnsigma_pt_2D",";pt (GeV/c);TPC Electron N#sigma",1000,0.3,30,1000,-15,10);
 492
 493
 494         fOutputList->Add(fPtElec_Inc);
 495         fOutputList->Add(fPtElec_ULS);
 496         fOutputList->Add(fPtElec_LS);
 497         fOutputList->Add(fPtTrigger_Inc);
 498         fOutputList->Add(fTPCnsigma_pt_2D);
 499
 500         //General Histograms
 501
 502         //Steps
 503         //Step 1: Before Track cuts
 504         //Step 2: Before PID
 505         //Step 3: After PID
 506
 507         fEoverP_pt = new TH2F *[3];
 508         fTPC_p = new TH2F *[3];
 509         fTPCnsigma_p = new TH2F *[3];
 510
 511         fECluster= new TH1F *[3];
 512         fEtaPhi= new TH2F *[3];
 513         fVtxZ= new  TH1F *[3];
 514         fNTracks= new  TH1F *[3];
 515         fNClusters= new TH1F *[3];
 516         fTPCNcls_EoverP= new TH2F *[3];
 517
 518         for(Int_t i = 0; i < 3; i++)
 519         {
 520           fEoverP_pt[i] = new TH2F(Form("fEoverP_pt%d",i),";p_{t} (GeV/c);E / p ",1000,0,30,500,0,2);
 521           fTPC_p[i] = new TH2F(Form("fTPC_p%d",i),";pt (GeV/c);TPC dE/dx (a. u.)",1000,0.3,15,1000,-20,200);
 522           fTPCnsigma_p[i] = new TH2F(Form("fTPCnsigma_p%d",i),";p (GeV/c);TPC Electron N#sigma",1000,0.3,15,1000,-15,10);
 523
 524
 525           fECluster[i]= new TH1F(Form("fECluster%d",i), ";ECluster",2000, 0,100);
 526           fEtaPhi[i]= new TH2F(Form("fEtaPhi%d",i),"#eta x #phi Clusters;#phi;#eta",200,0.,TMath::Pi(),50,-1.,1.);
 527       fVtxZ[i]= new  TH1F(Form("fVtxZ%d",i),"VtxZ",1000, -50,50);
 528           fNTracks[i]= new  TH1F(Form("fNTracks%d",i),"NTracks",1000, 0,1000);
 529       fNClusters[i]= new TH1F(Form("fNClusters%d",i),"fNClusters0",200, 0,100);
 530           fTPCNcls_EoverP[i]= new TH2F(Form("fTPCNcls_EoverP%d",i),"TPCNcls_EoverP",1000,0,200,200,0,2);
 531
 532
 533           fOutputList->Add(fEoverP_pt[i]);
 534           fOutputList->Add(fTPC_p[i]);
 535           fOutputList->Add(fTPCnsigma_p[i]);
 536
 537
 538           fOutputList->Add(fECluster[i]);
 539           fOutputList->Add(fEtaPhi[i]);
 540           fOutputList->Add(fVtxZ[i]);
 541           fOutputList->Add(fNTracks[i]);
 542           fOutputList->Add(fNClusters[i]);
 543           fOutputList->Add(fTPCNcls_EoverP[i]);
 544         }
 545
 546         //pt bin
 547         Int_t fPtBin[6] = {2,4,6,8,10,15};
 548
 549         fEoverP_tpc = new TH2F *[5];
 550         fTPC_pt = new TH1F *[5];
 551         fTPCnsigma_pt = new TH1F *[5];
 552
 553         fEta=new TH1F *[5];
 554         fPhi=new TH1F *[5];
 555         fR=new TH1F *[5];
 556         fR_EoverP=new TH2F *[5];
 557         fNcells=new TH1F *[5];
 558         fNcells_EoverP=new TH2F *[5];
 559         fM02_EoverP= new TH2F *[5];
 560         fM20_EoverP= new TH2F *[5];
 561         fEoverP_ptbins=new TH1F *[5];
 562         fECluster_ptbins=new TH1F *[5];
 563         fEoverP_wSSCut=new TH1F *[5];
 564         fNcells_electrons=new TH1F *[5];
 565         fNcells_hadrons=new TH1F *[5];
 566         fTPCnsigma_eta_electrons=new TH2F *[5];
 567         fTPCnsigma_eta_hadrons=new TH2F *[5];
 568
 569         if(fCorrelationFlag)
 570         {
 571                 fCEtaPhi_Inc = new TH2F *[5];
 572                 fCEtaPhi_Inc_DiHadron = new TH2F *[5];
 573
 574                 fCEtaPhi_ULS = new TH2F *[5];
 575                 fCEtaPhi_LS = new TH2F *[5];
 576                 fCEtaPhi_ULS_NoP = new TH2F *[5];
 577                 fCEtaPhi_LS_NoP = new TH2F *[5];
 578
 579                 fCEtaPhi_ULS_Weight = new TH2F *[5];
 580                 fCEtaPhi_LS_Weight = new TH2F *[5];
 581                 fCEtaPhi_ULS_NoP_Weight = new TH2F *[5];
 582                 fCEtaPhi_LS_NoP_Weight = new TH2F *[5];
 583
 584                 fCEtaPhi_Inc_EM = new TH2F *[5];
 585
 586                 fCEtaPhi_ULS_EM = new TH2F *[5];
 587                 fCEtaPhi_LS_EM = new TH2F *[5];
 588
 589                 fCEtaPhi_ULS_Weight_EM = new TH2F *[5];
 590                 fCEtaPhi_LS_Weight_EM = new TH2F *[5];
 591
 592                 fInvMass = new TH1F("fInvMass","",200,0,0.3);
 593                 fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
 594                 fDCA = new TH1F("fDCA","",200,0,1);
 595                 fDCABack = new TH1F("fDCABack","",200,0,1);
 596                 fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
 597                 fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
 598
 599                 fOutputList->Add(fInvMass);
 600                 fOutputList->Add(fInvMassBack);
 601                 fOutputList->Add(fDCA);
 602                 fOutputList->Add(fDCABack);
 603                 fOutputList->Add(fOpAngle);
 604                 fOutputList->Add(fOpAngleBack);
 605         }
 606
 607         for(Int_t i = 0; i < 5; i++)
 608         {
 609           fEoverP_tpc[i] = new TH2F(Form("fEoverP_tpc%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;E / p )",fPtBin[i],fPtBin[i+1]),1000,-15,15,100,0,2);
 610           fTPC_pt[i] = new TH1F(Form("fTPC_pt%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;Count",fPtBin[i],fPtBin[i+1]),200,20,200);
 611           fTPCnsigma_pt[i] = new TH1F(Form("fTPCnsigma_pt%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;Count",fPtBin[i],fPtBin[i+1]),200,-15,10);
 612
 613           fEta[i]=new TH1F(Form("fEta%d",i), Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma; dEta )",fPtBin[i],fPtBin[i+1]),100, -0.1,0.1);
 614           fPhi[i]=new TH1F(Form("fPhi%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma; dPhi )",fPtBin[i],fPtBin[i+1]), 100, -0.1,0.1);
 615           fR[i]=new TH1F(Form("fR%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma; dR )",fPtBin[i],fPtBin[i+1]), 100, -0.1,0.1);
 616           fR_EoverP[i]=new TH2F(Form("fR_EoverP%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;E / p )",fPtBin[i],fPtBin[i+1]),100, 0,0.1,1000,0,10);
 617           fNcells[i]=new TH1F(Form("fNcells%d",i), Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;E / p )",fPtBin[i],fPtBin[i+1]),100, 0, 30);
 618           fNcells_electrons[i]=new TH1F(Form("fNcells_electrons%d",i), Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma )",fPtBin[i],fPtBin[i+1]),100, 0, 30);
 619           fNcells_hadrons[i]=new TH1F(Form("fNcells_hadrons%d",i), Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;  Ncells; E / p )",fPtBin[i],fPtBin[i+1]),100, 0, 30);
 620           fNcells_EoverP[i]=new TH2F(Form("fNcells_EoverP%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma; Ncells; E / p )",fPtBin[i],fPtBin[i+1]),1000, 0,20,100,0,30);
 621           fM02_EoverP[i]= new TH2F(Form("fM02_EoverP%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma; M02; E / p )",fPtBin[i],fPtBin[i+1]),1000,0,1,100,0,2);
 622           fM20_EoverP[i]= new TH2F(Form("fM20_EoverP%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma; M20; E / p )",fPtBin[i],fPtBin[i+1]),1000,0,1,100,0,2);
 623           fEoverP_ptbins[i] = new TH1F(Form("fEoverP_ptbins%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;E / p )",fPtBin[i],fPtBin[i+1]),500,0,2);
 624           fECluster_ptbins[i]= new TH1F(Form("fECluster_ptbins%d",i), Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;E / p )",fPtBin[i],fPtBin[i+1]),2000, 0,100);
 625       fEoverP_wSSCut[i]=new TH1F(Form("fEoverP_wSSCut%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;E / p )",fPtBin[i],fPtBin[i+1]),500,0,2);
 626           fTPCnsigma_eta_electrons[i]=new TH2F(Form("fTPCnsigma_eta_electrons%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;Eta )",fPtBin[i],fPtBin[i+1]),1000,-15,15,100,-1,1);
 627           fTPCnsigma_eta_hadrons[i]=new TH2F(Form("fTPCnsigma_eta_hadrons%d",i),Form("%d < p_{t} < %d GeV/c;TPC Electron N#sigma;Eta )",fPtBin[i],fPtBin[i+1]),1000,-15,15,100,-1,1);
 628
 629           fOutputList->Add(fEoverP_tpc[i]);
 630           fOutputList->Add(fTPC_pt[i]);
 631           fOutputList->Add(fTPCnsigma_pt[i]);
 632
 633           fOutputList->Add(fEta[i]);
 634           fOutputList->Add(fPhi[i]);
 635           fOutputList->Add(fR[i]);
 636           fOutputList->Add(fR_EoverP[i]);
 637           fOutputList->Add(fNcells[i]);
 638           fOutputList->Add(fNcells_electrons[i]);
 639       fOutputList->Add(fNcells_hadrons[i]);
 640           fOutputList->Add(fNcells_EoverP[i]);
 641           fOutputList->Add(fECluster_ptbins[i]);
 642           fOutputList->Add(fEoverP_ptbins[i]);
 643           fOutputList->Add(fEoverP_wSSCut[i]);
 644           fOutputList->Add(fM02_EoverP[i]);
 645           fOutputList->Add(fM20_EoverP[i]);
 646           fOutputList->Add(fTPCnsigma_eta_electrons[i]);
 647           fOutputList->Add(fTPCnsigma_eta_hadrons[i]);
 648
 649
 650                 if(fCorrelationFlag)
 651                 {
 652                         fCEtaPhi_Inc[i] = new TH2F(Form("fCEtaPhi_Inc%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 653                         fCEtaPhi_Inc_DiHadron[i] = new TH2F(Form("fCEtaPhi_Inc_DiHadron%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 654
 655                         fCEtaPhi_ULS[i] = new TH2F(Form("fCEtaPhi_ULS%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 656                         fCEtaPhi_LS[i] = new TH2F(Form("fCEtaPhi_LS%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 657                         fCEtaPhi_ULS_NoP[i] = new TH2F(Form("fCEtaPhi_ULS_NoP%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 658                         fCEtaPhi_LS_NoP[i] = new TH2F(Form("fCEtaPhi_LS_NoP%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 659
 660                         fCEtaPhi_ULS_Weight[i] = new TH2F(Form("fCEtaPhi_ULS_Weight%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 661                         fCEtaPhi_LS_Weight[i] = new TH2F(Form("fCEtaPhi_LS_Weight%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 662                         fCEtaPhi_ULS_NoP_Weight[i] = new TH2F(Form("fCEtaPhi_ULS_NoP_Weight%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 663                         fCEtaPhi_LS_NoP_Weight[i] = new TH2F(Form("fCEtaPhi_LS_NoP_Weight%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 664
 665                         fOutputList->Add(fCEtaPhi_Inc[i]);
 666                         fOutputList->Add(fCEtaPhi_Inc_DiHadron[i]);
 667
 668                         fOutputList->Add(fCEtaPhi_ULS[i]);
 669                         fOutputList->Add(fCEtaPhi_LS[i]);
 670                         fOutputList->Add(fCEtaPhi_ULS_NoP[i]);
 671                         fOutputList->Add(fCEtaPhi_LS_NoP[i]);
 672
 673                         fOutputList->Add(fCEtaPhi_ULS_Weight[i]);
 674                         fOutputList->Add(fCEtaPhi_LS_Weight[i]);
 675                         fOutputList->Add(fCEtaPhi_ULS_NoP_Weight[i]);
 676                         fOutputList->Add(fCEtaPhi_LS_NoP_Weight[i]);
 677
 678                         if(fEventMixingFlag)
 679                         {
 680                                 fCEtaPhi_Inc_EM[i] = new TH2F(Form("fCEtaPhi_Inc_EM%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 681
 682                                 fCEtaPhi_ULS_EM[i] = new TH2F(Form("fCEtaPhi_ULS_EM%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 683                                 fCEtaPhi_LS_EM[i] = new TH2F(Form("fCEtaPhi_LS_EM%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 684
 685                                 fCEtaPhi_ULS_Weight_EM[i] = new TH2F(Form("fCEtaPhi_ULS_Weight_EM%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 686                                 fCEtaPhi_LS_Weight_EM[i] = new TH2F(Form("fCEtaPhi_LS_Weight_EM%d",i),Form("%d < p_{t} < %d GeV/c;#DeltaPhi (rad);#Delta#eta",fPtBin[i],fPtBin[i+1]),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-2,2);
 687
 688                                 fOutputList->Add(fCEtaPhi_Inc_EM[i]);
 689
 690                                 fOutputList->Add(fCEtaPhi_ULS_EM[i]);
 691                                 fOutputList->Add(fCEtaPhi_LS_EM[i]);
 692
 693                                 fOutputList->Add(fCEtaPhi_ULS_Weight_EM[i]);
 694                                 fOutputList->Add(fCEtaPhi_LS_Weight_EM[i]);
 695                         }
 696                 }
 697         }
 698
 699         //pt integrated
 700         fTPCnsigma_eta = new TH2F("fTPCnsigma_eta",";Pseudorapidity #eta; TPC signal - <TPC signal>_{elec} (#sigma)",200,-0.9,0.9,200,-15,15);
 701         fTPCnsigma_phi = new TH2F("fTPCnsigma_phi",";Azimuthal Angle #phi; TPC signal - <TPC signal>_{elec} (#sigma)",200,0,2*TMath::Pi(),200,-15,15);
 702
 703
 704         fNcells_pt=new TH2F("fNcells_pt","fNcells_pt",1000, 0,20,100,0,30);
 705         fEoverP_pt_pions= new TH2F("fEoverP_pt_pions","fEoverP_pt_pions",1000,0,30,500,0,2);
 706
 707         ftpc_p_EoverPcut= new TH2F("ftpc_p_EoverPcut","ftpc_p_EoverPcut",1000,0,30,200,20,200);
 708         fnsigma_p_EoverPcut= new TH2F("fnsigma_p_EoverPcut","fnsigma_p_EoverPcut",1000,0,30,500,-15,15);
 709
 710         fEoverP_pt_pions2= new TH2F("fEoverP_pt_pions2","fEoverP_pt_pions2",1000,0,30,500,0,2);
 711         fEoverP_pt_hadrons= new TH2F("fEoverP_pt_hadrons","fEoverP_pt_hadrons",1000,0,30,500,0,2);
 712
 713
 714         fOutputList->Add(fTPCnsigma_eta);
 715         fOutputList->Add(fTPCnsigma_phi);
 716
 717         fOutputList->Add(fNcells_pt);
 718         fOutputList->Add(fEoverP_pt_pions);
 719
 720         fOutputList->Add(ftpc_p_EoverPcut);
 721         fOutputList->Add(fnsigma_p_EoverPcut);
 722
 723         fOutputList->Add(fEoverP_pt_pions2);
 724         fOutputList->Add(fEoverP_pt_hadrons);
 725
 726         //__________________________________________________________________
 727         //Efficiency studies
 728         if(fIsMC)
 729         {
 730                 fPtBackgroundBeforeReco = new TH1F("fPtBackgroundBeforeReco",";p_{t} (GeV/c);Count",30,0,15);
 731                 fPtBackgroundAfterReco = new TH1F("fPtBackgroundAfterReco",";p_{t} (GeV/c);Count",30,0,15);
 732                 fPtMCparticleAll = new TH1F("fPtMCparticleAll",";p_{t} (GeV/c);Count",200,0,40);
 733                 fPtMCparticleReco = new TH1F("fPtMCparticleReco",";p_{t} (GeV/c);Count",200,0,40);
 734                 fPtMCparticleAllHfe1 = new TH1F("fPtMCparticleAllHfe1",";p_{t} (GeV/c);Count",200,0,40);
 735                 fPtMCparticleRecoHfe1 = new TH1F("fPtMCparticleRecoHfe1",";p_{t} (GeV/c);Count",200,0,40);
 736                 fPtMCparticleAllHfe2 = new TH1F("fPtMCparticleAllHfe2",";p_{t} (GeV/c);Count",200,0,40);
 737                 fPtMCparticleRecoHfe2 = new TH1F("fPtMCparticleRecoHfe2",";p_{t} (GeV/c);Count",200,0,40);
 738                 fPtMCelectronAfterAll = new TH1F("fPtMCelectronAfterAll",";p_{t} (GeV/c);Count",200,0,40);
 739                 fPtMCpi0 = new TH1F("fPtMCpi0",";p_{t} (GeV/c);Count",200,0,40);
 740                 fPtMC_EMCal_All= new TH1F("fPtMC_EMCal_All",";p_{t} (GeV/c);Count",200,0,40);
 741                 fPtMC_EMCal_Selected= new TH1F("fPtMC_EMCal_Selected",";p_{t} (GeV/c);Count",200,0,40);
 742                 fPtMC_TPC_All= new TH1F("fPtMC_TPC_All",";p_{t} (GeV/c);Count",200,0,40);
 743                 fPtMC_TPC_Selected = new TH1F("fPtMC_TPC_Selected",";p_{t} (GeV/c);Count",200,0,40);
 744                 fPtMCWithLabel = new TH1F("fPtMCWithLabel",";p_{t} (GeV/c);Count",200,0,40);
 745                 fPtMCWithoutLabel = new TH1F("fPtMCWithoutLabel",";p_{t} (GeV/c);Count",200,0,40);
 746                 fPtIsPhysicaPrimary = new TH1F("fPtIsPhysicaPrimary",";p_{t} (GeV/c);Count",200,0,40);
 747
 748                 fOutputList->Add(fPtBackgroundBeforeReco);
 749                 fOutputList->Add(fPtBackgroundAfterReco);
 750                 fOutputList->Add(fPtMCparticleAll);
 751                 fOutputList->Add(fPtMCparticleReco);
 752                 fOutputList->Add(fPtMCparticleAllHfe1);
 753                 fOutputList->Add(fPtMCparticleRecoHfe1);
 754                 fOutputList->Add(fPtMCparticleAllHfe2);
 755                 fOutputList->Add(fPtMCparticleRecoHfe2);
 756                 fOutputList->Add(fPtMCelectronAfterAll);
 757                 fOutputList->Add(fPtMCpi0);
 758                 fOutputList->Add(fPtMC_EMCal_All);
 759                 fOutputList->Add(fPtMC_EMCal_Selected);
 760                 fOutputList->Add(fPtMC_TPC_All);
 761                 fOutputList->Add(fPtMC_TPC_Selected);
 762                 fOutputList->Add(fPtMCWithLabel);
 763                 fOutputList->Add(fPtMCWithoutLabel);
 764                 fOutputList->Add(fPtIsPhysicaPrimary);
 765         }
 766
 767         fCentralityHist = new TH1F("fCentralityHist",";Centrality (%); Count",1000000,0,100);
 768         fCentralityHistPass = new TH1F("fCentralityHistPass",";Centrality (%); Count",1000000,0,100);
 769         fOutputList->Add(fCentralityHist);
 770         fOutputList->Add(fCentralityHistPass);
 771
 772         //______________________________________________________________________
 773         //Mixed event analysis
 774         if(fEventMixingFlag)
 775         {
 776                 fPoolNevents = new TH1F("fPoolNevents","Event Mixing Statistics; Number of events; Count",1000,0,1000);
 777                 fOutputList->Add(fPoolNevents);
 778
 779                 Int_t trackDepth = 2000; // number of objects (tracks) kept per event buffer bin. Once the number of stored objects (tracks) is above that limit, the oldest ones are removed.
 780                 Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 781
 782                 Int_t nCentralityBins  = 15;
 783                 Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
 784
 785                 Int_t nZvtxBins  = 9;
 786                 Double_t vertexBins[] = {-10, -7, -5, -3, -1, 1, 3, 5, 7, 10};
 787
 788                 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, (Double_t*) centralityBins, nZvtxBins, (Double_t*) vertexBins);
 789         }
 790         //______________________________________________________________________
 791
 792         PostData(1, fOutputList);
 793
 794 ///______________________________________________________________________
 795 }
 796
 797 //______________________________________________________________________
 798 //Main loop
 799 //Called for each event
 800 void AliAnalysisTaskEMCalHFEpA::UserExec(Option_t *)
 801 {
 802 //Check Event
 803         fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 804         fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 805
 806         if(!(fESD || fAOD))
 807         {
 808                 printf("ERROR: fESD & fAOD not available\n");
 809                 return;
 810         }
 811
 812         fVevent = dynamic_cast<AliVEvent*>(InputEvent());
 813
 814         if(!fVevent)
 815         {
 816                 printf("ERROR: fVEvent not available\n");
 817                 return;
 818         }
 819
 820 //Check Cuts
 821         if(!fCuts)
 822         {
 823                 AliError("HFE cuts not available");
 824                 return;
 825         }
 826 //Check PID
 827         if(!fPID->IsInitialized())
 828         {
 829     // Initialize PID with the given run number
 830                 AliWarning("PID not initialised, get from Run no");
 831
 832                 if(fIsAOD)
 833                 {
 834                         fPID->InitializePID(fAOD->GetRunNumber());
 835                 }
 836                 else
 837                 {
 838                         fPID->InitializePID(fESD->GetRunNumber());
 839                 }
 840         }
 841
 842 //PID response
 843         fPidResponse = fInputHandler->GetPIDResponse();
 844
 845
 846 //Check PID response
 847         if(!fPidResponse)
 848         {
 849                 AliDebug(1, "Using default PID Response");
 850                 fPidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 851         }
 852
 853         fPID->SetPIDResponse(fPidResponse);
 854
 855         fCFM->SetRecEventInfo(fVevent);
 856
 857         Double_t *fListOfmotherkink = 0;
 858         Int_t fNumberOfVertices = 0;
 859         Int_t fNumberOfMotherkink = 0;
 860
 861 //______________________________________________________________________
 862 //Vertex Selection
 863         if(fIsAOD)
 864         {
 865                 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
 866                 if(!trkVtx || trkVtx->GetNContributors()<=0) return;
 867                 TString vtxTtl = trkVtx->GetTitle();
 868                 if(!vtxTtl.Contains("VertexerTracks")) return;
 869                 Float_t zvtx = trkVtx->GetZ();
 870                 fZvtx = zvtx;
 871                 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
 872                 if(spdVtx->GetNContributors()<=0) return;
 873                 TString vtxTyp = spdVtx->GetTitle();
 874                 Double_t cov[6]={0};
 875                 spdVtx->GetCovarianceMatrix(cov);
 876                 Double_t zRes = TMath::Sqrt(cov[5]);
 877                 if(vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
 878                 if(TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
 879                 if(TMath::Abs(zvtx) > 10) return;
 880
 881         //Look for kink mother for AOD
 882
 883                 fNumberOfVertices = 0;
 884                 fNumberOfMotherkink = 0;
 885
 886                 if(fIsAOD)
 887                 {
 888                         fNumberOfVertices = fAOD->GetNumberOfVertices();
 889
 890                         fListOfmotherkink = new Double_t[fNumberOfVertices];
 891
 892                         for(Int_t ivertex=0; ivertex < fNumberOfVertices; ivertex++)
 893                         {
 894                                 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
 895                                 if(!aodvertex) continue;
 896                                 if(aodvertex->GetType()==AliAODVertex::kKink)
 897                                 {
 898                                         AliAODTrack *mother1 = (AliAODTrack *) aodvertex->GetParent();
 899                                         if(!mother1) continue;
 900                                         Int_t idmother = mother1->GetID();
 901                                         fListOfmotherkink[fNumberOfMotherkink] = idmother;
 902                                         fNumberOfMotherkink++;
 903                                 }
 904                         }
 905                 }
 906         }
 907         else
 908         {
 909                 /// (has to be checked)
 910
 911                 const AliESDVertex *trkVtx = fESD->GetPrimaryVertex();
 912                 Float_t zvtx = trkVtx->GetZ();
 913                 fZvtx = zvtx;
 914                 if(TMath::Abs(zvtx) > 10) return;
 915
 916                 //Double_t NTracks = fESD->GetNumberOfTracks();
 917
 918                 /*const AliVVertex* trkVtx = fVevent->GetPrimaryVertex();
 919                 if(!trkVtx || trkVtx->GetNContributors()<=0) return;
 920                 TString vtxTtl = trkVtx->GetTitle();
 921                 if(!vtxTtl.Contains("VertexerTracks")) return;
 922                 Float_t zvtx = trkVtx->GetZ();
 923                 const AliVVertex* spdVtx = fVevent->GetPrimaryVertexSPD();
 924                 if(spdVtx->GetNContributors()<=0) return;
 925                 TString vtxTyp = spdVtx->GetTitle();
 926                 Double_t cov[6]={0};
 927                 spdVtx->GetCovarianceMatrix(cov);
 928                 Double_t zRes = TMath::Sqrt(cov[5]);
 929                 if(vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
 930                 if(TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
 931                 if(TMath::Abs(zvtx) > 10) return;*/
 932         }
 933
 934 //______________________________________________________________________
 935
 936 //Only events with at least 2 tracks are accepted
 937         Int_t fNOtrks =  fVevent->GetNumberOfTracks();
 938         if(fNOtrks<2) return;
 939
 940 //______________________________________________________________________
 941 //Centrality Selection
 942         if(fHasCentralitySelection)
 943         {
 944                 Float_t centrality = -1;
 945
 946                 if(fIsAOD)
 947                 {
 948                         fCentrality = fAOD->GetHeader()->GetCentralityP();
 949                 }
 950                 else
 951                 {
 952                         fCentrality = fESD->GetCentrality();
 953                 }
 954
 955                 if(fEstimator==1) centrality = fCentrality->GetCentralityPercentile("ZDC");
 956                 else centrality = fCentrality->GetCentralityPercentile("V0A");
 957
 958                 //cout << "Centrality = " << centrality << " %" << endl;
 959
 960                 fCentralityHist->Fill(centrality);
 961
 962                 if(centrality<fCentralityMin || centrality>fCentralityMax) return;
 963
 964                 fCentralityHistPass->Fill(centrality);
 965         }
 966 //______________________________________________________________________
 967
 968 //______________________________________________________________________
 969
 970         if(fIsMC)
 971         {
 972                 if(fIsAOD)
 973                 {
 974                         fMCarray = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 975
 976                         if(!fMCarray)
 977                         {
 978                                 AliError("Array of MC particles not found");
 979                                 return;
 980                         }
 981
 982                         fMCheader = dynamic_cast<AliAODMCHeader*>(fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 983
 984                         if(!fMCheader)
 985                         {
 986                                 AliError("Could not find MC Header in AOD");
 987                                 return;
 988                         }
 989
 990                         for(Int_t iMC = 0; iMC < fMCarray->GetEntries(); iMC++)
 991                         {
 992                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(iMC);
 993
 994                                 Int_t pdg = fMCparticle->GetPdgCode();
 995                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax && fMCparticle->Charge()!=0)
 996                                 {
 997                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
 998                                         {
 999                                                 if(fMCparticle->IsPhysicalPrimary())
1000                                                 {
1001                                                         fPtMCparticleAll->Fill(fMCparticle->Pt());
1002
1003                                                         Bool_t MotherFound = FindMother(iMC);
1004                                                         if(MotherFound)
1005                                                         {
1006                                                                         if(fIsHFE1) fPtMCparticleAllHfe1->Fill(fMCparticle->Pt()); //denominator for total efficiency
1007                                                                         if(fIsHFE2) fPtMCparticleAllHfe2->Fill(fMCparticle->Pt());
1008                                                         }
1009                                                 }
1010                                         }
1011                                 }
1012                                 if(TMath::Abs(pdg)==111) fPtMCpi0->Fill(fMCparticle->Pt());
1013                         }
1014                 }
1015                 else
1016                 {
1017                         fEventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1018                         if (!fEventHandler) {
1019                                                         Printf("ERROR: Could not retrieve MC event handler");
1020                                                         return;
1021                         }
1022
1023                         fMCevent = fEventHandler->MCEvent();
1024                         if (!fMCevent) {
1025                                                         Printf("ERROR: Could not retrieve MC event");
1026                                                         return;
1027                         }
1028
1029                         fMCstack = fMCevent->Stack();
1030
1031                 for(Int_t iMC = 0; iMC < fMCstack->GetNtrack(); iMC++)
1032                 {
1033
1034                                 fMCtrack = fMCstack->Particle(iMC);
1035                                 Int_t pdg = fMCtrack->GetPdgCode();
1036                                 if(TMath::Abs(pdg)==111) fPtMCpi0->Fill(fMCtrack->Pt());
1037
1038                                 if(!fMCstack->IsPhysicalPrimary(iMC)) continue;
1039
1040
1041                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
1042                                 {
1043                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1044                                         {
1045                                                 fPtMCparticleAll->Fill(fMCtrack->Pt());
1046
1047                                                 Bool_t MotherFound = FindMother(iMC);
1048                                                 if(MotherFound)
1049                                                 {
1050                                                         if(fIsHFE1) fPtMCparticleAllHfe1->Fill(fMCtrack->Pt());//denominator for total efficiency
1051                                                         if(fIsHFE2) fPtMCparticleAllHfe2->Fill(fMCtrack->Pt());
1052                                                 }
1053                                         }
1054                                 }
1055                 }
1056                 }
1057         }
1058
1059 //______________________________________________________________________
1060 //EMCal Trigger Selection (Threshould selection)
1061         TString firedTrigger;
1062         TString TriggerEG1("EG1");
1063         TString TriggerEG2("EG2");
1064
1065         if(fAOD) firedTrigger = fAOD->GetFiredTriggerClasses();
1066         else if(fESD) firedTrigger = fESD->GetFiredTriggerClasses();
1067
1068         fNevent->Fill(0);
1069         if(firedTrigger.Contains(TriggerEG1)) fNevent->Fill(1);
1070         if(firedTrigger.Contains(TriggerEG2)) fNevent->Fill(2);
1071
1072         //EG1
1073         if(firedTrigger.Contains(TriggerEG1))
1074         {
1075                 fNevent->Fill(3);
1076         }
1077         else
1078         {
1079                 if(fEMCEG1) return;
1080         }
1081
1082         //EG2
1083         if(firedTrigger.Contains(TriggerEG2))
1084         {
1085                 fNevent->Fill(4);
1086         }
1087         else
1088         {
1089                 if(fEMCEG2) return;
1090         }
1091
1092 //______________________________________________________________________
1093
1094         Int_t ClsNo = -999;
1095         if(!fIsAOD) ClsNo = fESD->GetNumberOfCaloClusters();
1096         else ClsNo = fAOD->GetNumberOfCaloClusters();
1097
1098 //______________________________________________________________________
1099
1100 ///______________________________________________________________________
1101 ///Track loop
1102         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
1103         {
1104                 AliVParticle* Vtrack = fVevent->GetTrack(iTracks);
1105                 if (!Vtrack)
1106                 {
1107                         printf("ERROR: Could not receive track %d\n", iTracks);
1108                         continue;
1109                 }
1110
1111                 AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
1112                 AliESDtrack *etrack = dynamic_cast<AliESDtrack*>(Vtrack);
1113                 AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(Vtrack);
1114
1115                 Double_t fTPCnSigma = -999;
1116                 Double_t fTPCnSigma_pion = -999;
1117                 Double_t fTPCnSigma_proton = -999;
1118                 Double_t fTPCnSigma_kaon = -999;
1119                 Double_t fTPCsignal = -999;
1120                 Double_t fPt = -999;
1121                 Double_t fP = -999;
1122
1123                 ///_____________________________________________________________________________
1124                 ///Fill QA plots without track selection
1125                 fPt = track->Pt();
1126                 fP = TMath::Sqrt((track->Pt())*(track->Pt()) + (track->Pz())*(track->Pz()));
1127
1128                 fTPCsignal = track->GetTPCsignal();
1129                 fTPCnSigma = fPidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
1130                 fTPCnSigma_pion = fPidResponse->NumberOfSigmasTPC(track, AliPID::kPion);
1131                 fTPCnSigma_proton = fPidResponse->NumberOfSigmasTPC(track, AliPID::kProton);
1132                 fTPCnSigma_kaon = fPidResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
1133
1134                 fTPC_p[0]->Fill(fPt,fTPCsignal);
1135                 fTPCnsigma_p[0]->Fill(fP,fTPCnSigma);
1136
1137
1138                 Float_t TPCNcls = track->GetTPCNcls();
1139                 Float_t pos[3]={0,0,0};
1140
1141                 Double_t fEMCflag = kFALSE;
1142                 if(track->GetEMCALcluster()>0)
1143                 {
1144                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
1145                         if(fClus->IsEMCAL())
1146                         {
1147                           if(TMath::Abs(fClus->GetTrackDx())<=0.05 && TMath::Abs(fClus->GetTrackDz())<=0.05)
1148                           {
1149                               fEMCflag = kTRUE;
1150                               fEoverP_pt[0]->Fill(fPt,(fClus->E() / fP));
1151
1152
1153                                   Float_t Energy        = fClus->E();
1154                                   Float_t EoverP        = Energy/track->P();
1155                                   //Float_t M02 = fClus->GetM02();
1156                                   //Float_t M20 = fClus->GetM20();
1157
1158                                   /////////////// for Eta Phi distribution
1159                                   fClus->GetPosition(pos);
1160                                   TVector3 vpos(pos[0],pos[1],pos[2]);
1161                                   Double_t cphi = vpos.Phi();
1162                                   Double_t ceta = vpos.Eta();
1163                                   fEtaPhi[0]->Fill(cphi,ceta);
1164
1165                                   fECluster[0]->Fill(Energy);
1166                                   fTPCNcls_EoverP[0]->Fill(TPCNcls, EoverP);
1167                           }
1168                         }
1169                 }
1170
1171                 //______________________________________________________________
1172                 // Vertex
1173
1174                 fVtxZ[0]->Fill(fZvtx);
1175                 fNTracks[0]->Fill(fNOtrks);
1176                 fNClusters[0]->Fill(ClsNo);
1177
1178                 //______________________________________________________________
1179
1180                 ///Fill QA plots without track selection
1181                 ///_____________________________________________________________________________
1182 //______________________________________________________________________________________
1183 //Track Selection Cuts
1184
1185 //AOD (Test Filter Bit)
1186                 if(fIsAOD)
1187                 {
1188                         // standard cuts with very loose DCA - BIT(4)
1189                         // Description:
1190                         /*
1191                                 GetStandardITSTPCTrackCuts2011(kFALSE)
1192                                     SetMaxChi2PerClusterTPC(4);
1193                                     SetAcceptKinkDaughters(kFALSE);
1194                                     SetRequireTPCRefit(kTRUE);
1195                                     SetRequireITSRefit(kTRUE);
1196                                     SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
1197                                     SetMaxDCAToVertexZ(2);
1198                                     SetMaxDCAToVertex2D(kFALSE);
1199                                     SetRequireSigmaToVertex(kFALSE);
1200                                     SetMaxChi2PerClusterITS(36);
1201                                 SetMaxDCAToVertexXY(2.4)
1202                                 SetMaxDCAToVertexZ(3.2)
1203                                 SetDCaToVertex2D(kTRUE)
1204                         */
1205
1206                         if(!atrack->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;
1207                 }
1208
1209 //RecKine: ITSTPC cuts
1210                 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
1211 //RecKink
1212                 if(fRejectKinkMother)
1213                 {
1214                         if(fIsAOD)
1215                         {
1216                                 Bool_t kinkmotherpass = kTRUE;
1217                                 for(Int_t kinkmother = 0; kinkmother < fNumberOfMotherkink; kinkmother++)
1218                                 {
1219                                         if(track->GetID() == fListOfmotherkink[kinkmother])
1220                                         {
1221                                                 kinkmotherpass = kFALSE;
1222                                                 continue;
1223                                         }
1224                                 }
1225                                 if(!kinkmotherpass) continue;
1226                         }
1227                         else
1228                         {
1229                                 if(etrack->GetKinkIndex(0) != 0) continue;
1230                         }
1231                 }
1232
1233 //RecPrim
1234                 if(!fIsAOD)
1235                 {
1236                         if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
1237                 }
1238
1239 //HFEcuts: ITS layers cuts
1240                 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
1241
1242 //HFE cuts: TPC PID cleanup
1243                 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
1244 //______________________________________________________________________________________
1245
1246                 ///_____________________________________________________________
1247                 ///QA plots after track selection
1248                 if(fIsMC)
1249                 {
1250                         if(track->GetLabel()>=0) fPtMCWithLabel->Fill(fPt);
1251                         else fPtMCWithoutLabel->Fill(fPt);
1252                 }
1253
1254                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
1255                 {
1256                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1257
1258                         if(fMCparticle->IsPhysicalPrimary()) fPtIsPhysicaPrimary->Fill(fPt);
1259
1260                         Int_t pdg = fMCparticle->GetPdgCode();
1261                         if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax && fMCparticle->Charge()!=0)
1262                         {
1263                                 if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1264                                 {
1265                                         if(fMCparticle->IsPhysicalPrimary())
1266                                         {
1267                                                 fPtMCparticleReco->Fill(fMCparticle->Pt());
1268
1269                                                 Bool_t MotherFound = FindMother(track->GetLabel());
1270                                                 if(MotherFound)
1271                                                 {
1272                                                         if(fIsHFE1) fPtMCparticleRecoHfe1->Fill(fMCparticle->Pt());
1273                                                         if(fIsHFE2) fPtMCparticleRecoHfe2->Fill(fMCparticle->Pt());
1274                                                 }
1275                                         }
1276                                 }
1277                         }
1278                 }
1279                 else if(fIsMC && track->GetLabel()>=0)
1280                 {
1281                         if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
1282                         {
1283                                 fPtIsPhysicaPrimary->Fill(fPt);
1284                                 fMCtrack = fMCstack->Particle(track->GetLabel());
1285
1286                                 Int_t pdg = fMCtrack->GetPdgCode();
1287                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
1288                                 {
1289                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1290                                         {
1291                                                 fPtMCparticleReco->Fill(fMCtrack->Pt());
1292
1293                                                 Bool_t MotherFound = FindMother(track->GetLabel());
1294                                                 if(MotherFound)
1295                                                 {
1296                                                         if(fIsHFE1) fPtMCparticleRecoHfe1->Fill(fMCtrack->Pt());
1297                                                         if(fIsHFE2) fPtMCparticleRecoHfe2->Fill(fMCtrack->Pt());
1298                                                 }
1299                                         }
1300                                 }
1301                         }
1302                 }
1303
1304                 fTPC_p[1]->Fill(fPt,fTPCsignal);
1305                 fTPCnsigma_p[1]->Fill(fP,fTPCnSigma);
1306                 Double_t fPtBin[6] = {2,4,6,8,10,15};
1307
1308                 TPCNcls = track->GetTPCNcls();
1309                 Float_t pos2[3]={0,0,0};
1310
1311                 if(track->GetEMCALcluster()>0)
1312                 {
1313                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
1314                         if(fClus->IsEMCAL())
1315                         {
1316                           if(TMath::Abs(fClus->GetTrackDx())<=0.05 && TMath::Abs(fClus->GetTrackDz())<=0.05)
1317                           {
1318                               fEoverP_pt[1]->Fill(fPt,(fClus->E() / fP));
1319
1320                                   Float_t Energy        = fClus->E();
1321                                   Float_t EoverP        = Energy/track->P();
1322                                   //Float_t M02 = fClus->GetM02();
1323                                   //Float_t M20 = fClus->GetM20();
1324                                   Float_t ncells        = fClus->GetNCells();
1325
1326                                   /////////////// for Eta Phi distribution
1327                                   fClus->GetPosition(pos2);
1328                                   TVector3 vpos(pos2[0],pos2[1],pos2[2]);
1329                                   Double_t cphi = vpos.Phi();
1330                                   Double_t ceta = vpos.Eta();
1331                                   fEtaPhi[1]->Fill(cphi,ceta);
1332
1333                                   fECluster[1]->Fill(Energy);
1334                                   fTPCNcls_EoverP[1]->Fill(TPCNcls, EoverP);
1335
1336                                   // 31/05/2103
1337                                   //for EMCal triggger performance
1338                                   if(EoverP > 0.9){
1339                                           ftpc_p_EoverPcut->Fill(track->P(), fTPCsignal);
1340                                           fnsigma_p_EoverPcut->Fill(track->P(), fTPCnSigma);
1341
1342                                   }
1343
1344
1345                                   //for hadron contamination calculations
1346
1347                                   // 31/05/2103
1348                                   // EtaCut -> dados
1349                                   if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
1350                                          //main
1351                                           if(TMath::Abs(fTPCnSigma_pion)<3 || TMath::Abs(fTPCnSigma_proton)<3 || TMath::Abs(fTPCnSigma_kaon)<3 ){
1352
1353                                                   if(fTPCnSigma<-3.5){
1354                                                           fEoverP_pt_hadrons->Fill(fPt,EoverP);
1355                                                   }
1356                                           }
1357
1358                                           if(fTPCnSigma < -3.5){
1359                                                   fEoverP_pt_pions->Fill(fPt, EoverP);
1360
1361                                           }
1362
1363                                           if(fTPCnSigma < -3.5 && fTPCnSigma > -10){
1364                                                   fEoverP_pt_pions2->Fill(fPt, EoverP);
1365
1366                                           }
1367                                   }
1368
1369                                           for(Int_t i = 0; i < 5; i++)
1370                                           {
1371                                                   if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
1372                                                   {
1373
1374                                                            if(fTPCnSigma < -5 && fTPCnSigma > -10){
1375                                                                         fNcells_hadrons[i]->Fill(ncells);
1376                                                        }
1377                                                           //hadrons selection using E/p
1378                                                            if((fClus->E() / fP) >= 0.2 && (fClus->E() / fP) <= 0.7){
1379                                                                  fTPCnsigma_eta_hadrons[i]->Fill(fTPCnSigma, ceta);
1380                                                            }
1381                                                           //electrons selection using E/p
1382                                                            if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax) {
1383                                                                   fTPCnsigma_eta_electrons[i]->Fill(fTPCnSigma, ceta);
1384                                                            }
1385                                                }
1386                                      }
1387
1388                               if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax)
1389                               {
1390                                                 fTPCnsigma_eta->Fill(track->Eta(),fTPCnSigma);
1391                                                 fTPCnsigma_phi->Fill(track->Phi(),fTPCnSigma);
1392
1393                                                 if(fUseEMCal)
1394                                                 {
1395                                                         if(fTPCnSigma < 3.5 && fCorrelationFlag)
1396                                                         {
1397                                                                 fPtTrigger_Inc->Fill(fPt);
1398                                                                 DiHadronCorrelation(track, iTracks);
1399                                                         }
1400                                                 }
1401                                   }
1402
1403                           }
1404                         }
1405                 }
1406
1407                 //______________________________________________________________
1408                 // Vertex
1409
1410                 fVtxZ[1]->Fill(fZvtx);
1411                 fNTracks[1]->Fill(fNOtrks);
1412                 fNClusters[1]->Fill(ClsNo);
1413                 //______________________________________________________________
1414
1415                 ///______________________________________________________________________
1416                 ///Histograms for PID Studies
1417                 //Double_t fPtBin[6] = {2,4,6,8,10,15};
1418
1419                 for(Int_t i = 0; i < 5; i++)
1420                 {
1421                   if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
1422                   {
1423                           if(fEMCflag) fEoverP_tpc[i]->Fill(fTPCnSigma,(fClus->E() / fP));
1424
1425
1426                     fTPC_pt[i]->Fill(fTPCsignal);
1427                     fTPCnsigma_pt[i]->Fill(fTPCnSigma);
1428
1429                   }
1430                 }
1431
1432                 ///QA plots after track selection
1433                 ///_____________________________________________________________
1434
1435                 //_______________________________________________________
1436                 //Correlation Analysis - DiHadron
1437                 if(!fUseEMCal)
1438                 {
1439                         if(fTPCnSigma < 3.5 && fCorrelationFlag)
1440                         {
1441                                 fPtTrigger_Inc->Fill(fPt);
1442                                 DiHadronCorrelation(track, iTracks);
1443                         }
1444                 }
1445                 //_______________________________________________________
1446
1447
1448                 //!//////////////////////////////////////////////////////////////////
1449                 ///TPC - efficiency calculation //
1450
1451                 /// changing start here
1452                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
1453                 {
1454                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1455                         Int_t pdg = fMCparticle->GetPdgCode();
1456
1457                         //
1458                         if(fMCparticle->IsPhysicalPrimary()){
1459
1460
1461                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
1462
1463                                         Bool_t MotherFound = FindMother(track->GetLabel());
1464                                         if(MotherFound){
1465                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1466                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
1467                                                if(fIsHFE1) fPtMC_TPC_All->Fill(fMCparticle->Pt());
1468                                         }
1469                                         }
1470                              }
1471                         }
1472                 }///until here
1473
1474                 else if(fIsMC && track->GetLabel()>=0)//ESD
1475                 {
1476
1477                         if(fMCstack->IsPhysicalPrimary(track->GetLabel())){
1478                                 fMCtrack = fMCstack->Particle(track->GetLabel());
1479
1480                                 Int_t pdg = fMCtrack->GetPdgCode();
1481                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax ){
1482
1483                                         if(fMCtrack->GetFirstMother()>0){
1484                                             fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1485                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
1486                                                         if(fIsHFE1) fPtMC_TPC_All->Fill(fMCtrack->Pt());
1487                                         }
1488                                         }
1489                                 }
1490                         }
1491                 }
1492
1493
1494 ///________________________________________________________________________
1495 ///PID
1496 ///Here the PID cuts defined in the file "ConfigEMCalHFEpA.C" is applied
1497             Int_t pidpassed = 1;
1498             AliHFEpidObject hfetrack;
1499             hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1500             hfetrack.SetRecTrack(track);
1501             hfetrack.SetPP();   //proton-proton analysis
1502             if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
1503             fpid->Fill(pidpassed);
1504
1505             if(pidpassed==0) continue;
1506 ///________________________________________________________________________
1507
1508
1509                 //!//////////////////////////////////////////////////////////////////
1510                 ///TPC efficiency calculations
1511
1512                 /// changing start here
1513                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
1514                 {
1515                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1516                         Int_t pdg = fMCparticle->GetPdgCode();
1517
1518                         //
1519                         if(fMCparticle->IsPhysicalPrimary()){
1520
1521
1522                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
1523
1524                                         Bool_t MotherFound = FindMother(track->GetLabel());
1525                                         if(MotherFound){
1526                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1527                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
1528                                                         if(fIsHFE1) fPtMC_TPC_Selected->Fill(fMCparticle->Pt());
1529                                         }
1530                                         }
1531                                 }
1532                         }
1533                 }///until here
1534
1535                 else if(fIsMC && track->GetLabel()>=0)//ESD
1536                 {
1537
1538                         if(fMCstack->IsPhysicalPrimary(track->GetLabel())){
1539                                 fMCtrack = fMCstack->Particle(track->GetLabel());
1540
1541                                 Int_t pdg = fMCtrack->GetPdgCode();
1542                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax ){
1543
1544                                         if(fMCtrack->GetFirstMother()>0){
1545                                             fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1546                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
1547                                                         if(fIsHFE1) fPtMC_TPC_Selected->Fill(fMCtrack->Pt());
1548                                         }
1549                                         }
1550                                 }
1551                         }
1552                 }
1553
1554                 //Eta Cut for TPC only
1555                 if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
1556                         fTPCnsigma_pt_2D->Fill(fPt,fTPCnSigma);
1557                 }
1558
1559                 fTPCnsigma_p[2]->Fill(fP,fTPCnSigma);
1560                 fTPC_p[2]->Fill(fP,fTPCsignal);
1561                 TPCNcls = track->GetTPCNcls();
1562                 Float_t pos3[3]={0,0,0};
1563
1564                 if(track->GetEMCALcluster()>0)
1565                 {
1566                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
1567                         if(fClus->IsEMCAL())
1568                         {
1569
1570
1571                           /////////////// Residuals
1572                           Double_t Dx = fClus->GetTrackDx();
1573                           Double_t Dz = fClus->GetTrackDz();
1574                           Double_t R=TMath::Sqrt(Dx*Dx+Dz*Dz);
1575                           for(Int_t i = 0; i < 5; i++)
1576                                 {
1577                                         if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
1578                                         {
1579
1580                                                 fEta[i]->Fill(Dz);
1581                                                 fPhi[i]->Fill(Dx);
1582                                                 fR[i]->Fill(R);
1583                                         }
1584                           }
1585
1586                           if(TMath::Abs(fClus->GetTrackDx())<=0.05 && TMath::Abs(fClus->GetTrackDz())<=0.05)
1587                           {
1588
1589                                   // EtaCut
1590                                   if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
1591                                           fEoverP_pt[2]->Fill(fPt,(fClus->E() / fP));
1592                                   }
1593
1594                                   Float_t Energy        = fClus->E();
1595                                   Float_t EoverP        = Energy/track->P();
1596                                   Float_t M02   = fClus->GetM02();
1597                                   Float_t M20   = fClus->GetM20();
1598                                   Float_t ncells        = fClus->GetNCells();
1599
1600
1601                                   /////////////// for Eta Phi distribution
1602                                   fClus->GetPosition(pos3);
1603                                   TVector3 vpos(pos3[0],pos3[1],pos3[2]);
1604                                   Double_t cphi = vpos.Phi();
1605                                   Double_t ceta = vpos.Eta();
1606                                   fEtaPhi[2]->Fill(cphi,ceta);
1607
1608
1609
1610                                   fTPCNcls_EoverP[2]->Fill(TPCNcls, EoverP);
1611
1612                                   for(Int_t i = 0; i < 5; i++)
1613                                   {
1614                                           if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
1615                                           {
1616
1617                                                   fR_EoverP[i]->Fill(R, EoverP);
1618                                                   fNcells[i]->Fill(ncells);
1619                                                   fNcells_EoverP[i]->Fill(EoverP, ncells);
1620                                                   fM02_EoverP[i]->Fill(M02,EoverP);
1621                                                   fM20_EoverP[i]->Fill(M20,EoverP);
1622                                                   fECluster_ptbins[i]->Fill(Energy);
1623                                                   fEoverP_ptbins[i]->Fill(EoverP);
1624
1625                                                   if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax) {
1626                                                           fNcells_electrons[i]->Fill(ncells);
1627                                                   }
1628
1629                                                   if(M02<0.5 && M20<0.3) {
1630                                                           fEoverP_wSSCut[i]->Fill(EoverP);
1631                                                   }
1632                                           }
1633                                   }
1634
1635                                   fNcells_pt->Fill(fPt, ncells);
1636
1637
1638                                   //!//////////////////////////////////////////////////////////////////
1639                                   ///EMCal - Efficiency calculations
1640
1641                                   /// changing start here
1642                                   if(fIsMC && fIsAOD && track->GetLabel()>=0)
1643                                   {
1644                                           fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1645                                           Int_t pdg = fMCparticle->GetPdgCode();
1646
1647                                           //
1648                                           if(fMCparticle->IsPhysicalPrimary()){
1649
1650
1651                                                   if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
1652
1653                                                           Bool_t MotherFound = FindMother(track->GetLabel());
1654                                                           if(MotherFound){
1655                                                                   fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1656                                                                   if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
1657                                                                           fPtMC_EMCal_All->Fill(fMCparticle->Pt());
1658                                                                   }
1659                                                           }
1660                                                   }
1661                                           }
1662                                   }///until here
1663
1664                                   else if(fIsMC && track->GetLabel()>=0)//ESD
1665                                   {
1666
1667                                           if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
1668                                           {
1669
1670                                                   fMCtrack = fMCstack->Particle(track->GetLabel());
1671
1672                                                   Int_t pdg = fMCtrack->GetPdgCode();
1673                                                   if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax && fMCtrack->Phi()>=(TMath::Pi()*80/180) && fMCtrack->Phi()<=TMath::Pi())
1674                                                   {
1675                                                           if(fMCtrack->GetFirstMother()>0){
1676                                                                   fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1677                                                                   if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
1678
1679                                                                           fPtMC_EMCal_All->Fill(fMCtrack->Pt());
1680                                                                   }
1681                                                           }
1682                                                   }
1683                                           }
1684                                   }
1685
1686
1687                               if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax)
1688                               {
1689
1690                                             fECluster[2]->Fill(Energy);
1691                                                 //_______________________________________________________
1692                                                 //Correlation Analysis
1693                                                 if(fUseEMCal)
1694                                                 {
1695                                                         fPtElec_Inc->Fill(fPt);
1696
1697                                                         if(fCorrelationFlag)
1698                                                         {
1699                                                                 ElectronHadronCorrelation(track, iTracks, Vtrack);
1700                                                         }
1701                                                 }
1702                                                 //_______________________________________________________
1703
1704                                           ////////////////////////////////////////////////////////////////////
1705                                           ///EMCal - efficiency calculations
1706
1707                                           /// changing start here
1708                                           if(fIsMC && fIsAOD && track->GetLabel()>=0)
1709                                           {
1710                                                   fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1711                                                   Int_t pdg = fMCparticle->GetPdgCode();
1712
1713                                                   //
1714                                                   if(fMCparticle->IsPhysicalPrimary()){
1715
1716
1717                                                           if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
1718                                                                   //check
1719                                                                   if(fIsHFE1) fPtMCelectronAfterAll->Fill(fMCparticle->Pt()); //numerator for the total efficiency
1720
1721
1722                                                                   Bool_t MotherFound = FindMother(track->GetLabel());
1723                                                                   if(MotherFound){
1724                                                                           fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1725                                                                           if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
1726                                                                                   fPtMC_EMCal_Selected->Fill(fMCparticle->Pt());
1727                                                                           }
1728                                                                   }
1729                                                           }
1730                                                   }
1731                                           }///until here
1732
1733                                           else if(fIsMC && track->GetLabel()>=0)//ESD
1734                                           {
1735                                                   if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
1736                                                   {
1737                                                                 Bool_t MotherFound = FindMother(track->GetLabel());
1738                                                             fMCtrack = fMCstack->Particle(track->GetLabel());
1739                                                                 if(MotherFound)
1740                                                                 {
1741                                                                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax){
1742                                                                                 if(fIsHFE1) fPtMCelectronAfterAll->Fill(fMCtrack->Pt()); //numerator for the total efficiency
1743                                                                         }
1744                                                                 }
1745
1746
1747
1748                                                       Int_t pdg = fMCtrack->GetPdgCode();
1749                                                       if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax && fMCtrack->Phi()>=(TMath::Pi()*80/180) && fMCtrack->Phi()<=TMath::Pi())
1750                                                       {
1751                                                               if(fMCtrack->GetFirstMother()>0){
1752                                                                    fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1753                                                                    if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
1754
1755                                                                                   fPtMC_EMCal_Selected->Fill(fMCtrack->Pt());
1756                                                                    }
1757                                                                   }
1758                                                       }
1759                                                   }
1760                                           }
1761                                           ///////////////////////////////////////////////////////////////////
1762
1763
1764                                   }
1765                           }
1766                         }
1767                 }
1768
1769                 //______________________________________________________________
1770                 // Vertex
1771
1772                 fVtxZ[2]->Fill(fZvtx);
1773                 fNTracks[2]->Fill(fNOtrks);
1774                 fNClusters[2]->Fill(ClsNo);
1775
1776                 //______________________________________________________________
1777
1778                 //_______________________________________________________
1779                 //Correlation Analysis
1780                 if(!fUseEMCal)
1781                 {
1782                         fPtElec_Inc->Fill(fPt);
1783
1784                         if(fCorrelationFlag)
1785                         {
1786                                 ElectronHadronCorrelation(track, iTracks, Vtrack);
1787                         }
1788                 }
1789                 //_______________________________________________________
1790
1791 ///________________________________________________________________________
1792         }
1793
1794         //__________________________________________________________________
1795         //Event Mixing Analysis
1796         //Filling pool
1797         if(fEventMixingFlag)
1798         {
1799                 fPool = fPoolMgr->GetEventPool(fCentrality->GetCentralityPercentile("V0A"), fZvtx); // Get the buffer associated with the current centrality and z-vtx
1800
1801                 if(!fPool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality->GetCentralityPercentile("V0A"), fZvtx));
1802
1803                 fPool->UpdatePool(SelectedHadrons()); // fill the tracks in the event buffer. The ownership is handed over to the event mixing class. We are not allowed to delete tracksClone anymore!
1804
1805
1806         }
1807
1808         //__________________________________________________________________
1809
1810         delete fListOfmotherkink;
1811         PostData(1, fOutputList);
1812 }
1813
1814 //______________________________________________________________________
1815 void AliAnalysisTaskEMCalHFEpA::Terminate(Option_t *)
1816 {
1817 //Draw result to the screen
1818 //Called once at the end of the query
1819
1820         fOutputList = dynamic_cast<TList*> (GetOutputData(1));
1821
1822         if(!fOutputList)
1823         {
1824                 printf("ERROR: Output list not available\n");
1825                 return;
1826         }
1827 }
1828
1829 //______________________________________________________________________
1830 Bool_t AliAnalysisTaskEMCalHFEpA::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1831 {
1832 //Check single track cuts for a given cut step
1833 //Note this function is called inside the UserExec function
1834         const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1835         if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1836         return kTRUE;
1837 }
1838
1839 //______________________________________________________________________
1840 void AliAnalysisTaskEMCalHFEpA::ElectronHadronCorrelation(AliVTrack *track, Int_t trackIndex, AliVParticle *vtrack)
1841 {
1842         ///_________________________________________________________________
1843         ///MC analysis
1844         if(fIsMC)
1845         {
1846                  if(track->GetLabel() < 0)
1847         {
1848                 AliWarning(Form("The track %d does not have a valid MC label",trackIndex));
1849                 return;
1850         }
1851
1852                 if(fIsAOD)
1853                 {
1854                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1855
1856                         if(fMCparticle->GetMother()<0) return;
1857
1858                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1859
1860                 if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
1861                 {
1862                                 //Is Background
1863                                 fPtBackgroundBeforeReco->Fill(track->Pt());
1864                         }
1865                 }
1866                 else
1867                 {
1868                 fMCtrack = fMCstack->Particle(track->GetLabel());
1869
1870                 if(fMCtrack->GetFirstMother()<0) return;
1871
1872                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1873
1874                 if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
1875                 {
1876                                 //Is Background
1877                                 fPtBackgroundBeforeReco->Fill(track->Pt());
1878                         }
1879                 }
1880         }
1881         ///_________________________________________________________________
1882
1883         //________________________________________________
1884         //Associated particle cut
1885         fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
1886     fPartnerCuts->SetRequireITSRefit(kTRUE);
1887     fPartnerCuts->SetRequireTPCRefit(kTRUE);
1888     fPartnerCuts->SetEtaRange(-0.9,0.9);
1889     fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
1890     fPartnerCuts->SetMinNClustersTPC(80);
1891     fPartnerCuts->SetPtRange(0.3,1e10);
1892     //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
1893     //fPartnerCuts->SetMaxDCAToVertexXY(1);
1894     //fPartnerCuts->SetMaxDCAToVertexZ(3);
1895         //_________________________________________________
1896
1897         ///#################################################################
1898         //Non-HFE reconstruction
1899         fNonHFE = new AliSelectNonHFE();
1900         fNonHFE->SetAODanalysis(fIsAOD);
1901         if(fMassCutFlag) fNonHFE->SetInvariantMassCut(fMassCut);
1902         if(fAngleCutFlag) fNonHFE->SetOpeningAngleCut(fAngleCut);
1903         if(fChi2CutFlag) fNonHFE->SetChi2OverNDFCut(fChi2Cut);
1904         if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
1905         fNonHFE->SetAlgorithm("DCA"); //KF
1906         fNonHFE->SetPIDresponse(fPidResponse);
1907         fNonHFE->SetTrackCuts(-3.5,3.5,fPartnerCuts);
1908
1909         fNonHFE->SetHistAngleBack(fOpAngleBack);
1910         fNonHFE->SetHistAngle(fOpAngle);
1911         fNonHFE->SetHistDCABack(fDCABack);
1912         fNonHFE->SetHistDCA(fDCA);
1913         fNonHFE->SetHistMassBack(fInvMassBack);
1914         fNonHFE->SetHistMass(fInvMass);
1915
1916         fNonHFE->FindNonHFE(trackIndex,vtrack,fVevent);
1917
1918         Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
1919         Int_t *fLsPartner = fNonHFE->GetPartnersLS();
1920         Bool_t fUlsIsPartner = kFALSE;
1921         Bool_t fLsIsPartner = kFALSE;
1922         ///#################################################################
1923
1924         //Electron Information
1925         Double_t fPhiE = -999;
1926         Double_t fEtaE = -999;
1927         Double_t fPhiH = -999;
1928         Double_t fEtaH = -999;
1929         Double_t fDphi = -999;
1930         Double_t fDeta = -999;
1931         Double_t fPtE = -999;
1932         Double_t fPtH = -999;
1933
1934         Double_t pi = TMath::Pi();
1935
1936         fPhiE = track->Phi();
1937         fEtaE = track->Eta();
1938         fPtE = track->Pt();
1939
1940         ///_________________________________________________________________
1941         ///MC analysis
1942         if(fIsMC)
1943         {
1944                 if(fIsAOD)
1945                 {
1946                         if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
1947                 {
1948                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
1949                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
1950                         }
1951                 }
1952                 else
1953                 {
1954                         if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
1955                         {
1956                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
1957                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
1958                         }
1959                 }
1960         }
1961         ///_________________________________________________________________
1962         else
1963         {
1964                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
1965                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
1966         }
1967
1968         //__________________________________________________________________
1969         //Event Mixing Analysis - Hadron Loop
1970         //Retrieve
1971         if(fEventMixingFlag)
1972         {
1973                 fPool = fPoolMgr->GetEventPool(fCentrality->GetCentralityPercentile("V0A"), fZvtx); // Get the buffer associated with the current centrality and z-vtx
1974
1975                 if(!fPool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f",fCentrality->GetCentralityPercentile("V0A"), fZvtx));
1976
1977                 if(fPool->GetCurrentNEvents() >= 5) // start mixing when 5 events are in the buffer
1978                 {
1979                         fPoolNevents->Fill(fPool->GetCurrentNEvents());
1980
1981                         for (Int_t jMix = 0; jMix < fPool->GetCurrentNEvents(); jMix++)  // mix with each event in the buffer
1982                         {
1983                                 TObjArray* bgTracks = fPool->GetEvent(jMix);
1984
1985                                 for (Int_t kMix = 0; kMix < bgTracks->GetEntriesFast(); kMix++)  // mix with each track in the event
1986                                 {
1987                                         const AliEHCParticle* MixedTrack(dynamic_cast<AliEHCParticle*>(bgTracks->At(kMix)));
1988                                         if (NULL == MixedTrack) continue;
1989
1990                                         fPhiH = MixedTrack->Phi();
1991                                         fEtaH = MixedTrack->Eta();
1992                                         fPtH = MixedTrack->Pt();
1993
1994                                         if(fPtH>fPtE || fPtH<1) continue;
1995
1996                                         fDphi = fPhiE - fPhiH;
1997
1998                                         if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
1999                                         if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
2000
2001                                         fDeta = fEtaE - fEtaH;
2002
2003                                         Double_t fPtBin[6] = {2,4,6,8,10,15};
2004
2005                                         for(Int_t i = 0; i < 5; i++)
2006                                         {
2007                                             if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
2008                                             {
2009                                                         fCEtaPhi_Inc_EM[i]->Fill(fDphi,fDeta);
2010
2011                                                         if(fNonHFE->IsULS()) fCEtaPhi_ULS_EM[i]->Fill(fDphi,fDeta);
2012                                                         if(fNonHFE->IsLS()) fCEtaPhi_LS_EM[i]->Fill(fDphi,fDeta);
2013
2014                                                         if(fNonHFE->IsULS()) fCEtaPhi_ULS_Weight_EM[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
2015                                                         if(fNonHFE->IsLS()) fCEtaPhi_LS_Weight_EM[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
2016                                             }
2017                                         }
2018
2019                                         // TODO your code: do event mixing with current event and bgTracks
2020                                         // note that usually the content filled now is weighted by 1 / pool->GetCurrentNEvents()
2021                                 }
2022                         }
2023                 }
2024         }
2025         //__________________________________________________________________
2026
2027         //__________________________________________________________________
2028         //Same Event Analysis - Hadron Loop
2029         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
2030         {
2031                 if(trackIndex==iTracks) continue;
2032
2033                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
2034                 if (!Vtrack2)
2035                 {
2036                         printf("ERROR: Could not receive track %d\n", iTracks);
2037                         continue;
2038                 }
2039
2040                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
2041
2042                  if(fIsAOD)
2043                  {
2044                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
2045                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
2046                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
2047                         if(atrack2->GetTPCNcls() < 80) continue;
2048                 }
2049                 else
2050                 {
2051                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
2052                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
2053                 }
2054
2055                 fPhiH = track2->Phi();
2056                 fEtaH = track2->Eta();
2057                 fPtH = track2->Pt();
2058
2059                 if(fPtH>fPtE || fPtH<1) continue;
2060
2061                 fDphi = fPhiE - fPhiH;
2062
2063                 if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
2064                 if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
2065
2066                 fDeta = fEtaE - fEtaH;
2067
2068                 Double_t fPtBin[6] = {2,4,6,8,10,15};
2069
2070                 //______________________________________________________________
2071                 //Check if this track is a Non-HFE partner
2072                         for(Int_t i = 0; i < fNonHFE->GetNULS(); i++)
2073                         {
2074                                 if(fUlsPartner[i]==iTracks) fUlsIsPartner=kTRUE;
2075                         }
2076                         for(Int_t i = 0; i < fNonHFE->GetNLS(); i++)
2077                         {
2078                                 if(fLsPartner[i]==iTracks) fLsIsPartner=kTRUE;
2079                         }
2080                 //______________________________________________________________
2081
2082                   for(Int_t i = 0; i < 5; i++)
2083                   {
2084                     if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
2085                     {
2086                                 fCEtaPhi_Inc[i]->Fill(fDphi,fDeta);
2087
2088                                 if(fNonHFE->IsULS()) fCEtaPhi_ULS[i]->Fill(fDphi,fDeta);
2089                                 if(fNonHFE->IsLS()) fCEtaPhi_LS[i]->Fill(fDphi,fDeta);
2090                                 if(fNonHFE->IsULS() && !fUlsIsPartner) fCEtaPhi_ULS_NoP[i]->Fill(fDphi,fDeta);
2091                                 if(fNonHFE->IsLS() && !fLsIsPartner) fCEtaPhi_LS_NoP[i]->Fill(fDphi,fDeta);
2092
2093                                 if(fNonHFE->IsULS()) fCEtaPhi_ULS_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
2094                                 if(fNonHFE->IsLS()) fCEtaPhi_LS_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
2095                                 if(fNonHFE->IsULS() && !fUlsIsPartner) fCEtaPhi_ULS_NoP_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
2096                                 if(fNonHFE->IsLS() && !fLsIsPartner) fCEtaPhi_LS_NoP_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
2097                     }
2098                   }
2099         }
2100 }
2101
2102 //____________________________________________________________________________________________________________
2103 //Create a TObjArray with selected hadrons, for the mixed event analysis
2104 TObjArray* AliAnalysisTaskEMCalHFEpA::SelectedHadrons()
2105 {
2106         fTracksClone = new TObjArray;
2107         fTracksClone->SetOwner(kTRUE);
2108
2109         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
2110         {
2111                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
2112                 if (!Vtrack2)
2113                 {
2114                         printf("ERROR: Could not receive track %d\n", iTracks);
2115                         continue;
2116                 }
2117
2118                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
2119
2120                  if(fIsAOD)
2121                  {
2122                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
2123                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
2124                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
2125                         if(atrack2->GetTPCNcls() < 80) continue;
2126                 }
2127                 else
2128                 {
2129                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
2130                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
2131                 }
2132
2133                 fTracksClone->Add(new AliEHCParticle(track2->Eta(), track2->Phi(), track2->Pt()));
2134         }
2135         return fTracksClone;
2136 }
2137 //____________________________________________________________________________________________________________
2138
2139 //______________________________________________________________________
2140 void AliAnalysisTaskEMCalHFEpA::DiHadronCorrelation(AliVTrack *track, Int_t trackIndex)
2141 {
2142         //________________________________________________
2143         //Associated particle cut
2144         fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
2145     fPartnerCuts->SetRequireITSRefit(kTRUE);
2146     fPartnerCuts->SetRequireTPCRefit(kTRUE);
2147     fPartnerCuts->SetEtaRange(-0.9,0.9);
2148     fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
2149     fPartnerCuts->SetMinNClustersTPC(80);
2150     fPartnerCuts->SetPtRange(0.3,1e10);
2151     //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
2152     //fPartnerCuts->SetMaxDCAToVertexXY(1);
2153     //fPartnerCuts->SetMaxDCAToVertexZ(3);
2154         //_________________________________________________
2155
2156         //Electron Information
2157         Double_t fPhiE = -999;
2158         Double_t fEtaE = -999;
2159         Double_t fPhiH = -999;
2160         Double_t fEtaH = -999;
2161         Double_t fDphi = -999;
2162         Double_t fDeta = -999;
2163         Double_t fPtE = -999;
2164         Double_t fPtH = -999;
2165
2166         Double_t pi = TMath::Pi();
2167
2168         fPhiE = track->Phi();
2169         fEtaE = track->Eta();
2170         fPtE = track->Pt();
2171
2172         //__________________________________________________________________
2173         //Same Event Analysis - Hadron Loop
2174         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
2175         {
2176                 if(trackIndex==iTracks) continue;
2177
2178                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
2179                 if (!Vtrack2)
2180                 {
2181                         printf("ERROR: Could not receive track %d\n", iTracks);
2182                         continue;
2183                 }
2184
2185                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
2186
2187                  if(fIsAOD)
2188                  {
2189                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
2190                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
2191                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
2192                         if(atrack2->GetTPCNcls() < 80) continue;
2193                 }
2194                 else
2195                 {
2196                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
2197                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
2198                 }
2199
2200                 fPhiH = track2->Phi();
2201                 fEtaH = track2->Eta();
2202                 fPtH = track2->Pt();
2203
2204                 if(fPtH>fPtE || fPtH<1) continue;
2205
2206                 fDphi = fPhiE - fPhiH;
2207
2208                 if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
2209                 if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
2210
2211                 fDeta = fEtaE - fEtaH;
2212
2213                 Double_t fPtBin[6] = {2,4,6,8,10,15};
2214
2215                   for(Int_t i = 0; i < 5; i++)
2216                   {
2217                     if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
2218                     {
2219                                 fCEtaPhi_Inc_DiHadron[i]->Fill(fDphi,fDeta);
2220                     }
2221                   }
2222         }
2223 }
2224 //____________________________________________________________________________________________________________
2225
2226 //______________________________________________________________________
2227 Bool_t AliAnalysisTaskEMCalHFEpA::FindMother(Int_t mcIndex)
2228 {
2229         fIsHFE1 = kFALSE;
2230         fIsHFE2 = kFALSE;
2231         fIsNonHFE = kFALSE;
2232         fIsFromD = kFALSE;
2233         fIsFromB = kFALSE;
2234         fIsFromPi0 = kFALSE;
2235         fIsFromEta = kFALSE;
2236         fIsFromGamma = kFALSE;
2237
2238         if(mcIndex < 0 || !fIsMC)
2239         {
2240                 return kFALSE;
2241         }
2242
2243         Int_t pdg = -99999;
2244         Int_t mpdg = -99999;
2245         Int_t gmpdg = -99999;
2246         Int_t ggmpdg = -99999;
2247         Int_t gggmpdg = -99999;
2248
2249         if(fIsAOD)
2250         {
2251                 fMCparticle = (AliAODMCParticle*) fMCarray->At(mcIndex);
2252
2253                 pdg = TMath::Abs(fMCparticle->GetPdgCode());
2254
2255
2256                 if(pdg!=11)
2257                 {
2258                         fIsHFE1 = kFALSE;
2259                         fIsHFE2 = kFALSE;
2260                         fIsNonHFE = kFALSE;
2261                         fIsFromD = kFALSE;
2262                         fIsFromB = kFALSE;
2263                         fIsFromPi0 = kFALSE;
2264                         fIsFromEta = kFALSE;
2265                         fIsFromGamma = kFALSE;
2266                         return kFALSE;
2267                 }
2268
2269                 if(fMCparticle->GetMother()<0)
2270                 {
2271                         fIsHFE1 = kFALSE;
2272                         fIsHFE2 = kFALSE;
2273                         fIsNonHFE = kFALSE;
2274                         fIsFromD = kFALSE;
2275                         fIsFromB = kFALSE;
2276                         fIsFromPi0 = kFALSE;
2277                         fIsFromEta = kFALSE;
2278                         fIsFromGamma = kFALSE;
2279                         return kFALSE;
2280                 }
2281
2282                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2283                 mpdg = TMath::Abs(fMCparticleMother->GetPdgCode());
2284
2285                 if(fMCparticleMother->GetMother()<0)
2286                 {
2287                         gmpdg = 0;
2288                         ggmpdg = 0;
2289                         gggmpdg = 0;
2290                 }
2291                 else
2292                 {
2293                         fMCparticleGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleMother->GetMother());
2294                         gmpdg = TMath::Abs(fMCparticleGMother->GetPdgCode());
2295                         if(fMCparticleGMother->GetMother()<0)
2296                         {
2297                                 ggmpdg = 0;
2298                                 gggmpdg = 0;
2299                         }
2300                         else
2301                         {
2302                                 fMCparticleGGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleGMother->GetMother());
2303                                 ggmpdg = TMath::Abs(fMCparticleGGMother->GetPdgCode());
2304                                 if(fMCparticleGGMother->GetMother()<0)
2305                                 {
2306                                         gggmpdg = 0;
2307                                 }
2308                                 else
2309                                 {
2310                                         fMCparticleGGGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleGGMother->GetMother());
2311                                         gggmpdg = TMath::Abs(fMCparticleGGGMother->GetPdgCode());
2312                                 }
2313                         }
2314                 }
2315         }
2316         else
2317         {
2318                 fMCtrack = fMCstack->Particle(mcIndex);
2319
2320                 pdg = TMath::Abs(fMCtrack->GetPdgCode());
2321
2322                 if(pdg!=11)
2323                 {
2324                         fIsHFE1 = kFALSE;
2325                         fIsHFE2 = kFALSE;
2326                         fIsNonHFE = kFALSE;
2327                         fIsFromD = kFALSE;
2328                         fIsFromB = kFALSE;
2329                         fIsFromPi0 = kFALSE;
2330                         fIsFromEta = kFALSE;
2331                         fIsFromGamma = kFALSE;
2332                         return kFALSE;
2333                 }
2334
2335                 if(fMCtrack->GetFirstMother()<0)
2336                 {
2337                         fIsHFE1 = kFALSE;
2338                         fIsHFE2 = kFALSE;
2339                         fIsNonHFE = kFALSE;
2340                         fIsFromD = kFALSE;
2341                         fIsFromB = kFALSE;
2342                         fIsFromPi0 = kFALSE;
2343                         fIsFromEta = kFALSE;
2344                         fIsFromGamma = kFALSE;
2345                         return kFALSE;
2346                 }
2347
2348                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
2349                 mpdg = TMath::Abs(fMCtrackMother->GetPdgCode());
2350
2351                 if(fMCtrackMother->GetFirstMother()<0)
2352                 {
2353                         gmpdg = 0;
2354                         ggmpdg = 0;
2355                         gggmpdg = 0;
2356                 }
2357                 else
2358                 {
2359                         fMCtrackGMother = fMCstack->Particle(fMCtrackMother->GetFirstMother());
2360                         gmpdg = TMath::Abs(fMCtrackGMother->GetPdgCode());
2361
2362                         if(fMCtrackGMother->GetFirstMother()<0)
2363                         {
2364                                 ggmpdg = 0;
2365                                 gggmpdg = 0;
2366                         }
2367                         else
2368                         {
2369                                 fMCtrackGGMother = fMCstack->Particle(fMCtrackGMother->GetFirstMother());
2370                                 ggmpdg = TMath::Abs(fMCtrackGGMother->GetPdgCode());
2371
2372                                 if(fMCtrackGGMother->GetFirstMother()<0)
2373                                 {
2374                                         gggmpdg = 0;
2375                                 }
2376                                 else
2377                                 {
2378                                         fMCtrackGGGMother = fMCstack->Particle(fMCtrackGGMother->GetFirstMother());
2379                                         gggmpdg = TMath::Abs(fMCtrackGGGMother->GetPdgCode());
2380                                 }
2381                         }
2382                 }
2383         }
2384
2385         //Tag Electron Source
2386         if(mpdg==111 || mpdg==221 || mpdg==22)
2387         {
2388                 fIsHFE1 = kFALSE;
2389                 fIsHFE2 = kFALSE;
2390                 fIsNonHFE = kTRUE;
2391                 fIsFromD = kFALSE;
2392                 fIsFromB = kFALSE;
2393
2394                 fIsFromPi0 = kFALSE;
2395                 fIsFromEta = kFALSE;
2396                 fIsFromGamma = kFALSE;
2397
2398                 if(mpdg==111) fIsFromPi0 = kFALSE;
2399                 if(mpdg==221)fIsFromEta = kFALSE;
2400                 if(mpdg==22) fIsFromGamma = kFALSE;
2401
2402                 return kTRUE;
2403         }
2404         else
2405         {
2406                 fIsHFE1 = kFALSE;
2407                 fIsHFE2 = kTRUE;
2408
2409                 fIsFromPi0 = kFALSE;
2410                 fIsFromEta = kFALSE;
2411                 fIsFromGamma = kFALSE;
2412
2413                 fIsNonHFE = kFALSE;
2414
2415                 fIsFromD = kFALSE;
2416                 fIsFromB = kFALSE;
2417
2418                 if(mpdg>400 && mpdg<500)
2419                 {
2420                         if((gmpdg>500 && gmpdg<600) || (ggmpdg>500 && ggmpdg<600) || (gggmpdg>500 && gggmpdg<600))
2421                         {
2422                                 fIsHFE1 = kTRUE;
2423                                 fIsFromD = kFALSE;
2424                                 fIsFromB = kTRUE;
2425                                 return kTRUE;
2426                         }
2427                         else
2428                         {
2429                                 fIsHFE1 = kTRUE;
2430                                 fIsFromD = kTRUE;
2431                                 fIsFromB = kFALSE;
2432                                 return kTRUE;
2433                         }
2434                 }
2435                 else if(mpdg>500 && mpdg<600)
2436                 {
2437                         fIsHFE1 = kTRUE;
2438                         fIsFromD = kFALSE;
2439                         fIsFromB = kTRUE;
2440                         return kTRUE;
2441                 }
2442                 else
2443                 {
2444                         fIsHFE1 = kFALSE;
2445                         fIsFromD = kFALSE;
2446                         fIsFromB = kFALSE;
2447                         return kFALSE;
2448                 }
2449         }
2450 }