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         //              version: September 12th, 2013.                                                            //
  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 ,fUseShowerShapeCut(kFALSE)
 103 ,fFillBackground(kFALSE)
 104 ,fAssocWithSPD(kFALSE)
 105
 106 ,fEMCEG1(kFALSE)
 107 ,fEMCEG2(kFALSE)
 108 ,fIsHFE1(kFALSE)
 109 ,fIsHFE2(kFALSE)
 110 ,fIsNonHFE(kFALSE)
 111 ,fIsFromD(kFALSE)
 112 ,fIsFromB(kFALSE)
 113 ,fIsFromPi0(kFALSE)
 114 ,fIsFromEta(kFALSE)
 115 ,fIsFromGamma(kFALSE)
 116 ,fESD(0)
 117 ,fAOD(0)
 118 ,fVevent(0)
 119 ,fPartnerCuts(new AliESDtrackCuts())
 120 ,fOutputList(0)
 121 ,fPidResponse(0)
 122 ,fNonHFE(new AliSelectNonHFE())
 123 ,fIsAOD(kFALSE)
 124 ,fCentrality(0)
 125 ,fCentralityMin(0)
 126 ,fCentralityMax(100)
 127 ,fHasCentralitySelection(kFALSE)
 128 ,fCentralityHist(0)
 129 ,fCentralityHistPass(0)
 130 ,fZvtx(0)
 131 ,fEstimator(0)
 132 ,fClus(0)
 133         //,fClusESD(0)
 134 ,fNevent(0)
 135 ,fPtElec_Inc(0)
 136
 137 ,fCharge_n(0)
 138 ,fCharge_p(0)
 139
 140 ,fTime(0)
 141 ,fTime2(0)
 142 ,ftimingEle(0)
 143 ,ftimingEle2(0)
 144
 145
 146 ,fPtElec_ULS(0)
 147 ,fPtElec_LS(0)
 148 ,fPtElec_ULS2(0)
 149 ,fPtElec_LS2(0)
 150
 151 ,fPtElec_ULS_weight(0)
 152 ,fPtElec_LS_weight(0)
 153 ,fPtElec_ULS2_weight(0)
 154 ,fPtElec_LS2_weight(0)
 155
 156 ,fTOF01(0)
 157 ,fTOF02(0)
 158 ,fTOF03(0)
 159 ,fpid(0)
 160 ,fEoverP_pt(0)
 161 ,fEoverP_tpc(0)
 162 ,fTPC_pt(0)
 163 ,fTPC_p(0)
 164 ,fTPCnsigma_pt(0)
 165 ,fTPCnsigma_p(0)
 166 ,fTPCnsigma_pt_2D(0)
 167 ,fShowerShapeCut(0)
 168
 169 ,fShowerShapeM02_EoverP(0)
 170 ,fShowerShapeM20_EoverP(0)
 171
 172 ,fShowerShape_ha(0)
 173 ,fShowerShape_ele(0)
 174
 175 ,fTPCnsigma_eta(0)
 176 ,fTPCnsigma_phi(0)
 177 ,fECluster(0)
 178 ,fEtaPhi(0)
 179 ,fVtxZ(0)
 180 ,fNTracks(0)
 181 ,fNClusters(0)
 182 ,fTPCNcls_EoverP(0)
 183 ,fEta(0)
 184 ,fPhi(0)
 185 ,fR(0)
 186 ,fR_EoverP(0)
 187 ,fNcells(0)
 188 ,fNcells_EoverP(0)
 189 ,fNcells_electrons(0)
 190 ,fNcells_hadrons(0)
 191 ,fECluster_ptbins(0)
 192 ,fEoverP_ptbins(0)
 193 ,fEoverP_wSSCut(0)
 194 ,fM02_EoverP(0)
 195 ,fM20_EoverP(0)
 196 ,fTPCnsigma_eta_electrons(0)
 197 ,fTPCnsigma_eta_hadrons(0)
 198 ,fEoverP_pt_pions(0)
 199 ,ftpc_p_EoverPcut(0)
 200 ,fnsigma_p_EoverPcut(0)
 201 ,fEoverP_pt_pions2(0)
 202 ,fNcells_pt(0)
 203 ,fEoverP_pt_hadrons(0)
 204 ,fCEtaPhi_Inc(0)
 205 ,fCEtaPhi_ULS(0)
 206 ,fCEtaPhi_LS(0)
 207 ,fCEtaPhi_ULS_NoP(0)
 208 ,fCEtaPhi_LS_NoP(0)
 209 ,fCEtaPhi_ULS_Weight(0)
 210 ,fCEtaPhi_LS_Weight(0)
 211 ,fCEtaPhi_ULS_NoP_Weight(0)
 212 ,fCEtaPhi_LS_NoP_Weight(0)
 213
 214 ,fInvMass(0)
 215 ,fInvMassBack(0)
 216 ,fDCA(0)
 217 ,fDCABack(0)
 218 ,fOpAngle(0)
 219 ,fOpAngleBack(0)
 220
 221 ,fInvMass2(0)
 222 ,fInvMassBack2(0)
 223 ,fDCA2(0)
 224 ,fDCABack2(0)
 225 ,fOpAngle2(0)
 226 ,fOpAngleBack2(0)
 227
 228 ,fMassCut(0.1)
 229 ,fEtaCutMin(-0.9)
 230 ,fEtaCutMax(0.9)
 231
 232 ,fdPhiCut(0.05)
 233 ,fdEtaCut(0.05)
 234
 235
 236 ,fEoverPCutMin(0.8)
 237 ,fEoverPCutMax(1.2)
 238
 239 ,fM20CutMin(0.0)
 240 ,fM20CutMax(10)
 241 ,fM02CutMin(0.0)
 242 ,fM02CutMax(10)
 243
 244 ,fAngleCut(999)
 245 ,fChi2Cut(3.5)
 246 ,fDCAcut(999)
 247 ,fMassCutFlag(kTRUE)
 248 ,fAngleCutFlag(kFALSE)
 249 ,fChi2CutFlag(kFALSE)
 250 ,fDCAcutFlag(kFALSE)
 251 ,fPtBackgroundBeforeReco(0)
 252 ,fPtBackgroundBeforeReco2(0)
 253 ,fPtBackgroundBeforeReco_weight(0)
 254 ,fPtBackgroundBeforeReco2_weight(0)
 255 ,fPtBackgroundAfterReco(0)
 256
 257 ,fPtMinAsso(0.3)
 258 ,fTpcNclsAsso(80)
 259
 260 ,fPtMCparticleAll(0)
 261 ,fPtMCparticleAll_nonPrimary(0)
 262 ,fPtMCparticleAlle_nonPrimary(0)
 263 ,fPtMCparticleAlle_Primary(0)
 264 ,fPtMCparticleReco(0)
 265 ,fPtMCparticleReco_nonPrimary(0)
 266 ,fPtMCparticleAllHfe1(0)
 267 ,fPtMCparticleRecoHfe1(0)
 268 ,fPtMCparticleAllHfe2(0)
 269 ,fPtMCparticleRecoHfe2(0)
 270 ,fPtMCelectronAfterAll(0)
 271 ,fPtMCelectronAfterAll_nonPrimary(0)
 272 ,fPtMCelectronAfterAll_Primary(0)
 273 ,fPtMCpi0(0)
 274 ,fPtMCeta(0)
 275 ,fPtMC_EMCal_All(0)
 276 ,fPtMC_EMCal_Selected(0)
 277 ,fPtMC_TPC_All(0)
 278 ,fPtMC_TPC_Selected(0)
 279 ,fPtMCWithLabel(0)
 280 ,fPtMCWithoutLabel(0)
 281 ,fPtIsPhysicaPrimary(0)
 282 ,fCuts(0)
 283 ,fCFM(0)
 284 ,fPID(new AliHFEpid("hfePid"))
 285 ,fPIDqa(0)
 286 ,fMCstack(0)
 287 ,fRejectKinkMother(kFALSE)
 288 ,fMCtrack(0)
 289 ,fMCtrackMother(0)
 290 ,fMCtrackGMother(0)
 291 ,fMCtrackGGMother(0)
 292 ,fMCtrackGGGMother(0)
 293 ,fMCarray(0)
 294 ,fMCheader(0)
 295 ,fMCparticle(0)
 296 ,fMCparticleMother(0)
 297 ,fMCparticleGMother(0)
 298 ,fMCparticleGGMother(0)
 299 ,fMCparticleGGGMother(0)
 300 ,fEventHandler(0)
 301 ,fMCevent(0)
 302 ,fPoolMgr(0)
 303 ,fPool(0)
 304 ,fTracksClone(0)
 305 ,fTracks(0)
 306 ,fCEtaPhi_Inc_EM(0)
 307 ,fCEtaPhi_ULS_EM(0)
 308 ,fCEtaPhi_LS_EM(0)
 309 ,fCEtaPhi_ULS_Weight_EM(0)
 310 ,fCEtaPhi_LS_Weight_EM(0)
 311 ,fPoolNevents(0)
 312 ,fEventMixingFlag(0)
 313 ,fCEtaPhi_Inc_DiHadron(0)
 314 ,fPtTrigger_Inc(0)
 315 {
 316                 //Named constructor
 317                 // Define input and output slots here
 318                 // Input slot #0 works with a TChain
 319         DefineInput(0, TChain::Class());
 320                 // Output slot #0 id reserved by the base class for AOD
 321                 // Output slot #1 writes into a TH1 container
 322                 // DefineOutput(1, TH1I::Class());
 323         DefineOutput(1, TList::Class());
 324                 //  DefineOutput(3, TTree::Class());
 325 }
 326
 327         //________________________________________________________________________
 328 AliAnalysisTaskEMCalHFEpA::AliAnalysisTaskEMCalHFEpA()
 329 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCalHFEpA")
 330 ,fCorrelationFlag(0)
 331 ,fIsMC(0)
 332 ,fUseEMCal(kFALSE)
 333 ,fUseShowerShapeCut(kFALSE)
 334 ,fFillBackground(kFALSE)
 335 ,fAssocWithSPD(kFALSE)
 336 ,fEMCEG1(kFALSE)
 337 ,fEMCEG2(kFALSE)
 338 ,fIsHFE1(kFALSE)
 339 ,fIsHFE2(kFALSE)
 340 ,fIsNonHFE(kFALSE)
 341 ,fIsFromD(kFALSE)
 342 ,fIsFromB(kFALSE)
 343 ,fIsFromPi0(kFALSE)
 344 ,fIsFromEta(kFALSE)
 345 ,fIsFromGamma(kFALSE)
 346 ,fESD(0)
 347 ,fAOD(0)
 348 ,fVevent(0)
 349 ,fPartnerCuts(new AliESDtrackCuts())
 350 ,fOutputList(0)
 351 ,fPidResponse(0)
 352 ,fNonHFE(new AliSelectNonHFE())
 353 ,fIsAOD(kFALSE)
 354 ,fCentrality(0)
 355 ,fCentralityMin(0)
 356 ,fCentralityMax(100)
 357 ,fHasCentralitySelection(kFALSE)
 358 ,fCentralityHist(0)
 359 ,fCentralityHistPass(0)
 360 ,fZvtx(0)
 361 ,fEstimator(0)
 362 ,fClus(0)
 363         //,fClusESD(0)
 364 ,fNevent(0)
 365 ,fPtElec_Inc(0)
 366
 367 ,fCharge_n(0)
 368 ,fCharge_p(0)
 369
 370 ,fTime(0)
 371 ,fTime2(0)
 372 ,ftimingEle(0)
 373 ,ftimingEle2(0)
 374
 375
 376
 377 ,fPtElec_ULS(0)
 378 ,fPtElec_LS(0)
 379 ,fPtElec_ULS2(0)
 380 ,fPtElec_LS2(0)
 381
 382 ,fPtElec_ULS_weight(0)
 383 ,fPtElec_LS_weight(0)
 384 ,fPtElec_ULS2_weight(0)
 385 ,fPtElec_LS2_weight(0)
 386
 387 ,fTOF01(0)
 388 ,fTOF02(0)
 389 ,fTOF03(0)
 390 ,fpid(0)
 391 ,fEoverP_pt(0)
 392 ,fEoverP_tpc(0)
 393 ,fTPC_pt(0)
 394 ,fTPC_p(0)
 395 ,fTPCnsigma_pt(0)
 396 ,fTPCnsigma_p(0)
 397 ,fTPCnsigma_pt_2D(0)
 398 ,fShowerShapeCut(0)
 399 ,fShowerShapeM02_EoverP(0)
 400 ,fShowerShapeM20_EoverP(0)
 401
 402 ,fShowerShape_ha(0)
 403 ,fShowerShape_ele(0)
 404
 405 ,fTPCnsigma_eta(0)
 406 ,fTPCnsigma_phi(0)
 407 ,fECluster(0)
 408 ,fEtaPhi(0)
 409 ,fVtxZ(0)
 410 ,fNTracks(0)
 411 ,fNClusters(0)
 412 ,fTPCNcls_EoverP(0)
 413 ,fEta(0)
 414 ,fPhi(0)
 415 ,fR(0)
 416 ,fR_EoverP(0)
 417 ,fNcells(0)
 418 ,fNcells_EoverP(0)
 419 ,fNcells_electrons(0)
 420 ,fNcells_hadrons(0)
 421 ,fECluster_ptbins(0)
 422 ,fEoverP_ptbins(0)
 423 ,fEoverP_wSSCut(0)
 424 ,fM02_EoverP(0)
 425 ,fM20_EoverP(0)
 426 ,fTPCnsigma_eta_electrons(0)
 427 ,fTPCnsigma_eta_hadrons(0)
 428 ,fEoverP_pt_pions(0)
 429 ,ftpc_p_EoverPcut(0)
 430 ,fnsigma_p_EoverPcut(0)
 431 ,fEoverP_pt_pions2(0)
 432 ,fNcells_pt(0)
 433 ,fEoverP_pt_hadrons(0)
 434 ,fCEtaPhi_Inc(0)
 435 ,fCEtaPhi_ULS(0)
 436 ,fCEtaPhi_LS(0)
 437 ,fCEtaPhi_ULS_NoP(0)
 438 ,fCEtaPhi_LS_NoP(0)
 439 ,fCEtaPhi_ULS_Weight(0)
 440 ,fCEtaPhi_LS_Weight(0)
 441 ,fCEtaPhi_ULS_NoP_Weight(0)
 442 ,fCEtaPhi_LS_NoP_Weight(0)
 443
 444 ,fInvMass(0)
 445 ,fInvMassBack(0)
 446 ,fDCA(0)
 447 ,fDCABack(0)
 448 ,fOpAngle(0)
 449 ,fOpAngleBack(0)
 450
 451 ,fInvMass2(0)
 452 ,fInvMassBack2(0)
 453 ,fDCA2(0)
 454 ,fDCABack2(0)
 455 ,fOpAngle2(0)
 456 ,fOpAngleBack2(0)
 457
 458 ,fMassCut(0.1)
 459 ,fEtaCutMin(-0.9)
 460 ,fEtaCutMax(0.9)
 461
 462 ,fdPhiCut(0.05)
 463 ,fdEtaCut(0.05)
 464
 465 ,fEoverPCutMin(0.8)
 466 ,fEoverPCutMax(1.2)
 467
 468 ,fM20CutMin(0)
 469 ,fM20CutMax(10)
 470 ,fM02CutMin(0)
 471 ,fM02CutMax(10)
 472
 473 ,fAngleCut(999)
 474 ,fChi2Cut(3.5)
 475 ,fDCAcut(999)
 476 ,fMassCutFlag(kTRUE)
 477 ,fAngleCutFlag(kFALSE)
 478 ,fChi2CutFlag(kFALSE)
 479 ,fDCAcutFlag(kFALSE)
 480 ,fPtBackgroundBeforeReco(0)
 481 ,fPtBackgroundBeforeReco2(0)
 482 ,fPtBackgroundBeforeReco_weight(0)
 483 ,fPtBackgroundBeforeReco2_weight(0)
 484 ,fPtBackgroundAfterReco(0)
 485
 486 ,fPtMinAsso(0.3)
 487 ,fTpcNclsAsso(80)
 488
 489
 490 ,fPtMCparticleAll(0)
 491 ,fPtMCparticleAll_nonPrimary(0)
 492 ,fPtMCparticleAlle_nonPrimary(0)
 493 ,fPtMCparticleAlle_Primary(0)
 494 ,fPtMCparticleReco(0)
 495 ,fPtMCparticleReco_nonPrimary(0)
 496
 497 ,fPtMCparticleAllHfe1(0)
 498 ,fPtMCparticleRecoHfe1(0)
 499 ,fPtMCparticleAllHfe2(0)
 500 ,fPtMCparticleRecoHfe2(0)
 501 ,fPtMCelectronAfterAll(0)
 502 ,fPtMCelectronAfterAll_nonPrimary(0)
 503 ,fPtMCelectronAfterAll_Primary(0)
 504
 505 ,fPtMCpi0(0)
 506 ,fPtMCeta(0)
 507 ,fPtMC_EMCal_All(0)
 508 ,fPtMC_EMCal_Selected(0)
 509 ,fPtMC_TPC_All(0)
 510 ,fPtMC_TPC_Selected(0)
 511 ,fPtMCWithLabel(0)
 512 ,fPtMCWithoutLabel(0)
 513 ,fPtIsPhysicaPrimary(0)
 514 ,fCuts(0)
 515 ,fCFM(0)
 516 ,fPID(new AliHFEpid("hfePid"))
 517 ,fPIDqa(0)
 518 ,fMCstack(0)
 519 ,fRejectKinkMother(kFALSE)
 520 ,fMCtrack(0)
 521 ,fMCtrackMother(0)
 522 ,fMCtrackGMother(0)
 523 ,fMCtrackGGMother(0)
 524 ,fMCtrackGGGMother(0)
 525 ,fMCarray(0)
 526 ,fMCheader(0)
 527 ,fMCparticle(0)
 528 ,fMCparticleMother(0)
 529 ,fMCparticleGMother(0)
 530 ,fMCparticleGGMother(0)
 531 ,fMCparticleGGGMother(0)
 532 ,fEventHandler(0)
 533 ,fMCevent(0)
 534 ,fPoolMgr(0)
 535 ,fPool(0)
 536 ,fTracksClone(0)
 537 ,fTracks(0)
 538 ,fCEtaPhi_Inc_EM(0)
 539 ,fCEtaPhi_ULS_EM(0)
 540 ,fCEtaPhi_LS_EM(0)
 541 ,fCEtaPhi_ULS_Weight_EM(0)
 542 ,fCEtaPhi_LS_Weight_EM(0)
 543 ,fPoolNevents(0)
 544 ,fEventMixingFlag(0)
 545 ,fCEtaPhi_Inc_DiHadron(0)
 546 ,fPtTrigger_Inc(0)
 547 {
 548                 // Constructor
 549                 // Define input and output slots here
 550                 // Input slot #0 works with a TChain
 551         DefineInput(0, TChain::Class());
 552                 // Output slot #0 id reserved by the base class for AOD
 553                 // Output slot #1 writes into a TH1 container
 554                 // DefineOutput(1, TH1I::Class());
 555         DefineOutput(1, TList::Class());
 556                 //DefineOutput(3, TTree::Class());
 557 }
 558
 559         //______________________________________________________________________
 560 AliAnalysisTaskEMCalHFEpA::~AliAnalysisTaskEMCalHFEpA()
 561 {
 562                 //Destructor
 563         delete fOutputList;
 564         delete fPID;
 565         delete fCFM;
 566         delete fPIDqa;
 567 }
 568
 569         //______________________________________________________________________
 570         //Create Output Objects
 571         //Here we can define the histograms and others output files
 572         //Called once
 573 void AliAnalysisTaskEMCalHFEpA::UserCreateOutputObjects()
 574 {
 575                 //______________________________________________________________________
 576                 //Initialize PID
 577         if(!fPID->GetNumberOfPIDdetectors())
 578     {
 579                 fPID->AddDetector("TPC", 0);
 580     }
 581
 582         fPID->SortDetectors();
 583
 584         fPIDqa = new AliHFEpidQAmanager();
 585         fPIDqa->Initialize(fPID);
 586                 //______________________________________________________________________
 587
 588                 //______________________________________________________________________
 589                 //Initialize correction Framework and Cuts
 590         fCFM = new AliCFManager;
 591         const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 592         fCFM->SetNStepParticle(kNcutSteps);
 593         for(Int_t istep = 0; istep < kNcutSteps; istep++) fCFM->SetParticleCutsList(istep, NULL);
 594
 595         if(!fCuts)
 596         {
 597                 AliWarning("Cuts not available. Default cuts will be used");
 598                 fCuts = new AliHFEcuts;
 599                 fCuts->CreateStandardCuts();
 600         }
 601
 602         fCuts->Initialize(fCFM);
 603                 //______________________________________________________________________
 604
 605                 ///______________________________________________________________________
 606                 ///Output Tlist
 607                 //Create TList
 608         fOutputList = new TList();
 609         fOutputList->SetOwner();
 610
 611                 //PIDqa
 612         fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 613
 614                 //Store the number of events
 615                 //Define the histo
 616         fNevent = new TH1F("fNevent","Number of Events",5,-0.5,4.5);
 617                 //And then, add to the output list
 618         fOutputList->Add(fNevent);
 619
 620         fpid = new TH1F("fpid","PID flag",5,0,5);
 621         fOutputList->Add(fpid);
 622
 623                 //pt Distribution
 624         fPtElec_Inc = new TH1F("fPtElec_Inc","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 625
 626         fPtElec_ULS = new TH1F("fPtElec_ULS","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 627         fPtElec_LS = new TH1F("fPtElec_LS","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 628
 629         fPtElec_ULS_weight = new TH1F("fPtElec_ULS_weight","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 630         fPtElec_LS_weight = new TH1F("fPtElec_LS_weight","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 631
 632         fTOF01 = new TH2F("fTOF01","",200,-20,20,200,-20,20);
 633         fTOF02 = new TH2F("fTOF02","",200,-20,20,200,-20,20);
 634         fTOF03 = new TH2F("fTOF03","",200,-20,20,200,-20,20);
 635
 636         if(fFillBackground){
 637                 fPtElec_ULS2 = new TH1F("fPtElec_ULS2","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 638                 fPtElec_LS2 = new TH1F("fPtElec_LS2","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 639
 640                 fPtElec_ULS2_weight = new TH1F("fPtElec_ULS2_weight","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 641                 fPtElec_LS2_weight = new TH1F("fPtElec_LS2_weight","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 642
 643         }
 644
 645         fPtTrigger_Inc = new TH1F("fPtTrigger_Inc","pT dist for Hadron Contamination; p_{t} (GeV/c); Count",300,0,30);
 646         fTPCnsigma_pt_2D = new TH2F("fTPCnsigma_pt_2D",";pt (GeV/c);TPC Electron N#sigma",1000,0.3,30,1000,-15,10);
 647         fShowerShapeCut = new TH2F("fShowerShapeCut","Shower Shape;M02;M20",500,0,1.8,500,0,1.8);
 648
 649
 650
 651         fCharge_n = new TH1F("fCharge_n","Inclusive Electrons (Negative Charge); p_{t} (GeV/c); Count",200,0,30);
 652         fCharge_p = new TH1F("fCharge_p","Inclusive Positrons (Positive Charge); p_{t} (GeV/c); Count",200,0,30);
 653
 654         if(fUseEMCal){
 655
 656                 if(!fIsAOD){
 657                         fTime = new TH2D("fTime","Cells Cluster Time; p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 658                         fTime2 = new TH2D("fTime2","Cells Cluster Time;  p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 659                 }
 660
 661                 ftimingEle = new TH2D("ftimingEle","Cluster Time;  p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 662                 ftimingEle2 = new TH2D("ftimingEle2","Cluster Time;  p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 663
 664                 fShowerShape_ha = new TH2F("fShowerShape_ha","Shower Shape hadrons;M02;M20",500,0,1.8,500,0,1.8);
 665                 fShowerShape_ele = new TH2F("fShowerShape_ele","Shower Shape electrons;M02;M20",500,0,1.8,500,0,1.8);
 666
 667                 fShowerShapeM02_EoverP = new TH2F("fShowerShapeM02_EoverP","Shower Shape;M02;E/p",500,0,1.8,500,0,1.8);
 668                 fShowerShapeM20_EoverP = new TH2F("fShowerShapeM20_EoverP","Shower Shape;M20;E/p",500,0,1.8,500,0,1.8);
 669
 670         }
 671
 672         fOutputList->Add(fTOF01);
 673         fOutputList->Add(fTOF02);
 674         fOutputList->Add(fTOF03);
 675
 676         fOutputList->Add(fPtElec_Inc);
 677         fOutputList->Add(fPtElec_ULS);
 678         fOutputList->Add(fPtElec_LS);
 679         fOutputList->Add(fPtElec_ULS_weight);
 680         fOutputList->Add(fPtElec_LS_weight);
 681
 682         if(fFillBackground){
 683                 fOutputList->Add(fPtElec_ULS2);
 684                 fOutputList->Add(fPtElec_LS2);
 685                 fOutputList->Add(fPtElec_ULS2_weight);
 686                 fOutputList->Add(fPtElec_LS2_weight);
 687         }
 688
 689
 690         fOutputList->Add(fPtTrigger_Inc);
 691         fOutputList->Add(fTPCnsigma_pt_2D);
 692         fOutputList->Add(fShowerShapeCut);
 693
 694         fOutputList->Add(fCharge_n);
 695         fOutputList->Add(fCharge_p);
 696
 697         if(fUseEMCal){
 698
 699                 if(!fIsAOD){
 700                         fOutputList->Add(fTime);
 701                         fOutputList->Add(fTime2);
 702
 703                 }
 704
 705                 fOutputList->Add(ftimingEle);
 706                 fOutputList->Add(ftimingEle2);
 707
 708                 fOutputList->Add(fShowerShape_ha);
 709                 fOutputList->Add(fShowerShape_ele);
 710
 711                 fOutputList->Add(fShowerShapeM02_EoverP);
 712                 fOutputList->Add(fShowerShapeM20_EoverP);
 713
 714
 715
 716         }
 717
 718                 //General Histograms
 719
 720                 //Steps
 721                 //Step 1: Before Track cuts
 722                 //Step 2: Before PID
 723                 //Step 3: After PID
 724
 725         fEoverP_pt = new TH2F *[3];
 726         fTPC_p = new TH2F *[3];
 727         fTPCnsigma_p = new TH2F *[3];
 728
 729         fECluster= new TH1F *[3];
 730         fEtaPhi= new TH2F *[3];
 731         fVtxZ= new  TH1F *[3];
 732         fNTracks= new  TH1F *[3];
 733         fNClusters= new TH1F *[3];
 734         fTPCNcls_EoverP= new TH2F *[3];
 735
 736         for(Int_t i = 0; i < 3; i++)
 737         {
 738                 fEoverP_pt[i] = new TH2F(Form("fEoverP_pt%d",i),";p_{t} (GeV/c);E / p ",1000,0,30,500,0,2);
 739                 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);
 740                 fTPCnsigma_p[i] = new TH2F(Form("fTPCnsigma_p%d",i),";p (GeV/c);TPC Electron N#sigma",1000,0.3,15,1000,-15,10);
 741
 742
 743                 fECluster[i]= new TH1F(Form("fECluster%d",i), ";ECluster",2000, 0,100);
 744                 fEtaPhi[i]= new TH2F(Form("fEtaPhi%d",i),"#eta x #phi Clusters;#phi;#eta",200,0.,5,50,-1.,1.);
 745                 fVtxZ[i]= new  TH1F(Form("fVtxZ%d",i),"VtxZ",1000, -50,50);
 746                 fNTracks[i]= new  TH1F(Form("fNTracks%d",i),"NTracks",1000, 0,1000);
 747                 fNClusters[i]= new TH1F(Form("fNClusters%d",i),"fNClusters0",200, 0,100);
 748                 fTPCNcls_EoverP[i]= new TH2F(Form("fTPCNcls_EoverP%d",i),"TPCNcls_EoverP",1000,0,200,200,0,2);
 749
 750
 751                 fOutputList->Add(fEoverP_pt[i]);
 752                 fOutputList->Add(fTPC_p[i]);
 753                 fOutputList->Add(fTPCnsigma_p[i]);
 754
 755
 756                 fOutputList->Add(fECluster[i]);
 757                 fOutputList->Add(fEtaPhi[i]);
 758                 fOutputList->Add(fVtxZ[i]);
 759                 fOutputList->Add(fNTracks[i]);
 760                 fOutputList->Add(fNClusters[i]);
 761                 fOutputList->Add(fTPCNcls_EoverP[i]);
 762         }
 763
 764                 //pt bin
 765         Int_t fPtBin[7] = {1,2,4,6,8,10,15};
 766
 767         fEoverP_tpc = new TH2F *[6];
 768         fTPC_pt = new TH1F *[6];
 769         fTPCnsigma_pt = new TH1F *[6];
 770
 771         fEta=new TH1F *[6];
 772         fPhi=new TH1F *[6];
 773         fR=new TH1F *[6];
 774         fR_EoverP=new TH2F *[6];
 775         fNcells=new TH1F *[6];
 776         fNcells_EoverP=new TH2F *[6];
 777         fM02_EoverP= new TH2F *[6];
 778         fM20_EoverP= new TH2F *[6];
 779         fEoverP_ptbins=new TH1F *[6];
 780         fECluster_ptbins=new TH1F *[6];
 781         fEoverP_wSSCut=new TH1F *[6];
 782         fNcells_electrons=new TH1F *[6];
 783         fNcells_hadrons=new TH1F *[6];
 784         fTPCnsigma_eta_electrons=new TH2F *[6];
 785         fTPCnsigma_eta_hadrons=new TH2F *[6];
 786
 787         if(fCorrelationFlag)
 788         {
 789                 fCEtaPhi_Inc = new TH2F *[6];
 790                 fCEtaPhi_Inc_DiHadron = new TH2F *[6];
 791
 792                 fCEtaPhi_ULS = new TH2F *[6];
 793                 fCEtaPhi_LS = new TH2F *[6];
 794                 fCEtaPhi_ULS_NoP = new TH2F *[6];
 795                 fCEtaPhi_LS_NoP = new TH2F *[6];
 796
 797                 fCEtaPhi_ULS_Weight = new TH2F *[6];
 798                 fCEtaPhi_LS_Weight = new TH2F *[6];
 799                 fCEtaPhi_ULS_NoP_Weight = new TH2F *[6];
 800                 fCEtaPhi_LS_NoP_Weight = new TH2F *[6];
 801
 802                 fCEtaPhi_Inc_EM = new TH2F *[6];
 803
 804                 fCEtaPhi_ULS_EM = new TH2F *[6];
 805                 fCEtaPhi_LS_EM = new TH2F *[6];
 806
 807                 fCEtaPhi_ULS_Weight_EM = new TH2F *[6];
 808                 fCEtaPhi_LS_Weight_EM = new TH2F *[6];
 809
 810                 fInvMass = new TH1F("fInvMass","",200,0,0.3);
 811                 fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
 812                 fDCA = new TH1F("fDCA","",200,0,1);
 813                 fDCABack = new TH1F("fDCABack","",200,0,1);
 814                 fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
 815                 fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
 816
 817                 fOutputList->Add(fInvMass);
 818                 fOutputList->Add(fInvMassBack);
 819                 fOutputList->Add(fDCA);
 820                 fOutputList->Add(fDCABack);
 821                 fOutputList->Add(fOpAngle);
 822                 fOutputList->Add(fOpAngleBack);
 823         }
 824
 825         if(fFillBackground){
 826
 827                 fInvMass = new TH1F("fInvMass","",200,0,0.3);
 828                 fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
 829                 fDCA = new TH1F("fDCA","",200,0,1);
 830                 fDCABack = new TH1F("fDCABack","",200,0,1);
 831                 fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
 832                 fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
 833
 834                 fOutputList->Add(fInvMass);
 835                 fOutputList->Add(fInvMassBack);
 836                 fOutputList->Add(fDCA);
 837                 fOutputList->Add(fDCABack);
 838                 fOutputList->Add(fOpAngle);
 839                 fOutputList->Add(fOpAngleBack);
 840
 841                         //histos for TPC-only
 842                 fInvMass2 = new TH1F("fInvMass2","",200,0,0.3);
 843                 fInvMassBack2 = new TH1F("fInvMassBack2","",200,0,0.3);
 844                 fDCA2 = new TH1F("fDCA2","",200,0,1);
 845                 fDCABack2 = new TH1F("fDCABack2","",200,0,1);
 846                 fOpAngle2 = new TH1F("fOpAngle2","",200,0,0.5);
 847                 fOpAngleBack2 = new TH1F("fOpAngleBack2","",200,0,0.5);
 848
 849                 fOutputList->Add(fInvMass2);
 850                 fOutputList->Add(fInvMassBack2);
 851                 fOutputList->Add(fDCA2);
 852                 fOutputList->Add(fDCABack2);
 853                 fOutputList->Add(fOpAngle2);
 854                 fOutputList->Add(fOpAngleBack2);
 855
 856         }
 857
 858         for(Int_t i = 0; i < 6; i++)
 859         {
 860                 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);
 861                 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);
 862                 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);
 863
 864                 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);
 865                 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);
 866                 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);
 867                 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);
 868                 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);
 869                 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);
 870                 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);
 871                 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);
 872                 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,100,100,0,2);
 873                 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,100,100,0,2);
 874                 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);
 875                 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);
 876                 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);
 877                 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);
 878                 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);
 879
 880                 fOutputList->Add(fEoverP_tpc[i]);
 881                 fOutputList->Add(fTPC_pt[i]);
 882                 fOutputList->Add(fTPCnsigma_pt[i]);
 883
 884                 fOutputList->Add(fEta[i]);
 885                 fOutputList->Add(fPhi[i]);
 886                 fOutputList->Add(fR[i]);
 887                 fOutputList->Add(fR_EoverP[i]);
 888                 fOutputList->Add(fNcells[i]);
 889                 fOutputList->Add(fNcells_electrons[i]);
 890                 fOutputList->Add(fNcells_hadrons[i]);
 891                 fOutputList->Add(fNcells_EoverP[i]);
 892                 fOutputList->Add(fECluster_ptbins[i]);
 893                 fOutputList->Add(fEoverP_ptbins[i]);
 894                 fOutputList->Add(fEoverP_wSSCut[i]);
 895                 fOutputList->Add(fM02_EoverP[i]);
 896                 fOutputList->Add(fM20_EoverP[i]);
 897                 fOutputList->Add(fTPCnsigma_eta_electrons[i]);
 898                 fOutputList->Add(fTPCnsigma_eta_hadrons[i]);
 899
 900
 901                 if(fCorrelationFlag)
 902                 {
 903                         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);
 904                         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);
 905
 906                         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);
 907                         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);
 908                         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);
 909                         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);
 910
 911                         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);
 912                         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);
 913                         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);
 914                         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);
 915
 916                         fOutputList->Add(fCEtaPhi_Inc[i]);
 917                         fOutputList->Add(fCEtaPhi_Inc_DiHadron[i]);
 918
 919                         fOutputList->Add(fCEtaPhi_ULS[i]);
 920                         fOutputList->Add(fCEtaPhi_LS[i]);
 921                         fOutputList->Add(fCEtaPhi_ULS_NoP[i]);
 922                         fOutputList->Add(fCEtaPhi_LS_NoP[i]);
 923
 924                         fOutputList->Add(fCEtaPhi_ULS_Weight[i]);
 925                         fOutputList->Add(fCEtaPhi_LS_Weight[i]);
 926                         fOutputList->Add(fCEtaPhi_ULS_NoP_Weight[i]);
 927                         fOutputList->Add(fCEtaPhi_LS_NoP_Weight[i]);
 928
 929                         if(fEventMixingFlag)
 930                         {
 931                                 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);
 932
 933                                 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);
 934                                 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);
 935
 936                                 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);
 937                                 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);
 938
 939                                 fOutputList->Add(fCEtaPhi_Inc_EM[i]);
 940
 941                                 fOutputList->Add(fCEtaPhi_ULS_EM[i]);
 942                                 fOutputList->Add(fCEtaPhi_LS_EM[i]);
 943
 944                                 fOutputList->Add(fCEtaPhi_ULS_Weight_EM[i]);
 945                                 fOutputList->Add(fCEtaPhi_LS_Weight_EM[i]);
 946                         }
 947                 }
 948         }
 949
 950                 //pt integrated
 951         fTPCnsigma_eta = new TH2F("fTPCnsigma_eta",";Pseudorapidity #eta; TPC signal - <TPC signal>_{elec} (#sigma)",200,-0.9,0.9,200,-15,15);
 952         fTPCnsigma_phi = new TH2F("fTPCnsigma_phi",";Azimuthal Angle #phi; TPC signal - <TPC signal>_{elec} (#sigma)",200,0,2*TMath::Pi(),200,-15,15);
 953
 954
 955         fNcells_pt=new TH2F("fNcells_pt","fNcells_pt",1000, 0,20,100,0,30);
 956         fEoverP_pt_pions= new TH2F("fEoverP_pt_pions","fEoverP_pt_pions",1000,0,30,500,0,2);
 957
 958         ftpc_p_EoverPcut= new TH2F("ftpc_p_EoverPcut","ftpc_p_EoverPcut",1000,0,30,200,20,200);
 959         fnsigma_p_EoverPcut= new TH2F("fnsigma_p_EoverPcut","fnsigma_p_EoverPcut",1000,0,30,500,-15,15);
 960
 961         fEoverP_pt_pions2= new TH2F("fEoverP_pt_pions2","fEoverP_pt_pions2",1000,0,30,500,0,2);
 962         fEoverP_pt_hadrons= new TH2F("fEoverP_pt_hadrons","fEoverP_pt_hadrons",1000,0,30,500,0,2);
 963
 964
 965         fOutputList->Add(fTPCnsigma_eta);
 966         fOutputList->Add(fTPCnsigma_phi);
 967
 968         fOutputList->Add(fNcells_pt);
 969         fOutputList->Add(fEoverP_pt_pions);
 970
 971         fOutputList->Add(ftpc_p_EoverPcut);
 972         fOutputList->Add(fnsigma_p_EoverPcut);
 973
 974         fOutputList->Add(fEoverP_pt_pions2);
 975         fOutputList->Add(fEoverP_pt_hadrons);
 976
 977                 //__________________________________________________________________
 978                 //Efficiency studies
 979         if(fIsMC)
 980         {
 981                 fPtBackgroundBeforeReco = new TH1F("fPtBackgroundBeforeReco",";p_{T} (GeV/c);Count",300,0,30);
 982                 fPtBackgroundBeforeReco_weight = new TH1F("fPtBackgroundBeforeReco_weight",";p_{T} (GeV/c);Count",300,0,30);
 983                 if(fFillBackground)fPtBackgroundBeforeReco2 = new TH1F("fPtBackgroundBeforeReco2",";p_{T} (GeV/c);Count",300,0,30);
 984                 if(fFillBackground)fPtBackgroundBeforeReco2_weight = new TH1F("fPtBackgroundBeforeReco2_weight",";p_{T} (GeV/c);Count",300,0,30);
 985
 986                 fPtBackgroundAfterReco = new TH1F("fPtBackgroundAfterReco",";p_{T} (GeV/c);Count",300,0,30);
 987                 fPtMCparticleAll = new TH1F("fPtMCparticleAll",";p_{T} (GeV/c);Count",200,0,40);
 988                 fPtMCparticleReco = new TH1F("fPtMCparticleReco",";p_{T} (GeV/c);Count",200,0,40);
 989
 990                 fPtMCparticleAll_nonPrimary = new TH1F("fPtMCparticleAll_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
 991                 fPtMCparticleAlle_nonPrimary = new TH1F("fPtMCparticleAlle_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
 992                 fPtMCparticleAlle_Primary = new TH1F("fPtMCparticleAlle_Primary",";p_{T} (GeV/c);Count",200,0,40);
 993
 994                 fPtMCparticleReco_nonPrimary = new TH1F("fPtMCparticleReco_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
 995
 996                 fPtMCparticleAllHfe1 = new TH1F("fPtMCparticleAllHfe1",";p_{t} (GeV/c);Count",200,0,40);
 997                 fPtMCparticleRecoHfe1 = new TH1F("fPtMCparticleRecoHfe1",";p_{t} (GeV/c);Count",200,0,40);
 998                 fPtMCparticleAllHfe2 = new TH1F("fPtMCparticleAllHfe2",";p_{t} (GeV/c);Count",200,0,40);
 999                 fPtMCparticleRecoHfe2 = new TH1F("fPtMCparticleRecoHfe2",";p_{t} (GeV/c);Count",200,0,40);
1000
1001                 fPtMCelectronAfterAll = new TH1F("fPtMCelectronAfterAll",";p_{T} (GeV/c);Count",200,0,40);
1002                 fPtMCelectronAfterAll_nonPrimary = new TH1F("fPtMCelectronAfterAll_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
1003                 fPtMCelectronAfterAll_Primary = new TH1F("fPtMCelectronAfterAll_Primary",";p_{T} (GeV/c);Count",200,0,40);
1004
1005
1006
1007                 fPtMCpi0 = new TH1F("fPtMCpi0",";p_{t} (GeV/c);Count",200,0,30);
1008                 fPtMCeta = new TH1F("fPtMCeta",";p_{T} (GeV/c);Count",200,0,30);
1009                 fPtMC_EMCal_All= new TH1F("fPtMC_EMCal_All",";p_{t} (GeV/c);Count",200,0,40);
1010                 fPtMC_EMCal_Selected= new TH1F("fPtMC_EMCal_Selected",";p_{t} (GeV/c);Count",200,0,40);
1011                 fPtMC_TPC_All= new TH1F("fPtMC_TPC_All",";p_{t} (GeV/c);Count",200,0,40);
1012                 fPtMC_TPC_Selected = new TH1F("fPtMC_TPC_Selected",";p_{t} (GeV/c);Count",200,0,40);
1013                 fPtMCWithLabel = new TH1F("fPtMCWithLabel",";p_{t} (GeV/c);Count",200,0,40);
1014                 fPtMCWithoutLabel = new TH1F("fPtMCWithoutLabel",";p_{t} (GeV/c);Count",200,0,40);
1015                 fPtIsPhysicaPrimary = new TH1F("fPtIsPhysicaPrimary",";p_{t} (GeV/c);Count",200,0,40);
1016
1017                 fOutputList->Add(fPtBackgroundBeforeReco);
1018                 fOutputList->Add(fPtBackgroundBeforeReco_weight);
1019
1020                 if(fFillBackground) fOutputList->Add(fPtBackgroundBeforeReco2);
1021                 if(fFillBackground) fOutputList->Add(fPtBackgroundBeforeReco2_weight);
1022                 fOutputList->Add(fPtBackgroundAfterReco);
1023                 fOutputList->Add(fPtMCparticleAll);
1024                 fOutputList->Add(fPtMCparticleReco);
1025
1026                 fOutputList->Add(fPtMCparticleAll_nonPrimary);
1027                 fOutputList->Add(fPtMCparticleAlle_nonPrimary);
1028
1029                 fOutputList->Add(fPtMCparticleAlle_Primary);
1030                 fOutputList->Add(fPtMCparticleReco_nonPrimary);
1031
1032                 fOutputList->Add(fPtMCparticleAllHfe1);
1033                 fOutputList->Add(fPtMCparticleRecoHfe1);
1034                 fOutputList->Add(fPtMCparticleAllHfe2);
1035                 fOutputList->Add(fPtMCparticleRecoHfe2);
1036                 fOutputList->Add(fPtMCelectronAfterAll);
1037
1038                 fOutputList->Add(fPtMCelectronAfterAll_nonPrimary);
1039                 fOutputList->Add(fPtMCelectronAfterAll_Primary);
1040
1041
1042
1043                 fOutputList->Add(fPtMCpi0);
1044                 fOutputList->Add(fPtMCeta);
1045                 fOutputList->Add(fPtMC_EMCal_All);
1046                 fOutputList->Add(fPtMC_EMCal_Selected);
1047                 fOutputList->Add(fPtMC_TPC_All);
1048                 fOutputList->Add(fPtMC_TPC_Selected);
1049                 fOutputList->Add(fPtMCWithLabel);
1050                 fOutputList->Add(fPtMCWithoutLabel);
1051                 fOutputList->Add(fPtIsPhysicaPrimary);
1052         }
1053
1054         fCentralityHist = new TH1F("fCentralityHist",";Centrality (%); Count",1000000,0,100);
1055         fCentralityHistPass = new TH1F("fCentralityHistPass",";Centrality (%); Count",1000000,0,100);
1056         fOutputList->Add(fCentralityHist);
1057         fOutputList->Add(fCentralityHistPass);
1058
1059                 //______________________________________________________________________
1060                 //Mixed event analysis
1061         if(fEventMixingFlag)
1062         {
1063                 fPoolNevents = new TH1F("fPoolNevents","Event Mixing Statistics; Number of events; Count",1000,0,1000);
1064                 fOutputList->Add(fPoolNevents);
1065
1066                 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.
1067                 Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
1068
1069                 Int_t nCentralityBins  = 15;
1070                 Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
1071
1072                 Int_t nZvtxBins  = 9;
1073                 Double_t vertexBins[] = {-10, -7, -5, -3, -1, 1, 3, 5, 7, 10};
1074
1075                 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, (Double_t*) centralityBins, nZvtxBins, (Double_t*) vertexBins);
1076         }
1077                 //______________________________________________________________________
1078
1079         PostData(1, fOutputList);
1080
1081                 ///______________________________________________________________________
1082 }
1083
1084         //______________________________________________________________________
1085         //Main loop
1086         //Called for each event
1087 void AliAnalysisTaskEMCalHFEpA::UserExec(Option_t *)
1088 {
1089                 //Check Event
1090         fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1091         fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1092
1093         if(!(fESD || fAOD))
1094         {
1095                 printf("ERROR: fESD & fAOD not available\n");
1096                 return;
1097         }
1098
1099         fVevent = dynamic_cast<AliVEvent*>(InputEvent());
1100
1101         if(!fVevent)
1102         {
1103                 printf("ERROR: fVEvent not available\n");
1104                 return;
1105         }
1106
1107                 //Check Cuts
1108         if(!fCuts)
1109         {
1110                 AliError("HFE cuts not available");
1111                 return;
1112         }
1113                 //Check PID
1114         if(!fPID->IsInitialized())
1115         {
1116                         // Initialize PID with the given run number
1117                 AliWarning("PID not initialised, get from Run no");
1118
1119                 if(fIsAOD)
1120                 {
1121                         fPID->InitializePID(fAOD->GetRunNumber());
1122                 }
1123                 else
1124                 {
1125                         fPID->InitializePID(fESD->GetRunNumber());
1126                 }
1127         }
1128
1129                 //PID response
1130         fPidResponse = fInputHandler->GetPIDResponse();
1131
1132
1133                 //Check PID response
1134         if(!fPidResponse)
1135         {
1136                 AliDebug(1, "Using default PID Response");
1137                 fPidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
1138         }
1139
1140         fPID->SetPIDResponse(fPidResponse);
1141
1142         fCFM->SetRecEventInfo(fVevent);
1143
1144         Double_t *fListOfmotherkink = 0;
1145         Int_t fNumberOfVertices = 0;
1146         Int_t fNumberOfMotherkink = 0;
1147
1148                 //______________________________________________________________________
1149                 //Vertex Selection
1150         if(fIsAOD)
1151         {
1152                 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
1153                 if(!trkVtx || trkVtx->GetNContributors()<=0) return;
1154                 TString vtxTtl = trkVtx->GetTitle();
1155                 if(!vtxTtl.Contains("VertexerTracks")) return;
1156                 Float_t zvtx = trkVtx->GetZ();
1157                 fZvtx = zvtx;
1158                 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
1159                 if(spdVtx->GetNContributors()<=0) return;
1160                 TString vtxTyp = spdVtx->GetTitle();
1161                 Double_t cov[6]={0};
1162                 spdVtx->GetCovarianceMatrix(cov);
1163                 Double_t zRes = TMath::Sqrt(cov[5]);
1164                 if(vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1165                 if(TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1166                 if(TMath::Abs(zvtx) > 10) return;
1167
1168                         //Look for kink mother for AOD
1169
1170                 fNumberOfVertices = 0;
1171                 fNumberOfMotherkink = 0;
1172
1173                 if(fIsAOD)
1174                 {
1175                         fNumberOfVertices = fAOD->GetNumberOfVertices();
1176
1177                         fListOfmotherkink = new Double_t[fNumberOfVertices];
1178
1179                         for(Int_t ivertex=0; ivertex < fNumberOfVertices; ivertex++)
1180                         {
1181                                 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
1182                                 if(!aodvertex) continue;
1183                                 if(aodvertex->GetType()==AliAODVertex::kKink)
1184                                 {
1185                                         AliAODTrack *mother1 = (AliAODTrack *) aodvertex->GetParent();
1186                                         if(!mother1) continue;
1187                                         Int_t idmother = mother1->GetID();
1188                                         fListOfmotherkink[fNumberOfMotherkink] = idmother;
1189                                         fNumberOfMotherkink++;
1190                                 }
1191                         }
1192                 }
1193         }
1194         else
1195         {
1196
1197
1198
1199                         /// ESD
1200                 const AliESDVertex* trkVtx = fESD->GetPrimaryVertex();
1201                 if(!trkVtx || trkVtx->GetNContributors()<=0) return;
1202                 TString vtxTtl = trkVtx->GetTitle();
1203                 if(!vtxTtl.Contains("VertexerTracks")) return;
1204                 Float_t zvtx = trkVtx->GetZ();
1205
1206                 const AliESDVertex* spdVtx = fESD->GetPrimaryVertexSPD();
1207                 if(spdVtx->GetNContributors()<=0) return;
1208                 TString vtxTyp = spdVtx->GetTitle();
1209                 Double_t cov[6]={0};
1210                 spdVtx->GetCovarianceMatrix(cov);
1211                 Double_t zRes = TMath::Sqrt(cov[5]);
1212                 if(vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1213                 if(TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1214                 if(TMath::Abs(zvtx) > 10) return;
1215         }
1216
1217                 //______________________________________________________________________
1218
1219                 //Only events with at least 2 tracks are accepted
1220         Int_t fNOtrks =  fVevent->GetNumberOfTracks();
1221         if(fNOtrks<2) return;
1222
1223                 //______________________________________________________________________
1224                 //Centrality Selection
1225         if(fHasCentralitySelection)
1226         {
1227                 Float_t centrality = -1;
1228
1229                 if(fIsAOD)
1230                 {
1231                         fCentrality = fAOD->GetHeader()->GetCentralityP();
1232                 }
1233                 else
1234                 {
1235                         fCentrality = fESD->GetCentrality();
1236                 }
1237
1238                 if(fEstimator==1) centrality = fCentrality->GetCentralityPercentile("ZDC");
1239                 else centrality = fCentrality->GetCentralityPercentile("V0A");
1240
1241                         //cout << "Centrality = " << centrality << " %" << endl;
1242
1243                 fCentralityHist->Fill(centrality);
1244
1245                 if(centrality<fCentralityMin || centrality>fCentralityMax) return;
1246
1247                 fCentralityHistPass->Fill(centrality);
1248         }
1249                 //______________________________________________________________________
1250
1251                 //______________________________________________________________________
1252
1253         if(fIsMC)
1254         {
1255                 if(fIsAOD)
1256                 {
1257                         fMCarray = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1258
1259                         if(!fMCarray)
1260                         {
1261                                 AliError("Array of MC particles not found");
1262                                 return;
1263                         }
1264
1265                         fMCheader = dynamic_cast<AliAODMCHeader*>(fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
1266
1267                         if(!fMCheader)
1268                         {
1269                                 AliError("Could not find MC Header in AOD");
1270                                 return;
1271                         }
1272
1273                         for(Int_t iMC = 0; iMC < fMCarray->GetEntries(); iMC++)
1274                         {
1275                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(iMC);
1276                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1277
1278                                 Int_t pdg = fMCparticle->GetPdgCode();
1279                                 Int_t mpdg = fMCparticleMother->GetPdgCode();
1280
1281                                 double proX = fMCparticle->Xv();
1282                                 double proY = fMCparticle->Yv();
1283                                 double proR = sqrt(pow(proX,2)+pow(proY,2));
1284
1285
1286                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax && fMCparticle->Charge()!=0)
1287                                 {
1288                                                 //to correct background
1289                                         if (TMath::Abs(pdg) == 11){
1290                                                 if(TMath::Abs(mpdg) == 221 || TMath::Abs(mpdg) == 22 || TMath::Abs(mpdg) == 111){
1291
1292                                                         if(proR<7){
1293                                                                 fPtMCparticleAlle_nonPrimary->Fill(fMCparticle->Pt()); //denominator for total efficiency for all electrons, and not primary
1294
1295                                                         }
1296                                                 }
1297                                         }
1298
1299                                         if (TMath::Abs(pdg) == 11 && fMCparticle->IsPhysicalPrimary()) fPtMCparticleAlle_Primary->Fill(fMCparticle->Pt()); //denominator for total efficiency for all electrons primary
1300
1301                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1302                                         {
1303
1304                                                 fPtMCparticleAll_nonPrimary->Fill(fMCparticle->Pt()); //denominator for total efficiency for all particles, and not primary
1305                                                 if(fMCparticle->IsPhysicalPrimary())
1306                                                 {
1307                                                         fPtMCparticleAll->Fill(fMCparticle->Pt());
1308
1309                                                         Bool_t MotherFound = FindMother(iMC);
1310                                                         if(MotherFound)
1311                                                         {
1312                                                                 if(fIsHFE1) fPtMCparticleAllHfe1->Fill(fMCparticle->Pt()); //denominator for total efficiency
1313                                                                 if(fIsHFE2) fPtMCparticleAllHfe2->Fill(fMCparticle->Pt());
1314                                                         }
1315                                                 }
1316                                         }
1317                                 }//eta cut
1318
1319                                 if(TMath::Abs(pdg)==111) fPtMCpi0->Fill(fMCparticle->Pt());
1320                                 if(TMath::Abs(pdg)==221) fPtMCeta->Fill(fMCparticle->Pt());
1321                         }//loop tracks
1322                 }//AOD
1323                 else
1324                 {
1325                         fEventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1326                         if (!fEventHandler) {
1327                                 Printf("ERROR: Could not retrieve MC event handler");
1328                                 return;
1329                         }
1330
1331                         fMCevent = fEventHandler->MCEvent();
1332                         if (!fMCevent) {
1333                                 Printf("ERROR: Could not retrieve MC event");
1334                                 return;
1335                         }
1336
1337                         fMCstack = fMCevent->Stack();
1338
1339                 for(Int_t iMC = 0; iMC < fMCstack->GetNtrack(); iMC++)
1340                 {
1341
1342                                 fMCtrack = fMCstack->Particle(iMC);
1343                                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1344                                 TParticle *particle=fMCstack->Particle(iMC);
1345
1346                                 Int_t pdg = fMCtrack->GetPdgCode();
1347                                 Int_t mpdg = fMCtrackMother->GetPdgCode();
1348
1349                                 if(TMath::Abs(pdg)==111) fPtMCpi0->Fill(fMCtrack->Pt());
1350                                 if(TMath::Abs(pdg)==221) fPtMCeta->Fill(fMCtrack->Pt());
1351
1352
1353                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
1354                                 {
1355
1356                                                 //to correct background
1357                                         if (TMath::Abs(pdg) == 11){
1358                                                 if(TMath::Abs(mpdg) == 221 || TMath::Abs(mpdg) == 22 || TMath::Abs(mpdg) == 111){
1359                                                         Double_t proR=particle->R();
1360                                                         if(proR<7){
1361                                                                 fPtMCparticleAlle_nonPrimary->Fill(fMCtrack->Pt()); //denominator for total efficiency for all electrons, and not primary
1362                                                         }
1363                                                 }
1364                                         }
1365
1366                                         if (TMath::Abs(pdg) == 11 && fMCstack->IsPhysicalPrimary(iMC))  fPtMCparticleAlle_Primary->Fill(fMCtrack->Pt());
1367
1368
1369                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1370                                         {
1371                                                 fPtMCparticleAll_nonPrimary->Fill(fMCtrack->Pt());//denominator for total efficiency for all particle, non Primary track
1372
1373                                                 if(!fMCstack->IsPhysicalPrimary(iMC)) continue;
1374                                                 fPtMCparticleAll->Fill(fMCtrack->Pt());
1375
1376                                                 Bool_t MotherFound = FindMother(iMC);
1377                                                 if(MotherFound)
1378                                                 {
1379                                                         if(fIsHFE1) fPtMCparticleAllHfe1->Fill(fMCtrack->Pt());//denominator for total efficiency
1380                                                         if(fIsHFE2) fPtMCparticleAllHfe2->Fill(fMCtrack->Pt());
1381                                                 }
1382                                         }
1383                                 }//particle kind
1384                 }//loop tracks
1385                 }//ESD
1386         }//Is MC
1387
1388                 //______________________________________________________________________
1389                 //EMCal Trigger Selection (Threshould selection)
1390         TString firedTrigger;
1391         TString TriggerEG1("EG1");
1392         TString TriggerEG2("EG2");
1393
1394         if(fAOD) firedTrigger = fAOD->GetFiredTriggerClasses();
1395         else if(fESD) firedTrigger = fESD->GetFiredTriggerClasses();
1396
1397         fNevent->Fill(0);
1398         if(firedTrigger.Contains(TriggerEG1)) fNevent->Fill(1);
1399         if(firedTrigger.Contains(TriggerEG2)) fNevent->Fill(2);
1400
1401                 //EG1
1402         if(firedTrigger.Contains(TriggerEG1))
1403         {
1404                 fNevent->Fill(3);
1405         }
1406         else
1407         {
1408                 if(fEMCEG1) return;
1409         }
1410
1411                 //EG2
1412         if(firedTrigger.Contains(TriggerEG2))
1413         {
1414                 fNevent->Fill(4);
1415         }
1416         else
1417         {
1418                 if(fEMCEG2) return;
1419         }
1420
1421                 //______________________________________________________________________
1422
1423         Int_t ClsNo = -999;
1424         if(!fIsAOD) ClsNo = fESD->GetNumberOfCaloClusters();
1425         else ClsNo = fAOD->GetNumberOfCaloClusters();
1426
1427                 //______________________________________________________________________
1428
1429                 ///______________________________________________________________________
1430                 ///Track loop
1431         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
1432         {
1433                 AliVParticle* Vtrack = fVevent->GetTrack(iTracks);
1434                 if (!Vtrack)
1435                 {
1436                         printf("ERROR: Could not receive track %d\n", iTracks);
1437                         continue;
1438                 }
1439
1440                 AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
1441                 AliESDtrack *etrack = dynamic_cast<AliESDtrack*>(Vtrack);
1442                 AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(Vtrack);
1443
1444                 Double_t fTPCnSigma = -999;
1445                 Double_t fTOFnSigma = -999;
1446                 Double_t fTPCnSigma_pion = -999;
1447                 Double_t fTPCnSigma_proton = -999;
1448                 Double_t fTPCnSigma_kaon = -999;
1449                 Double_t fTPCsignal = -999;
1450                 Double_t fPt = -999;
1451                 Double_t fP = -999;
1452
1453                         ///_____________________________________________________________________________
1454                         ///Fill QA plots without track selection
1455                 fPt = track->Pt();
1456                 fP = TMath::Sqrt((track->Pt())*(track->Pt()) + (track->Pz())*(track->Pz()));
1457
1458                 fTPCsignal = track->GetTPCsignal();
1459                 fTPCnSigma = fPidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
1460                 fTOFnSigma = fPidResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
1461                 fTPCnSigma_pion = fPidResponse->NumberOfSigmasTPC(track, AliPID::kPion);
1462                 fTPCnSigma_proton = fPidResponse->NumberOfSigmasTPC(track, AliPID::kProton);
1463                 fTPCnSigma_kaon = fPidResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
1464
1465                 fTPC_p[0]->Fill(fPt,fTPCsignal);
1466                 fTPCnsigma_p[0]->Fill(fP,fTPCnSigma);
1467
1468
1469                 Float_t TPCNcls = track->GetTPCNcls();
1470                 Float_t pos[3]={0,0,0};
1471
1472                 Double_t fEMCflag = kFALSE;
1473                 if(track->GetEMCALcluster()>0)
1474                 {
1475                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
1476                         if(fClus->IsEMCAL())
1477                         {
1478                                 if(TMath::Abs(fClus->GetTrackDx())<=fdPhiCut && TMath::Abs(fClus->GetTrackDz())<=fdEtaCut)
1479                                 {
1480                                         fEMCflag = kTRUE;
1481                                         fEoverP_pt[0]->Fill(fPt,(fClus->E() / fP));
1482
1483
1484                                         Float_t Energy  = fClus->E();
1485                                         Float_t EoverP  = Energy/track->P();
1486                                                 //Float_t M02   = fClus->GetM02();
1487                                                 //Float_t M20   = fClus->GetM20();
1488
1489                                                 /////////////// for Eta Phi distribution
1490                                         fClus->GetPosition(pos);
1491                                         TVector3 vpos(pos[0],pos[1],pos[2]);
1492                                         Double_t cphi = vpos.Phi();
1493                                         Double_t ceta = vpos.Eta();
1494                                         fEtaPhi[0]->Fill(cphi,ceta);
1495
1496                                         fECluster[0]->Fill(Energy);
1497                                         fTPCNcls_EoverP[0]->Fill(TPCNcls, EoverP);
1498                                 }
1499                         }
1500                 }
1501
1502                         //______________________________________________________________
1503                         // Vertex
1504
1505                 fVtxZ[0]->Fill(fZvtx);
1506                 fNTracks[0]->Fill(fNOtrks);
1507                 fNClusters[0]->Fill(ClsNo);
1508
1509                         //______________________________________________________________
1510
1511                         ///Fill QA plots without track selection
1512                         ///_____________________________________________________________________________
1513                         //______________________________________________________________________________________
1514                         //Track Selection Cuts
1515
1516                         //AOD (Test Filter Bit)
1517                 if(fIsAOD)
1518                 {
1519                                 // standard cuts with very loose DCA - BIT(4)
1520                                 // Description:
1521                         /*
1522                          GetStandardITSTPCTrackCuts2011(kFALSE)
1523                          SetMaxChi2PerClusterTPC(4);
1524                          SetAcceptKinkDaughters(kFALSE);
1525                          SetRequireTPCRefit(kTRUE);
1526                          SetRequireITSRefit(kTRUE);
1527                          SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
1528                          SetMaxDCAToVertexZ(2);
1529                          SetMaxDCAToVertex2D(kFALSE);
1530                          SetRequireSigmaToVertex(kFALSE);
1531                          SetMaxChi2PerClusterITS(36);
1532                          SetMaxDCAToVertexXY(2.4)
1533                          SetMaxDCAToVertexZ(3.2)
1534                          SetDCaToVertex2D(kTRUE)
1535                          */
1536
1537                         if(!atrack->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;
1538                 }
1539
1540                         //RecKine: ITSTPC cuts
1541                 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
1542                         //RecKink
1543                 if(fRejectKinkMother)
1544                 {
1545                         if(fIsAOD)
1546                         {
1547                                 Bool_t kinkmotherpass = kTRUE;
1548                                 for(Int_t kinkmother = 0; kinkmother < fNumberOfMotherkink; kinkmother++)
1549                                 {
1550                                         if(track->GetID() == fListOfmotherkink[kinkmother])
1551                                         {
1552                                                 kinkmotherpass = kFALSE;
1553                                                 continue;
1554                                         }
1555                                 }
1556                                 if(!kinkmotherpass) continue;
1557                         }
1558                         else
1559                         {
1560                                 if(etrack->GetKinkIndex(0) != 0) continue;
1561                         }
1562                 }
1563
1564                         //RecPrim
1565                 if(!fIsAOD)
1566                 {
1567                         if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
1568                 }
1569
1570                         //HFEcuts: ITS layers cuts
1571                 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
1572
1573                         //HFE cuts: TPC PID cleanup
1574                 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
1575                         //______________________________________________________________________________________
1576
1577                         ///_____________________________________________________________
1578                         ///QA plots after track selection
1579                 if(fIsMC)
1580                 {
1581                         if(track->GetLabel()>=0) fPtMCWithLabel->Fill(fPt);
1582                         else fPtMCWithoutLabel->Fill(fPt);
1583                 }
1584
1585                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
1586                 {
1587                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1588
1589                         if(fMCparticle->IsPhysicalPrimary()) fPtIsPhysicaPrimary->Fill(fPt);
1590
1591                         Int_t pdg = fMCparticle->GetPdgCode();
1592                         if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax && fMCparticle->Charge()!=0)
1593                         {
1594
1595
1596                                 if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1597                                 {
1598                                         fPtMCparticleReco_nonPrimary->Fill(fMCparticle->Pt()); //not Primary track
1599
1600                                         if(fMCparticle->IsPhysicalPrimary())
1601                                         {
1602                                                 fPtMCparticleReco->Fill(fMCparticle->Pt());
1603
1604                                                 Bool_t MotherFound = FindMother(track->GetLabel());
1605                                                 if(MotherFound)
1606                                                 {
1607                                                         if(fIsHFE1) fPtMCparticleRecoHfe1->Fill(fMCparticle->Pt());
1608                                                         if(fIsHFE2) fPtMCparticleRecoHfe2->Fill(fMCparticle->Pt());
1609                                                 }
1610                                         }
1611                                 }
1612                         }
1613                 }
1614                 else if(fIsMC && track->GetLabel()>=0)
1615                 {
1616                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
1617                         {
1618
1619                                 fMCtrack = fMCstack->Particle(track->GetLabel());
1620                                 Int_t pdg = fMCtrack->GetPdgCode();
1621
1622                                 if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1623                                 {
1624                                         fPtMCparticleReco_nonPrimary->Fill(fMCtrack->Pt());//not Primary track
1625                                 }
1626
1627
1628                                 if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
1629                                 {
1630                                                 fPtIsPhysicaPrimary->Fill(fPt);
1631
1632
1633
1634
1635                                                 if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1636                                                 {
1637                                                         fPtMCparticleReco->Fill(fMCtrack->Pt());
1638
1639                                                         Bool_t MotherFound = FindMother(track->GetLabel());
1640                                                         if(MotherFound)
1641                                                         {
1642                                                                 if(fIsHFE1) fPtMCparticleRecoHfe1->Fill(fMCtrack->Pt());
1643                                                                 if(fIsHFE2) fPtMCparticleRecoHfe2->Fill(fMCtrack->Pt());
1644                                                         }
1645                                                 }
1646                                 }
1647                         }
1648                 }
1649
1650                 fTPC_p[1]->Fill(fPt,fTPCsignal);
1651                 fTPCnsigma_p[1]->Fill(fP,fTPCnSigma);
1652                 Double_t fPtBin[7] = {1,2,4,6,8,10,15};
1653
1654                 TPCNcls = track->GetTPCNcls();
1655                 Float_t pos2[3]={0,0,0};
1656
1657                 if(track->GetEMCALcluster()>0)
1658                 {
1659                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
1660                         if(fClus->IsEMCAL())
1661                         {
1662                                 if(TMath::Abs(fClus->GetTrackDx())<=fdPhiCut && TMath::Abs(fClus->GetTrackDz())<=fdEtaCut)
1663                                 {
1664                                         fEoverP_pt[1]->Fill(fPt,(fClus->E() / fP));
1665
1666                                         Float_t Energy  = fClus->E();
1667                                         Float_t EoverP  = Energy/track->P();
1668                                         Float_t M02     = fClus->GetM02();
1669                                         Float_t M20     = fClus->GetM20();
1670                                         Float_t ncells  = fClus->GetNCells();
1671
1672                                                 /////////////// for Eta Phi distribution
1673                                         fClus->GetPosition(pos2);
1674                                         TVector3 vpos(pos2[0],pos2[1],pos2[2]);
1675                                         Double_t cphi = vpos.Phi();
1676                                         Double_t ceta = vpos.Eta();
1677                                         fEtaPhi[1]->Fill(cphi,ceta);
1678
1679                                         fECluster[1]->Fill(Energy);
1680                                         fTPCNcls_EoverP[1]->Fill(TPCNcls, EoverP);
1681
1682
1683                                                 //for EMCal trigger performance
1684                                         if(EoverP > 0.9){
1685                                                 ftpc_p_EoverPcut->Fill(track->P(), fTPCsignal);
1686                                                 fnsigma_p_EoverPcut->Fill(track->P(), fTPCnSigma);
1687
1688                                         }
1689
1690
1691                                                 //for hadron contamination calculations
1692
1693
1694                                                 // EtaCut -> dados
1695                                         if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
1696                                                         //main
1697                                                 if(TMath::Abs(fTPCnSigma_pion)<3 || TMath::Abs(fTPCnSigma_proton)<3 || TMath::Abs(fTPCnSigma_kaon)<3 ){
1698
1699                                                         if(fTPCnSigma<-3.5){
1700                                                                 fEoverP_pt_hadrons->Fill(fPt,EoverP);
1701                                                                 if(fUseEMCal) fShowerShape_ha->Fill(M02,M20);
1702                                                         }
1703                                                 }
1704
1705                                                 if(fTPCnSigma < -3.5){
1706                                                         fEoverP_pt_pions->Fill(fPt, EoverP);
1707
1708                                                 }
1709
1710                                                 if(fTPCnSigma < -3.5 && fTPCnSigma > -10){
1711                                                         fEoverP_pt_pions2->Fill(fPt, EoverP);
1712
1713                                                 }
1714
1715
1716                                         }
1717
1718
1719
1720
1721                                         for(Int_t i = 0; i < 6; i++)
1722                                         {
1723                                                 if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
1724                                                 {
1725
1726                                                         if(fTPCnSigma < -5 && fTPCnSigma > -10){
1727                                                                 fNcells_hadrons[i]->Fill(ncells);
1728                                                         }
1729                                                                 //hadrons selection using E/p
1730                                                         if((fClus->E() / fP) >= 0.2 && (fClus->E() / fP) <= 0.7){
1731                                                                 fTPCnsigma_eta_hadrons[i]->Fill(fTPCnSigma, ceta);
1732                                                         }
1733                                                                 //electrons selection using E/p
1734                                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax) {
1735                                                                 fTPCnsigma_eta_electrons[i]->Fill(fTPCnSigma, ceta);
1736                                                         }
1737                                                 }
1738                                         }
1739
1740                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax)
1741                                         {
1742                                                 fTPCnsigma_eta->Fill(track->Eta(),fTPCnSigma);
1743                                                 fTPCnsigma_phi->Fill(track->Phi(),fTPCnSigma);
1744
1745                                                 if(fUseEMCal)
1746                                                 {
1747                                                         if(fTPCnSigma < 3.5 && fCorrelationFlag)
1748                                                         {
1749                                                                 fPtTrigger_Inc->Fill(fPt);
1750                                                                 DiHadronCorrelation(track, iTracks);
1751                                                         }
1752                                                 }
1753                                         }
1754
1755                                 }
1756                         }
1757                 }
1758
1759                         //______________________________________________________________
1760                         // Vertex
1761
1762                 fVtxZ[1]->Fill(fZvtx);
1763                 fNTracks[1]->Fill(fNOtrks);
1764                 fNClusters[1]->Fill(ClsNo);
1765                         //______________________________________________________________
1766
1767                         ///______________________________________________________________________
1768                         ///Histograms for PID Studies
1769                         //Double_t fPtBin[6] = {2,4,6,8,10,15};
1770
1771                 for(Int_t i = 0; i < 6; i++)
1772                 {
1773                         if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
1774                         {
1775                                 if(fEMCflag) fEoverP_tpc[i]->Fill(fTPCnSigma,(fClus->E() / fP));
1776
1777
1778                                 fTPC_pt[i]->Fill(fTPCsignal);
1779                                 fTPCnsigma_pt[i]->Fill(fTPCnSigma);
1780
1781                         }
1782                 }
1783
1784                         ///QA plots after track selection
1785                         ///_____________________________________________________________
1786
1787                         //_______________________________________________________
1788                         //Correlation Analysis - DiHadron
1789                 if(!fUseEMCal)
1790                 {
1791                         if(fTPCnSigma < 3.5 && fCorrelationFlag)
1792                         {
1793                                 fPtTrigger_Inc->Fill(fPt);
1794                                 DiHadronCorrelation(track, iTracks);
1795                         }
1796                 }
1797                         //_______________________________________________________
1798
1799
1800                         ///////////////////////////////////////////////////////////////////
1801                         ///TPC - efficiency calculation //
1802
1803                         /// changing start here
1804                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
1805                 {
1806                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1807                         Int_t pdg = fMCparticle->GetPdgCode();
1808
1809                                 //
1810                         if(fMCparticle->IsPhysicalPrimary()){
1811
1812
1813                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
1814
1815                                         Bool_t MotherFound = FindMother(track->GetLabel());
1816                                         if(MotherFound){
1817                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1818                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
1819                                                         if(fIsHFE1) fPtMC_TPC_All->Fill(fMCparticle->Pt());
1820                                         }
1821                                         }
1822                                 }
1823                         }
1824                 }///until here
1825
1826                 else if(fIsMC && track->GetLabel()>=0)//ESD
1827                 {
1828
1829                         if(fMCstack->IsPhysicalPrimary(track->GetLabel())){
1830                                 fMCtrack = fMCstack->Particle(track->GetLabel());
1831
1832                                 Int_t pdg = fMCtrack->GetPdgCode();
1833                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax ){
1834
1835                                         if(fMCtrack->GetFirstMother()>0){
1836                                             fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1837                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
1838                                                         if(fIsHFE1) fPtMC_TPC_All->Fill(fMCtrack->Pt());
1839                                         }
1840                                         }
1841                                 }
1842                         }
1843                 }
1844
1845
1846                         if(fPt>1 && fPt<2) fTOF01->Fill(fTOFnSigma,fTPCnSigma);
1847                         if(fPt>2 && fPt<4) fTOF02->Fill(fTOFnSigma,fTPCnSigma);
1848                         if(fPt>4 && fPt<6) fTOF03->Fill(fTOFnSigma,fTPCnSigma);
1849
1850                         ///________________________________________________________________________
1851                         ///PID
1852                         ///Here the PID cuts defined in the file "ConfigEMCalHFEpA.C" is applied
1853             Int_t pidpassed = 1;
1854             AliHFEpidObject hfetrack;
1855             hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1856             hfetrack.SetRecTrack(track);
1857             hfetrack.SetPP();   //proton-proton analysis
1858             if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
1859             fpid->Fill(pidpassed);
1860
1861             if(pidpassed==0) continue;
1862                         ///________________________________________________________________________
1863
1864
1865                         ////////////////////////////////////////////////////////////////////
1866                         ///TPC efficiency calculations
1867
1868                         /// changing start here
1869                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
1870                 {
1871                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
1872                         Int_t pdg = fMCparticle->GetPdgCode();
1873
1874                                 //
1875                         if(fMCparticle->IsPhysicalPrimary()){
1876
1877
1878                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
1879
1880                                         Bool_t MotherFound = FindMother(track->GetLabel());
1881                                         if(MotherFound){
1882                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1883                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
1884                                                         if(fIsHFE1) fPtMC_TPC_Selected->Fill(fMCparticle->Pt());
1885                                         }
1886                                         }
1887                                 }
1888                         }
1889                 }///until here
1890
1891                 else if(fIsMC && track->GetLabel()>=0)//ESD
1892                 {
1893
1894                         if(fMCstack->IsPhysicalPrimary(track->GetLabel())){
1895                                 fMCtrack = fMCstack->Particle(track->GetLabel());
1896
1897                                 Int_t pdg = fMCtrack->GetPdgCode();
1898                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax ){
1899
1900                                         if(fMCtrack->GetFirstMother()>0){
1901                                             fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1902                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
1903                                                         if(fIsHFE1) fPtMC_TPC_Selected->Fill(fMCtrack->Pt());
1904                                         }
1905                                         }
1906                                 }
1907                         }
1908                 }
1909
1910                         //Eta Cut for TPC only
1911                 if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
1912                         fTPCnsigma_pt_2D->Fill(fPt,fTPCnSigma);
1913                 }
1914
1915                         //Background for TPC only
1916                 if(fFillBackground){
1917                         Background(track, iTracks, Vtrack, kTRUE); //IsTPConly=kTRUE
1918                 }
1919
1920
1921                 fTPCnsigma_p[2]->Fill(fP,fTPCnSigma);
1922                 fTPC_p[2]->Fill(fP,fTPCsignal);
1923                 TPCNcls = track->GetTPCNcls();
1924                 Float_t pos3[3]={0,0,0};
1925
1926                 if(track->GetEMCALcluster()>0)
1927                 {
1928                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
1929                         if(fClus->IsEMCAL())
1930                         {
1931
1932                                         //________________________________________________________________________
1933
1934
1935                                         //Cluster timing distribution
1936                                 if(fUseEMCal && !fIsAOD ){
1937                                         AliESDCaloCells &cells=*(fESD->GetEMCALCells());
1938                                                 //      Int_t nTotalCells = cells.GetNumberOfCells() ;
1939                                                 //Int_t type        = cells.GetType();
1940                                                 //for (Int_t icell=  0; icell <  nTotalCells; icell++) {
1941                                                 //fTime->Fill(cells.GetTime(icell));
1942                                                 //}
1943
1944                                         TRefArray* caloClusters = new TRefArray();
1945                                         fESD->GetEMCALClusters(caloClusters);
1946
1947                                         Int_t nclus = caloClusters->GetEntries();
1948
1949                                                 //fClusESD = fESD->GetCaloCluster(track->GetEMCALcluster());
1950
1951                                         for (Int_t icl = 0; icl < nclus; icl++) {
1952
1953                                                 AliESDCaloCluster* clus = (AliESDCaloCluster*)caloClusters->At(icl);
1954
1955                                                 if(fClus->IsEMCAL()){
1956                                                         Float_t ncells  = fClus->GetNCells();
1957                                                         UShort_t * index = clus->GetCellsAbsId() ;
1958                                                         UShort_t * index2 = fClus->GetCellsAbsId() ;
1959
1960
1961                                                         for(Int_t i = 0; i < ncells ; i++){
1962
1963                                                                 Int_t absId =   index[i];
1964                                                                 fTime->Fill(fPt,cells.GetCellTime(absId));
1965
1966                                                                 Int_t absId2 =   index2[i];
1967                                                                 fTime2->Fill(fPt,cells.GetCellTime(absId2));
1968                                                         }
1969
1970                                                 }
1971                                         }
1972
1973
1974
1975
1976                                 }
1977
1978                                 if(fUseEMCal){
1979                                         double emctof = fClus->GetTOF();
1980                                         ftimingEle->Fill(fPt,emctof);
1981                                 }
1982                                         //________________________________________________________________________
1983
1984
1985
1986
1987                                         /////////////// Residuals
1988                                 Double_t Dx = fClus->GetTrackDx();
1989                                 Double_t Dz = fClus->GetTrackDz();
1990                                 Double_t R=TMath::Sqrt(Dx*Dx+Dz*Dz);
1991                                 for(Int_t i = 0; i < 6; i++)
1992                                 {
1993                                         if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
1994                                         {
1995
1996                                                 fEta[i]->Fill(Dz);
1997                                                 fPhi[i]->Fill(Dx);
1998                                                 fR[i]->Fill(R);
1999                                         }
2000                                 }
2001
2002                                 if(TMath::Abs(fClus->GetTrackDx())<=fdPhiCut && TMath::Abs(fClus->GetTrackDz())<=fdEtaCut)
2003                                 {
2004                                         Float_t Energy  = fClus->E();
2005                                         Float_t EoverP  = Energy/track->P();
2006                                         Float_t M02     = fClus->GetM02();
2007                                         Float_t M20     = fClus->GetM20();
2008                                         Float_t ncells  = fClus->GetNCells();
2009                                                 //----------------------------------------------------------------------------------------
2010                                                 // EtaCut electrons histogram
2011                                                 //Shower Shape Cut
2012                                         if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
2013
2014                                                 if(fUseShowerShapeCut){
2015                                                         if(M02 >= fM02CutMin && M02<=fM02CutMax && M20>=fM20CutMin && M20<=fM20CutMax){
2016                                                                 fEoverP_pt[2]->Fill(fPt,(fClus->E() / fP));
2017                                                                 fShowerShapeCut->Fill(M02,M20);
2018
2019                                                         }
2020
2021                                                 }
2022                                                 if(!fUseShowerShapeCut){
2023                                                         fEoverP_pt[2]->Fill(fPt,(fClus->E() / fP));
2024                                                         fShowerShapeCut->Fill(M02,M20);
2025
2026                                                 }
2027                                                 if(fUseEMCal) fShowerShape_ele->Fill(M02,M20);
2028
2029                                                         //for shower shape cut studies - now with TPC PID Cut
2030                                                 if(fUseEMCal){
2031                                                         fShowerShapeM02_EoverP->Fill(M02,EoverP);
2032                                                         fShowerShapeM20_EoverP->Fill(M20,EoverP);
2033                                                 }
2034
2035                                         }
2036
2037                                                 //----------------------------------------------------------------------------------------
2038
2039
2040
2041                                                 /////////////// for Eta Phi distribution
2042                                         fClus->GetPosition(pos3);
2043                                         TVector3 vpos(pos3[0],pos3[1],pos3[2]);
2044                                         Double_t cphi = vpos.Phi();
2045                                         Double_t ceta = vpos.Eta();
2046                                         fEtaPhi[2]->Fill(cphi,ceta);
2047
2048
2049
2050                                         fTPCNcls_EoverP[2]->Fill(TPCNcls, EoverP);
2051
2052                                         for(Int_t i = 0; i < 6; i++)
2053                                         {
2054                                                 if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
2055                                                 {
2056
2057                                                         fR_EoverP[i]->Fill(R, EoverP);
2058                                                         fNcells[i]->Fill(ncells);
2059                                                         fNcells_EoverP[i]->Fill(EoverP, ncells);
2060                                                         fM02_EoverP[i]->Fill(M02,EoverP);
2061                                                         fM20_EoverP[i]->Fill(M20,EoverP);
2062                                                         fECluster_ptbins[i]->Fill(Energy);
2063                                                         fEoverP_ptbins[i]->Fill(EoverP);
2064
2065                                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax) {
2066                                                                 fNcells_electrons[i]->Fill(ncells);
2067                                                         }
2068
2069                                                         if(M02<0.5 && M20<0.3) {
2070                                                                 fEoverP_wSSCut[i]->Fill(EoverP);
2071                                                         }
2072                                                 }
2073                                         }
2074
2075                                         fNcells_pt->Fill(fPt, ncells);
2076
2077
2078                                                 ////////////////////////////////////////////////////////////////////
2079                                                 ///EMCal - Efficiency calculations
2080
2081                                                 /// changing start here
2082                                         if(fIsMC && fIsAOD && track->GetLabel()>=0)
2083                                         {
2084                                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2085                                                 Int_t pdg = fMCparticle->GetPdgCode();
2086
2087                                                         //
2088                                                 if(fMCparticle->IsPhysicalPrimary()){
2089
2090
2091                                                         if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
2092
2093                                                                 Bool_t MotherFound = FindMother(track->GetLabel());
2094                                                                 if(MotherFound){
2095                                                                         fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2096                                                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
2097                                                                                 fPtMC_EMCal_All->Fill(fMCparticle->Pt());
2098                                                                         }
2099                                                                 }
2100                                                         }
2101                                                 }
2102                                         }///until here
2103
2104                                         else if(fIsMC && track->GetLabel()>=0)//ESD
2105                                         {
2106
2107                                                 if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
2108                                                 {
2109
2110                                                         fMCtrack = fMCstack->Particle(track->GetLabel());
2111
2112                                                         Int_t pdg = fMCtrack->GetPdgCode();
2113                                                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax && fMCtrack->Phi()>=(TMath::Pi()*80/180) && fMCtrack->Phi()<=TMath::Pi())
2114                                                         {
2115                                                                 if(fMCtrack->GetFirstMother()>0){
2116                                                                         fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
2117                                                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
2118
2119                                                                                 fPtMC_EMCal_All->Fill(fMCtrack->Pt());
2120                                                                         }
2121                                                                 }
2122                                                         }
2123                                                 }
2124                                         }
2125
2126
2127                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax)
2128                                         {
2129
2130                                             fECluster[2]->Fill(Energy);
2131                                                         //_______________________________________________________
2132                                                         //Correlation Analysis
2133                                                 if(fUseEMCal)
2134                                                 {
2135                                                         fPtElec_Inc->Fill(fPt);
2136                                                                 //new function to fill non-HFE histos
2137                                                         if(fFillBackground){
2138                                                                 Background(track, iTracks, Vtrack, kFALSE);
2139                                                         }
2140
2141                                                         double emctof2 = fClus->GetTOF();
2142                                                         ftimingEle2->Fill(fPt,emctof2);
2143
2144                                                         if(fCorrelationFlag)
2145                                                         {
2146                                                                 ElectronHadronCorrelation(track, iTracks, Vtrack);
2147                                                         }
2148                                                 }
2149                                                         //_______________________________________________________
2150
2151                                                         ////////////////////////////////////////////////////////////////////
2152                                                         ///EMCal - efficiency calculations
2153
2154                                                 if(track->Charge()<0)  fCharge_n->Fill(fPt);
2155                                                 if(track->Charge()>0)  fCharge_p->Fill(fPt);
2156
2157
2158                                                         /// changing start here
2159                                                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
2160                                                 {
2161                                                         if(track->Charge()<0)  fCharge_n->Fill(fPt);
2162                                                         if(track->Charge()>0)  fCharge_p->Fill(fPt);
2163
2164                                                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2165                                                         fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2166                                                         Int_t pdg = fMCparticle->GetPdgCode();
2167                                                         Int_t mpdg = fMCparticleMother->GetPdgCode();
2168
2169
2170
2171                                                         if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
2172
2173                                                                 if( TMath::Abs(pdg) == 11 ){
2174                                                                         if(TMath::Abs(mpdg) == 221 || TMath::Abs(mpdg) == 22 || TMath::Abs(mpdg) == 111){
2175                                                                                 fPtMCelectronAfterAll_nonPrimary->Fill(fMCparticle->Pt()); //numerator for the total efficiency, non Primary track
2176                                                                         }
2177                                                                 }
2178                                                                 if( TMath::Abs(pdg) == 11 && fMCparticle->IsPhysicalPrimary()) fPtMCelectronAfterAll_Primary->Fill(fMCparticle->Pt());
2179                                                         }
2180
2181
2182                                                         //
2183                                                         if(fMCparticle->IsPhysicalPrimary()){
2184
2185
2186                                                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
2187                                                                                 //check
2188                                                                         if(fIsHFE1) fPtMCelectronAfterAll->Fill(fMCparticle->Pt()); //numerator for the total efficiency
2189
2190
2191                                                                         Bool_t MotherFound = FindMother(track->GetLabel());
2192                                                                         if(MotherFound){
2193                                                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2194                                                                                 if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
2195                                                                                         fPtMC_EMCal_Selected->Fill(fMCparticle->Pt());
2196                                                                                 }
2197                                                                         }
2198                                                                 }
2199                                                         }
2200                                                 }///until here
2201
2202                                                 else if(fIsMC && track->GetLabel()>=0)//ESD
2203                                                 {
2204                                                         if(track->Charge()<0)  fCharge_n->Fill(fPt);
2205                                                         if(track->Charge()>0)  fCharge_p->Fill(fPt);
2206
2207                                                         fMCtrack = fMCstack->Particle(track->GetLabel());
2208                                                         Int_t pdg = fMCtrack->GetPdgCode();
2209
2210                                                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax){
2211                                                                 if( TMath::Abs(pdg) == 11) fPtMCelectronAfterAll_nonPrimary->Fill(fMCtrack->Pt()); //numerator for the total efficiency, non Primary track
2212                                                                 if( TMath::Abs(pdg) == 11 && fMCstack->IsPhysicalPrimary(track->GetLabel())) fPtMCelectronAfterAll_Primary->Fill(fMCtrack->Pt());
2213                                                         }
2214
2215                                                         if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
2216                                                         {
2217                                                                 Bool_t MotherFound = FindMother(track->GetLabel());
2218
2219                                                                 if(MotherFound)
2220                                                                 {
2221                                                                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax){
2222                                                                                 if(fIsHFE1) fPtMCelectronAfterAll->Fill(fMCtrack->Pt()); //numerator for the total efficiency
2223                                                                         }
2224                                                                 }
2225
2226
2227
2228
2229                                                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax && fMCtrack->Phi()>=(TMath::Pi()*80/180) && fMCtrack->Phi()<=TMath::Pi())
2230                                                                 {
2231                                                                         if(fMCtrack->GetFirstMother()>0){
2232                                                                                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
2233                                                                                 if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
2234
2235                                                                                         fPtMC_EMCal_Selected->Fill(fMCtrack->Pt());
2236                                                                                 }
2237                                                                         }
2238                                                                 }
2239                                                         }
2240                                                 }
2241                                                         ///////////////////////////////////////////////////////////////////
2242
2243
2244                                         }
2245                                 }
2246                         }
2247                 }
2248
2249                         //______________________________________________________________
2250                         // Vertex
2251
2252                 fVtxZ[2]->Fill(fZvtx);
2253                 fNTracks[2]->Fill(fNOtrks);
2254                 fNClusters[2]->Fill(ClsNo);
2255
2256                         //______________________________________________________________
2257
2258                         //_______________________________________________________
2259                         //Correlation Analysis
2260                 if(!fUseEMCal)
2261                 {
2262                         fPtElec_Inc->Fill(fPt);
2263
2264                         if(fCorrelationFlag)
2265                         {
2266                                 ElectronHadronCorrelation(track, iTracks, Vtrack);
2267                         }
2268                 }
2269                         //_______________________________________________________
2270
2271                         ///________________________________________________________________________
2272         }
2273
2274                 //__________________________________________________________________
2275                 //Event Mixing Analysis
2276                 //Filling pool
2277         if(fEventMixingFlag)
2278         {
2279                 fPool = fPoolMgr->GetEventPool(fCentrality->GetCentralityPercentile("V0A"), fZvtx); // Get the buffer associated with the current centrality and z-vtx
2280
2281                 if(!fPool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality->GetCentralityPercentile("V0A"), fZvtx));
2282
2283                 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!
2284
2285
2286         }
2287
2288                 //__________________________________________________________________
2289
2290         delete fListOfmotherkink;
2291         PostData(1, fOutputList);
2292 }
2293
2294         //______________________________________________________________________
2295 void AliAnalysisTaskEMCalHFEpA::Terminate(Option_t *)
2296 {
2297                 //Draw result to the screen
2298                 //Called once at the end of the query
2299
2300         fOutputList = dynamic_cast<TList*> (GetOutputData(1));
2301
2302         if(!fOutputList)
2303         {
2304                 printf("ERROR: Output list not available\n");
2305                 return;
2306         }
2307 }
2308
2309         //______________________________________________________________________
2310 Bool_t AliAnalysisTaskEMCalHFEpA::ProcessCutStep(Int_t cutStep, AliVParticle *track)
2311 {
2312                 //Check single track cuts for a given cut step
2313                 //Note this function is called inside the UserExec function
2314         const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
2315         if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
2316         return kTRUE;
2317 }
2318
2319
2320         //______________________________________________________________________
2321
2322
2323 void AliAnalysisTaskEMCalHFEpA::Background(AliVTrack *track, Int_t trackIndex, AliVParticle *vtrack, Bool_t IsTPConly)
2324 {
2325                 ///_________________________________________________________________
2326                 ///MC analysis
2327         if(fIsMC)
2328         {
2329                 if(track->GetLabel() < 0)
2330         {
2331                         AliWarning(Form("The track %d does not have a valid MC label",trackIndex));
2332                         return;
2333         }
2334
2335                 if(fIsAOD)
2336                 {
2337                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2338
2339                         if(fMCparticle->GetMother()<0) return;
2340
2341                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2342
2343                 if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
2344                 {
2345                                         //Is Background
2346                                 if(!IsTPConly)fPtBackgroundBeforeReco->Fill(track->Pt());
2347                                 if(IsTPConly)fPtBackgroundBeforeReco2->Fill(track->Pt());
2348
2349
2350                                 //new 23 September      //weighted histograms //test
2351                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111 ){
2352                                                 //Double_t mPt=fMCparticleMother->Pt();
2353                                                 //Double_t mweight=3*mPt;
2354                                         Double_t mweight=1;
2355                                         if(!IsTPConly)fPtBackgroundBeforeReco_weight->Fill(track->Pt(), mweight);
2356                                         if(IsTPConly)fPtBackgroundBeforeReco2_weight->Fill(track->Pt(), mweight);
2357                                 }
2358                                 else if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111 ){
2359                                                 //Double_t gmPt=fMCparticleGMother->Pt();
2360                                                 //Double_t gmweight=3*gmPt;
2361                                         Double_t gmweight=1;
2362                                         if(!IsTPConly)fPtBackgroundBeforeReco_weight->Fill(track->Pt(), gmweight);
2363                                         if(IsTPConly)fPtBackgroundBeforeReco2_weight->Fill(track->Pt(), gmweight);
2364                                 }
2365                                 else{
2366                                         if(!IsTPConly)fPtBackgroundBeforeReco_weight->Fill(track->Pt());
2367                                         if(IsTPConly)fPtBackgroundBeforeReco2_weight->Fill(track->Pt());
2368                                 }
2369                         }//particle kind
2370                 }//IsAOD
2371                 //ESD
2372                 else
2373                 {
2374                 fMCtrack = fMCstack->Particle(track->GetLabel());
2375
2376                 if(fMCtrack->GetFirstMother()<0) return;
2377
2378                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
2379
2380                 if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
2381                 {
2382                                         //Is Background
2383                                 if(!IsTPConly)fPtBackgroundBeforeReco->Fill(track->Pt());
2384                                 if(IsTPConly)fPtBackgroundBeforeReco2->Fill(track->Pt());
2385                         }
2386                 }
2387         }//IsMC
2388
2389                 ///_________________________________________________________________
2390
2391                 //________________________________________________
2392                 //Associated particle cut
2393         fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
2394     fPartnerCuts->SetRequireITSRefit(kTRUE);
2395     fPartnerCuts->SetRequireTPCRefit(kTRUE);
2396     fPartnerCuts->SetEtaRange(-0.9,0.9);
2397     fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
2398     fPartnerCuts->SetMinNClustersTPC(80);
2399     fPartnerCuts->SetPtRange(0.3,1e10);
2400                 //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
2401                 //fPartnerCuts->SetMaxDCAToVertexXY(1);
2402                 //fPartnerCuts->SetMaxDCAToVertexZ(3);
2403                 //_________________________________________________
2404
2405                 ///#################################################################
2406                 //Non-HFE reconstruction
2407         fNonHFE = new AliSelectNonHFE();
2408         fNonHFE->SetAODanalysis(fIsAOD);
2409         if(fMassCutFlag) fNonHFE->SetInvariantMassCut(fMassCut);
2410         if(fAngleCutFlag) fNonHFE->SetOpeningAngleCut(fAngleCut);
2411         if(fChi2CutFlag) fNonHFE->SetChi2OverNDFCut(fChi2Cut);
2412         if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
2413         fNonHFE->SetAlgorithm("DCA"); //KF
2414         fNonHFE->SetPIDresponse(fPidResponse);
2415         fNonHFE->SetTrackCuts(-3.5,3.5,fPartnerCuts);
2416         fNonHFE->SetAdditionalCuts(fPtMinAsso,fTpcNclsAsso);
2417
2418         if(!IsTPConly){
2419                 fNonHFE->SetHistAngleBack(fOpAngleBack);
2420                 fNonHFE->SetHistAngle(fOpAngle);
2421                 fNonHFE->SetHistDCABack(fDCABack);
2422                 fNonHFE->SetHistDCA(fDCA);
2423                 fNonHFE->SetHistMassBack(fInvMassBack);
2424                 fNonHFE->SetHistMass(fInvMass);
2425         }
2426         if(IsTPConly){
2427                 fNonHFE->SetHistAngleBack(fOpAngleBack2);
2428                 fNonHFE->SetHistAngle(fOpAngle2);
2429                 fNonHFE->SetHistDCABack(fDCABack2);
2430                 fNonHFE->SetHistDCA(fDCA2);
2431                 fNonHFE->SetHistMassBack(fInvMassBack2);
2432                 fNonHFE->SetHistMass(fInvMass2);
2433         }
2434
2435         fNonHFE->FindNonHFE(trackIndex,vtrack,fVevent);
2436
2437
2438
2439                 //Electron Information
2440         Double_t fPhiE = -999;
2441         Double_t fEtaE = -999;
2442         Double_t fPtE = -999;
2443         fPhiE = track->Phi();
2444         fEtaE = track->Eta();
2445         fPtE = track->Pt();
2446
2447                 ///_________________________________________________________________
2448                 ///MC analysis
2449         if(fIsMC)
2450         {
2451                 if(fIsAOD)
2452                 {
2453                         if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
2454                 {
2455
2456                                 Double_t weight=1;
2457
2458                                 if(!IsTPConly){
2459                                         if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
2460                                         if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
2461
2462
2463
2464                                         //new 26 September      //weighted histograms
2465                                         if(TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221){
2466                                                 Double_t mPt=fMCparticleMother->Pt();
2467                                                 Double_t mweight1=1;
2468                                                 Double_t mweight2=1;
2469
2470
2471                                                 //for pions
2472                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111){
2473                                                         Double_t x=mPt;
2474                                                         if(mPt<=4.5) weight=x*x*0.089-0.277*x+1.46;
2475                                                         if(mPt>4.5)  weight=TMath::Erf((x-0.425)/13.05)*5.94;
2476                                                 }
2477                                                 //for eta
2478                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==221){
2479                                                         Double_t x=mPt;
2480                                                         if(mPt<=4.5)  weight=x*x*0.071-0.295*x+1.36;
2481                                                         if(mPt>4.5)  weight=TMath::Erf((x-0.341)/13.31)*4.32;
2482
2483                                                 }
2484
2485                                                         //check this
2486                                                 if(fNonHFE->IsULS()) mweight1=(fNonHFE->GetNULS())/weight;
2487                                                 if(fNonHFE->IsLS())  mweight2=(fNonHFE->GetNLS())/weight;
2488
2489                                                         //fill histos
2490                                                 if(fNonHFE->IsULS())fPtElec_ULS_weight->Fill(fPtE, mweight1);
2491                                                 if(fNonHFE->IsLS())fPtElec_LS_weight->Fill(fPtE, mweight2);
2492                                         }
2493                                         else if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111 || TMath::Abs(fMCparticleGMother->GetPdgCode())==221 ){
2494                                                 Double_t gmPt=fMCparticleGMother->Pt();
2495                                                 Double_t gmweight1=1;
2496                                                 Double_t gmweight2=1;
2497
2498
2499                                                 //for pions
2500                                                 if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111){
2501                                                         Double_t x=gmPt;
2502                                                         if(gmPt<=4.5)  weight=x*x*0.089-0.277*x+1.46;
2503                                                         if(gmPt>4.5)  weight=TMath::Erf((x-0.425)/13.05)*5.94;
2504                                                 }
2505                                                 //for eta
2506                                                 if(TMath::Abs(fMCparticleGMother->GetPdgCode())==221){
2507                                                         Double_t x=gmPt;
2508                                                         if(gmPt<=4.5) weight=x*x*0.071-0.295*x+1.36;
2509                                                         if(gmPt>4.5)  weight=TMath::Erf((x-0.341)/13.31)*4.32;
2510
2511                                                 }
2512
2513                                                         //check this
2514                                                 if(fNonHFE->IsULS()) gmweight1=(fNonHFE->GetNULS())/weight;
2515                                                 if(fNonHFE->IsLS())  gmweight2=(fNonHFE->GetNLS())/weight;
2516
2517                                                 //fill histos
2518                                                 if(fNonHFE->IsULS())fPtElec_ULS_weight->Fill(fPtE, gmweight1);
2519                                                 if(fNonHFE->IsLS())fPtElec_LS_weight->Fill(fPtE, gmweight2);
2520                                         }
2521                                         else{
2522                                                 if(fNonHFE->IsULS()) fPtElec_ULS_weight->Fill(fPtE,fNonHFE->GetNULS());
2523                                                 if(fNonHFE->IsLS()) fPtElec_LS_weight->Fill(fPtE,fNonHFE->GetNLS());
2524                                         }
2525
2526
2527                                 }//!IsTPConly
2528
2529                                 if(IsTPConly){
2530                                         if(fNonHFE->IsULS()) fPtElec_ULS2->Fill(fPtE,fNonHFE->GetNULS());
2531                                         if(fNonHFE->IsLS()) fPtElec_LS2->Fill(fPtE,fNonHFE->GetNLS());
2532
2533
2534
2535
2536                                                 //new 26 September      //weighted histograms
2537                                         if(TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221){
2538                                                 Double_t mPt=fMCparticleMother->Pt();
2539
2540                                                 Double_t mweight1=1;
2541                                                 Double_t mweight2=1;
2542
2543
2544                                                         //for pions
2545                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111){
2546                                                         Double_t x=mPt;
2547                                                         if(mPt<=4.5)  weight=x*x*0.089-0.277*x+1.46;
2548                                                         if(mPt>4.5) weight=TMath::Erf((x-0.425)/13.05)*5.94;
2549                                                 }
2550                                                         //for eta
2551                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==221){
2552                                                         Double_t x=mPt;
2553                                                         if(mPt<=4.5)  weight=x*x*0.071-0.295*x+1.36;
2554                                                         if(mPt>4.5)  weight=TMath::Erf((x-0.341)/13.31)*4.32;
2555
2556                                                 }
2557
2558
2559                                                 //check this
2560                                                 if(fNonHFE->IsULS()) mweight1=(fNonHFE->GetNULS())/weight;
2561                                                 if(fNonHFE->IsLS())  mweight2=(fNonHFE->GetNLS())/weight;
2562
2563                                                         //fill histos
2564                                                 if(fNonHFE->IsULS())fPtElec_ULS2_weight->Fill(fPtE, mweight1);
2565                                                 if(fNonHFE->IsLS())fPtElec_LS2_weight->Fill(fPtE, mweight2);
2566                                         }
2567                                         else if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111 || TMath::Abs(fMCparticleGMother->GetPdgCode())==221 ){
2568                                                 Double_t gmPt=fMCparticleGMother->Pt();
2569                                                 Double_t gmweight1=1;
2570                                                 Double_t gmweight2=1;
2571
2572
2573                                                 //for pions
2574                                                 if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111){
2575                                                         Double_t x=gmPt;
2576                                                         if(gmPt<=4.5)  weight=x*x*0.089-0.277*x+1.46;
2577                                                         if(gmPt>4.5)  weight=TMath::Erf((x-0.425)/13.05)*5.94;
2578                                                 }
2579                                                         //for eta
2580                                                 if(TMath::Abs(fMCparticleGMother->GetPdgCode())==221){
2581                                                         Double_t x=gmPt;
2582                                                         if(gmPt<=4.5)  weight=x*x*0.071-0.295*x+1.36;
2583                                                         if(gmPt>4.5)  weight=TMath::Erf((x-0.341)/13.31)*4.32;
2584
2585                                                 }
2586
2587                                                         //check this
2588                                                 if(fNonHFE->IsULS()) gmweight1=(fNonHFE->GetNULS())/weight;
2589                                                 if(fNonHFE->IsLS())  gmweight2=(fNonHFE->GetNLS())/weight;
2590
2591                                                         //fill histos
2592                                                 if(fNonHFE->IsULS())fPtElec_ULS2_weight->Fill(fPtE, gmweight1);
2593                                                 if(fNonHFE->IsLS())fPtElec_LS2_weight->Fill(fPtE, gmweight2);
2594                                         }
2595                                         else{
2596                                                 if(fNonHFE->IsULS()) fPtElec_ULS2_weight->Fill(fPtE,fNonHFE->GetNULS());
2597                                                 if(fNonHFE->IsLS()) fPtElec_LS2_weight->Fill(fPtE,fNonHFE->GetNLS());
2598                                         }
2599
2600                                 }//IsTPConly
2601
2602                         }//particle kind
2603                 }//close IsAOD
2604                  //It is ESD
2605                 else
2606                 {
2607                         if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
2608                         {
2609                                 if(!IsTPConly){
2610                                         if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
2611                                         if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
2612                                 }
2613
2614                                 if(IsTPConly){
2615                                         if(fNonHFE->IsULS()) fPtElec_ULS2->Fill(fPtE,fNonHFE->GetNULS());
2616                                         if(fNonHFE->IsLS()) fPtElec_LS2->Fill(fPtE,fNonHFE->GetNLS());
2617                                 }
2618                         }
2619                 }
2620         }//close IsMC
2621         ///_________________________________________________________________
2622         //not MC
2623         else
2624         {
2625                 if(!IsTPConly){
2626                         if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
2627                         if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
2628                 }
2629
2630                 if(IsTPConly){
2631                         if(fNonHFE->IsULS()) fPtElec_ULS2->Fill(fPtE,fNonHFE->GetNULS());
2632                         if(fNonHFE->IsLS()) fPtElec_LS2->Fill(fPtE,fNonHFE->GetNLS());
2633                 }
2634         }
2635
2636
2637
2638 }
2639
2640
2641         //______________________________________________________________________
2642 void AliAnalysisTaskEMCalHFEpA::ElectronHadronCorrelation(AliVTrack *track, Int_t trackIndex, AliVParticle *vtrack)
2643 {
2644
2645                 ///_________________________________________________________________
2646                 ///MC analysis
2647         if(fIsMC)
2648         {
2649                 if(track->GetLabel() < 0)
2650         {
2651                         AliWarning(Form("The track %d does not have a valid MC label",trackIndex));
2652                         return;
2653         }
2654
2655                 if(fIsAOD)
2656                 {
2657                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2658
2659                         if(fMCparticle->GetMother()<0) return;
2660
2661                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2662
2663                 if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
2664                 {
2665                                         //Is Background
2666                                 fPtBackgroundBeforeReco->Fill(track->Pt());
2667                         }
2668                 }
2669                 else
2670                 {
2671                 fMCtrack = fMCstack->Particle(track->GetLabel());
2672
2673                 if(fMCtrack->GetFirstMother()<0) return;
2674
2675                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
2676
2677                 if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
2678                 {
2679                                         //Is Background
2680                                 fPtBackgroundBeforeReco->Fill(track->Pt());
2681                         }
2682                 }
2683         }
2684                 ///_________________________________________________________________
2685
2686                 //________________________________________________
2687                 //Associated particle cut
2688         fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
2689     fPartnerCuts->SetRequireITSRefit(kTRUE);
2690     fPartnerCuts->SetRequireTPCRefit(kTRUE);
2691     fPartnerCuts->SetEtaRange(-0.9,0.9);
2692     fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
2693     fPartnerCuts->SetMinNClustersTPC(80);
2694     fPartnerCuts->SetPtRange(0.3,1e10);
2695                 //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
2696                 //fPartnerCuts->SetMaxDCAToVertexXY(1);
2697                 //fPartnerCuts->SetMaxDCAToVertexZ(3);
2698                 //_________________________________________________
2699
2700                 ///#################################################################
2701                 //Non-HFE reconstruction
2702         fNonHFE = new AliSelectNonHFE();
2703         fNonHFE->SetAODanalysis(fIsAOD);
2704         if(fMassCutFlag) fNonHFE->SetInvariantMassCut(fMassCut);
2705         if(fAngleCutFlag) fNonHFE->SetOpeningAngleCut(fAngleCut);
2706         if(fChi2CutFlag) fNonHFE->SetChi2OverNDFCut(fChi2Cut);
2707         if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
2708         fNonHFE->SetAlgorithm("DCA"); //KF
2709         fNonHFE->SetPIDresponse(fPidResponse);
2710         fNonHFE->SetTrackCuts(-3.5,3.5,fPartnerCuts);
2711         fNonHFE->SetAdditionalCuts(fPtMinAsso,fTpcNclsAsso);
2712
2713
2714         fNonHFE->SetHistAngleBack(fOpAngleBack);
2715         fNonHFE->SetHistAngle(fOpAngle);
2716         fNonHFE->SetHistDCABack(fDCABack);
2717         fNonHFE->SetHistDCA(fDCA);
2718         fNonHFE->SetHistMassBack(fInvMassBack);
2719         fNonHFE->SetHistMass(fInvMass);
2720
2721         fNonHFE->FindNonHFE(trackIndex,vtrack,fVevent);
2722
2723         Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
2724         Int_t *fLsPartner = fNonHFE->GetPartnersLS();
2725         Bool_t fUlsIsPartner = kFALSE;
2726         Bool_t fLsIsPartner = kFALSE;
2727                 ///#################################################################
2728
2729
2730                 //Electron Information
2731         Double_t fPhiE = -999;
2732         Double_t fEtaE = -999;
2733         Double_t fPhiH = -999;
2734         Double_t fEtaH = -999;
2735         Double_t fDphi = -999;
2736         Double_t fDeta = -999;
2737         Double_t fPtE = -999;
2738         Double_t fPtH = -999;
2739
2740         Double_t pi = TMath::Pi();
2741
2742         fPhiE = track->Phi();
2743         fEtaE = track->Eta();
2744         fPtE = track->Pt();
2745
2746
2747                 ///_________________________________________________________________
2748                 ///MC analysis
2749         if(fIsMC)
2750         {
2751                 if(fIsAOD)
2752                 {
2753                         if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
2754                 {
2755                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
2756                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
2757                         }
2758                 }
2759                 else
2760                 {
2761                         if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
2762                         {
2763                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
2764                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
2765                         }
2766                 }
2767         }
2768                 ///_________________________________________________________________
2769         else
2770         {
2771                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
2772                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
2773         }
2774
2775
2776
2777
2778                 //__________________________________________________________________
2779                 //Event Mixing Analysis - Hadron Loop
2780                 //Retrieve
2781         if(fEventMixingFlag)
2782         {
2783                 fPool = fPoolMgr->GetEventPool(fCentrality->GetCentralityPercentile("V0A"), fZvtx); // Get the buffer associated with the current centrality and z-vtx
2784
2785                 if(!fPool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f",fCentrality->GetCentralityPercentile("V0A"), fZvtx));
2786
2787                 if(fPool->GetCurrentNEvents() >= 5) // start mixing when 5 events are in the buffer
2788                 {
2789                         fPoolNevents->Fill(fPool->GetCurrentNEvents());
2790
2791                         for (Int_t jMix = 0; jMix < fPool->GetCurrentNEvents(); jMix++)  // mix with each event in the buffer
2792                         {
2793                                 TObjArray* bgTracks = fPool->GetEvent(jMix);
2794
2795                                 for (Int_t kMix = 0; kMix < bgTracks->GetEntriesFast(); kMix++)  // mix with each track in the event
2796                                 {
2797                                         const AliEHCParticle* MixedTrack(dynamic_cast<AliEHCParticle*>(bgTracks->At(kMix)));
2798                                         if (NULL == MixedTrack) continue;
2799
2800                                         fPhiH = MixedTrack->Phi();
2801                                         fEtaH = MixedTrack->Eta();
2802                                         fPtH = MixedTrack->Pt();
2803
2804                                         if(fPtH<0.5 || fPtH>2.0) continue;
2805
2806                                         fDphi = fPhiE - fPhiH;
2807
2808                                         if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
2809                                         if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
2810
2811                                         fDeta = fEtaE - fEtaH;
2812
2813                                         Double_t fPtBin[7] = {1,2,4,6,8,10,15};
2814
2815                                         for(Int_t i = 0; i < 6; i++)
2816                                         {
2817                                             if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
2818                                             {
2819                                                         fCEtaPhi_Inc_EM[i]->Fill(fDphi,fDeta);
2820
2821                                                         if(fNonHFE->IsULS()) fCEtaPhi_ULS_EM[i]->Fill(fDphi,fDeta);
2822                                                         if(fNonHFE->IsLS()) fCEtaPhi_LS_EM[i]->Fill(fDphi,fDeta);
2823
2824                                                         if(fNonHFE->IsULS()) fCEtaPhi_ULS_Weight_EM[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
2825                                                         if(fNonHFE->IsLS()) fCEtaPhi_LS_Weight_EM[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
2826                                             }
2827                                         }
2828
2829                                                 // TODO your code: do event mixing with current event and bgTracks
2830                                                 // note that usually the content filled now is weighted by 1 / pool->GetCurrentNEvents()
2831                                 }
2832                         }
2833                 }
2834         }
2835                 //__________________________________________________________________
2836
2837                 //__________________________________________________________________
2838                 //Same Event Analysis - Hadron Loop
2839         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
2840         {
2841                 if(trackIndex==iTracks) continue;
2842
2843                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
2844                 if (!Vtrack2)
2845                 {
2846                         printf("ERROR: Could not receive track %d\n", iTracks);
2847                         continue;
2848                 }
2849
2850                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
2851
2852                 if(fIsAOD)
2853                 {
2854                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
2855                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
2856                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
2857                         if(atrack2->GetTPCNcls() < 80) continue;
2858                         if(fAssocWithSPD && ((!(atrack2->HasPointOnITSLayer(0))) && (!(atrack2->HasPointOnITSLayer(1))))) continue;
2859                 }
2860                 else
2861                 {
2862                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
2863                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
2864                 }
2865
2866                 fPhiH = track2->Phi();
2867                 fEtaH = track2->Eta();
2868                 fPtH = track2->Pt();
2869
2870                 if(fPtH<0.5 || fPtH>2.0) continue;
2871
2872                 fDphi = fPhiE - fPhiH;
2873
2874                 if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
2875                 if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
2876
2877                 fDeta = fEtaE - fEtaH;
2878
2879                 Double_t fPtBin[7] = {1,2,4,6,8,10,15};
2880
2881                         //______________________________________________________________
2882                         //Check if this track is a Non-HFE partner
2883                 for(Int_t i = 0; i < fNonHFE->GetNULS(); i++)
2884                 {
2885                         if(fUlsPartner[i]==iTracks) fUlsIsPartner=kTRUE;
2886                 }
2887                 for(Int_t i = 0; i < fNonHFE->GetNLS(); i++)
2888                 {
2889                         if(fLsPartner[i]==iTracks) fLsIsPartner=kTRUE;
2890                 }
2891                         //______________________________________________________________
2892
2893                 for(Int_t i = 0; i < 6; i++)
2894                 {
2895                     if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
2896                     {
2897                                 fCEtaPhi_Inc[i]->Fill(fDphi,fDeta);
2898
2899                                 if(fNonHFE->IsULS()) fCEtaPhi_ULS[i]->Fill(fDphi,fDeta);
2900                                 if(fNonHFE->IsLS()) fCEtaPhi_LS[i]->Fill(fDphi,fDeta);
2901                                 if(fNonHFE->IsULS() && !fUlsIsPartner) fCEtaPhi_ULS_NoP[i]->Fill(fDphi,fDeta);
2902                                 if(fNonHFE->IsLS() && !fLsIsPartner) fCEtaPhi_LS_NoP[i]->Fill(fDphi,fDeta);
2903
2904                                 if(fNonHFE->IsULS()) fCEtaPhi_ULS_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
2905                                 if(fNonHFE->IsLS()) fCEtaPhi_LS_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
2906                                 if(fNonHFE->IsULS() && !fUlsIsPartner) fCEtaPhi_ULS_NoP_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
2907                                 if(fNonHFE->IsLS() && !fLsIsPartner) fCEtaPhi_LS_NoP_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
2908                     }
2909                 }
2910         }
2911 }
2912
2913         //____________________________________________________________________________________________________________
2914         //Create a TObjArray with selected hadrons, for the mixed event analysis
2915 TObjArray* AliAnalysisTaskEMCalHFEpA::SelectedHadrons()
2916 {
2917         fTracksClone = new TObjArray;
2918         fTracksClone->SetOwner(kTRUE);
2919
2920         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
2921         {
2922                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
2923                 if (!Vtrack2)
2924                 {
2925                         printf("ERROR: Could not receive track %d\n", iTracks);
2926                         continue;
2927                 }
2928
2929                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
2930
2931                 if(fIsAOD)
2932                 {
2933                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
2934                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
2935                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
2936                         if(atrack2->GetTPCNcls() < 80) continue;
2937                         if(fAssocWithSPD && ((!(atrack2->HasPointOnITSLayer(0))) && (!(atrack2->HasPointOnITSLayer(1))))) continue;
2938                 }
2939                 else
2940                 {
2941                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
2942                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
2943                 }
2944
2945                 fTracksClone->Add(new AliEHCParticle(track2->Eta(), track2->Phi(), track2->Pt()));
2946         }
2947         return fTracksClone;
2948 }
2949         //____________________________________________________________________________________________________________
2950
2951         //______________________________________________________________________
2952 void AliAnalysisTaskEMCalHFEpA::DiHadronCorrelation(AliVTrack *track, Int_t trackIndex)
2953 {
2954                 //________________________________________________
2955                 //Associated particle cut
2956         fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
2957     fPartnerCuts->SetRequireITSRefit(kTRUE);
2958     fPartnerCuts->SetRequireTPCRefit(kTRUE);
2959     fPartnerCuts->SetEtaRange(-0.9,0.9);
2960     fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
2961     fPartnerCuts->SetMinNClustersTPC(80);
2962     fPartnerCuts->SetPtRange(0.3,1e10);
2963                 //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
2964                 //fPartnerCuts->SetMaxDCAToVertexXY(1);
2965                 //fPartnerCuts->SetMaxDCAToVertexZ(3);
2966                 //_________________________________________________
2967
2968                 //Electron Information
2969         Double_t fPhiE = -999;
2970         Double_t fEtaE = -999;
2971         Double_t fPhiH = -999;
2972         Double_t fEtaH = -999;
2973         Double_t fDphi = -999;
2974         Double_t fDeta = -999;
2975         Double_t fPtE = -999;
2976         Double_t fPtH = -999;
2977
2978         Double_t pi = TMath::Pi();
2979
2980         fPhiE = track->Phi();
2981         fEtaE = track->Eta();
2982         fPtE = track->Pt();
2983
2984                 //__________________________________________________________________
2985                 //Same Event Analysis - Hadron Loop
2986         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
2987         {
2988                 if(trackIndex==iTracks) continue;
2989
2990                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
2991                 if (!Vtrack2)
2992                 {
2993                         printf("ERROR: Could not receive track %d\n", iTracks);
2994                         continue;
2995                 }
2996
2997                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
2998
2999                 if(fIsAOD)
3000                 {
3001                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
3002                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
3003                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
3004                         if(atrack2->GetTPCNcls() < 80) continue;
3005                         if(fAssocWithSPD && ((!(atrack2->HasPointOnITSLayer(0))) && (!(atrack2->HasPointOnITSLayer(1))))) continue;
3006                 }
3007                 else
3008                 {
3009                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
3010                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
3011                 }
3012
3013                 fPhiH = track2->Phi();
3014                 fEtaH = track2->Eta();
3015                 fPtH = track2->Pt();
3016
3017                 if(fPtH<0.5 || fPtH>2.0) continue;
3018
3019                 fDphi = fPhiE - fPhiH;
3020
3021                 if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
3022                 if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
3023
3024                 fDeta = fEtaE - fEtaH;
3025
3026                 Double_t fPtBin[7] = {1,2,4,6,8,10,15};
3027
3028                 for(Int_t i = 0; i < 6; i++)
3029                 {
3030                     if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
3031                     {
3032                                 fCEtaPhi_Inc_DiHadron[i]->Fill(fDphi,fDeta);
3033                     }
3034                 }
3035         }
3036 }
3037         //____________________________________________________________________________________________________________
3038
3039         //______________________________________________________________________
3040 Bool_t AliAnalysisTaskEMCalHFEpA::FindMother(Int_t mcIndex)
3041 {
3042         fIsHFE1 = kFALSE;
3043         fIsHFE2 = kFALSE;
3044         fIsNonHFE = kFALSE;
3045         fIsFromD = kFALSE;
3046         fIsFromB = kFALSE;
3047         fIsFromPi0 = kFALSE;
3048         fIsFromEta = kFALSE;
3049         fIsFromGamma = kFALSE;
3050
3051         if(mcIndex < 0 || !fIsMC)
3052         {
3053                 return kFALSE;
3054         }
3055
3056         Int_t pdg = -99999;
3057         Int_t mpdg = -99999;
3058         Int_t gmpdg = -99999;
3059         Int_t ggmpdg = -99999;
3060         Int_t gggmpdg = -99999;
3061
3062         if(fIsAOD)
3063         {
3064                 fMCparticle = (AliAODMCParticle*) fMCarray->At(mcIndex);
3065
3066                 pdg = TMath::Abs(fMCparticle->GetPdgCode());
3067
3068
3069                 if(pdg!=11)
3070                 {
3071                         fIsHFE1 = kFALSE;
3072                         fIsHFE2 = kFALSE;
3073                         fIsNonHFE = kFALSE;
3074                         fIsFromD = kFALSE;
3075                         fIsFromB = kFALSE;
3076                         fIsFromPi0 = kFALSE;
3077                         fIsFromEta = kFALSE;
3078                         fIsFromGamma = kFALSE;
3079                         return kFALSE;
3080                 }
3081
3082                 if(fMCparticle->GetMother()<0)
3083                 {
3084                         fIsHFE1 = kFALSE;
3085                         fIsHFE2 = kFALSE;
3086                         fIsNonHFE = kFALSE;
3087                         fIsFromD = kFALSE;
3088                         fIsFromB = kFALSE;
3089                         fIsFromPi0 = kFALSE;
3090                         fIsFromEta = kFALSE;
3091                         fIsFromGamma = kFALSE;
3092                         return kFALSE;
3093                 }
3094
3095                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
3096                 mpdg = TMath::Abs(fMCparticleMother->GetPdgCode());
3097
3098                 if(fMCparticleMother->GetMother()<0)
3099                 {
3100                         gmpdg = 0;
3101                         ggmpdg = 0;
3102                         gggmpdg = 0;
3103                 }
3104                 else
3105                 {
3106                         fMCparticleGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleMother->GetMother());
3107                         gmpdg = TMath::Abs(fMCparticleGMother->GetPdgCode());
3108                         if(fMCparticleGMother->GetMother()<0)
3109                         {
3110                                 ggmpdg = 0;
3111                                 gggmpdg = 0;
3112                         }
3113                         else
3114                         {
3115                                 fMCparticleGGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleGMother->GetMother());
3116                                 ggmpdg = TMath::Abs(fMCparticleGGMother->GetPdgCode());
3117                                 if(fMCparticleGGMother->GetMother()<0)
3118                                 {
3119                                         gggmpdg = 0;
3120                                 }
3121                                 else
3122                                 {
3123                                         fMCparticleGGGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleGGMother->GetMother());
3124                                         gggmpdg = TMath::Abs(fMCparticleGGGMother->GetPdgCode());
3125                                 }
3126                         }
3127                 }
3128         }
3129         else
3130         {
3131                 fMCtrack = fMCstack->Particle(mcIndex);
3132
3133                 pdg = TMath::Abs(fMCtrack->GetPdgCode());
3134
3135                 if(pdg!=11)
3136                 {
3137                         fIsHFE1 = kFALSE;
3138                         fIsHFE2 = kFALSE;
3139                         fIsNonHFE = kFALSE;
3140                         fIsFromD = kFALSE;
3141                         fIsFromB = kFALSE;
3142                         fIsFromPi0 = kFALSE;
3143                         fIsFromEta = kFALSE;
3144                         fIsFromGamma = kFALSE;
3145                         return kFALSE;
3146                 }
3147
3148                 if(fMCtrack->GetFirstMother()<0)
3149                 {
3150                         fIsHFE1 = kFALSE;
3151                         fIsHFE2 = kFALSE;
3152                         fIsNonHFE = kFALSE;
3153                         fIsFromD = kFALSE;
3154                         fIsFromB = kFALSE;
3155                         fIsFromPi0 = kFALSE;
3156                         fIsFromEta = kFALSE;
3157                         fIsFromGamma = kFALSE;
3158                         return kFALSE;
3159                 }
3160
3161                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
3162                 mpdg = TMath::Abs(fMCtrackMother->GetPdgCode());
3163
3164                 if(fMCtrackMother->GetFirstMother()<0)
3165                 {
3166                         gmpdg = 0;
3167                         ggmpdg = 0;
3168                         gggmpdg = 0;
3169                 }
3170                 else
3171                 {
3172                         fMCtrackGMother = fMCstack->Particle(fMCtrackMother->GetFirstMother());
3173                         gmpdg = TMath::Abs(fMCtrackGMother->GetPdgCode());
3174
3175                         if(fMCtrackGMother->GetFirstMother()<0)
3176                         {
3177                                 ggmpdg = 0;
3178                                 gggmpdg = 0;
3179                         }
3180                         else
3181                         {
3182                                 fMCtrackGGMother = fMCstack->Particle(fMCtrackGMother->GetFirstMother());
3183                                 ggmpdg = TMath::Abs(fMCtrackGGMother->GetPdgCode());
3184
3185                                 if(fMCtrackGGMother->GetFirstMother()<0)
3186                                 {
3187                                         gggmpdg = 0;
3188                                 }
3189                                 else
3190                                 {
3191                                         fMCtrackGGGMother = fMCstack->Particle(fMCtrackGGMother->GetFirstMother());
3192                                         gggmpdg = TMath::Abs(fMCtrackGGGMother->GetPdgCode());
3193                                 }
3194                         }
3195                 }
3196         }
3197
3198         //Tag Electron Source
3199         if(mpdg==111 || mpdg==221 || mpdg==22)
3200         {
3201                 fIsHFE1 = kFALSE;
3202                 fIsHFE2 = kFALSE;
3203                 fIsNonHFE = kTRUE;
3204                 fIsFromD = kFALSE;
3205                 fIsFromB = kFALSE;
3206
3207                 fIsFromPi0 = kFALSE;
3208                 fIsFromEta = kFALSE;
3209                 fIsFromGamma = kFALSE;
3210
3211                 if(mpdg==111) fIsFromPi0 = kFALSE;
3212                 if(mpdg==221)fIsFromEta = kFALSE;
3213                 if(mpdg==22) fIsFromGamma = kFALSE;
3214
3215                 return kTRUE;
3216         }
3217         else
3218         {
3219                 fIsHFE1 = kFALSE;
3220                 fIsHFE2 = kTRUE;
3221
3222                 fIsFromPi0 = kFALSE;
3223                 fIsFromEta = kFALSE;
3224                 fIsFromGamma = kFALSE;
3225
3226                 fIsNonHFE = kFALSE;
3227
3228                 fIsFromD = kFALSE;
3229                 fIsFromB = kFALSE;
3230
3231                 if(mpdg>400 && mpdg<500)
3232                 {
3233                         if((gmpdg>500 && gmpdg<600) || (ggmpdg>500 && ggmpdg<600) || (gggmpdg>500 && gggmpdg<600))
3234                         {
3235                                 fIsHFE1 = kTRUE;
3236                                 fIsFromD = kFALSE;
3237                                 fIsFromB = kTRUE;
3238                                 return kTRUE;
3239                         }
3240                         else
3241                         {
3242                                 fIsHFE1 = kTRUE;
3243                                 fIsFromD = kTRUE;
3244                                 fIsFromB = kFALSE;
3245                                 return kTRUE;
3246                         }
3247                 }
3248                 else if(mpdg>500 && mpdg<600)
3249                 {
3250                         fIsHFE1 = kTRUE;
3251                         fIsFromD = kFALSE;
3252                         fIsFromB = kTRUE;
3253                         return kTRUE;
3254                 }
3255                 else
3256                 {
3257                         fIsHFE1 = kFALSE;
3258                         fIsFromD = kFALSE;
3259                         fIsFromB = kFALSE;
3260                         return kFALSE;
3261                 }
3262         }
3263 }