PWGHF/hfe/AliAnalysisTaskEMCalHFEpA.cxx

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