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