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