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