PWGHF/hfe/AliAnalysisTaskEMCalHFEpA.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17         ////////////////////////////////////////////////////////////////////////
  18         //                                                                    //
  19         //      Task for Heavy-flavour electron analysis in pPb collisions    //
  20         //      (+ Electron-Hadron Jetlike Azimuthal Correlation)             //
  21         //                                                                                                                                        //
  22         //              version: May 30, 2014.                                                                //
  23         //                                                                    //
  24         //          Authors                                                                               //
  25         //              Elienos Pereira de Oliveira Filho (epereira@cern.ch)          //
  26         //              Cristiane Jahnke                (cristiane.jahnke@cern.ch)                    //
  27         //                                                                    //
  28         ////////////////////////////////////////////////////////////////////////
  29
  30 #include "TChain.h"
  31 #include "TTree.h"
  32 #include "TNtuple.h"
  33 #include "TH1F.h"
  34 #include "TH2F.h"
  35 #include "TCanvas.h"
  36 #include "AliAnalysisTask.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliESDEvent.h"
  39 #include "AliAODEvent.h"
  40 #include "AliVEvent.h"
  41 #include "AliESDInputHandler.h"
  42 #include "AliESDtrackCuts.h"
  43 #include "AliESDCaloCluster.h"
  44 #include "AliESDCaloCells.h"
  45 #include "AliEMCALTrack.h"
  46 #include "AliExternalTrackParam.h"
  47 #include "AliPhysicsSelection.h"
  48 #include "TGeoGlobalMagField.h"
  49 #include "AliMagF.h"
  50 #include "AliLog.h"
  51 #include "AliStack.h"
  52 #include "AliCentrality.h"
  53 #include "AliAODMCParticle.h"
  54 #include "AliAODMCHeader.h"
  55 #include "AliPID.h"
  56 #include "AliPIDResponse.h"
  57 #include "AliHFEcontainer.h"
  58 #include "AliHFEcuts.h"
  59 #include "AliHFEpid.h"
  60 #include "AliHFEpidBase.h"
  61 #include "AliHFEpidQAmanager.h"
  62 #include "AliHFEtools.h"
  63 #include "AliCFContainer.h"
  64 #include "AliCFManager.h"
  65 #include "AliSelectNonHFE.h"
  66 #include "AliHFEpidTPC.h"
  67 #include "AliAnalysisTaskEMCalHFEpA.h"
  68 #include "TMath.h"
  69 #include "THnSparse.h"
  70 #include "TLorentzVector.h"
  71 #include "TString.h"
  72 #include "TFile.h"
  73 #include "AliESDHandler.h"
  74 #include "AliMCEventHandler.h"
  75 #include "AliMCEvent.h"
  76 #include "AliStack.h"
  77 #include "TParticle.h"
  78 #include "AliLog.h"
  79 #include "AliAnalysisTaskSE.h"
  80 #include "TRefArray.h"
  81 #include "TVector.h"
  82 #include "stdio.h"
  83 #include "TGeoManager.h"
  84 #include "iostream"
  85 #include "fstream"
  86 #include "AliKFParticle.h"
  87 #include "AliKFVertex.h"
  88 #include "AliVParticle.h"
  89 #include "AliVTrack.h"
  90 #include "AliEventPoolManager.h"
  91 #include "TObjArray.h"
  92         //include to use reader as Lucile does
  93 #include "AliCaloTrackAODReader.h"
  94 #include "AliCaloTrackReader.h"
  95 #include "AliEMCALRecoUtils.h" //to remove exotics
  96 #include "AliAODHeader.h"
  97 #include "AliEMCALGeometry.h"
  98
  99
 100
 101         // --- ANALYSIS system ---
 102 #include "AliCalorimeterUtils.h"
 103 #include "AliESDEvent.h"
 104 #include "AliMCEvent.h"
 105 #include "AliStack.h"
 106 #include "AliAODPWG4Particle.h"
 107 #include "AliVCluster.h"
 108 #include "AliVCaloCells.h"
 109 #include "AliMixedEvent.h"
 110 #include "AliAODCaloCluster.h"
 111 #include "AliOADBContainer.h"
 112 #include "AliAnalysisManager.h"
 113
 114         // --- Detector ---
 115 #include "AliEMCALGeometry.h"
 116 #include "AliPHOSGeoUtils.h"
 117
 118         //______________________________________________________________________
 119
 120         //______________________________________________________________________
 121 ClassImp(AliAnalysisTaskEMCalHFEpA)
 122
 123         //______________________________________________________________________
 124 AliAnalysisTaskEMCalHFEpA::AliAnalysisTaskEMCalHFEpA(const char *name)
 125 : AliAnalysisTaskSE(name)
 126 ,fCorrelationFlag(0)
 127 ,fIsMC(0)
 128 ,fUseEMCal(kFALSE)
 129
 130 ,fUseTrigger(kFALSE)
 131 ,fUseShowerShapeCut(kFALSE)
 132 ,fFillBackground(kFALSE)
 133 ,fAssocWithSPD(kFALSE)
 134
 135 ,fEMCEG1(kFALSE)
 136 ,fEMCEG2(kFALSE)
 137 ,fIsHFE1(kFALSE)
 138 ,fIsHFE2(kFALSE)
 139 ,fIsNonHFE(kFALSE)
 140 ,fIsFromD(kFALSE)
 141 ,fIsFromB(kFALSE)
 142 ,fIsFromPi0(kFALSE)
 143 ,fIsFromEta(kFALSE)
 144 ,fIsFromGamma(kFALSE)
 145 ,fESD(0)
 146 ,fAOD(0)
 147 ,fVevent(0)
 148 ,fPartnerCuts(new AliESDtrackCuts())
 149 ,fOutputList(0)
 150 ,fPidResponse(0)
 151 ,fNonHFE(new AliSelectNonHFE())
 152 ,fIsAOD(kFALSE)
 153 ,fCentrality(0)
 154 ,fCentralityMin(0)
 155 ,fCentralityMax(100)
 156 ,fHasCentralitySelection(kFALSE)
 157 ,fCentralityHist(0)
 158 ,fCentralityHistPass(0)
 159 ,fZvtx(0)
 160 ,fEstimator(0)
 161 ,fClus(0)
 162 //,fClusESD(0)
 163 ,fNevent(0)
 164 ,fNevent2(0)
 165 ,fPtElec_Inc(0)
 166 ,fPtPrim(0)
 167 ,fPtSec(0)
 168 ,fPtPrim2(0)
 169 ,fPtSec2(0)
 170 ,fCharge_n(0)
 171 ,fCharge_p(0)
 172 ,fTime(0)
 173 ,fTime2(0)
 174 ,ftimingEle(0)
 175 ,ftimingEle2(0)
 176 ,fPtElec_ULS(0)
 177 ,fPtElec_LS(0)
 178 ,fPtElec_ULS_NoPid(0)
 179 ,fPtElec_LS_NoPid(0)
 180 ,fPtElec_ULS_MC(0)
 181 ,fPtElec_ULS_MC_weight(0)
 182 ,fPtElec_ULS2(0)
 183 ,fPtElec_LS2(0)
 184 ,fPtElec_ULS_weight(0)
 185 ,fPtElec_LS_weight(0)
 186 ,fPtElec_ULS2_weight(0)
 187 ,fPtElec_LS2_weight(0)
 188 ,fTOF01(0)
 189 ,fTOF02(0)
 190 ,fTOF03(0)
 191 ,fpid(0)
 192 ,fEoverP_pt(0)
 193 ,fEoverP_tpc(0)
 194 ,fTPC_pt(0)
 195 ,fTPC_p(0)
 196 ,fTPCnsigma_pt(0)
 197 ,fTPCnsigma_p(0)
 198 ,fTPCnsigma_pt_2D(0)
 199 ,fShowerShapeCut(0)
 200 ,fShowerShapeM02_EoverP(0)
 201 ,fShowerShapeM20_EoverP(0)
 202 ,fShowerShape_ha(0)
 203 ,fShowerShape_ele(0)
 204 ,fTPCnsigma_eta(0)
 205 ,fTPCnsigma_phi(0)
 206 ,fECluster(0)
 207 ,fECluster_pure(0)
 208 ,fECluster_not_exotic(0)
 209 ,fECluster_not_exotic1(0)
 210 ,fECluster_not_exotic2(0)
 211 ,fECluster_exotic(0)
 212 ,fNCluster_pure(0)
 213 ,fNCluster_pure_aod(0)
 214 ,fNCluster_ECluster(0)
 215 ,fNcells_energy(0)
 216 ,fNcells_energy_elec_selected(0)
 217 ,fNcells_energy_not_exotic(0)
 218
 219 ,fEtaPhi(0)
 220 ,fEtaPhi_num(0)
 221 ,fEtaPhi_den(0)
 222 ,fEtaPhi_data(0)
 223 ,fpt_reco_pt_MC_num(0)
 224 ,fpt_reco_pt_MC_den(0)
 225 ,fVtxZ(0)
 226
 227 ,fVtxZ_new1(0)
 228 ,fVtxZ_new2(0)
 229 ,fVtxZ_new3(0)
 230 ,fVtxZ_new4(0)
 231
 232 ,fzRes1(0)
 233 ,fzRes2(0)
 234 ,fSPD_track_vtx1(0)
 235 ,fSPD_track_vtx2(0)
 236
 237
 238 ,fEtad(0)
 239 ,fNTracks(0)
 240 ,fTrack_Multi(0)
 241 ,fNTracks_pt(0)
 242 ,fNTracks_eta(0)
 243 ,fNTracks_phi(0)
 244 ,fNClusters(0)
 245 ,fTPCNcls_EoverP(0)
 246 ,fTPCNcls_pid(0)
 247 ,fEta(0)
 248 ,fPhi(0)
 249 ,fR(0)
 250 ,fR_EoverP(0)
 251 ,fNcells(0)
 252 ,fNcells_EoverP(0)
 253 ,fNcells_electrons(0)
 254 ,fNcells_hadrons(0)
 255 ,fECluster_ptbins(0)
 256 ,fEoverP_ptbins(0)
 257 ,fEoverP_wSSCut(0)
 258 ,fM02_EoverP(0)
 259 ,fM20_EoverP(0)
 260 ,fTPCnsigma_eta_electrons(0)
 261 ,fTPCnsigma_eta_hadrons(0)
 262 ,fEoverP_pt_pions(0)
 263 ,ftpc_p_EoverPcut(0)
 264 ,fnsigma_p_EoverPcut(0)
 265 ,fEoverP_pt_pions2(0)
 266 ,fNcells_pt(0)
 267 ,fEoverP_pt_hadrons(0)
 268 ,fCEtaPhi_Inc(0)
 269 ,fCEtaPhi_ULS(0)
 270 ,fCEtaPhi_LS(0)
 271 ,fCEtaPhi_ULS_NoP(0)
 272 ,fCEtaPhi_LS_NoP(0)
 273 ,fCEtaPhi_ULS_Weight(0)
 274 ,fCEtaPhi_LS_Weight(0)
 275 ,fCEtaPhi_ULS_NoP_Weight(0)
 276 ,fCEtaPhi_LS_NoP_Weight(0)
 277 ,fInvMass(0)
 278 ,fInvMassBack(0)
 279 ,fDCA(0)
 280 ,fDCABack(0)
 281 ,fOpAngle(0)
 282 ,fOpAngleBack(0)
 283 ,fInvMass2(0)
 284 ,fInvMassBack2(0)
 285 ,fDCA2(0)
 286 ,fDCABack2(0)
 287 ,fOpAngle2(0)
 288 ,fOpAngleBack2(0)
 289 ,fMassCut(0.1)
 290 ,fEtaCutMin(-0.9)
 291 ,fEtaCutMax(0.9)
 292 ,fdPhiCut(0.05)
 293 ,fdEtaCut(0.05)
 294 ,fEoverPCutMin(0.8)
 295 ,fEoverPCutMax(1.2)
 296 ,fM20CutMin(0.0)
 297 ,fM20CutMax(10)
 298 ,fM02CutMin(0.0)
 299 ,fM02CutMax(10)
 300 ,fAngleCut(999)
 301 ,fChi2Cut(3.5)
 302 ,fDCAcut(999)
 303 ,fMassCutFlag(kTRUE)
 304 ,fAngleCutFlag(kFALSE)
 305 ,fChi2CutFlag(kFALSE)
 306 ,fDCAcutFlag(kFALSE)
 307 ,fAssHadronPtMin(0.5)
 308 ,fAssHadronPtMax(2.0)
 309 ,fPtBackgroundBeforeReco(0)
 310 ,fPtBackgroundBeforeReco2(0)
 311 ,fPtBackgroundBeforeReco_weight(0)
 312 ,fPtBackgroundBeforeReco2_weight(0)
 313 ,fpT_m_electron(0)
 314 ,fpT_gm_electron(0)
 315 ,fPtBackgroundAfterReco(0)
 316 ,fPtMinAsso(0.3)
 317 ,fTpcNclsAsso(80)
 318 ,fPtMCparticleAll(0)
 319 ,fPtMCparticleAll_nonPrimary(0)
 320 ,fPtMCparticleAlle_nonPrimary(0)
 321 ,fPtMCparticleAlle_Primary(0)
 322 ,fPtMCparticleReco(0)
 323 ,fPtMCparticleReco_nonPrimary(0)
 324 ,fPtMCparticleAllHfe1(0)
 325 ,fPtMCparticleRecoHfe1(0)
 326 ,fPtMCparticleAllHfe2(0)
 327 ,fPtMCparticleRecoHfe2(0)
 328 ,fPtMCelectronAfterAll(0)
 329 ,fPtMCelectronAfterAll_unfolding(0)
 330 ,fPtMCelectronAfterAll_nonPrimary(0)
 331 ,fPtMCelectronAfterAll_Primary(0)
 332 ,fPtMCpi0(0)
 333 ,fPtMCeta(0)
 334 ,fPtMCpi02(0)
 335 ,fPtMCeta2(0)
 336 ,fPtMCpi03(0)
 337 ,fPtMCeta3(0)
 338 ,fPtMC_EMCal_All(0)
 339 ,fPtMC_EMCal_Selected(0)
 340 ,fPtMC_TPC_All(0)
 341 ,fPtMC_TPC_Selected(0)
 342 ,fPt_track_match_den(0)
 343 ,fPt_track_match_num(0)
 344 ,fPtMCWithLabel(0)
 345 ,fPtMCWithoutLabel(0)
 346 ,fPtIsPhysicaPrimary(0)
 347 ,fCuts(0)
 348         //,reader(0)
 349 ,fCFM(0)
 350 ,fPID(new AliHFEpid("hfePid"))
 351         //Lucile
 352 ,fPIDqa(0)
 353 ,fMCstack(0)
 354 ,fRejectKinkMother(kFALSE)
 355 ,fMCtrack(0)
 356 ,fMCtrackMother(0)
 357 ,fMCtrackGMother(0)
 358 ,fMCtrackGGMother(0)
 359 ,fMCtrackGGGMother(0)
 360 ,fMCarray(0)
 361 ,fMCheader(0)
 362 ,fMCparticle(0)
 363 ,fMCparticle2(0)
 364 ,fMCparticleMother(0)
 365 ,fMCparticleGMother(0)
 366 ,fMCparticleGGMother(0)
 367 ,fMCparticleGGGMother(0)
 368 ,fEventHandler(0)
 369 ,fMCevent(0)
 370 ,fPoolMgr(0)
 371 ,fPool(0)
 372 ,fTracksClone(0)
 373 ,fTracks(0)
 374 ,fCEtaPhi_Inc_EM(0)
 375 ,fCEtaPhi_ULS_EM(0)
 376 ,fCEtaPhi_LS_EM(0)
 377 ,fCEtaPhi_ULS_Weight_EM(0)
 378 ,fCEtaPhi_LS_Weight_EM(0)
 379 ,fPoolNevents(0)
 380 ,fEventMixingFlag(0)
 381 ,fCEtaPhi_Inc_DiHadron(0)
 382 ,fPtTrigger_Inc(0)
 383         //,fEMCALRecoUtils(new AliEMCALRecoUtils)
 384         //,fEMCALGeo(0x0)
 385         //,fCaloUtils(0x0)
 386
 387 ,fBitEGA(0)
 388 //,fEMCALRecoUtils(0)//exotic
 389
 390 {
 391                 //Named constructor
 392                 // Define input and output slots here
 393                 // Input slot #0 works with a TChain
 394
 395                 //exotic
 396                 //fEMCALRecoUtils  = new AliEMCALRecoUtils();
 397
 398         DefineInput(0, TChain::Class());
 399                 // Output slot #0 id reserved by the base class for AOD
 400                 // Output slot #1 writes into a TH1 container
 401                 // DefineOutput(1, TH1I::Class());
 402         DefineOutput(1, TList::Class());
 403                 //  DefineOutput(3, TTree::Class());
 404 }
 405
 406         //________________________________________________________________________
 407 AliAnalysisTaskEMCalHFEpA::AliAnalysisTaskEMCalHFEpA()
 408 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCalHFEpA")
 409 ,fCorrelationFlag(0)
 410 ,fIsMC(0)
 411 ,fUseEMCal(kFALSE)
 412 ,fUseTrigger(kFALSE)
 413 ,fUseShowerShapeCut(kFALSE)
 414 ,fFillBackground(kFALSE)
 415 ,fAssocWithSPD(kFALSE)
 416 ,fEMCEG1(kFALSE)
 417 ,fEMCEG2(kFALSE)
 418 ,fIsHFE1(kFALSE)
 419 ,fIsHFE2(kFALSE)
 420 ,fIsNonHFE(kFALSE)
 421 ,fIsFromD(kFALSE)
 422 ,fIsFromB(kFALSE)
 423 ,fIsFromPi0(kFALSE)
 424 ,fIsFromEta(kFALSE)
 425 ,fIsFromGamma(kFALSE)
 426 ,fESD(0)
 427 ,fAOD(0)
 428 ,fVevent(0)
 429 ,fPartnerCuts(new AliESDtrackCuts())
 430 ,fOutputList(0)
 431 ,fPidResponse(0)
 432 ,fNonHFE(new AliSelectNonHFE())
 433 ,fIsAOD(kFALSE)
 434 ,fCentrality(0)
 435 ,fCentralityMin(0)
 436 ,fCentralityMax(100)
 437 ,fHasCentralitySelection(kFALSE)
 438 ,fCentralityHist(0)
 439 ,fCentralityHistPass(0)
 440 ,fZvtx(0)
 441 ,fEstimator(0)
 442 ,fClus(0)
 443         //,fClusESD(0)
 444 ,fNevent(0)
 445 ,fNevent2(0)
 446 ,fPtElec_Inc(0)
 447 ,fPtPrim(0)
 448 ,fPtSec(0)
 449 ,fPtPrim2(0)
 450 ,fPtSec2(0)
 451 ,fCharge_n(0)
 452 ,fCharge_p(0)
 453 ,fTime(0)
 454 ,fTime2(0)
 455 ,ftimingEle(0)
 456 ,ftimingEle2(0)
 457 ,fPtElec_ULS(0)
 458 ,fPtElec_LS(0)
 459 ,fPtElec_ULS_NoPid(0)
 460 ,fPtElec_LS_NoPid(0)
 461 ,fPtElec_ULS_MC(0)
 462 ,fPtElec_ULS_MC_weight(0)
 463 ,fPtElec_ULS2(0)
 464 ,fPtElec_LS2(0)
 465 ,fPtElec_ULS_weight(0)
 466 ,fPtElec_LS_weight(0)
 467 ,fPtElec_ULS2_weight(0)
 468 ,fPtElec_LS2_weight(0)
 469 ,fTOF01(0)
 470 ,fTOF02(0)
 471 ,fTOF03(0)
 472 ,fpid(0)
 473 ,fEoverP_pt(0)
 474 ,fEoverP_tpc(0)
 475 ,fTPC_pt(0)
 476 ,fTPC_p(0)
 477 ,fTPCnsigma_pt(0)
 478 ,fTPCnsigma_p(0)
 479 ,fTPCnsigma_pt_2D(0)
 480 ,fShowerShapeCut(0)
 481 ,fShowerShapeM02_EoverP(0)
 482 ,fShowerShapeM20_EoverP(0)
 483 ,fShowerShape_ha(0)
 484 ,fShowerShape_ele(0)
 485 ,fTPCnsigma_eta(0)
 486 ,fTPCnsigma_phi(0)
 487 ,fECluster(0)
 488 ,fECluster_pure(0)
 489 ,fECluster_not_exotic(0)
 490 ,fECluster_not_exotic1(0)
 491 ,fECluster_not_exotic2(0)
 492 ,fECluster_exotic(0)
 493 ,fNCluster_pure(0)
 494 ,fNCluster_pure_aod(0)
 495 ,fNCluster_ECluster(0)
 496 ,fNcells_energy(0)
 497 ,fNcells_energy_elec_selected(0)
 498 ,fNcells_energy_not_exotic(0)
 499 ,fEtaPhi(0)
 500 ,fEtaPhi_num(0)
 501 ,fEtaPhi_den(0)
 502 ,fEtaPhi_data(0)
 503 ,fpt_reco_pt_MC_num(0)
 504 ,fpt_reco_pt_MC_den(0)
 505 ,fVtxZ(0)
 506
 507 ,fVtxZ_new1(0)
 508 ,fVtxZ_new2(0)
 509 ,fVtxZ_new3(0)
 510 ,fVtxZ_new4(0)
 511
 512 ,fzRes1(0)
 513 ,fzRes2(0)
 514 ,fSPD_track_vtx1(0)
 515 ,fSPD_track_vtx2(0)
 516
 517
 518
 519 ,fEtad(0)
 520 ,fNTracks(0)
 521 ,fTrack_Multi(0)
 522 ,fNTracks_pt(0)
 523 ,fNTracks_eta(0)
 524 ,fNTracks_phi(0)
 525 ,fNClusters(0)
 526 ,fTPCNcls_EoverP(0)
 527 ,fTPCNcls_pid(0)
 528 ,fEta(0)
 529 ,fPhi(0)
 530 ,fR(0)
 531 ,fR_EoverP(0)
 532 ,fNcells(0)
 533 ,fNcells_EoverP(0)
 534 ,fNcells_electrons(0)
 535 ,fNcells_hadrons(0)
 536 ,fECluster_ptbins(0)
 537 ,fEoverP_ptbins(0)
 538 ,fEoverP_wSSCut(0)
 539 ,fM02_EoverP(0)
 540 ,fM20_EoverP(0)
 541 ,fTPCnsigma_eta_electrons(0)
 542 ,fTPCnsigma_eta_hadrons(0)
 543 ,fEoverP_pt_pions(0)
 544 ,ftpc_p_EoverPcut(0)
 545 ,fnsigma_p_EoverPcut(0)
 546 ,fEoverP_pt_pions2(0)
 547 ,fNcells_pt(0)
 548 ,fEoverP_pt_hadrons(0)
 549 ,fCEtaPhi_Inc(0)
 550 ,fCEtaPhi_ULS(0)
 551 ,fCEtaPhi_LS(0)
 552 ,fCEtaPhi_ULS_NoP(0)
 553 ,fCEtaPhi_LS_NoP(0)
 554 ,fCEtaPhi_ULS_Weight(0)
 555 ,fCEtaPhi_LS_Weight(0)
 556 ,fCEtaPhi_ULS_NoP_Weight(0)
 557 ,fCEtaPhi_LS_NoP_Weight(0)
 558 ,fInvMass(0)
 559 ,fInvMassBack(0)
 560 ,fDCA(0)
 561 ,fDCABack(0)
 562 ,fOpAngle(0)
 563 ,fOpAngleBack(0)
 564 ,fInvMass2(0)
 565 ,fInvMassBack2(0)
 566 ,fDCA2(0)
 567 ,fDCABack2(0)
 568 ,fOpAngle2(0)
 569 ,fOpAngleBack2(0)
 570 ,fMassCut(0.1)
 571 ,fEtaCutMin(-0.9)
 572 ,fEtaCutMax(0.9)
 573 ,fdPhiCut(0.05)
 574 ,fdEtaCut(0.05)
 575 ,fEoverPCutMin(0.8)
 576 ,fEoverPCutMax(1.2)
 577 ,fM20CutMin(0)
 578 ,fM20CutMax(10)
 579 ,fM02CutMin(0)
 580 ,fM02CutMax(10)
 581 ,fAngleCut(999)
 582 ,fChi2Cut(3.5)
 583 ,fDCAcut(999)
 584 ,fMassCutFlag(kTRUE)
 585 ,fAngleCutFlag(kFALSE)
 586 ,fChi2CutFlag(kFALSE)
 587 ,fDCAcutFlag(kFALSE)
 588 ,fAssHadronPtMin(0.5)
 589 ,fAssHadronPtMax(2.0)
 590 ,fPtBackgroundBeforeReco(0)
 591 ,fPtBackgroundBeforeReco2(0)
 592 ,fPtBackgroundBeforeReco_weight(0)
 593 ,fPtBackgroundBeforeReco2_weight(0)
 594 ,fpT_m_electron(0)
 595 ,fpT_gm_electron(0)
 596 ,fPtBackgroundAfterReco(0)
 597 ,fPtMinAsso(0.3)
 598 ,fTpcNclsAsso(80)
 599 ,fPtMCparticleAll(0)
 600 ,fPtMCparticleAll_nonPrimary(0)
 601 ,fPtMCparticleAlle_nonPrimary(0)
 602 ,fPtMCparticleAlle_Primary(0)
 603 ,fPtMCparticleReco(0)
 604 ,fPtMCparticleReco_nonPrimary(0)
 605 ,fPtMCparticleAllHfe1(0)
 606 ,fPtMCparticleRecoHfe1(0)
 607 ,fPtMCparticleAllHfe2(0)
 608 ,fPtMCparticleRecoHfe2(0)
 609 ,fPtMCelectronAfterAll(0)
 610 ,fPtMCelectronAfterAll_unfolding(0)
 611 ,fPtMCelectronAfterAll_nonPrimary(0)
 612 ,fPtMCelectronAfterAll_Primary(0)
 613 ,fPtMCpi0(0)
 614 ,fPtMCeta(0)
 615 ,fPtMCpi02(0)
 616 ,fPtMCeta2(0)
 617 ,fPtMCpi03(0)
 618 ,fPtMCeta3(0)
 619 ,fPtMC_EMCal_All(0)
 620 ,fPtMC_EMCal_Selected(0)
 621 ,fPtMC_TPC_All(0)
 622 ,fPtMC_TPC_Selected(0)
 623 ,fPt_track_match_den(0)
 624 ,fPt_track_match_num(0)
 625 ,fPtMCWithLabel(0)
 626 ,fPtMCWithoutLabel(0)
 627 ,fPtIsPhysicaPrimary(0)
 628 ,fCuts(0)
 629         //Lucile
 630         //,reader(0)
 631 ,fCFM(0)
 632 ,fPID(new AliHFEpid("hfePid"))
 633 ,fPIDqa(0)
 634 ,fMCstack(0)
 635 ,fRejectKinkMother(kFALSE)
 636 ,fMCtrack(0)
 637 ,fMCtrackMother(0)
 638 ,fMCtrackGMother(0)
 639 ,fMCtrackGGMother(0)
 640 ,fMCtrackGGGMother(0)
 641 ,fMCarray(0)
 642 ,fMCheader(0)
 643 ,fMCparticle(0)
 644 ,fMCparticle2(0)
 645 ,fMCparticleMother(0)
 646 ,fMCparticleGMother(0)
 647 ,fMCparticleGGMother(0)
 648 ,fMCparticleGGGMother(0)
 649 ,fEventHandler(0)
 650 ,fMCevent(0)
 651 ,fPoolMgr(0)
 652 ,fPool(0)
 653 ,fTracksClone(0)
 654 ,fTracks(0)
 655 ,fCEtaPhi_Inc_EM(0)
 656 ,fCEtaPhi_ULS_EM(0)
 657 ,fCEtaPhi_LS_EM(0)
 658 ,fCEtaPhi_ULS_Weight_EM(0)
 659 ,fCEtaPhi_LS_Weight_EM(0)
 660 ,fPoolNevents(0)
 661 ,fEventMixingFlag(0)
 662 ,fCEtaPhi_Inc_DiHadron(0)
 663 ,fPtTrigger_Inc(0)
 664         //,fEMCALRecoUtils(new AliEMCALRecoUtils)
 665         //,fEMCALGeo(0x0)
 666         //,fCaloUtils(0x0)
 667 ,fBitEGA(0)
 668         //,fEMCALRecoUtils(0)//exotic
 669 {
 670                 // Constructor
 671                 // Define input and output slots here
 672                 // Input slot #0 works with a TChain
 673
 674                 //exotic
 675                 // fEMCALRecoUtils  = new AliEMCALRecoUtils();
 676
 677         DefineInput(0, TChain::Class());
 678                 // Output slot #0 id reserved by the base class for AOD
 679                 // Output slot #1 writes into a TH1 container
 680                 // DefineOutput(1, TH1I::Class());
 681         DefineOutput(1, TList::Class());
 682                 //DefineOutput(3, TTree::Class());
 683 }
 684
 685         //______________________________________________________________________
 686 AliAnalysisTaskEMCalHFEpA::~AliAnalysisTaskEMCalHFEpA()
 687 {
 688         //Destructor
 689         delete fOutputList;
 690         delete fPID;
 691         delete fCFM;
 692         delete fPIDqa;
 693         //Lucile
 694         //delete reader;
 695         //if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
 696 }
 697
 698         //______________________________________________________________________
 699         //Create Output Objects
 700         //Here we can define the histograms and others output files
 701         //Called once
 702 void AliAnalysisTaskEMCalHFEpA::UserCreateOutputObjects()
 703 {
 704                 //______________________________________________________________________
 705                 //Initialize PID
 706         if(!fPID->GetNumberOfPIDdetectors())
 707     {
 708                 fPID->AddDetector("TPC", 0);
 709     }
 710
 711         fPID->SortDetectors();
 712
 713         fPIDqa = new AliHFEpidQAmanager();
 714         fPIDqa->Initialize(fPID);
 715                 //______________________________________________________________________
 716
 717                 //______________________________________________________________________
 718                 //Initialize correction Framework and Cuts
 719         fCFM = new AliCFManager;
 720         const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 721         fCFM->SetNStepParticle(kNcutSteps);
 722         for(Int_t istep = 0; istep < kNcutSteps; istep++) fCFM->SetParticleCutsList(istep, NULL);
 723
 724         if(!fCuts)
 725         {
 726                 AliWarning("Cuts not available. Default cuts will be used");
 727                 fCuts = new AliHFEcuts;
 728                 fCuts->CreateStandardCuts();
 729         }
 730
 731         fCuts->Initialize(fCFM);
 732                 //______________________________________________________________________
 733
 734                 ///___________________//Lucile
 735
 736                 //if(fIsAOD) {
 737                 //      reader = new AliCaloTrackAODReader();
 738                 //}
 739                         //___________________________________________________
 740                 ///Output Tlist
 741                 //Create TList
 742         fOutputList = new TList();
 743         fOutputList->SetOwner();
 744
 745                 //PIDqa
 746         fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 747
 748                 //Store the number of events
 749                 //Define the histo
 750         fNevent = new TH1F("fNevent","Number of Events",30,0,30);
 751         fNevent2 = new TH1F("fNevent2","Number of Events 2",30,0,30);
 752                 //And then, add to the output list
 753         fOutputList->Add(fNevent);
 754         fOutputList->Add(fNevent2);
 755
 756         fpid = new TH1F("fpid","PID flag",5,0,5);
 757         fOutputList->Add(fpid);
 758
 759         //pt Distribution
 760         fPtElec_Inc = new TH1F("fPtElec_Inc","Inclusive Electrons; p_{T} (GeV/c); Count",300,0,30);
 761
 762         fPtPrim = new TH1F("fPtPrim","Primary Electrons aod track; p_{T} (GeV/c); Count",300,0,30);
 763         fPtSec = new TH1F("fPtSec","Secundary Electrons aod track; p_{T} (GeV/c); Count",300,0,30);
 764         fPtPrim2 = new TH1F("fPtPrim2","Primary Electrons vtrack; p_{T} (GeV/c); Count",300,0,30);
 765         fPtSec2 = new TH1F("fPtSec2","Secundary Electrons vtrack; p_{T} (GeV/c); Count",300,0,30);
 766
 767
 768         fPtElec_ULS = new TH1F("fPtElec_ULS","ULS; p_{T} (GeV/c); Count",300,0,30);
 769         fPtElec_LS = new TH1F("fPtElec_LS","LS; p_{T} (GeV/c); Count",300,0,30);
 770
 771         fPtElec_ULS_NoPid = new TH1F("fPtElec_ULS_NoPid","ULS; p_{T} (GeV/c); Count",300,0,30);
 772         fPtElec_LS_NoPid = new TH1F("fPtElec_LS_NoPid","LS; p_{T} (GeV/c); Count",300,0,30);
 773
 774         fPtElec_ULS_MC = new TH1F("fPtElec_ULS_MC","ULS; p_{T} (GeV/c); Count",300,0,30);
 775         fPtElec_ULS_MC_weight = new TH1F("fPtElec_ULS_MC_weight","ULS; p_{T} (GeV/c); Count",300,0,30);
 776
 777         fPtElec_ULS_weight = new TH1F("fPtElec_ULS_weight","ULS; p_{T} (GeV/c); Count",300,0,30);
 778         fPtElec_LS_weight = new TH1F("fPtElec_LS_weight","LS; p_{T} (GeV/c); Count",300,0,30);
 779
 780         fTOF01 = new TH2F("fTOF01","",200,-20,20,200,-20,20);
 781         fTOF02 = new TH2F("fTOF02","",200,-20,20,200,-20,20);
 782         fTOF03 = new TH2F("fTOF03","",200,-20,20,200,-20,20);
 783
 784         if(fFillBackground){
 785                 fPtElec_ULS2 = new TH1F("fPtElec_ULS2","ULS; p_{T} (GeV/c); Count",300,0,30);
 786                 fPtElec_LS2 = new TH1F("fPtElec_LS2","LS; p_{T} (GeV/c); Count",300,0,30);
 787
 788                 fPtElec_ULS2_weight = new TH1F("fPtElec_ULS2_weight","ULS; p_{T} (GeV/c); Count",300,0,30);
 789                 fPtElec_LS2_weight = new TH1F("fPtElec_LS2_weight","LS; p_{T} (GeV/c); Count",300,0,30);
 790
 791         }
 792
 793         fPtTrigger_Inc = new TH1F("fPtTrigger_Inc","pT dist for Hadron Contamination; p_{t} (GeV/c); Count",300,0,30);
 794         fTPCnsigma_pt_2D = new TH2F("fTPCnsigma_pt_2D",";pt (GeV/c);TPC Electron N#sigma",1000,0.3,30,1000,-15,10);
 795         fShowerShapeCut = new TH2F("fShowerShapeCut","Shower Shape;M02;M20",500,0,1.8,500,0,1.8);
 796         fEtaPhi_num=new TH2F("fEtaPhi_num","#eta x #phi track;#phi;#eta",200,0.,5,50,-1.,1.);
 797         fEtaPhi_den=new TH2F("fEtaPhi_den","#eta x #phi track;#phi;#eta",200,0.,5,50,-1.,1.);
 798         fEtaPhi_data=new TH2F("fEtaPhi_data","#eta x #phi track;#phi;#eta",200,0.,5,50,-1.,1.);
 799
 800         fpt_reco_pt_MC_num=new TH2F("fpt_reco_pt_MC_num","pt reco x pt MC;pt reco; pt MC",300,0.,30,300,0.,30);
 801         fpt_reco_pt_MC_den=new TH2F("fpt_reco_pt_MC_den","pt reco x pt MC;pt reco; pt MC",300,0.,30,300,0.,30);
 802
 803
 804         fCharge_n = new TH1F("fCharge_n","Inclusive Electrons (Negative Charge); p_{t} (GeV/c); Count",200,0,30);
 805         fCharge_p = new TH1F("fCharge_p","Inclusive Positrons (Positive Charge); p_{t} (GeV/c); Count",200,0,30);
 806
 807         fECluster_pure= new TH1F("fECluster_pure", ";ECluster pure",2000,0,100);
 808         fECluster_not_exotic= new TH1F("fECluster_not_exotic", ";ECluster not exotic - function ",2000,0,100);
 809
 810         fECluster_not_exotic1= new TH1F("fECluster_not_exotic1", ";ECluster not exotic Ncells > E/3+1",2000,0,100);
 811
 812         fECluster_not_exotic2= new TH1F("fECluster_not_exotic2", ";ECluster not exotic 2",2000,0,100);
 813         fECluster_exotic= new TH1F("fECluster_exotic", ";ECluster exotic",2000,0,100);
 814
 815         //not associated with tracks
 816         fNCluster_pure= new TH1F("fNCluster_pure", ";Number of clusters - pure",2000,-1,1999);
 817         fNCluster_pure_aod= new TH1F("fNCluster_pure_aod", ";Number of clusters - pure -aod",2000,-1,1999);
 818         fNCluster_ECluster= new TH2F("fNCluster_ECluster", ";Number of clusters vs. Energy of Cluster",2000,-1,1999, 4000, -1, 1999);
 819
 820         fNcells_energy= new TH2F("fNcells_energy", "all clusters;Number of cells;Energy of Cluster",100,0,100, 2000, 0, 100);
 821         fNcells_energy_elec_selected= new TH2F("fNcells_energy_elec_selected", "clusters for electrons on TPC;Number of cells;Energy of Cluster",100,0,100, 2000, 0, 100);
 822         fNcells_energy_not_exotic= new TH2F("fNcells_energy_not_exotic", "not exotic cluster;Number of cells;Energy of Cluster ",100,0,100, 2000, 0, 100);
 823
 824         if(fUseEMCal){
 825
 826                 if(!fIsAOD){
 827                         fTime = new TH2D("fTime","Cells Cluster Time; p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 828                         fTime2 = new TH2D("fTime2","Cells Cluster Time;  p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 829                 }
 830
 831                 ftimingEle = new TH2D("ftimingEle","Cluster Time;  p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 832                 ftimingEle2 = new TH2D("ftimingEle2","Cluster Time;  p_{T} (GeV/c); Time (s)",300,0,30,1000,1e-8,1e-5);
 833
 834                 fShowerShape_ha = new TH2F("fShowerShape_ha","Shower Shape hadrons;M02;M20",500,0,1.8,500,0,1.8);
 835                 fShowerShape_ele = new TH2F("fShowerShape_ele","Shower Shape electrons;M02;M20",500,0,1.8,500,0,1.8);
 836
 837                 fShowerShapeM02_EoverP = new TH2F("fShowerShapeM02_EoverP","Shower Shape;M02;E/p",500,0,1.8,500,0,1.8);
 838                 fShowerShapeM20_EoverP = new TH2F("fShowerShapeM20_EoverP","Shower Shape;M20;E/p",500,0,1.8,500,0,1.8);
 839
 840         }
 841
 842         fOutputList->Add(fTOF01);
 843         fOutputList->Add(fTOF02);
 844         fOutputList->Add(fTOF03);
 845
 846         fOutputList->Add(fEtaPhi_num);
 847         fOutputList->Add(fEtaPhi_den);
 848         fOutputList->Add(fEtaPhi_data);
 849
 850         fOutputList->Add(fpt_reco_pt_MC_num);
 851         fOutputList->Add(fpt_reco_pt_MC_den);
 852
 853
 854         fOutputList->Add(fPtElec_Inc);
 855         fOutputList->Add(fPtElec_ULS);
 856         fOutputList->Add(fPtElec_LS);
 857         fOutputList->Add(fPtElec_ULS_NoPid);
 858         fOutputList->Add(fPtElec_LS_NoPid);
 859         fOutputList->Add(fPtElec_ULS_MC);
 860         fOutputList->Add(fPtElec_ULS_MC_weight);
 861
 862         fOutputList->Add(fPtPrim);
 863         fOutputList->Add(fPtSec);
 864         fOutputList->Add(fPtPrim2);
 865         fOutputList->Add(fPtSec2);
 866
 867
 868
 869         fOutputList->Add(fPtElec_ULS_weight);
 870         fOutputList->Add(fPtElec_LS_weight);
 871
 872         if(fFillBackground){
 873                 fOutputList->Add(fPtElec_ULS2);
 874                 fOutputList->Add(fPtElec_LS2);
 875                 fOutputList->Add(fPtElec_ULS2_weight);
 876                 fOutputList->Add(fPtElec_LS2_weight);
 877         }
 878
 879
 880         fOutputList->Add(fPtTrigger_Inc);
 881         fOutputList->Add(fTPCnsigma_pt_2D);
 882         fOutputList->Add(fShowerShapeCut);
 883
 884         fOutputList->Add(fCharge_n);
 885         fOutputList->Add(fCharge_p);
 886
 887         fOutputList->Add(fECluster_pure);
 888         fOutputList->Add(fECluster_not_exotic);
 889         fOutputList->Add(fECluster_not_exotic1);
 890         fOutputList->Add(fECluster_not_exotic2);
 891         fOutputList->Add(fECluster_exotic);
 892
 893         fOutputList->Add(fNCluster_pure);
 894         fOutputList->Add(fNCluster_pure_aod);
 895
 896         fOutputList->Add(fNCluster_ECluster);
 897         fOutputList->Add(fNcells_energy);
 898         fOutputList->Add(fNcells_energy_elec_selected);
 899         fOutputList->Add(fNcells_energy_not_exotic);
 900
 901
 902         if(fUseEMCal){
 903
 904                 if(!fIsAOD){
 905                         fOutputList->Add(fTime);
 906                         fOutputList->Add(fTime2);
 907
 908                 }
 909
 910                 fOutputList->Add(ftimingEle);
 911                 fOutputList->Add(ftimingEle2);
 912
 913                 fOutputList->Add(fShowerShape_ha);
 914                 fOutputList->Add(fShowerShape_ele);
 915
 916                 fOutputList->Add(fShowerShapeM02_EoverP);
 917                 fOutputList->Add(fShowerShapeM20_EoverP);
 918
 919
 920
 921         }
 922
 923         fVtxZ_new1= new  TH1F("fVtxZ_new1","fVtxZ_new1",4000, -50,50);
 924         fVtxZ_new2= new  TH1F("fVtxZ_new2","fVtxZ_new2",4000, -50,50);
 925         fVtxZ_new3= new  TH1F("fVtxZ_new3","fVtxZ_new3",4000, -50,50);
 926         fVtxZ_new4= new  TH1F("fVtxZ_new4","fVtxZ_new4",4000, -50,50);
 927
 928         fzRes1= new  TH1F("fzRes1","fzRes1",4000, 0,1);
 929         fzRes2= new  TH1F("fzRes2","fzRes2",4000, 0,1);
 930         fSPD_track_vtx1= new  TH1F("fSPD_track_vtx1","fSPD_track_vtx1",4000, -5,5);
 931         fSPD_track_vtx2= new  TH1F("fSPD_track_vtx2","fSPD_track_vtx2",4000, -5,5);
 932
 933
 934                 //General Histograms
 935
 936                 //Steps
 937                 //Step 1: Before Track cuts
 938                 //Step 2: Before PID
 939                 //Step 3: After PID
 940
 941         fEoverP_pt = new TH2F *[3];
 942         fTPC_p = new TH2F *[3];
 943         fTPCnsigma_p = new TH2F *[3];
 944
 945         fECluster= new TH1F *[3];
 946         fEtaPhi= new TH2F *[3];
 947         fVtxZ= new  TH1F *[3];
 948         fEtad= new  TH1F *[3];
 949         fNTracks= new  TH1F *[3];
 950
 951         fNTracks_pt= new  TH2F *[3];
 952         fNTracks_eta= new  TH2F *[3];
 953         fNTracks_phi= new  TH2F *[3];
 954
 955         fNClusters= new TH1F *[3];
 956         fTPCNcls_EoverP= new TH2F *[3];
 957         fTPCNcls_pid=new TH2F *[4];
 958
 959
 960
 961         for(Int_t i = 0; i < 3; i++)
 962         {
 963                 fEoverP_pt[i] = new TH2F(Form("fEoverP_pt%d",i),";p_{t} (GeV/c);E / p ",1000,0,30,500,0,2);
 964                 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);
 965                 fTPCnsigma_p[i] = new TH2F(Form("fTPCnsigma_p%d",i),";p (GeV/c);TPC Electron N#sigma",1000,0.3,15,1000,-15,10);
 966
 967
 968                 fECluster[i]= new TH1F(Form("fECluster%d",i), ";ECluster",2000, 0,100);
 969                 fEtaPhi[i]= new TH2F(Form("fEtaPhi%d",i),"#eta x #phi Clusters;#phi;#eta",200,0.,5,50,-1.,1.);
 970                 fVtxZ[i]= new  TH1F(Form("fVtxZ%d",i),"VtxZ",1000, -50,50);
 971                 fEtad[i]= new  TH1F(Form("fEtad%d",i),"Eta distribution",200, -1.2,1.2);
 972                 fNTracks[i]= new  TH1F(Form("fNTracks%d",i),"NTracks",1000, 0,5000);
 973
 974                 fNTracks_pt[i]= new  TH2F(Form("fNTracks_pt%d",i),"NTracks vs. pt",1000, 0,5000, 1000, 0, 100);
 975                 fNTracks_eta[i]= new  TH2F(Form("fNTracks_eta%d",i),"NTracks vs. pt",1000, 0,5000, 500, -1.0, 1.0);
 976                 fNTracks_phi[i]= new  TH2F(Form("fNTracks_phi%d",i),"NTracks vs. pt",1000, 0,5000, 500, 0, 5.0);
 977
 978
 979
 980                 fNClusters[i]= new TH1F(Form("fNClusters%d",i),"fNClusters0",200, 0,100);
 981                 fTPCNcls_EoverP[i]= new TH2F(Form("fTPCNcls_EoverP%d",i),"TPCNcls_EoverP",1000,0,200,200,0,2);
 982
 983
 984
 985                 fOutputList->Add(fEoverP_pt[i]);
 986                 fOutputList->Add(fTPC_p[i]);
 987                 fOutputList->Add(fTPCnsigma_p[i]);
 988
 989
 990                 fOutputList->Add(fECluster[i]);
 991                 fOutputList->Add(fEtaPhi[i]);
 992                 fOutputList->Add(fVtxZ[i]);
 993                 fOutputList->Add(fEtad[i]);
 994                 fOutputList->Add(fNTracks[i]);
 995
 996                 fOutputList->Add(fNTracks_pt[i]);
 997                 fOutputList->Add(fNTracks_eta[i]);
 998                 fOutputList->Add(fNTracks_phi[i]);
 999
1000                 fOutputList->Add(fNClusters[i]);
1001                 fOutputList->Add(fTPCNcls_EoverP[i]);
1002         }
1003
1004         fTrack_Multi= new  TH1F("fTrack_Multi","fTrack_Multi",1000, 0,1000);
1005
1006         for(Int_t i = 0; i < 4; i++)
1007         {
1008                 fTPCNcls_pid[i]= new TH2F(Form("fTPCNcls_pid%d",i),"fTPCNcls_pid;NCls;NCls for PID",200,0,200,200,0,200);
1009                 fOutputList->Add(fTPCNcls_pid[i]);
1010         }
1011
1012                 //pt bin
1013         Int_t fPtBin[7] = {1,2,4,6,8,10,15};
1014
1015         fEoverP_tpc = new TH2F *[6];
1016         fTPC_pt = new TH1F *[6];
1017         fTPCnsigma_pt = new TH1F *[6];
1018
1019         fEta=new TH1F *[6];
1020         fPhi=new TH1F *[6];
1021         fR=new TH1F *[6];
1022         fR_EoverP=new TH2F *[6];
1023         fNcells=new TH1F *[6];
1024         fNcells_EoverP=new TH2F *[6];
1025         fM02_EoverP= new TH2F *[6];
1026         fM20_EoverP= new TH2F *[6];
1027         fEoverP_ptbins=new TH1F *[6];
1028         fECluster_ptbins=new TH1F *[6];
1029         fEoverP_wSSCut=new TH1F *[6];
1030         fNcells_electrons=new TH1F *[6];
1031         fNcells_hadrons=new TH1F *[6];
1032         fTPCnsigma_eta_electrons=new TH2F *[6];
1033         fTPCnsigma_eta_hadrons=new TH2F *[6];
1034
1035         if(fCorrelationFlag)
1036         {
1037                 fCEtaPhi_Inc = new TH2F *[6];
1038                 fCEtaPhi_Inc_DiHadron = new TH2F *[6];
1039
1040                 fCEtaPhi_ULS = new TH2F *[6];
1041                 fCEtaPhi_LS = new TH2F *[6];
1042                 fCEtaPhi_ULS_NoP = new TH2F *[6];
1043                 fCEtaPhi_LS_NoP = new TH2F *[6];
1044
1045                 fCEtaPhi_ULS_Weight = new TH2F *[6];
1046                 fCEtaPhi_LS_Weight = new TH2F *[6];
1047                 fCEtaPhi_ULS_NoP_Weight = new TH2F *[6];
1048                 fCEtaPhi_LS_NoP_Weight = new TH2F *[6];
1049
1050                 fCEtaPhi_Inc_EM = new TH2F *[6];
1051
1052                 fCEtaPhi_ULS_EM = new TH2F *[6];
1053                 fCEtaPhi_LS_EM = new TH2F *[6];
1054
1055                 fCEtaPhi_ULS_Weight_EM = new TH2F *[6];
1056                 fCEtaPhi_LS_Weight_EM = new TH2F *[6];
1057
1058                 fInvMass = new TH1F("fInvMass","",200,0,0.3);
1059                 fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
1060                 fDCA = new TH1F("fDCA","",200,0,1);
1061                 fDCABack = new TH1F("fDCABack","",200,0,1);
1062                 fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
1063                 fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
1064
1065                 fOutputList->Add(fInvMass);
1066                 fOutputList->Add(fInvMassBack);
1067                 fOutputList->Add(fDCA);
1068                 fOutputList->Add(fDCABack);
1069                 fOutputList->Add(fOpAngle);
1070                 fOutputList->Add(fOpAngleBack);
1071         }
1072
1073         if(fFillBackground){
1074
1075                 fInvMass = new TH1F("fInvMass","",200,0,0.3);
1076                 fInvMassBack = new TH1F("fInvMassBack","",200,0,0.3);
1077                 fDCA = new TH1F("fDCA","",200,0,1);
1078                 fDCABack = new TH1F("fDCABack","",200,0,1);
1079                 fOpAngle = new TH1F("fOpAngle","",200,0,0.5);
1080                 fOpAngleBack = new TH1F("fOpAngleBack","",200,0,0.5);
1081
1082                 fOutputList->Add(fInvMass);
1083                 fOutputList->Add(fInvMassBack);
1084                 fOutputList->Add(fDCA);
1085                 fOutputList->Add(fDCABack);
1086                 fOutputList->Add(fOpAngle);
1087                 fOutputList->Add(fOpAngleBack);
1088
1089                         //histos for TPC-only
1090                 fInvMass2 = new TH1F("fInvMass2","",200,0,0.3);
1091                 fInvMassBack2 = new TH1F("fInvMassBack2","",200,0,0.3);
1092                 fDCA2 = new TH1F("fDCA2","",200,0,1);
1093                 fDCABack2 = new TH1F("fDCABack2","",200,0,1);
1094                 fOpAngle2 = new TH1F("fOpAngle2","",200,0,0.5);
1095                 fOpAngleBack2 = new TH1F("fOpAngleBack2","",200,0,0.5);
1096
1097                 fOutputList->Add(fInvMass2);
1098                 fOutputList->Add(fInvMassBack2);
1099                 fOutputList->Add(fDCA2);
1100                 fOutputList->Add(fDCABack2);
1101                 fOutputList->Add(fOpAngle2);
1102                 fOutputList->Add(fOpAngleBack2);
1103
1104         }
1105
1106         for(Int_t i = 0; i < 6; i++)
1107         {
1108                 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);
1109                 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);
1110                 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);
1111
1112                 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);
1113                 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);
1114                 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);
1115                 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);
1116                 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);
1117                 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);
1118                 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);
1119                 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);
1120                 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);
1121                 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);
1122                 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);
1123                 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);
1124                 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);
1125                 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);
1126                 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);
1127
1128                 fOutputList->Add(fEoverP_tpc[i]);
1129                 fOutputList->Add(fTPC_pt[i]);
1130                 fOutputList->Add(fTPCnsigma_pt[i]);
1131
1132                 fOutputList->Add(fEta[i]);
1133                 fOutputList->Add(fPhi[i]);
1134                 fOutputList->Add(fR[i]);
1135                 fOutputList->Add(fR_EoverP[i]);
1136                 fOutputList->Add(fNcells[i]);
1137                 fOutputList->Add(fNcells_electrons[i]);
1138                 fOutputList->Add(fNcells_hadrons[i]);
1139                 fOutputList->Add(fNcells_EoverP[i]);
1140                 fOutputList->Add(fECluster_ptbins[i]);
1141                 fOutputList->Add(fEoverP_ptbins[i]);
1142                 fOutputList->Add(fEoverP_wSSCut[i]);
1143                 fOutputList->Add(fM02_EoverP[i]);
1144                 fOutputList->Add(fM20_EoverP[i]);
1145                 fOutputList->Add(fTPCnsigma_eta_electrons[i]);
1146                 fOutputList->Add(fTPCnsigma_eta_hadrons[i]);
1147
1148
1149                 if(fCorrelationFlag)
1150                 {
1151                         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);
1152                         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);
1153
1154                         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);
1155                         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);
1156                         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);
1157                         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);
1158
1159                         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);
1160                         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);
1161                         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);
1162                         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);
1163
1164                         fOutputList->Add(fCEtaPhi_Inc[i]);
1165                         fOutputList->Add(fCEtaPhi_Inc_DiHadron[i]);
1166
1167                         fOutputList->Add(fCEtaPhi_ULS[i]);
1168                         fOutputList->Add(fCEtaPhi_LS[i]);
1169                         fOutputList->Add(fCEtaPhi_ULS_NoP[i]);
1170                         fOutputList->Add(fCEtaPhi_LS_NoP[i]);
1171
1172                         fOutputList->Add(fCEtaPhi_ULS_Weight[i]);
1173                         fOutputList->Add(fCEtaPhi_LS_Weight[i]);
1174                         fOutputList->Add(fCEtaPhi_ULS_NoP_Weight[i]);
1175                         fOutputList->Add(fCEtaPhi_LS_NoP_Weight[i]);
1176
1177                         if(fEventMixingFlag)
1178                         {
1179                                 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);
1180
1181                                 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);
1182                                 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);
1183
1184                                 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);
1185                                 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);
1186
1187                                 fOutputList->Add(fCEtaPhi_Inc_EM[i]);
1188
1189                                 fOutputList->Add(fCEtaPhi_ULS_EM[i]);
1190                                 fOutputList->Add(fCEtaPhi_LS_EM[i]);
1191
1192                                 fOutputList->Add(fCEtaPhi_ULS_Weight_EM[i]);
1193                                 fOutputList->Add(fCEtaPhi_LS_Weight_EM[i]);
1194                         }
1195                 }
1196         }
1197
1198                 //pt integrated
1199         fTPCnsigma_eta = new TH2F("fTPCnsigma_eta",";Pseudorapidity #eta; TPC signal - <TPC signal>_{elec} (#sigma)",200,-0.9,0.9,200,-15,15);
1200         fTPCnsigma_phi = new TH2F("fTPCnsigma_phi",";Azimuthal Angle #phi; TPC signal - <TPC signal>_{elec} (#sigma)",200,0,2*TMath::Pi(),200,-15,15);
1201
1202
1203
1204         fNcells_pt=new TH2F("fNcells_pt","fNcells_pt",1000, 0,20,100,0,30);
1205         fEoverP_pt_pions= new TH2F("fEoverP_pt_pions","fEoverP_pt_pions",1000,0,30,500,0,2);
1206
1207         ftpc_p_EoverPcut= new TH2F("ftpc_p_EoverPcut","ftpc_p_EoverPcut",1000,0,30,200,20,200);
1208         fnsigma_p_EoverPcut= new TH2F("fnsigma_p_EoverPcut","fnsigma_p_EoverPcut",1000,0,30,500,-15,15);
1209
1210         fEoverP_pt_pions2= new TH2F("fEoverP_pt_pions2","fEoverP_pt_pions2",1000,0,30,500,0,2);
1211         fEoverP_pt_hadrons= new TH2F("fEoverP_pt_hadrons","fEoverP_pt_hadrons",1000,0,30,500,0,2);
1212
1213
1214         fOutputList->Add(fTPCnsigma_eta);
1215         fOutputList->Add(fTPCnsigma_phi);
1216
1217         fOutputList->Add(fNcells_pt);
1218         fOutputList->Add(fEoverP_pt_pions);
1219
1220         fOutputList->Add(ftpc_p_EoverPcut);
1221         fOutputList->Add(fnsigma_p_EoverPcut);
1222
1223         fOutputList->Add(fEoverP_pt_pions2);
1224         fOutputList->Add(fEoverP_pt_hadrons);
1225
1226         fOutputList->Add(fVtxZ_new1);
1227         fOutputList->Add(fVtxZ_new2);
1228         fOutputList->Add(fVtxZ_new3);
1229         fOutputList->Add(fVtxZ_new4);
1230
1231         fOutputList->Add(fzRes1);
1232         fOutputList->Add(fzRes2);
1233         fOutputList->Add(fSPD_track_vtx1);
1234         fOutputList->Add(fSPD_track_vtx2);
1235
1236
1237
1238                 //__________________________________________________________________
1239                 //Efficiency studies
1240         if(fIsMC)
1241         {
1242                 fPtBackgroundBeforeReco = new TH1F("fPtBackgroundBeforeReco",";p_{T} (GeV/c);Count",300,0,30);
1243                 fPtBackgroundBeforeReco_weight = new TH1F("fPtBackgroundBeforeReco_weight",";p_{T} (GeV/c);Count",300,0,30);
1244                 if(fFillBackground)fPtBackgroundBeforeReco2 = new TH1F("fPtBackgroundBeforeReco2",";p_{T} (GeV/c);Count",300,0,30);
1245                 if(fFillBackground)fPtBackgroundBeforeReco2_weight = new TH1F("fPtBackgroundBeforeReco2_weight",";p_{T} (GeV/c);Count",300,0,30);
1246                 fpT_m_electron= new TH2F("fpT_m_electron","fpT_m_electron",300,0,30,300,0,30);
1247                 fpT_gm_electron= new TH2F("fpT_gm_electron","fpT_gm_electron",300,0,30,300,0,30);
1248
1249                 fPtBackgroundAfterReco = new TH1F("fPtBackgroundAfterReco",";p_{T} (GeV/c);Count",300,0,30);
1250                 fPtMCparticleAll = new TH1F("fPtMCparticleAll",";p_{T} (GeV/c);Count",200,0,40);
1251                 fPtMCparticleReco = new TH1F("fPtMCparticleReco",";p_{T} (GeV/c);Count",200,0,40);
1252
1253                 fPtMCparticleAll_nonPrimary = new TH1F("fPtMCparticleAll_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
1254                 fPtMCparticleAlle_nonPrimary = new TH1F("fPtMCparticleAlle_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
1255                 fPtMCparticleAlle_Primary = new TH1F("fPtMCparticleAlle_Primary",";p_{T} (GeV/c);Count",200,0,40);
1256
1257                 fPtMCparticleReco_nonPrimary = new TH1F("fPtMCparticleReco_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
1258
1259                 fPtMCparticleAllHfe1 = new TH1F("fPtMCparticleAllHfe1",";p_{t} (GeV/c);Count",200,0,40);
1260                 fPtMCparticleRecoHfe1 = new TH1F("fPtMCparticleRecoHfe1",";p_{t} (GeV/c);Count",200,0,40);
1261                 fPtMCparticleAllHfe2 = new TH1F("fPtMCparticleAllHfe2",";p_{t} (GeV/c);Count",200,0,40);
1262                 fPtMCparticleRecoHfe2 = new TH1F("fPtMCparticleRecoHfe2",";p_{t} (GeV/c);Count",200,0,40);
1263
1264                 fPtMCelectronAfterAll = new TH1F("fPtMCelectronAfterAll",";p_{T} (GeV/c);Count",200,0,40);
1265                 fPtMCelectronAfterAll_unfolding = new TH1F("fPtMCelectronAfterAll_unfolding",";p_{T} (GeV/c);Count",200,0,40);
1266                 fPtMCelectronAfterAll_nonPrimary = new TH1F("fPtMCelectronAfterAll_nonPrimary",";p_{T} (GeV/c);Count",200,0,40);
1267                 fPtMCelectronAfterAll_Primary = new TH1F("fPtMCelectronAfterAll_Primary",";p_{T} (GeV/c);Count",200,0,40);
1268
1269
1270
1271                 fPtMCpi0 = new TH1F("fPtMCpi0",";p_{t} (GeV/c);Count",200,0,30);
1272                 fPtMCeta = new TH1F("fPtMCeta",";p_{T} (GeV/c);Count",200,0,30);
1273                 fPtMCpi02 = new TH1F("fPtMCpi02",";p_{t} (GeV/c);Count",200,0,30);
1274                 fPtMCeta2 = new TH1F("fPtMCeta2",";p_{T} (GeV/c);Count",200,0,30);
1275                 fPtMCpi03 = new TH1F("fPtMCpi03",";p_{t} (GeV/c);Count",200,0,30);
1276                 fPtMCeta3 = new TH1F("fPtMCeta3",";p_{T} (GeV/c);Count",200,0,30);
1277
1278                 fPtMC_EMCal_All= new TH1F("fPtMC_EMCal_All",";p_{t} (GeV/c);Count",200,0,40);
1279                 fPtMC_EMCal_Selected= new TH1F("fPtMC_EMCal_Selected",";p_{t} (GeV/c);Count",200,0,40);
1280                 fPtMC_TPC_All= new TH1F("fPtMC_TPC_All",";p_{T} (GeV/c);Count",200,0,40);
1281                 fPtMC_TPC_Selected = new TH1F("fPtMC_TPC_Selected",";p_{T} (GeV/c);Count",200,0,40);
1282
1283                 fPt_track_match_den = new TH1F("fPt_track_match_den",";p_{T} (GeV/c);Count",200,0,40);
1284                 fPt_track_match_num = new TH1F("fPt_track_match_num",";p_{T} (GeV/c);Count",200,0,40);
1285
1286                 fPtMCWithLabel = new TH1F("fPtMCWithLabel",";p_{t} (GeV/c);Count",200,0,40);
1287                 fPtMCWithoutLabel = new TH1F("fPtMCWithoutLabel",";p_{t} (GeV/c);Count",200,0,40);
1288                 fPtIsPhysicaPrimary = new TH1F("fPtIsPhysicaPrimary",";p_{t} (GeV/c);Count",200,0,40);
1289
1290                 fOutputList->Add(fPtBackgroundBeforeReco);
1291                 fOutputList->Add(fPtBackgroundBeforeReco_weight);
1292
1293                 if(fFillBackground) fOutputList->Add(fPtBackgroundBeforeReco2);
1294                 if(fFillBackground) fOutputList->Add(fPtBackgroundBeforeReco2_weight);
1295
1296                 fOutputList->Add(fpT_m_electron);
1297                 fOutputList->Add(fpT_gm_electron);
1298
1299                 fOutputList->Add(fPtBackgroundAfterReco);
1300                 fOutputList->Add(fPtMCparticleAll);
1301                 fOutputList->Add(fPtMCparticleReco);
1302
1303                 fOutputList->Add(fPtMCparticleAll_nonPrimary);
1304                 fOutputList->Add(fPtMCparticleAlle_nonPrimary);
1305
1306                 fOutputList->Add(fPtMCparticleAlle_Primary);
1307                 fOutputList->Add(fPtMCparticleReco_nonPrimary);
1308
1309                 fOutputList->Add(fPtMCparticleAllHfe1);
1310                 fOutputList->Add(fPtMCparticleRecoHfe1);
1311                 fOutputList->Add(fPtMCparticleAllHfe2);
1312                 fOutputList->Add(fPtMCparticleRecoHfe2);
1313                 fOutputList->Add(fPtMCelectronAfterAll);
1314                 fOutputList->Add(fPtMCelectronAfterAll_unfolding);
1315
1316                 fOutputList->Add(fPtMCelectronAfterAll_nonPrimary);
1317                 fOutputList->Add(fPtMCelectronAfterAll_Primary);
1318
1319
1320
1321                 fOutputList->Add(fPtMCpi0);
1322                 fOutputList->Add(fPtMCeta);
1323                 fOutputList->Add(fPtMCpi02);
1324                 fOutputList->Add(fPtMCeta2);
1325                 fOutputList->Add(fPtMCpi03);
1326                 fOutputList->Add(fPtMCeta3);
1327                 fOutputList->Add(fPtMC_EMCal_All);
1328                 fOutputList->Add(fPtMC_EMCal_Selected);
1329                 fOutputList->Add(fPtMC_TPC_All);
1330                 fOutputList->Add(fPtMC_TPC_Selected);
1331
1332                 fOutputList->Add(fPt_track_match_den);
1333                 fOutputList->Add(fPt_track_match_num);
1334
1335                 fOutputList->Add(fPtMCWithLabel);
1336                 fOutputList->Add(fPtMCWithoutLabel);
1337                 fOutputList->Add(fPtIsPhysicaPrimary);
1338         }
1339
1340         fCentralityHist = new TH1F("fCentralityHist",";Centrality (%); Count",1000000,0,100);
1341         fCentralityHistPass = new TH1F("fCentralityHistPass",";Centrality (%); Count",1000000,0,100);
1342         fOutputList->Add(fCentralityHist);
1343         fOutputList->Add(fCentralityHistPass);
1344
1345         //______________________________________________________________________
1346         //Mixed event analysis
1347         if(fEventMixingFlag)
1348         {
1349                 fPoolNevents = new TH1F("fPoolNevents","Event Mixing Statistics; Number of events; Count",1000,0,1000);
1350                 fOutputList->Add(fPoolNevents);
1351
1352                 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.
1353                 Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
1354
1355                 Int_t nCentralityBins  = 15;
1356                 Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
1357
1358                 Int_t nZvtxBins  = 9;
1359                 Double_t vertexBins[] = {-10, -7, -5, -3, -1, 1, 3, 5, 7, 10};
1360
1361                 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, (Double_t*) centralityBins, nZvtxBins, (Double_t*) vertexBins);
1362         }
1363                 //______________________________________________________________________
1364
1365         PostData(1, fOutputList);
1366
1367                 ///______________________________________________________________________
1368 }
1369
1370 //______________________________________________________________________
1371 //Main loop
1372 //Called for each event
1373 void AliAnalysisTaskEMCalHFEpA::UserExec(Option_t *)
1374 {
1375                 //Check Event
1376         fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1377         fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1378
1379         if(!(fESD || fAOD))
1380         {
1381                 printf("ERROR: fESD & fAOD not available\n");
1382                 return;
1383         }
1384
1385         fVevent = dynamic_cast<AliVEvent*>(InputEvent());
1386
1387         if(!fVevent)
1388         {
1389                 printf("ERROR: fVEvent not available\n");
1390                 return;
1391         }
1392
1393                 //Check Cuts
1394         if(!fCuts)
1395         {
1396                 AliError("HFE cuts not available");
1397                 return;
1398         }
1399                 //Check PID
1400         if(!fPID->IsInitialized())
1401         {
1402                         // Initialize PID with the given run number
1403                 AliWarning("PID not initialised, get from Run no");
1404
1405                 if(fIsAOD)
1406                 {
1407                         fPID->InitializePID(fAOD->GetRunNumber());
1408                 }
1409                 else
1410                 {
1411                         fPID->InitializePID(fESD->GetRunNumber());
1412                 }
1413         }
1414
1415                 //PID response
1416         fPidResponse = fInputHandler->GetPIDResponse();
1417
1418
1419                 //Check PID response
1420         if(!fPidResponse)
1421         {
1422                 AliDebug(1, "Using default PID Response");
1423                 fPidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
1424         }
1425
1426         fPID->SetPIDResponse(fPidResponse);
1427
1428         fCFM->SetRecEventInfo(fVevent);
1429
1430         Double_t *fListOfmotherkink = 0;
1431         Int_t fNumberOfVertices = 0;
1432         Int_t fNumberOfMotherkink = 0;
1433
1434
1435         //total event before event selection
1436         fNevent->Fill(1);
1437
1438         //______________________________________________________________________
1439         //Vertex Selection
1440         if(fIsAOD)
1441         {
1442                 const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
1443                 if(!trkVtx || trkVtx->GetNContributors()<=0) return;
1444                 TString vtxTtl = trkVtx->GetTitle();
1445                 if(!vtxTtl.Contains("VertexerTracks")) return;
1446                         //Float_t zvtx = trkVtx->GetZ();
1447                 Float_t zvtx = -100;
1448                 zvtx=trkVtx->GetZ();
1449                 fZvtx = zvtx;
1450
1451                 fVtxZ_new1->Fill(fZvtx);
1452
1453                 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
1454                 if(spdVtx->GetNContributors()<=0) return;
1455                 TString vtxTyp = spdVtx->GetTitle();
1456                 Double_t cov[6]={0};
1457                 spdVtx->GetCovarianceMatrix(cov);
1458                 Double_t zRes = TMath::Sqrt(cov[5]);
1459
1460                 fzRes1->Fill(zRes);
1461                 if(vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1462                 fzRes2->Fill(zRes);
1463
1464                 fSPD_track_vtx1->Fill(spdVtx->GetZ() - trkVtx->GetZ());
1465                 if(TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1466                 fSPD_track_vtx2->Fill(spdVtx->GetZ() - trkVtx->GetZ());
1467
1468
1469                 if(TMath::Abs(zvtx) > 10) return;
1470                 fVtxZ_new2->Fill(fZvtx);
1471
1472                 if(fabs(zvtx>10.0))return;
1473                 fVtxZ_new3->Fill(fZvtx);
1474
1475
1476                 //Look for kink mother for AOD
1477
1478                 fNumberOfVertices = 0;
1479                 fNumberOfMotherkink = 0;
1480
1481                 if(fIsAOD)
1482                 {
1483                         fNumberOfVertices = fAOD->GetNumberOfVertices();
1484
1485                         fListOfmotherkink = new Double_t[fNumberOfVertices];
1486
1487                         for(Int_t ivertex=0; ivertex < fNumberOfVertices; ivertex++)
1488                         {
1489                                 AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
1490                                 if(!aodvertex) continue;
1491                                 if(aodvertex->GetType()==AliAODVertex::kKink)
1492                                 {
1493                                         AliAODTrack *mother1 = (AliAODTrack *) aodvertex->GetParent();
1494                                         if(!mother1) continue;
1495                                         Int_t idmother = mother1->GetID();
1496                                         fListOfmotherkink[fNumberOfMotherkink] = idmother;
1497                                         fNumberOfMotherkink++;
1498                                 }
1499                         }
1500                 }
1501         }
1502         else
1503         {
1504
1505
1506
1507                 /// ESD
1508                 const AliESDVertex* trkVtx = fESD->GetPrimaryVertex();
1509                 if(!trkVtx || trkVtx->GetNContributors()<=0) return;
1510                 TString vtxTtl = trkVtx->GetTitle();
1511                 if(!vtxTtl.Contains("VertexerTracks")) return;
1512                 Float_t zvtx = -100;
1513                 zvtx=trkVtx->GetZ();
1514
1515
1516                 const AliESDVertex* spdVtx = fESD->GetPrimaryVertexSPD();
1517                 if(spdVtx->GetNContributors()<=0) return;
1518                 TString vtxTyp = spdVtx->GetTitle();
1519                 Double_t cov[6]={0};
1520                 spdVtx->GetCovarianceMatrix(cov);
1521                 Double_t zRes = TMath::Sqrt(cov[5]);
1522                 if(vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1523                 if(TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1524                 if(TMath::Abs(zvtx) > 10) return;
1525         }
1526
1527         //______________________________________________________________________
1528         //after vertex selection
1529         fNevent->Fill(10);
1530
1531         //______________________________________________________________________
1532         //EMCal Trigger Selection (Threshold selection)
1533
1534         TString firedTrigger;
1535         TString TriggerEG1("EG1"); //takes trigger with name with EG1, ex: CEMC7EG1-B-NOPF-CENTNOTRD
1536         TString TriggerEG2("EG2");
1537         //Jan 17, 2014
1538         //TString TriggerEJE("EJE");
1539
1540         if(fAOD) firedTrigger = fAOD->GetFiredTriggerClasses();
1541         else if(fESD) firedTrigger = fESD->GetFiredTriggerClasses();
1542
1543                 //Bool_t IsEventEMCALL0=kTRUE;
1544         Bool_t IsEventEMCALL1=kFALSE;
1545
1546         if(firedTrigger.Contains(TriggerEG1)){
1547                 fNevent->Fill(2);
1548                 IsEventEMCALL1=kTRUE;
1549         }
1550         if(firedTrigger.Contains(TriggerEG2)){
1551                 fNevent->Fill(3);
1552                 IsEventEMCALL1=kTRUE;
1553         }
1554
1555         //if the flag is for a given threshold and it was not fired, return.
1556
1557         if(fEMCEG1){
1558                 if(!firedTrigger.Contains(TriggerEG1))return;
1559                 if(firedTrigger.Contains(TriggerEG2)){
1560                         fNevent->Fill(4);
1561
1562                 }
1563
1564         }
1565
1566
1567         if(fEMCEG2){
1568                 if(!firedTrigger.Contains(TriggerEG2))return;
1569                 if(firedTrigger.Contains(TriggerEG1)){
1570                         fNevent->Fill(5);
1571                 }
1572
1573         }
1574
1575
1576
1577         //______________________________________________________________________
1578         //Testing if there is an overlap EGA and EJE
1579         //none
1580         /*
1581         if(!(firedTrigger.Contains(TriggerEG1) && firedTrigger.Contains(TriggerEG2) ) && !firedTrigger.Contains(TriggerEJE))
1582         {
1583                 fNevent->Fill(6);
1584         }
1585                 //only GA
1586         if((firedTrigger.Contains(TriggerEG1) || firedTrigger.Contains(TriggerEG2)) && !firedTrigger.Contains(TriggerEJE))
1587         {
1588                 fNevent->Fill(7);
1589         }
1590                 //only JE
1591         if(!(firedTrigger.Contains(TriggerEG1) && firedTrigger.Contains(TriggerEG2)) && firedTrigger.Contains(TriggerEJE))
1592         {
1593                 fNevent->Fill(8);
1594         }
1595                 //both
1596         if((firedTrigger.Contains(TriggerEG1) || firedTrigger.Contains(TriggerEG2)) && firedTrigger.Contains(TriggerEJE))
1597         {
1598                 fNevent->Fill(9);
1599         }
1600         */
1601
1602
1603
1604
1605         //Only events with at least 2 tracks are accepted
1606         Int_t fNOtrks =  fVevent->GetNumberOfTracks();
1607         if(fNOtrks<2) return;
1608
1609         fNevent->Fill(11);
1610
1611         if(fIsAOD){
1612                 Int_t fNOtrks2 =  fAOD->GetNumberOfTracks();
1613                 if(fNOtrks2<2) return;
1614         }
1615         fNevent->Fill(12);
1616
1617         //______________________________________________________________________
1618         //new track loop to select events
1619         //track pt cut (at least 2)
1620         /*
1621         if(fUseTrigger){
1622                 if(fIsAOD){
1623                         double fTrackMulti=0;
1624                         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
1625                         {
1626                                 AliVParticle* Vtrack = fVevent->GetTrack(iTracks);
1627                                 if (!Vtrack) continue;
1628
1629                                 AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
1630                                         //AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(Vtrack);
1631
1632                                 if((track->Pt())<0.2 || (track->Pt())>1000.0) continue;
1633                                 else fTrackMulti=fTrackMulti+1;
1634
1635                         }
1636                                 //Only take event if track multiplicity is bigger than 2.
1637                         if(fTrackMulti<2) return;
1638                 }
1639         }
1640         fNevent->Fill(13);
1641         //______________________________________________________________________
1642         //Using more cuts than I have beeing using
1643         //eta cut and primary (at least 2)
1644         if(fUseTrigger){
1645                 if(fIsAOD){
1646                         double fTrackMulti2=0;
1647                         for(Int_t i = 0; i < fVevent->GetNumberOfTracks(); i++)
1648                         {
1649                                 AliVParticle* Vtrack2 = fVevent->GetTrack(i);
1650                                 if (!Vtrack2) continue;
1651
1652
1653                                 AliVTrack *track_new = dynamic_cast<AliVTrack*>(Vtrack2);
1654                                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(Vtrack2);
1655
1656
1657                                 if(aodtrack)
1658                                 {
1659
1660
1661                                     if(TMath::Abs(track_new->Eta())> 0.9) continue;
1662                                         if (aodtrack->GetType()!= AliAODTrack::kPrimary) continue ;
1663                                     else fTrackMulti2=fTrackMulti2+1;
1664                                 }
1665                         }
1666                                 //Only take event if track multiplicity is bigger than 2.
1667                         if(fTrackMulti2<2) return;
1668
1669
1670                 }
1671         }
1672         fNevent->Fill(14);
1673 //______________________________________________________________________
1674 //Using more cuts than I have beeing using
1675 //hybrid (at least2)
1676         if(fUseTrigger){
1677                 if(fIsAOD){
1678                         double fTrackMulti3=0;
1679                         for(Int_t i = 0; i < fVevent->GetNumberOfTracks(); i++)
1680                         {
1681                                 AliVParticle* Vtrack3 = fVevent->GetTrack(i);
1682                                 if (!Vtrack3) continue;
1683
1684                                         //AliVTrack *track_new = dynamic_cast<AliVTrack*>(Vtrack3);
1685                                 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(Vtrack3);
1686
1687
1688                                 if(aodtrack)
1689                                 {
1690
1691                                         if (!aodtrack->IsHybridGlobalConstrainedGlobal()) continue ;
1692                                                 //another option if I don't want to use hybrid
1693                                                 //if ( aodtrack->TestFilterBit(128)==kFALSE) continue ;
1694                                     else fTrackMulti3=fTrackMulti3+1;
1695                                 }
1696                         }
1697                         //Only take event if track multiplicity is bigger than 2.
1698                         if(fTrackMulti3<2) return;
1699
1700                 }
1701         }
1702         fNevent->Fill(15);
1703 //______________________________________________________________________
1704
1705
1706         if(fUseTrigger){
1707                 if(fIsAOD){
1708                         double fTrackMulti4=0;
1709                         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
1710                         {
1711                                 AliVParticle* Vtrack4 = fVevent->GetTrack(iTracks);
1712                                 if (!Vtrack4) continue;
1713
1714
1715                                         //AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack4);
1716                                 AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(Vtrack4);
1717
1718                                 if(!atrack->TestFilterBit(768)) continue;
1719                                 if(!atrack->IsHybridGlobalConstrainedGlobal()) continue ;
1720
1721
1722                                 else fTrackMulti4=fTrackMulti4+1;
1723
1724                         }
1725                         //Only take event if track multiplicity is bigger than 2.
1726                         if(fTrackMulti4<2) return;
1727                         fTrack_Multi->Fill(fTrackMulti4);
1728                 }
1729         }
1730         fNevent->Fill(16);
1731 //______________________________________________________________________
1732         */
1733
1734 //______________________________________________________________________
1735 //Centrality Selection
1736         if(fHasCentralitySelection)
1737         {
1738                 Float_t centrality = -1;
1739
1740                 if(fIsAOD)
1741                 {
1742                         fCentrality = fAOD->GetHeader()->GetCentralityP();
1743                 }
1744                 else
1745                 {
1746                         fCentrality = fESD->GetCentrality();
1747                 }
1748
1749                 if(fEstimator==1) centrality = fCentrality->GetCentralityPercentile("ZDC");
1750                 else centrality = fCentrality->GetCentralityPercentile("V0A");
1751
1752                         //cout << "Centrality = " << centrality << " %" << endl;
1753
1754                 fCentralityHist->Fill(centrality);
1755
1756                 if(centrality<fCentralityMin || centrality>fCentralityMax) return;
1757
1758                 fCentralityHistPass->Fill(centrality);
1759         }
1760         //______________________________________________________________________
1761
1762
1763         fNevent->Fill(17);
1764
1765         //______________________________________________________________________
1766
1767         if(fIsMC)
1768         {
1769                 if(fIsAOD)
1770                 {
1771                         fMCarray = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
1772
1773                         if(!fMCarray)
1774                         {
1775                                 AliError("Array of MC particles not found");
1776                                 return;
1777                         }
1778
1779                         fMCheader = dynamic_cast<AliAODMCHeader*>(fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
1780
1781                         if(!fMCheader)
1782                         {
1783                                 AliError("Could not find MC Header in AOD");
1784                                 return;
1785                         }
1786
1787                         for(Int_t iMC = 0; iMC < fMCarray->GetEntries(); iMC++)
1788                         {
1789                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(iMC);
1790                                 if(fMCparticle->GetMother()>0) fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
1791
1792                                 Int_t pdg = fMCparticle->GetPdgCode();
1793
1794                                 //====================================================================
1795                                 //trying take pions spectra 27/May/2014
1796                                 //IsPrimary only take events from pythia
1797                                 //IsPhysicalPrimariee: (all prompt particles, including strong decay products plus weak decay product from heavy-flavor).
1798                                 //eta cut same as MinJung
1799
1800                                 if(fMCparticle->Eta()>=-0.8 && fMCparticle->Eta()<=0.8)
1801                                 {
1802                                         if(fMCparticle->IsPrimary()){
1803
1804                                                 if(TMath::Abs(pdg)==111) fPtMCpi0->Fill(fMCparticle->Pt());
1805                                                 if(TMath::Abs(pdg)==221) fPtMCeta->Fill(fMCparticle->Pt());
1806                                                 //eta cut same as MinJung
1807                                         }
1808
1809                                         if(fMCparticle->IsPhysicalPrimary()){
1810                                                 if(TMath::Abs(pdg)==111) fPtMCpi02->Fill(fMCparticle->Pt());
1811                                                 if(TMath::Abs(pdg)==221) fPtMCeta2->Fill(fMCparticle->Pt());
1812
1813                                         }
1814
1815                                         if(TMath::Abs(pdg)==111) fPtMCpi03->Fill(fMCparticle->Pt());
1816                                         if(TMath::Abs(pdg)==221) fPtMCeta3->Fill(fMCparticle->Pt());
1817                                 }
1818                                 //====================================================================
1819
1820                                 double proX = fMCparticle->Xv();
1821                                 double proY = fMCparticle->Yv();
1822                                 double proR = sqrt(pow(proX,2)+pow(proY,2));
1823
1824
1825                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax && fMCparticle->Charge()!=0)
1826                                 {
1827                                                 //to correct background
1828                                         if (TMath::Abs(pdg) == 11 && fMCparticle->GetMother()>0){
1829                                                 Int_t mpdg = fMCparticleMother->GetPdgCode();
1830
1831                                                 if(TMath::Abs(mpdg) == 221 || TMath::Abs(mpdg) == 22 || TMath::Abs(mpdg) == 111){
1832
1833                                                         if(proR<7){
1834                                                                 fPtMCparticleAlle_nonPrimary->Fill(fMCparticle->Pt()); //denominator for total efficiency for all electrons, and not primary
1835
1836                                                         }
1837                                                 }
1838                                         }
1839
1840                                         if (TMath::Abs(pdg) == 11 && fMCparticle->IsPhysicalPrimary()){
1841                                                 fPtMCparticleAlle_Primary->Fill(fMCparticle->Pt()); //denominator for total efficiency for all electrons primary
1842                                         }
1843
1844                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1845                                         {
1846
1847                                                 fPtMCparticleAll_nonPrimary->Fill(fMCparticle->Pt()); //denominator for total efficiency for all particles, and not primary
1848                                                 if(fMCparticle->IsPhysicalPrimary())
1849                                                 {
1850                                                         fPtMCparticleAll->Fill(fMCparticle->Pt());
1851
1852                                                         Bool_t MotherFound = FindMother(iMC);
1853                                                         //Bool_t MotherFound = FindMother(track->GetLabel());
1854                                                         if(MotherFound)
1855                                                         {
1856                                                                 if(fIsHFE1){
1857                                                                         //denominator for total efficiency and tracking
1858                                                                         //unfolding: denominator is pt_MC and numerator is pt_reco
1859                                                                         fPtMCparticleAllHfe1->Fill(fMCparticle->Pt());
1860                                                                         fEtaPhi_den->Fill(fMCparticle->Phi(),fMCparticle->Eta());
1861
1862
1863                                                                 } //denominator for total efficiency and tracking
1864                                                                 if(fIsHFE2){
1865                                                                         fPtMCparticleAllHfe2->Fill(fMCparticle->Pt());
1866                                                                 }
1867                                                         }
1868                                                 }
1869                                         }
1870                                 }//eta cut
1871
1872
1873
1874                         }//loop tracks
1875
1876
1877
1878                         //second loop over track, but only the primaries ones
1879                         //only primary pions --> how to take the primaries ones in AOD?
1880                         /*
1881                         for(Int_t iMC = 0; iMC < fMCarray->GetNPrimary(); iMC++){
1882                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(iMC);
1883                                 pdg = fMCparticle->GetPdgCode();
1884
1885                                 if(TMath::Abs(pdg)==111) fPtMCpi0->Fill(fMCparticle->Pt());
1886                                 if(TMath::Abs(pdg)==221) fPtMCeta->Fill(fMCparticle->Pt());
1887
1888                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax)
1889                                 {
1890
1891                                         if(TMath::Abs(pdg)==111) fPtMCpi02->Fill(fMCparticle->Pt());
1892                                         if(TMath::Abs(pdg)==221) fPtMCeta2->Fill(fMCparticle->Pt());
1893
1894                                 }
1895                         }
1896                          */
1897
1898
1899                 }//AOD
1900                 else
1901                 {
1902                         fEventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1903                         if (!fEventHandler) {
1904                                 Printf("ERROR: Could not retrieve MC event handler");
1905                                 return;
1906                         }
1907
1908                         fMCevent = fEventHandler->MCEvent();
1909                         if (!fMCevent) {
1910                                 Printf("ERROR: Could not retrieve MC event");
1911                                 return;
1912                         }
1913
1914                         fMCstack = fMCevent->Stack();
1915
1916                         //pion and eta spectrum
1917                         //MinJung code
1918
1919                         //----------------------------------------------------------------------------------------------------
1920                         AliVParticle *mctrack2 = NULL;
1921                         AliMCParticle *mctrack0 = NULL;
1922
1923
1924                         for(Int_t imc = 0; imc <fMCEvent->GetNumberOfPrimaries(); imc++){
1925                                 if(!(mctrack2 = fMCEvent->GetTrack(imc))) continue;
1926                                 TParticle* mcpart0 = fMCEvent->Stack()->Particle(imc);
1927                                 if(!mcpart0) continue;
1928                                 mctrack0 = dynamic_cast<AliMCParticle *>(mctrack2);
1929                                 if(!mctrack0) continue;
1930
1931                                 if(TMath::Abs(AliHFEtools::GetRapidity(mcpart0))<0.8){
1932
1933                                         if(TMath::Abs(mctrack0->PdgCode()) == 111) // pi0
1934                                         {
1935                                                 fPtMCpi0->Fill(mctrack0->Pt());
1936                                         }
1937
1938                                         if(TMath::Abs(mctrack0->PdgCode()) == 221) // eta
1939                                         {
1940                                                 fPtMCeta->Fill(mctrack0->Pt());
1941                                         }
1942
1943                                 }
1944
1945                         }
1946                         // end of MinJung
1947                         //----------------------------------------------------------------------------------------------------
1948
1949
1950                 for(Int_t iMC = 0; iMC < fMCstack->GetNtrack(); iMC++)
1951                 {
1952
1953                                 fMCtrack = fMCstack->Particle(iMC);
1954                                 if(fMCtrack->GetFirstMother()>0) fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
1955                                 TParticle *particle=fMCstack->Particle(iMC);
1956
1957                                 Int_t pdg = fMCtrack->GetPdgCode();
1958
1959
1960                                 if(TMath::Abs(pdg)==111) fPtMCpi0->Fill(fMCtrack->Pt());
1961                                 if(TMath::Abs(pdg)==221) fPtMCeta->Fill(fMCtrack->Pt());
1962
1963
1964                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
1965                                 {
1966
1967                                                 //to correct background
1968                                         if (TMath::Abs(pdg) == 11 && fMCtrack->GetFirstMother()>0){
1969                                                 Int_t mpdg = fMCtrackMother->GetPdgCode();
1970                                                 if(TMath::Abs(mpdg) == 221 || TMath::Abs(mpdg) == 22 || TMath::Abs(mpdg) == 111){
1971                                                                 Double_t proR=particle->R();
1972                                                                 if(proR<7){
1973                                                                         fPtMCparticleAlle_nonPrimary->Fill(fMCtrack->Pt()); //denominator for total efficiency for all electrons, and not primary
1974                                                                 }
1975                                                 }
1976                                         }
1977
1978                                         if(TMath::Abs(pdg) == 11 && fMCstack->IsPhysicalPrimary(iMC)){
1979
1980                                                 fPtMCparticleAlle_Primary->Fill(fMCtrack->Pt());
1981                                         }
1982
1983                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
1984                                         {
1985                                                 fPtMCparticleAll_nonPrimary->Fill(fMCtrack->Pt());//denominator for total efficiency for all particle, non Primary track
1986
1987                                                 if(fMCstack->IsPhysicalPrimary(iMC))
1988                                                 {
1989                                                         fPtMCparticleAll->Fill(fMCtrack->Pt());
1990
1991                                                         Bool_t MotherFound = FindMother(iMC);
1992                                                         //Bool_t MotherFound = FindMother(track->GetLabel());
1993                                                         if(MotherFound)
1994                                                         {
1995                                                                 if(fIsHFE1){
1996                                                                         fPtMCparticleAllHfe1->Fill(fMCtrack->Pt());//denominator for total efficiency and tracking
1997                                                                         fEtaPhi_den->Fill(fMCtrack->Phi(),fMCtrack->Eta());
1998                                                                 }
1999                                                                 if(fIsHFE2){
2000                                                                         fPtMCparticleAllHfe2->Fill(fMCtrack->Pt());
2001                                                                 }
2002                                                         }
2003                                                 }//Is Physical primary
2004                                         }
2005                                 }//eta cut
2006                 }//loop tracks
2007                 }//ESD
2008         }//Is MC
2009
2010 //______________________________________________________________________
2011 //threshold selection was here
2012 //______________________________________________________________________
2013 //all events selected
2014
2015         fNevent->Fill(0);
2016
2017
2018 //______________________________________________________________________
2019 //events in the threshold
2020
2021         if(firedTrigger.Contains(TriggerEG1))
2022         {
2023                 if(fEMCEG1){
2024                         fNevent->Fill(18);
2025                     if(!firedTrigger.Contains(TriggerEG2)) fNevent->Fill(19);
2026                                 //if(firedTrigger.Contains(TriggerEG2)) return;
2027                 }
2028         }
2029
2030
2031         //EG2
2032         if(firedTrigger.Contains(TriggerEG2))
2033         {
2034                 if(fEMCEG2){
2035                         fNevent->Fill(20);
2036                         if(!firedTrigger.Contains(TriggerEG1)) fNevent->Fill(21);
2037                                 //if(firedTrigger.Contains(TriggerEG1)) return;
2038                 }
2039         }
2040
2041
2042
2043         //New cluster information
2044         //after trigger threshold selection
2045         Int_t ClsNo2 = -999;
2046         ClsNo2 = fVevent->GetNumberOfCaloClusters();
2047         fNCluster_pure->Fill(ClsNo2);
2048
2049
2050
2051         if(ClsNo2<=0){
2052                 fNevent->Fill(22); //events with no cluster
2053                 return;
2054         }
2055
2056         //in order to include time cut
2057         //fEMCALCells = fAOD->GetEMCALCells();
2058         //Double_t tof = clus->GetTOF();
2059         //clus->GetNCells()
2060         //if ( clus->E() < minE ) continue ;
2061
2062
2063
2064         if(fUseTrigger && fIsAOD){
2065
2066           //AliAODHeader * aodh = fAOD->GetHeader();
2067           //Int_t bc= aodh->GetBunchCrossNumber();
2068
2069
2070                 Int_t ClsNo_aod = -999;
2071                 ClsNo_aod = fAOD->GetNumberOfCaloClusters();
2072                 fNCluster_pure_aod->Fill(ClsNo_aod);
2073                         //Bool_t exotic=kTRUE;
2074
2075                 for (Int_t i=0; i< ClsNo_aod; i++ ){
2076
2077                         //fClus = fVevent->GetCaloCluster(i);
2078                         //to be compatible with Shingo
2079                     AliVCluster *clust = 0x0;
2080                         clust = (AliVCluster*) fAOD->GetCaloCluster(i);
2081
2082                         if(clust && clust->IsEMCAL())
2083                         {
2084                                 //pure cluster information
2085                                 fECluster_pure->Fill(clust->E());
2086
2087                                 fNcells_energy->Fill(clust->GetNCells(),clust->E());
2088                                 fNCluster_ECluster->Fill(ClsNo_aod,clust->E());
2089
2090                                 if(clust->E()>1000) fNevent->Fill(23);
2091
2092                                 //exotic
2093                                 /*
2094                                 exotic   = fEMCALRecoUtils->IsExoticCluster(clust, (AliVCaloCells*)fAOD->GetEMCALCells(), bc);
2095                                 if(exotic == kFALSE){
2096                                         fECluster_not_exotic->Fill(clust->E());
2097                                         fNcells_energy_not_exotic->Fill(clust->GetNCells(),clust->E());
2098                                 }
2099                                 */
2100
2101                                 //approximation to remove exotics
2102                                 if(clust->GetNCells()<5 && clust->E()>15.0){
2103                                         fECluster_exotic->Fill(clust->E());
2104                                 }
2105                                 //Marcel cut
2106                                 else if((clust->GetNCells())> ((clust->E())/3+1)){
2107                                         fECluster_not_exotic1->Fill(clust->E());
2108                                 }
2109                                 else{
2110                                         fECluster_not_exotic2->Fill(clust->E());
2111                                 }
2112
2113
2114                         }
2115                         /*
2116                         //______________________________________________________________________
2117                         //Trying to remove events with bad cells and find patches
2118                         //First, I will try to count them
2119                         //______________________________________________________________________
2120
2121                         if(clust && clust->IsEMCAL())
2122                         {
2123                                 Bool_t badchannel = ContainsBadChannel("EMCAL", clust->GetCellsAbsId(),clust->GetNCells() );
2124                                 printf("Contém bad channel? %d ", badchannel);
2125                                 if(badchannel)fNevent2->Fill(0);
2126
2127                                 //trying to find patches
2128                                 TArrayI patches_found=GetTriggerPatches(IsEventEMCALL0, IsEventEMCALL1);
2129                                 printf("N patches %d, first %d, last %d\n",patches_found.GetSize(),  patches_found.At(0), patches_found.At(patches_found.GetSize()-1));
2130
2131                         }
2132
2133                         //end of bad cells
2134                         //______________________________________________________________________
2135 */
2136
2137                 }
2138         }
2139
2140
2141         fNevent->Fill(24); //events with cluster
2142
2143
2144         fVtxZ_new4->Fill(fZvtx);
2145
2146
2147         Int_t ClsNo = -999;
2148         if(!fIsAOD) ClsNo = fESD->GetNumberOfCaloClusters();
2149         else ClsNo = fAOD->GetNumberOfCaloClusters();
2150
2151         //______________________________________________________________________
2152
2153         ///_____________________________________________________________________
2154         ///Track loop
2155         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
2156         {
2157                 AliVParticle* Vtrack = fVevent->GetTrack(iTracks);
2158                 if (!Vtrack)
2159                 {
2160                         printf("ERROR: Could not receive track %d\n", iTracks);
2161                         continue;
2162                 }
2163
2164                 AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
2165                 AliESDtrack *etrack = dynamic_cast<AliESDtrack*>(Vtrack);
2166                 AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(Vtrack);
2167
2168
2169                 Double_t fTPCnSigma = -999;
2170                 Double_t fTOFnSigma = -999;
2171                 Double_t fTPCnSigma_pion = -999;
2172                 Double_t fTPCnSigma_proton = -999;
2173                 Double_t fTPCnSigma_kaon = -999;
2174                 Double_t fTPCsignal = -999;
2175                 Double_t fPt = -999;
2176                 Double_t fP = -999;
2177
2178                 //December 9th 2013
2179                 //Etacut test on the begging
2180                 fEtad[0]->Fill(track->Eta());
2181                         //if(track->Eta()<fEtaCutMin || track->Eta()>fEtaCutMax) continue;
2182                 fEtad[1]->Fill(track->Eta());
2183
2184
2185
2186
2187                 ///_____________________________________________________________________________
2188                 ///Fill QA plots without track selection
2189                 fPt = track->Pt();
2190                 fP = TMath::Sqrt((track->Pt())*(track->Pt()) + (track->Pz())*(track->Pz()));
2191
2192                 fTPCsignal = track->GetTPCsignal();
2193                 fTPCnSigma = fPidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2194                 fTOFnSigma = fPidResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
2195                 fTPCnSigma_pion = fPidResponse->NumberOfSigmasTPC(track, AliPID::kPion);
2196                 fTPCnSigma_proton = fPidResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2197                 fTPCnSigma_kaon = fPidResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
2198
2199                 fTPC_p[0]->Fill(fPt,fTPCsignal);
2200                 fTPCnsigma_p[0]->Fill(fP,fTPCnSigma);
2201
2202
2203                 Float_t TPCNcls = track->GetTPCNcls();
2204                 //TPC Ncls for pid
2205                 Float_t TPCNcls_pid = track->GetTPCsignalN();
2206
2207
2208
2209                 Float_t pos[3]={0,0,0};
2210
2211                 Double_t fEMCflag = kFALSE;
2212                 if(track->GetEMCALcluster()>0)
2213                 {
2214                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
2215                         if(fClus->IsEMCAL())
2216                         {
2217
2218                                 //only for charged tracks
2219                                 fECluster[0]->Fill(fClus->E());
2220
2221
2222                                 if(TMath::Abs(fClus->GetTrackDx())<=fdPhiCut && TMath::Abs(fClus->GetTrackDz())<=fdEtaCut)
2223                                 {
2224                                         fEMCflag = kTRUE;
2225                                         fEoverP_pt[0]->Fill(fPt,(fClus->E() / fP));
2226
2227
2228                                         Float_t Energy  = fClus->E();
2229                                         Float_t EoverP  = Energy/track->P();
2230                                                 //Float_t M02   = fClus->GetM02();
2231                                                 //Float_t M20   = fClus->GetM20();
2232
2233                                                 /////////////// for Eta Phi distribution
2234                                         fClus->GetPosition(pos);
2235                                         TVector3 vpos(pos[0],pos[1],pos[2]);
2236                                         Double_t cphi = vpos.Phi();
2237                                         Double_t ceta = vpos.Eta();
2238                                         fEtaPhi[0]->Fill(cphi,ceta);
2239
2240
2241                                         fTPCNcls_EoverP[0]->Fill(TPCNcls, EoverP);
2242                                 }
2243                         }
2244                 }
2245
2246                         //______________________________________________________________
2247                         // Vertex
2248
2249                 fVtxZ[0]->Fill(fZvtx);
2250                 if(iTracks == 0)fNTracks[0]->Fill(fNOtrks);
2251                 fNTracks_pt[0]->Fill(fNOtrks, fPt);
2252                 fNTracks_eta[0]->Fill(fNOtrks, track->Eta());
2253                 fNTracks_phi[0]->Fill(fNOtrks, track->Phi());
2254
2255
2256                 fNClusters[0]->Fill(ClsNo);
2257                 fTPCNcls_pid[0]->Fill(TPCNcls, TPCNcls_pid);
2258                         //______________________________________________________________
2259
2260 ///Fill QA plots without track selection
2261 ///_____________________________________________________________________________
2262 //______________________________________________________________________________________
2263 //Track Selection Cuts
2264
2265 //AOD (Test Filter Bit)
2266                 if(fIsAOD)
2267                 {
2268                                 // standard cuts with very loose DCA - BIT(4)
2269                                 // Description:
2270                         /*
2271                          GetStandardITSTPCTrackCuts2011(kFALSE)
2272                          SetMaxChi2PerClusterTPC(4);
2273                          SetAcceptKinkDaughters(kFALSE);
2274                          SetRequireTPCRefit(kTRUE);
2275                          SetRequireITSRefit(kTRUE);
2276                          SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
2277                          SetMaxDCAToVertexZ(2);
2278                          SetMaxDCAToVertex2D(kFALSE);
2279                          SetRequireSigmaToVertex(kFALSE);
2280                          SetMaxChi2PerClusterITS(36);
2281                          SetMaxDCAToVertexXY(2.4)
2282                          SetMaxDCAToVertexZ(3.2)
2283                          SetDCaToVertex2D(kTRUE)
2284                          */
2285
2286                         if(!atrack->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;
2287                 }
2288
2289 //RecKine: ITSTPC cuts
2290                 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
2291 //RecKink
2292                 if(fRejectKinkMother)
2293                 {
2294                         if(fIsAOD)
2295                         {
2296                                 Bool_t kinkmotherpass = kTRUE;
2297                                 for(Int_t kinkmother = 0; kinkmother < fNumberOfMotherkink; kinkmother++)
2298                                 {
2299                                         if(track->GetID() == fListOfmotherkink[kinkmother])
2300                                         {
2301                                                 kinkmotherpass = kFALSE;
2302                                                 continue;
2303                                         }
2304                                 }
2305                                 if(!kinkmotherpass) continue;
2306                         }
2307                         else
2308                         {
2309                                 if(etrack->GetKinkIndex(0) != 0) continue;
2310                         }
2311                 }
2312
2313 //RecPrim
2314                 if(!fIsAOD)
2315                 {
2316                         if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
2317                 }
2318
2319 //HFEcuts: ITS layers cuts
2320                 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
2321
2322 //HFE cuts: TPC PID cleanup
2323                 if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
2324 //______________________________________________________________________________________
2325
2326                 if(fIsAOD){
2327                                 //AOD test -- Fancesco suggestion
2328                                 //aod test -- Francesco suggestion
2329                         AliAODTrack *aod_track=fAOD->GetTrack(iTracks);
2330
2331                         Int_t type=aod_track->GetType();
2332                         if(type==AliAODTrack::kPrimary) fPtPrim->Fill(aod_track->Pt());
2333                         if(type==AliAODTrack::kSecondary) fPtSec->Fill(aod_track->Pt());
2334
2335                                 //Int_t type2=track->GetType();
2336                         if(type==AliAODTrack::kPrimary) fPtPrim2->Fill(track->Pt());
2337                         if(type==AliAODTrack::kSecondary) fPtSec2->Fill(track->Pt());
2338                 }
2339
2340
2341 ///_____________________________________________________________
2342 ///QA plots after track selection
2343                 if(fIsMC)
2344                 {
2345                         if(track->GetLabel()>=0) fPtMCWithLabel->Fill(fPt);
2346                         else fPtMCWithoutLabel->Fill(fPt);
2347                 }
2348
2349                 if(fIsMC  && track->GetLabel()>=0)
2350                 {
2351                         if(fIsAOD)
2352                         {
2353                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2354
2355
2356                                 if(fMCparticle->IsPhysicalPrimary())
2357                                 {
2358                                         fPtIsPhysicaPrimary->Fill(fPt);
2359                                 }
2360
2361                                 Int_t pdg = fMCparticle->GetPdgCode();
2362                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax && fMCparticle->Charge()!=0)
2363                                 {
2364
2365
2366                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
2367                                         {
2368                                                 fPtMCparticleReco_nonPrimary->Fill(fMCparticle->Pt()); //not Primary track
2369
2370                                                 if(fMCparticle->IsPhysicalPrimary())
2371                                                 {
2372                                                         fPtMCparticleReco->Fill(fMCparticle->Pt());
2373
2374                                                         Bool_t MotherFound = FindMother(track->GetLabel());
2375                                                         if(MotherFound)
2376                                                         {
2377                                                                 if(fIsHFE1)
2378                                                                 {
2379                                                                         fPtMCparticleRecoHfe1->Fill(fMCparticle->Pt());//numerator tracking
2380                                                                         //unfolding
2381                                                                         fpt_reco_pt_MC_den->Fill(track->Pt(),fMCparticle->Pt());
2382
2383                                                                 }
2384                                                                 if(fIsHFE2){
2385                                                                         fPtMCparticleRecoHfe2->Fill(fMCparticle->Pt());
2386                                                                 }
2387                                                         }
2388                                                 }
2389                                         }
2390                                 }
2391
2392
2393
2394                         }//close AOD
2395                          //ESD
2396                         else
2397                         {
2398
2399
2400                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
2401                                 {
2402
2403                                         fMCtrack = fMCstack->Particle(track->GetLabel());
2404                                         Int_t pdg = fMCtrack->GetPdgCode();
2405
2406                                         if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
2407                                         {
2408                                                 fPtMCparticleReco_nonPrimary->Fill(fMCtrack->Pt());//not Primary track
2409                                         }
2410
2411
2412                                         if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
2413                                         {
2414                                                 fPtIsPhysicaPrimary->Fill(fPt);
2415
2416
2417                                                 if( TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13 )
2418                                                 {
2419                                                         fPtMCparticleReco->Fill(fMCtrack->Pt());
2420
2421                                                         Bool_t MotherFound = FindMother(track->GetLabel());
2422                                                         if(MotherFound)
2423                                                         {
2424                                                                 if(fIsHFE1) fPtMCparticleRecoHfe1->Fill(fMCtrack->Pt());//numerator tracking
2425                                                                 if(fIsHFE2) fPtMCparticleRecoHfe2->Fill(fMCtrack->Pt());
2426                                                         }
2427                                                 }
2428                                         }
2429                                 }
2430                         }//close ESD
2431                 }//close IsMC
2432
2433                 fTPC_p[1]->Fill(fPt,fTPCsignal);
2434                 fTPCnsigma_p[1]->Fill(fP,fTPCnSigma);
2435                 Double_t fPtBin[7] = {1,2,4,6,8,10,15};
2436
2437                 TPCNcls = track->GetTPCNcls();
2438                 Float_t pos2[3]={0,0,0};
2439
2440                 if(track->GetEMCALcluster()>0)
2441                 {
2442                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
2443                         if(fClus->IsEMCAL())
2444                         {
2445                                 if(TMath::Abs(fClus->GetTrackDx())<=fdPhiCut && TMath::Abs(fClus->GetTrackDz())<=fdEtaCut)
2446                                 {
2447                                         fEoverP_pt[1]->Fill(fPt,(fClus->E() / fP));
2448
2449                                         Float_t Energy  = fClus->E();
2450                                         Float_t EoverP  = Energy/track->P();
2451                                         Float_t M02     = fClus->GetM02();
2452                                         Float_t M20     = fClus->GetM20();
2453                                         Float_t ncells  = fClus->GetNCells();
2454
2455                                                 /////////////// for Eta Phi distribution
2456                                         fClus->GetPosition(pos2);
2457                                         TVector3 vpos(pos2[0],pos2[1],pos2[2]);
2458                                         Double_t cphi = vpos.Phi();
2459                                         Double_t ceta = vpos.Eta();
2460                                         fEtaPhi[1]->Fill(cphi,ceta);
2461
2462                                         fECluster[1]->Fill(Energy);
2463                                         fTPCNcls_EoverP[1]->Fill(TPCNcls, EoverP);
2464
2465
2466                                         //for EMCal trigger performance
2467                                         if(EoverP > 0.9){
2468                                                 ftpc_p_EoverPcut->Fill(track->P(), fTPCsignal);
2469                                                 fnsigma_p_EoverPcut->Fill(track->P(), fTPCnSigma);
2470
2471                                         }
2472
2473
2474                                         //for hadron contamination calculations
2475
2476
2477                                         // EtaCut -> dados
2478                                         if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
2479                                                         //main
2480                                                 if(TMath::Abs(fTPCnSigma_pion)<3 || TMath::Abs(fTPCnSigma_proton)<3 || TMath::Abs(fTPCnSigma_kaon)<3 ){
2481
2482                                                         if(fTPCnSigma<-3.5){
2483                                                                 fEoverP_pt_hadrons->Fill(fPt,EoverP);
2484                                                                 if(fUseEMCal) fShowerShape_ha->Fill(M02,M20);
2485                                                         }
2486                                                 }
2487                                                         //for systematic studies of hadron contamination
2488                                                 if(fTPCnSigma < -4){
2489                                                         fEoverP_pt_pions->Fill(fPt, EoverP);
2490
2491                                                 }
2492
2493                                                 if(fTPCnSigma < -5){
2494                                                         fEoverP_pt_pions2->Fill(fPt, EoverP);
2495
2496                                                 }
2497
2498
2499                                         }
2500
2501
2502
2503
2504                                         for(Int_t i = 0; i < 6; i++)
2505                                         {
2506                                                 if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
2507                                                 {
2508
2509                                                         if(fTPCnSigma < -5 && fTPCnSigma > -10){
2510                                                                 fNcells_hadrons[i]->Fill(ncells);
2511                                                         }
2512                                                                 //hadrons selection using E/p
2513                                                         if((fClus->E() / fP) >= 0.2 && (fClus->E() / fP) <= 0.7){
2514                                                                 fTPCnsigma_eta_hadrons[i]->Fill(fTPCnSigma, ceta);
2515                                                         }
2516                                                                 //electrons selection using E/p
2517                                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax) {
2518                                                                 fTPCnsigma_eta_electrons[i]->Fill(fTPCnSigma, ceta);
2519                                                         }
2520                                                 }
2521                                         }
2522
2523                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax)
2524                                         {
2525                                                 fTPCnsigma_eta->Fill(track->Eta(),fTPCnSigma);
2526                                                 fTPCnsigma_phi->Fill(track->Phi(),fTPCnSigma);
2527
2528                                                 if(fUseEMCal)
2529                                                 {
2530                                                         if(fTPCnSigma < 3.5 && fCorrelationFlag)
2531                                                         {
2532                                                                 fPtTrigger_Inc->Fill(fPt);
2533                                                                 DiHadronCorrelation(track, iTracks);
2534                                                         }
2535                                                 }
2536                                         }
2537
2538                                 }
2539                         }
2540                 }
2541
2542                 //______________________________________________________________
2543                 // Vertex
2544
2545                 fVtxZ[1]->Fill(fZvtx);
2546                 if(iTracks == 0)fNTracks[1]->Fill(fNOtrks);
2547                 fNTracks_pt[1]->Fill(fNOtrks, fPt);
2548                 fNTracks_eta[1]->Fill(fNOtrks, track->Eta());
2549                 fNTracks_phi[1]->Fill(fNOtrks, track->Phi());
2550                 fNClusters[1]->Fill(ClsNo);
2551                 fTPCNcls_pid[1]->Fill(TPCNcls, TPCNcls_pid);
2552
2553                 //______________________________________________________________
2554
2555                 ///______________________________________________________________________
2556                 ///Histograms for PID Studies
2557                 //Double_t fPtBin[6] = {2,4,6,8,10,15};
2558
2559                 for(Int_t i = 0; i < 6; i++)
2560                 {
2561                         if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
2562                         {
2563                                 if(fEMCflag) fEoverP_tpc[i]->Fill(fTPCnSigma,(fClus->E() / fP));
2564
2565
2566                                 fTPC_pt[i]->Fill(fTPCsignal);
2567                                 fTPCnsigma_pt[i]->Fill(fTPCnSigma);
2568
2569                         }
2570                 }
2571
2572                         ///QA plots after track selection
2573                         ///_____________________________________________________________
2574
2575                         //_______________________________________________________
2576                         //Correlation Analysis - DiHadron
2577                 if(!fUseEMCal)
2578                 {
2579                         if(fTPCnSigma < 3.5 && fCorrelationFlag)
2580                         {
2581                                 fPtTrigger_Inc->Fill(fPt);
2582                                 DiHadronCorrelation(track, iTracks);
2583                         }
2584                 }
2585                         //_______________________________________________________
2586
2587
2588                         ///////////////////////////////////////////////////////////////////
2589                         ///TPC - efficiency calculation //
2590
2591                         /// changing start here
2592                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
2593                 {
2594                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2595                         Int_t pdg = fMCparticle->GetPdgCode();
2596
2597                                 //
2598                         if(fMCparticle->IsPhysicalPrimary()){
2599
2600
2601                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
2602
2603                                         Bool_t MotherFound = FindMother(track->GetLabel());
2604                                         if(MotherFound){
2605                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2606                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
2607                                                         if(fIsHFE1) fPtMC_TPC_All->Fill(fMCparticle->Pt());
2608                                         }
2609                                         }
2610                                 }
2611                         }
2612                 }///until here
2613
2614                 else if(fIsMC && track->GetLabel()>=0)//ESD
2615                 {
2616
2617                         if(fMCstack->IsPhysicalPrimary(track->GetLabel())){
2618                                 fMCtrack = fMCstack->Particle(track->GetLabel());
2619
2620                                 Int_t pdg = fMCtrack->GetPdgCode();
2621                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax ){
2622
2623                                         if(fMCtrack->GetFirstMother()>0){
2624                                             fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
2625                                                 Bool_t MotherFound = FindMother(track->GetLabel());
2626                                                 if(MotherFound){
2627                                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
2628                                                                 if(fIsHFE1) fPtMC_TPC_All->Fill(fMCtrack->Pt());
2629                                                         }
2630                                                 }
2631                                         }
2632                                 }
2633                         }
2634                 }
2635
2636
2637                         if(fPt>1 && fPt<2) fTOF01->Fill(fTOFnSigma,fTPCnSigma);
2638                         if(fPt>2 && fPt<4) fTOF02->Fill(fTOFnSigma,fTPCnSigma);
2639                         if(fPt>4 && fPt<6) fTOF03->Fill(fTOFnSigma,fTPCnSigma);
2640
2641 ///________________________________________________________________________
2642 ///PID
2643 ///Here the PID cuts defined in the file "ConfigEMCalHFEpA.C" is applied
2644             Int_t pidpassed = 1;
2645             AliHFEpidObject hfetrack;
2646             hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
2647             hfetrack.SetRecTrack(track);
2648             hfetrack.SetPP();   //proton-proton analysis
2649             if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
2650             fpid->Fill(pidpassed);
2651
2652             if(pidpassed==0) continue;
2653 ///________________________________________________________________________
2654
2655
2656 ////////////////////////////////////////////////////////////////////
2657 ///TPC efficiency calculations
2658
2659
2660                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
2661                 {
2662                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2663                         Int_t pdg = fMCparticle->GetPdgCode();
2664
2665                                 //
2666                         if(fMCparticle->IsPhysicalPrimary()){
2667
2668
2669                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
2670
2671                                         Bool_t MotherFound = FindMother(track->GetLabel());
2672                                         if(MotherFound){
2673                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2674                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
2675                                                         if(fIsHFE1) fPtMC_TPC_Selected->Fill(fMCparticle->Pt());
2676                                         }
2677                                         }
2678                                 }
2679                         }
2680                 }///until here
2681
2682                 else if(fIsMC && track->GetLabel()>=0)//ESD
2683                 {
2684
2685                         if(fMCstack->IsPhysicalPrimary(track->GetLabel())){
2686                                 fMCtrack = fMCstack->Particle(track->GetLabel());
2687
2688                                 Int_t pdg = fMCtrack->GetPdgCode();
2689                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax ){
2690
2691                                         if(fMCtrack->GetFirstMother()>0){
2692                                             fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
2693                                                 Bool_t MotherFound = FindMother(track->GetLabel());
2694                                                 if(MotherFound){
2695                                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
2696                                                                 if(fIsHFE1) fPtMC_TPC_Selected->Fill(fMCtrack->Pt());
2697                                                         }
2698                                                 }
2699                                         }
2700                                 }
2701                         }
2702                 }
2703
2704                 //Eta Cut for TPC only
2705                 if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
2706                         fTPCnsigma_pt_2D->Fill(fPt,fTPCnSigma);
2707                 }
2708
2709                 //Background for TPC only
2710                 if(fFillBackground){
2711
2712                         //efficiency without SS cut for TPC only
2713                         if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax){
2714                                 Background(track, iTracks, Vtrack, kTRUE); //IsTPConly=kTRUE
2715                         } //Eta cut to be consistent with other efficiency
2716                 }
2717
2718
2719                 fTPCnsigma_p[2]->Fill(fP,fTPCnSigma);
2720                 fTPC_p[2]->Fill(fP,fTPCsignal);
2721                 TPCNcls = track->GetTPCNcls();
2722                 Float_t pos3[3]={0,0,0};
2723
2724
2725                 //here denominator for track-matching efficiency
2726                 if(fIsMC && fIsAOD && track->GetLabel()>=0)
2727                 {
2728                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2729                         Int_t pdg = fMCparticle->GetPdgCode();
2730
2731                                 //
2732                         if(fMCparticle->IsPhysicalPrimary()){
2733
2734
2735                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
2736
2737                                         Bool_t MotherFound = FindMother(track->GetLabel());
2738                                         if(MotherFound){
2739                                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2740                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
2741                                                         if(fIsHFE1) fPt_track_match_den->Fill(fMCparticle->Pt());
2742                                         }
2743                                         }
2744                                 }
2745                         }
2746                 }///until here
2747
2748
2749                 if(track->GetEMCALcluster()>0)
2750                 {
2751                         fClus = fVevent->GetCaloCluster(track->GetEMCALcluster());
2752                         if(fClus->IsEMCAL())
2753                         {
2754
2755                 //________________________________________________________________________
2756
2757
2758                                 //Cluster timing distribution -- for ESD
2759                                 if(fUseEMCal && !fIsAOD){
2760
2761                                         AliESDCaloCells &cells_esd=*(fESD->GetEMCALCells());
2762                                         TRefArray* caloClusters_esd = new TRefArray();
2763                                     fESD->GetEMCALClusters(caloClusters_esd);
2764                                         Int_t nclus_esd = caloClusters_esd->GetEntries();
2765
2766
2767                                         for (Int_t icl = 0; icl < nclus_esd; icl++) {
2768
2769                                                 AliESDCaloCluster* clus_esd = (AliESDCaloCluster*)caloClusters_esd->At(icl);
2770
2771                                                 if(clus_esd->IsEMCAL()){
2772                                                         Float_t ncells_esd      = fClus->GetNCells();
2773                                                         UShort_t * index_esd = clus_esd->GetCellsAbsId() ;
2774                                                         UShort_t * index2_esd = fClus->GetCellsAbsId() ;
2775
2776
2777                                                         for(Int_t i = 0; i < ncells_esd ; i++){
2778
2779                                                                 Int_t absId_esd =   index_esd[i];
2780                                                                 fTime->Fill(fPt,cells_esd.GetCellTime(absId_esd));
2781
2782                                                                 Int_t absId2_esd =   index2_esd[i];
2783                                                                 fTime2->Fill(fPt,cells_esd.GetCellTime(absId2_esd));
2784                                                         }
2785
2786                                                 }
2787                                         }
2788
2789                                 }
2790                                 /* not working!
2791                                 //Cluster timing distribution -- for AOD
2792                                 if(fUseEMCal && fIsAOD){
2793
2794                                         AliAODCaloCells &cells_aod=*(fAOD->GetEMCALCells());
2795
2796                                         TRefArray* caloClusters_aod = new TRefArray();
2797                                         fAOD->GetEMCALClusters(caloClusters_aod);
2798
2799                                         Int_t nclus_aod = caloClusters_aod->GetEntries();
2800
2801                                         for (Int_t icl = 0; icl < nclus_aod; icl++) {
2802
2803                                                 AliAODCaloCluster* clus_aod = (AliAODCaloCluster*)caloClusters_aod->At(icl);
2804
2805
2806                                                 if(clus_aod->IsEMCAL()){
2807                                                         Float_t ncells_aod      = fClus->GetNCells();
2808                                                         UShort_t * index_aod = clus_aod->GetCellsAbsId() ;
2809                                                         UShort_t * index2_aod = fClus->GetCellsAbsId() ;
2810
2811
2812                                                         for(Int_t i = 0; i < ncells_aod ; i++){
2813
2814                                                                 Int_t absId_aod =   index_aod[i];
2815                                                                 fTime->Fill(fPt,cells_aod.GetCellTime(absId_aod));
2816
2817                                                                 Int_t absId2_aod =   index2_aod[i];
2818                                                                 fTime2->Fill(fPt,cells_aod.GetCellTime(absId2_aod));
2819                                                         }
2820
2821                                                 }
2822                                         }
2823
2824                                 }
2825                                  */
2826
2827
2828                                 if(fUseEMCal){
2829                                         double emctof = fClus->GetTOF();
2830                                         ftimingEle->Fill(fPt,emctof);
2831                                 }
2832                                         //________________________________________________________________________
2833
2834
2835
2836
2837                                 // Residuals
2838                                 Double_t Dx = fClus->GetTrackDx();
2839                                 Double_t Dz = fClus->GetTrackDz();
2840                                 Double_t R=TMath::Sqrt(Dx*Dx+Dz*Dz);
2841
2842                                 for(Int_t i = 0; i < 6; i++)
2843                                 {
2844                                         if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
2845                                         {
2846
2847                                                 fEta[i]->Fill(Dz);
2848                                                 fPhi[i]->Fill(Dx);
2849                                                 fR[i]->Fill(R);
2850                                         }
2851                                 }
2852
2853                                 if(TMath::Abs(fClus->GetTrackDx())<=fdPhiCut && TMath::Abs(fClus->GetTrackDz())<=fdEtaCut)
2854                                 {
2855                                         Float_t Energy  = fClus->E();
2856                                         Float_t EoverP  = Energy/track->P();
2857                                         Float_t M02     = fClus->GetM02();
2858                                         Float_t M20     = fClus->GetM20();
2859                                         Float_t ncells  = fClus->GetNCells();
2860
2861                                         //here numerator for track-matching efficiency
2862                                         if(fIsMC && fIsAOD && track->GetLabel()>=0)
2863                                         {
2864                                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2865                                                 Int_t pdg = fMCparticle->GetPdgCode();
2866
2867                                                         //
2868                                                 if(fMCparticle->IsPhysicalPrimary()){
2869
2870
2871                                                         if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
2872
2873                                                                 Bool_t MotherFound = FindMother(track->GetLabel());
2874                                                                 if(MotherFound){
2875                                                                         fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2876                                                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
2877                                                                                 if(fIsHFE1) fPt_track_match_num->Fill(fMCparticle->Pt());
2878                                                                         }
2879                                                                 }
2880                                                         }
2881                                                 }
2882                                         }///until here
2883
2884
2885
2886                                         //----------------------------------------------------------------------------------------
2887                                         //
2888                                         //EtaCut electrons histogram
2889                                         //Shower Shape Cut
2890                                         if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
2891
2892                                                 if(fUseShowerShapeCut){
2893                                                         if(M02 >= fM02CutMin && M02<=fM02CutMax && M20>=fM20CutMin && M20<=fM20CutMax){
2894                                                                 fEoverP_pt[2]->Fill(fPt,(fClus->E() / fP));
2895                                                                 fShowerShapeCut->Fill(M02,M20);
2896                                                                         //in order to check if there are exotic cluster in this selected cluster (27 may 2014)
2897                                                                 fNcells_energy_elec_selected->Fill(ncells,Energy);
2898
2899                                                         }
2900
2901                                                 }
2902                                                 if(!fUseShowerShapeCut){
2903                                                         fEoverP_pt[2]->Fill(fPt,(fClus->E() / fP));
2904                                                         fShowerShapeCut->Fill(M02,M20);
2905                                                         fNcells_energy_elec_selected->Fill(ncells,Energy);
2906
2907
2908                                                 }
2909                                                 if(fUseEMCal) fShowerShape_ele->Fill(M02,M20);
2910
2911                                                 //for shower shape cut studies - now with TPC PID Cut
2912                                                 if(fUseEMCal){
2913                                                         fShowerShapeM02_EoverP->Fill(M02,EoverP);
2914                                                         fShowerShapeM20_EoverP->Fill(M20,EoverP);
2915                                                 }
2916
2917                                         }
2918
2919                                         //----------------------------------------------------------------------------------------
2920
2921
2922
2923                                         // for Eta Phi distribution
2924                                         fClus->GetPosition(pos3);
2925                                         TVector3 vpos(pos3[0],pos3[1],pos3[2]);
2926                                         Double_t cphi = vpos.Phi();
2927                                         Double_t ceta = vpos.Eta();
2928                                         fEtaPhi[2]->Fill(cphi,ceta);
2929
2930
2931
2932                                         fTPCNcls_EoverP[2]->Fill(TPCNcls, EoverP);
2933
2934                                         for(Int_t i = 0; i < 6; i++)
2935                                         {
2936                                                 if(fPt>=fPtBin[i] && fPt<fPtBin[i+1])
2937                                                 {
2938
2939                                                         fR_EoverP[i]->Fill(R, EoverP);
2940                                                         fNcells[i]->Fill(ncells);
2941                                                         fNcells_EoverP[i]->Fill(EoverP, ncells);
2942                                                         fM02_EoverP[i]->Fill(M02,EoverP);
2943                                                         fM20_EoverP[i]->Fill(M20,EoverP);
2944                                                         fECluster_ptbins[i]->Fill(Energy);
2945                                                         fEoverP_ptbins[i]->Fill(EoverP);
2946
2947                                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax) {
2948                                                                 fNcells_electrons[i]->Fill(ncells);
2949                                                         }
2950
2951                                                         if(M02<0.5 && M20<0.3) {
2952                                                                 fEoverP_wSSCut[i]->Fill(EoverP);
2953                                                         }
2954                                                 }
2955                                         }
2956
2957                                         fNcells_pt->Fill(fPt, ncells);
2958
2959
2960                                         ////////////////////////////////////////////////////////////////////
2961                                         ///EMCal - Efficiency calculations
2962
2963
2964                                         if(fIsMC && fIsAOD && track->GetLabel()>=0)
2965                                         {
2966                                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
2967                                                 Int_t pdg = fMCparticle->GetPdgCode();
2968
2969                                                 //
2970                                                 if(fMCparticle->IsPhysicalPrimary()){
2971
2972                                                         //Phi cut && fMCparticle->Phi()>=(TMath::Pi()*80/180) && fMCparticle->Phi()<=TMath::Pi()
2973                                                         if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax  ){
2974
2975                                                                 Bool_t MotherFound = FindMother(track->GetLabel());
2976                                                                 if(MotherFound){
2977                                                                         fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
2978                                                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
2979                                                                                 if(fIsHFE1)fPtMC_EMCal_All->Fill(fMCparticle->Pt());
2980                                                                         }
2981                                                                 }
2982                                                         }
2983                                                 }
2984                                         }
2985
2986                                         else if(fIsMC && track->GetLabel()>=0)//ESD
2987                                         {
2988
2989                                                 if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
2990                                                 {
2991
2992                                                         fMCtrack = fMCstack->Particle(track->GetLabel());
2993
2994                                                         Int_t pdg = fMCtrack->GetPdgCode();
2995                                                         //Phi cut && fMCtrack->Phi()>=(TMath::Pi()*80/180) && fMCtrack->Phi()<=TMath::Pi()
2996                                                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax )
2997                                                         {
2998                                                                 Bool_t MotherFound = FindMother(track->GetLabel());
2999                                                                         //MotherFound included
3000                                                                 if(MotherFound){
3001                                                                         if(fMCtrack->GetFirstMother()>0){
3002                                                                                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
3003                                                                                 if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
3004                                                                                         if(fIsHFE1)fPtMC_EMCal_All->Fill(fMCtrack->Pt());
3005                                                                                 }
3006                                                                         }
3007                                                                 }
3008                                                         }
3009                                                 }
3010                                         }
3011
3012                                         //_______________________________________________________
3013                                         //PID using EMCal
3014                                         if((fClus->E() / fP) >= fEoverPCutMin && (fClus->E() / fP) <= fEoverPCutMax)
3015                                         {
3016
3017
3018                                                         fECluster[2]->Fill(Energy);
3019                                                         fTPCNcls_pid[3]->Fill(TPCNcls, TPCNcls_pid);
3020
3021
3022                                                 if(fUseEMCal)
3023                                                 {
3024                                                         if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax ){
3025                                                                 fPtElec_Inc->Fill(fPt);
3026                                                                 //eta phi distribution for data
3027                                                                 fEtaPhi_data->Fill(track->Phi(),track->Eta());
3028                                                         }
3029
3030                                                         //Eta cut for background
3031                                                         if(fFillBackground){
3032                                                                 fEtad[2]->Fill(track->Eta());
3033
3034                                                                 //background for triggered data: trigger electron must have same cuts on shower shape  06/Jan/2014
3035                                                                 if(fUseShowerShapeCut){
3036                                                                         if(M02 >= fM02CutMin && M02<=fM02CutMax && M20>=fM20CutMin && M20<=fM20CutMax){
3037                                                                                 if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax){
3038                                                                                         Background(track, iTracks, Vtrack, kFALSE);
3039                                                                                 }
3040                                                                         }
3041                                                                 }
3042                                                                 else{
3043                                                                         if(track->Eta()>=fEtaCutMin && track->Eta()<=fEtaCutMax){
3044                                                                                 Background(track, iTracks, Vtrack, kFALSE);
3045                                                                         }
3046                                                                 }
3047
3048                                                         }
3049
3050                                                         double emctof2 = fClus->GetTOF();
3051                                                         ftimingEle2->Fill(fPt,emctof2);
3052                                                         //Correlation Analysis
3053                                                         if(fCorrelationFlag)
3054                                                         {
3055                                                                 ElectronHadronCorrelation(track, iTracks, Vtrack);
3056                                                         }
3057                                                 }
3058                                                 //_______________________________________________________
3059
3060                                                 ////////////////////////////////////////////////////////////////////
3061                                                 ///EMCal - efficiency calculations
3062
3063                                                 if(track->Charge()<0)  fCharge_n->Fill(fPt);
3064                                                 if(track->Charge()>0)  fCharge_p->Fill(fPt);
3065
3066
3067
3068                                                 if(fIsMC && fIsAOD && track->GetLabel()>=0)//AOD
3069                                                 {
3070                                                         if(track->Charge()<0)  fCharge_n->Fill(fPt);
3071                                                         if(track->Charge()>0)  fCharge_p->Fill(fPt);
3072
3073                                                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
3074                                                         if(fMCparticle->GetMother()>0) fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
3075                                                         Int_t pdg = fMCparticle->GetPdgCode();
3076
3077                                                         double proX = fMCparticle->Xv();
3078                                                         double proY = fMCparticle->Yv();
3079                                                         double proR = sqrt(pow(proX,2)+pow(proY,2));
3080
3081
3082                                                         if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax && fMCparticle->Phi()>=(TMath::Pi()*80/180) && fMCparticle->Phi()<=TMath::Pi()  ){
3083
3084                                                                 if( TMath::Abs(pdg) == 11 && fMCparticle->GetMother()>0 ){
3085                                                                         Int_t mpdg = fMCparticleMother->GetPdgCode();
3086                                                                         if(TMath::Abs(mpdg) == 221 || TMath::Abs(mpdg) == 22 || TMath::Abs(mpdg) == 111){
3087                                                                                 if(proR<7)fPtMCelectronAfterAll_nonPrimary->Fill(fMCparticle->Pt()); //numerator for the total efficiency, non Primary track
3088                                                                         }
3089                                                                 }
3090                                                                 if( TMath::Abs(pdg) == 11 && fMCparticle->IsPhysicalPrimary()) fPtMCelectronAfterAll_Primary->Fill(fMCparticle->Pt());
3091                                                         }
3092
3093
3094                                                         //
3095                                                         if(fMCparticle->IsPhysicalPrimary()){
3096
3097
3098                                                                 if(fMCparticle->Eta()>=fEtaCutMin && fMCparticle->Eta()<=fEtaCutMax ){
3099
3100                                                                         Bool_t MotherFound = FindMother(track->GetLabel());
3101                                                                         if(MotherFound){
3102
3103                                                                                 if(!fUseShowerShapeCut){
3104                                                                                         if(fIsHFE1){
3105                                                                                                 //Unfolding   pt_reco/pt_MC  in the efficiency
3106                                                                                                 fPtMCelectronAfterAll->Fill(fMCparticle->Pt());
3107                                                                                                 fPtMCelectronAfterAll_unfolding->Fill(track->Pt());
3108                                                                                                 fEtaPhi_num->Fill(fMCparticle->Phi(),fMCparticle->Eta());
3109
3110                                                                                                 //new histo to estimate how different is pt reco from pt MC
3111                                                                                                 fpt_reco_pt_MC_num->Fill(track->Pt(),fMCparticle->Pt());
3112                                                                                         }//numerator for the total efficiency AOD
3113                                                                                 }
3114                                                                                         //November 11 for efficiency of triggered data
3115                                                                                 if(fUseShowerShapeCut){
3116                                                                                         if(M02 >= fM02CutMin && M02<=fM02CutMax && M20>=fM20CutMin && M20<=fM20CutMax){
3117                                                                                                 if(fIsHFE1){
3118                                                                                                         //Unfolding   pt_reco/pt_MC  in the efficiency
3119                                                                                                         fPtMCelectronAfterAll->Fill(fMCparticle->Pt());
3120                                                                                                         fPtMCelectronAfterAll_unfolding->Fill(track->Pt());
3121                                                                                                         fEtaPhi_num->Fill(fMCparticle->Phi(),fMCparticle->Eta());
3122
3123                                                                                                         //new histo to estimate how different is pt reco from pt MC
3124                                                                                                         fpt_reco_pt_MC_num->Fill(track->Pt(),fMCparticle->Pt());
3125                                                                                                 }//numerator for the total efficiency AOD
3126                                                                                         }
3127                                                                                 }
3128
3129
3130
3131                                                                                         fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
3132                                                                                         if( TMath::Abs(pdg) == 11 && fMCparticleMother->GetPdgCode()!=22 ){
3133                                                                                                 if(fIsHFE1)fPtMC_EMCal_Selected->Fill(fMCparticle->Pt());
3134                                                                                         }
3135                                                                         }//if MotherFound
3136                                                                 }//eta cut
3137                                                         }
3138                                                 }///close AOD
3139
3140                                                 else if(fIsMC && track->GetLabel()>=0)//ESD
3141                                                 {
3142                                                         if(track->Charge()<0)  fCharge_n->Fill(fPt);
3143                                                         if(track->Charge()>0)  fCharge_p->Fill(fPt);
3144
3145                                                         fMCtrack = fMCstack->Particle(track->GetLabel());
3146                                                         if(fMCtrack->GetFirstMother()>0)
3147                                                         {
3148                                                                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
3149                                                         }
3150                                                         TParticle *particle=fMCstack->Particle(track->GetLabel());
3151
3152                                                         Int_t pdg = fMCtrack->GetPdgCode();
3153
3154
3155                                                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
3156                                                         {
3157                                                                 if( TMath::Abs(pdg) == 11 && fMCtrack->GetFirstMother()>0 )
3158                                                                 {
3159                                                                         Int_t mpdg = fMCtrackMother->GetPdgCode();
3160                                                                         if(TMath::Abs(mpdg) == 221 || TMath::Abs(mpdg) == 22 || TMath::Abs(mpdg) == 111)
3161                                                                         {
3162                                                                                 Double_t proR=particle->R();
3163                                                                                 if(proR<7)
3164                                                                                 {
3165                                                                                   fPtMCelectronAfterAll_nonPrimary->Fill(fMCtrack->Pt()); //numerator for the total efficiency, non Primary track
3166                                                                                 }
3167                                                                         }
3168                                                                 }
3169                                                                 if( TMath::Abs(pdg) == 11 && fMCstack->IsPhysicalPrimary(track->GetLabel()))
3170                                                                 {
3171                                                                         fPtMCelectronAfterAll_Primary->Fill(fMCtrack->Pt());
3172                                                                 }
3173                                                         }
3174
3175                                                         if(fMCstack->IsPhysicalPrimary(track->GetLabel()))
3176                                                         {
3177                                                                 Bool_t MotherFound = FindMother(track->GetLabel());
3178
3179                                                                 if(MotherFound)
3180                                                                 {
3181                                                                         if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
3182                                                                         {
3183                                                                                 if(!fUseShowerShapeCut){
3184                                                                                         if(fIsHFE1){
3185                                                                                                 fPtMCelectronAfterAll->Fill(fMCtrack->Pt()); //numerator for the total efficiency ESD
3186                                                                                                 fEtaPhi_num->Fill(fMCtrack->Phi(),fMCtrack->Eta());
3187                                                                                         }
3188                                                                                 }
3189                                                                                 //November 11 for efficiency of triggered data
3190                                                                                 if(fUseShowerShapeCut){
3191                                                                                         if(M02 >= fM02CutMin && M02<=fM02CutMax && M20>=fM20CutMin && M20<=fM20CutMax){
3192                                                                                                 if(fIsHFE1){
3193                                                                                                         fPtMCelectronAfterAll->Fill(fMCtrack->Pt()); //numerator for the total efficiency ESD
3194                                                                                                         fEtaPhi_num->Fill(fMCtrack->Phi(),fMCtrack->Eta());
3195                                                                                                 }
3196                                                                                         }
3197                                                                                 }
3198
3199
3200
3201                                                                         }
3202                                                                 }
3203
3204
3205
3206                                                                 // Phi cut && fMCtrack->Phi()>=(TMath::Pi()*80/180) && fMCtrack->Phi()<=TMath::Pi()
3207                                                                 if(fMCtrack->Eta()>=fEtaCutMin && fMCtrack->Eta()<=fEtaCutMax)
3208                                                                 {
3209                                                                         //included MotherFound
3210
3211                                                                         if(MotherFound)
3212                                                                         {
3213                                                                                 if(fMCtrack->GetFirstMother()>0){
3214                                                                                         fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
3215                                                                                         if( TMath::Abs(pdg) == 11 && fMCtrackMother->GetPdgCode()!=22 ){
3216
3217                                                                                                 if(fIsHFE1){fPtMC_EMCal_Selected->Fill(fMCtrack->Pt());}
3218                                                                                         }
3219                                                                                 }
3220                                                                         }
3221                                                                 }
3222                                                         }
3223                                                 }//close ESD
3224                                                 ///////////////////////////////////////////////////////////////////
3225
3226
3227                                         }
3228                                 }
3229                         }
3230                 }
3231
3232                         //______________________________________________________________
3233                         // Vertex
3234
3235                 fVtxZ[2]->Fill(fZvtx);
3236                 if(iTracks == 0)fNTracks[2]->Fill(fNOtrks);
3237                 fNTracks_pt[2]->Fill(fNOtrks, fPt);
3238                 fNTracks_eta[2]->Fill(fNOtrks, track->Eta());
3239                 fNTracks_phi[2]->Fill(fNOtrks, track->Phi());
3240                 fNClusters[2]->Fill(ClsNo);
3241                 fTPCNcls_pid[2]->Fill(TPCNcls, TPCNcls_pid);
3242
3243                         //______________________________________________________________
3244
3245                         //_______________________________________________________
3246                         //Correlation Analysis
3247                 if(!fUseEMCal)
3248                 {
3249                         fPtElec_Inc->Fill(fPt);
3250
3251                         if(fCorrelationFlag)
3252                         {
3253                                 ElectronHadronCorrelation(track, iTracks, Vtrack);
3254                         }
3255                 }
3256                         //_______________________________________________________
3257
3258                         ///________________________________________________________________________
3259         }
3260
3261                 //__________________________________________________________________
3262                 //Event Mixing Analysis
3263                 //Filling pool
3264         if(fEventMixingFlag)
3265         {
3266                 fPool = fPoolMgr->GetEventPool(fCentrality->GetCentralityPercentile("V0A"), fZvtx); // Get the buffer associated with the current centrality and z-vtx
3267
3268                 if(!fPool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality->GetCentralityPercentile("V0A"), fZvtx));
3269
3270                 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!
3271
3272
3273         }
3274
3275                 //__________________________________________________________________
3276
3277         delete fListOfmotherkink;
3278         PostData(1, fOutputList);
3279 }
3280
3281         //______________________________________________________________________
3282 void AliAnalysisTaskEMCalHFEpA::Terminate(Option_t *)
3283 {
3284                 //Draw result to the screen
3285                 //Called once at the end of the query
3286
3287         fOutputList = dynamic_cast<TList*> (GetOutputData(1));
3288
3289         if(!fOutputList)
3290         {
3291                 printf("ERROR: Output list not available\n");
3292                 return;
3293         }
3294 }
3295
3296         //______________________________________________________________________
3297 Bool_t AliAnalysisTaskEMCalHFEpA::ProcessCutStep(Int_t cutStep, AliVParticle *track)
3298 {
3299                 //Check single track cuts for a given cut step
3300                 //Note this function is called inside the UserExec function
3301         const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
3302         if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
3303         return kTRUE;
3304 }
3305
3306
3307         //______________________________________________________________________
3308
3309
3310 void AliAnalysisTaskEMCalHFEpA::Background(AliVTrack *track, Int_t trackIndex, AliVParticle *vtrack, Bool_t IsTPConly)
3311 {
3312         ///_________________________________________________________________
3313         ///MC analysis
3314         //Bool_t IsMCefix=kFALSE; //to make correction on efix, use kTRUE (do not change the efficiency, so I will keep the correction only for d3)
3315
3316                 if(fIsMC)
3317                 {
3318                         if(track->GetLabel() < 0)
3319                         {
3320                                 AliWarning(Form("The track %d does not have a valid MC label",trackIndex));
3321                                 return;
3322                         }
3323
3324                         if(fIsAOD)
3325                         {
3326                                 fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
3327
3328                                 if(fMCparticle->GetMother()<0) return;
3329
3330                                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
3331                                 if(fMCparticleMother->GetMother()>0)fMCparticleGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleMother->GetMother());
3332
3333                                 if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
3334                                 {
3335                                                 //Is Background
3336                                         if(!IsTPConly)fPtBackgroundBeforeReco->Fill(track->Pt());
3337                                         if(IsTPConly)fPtBackgroundBeforeReco2->Fill(track->Pt());
3338
3339
3340                                                 //October 08th weighted histos
3341                                         if(TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221 ){
3342
3343                                                 Double_t mPt=fMCparticleMother->Pt();
3344                                                 Double_t mweight=1;
3345
3346
3347                                                 //________________________________________________________________
3348                                                 //correction for d3 based on data //from Jan
3349
3350                                                         if(TMath::Abs(fMCparticleMother->GetPdgCode())==111){
3351                                                                 Double_t x=mPt;
3352
3353                                                                 mweight=CalculateWeight(111, x);
3354
3355                                                         }
3356
3357                                                 //________________________________________________________________
3358
3359                                                 //Histo pT mother versus pT electron
3360                                                 fpT_m_electron->Fill(mPt, track->Pt());
3361
3362                                                 if(!IsTPConly)fPtBackgroundBeforeReco_weight->Fill(track->Pt(), 1./mweight);
3363                                                 if(IsTPConly)fPtBackgroundBeforeReco2_weight->Fill(track->Pt(), 1./mweight);
3364                                         }
3365                                         else if(fMCparticleMother->GetMother()>0 && (TMath::Abs(fMCparticleGMother->GetPdgCode())==111 || TMath::Abs(fMCparticleGMother->GetPdgCode())==221 )){
3366
3367                                                 Double_t gmPt=fMCparticleGMother->Pt();
3368                                                 Double_t gmweight=1;
3369
3370                                                 //________________________________________________________________
3371                                                 //correction for d3 based on data
3372
3373                                                         if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111){
3374                                                                 Double_t x=gmPt;
3375                                                                 gmweight=CalculateWeight(111, x);
3376                                                         }
3377
3378
3379                                                 //________________________________________________________________
3380                                                 //Histo pT gmother versus pT electron
3381
3382                                                 fpT_gm_electron->Fill(gmPt, track->Pt());
3383
3384                                                 if(!IsTPConly)fPtBackgroundBeforeReco_weight->Fill(track->Pt(), 1./gmweight);
3385                                                 if(IsTPConly)fPtBackgroundBeforeReco2_weight->Fill(track->Pt(), 1./gmweight);
3386                                         }
3387                                         else{
3388                                                 if(!IsTPConly)fPtBackgroundBeforeReco_weight->Fill(track->Pt());
3389                                                 if(IsTPConly)fPtBackgroundBeforeReco2_weight->Fill(track->Pt());
3390                                         }
3391                                 }//particle kind
3392                         }//IsAOD
3393                          //ESD
3394                         else
3395                         {
3396                                 fMCtrack = fMCstack->Particle(track->GetLabel());
3397
3398                                 if(fMCtrack->GetFirstMother()<0) return;
3399
3400                                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
3401
3402                                 if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
3403                                 {
3404                                                 //Is Background
3405                                         if(!IsTPConly)fPtBackgroundBeforeReco->Fill(track->Pt());
3406                                         if(IsTPConly)fPtBackgroundBeforeReco2->Fill(track->Pt());
3407                                 }
3408                         }
3409                 }//IsMC
3410
3411                 ///_________________________________________________________________
3412
3413                 //________________________________________________
3414                 //Associated particle cut
3415                 fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
3416                 fPartnerCuts->SetRequireITSRefit(kTRUE);
3417                 fPartnerCuts->SetRequireTPCRefit(kTRUE);
3418                 fPartnerCuts->SetEtaRange(-0.9,0.9);
3419                 fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
3420                 fPartnerCuts->SetMinNClustersTPC(80);
3421                 fPartnerCuts->SetPtRange(0,1e10);
3422                 //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
3423                 //fPartnerCuts->SetMaxDCAToVertexXY(1);
3424                 //fPartnerCuts->SetMaxDCAToVertexZ(3);
3425                 //_________________________________________________
3426
3427                 ///#################################################################
3428                 //Non-HFE reconstruction
3429                 fNonHFE = new AliSelectNonHFE();
3430                 fNonHFE->SetAODanalysis(fIsAOD);
3431                 if(fMassCutFlag) fNonHFE->SetInvariantMassCut(fMassCut);
3432                 if(fAngleCutFlag) fNonHFE->SetOpeningAngleCut(fAngleCut);
3433                 if(fChi2CutFlag) fNonHFE->SetChi2OverNDFCut(fChi2Cut);
3434                 if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
3435                 fNonHFE->SetAlgorithm("DCA"); //KF
3436                 fNonHFE->SetPIDresponse(fPidResponse);
3437                 fNonHFE->SetTrackCuts(-3.5,3.5,fPartnerCuts);
3438                 fNonHFE->SetAdditionalCuts(fPtMinAsso,fTpcNclsAsso);
3439
3440                 if(!IsTPConly){
3441                         fNonHFE->SetHistAngleBack(fOpAngleBack);
3442                         fNonHFE->SetHistAngle(fOpAngle);
3443                         fNonHFE->SetHistDCABack(fDCABack);
3444                         fNonHFE->SetHistDCA(fDCA);
3445                         fNonHFE->SetHistMassBack(fInvMassBack);
3446                         fNonHFE->SetHistMass(fInvMass);
3447                 }
3448                 if(IsTPConly){
3449                         fNonHFE->SetHistAngleBack(fOpAngleBack2);
3450                         fNonHFE->SetHistAngle(fOpAngle2);
3451                         fNonHFE->SetHistDCABack(fDCABack2);
3452                         fNonHFE->SetHistDCA(fDCA2);
3453                         fNonHFE->SetHistMassBack(fInvMassBack2);
3454                         fNonHFE->SetHistMass(fInvMass2);
3455                 }
3456
3457                 fNonHFE->FindNonHFE(trackIndex,vtrack,fVevent);
3458
3459                         //index of track selected as partner
3460                 Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
3461
3462                         //Electron Information
3463                 Double_t fPhiE = -999;
3464                 Double_t fEtaE = -999;
3465                 Double_t fPtE = -999;
3466                 fPhiE = track->Phi();
3467                 fEtaE = track->Eta();
3468                 fPtE = track->Pt();
3469
3470                         ///_________________________________________________________________
3471                         ///MC analysis
3472                 if(fIsMC)
3473                 {
3474                         if(fIsAOD)
3475                         {
3476                                 if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
3477                                 {
3478
3479                                         Double_t weight=1;
3480
3481                                         if(!IsTPConly){
3482                                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
3483                                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
3484
3485
3486
3487
3488                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221){
3489                                                         Double_t mPt=fMCparticleMother->Pt();
3490                                                                 Double_t mweight1=1;
3491                                                                 Double_t mweight2=1;
3492                                                                 //Double_t weight=1;
3493
3494
3495                                                          //----------------------------------------------------------------------------
3496                                                          //correction based on data only for pi0
3497                                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111){
3498                                                                         Double_t x=mPt;
3499                                                                         weight=CalculateWeight(111, x);
3500                                                                 }
3501
3502
3503                                                                 //----------------------------------------------------------------------------
3504
3505                                                                 //check this
3506                                                         if(fNonHFE->IsULS()) mweight1=(fNonHFE->GetNULS())/weight;
3507                                                         if(fNonHFE->IsLS())  mweight2=(fNonHFE->GetNLS())/weight;
3508
3509                                                                 //fill histos
3510                                                         if(fNonHFE->IsULS())fPtElec_ULS_weight->Fill(fPtE, mweight1);
3511                                                         if(fNonHFE->IsLS())fPtElec_LS_weight->Fill(fPtE, mweight2);
3512                                                 }
3513                                                 else if(fMCparticleMother->GetMother()>0 && (TMath::Abs(fMCparticleGMother->GetPdgCode())==111 || TMath::Abs(fMCparticleGMother->GetPdgCode())==221 )){
3514                                                         Double_t gmPt=fMCparticleGMother->Pt();
3515                                                                 Double_t gmweight1=1;
3516                                                                 Double_t gmweight2=1;
3517                                                                 //Double_t weight=1;
3518
3519                                                         //----------------------------------------------------------------------------
3520
3521                                                         //----------------------------------------------------------------------------
3522
3523                                                         //correction based on data only for pi0
3524
3525                                                                 if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111){
3526                                                                         Double_t x=gmPt;
3527                                                                         weight=CalculateWeight(111, x);
3528                                                                 }
3529
3530
3531
3532
3533
3534                                                                 //check this
3535                                                         if(fNonHFE->IsULS()) gmweight1=(fNonHFE->GetNULS())/weight;
3536                                                         if(fNonHFE->IsLS())  gmweight2=(fNonHFE->GetNLS())/weight;
3537
3538                                                                 //fill histos
3539                                                         if(fNonHFE->IsULS())fPtElec_ULS_weight->Fill(fPtE, gmweight1);
3540                                                         if(fNonHFE->IsLS())fPtElec_LS_weight->Fill(fPtE, gmweight2);
3541                                                 }
3542                                                 else{
3543                                                         if(fNonHFE->IsULS()) fPtElec_ULS_weight->Fill(fPtE,fNonHFE->GetNULS());
3544                                                         if(fNonHFE->IsLS()) fPtElec_LS_weight->Fill(fPtE,fNonHFE->GetNLS());
3545                                                 }
3546
3547
3548                                         }//!IsTPConly
3549
3550                                         if(IsTPConly){
3551                                                 if(fNonHFE->IsULS()) fPtElec_ULS2->Fill(fPtE,fNonHFE->GetNULS());
3552                                                 if(fNonHFE->IsLS()) fPtElec_LS2->Fill(fPtE,fNonHFE->GetNLS());
3553
3554                                                         //new 08 October        //weighted histograms
3555                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221){
3556                                                         Double_t mPt=fMCparticleMother->Pt();
3557
3558                                                                 Double_t mweight1=1;
3559                                                                 Double_t mweight2=1;
3560                                                                 //Double_t weight=1;
3561
3562
3563                                                          //----------------------------------------------------------------------------
3564
3565                                                         //correction based on data only for pi0 for d3
3566
3567                                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111){
3568                                                                         Double_t x=mPt;
3569                                                                         weight=CalculateWeight(111, x);
3570
3571                                                                 }
3572
3573
3574                                                                 //check this
3575                                                         if(fNonHFE->IsULS()) mweight1=(fNonHFE->GetNULS())/weight;
3576                                                         if(fNonHFE->IsLS())  mweight2=(fNonHFE->GetNLS())/weight;
3577
3578                                                                 //fill histos
3579                                                         if(fNonHFE->IsULS())fPtElec_ULS2_weight->Fill(fPtE, mweight1);
3580                                                         if(fNonHFE->IsLS())fPtElec_LS2_weight->Fill(fPtE, mweight2);
3581                                                 }
3582                                                 else if(fMCparticleMother->GetMother()>0 && (TMath::Abs(fMCparticleGMother->GetPdgCode())==111 || TMath::Abs(fMCparticleGMother->GetPdgCode())==221 )){
3583                                                         Double_t gmPt=fMCparticleGMother->Pt();
3584                                                         Double_t gmweight1=1;
3585                                                         Double_t gmweight2=1;
3586                                                                 //Double_t weight=1;
3587
3588
3589                                                         //----------------------------------------------------------------------------
3590                                                          //correction based on data only for pi0
3591
3592                                                                 if(TMath::Abs(fMCparticleGMother->GetPdgCode())==111){
3593                                                                         Double_t x=gmPt;
3594
3595                                                                         weight=CalculateWeight(111, x);
3596
3597                                                                 }
3598
3599                                                         //----------------------------------------------------------------------------
3600
3601
3602
3603                                                                 //check this
3604                                                         if(fNonHFE->IsULS()) gmweight1=(fNonHFE->GetNULS())/weight;
3605                                                         if(fNonHFE->IsLS())  gmweight2=(fNonHFE->GetNLS())/weight;
3606
3607                                                                 //fill histos
3608                                                         if(fNonHFE->IsULS())fPtElec_ULS2_weight->Fill(fPtE, gmweight1);
3609                                                         if(fNonHFE->IsLS())fPtElec_LS2_weight->Fill(fPtE, gmweight2);
3610                                                 }
3611                                                 else{
3612                                                         if(fNonHFE->IsULS()) fPtElec_ULS2_weight->Fill(fPtE,fNonHFE->GetNULS());
3613                                                         if(fNonHFE->IsLS()) fPtElec_LS2_weight->Fill(fPtE,fNonHFE->GetNLS());
3614                                                 }
3615
3616
3617                                                 //----------------------------------------------------------------------------
3618                                                 //to check other way to calculate efficiency
3619                                                 //ULS with no weight from ULS-LS original
3620                                                 //we have to know if track2 comes from same mother!!!
3621                                                 //----------------------------------------------------------------------------
3622                                                 if(fNonHFE->IsULS()){
3623
3624                                                         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
3625                                                         {
3626
3627                                                                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
3628                                                                 if (!Vtrack2)
3629                                                                 {
3630                                                                         printf("ERROR: Could not receive track %d\n", iTracks);
3631                                                                         continue;
3632                                                                 }
3633                                                                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
3634                                                                 if(track2->GetLabel()<0) continue;
3635                                                                 fMCparticle2 = (AliAODMCParticle*) fMCarray->At(track2->GetLabel());
3636                                                                 if(fMCparticle2->GetMother()<0) continue;
3637
3638                                                                 for(Int_t i = 0; i < fNonHFE->GetNULS(); i++)
3639                                                                 {
3640                                                                         if(fUlsPartner[i]==iTracks){
3641                                                                                         //only fill if it has same mother
3642                                                                                         //with weight to take into account the number of partners
3643                                                                                 if(fMCparticle2->GetMother()==fMCparticle->GetMother()) fPtElec_ULS_MC->Fill(fPtE, fNonHFE->GetNULS());
3644
3645                                                                                 //-----------------------------------------------------------------------------------------------------------
3646                                                                                 //weight for mother
3647                                                                                 //Double_t weight2=1;
3648                                                                                 Double_t mPt=fMCparticleMother->Pt();
3649
3650
3651                                                                                 if(TMath::Abs(fMCparticleMother->GetPdgCode())==111){
3652                                                                                         Double_t x=mPt;
3653                                                                                         weight=CalculateWeight(111, x);
3654
3655
3656                                                                                 }
3657
3658                                                                                 //weight for grandmother
3659                                                                                 Double_t gmPt=fMCparticleGMother->Pt();
3660                                                                                 if(TMath::Abs((fMCparticleMother->GetMother()>0) && ((fMCparticleGMother->GetPdgCode())==111))){
3661                                                                                         Double_t x=gmPt;
3662                                                                                         weight=CalculateWeight(111, x);
3663
3664                                                                                 }
3665
3666
3667                                                                                 if(fMCparticle2->GetMother()==fMCparticle->GetMother()) fPtElec_ULS_MC_weight->Fill(fPtE, (fNonHFE->GetNULS())*1./weight);
3668
3669                                                                                 //-----------------------------------------------------------------------------------------------------------
3670                                                                                 //end of weight
3671
3672                                                                         }//partner found same as track
3673                                                                 }//loop in all partner
3674
3675                                                         }//track
3676                                                 }//is ULS
3677                                                 //----------------------------------------------------------------------------
3678                                                 //end of part to check other way to calculate efficiency
3679                                                 //----------------------------------------------------------------------------
3680
3681                                         }//IsTPConly
3682
3683                                 }//particle kind
3684
3685                                 if(IsTPConly){
3686                                                 //ULS-LS with no pid AOD
3687                                         if(fNonHFE->IsULS()) fPtElec_ULS_NoPid->Fill(fPtE,fNonHFE->GetNULS());
3688                                         if(fNonHFE->IsLS()) fPtElec_LS_NoPid->Fill(fPtE,fNonHFE->GetNLS());
3689                                 }
3690
3691                         }//close IsAOD
3692                          //It is ESD
3693                         else
3694                         {
3695                                 if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
3696                                 {
3697                                         if(!IsTPConly){
3698                                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
3699                                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
3700                                         }
3701
3702                                         if(IsTPConly){
3703                                                 if(fNonHFE->IsULS()) fPtElec_ULS2->Fill(fPtE,fNonHFE->GetNULS());
3704                                                 if(fNonHFE->IsLS()) fPtElec_LS2->Fill(fPtE,fNonHFE->GetNLS());
3705                                         }
3706                                 }
3707
3708
3709                         }
3710                 }//close IsMC
3711                  ///_________________________________________________________________
3712                  //not MC
3713                 else
3714                 {
3715                         if(!IsTPConly){
3716                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
3717                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
3718                         }
3719
3720                         if(IsTPConly){
3721                                 if(fNonHFE->IsULS()) fPtElec_ULS2->Fill(fPtE,fNonHFE->GetNULS());
3722                                 if(fNonHFE->IsLS()) fPtElec_LS2->Fill(fPtE,fNonHFE->GetNLS());
3723                         }
3724                 }
3725
3726
3727
3728 }//end of Background function
3729
3730         //______________________________________________________________________
3731 void AliAnalysisTaskEMCalHFEpA::ElectronHadronCorrelation(AliVTrack *track, Int_t trackIndex, AliVParticle *vtrack)
3732 {
3733
3734         ///_________________________________________________________________
3735         ///MC analysis
3736         if(fIsMC)
3737         {
3738                 if(track->GetLabel() < 0)
3739         {
3740                         AliWarning(Form("The track %d does not have a valid MC label",trackIndex));
3741                         return;
3742         }
3743
3744                 if(fIsAOD)
3745                 {
3746                         fMCparticle = (AliAODMCParticle*) fMCarray->At(track->GetLabel());
3747
3748                         if(fMCparticle->GetMother()<0) return;
3749
3750                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
3751
3752                 if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
3753                 {
3754                                         //Is Background
3755                                 fPtBackgroundBeforeReco->Fill(track->Pt());
3756                         }
3757                 }
3758                 else
3759                 {
3760                 fMCtrack = fMCstack->Particle(track->GetLabel());
3761
3762                 if(fMCtrack->GetFirstMother()<0) return;
3763
3764                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
3765
3766                 if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
3767                 {
3768                                         //Is Background
3769                                 fPtBackgroundBeforeReco->Fill(track->Pt());
3770                         }
3771                 }
3772         }
3773                 ///_________________________________________________________________
3774
3775                 //________________________________________________
3776                 //Associated particle cut
3777         fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
3778     fPartnerCuts->SetRequireITSRefit(kTRUE);
3779     fPartnerCuts->SetRequireTPCRefit(kTRUE);
3780     fPartnerCuts->SetEtaRange(-0.9,0.9);
3781     fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
3782     fPartnerCuts->SetMinNClustersTPC(80);
3783     fPartnerCuts->SetPtRange(0.3,1e10);
3784                 //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
3785                 //fPartnerCuts->SetMaxDCAToVertexXY(1);
3786                 //fPartnerCuts->SetMaxDCAToVertexZ(3);
3787                 //_________________________________________________
3788
3789                 ///#################################################################
3790                 //Non-HFE reconstruction
3791         fNonHFE = new AliSelectNonHFE();
3792         fNonHFE->SetAODanalysis(fIsAOD);
3793         if(fMassCutFlag) fNonHFE->SetInvariantMassCut(fMassCut);
3794         if(fAngleCutFlag) fNonHFE->SetOpeningAngleCut(fAngleCut);
3795         if(fChi2CutFlag) fNonHFE->SetChi2OverNDFCut(fChi2Cut);
3796         if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
3797         fNonHFE->SetAlgorithm("DCA"); //KF
3798         fNonHFE->SetPIDresponse(fPidResponse);
3799         fNonHFE->SetTrackCuts(-3.5,3.5,fPartnerCuts);
3800         fNonHFE->SetAdditionalCuts(fPtMinAsso,fTpcNclsAsso);
3801
3802
3803         fNonHFE->SetHistAngleBack(fOpAngleBack);
3804         fNonHFE->SetHistAngle(fOpAngle);
3805         fNonHFE->SetHistDCABack(fDCABack);
3806         fNonHFE->SetHistDCA(fDCA);
3807         fNonHFE->SetHistMassBack(fInvMassBack);
3808         fNonHFE->SetHistMass(fInvMass);
3809
3810         fNonHFE->FindNonHFE(trackIndex,vtrack,fVevent);
3811
3812         Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
3813         Int_t *fLsPartner = fNonHFE->GetPartnersLS();
3814         Bool_t fUlsIsPartner = kFALSE;
3815         Bool_t fLsIsPartner = kFALSE;
3816                 ///#################################################################
3817
3818
3819                 //Electron Information
3820         Double_t fPhiE = -999;
3821         Double_t fEtaE = -999;
3822         Double_t fPhiH = -999;
3823         Double_t fEtaH = -999;
3824         Double_t fDphi = -999;
3825         Double_t fDeta = -999;
3826         Double_t fPtE = -999;
3827         Double_t fPtH = -999;
3828
3829         Double_t pi = TMath::Pi();
3830
3831         fPhiE = track->Phi();
3832         fEtaE = track->Eta();
3833         fPtE = track->Pt();
3834
3835
3836                 ///_________________________________________________________________
3837                 ///MC analysis
3838         if(fIsMC)
3839         {
3840                 if(fIsAOD)
3841                 {
3842                         if(TMath::Abs(fMCparticle->GetPdgCode())==11 && (TMath::Abs(fMCparticleMother->GetPdgCode())==22 || TMath::Abs(fMCparticleMother->GetPdgCode())==111 || TMath::Abs(fMCparticleMother->GetPdgCode())==221))
3843                 {
3844                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
3845                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
3846
3847
3848                         }
3849
3850                         if(fNonHFE->IsULS()) fPtElec_ULS_NoPid->Fill(fPtE,fNonHFE->GetNULS());
3851                         if(fNonHFE->IsLS()) fPtElec_LS_NoPid->Fill(fPtE,fNonHFE->GetNLS());
3852                 }
3853                 else
3854                 {
3855                         if(TMath::Abs(fMCtrack->GetPdgCode())==11 && (TMath::Abs(fMCtrackMother->GetPdgCode())==22 || TMath::Abs(fMCtrackMother->GetPdgCode())==111 || TMath::Abs(fMCtrackMother->GetPdgCode())==221))
3856                         {
3857                                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
3858                                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
3859
3860                         }
3861
3862
3863                 }
3864         }
3865                 ///_________________________________________________________________
3866         else
3867         {
3868                 if(fNonHFE->IsULS()) fPtElec_ULS->Fill(fPtE,fNonHFE->GetNULS());
3869                 if(fNonHFE->IsLS()) fPtElec_LS->Fill(fPtE,fNonHFE->GetNLS());
3870         }
3871
3872
3873
3874
3875                 //__________________________________________________________________
3876                 //Event Mixing Analysis - Hadron Loop
3877                 //Retrieve
3878         if(fEventMixingFlag)
3879         {
3880                 fPool = fPoolMgr->GetEventPool(fCentrality->GetCentralityPercentile("V0A"), fZvtx); // Get the buffer associated with the current centrality and z-vtx
3881
3882                 if(!fPool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f",fCentrality->GetCentralityPercentile("V0A"), fZvtx));
3883
3884                 if(fPool->GetCurrentNEvents() >= 5) // start mixing when 5 events are in the buffer
3885                 {
3886                         fPoolNevents->Fill(fPool->GetCurrentNEvents());
3887
3888                         for (Int_t jMix = 0; jMix < fPool->GetCurrentNEvents(); jMix++)  // mix with each event in the buffer
3889                         {
3890                                 TObjArray* bgTracks = fPool->GetEvent(jMix);
3891
3892                                 for (Int_t kMix = 0; kMix < bgTracks->GetEntriesFast(); kMix++)  // mix with each track in the event
3893                                 {
3894                                         const AliEHCParticle* MixedTrack(dynamic_cast<AliEHCParticle*>(bgTracks->At(kMix)));
3895                                         if (NULL == MixedTrack) continue;
3896
3897                                         fPhiH = MixedTrack->Phi();
3898                                         fEtaH = MixedTrack->Eta();
3899                                         fPtH = MixedTrack->Pt();
3900
3901                                         if(fPtH<fAssHadronPtMin || fPtH>fAssHadronPtMax) continue;
3902
3903                                         fDphi = fPhiE - fPhiH;
3904
3905                                         if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
3906                                         if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
3907
3908                                         fDeta = fEtaE - fEtaH;
3909
3910                                         Double_t fPtBin[7] = {1,2,4,6,8,10,15};
3911
3912                                         for(Int_t i = 0; i < 6; i++)
3913                                         {
3914                                             if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
3915                                             {
3916                                                         fCEtaPhi_Inc_EM[i]->Fill(fDphi,fDeta);
3917
3918                                                         if(fNonHFE->IsULS()) fCEtaPhi_ULS_EM[i]->Fill(fDphi,fDeta);
3919                                                         if(fNonHFE->IsLS()) fCEtaPhi_LS_EM[i]->Fill(fDphi,fDeta);
3920
3921                                                         if(fNonHFE->IsULS()) fCEtaPhi_ULS_Weight_EM[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
3922                                                         if(fNonHFE->IsLS()) fCEtaPhi_LS_Weight_EM[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
3923                                             }
3924                                         }
3925
3926                                                 // TODO your code: do event mixing with current event and bgTracks
3927                                                 // note that usually the content filled now is weighted by 1 / pool->GetCurrentNEvents()
3928                                 }
3929                         }
3930                 }
3931         }
3932                 //__________________________________________________________________
3933
3934                 //__________________________________________________________________
3935                 //Same Event Analysis - Hadron Loop
3936         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
3937         {
3938                 if(trackIndex==iTracks) continue;
3939
3940                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
3941                 if (!Vtrack2)
3942                 {
3943                         printf("ERROR: Could not receive track %d\n", iTracks);
3944                         continue;
3945                 }
3946
3947                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
3948
3949                 if(track2->Eta()<fEtaCutMin || track2->Eta()>fEtaCutMax ) continue;
3950
3951                 if(fIsAOD)
3952                 {
3953                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
3954                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
3955                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
3956                         if(atrack2->GetTPCNcls() < 80) continue;
3957                         if(fAssocWithSPD && ((!(atrack2->HasPointOnITSLayer(0))) && (!(atrack2->HasPointOnITSLayer(1))))) continue;
3958                 }
3959                 else
3960                 {
3961                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
3962                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
3963                 }
3964
3965                 fPhiH = track2->Phi();
3966                 fEtaH = track2->Eta();
3967                 fPtH = track2->Pt();
3968
3969                 if(fPtH<fAssHadronPtMin || fPtH>fAssHadronPtMax) continue;
3970
3971                 fDphi = fPhiE - fPhiH;
3972
3973                 if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
3974                 if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
3975
3976                 fDeta = fEtaE - fEtaH;
3977
3978                 Double_t fPtBin[7] = {1,2,4,6,8,10,15};
3979
3980                 //______________________________________________________________
3981                 //Check if this track is a Non-HFE partner
3982                 for(Int_t i = 0; i < fNonHFE->GetNULS(); i++)
3983                 {
3984                         if(fUlsPartner[i]==iTracks) fUlsIsPartner=kTRUE;
3985                 }
3986                 for(Int_t i = 0; i < fNonHFE->GetNLS(); i++)
3987                 {
3988                         if(fLsPartner[i]==iTracks) fLsIsPartner=kTRUE;
3989                 }
3990                 //______________________________________________________________
3991
3992                 for(Int_t i = 0; i < 6; i++)
3993                 {
3994                     if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
3995                     {
3996                                 fCEtaPhi_Inc[i]->Fill(fDphi,fDeta);
3997
3998                                 if(fNonHFE->IsULS()) fCEtaPhi_ULS[i]->Fill(fDphi,fDeta);
3999                                 if(fNonHFE->IsLS()) fCEtaPhi_LS[i]->Fill(fDphi,fDeta);
4000                                 if(fNonHFE->IsULS() && !fUlsIsPartner) fCEtaPhi_ULS_NoP[i]->Fill(fDphi,fDeta);
4001                                 if(fNonHFE->IsLS() && !fLsIsPartner) fCEtaPhi_LS_NoP[i]->Fill(fDphi,fDeta);
4002
4003                                 if(fNonHFE->IsULS()) fCEtaPhi_ULS_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
4004                                 if(fNonHFE->IsLS()) fCEtaPhi_LS_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
4005                                 if(fNonHFE->IsULS() && !fUlsIsPartner) fCEtaPhi_ULS_NoP_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNULS());
4006                                 if(fNonHFE->IsLS() && !fLsIsPartner) fCEtaPhi_LS_NoP_Weight[i]->Fill(fDphi,fDeta,fNonHFE->GetNLS());
4007                     }
4008                 }
4009         }
4010 }
4011
4012         //____________________________________________________________________________________________________________
4013         //Create a TObjArray with selected hadrons, for the mixed event analysis
4014 TObjArray* AliAnalysisTaskEMCalHFEpA::SelectedHadrons()
4015 {
4016         fTracksClone = new TObjArray;
4017         fTracksClone->SetOwner(kTRUE);
4018
4019         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
4020         {
4021                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
4022                 if (!Vtrack2)
4023                 {
4024                         printf("ERROR: Could not receive track %d\n", iTracks);
4025                         continue;
4026                 }
4027
4028                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
4029
4030                 if(track2->Eta()<fEtaCutMin || track2->Eta()>fEtaCutMax ) continue;
4031
4032                 if(fIsAOD)
4033                 {
4034                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
4035                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
4036                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
4037                         if(atrack2->GetTPCNcls() < 80) continue;
4038                         if(fAssocWithSPD && ((!(atrack2->HasPointOnITSLayer(0))) && (!(atrack2->HasPointOnITSLayer(1))))) continue;
4039                 }
4040                 else
4041                 {
4042                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
4043                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
4044                 }
4045
4046                 fTracksClone->Add(new AliEHCParticle(track2->Eta(), track2->Phi(), track2->Pt()));
4047         }
4048         return fTracksClone;
4049 }
4050         //____________________________________________________________________________________________________________
4051
4052         //______________________________________________________________________
4053 void AliAnalysisTaskEMCalHFEpA::DiHadronCorrelation(AliVTrack *track, Int_t trackIndex)
4054 {
4055                 //________________________________________________
4056                 //Associated particle cut
4057         fPartnerCuts->SetAcceptKinkDaughters(kFALSE);
4058     fPartnerCuts->SetRequireITSRefit(kTRUE);
4059     fPartnerCuts->SetRequireTPCRefit(kTRUE);
4060     fPartnerCuts->SetEtaRange(-0.9,0.9);
4061     fPartnerCuts->SetMaxChi2PerClusterTPC(4.0);
4062     fPartnerCuts->SetMinNClustersTPC(80);
4063     fPartnerCuts->SetPtRange(0.3,1e10);
4064                 //fPartnerCuts->SetRequireSigmaToVertex(kTRUE);
4065                 //fPartnerCuts->SetMaxDCAToVertexXY(1);
4066                 //fPartnerCuts->SetMaxDCAToVertexZ(3);
4067                 //_________________________________________________
4068
4069                 //Electron Information
4070         Double_t fPhiE = -999;
4071         Double_t fEtaE = -999;
4072         Double_t fPhiH = -999;
4073         Double_t fEtaH = -999;
4074         Double_t fDphi = -999;
4075         Double_t fDeta = -999;
4076         Double_t fPtE = -999;
4077         Double_t fPtH = -999;
4078
4079         Double_t pi = TMath::Pi();
4080
4081         fPhiE = track->Phi();
4082         fEtaE = track->Eta();
4083         fPtE = track->Pt();
4084
4085                 //__________________________________________________________________
4086                 //Same Event Analysis - Hadron Loop
4087         for(Int_t iTracks = 0; iTracks < fVevent->GetNumberOfTracks(); iTracks++)
4088         {
4089                 if(trackIndex==iTracks) continue;
4090
4091                 AliVParticle* Vtrack2 = fVevent->GetTrack(iTracks);
4092                 if (!Vtrack2)
4093                 {
4094                         printf("ERROR: Could not receive track %d\n", iTracks);
4095                         continue;
4096                 }
4097
4098                 AliVTrack *track2 = dynamic_cast<AliVTrack*>(Vtrack2);
4099                 if(track2->Eta()<fEtaCutMin || track2->Eta()>fEtaCutMax ) continue;
4100
4101                 if(fIsAOD)
4102                 {
4103                         AliAODTrack *atrack2 = dynamic_cast<AliAODTrack*>(Vtrack2);
4104                         if(!atrack2->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
4105                         if((!(atrack2->GetStatus()&AliESDtrack::kITSrefit)|| (!(atrack2->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
4106                         if(atrack2->GetTPCNcls() < 80) continue;
4107                         if(fAssocWithSPD && ((!(atrack2->HasPointOnITSLayer(0))) && (!(atrack2->HasPointOnITSLayer(1))))) continue;
4108                 }
4109                 else
4110                 {
4111                         AliESDtrack *etrack2 = dynamic_cast<AliESDtrack*>(Vtrack2);
4112                         if(!fPartnerCuts->AcceptTrack(etrack2)) continue;
4113                 }
4114
4115                 fPhiH = track2->Phi();
4116                 fEtaH = track2->Eta();
4117                 fPtH = track2->Pt();
4118
4119                 if(fPtH<fAssHadronPtMin || fPtH>fAssHadronPtMax) continue;
4120
4121                 fDphi = fPhiE - fPhiH;
4122
4123                 if (fDphi > 3*pi/2) fDphi = fDphi - 2*pi;
4124                 if (fDphi < -pi/2)  fDphi = fDphi + 2*pi;
4125
4126                 fDeta = fEtaE - fEtaH;
4127
4128                 Double_t fPtBin[7] = {1,2,4,6,8,10,15};
4129
4130                 for(Int_t i = 0; i < 6; i++)
4131                 {
4132                     if(fPtE>=fPtBin[i] && fPtE<fPtBin[i+1])
4133                     {
4134                                 fCEtaPhi_Inc_DiHadron[i]->Fill(fDphi,fDeta);
4135                     }
4136                 }
4137         }
4138 }
4139         //____________________________________________________________________________________________________________
4140
4141         //______________________________________________________________________
4142 Bool_t AliAnalysisTaskEMCalHFEpA::FindMother(Int_t mcIndex)
4143 {
4144         fIsHFE1 = kFALSE;
4145         fIsHFE2 = kFALSE;
4146         fIsNonHFE = kFALSE;
4147         fIsFromD = kFALSE;
4148         fIsFromB = kFALSE;
4149         fIsFromPi0 = kFALSE;
4150         fIsFromEta = kFALSE;
4151         fIsFromGamma = kFALSE;
4152
4153         if(mcIndex < 0 || !fIsMC)
4154         {
4155                 return kFALSE;
4156         }
4157
4158         Int_t pdg = -99999;
4159         Int_t mpdg = -99999;
4160         Int_t gmpdg = -99999;
4161         Int_t ggmpdg = -99999;
4162         Int_t gggmpdg = -99999;
4163
4164         if(fIsAOD)
4165         {
4166                 fMCparticle = (AliAODMCParticle*) fMCarray->At(mcIndex);
4167
4168                 pdg = TMath::Abs(fMCparticle->GetPdgCode());
4169
4170
4171                 if(pdg!=11)
4172                 {
4173                         fIsHFE1 = kFALSE;
4174                         fIsHFE2 = kFALSE;
4175                         fIsNonHFE = kFALSE;
4176                         fIsFromD = kFALSE;
4177                         fIsFromB = kFALSE;
4178                         fIsFromPi0 = kFALSE;
4179                         fIsFromEta = kFALSE;
4180                         fIsFromGamma = kFALSE;
4181                         return kFALSE;
4182                 }
4183
4184                 if(fMCparticle->GetMother()<0)
4185                 {
4186                         fIsHFE1 = kFALSE;
4187                         fIsHFE2 = kFALSE;
4188                         fIsNonHFE = kFALSE;
4189                         fIsFromD = kFALSE;
4190                         fIsFromB = kFALSE;
4191                         fIsFromPi0 = kFALSE;
4192                         fIsFromEta = kFALSE;
4193                         fIsFromGamma = kFALSE;
4194                         return kFALSE;
4195                 }
4196
4197                 fMCparticleMother = (AliAODMCParticle*) fMCarray->At(fMCparticle->GetMother());
4198                 mpdg = TMath::Abs(fMCparticleMother->GetPdgCode());
4199
4200                 if(fMCparticleMother->GetMother()<0)
4201                 {
4202                         gmpdg = 0;
4203                         ggmpdg = 0;
4204                         gggmpdg = 0;
4205                 }
4206                 else
4207                 {
4208                         fMCparticleGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleMother->GetMother());
4209                         gmpdg = TMath::Abs(fMCparticleGMother->GetPdgCode());
4210                         if(fMCparticleGMother->GetMother()<0)
4211                         {
4212                                 ggmpdg = 0;
4213                                 gggmpdg = 0;
4214                         }
4215                         else
4216                         {
4217                                 fMCparticleGGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleGMother->GetMother());
4218                                 ggmpdg = TMath::Abs(fMCparticleGGMother->GetPdgCode());
4219                                 if(fMCparticleGGMother->GetMother()<0)
4220                                 {
4221                                         gggmpdg = 0;
4222                                 }
4223                                 else
4224                                 {
4225                                         fMCparticleGGGMother = (AliAODMCParticle*) fMCarray->At(fMCparticleGGMother->GetMother());
4226                                         gggmpdg = TMath::Abs(fMCparticleGGGMother->GetPdgCode());
4227                                 }
4228                         }
4229                 }
4230         }
4231         else
4232         {
4233                 fMCtrack = fMCstack->Particle(mcIndex);
4234
4235                 pdg = TMath::Abs(fMCtrack->GetPdgCode());
4236
4237                 if(pdg!=11)
4238                 {
4239                         fIsHFE1 = kFALSE;
4240                         fIsHFE2 = kFALSE;
4241                         fIsNonHFE = kFALSE;
4242                         fIsFromD = kFALSE;
4243                         fIsFromB = kFALSE;
4244                         fIsFromPi0 = kFALSE;
4245                         fIsFromEta = kFALSE;
4246                         fIsFromGamma = kFALSE;
4247                         return kFALSE;
4248                 }
4249
4250                 if(fMCtrack->GetFirstMother()<0)
4251                 {
4252                         fIsHFE1 = kFALSE;
4253                         fIsHFE2 = kFALSE;
4254                         fIsNonHFE = kFALSE;
4255                         fIsFromD = kFALSE;
4256                         fIsFromB = kFALSE;
4257                         fIsFromPi0 = kFALSE;
4258                         fIsFromEta = kFALSE;
4259                         fIsFromGamma = kFALSE;
4260                         return kFALSE;
4261                 }
4262
4263                 fMCtrackMother = fMCstack->Particle(fMCtrack->GetFirstMother());
4264                 mpdg = TMath::Abs(fMCtrackMother->GetPdgCode());
4265
4266                 if(fMCtrackMother->GetFirstMother()<0)
4267                 {
4268                         gmpdg = 0;
4269                         ggmpdg = 0;
4270                         gggmpdg = 0;
4271                 }
4272                 else
4273                 {
4274                         fMCtrackGMother = fMCstack->Particle(fMCtrackMother->GetFirstMother());
4275                         gmpdg = TMath::Abs(fMCtrackGMother->GetPdgCode());
4276
4277                         if(fMCtrackGMother->GetFirstMother()<0)
4278                         {
4279                                 ggmpdg = 0;
4280                                 gggmpdg = 0;
4281                         }
4282                         else
4283                         {
4284                                 fMCtrackGGMother = fMCstack->Particle(fMCtrackGMother->GetFirstMother());
4285                                 ggmpdg = TMath::Abs(fMCtrackGGMother->GetPdgCode());
4286
4287                                 if(fMCtrackGGMother->GetFirstMother()<0)
4288                                 {
4289                                         gggmpdg = 0;
4290                                 }
4291                                 else
4292                                 {
4293                                         fMCtrackGGGMother = fMCstack->Particle(fMCtrackGGMother->GetFirstMother());
4294                                         gggmpdg = TMath::Abs(fMCtrackGGGMother->GetPdgCode());
4295                                 }
4296                         }
4297                 }
4298         }
4299
4300         //Tag Electron Source
4301         if(mpdg==111 || mpdg==221 || mpdg==22)
4302         {
4303                 fIsHFE1 = kFALSE;
4304                 fIsHFE2 = kFALSE;
4305                 fIsNonHFE = kTRUE;
4306                 fIsFromD = kFALSE;
4307                 fIsFromB = kFALSE;
4308
4309                 fIsFromPi0 = kFALSE;
4310                 fIsFromEta = kFALSE;
4311                 fIsFromGamma = kFALSE;
4312
4313                 if(mpdg==111) fIsFromPi0 = kFALSE;
4314                 if(mpdg==221)fIsFromEta = kFALSE;
4315                 if(mpdg==22) fIsFromGamma = kFALSE;
4316
4317                 return kTRUE;
4318         }
4319         else
4320         {
4321                 fIsHFE1 = kFALSE;
4322                 fIsHFE2 = kTRUE;
4323
4324                 fIsFromPi0 = kFALSE;
4325                 fIsFromEta = kFALSE;
4326                 fIsFromGamma = kFALSE;
4327
4328                 fIsNonHFE = kFALSE;
4329
4330                 fIsFromD = kFALSE;
4331                 fIsFromB = kFALSE;
4332
4333                 if(mpdg>400 && mpdg<500)
4334                 {
4335                         if((gmpdg>500 && gmpdg<600) || (ggmpdg>500 && ggmpdg<600) || (gggmpdg>500 && gggmpdg<600))
4336                         {
4337                                 fIsHFE1 = kTRUE;
4338                                 fIsFromD = kFALSE;
4339                                 fIsFromB = kTRUE;
4340                                 return kTRUE;
4341                         }
4342                         else
4343                         {
4344                                 fIsHFE1 = kTRUE;
4345                                 fIsFromD = kTRUE;
4346                                 fIsFromB = kFALSE;
4347                                 return kTRUE;
4348                         }
4349                 }
4350                 else if(mpdg>500 && mpdg<600)
4351                 {
4352                         fIsHFE1 = kTRUE;
4353                         fIsFromD = kFALSE;
4354                         fIsFromB = kTRUE;
4355                         return kTRUE;
4356                 }
4357                 else
4358                 {
4359                         fIsHFE1 = kFALSE;
4360                         fIsFromD = kFALSE;
4361                         fIsFromB = kFALSE;
4362                         return kFALSE;
4363                 }
4364         }
4365 }
4366 /*
4367 Bool_t AliAnalysisTaskEMCalHFEpA::ContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
4368 {
4369                 // Check that in the cluster cells, there is no bad channel of those stored
4370                 // in fEMCALBadChannelMap
4371                 // adapted from AliCalorimeterUtils
4372
4373                 //if (!fRemoveBadChannels) return kFALSE;
4374                 //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
4375         if( !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
4376
4377                 //Int_t icol = -1;
4378                 //Int_t irow = -1;
4379                 //Int_t imod = -1;
4380         for(Int_t iCell = 0; iCell<nCells; iCell++){
4381
4382                         //Get the column and row
4383                 if(calorimeter == "EMCAL"){
4384                         return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
4385                 }
4386                 else return kFALSE;
4387
4388         }// cell cluster loop
4389
4390         return kFALSE;
4391
4392 }
4393 */
4394 /*
4395
4396 //________________________________________________________________________________
4397 TArrayI AliAnalysisTaskEMCalHFEpA::GetTriggerPatches(Bool_t IsEventEMCALL0, Bool_t IsEventEMCALL1)
4398 {
4399         // Select the patches that triggered
4400         // Depend on L0 or L1
4401
4402         // some variables
4403         //Int_t  trigtimes[30], globCol, globRow,ntimes, i;
4404         Int_t   globCol, globRow;
4405
4406         Int_t  absId  = -1; //[100];
4407         Int_t  nPatch = 0;
4408
4409         TArrayI patches(0);
4410
4411                 // get object pointer
4412         AliVCaloTrigger *caloTrigger = InputEvent()->GetCaloTrigger( "EMCAL" );
4413
4414                 // class is not empty
4415         if( caloTrigger->GetEntries() > 0 )
4416         {
4417                         // must reset before usage, or the class will fail
4418                 caloTrigger->Reset();
4419
4420                         // go throuth the trigger channels
4421                 while( caloTrigger->Next() )
4422                 {
4423                                 // get position in global 2x2 tower coordinates
4424                         caloTrigger->GetPosition( globCol, globRow );
4425
4426                         //L0
4427                         if(IsEventEMCALL0)
4428                         {
4429                             //not implemented
4430                         }
4431
4432
4433                         else if(IsEventEMCALL1) // L1
4434                         {
4435                                         Int_t bit = 0;
4436                                         caloTrigger->GetTriggerBits(bit);
4437
4438                                         Bool_t isEGA = ((bit >> fBitEGA) & 0x1);
4439                                         //Bool_t isEJE = ((bit >> fBitEJE) & 0x1) && IsEventEMCALL1Jet  () ;
4440
4441                                         if(!isEGA) continue;
4442
4443                                         Int_t patchsize = -1;
4444                                         if(isEGA) patchsize =  2;
4445                                         //else if (isEJE) patchsize = 16;
4446
4447                                         // add 2x2 (EGA) or 16x16 (EJE) patches
4448                                         for(Int_t irow=0; irow < patchsize; irow++)
4449                                         {
4450                                                 for(Int_t icol=0; icol < patchsize; icol++)
4451                                                 {
4452                                                         GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol+icol,globRow+irow, absId);
4453                                                         //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absIDTrig[nTrig]);
4454                                                         patches.Set(nPatch+1); // Set size of this array to nPatch+1 ints.
4455                                                         patches.AddAt(absId,nPatch++); //Add Int_t absId at position nPatch++
4456                                                 }
4457                                         }
4458
4459                         } // L1
4460
4461                 } // trigger iterator
4462         } // go thorough triggers
4463
4464         printf("N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nPatch, patches.At(0), patches.At(patches.GetSize()-1));
4465
4466         return patches;
4467 }
4468  */
4469 Double_t AliAnalysisTaskEMCalHFEpA::CalculateWeight(Int_t pdg_particle, Double_t x)
4470 {
4471                 //weight for d3 based on MinJung parametrization //sent by Jan
4472             Double_t weight=1;
4473                 if(pdg_particle==111){
4474                         if(x>= 0.100000 &&  x < 0.112797 ) weight=1.262120;
4475                         if(x>= 0.112797 &&  x < 0.127231 ) weight=1.277765;
4476                         if(x>= 0.127231 &&  x < 0.143512 ) weight=1.295605;
4477                         if(x>= 0.143512 &&  x < 0.161877 ) weight=1.318155;
4478                         if(x>= 0.161877 &&  x < 0.182592 ) weight=1.348693;
4479                         if(x>= 0.182592 &&  x < 0.205957 ) weight=1.388636;
4480                         if(x>= 0.205957 &&  x < 0.232313 ) weight=1.439122;
4481                         if(x>= 0.232313 &&  x < 0.262041 ) weight=1.497452;
4482                         if(x>= 0.262041 &&  x < 0.295573 ) weight=1.559409;
4483                         if(x>= 0.295573 &&  x < 0.333397 ) weight=1.615169;
4484                         if(x>= 0.333397 &&  x < 0.376060 ) weight=1.654954;
4485                         if(x>= 0.376060 &&  x < 0.424183 ) weight=1.668753;
4486                         if(x>= 0.424183 &&  x < 0.478465 ) weight=1.652225;
4487                         if(x>= 0.478465 &&  x < 0.539692 ) weight=1.603119;
4488                         if(x>= 0.539692 &&  x < 0.608754 ) weight=1.526049;
4489                         if(x>= 0.608754 &&  x < 0.686654 ) weight=1.426724;
4490                         if(x>= 0.686654 &&  x < 0.774523 ) weight=1.312684;
4491                         if(x>= 0.774523 &&  x < 0.873636 ) weight=1.195395;
4492                         if(x>= 0.873636 &&  x < 0.985432 ) weight=1.086264;
4493                         if(x>= 0.985432 &&  x < 1.111534 ) weight=0.993666;
4494                         if(x>= 1.111534 &&  x < 1.253773 ) weight=0.922587;
4495                         if(x>= 1.253773 &&  x < 1.414214 ) weight=0.875739;
4496                         if(x>= 1.414214 &&  x < 1.595185 ) weight=0.852181;
4497                         if(x>= 1.595185 &&  x < 1.799315 ) weight=0.847828;
4498                         if(x>= 1.799315 &&  x < 2.029567 ) weight=0.863875;
4499                         if(x>= 2.029567 &&  x < 2.289283 ) weight=0.899112;
4500                         if(x>= 2.289283 &&  x < 2.582235 ) weight=0.955194;
4501                         if(x>= 2.582235 &&  x < 2.912674 ) weight=1.033824;
4502                         if(x>= 2.912674 &&  x < 3.285398 ) weight=1.133714;
4503                         if(x>= 3.285398 &&  x < 3.705818 ) weight=1.259471;
4504                         if(x>= 3.705818 &&  x < 4.180038 ) weight=1.406883;
4505                         if(x>= 4.180038 &&  x < 4.714942 ) weight=1.578923;
4506                         if(x>= 4.714942 &&  x < 5.318296 ) weight=1.778513;
4507                         if(x>= 5.318296 &&  x < 5.998859 ) weight=2.001171;
4508                         if(x>= 5.998859 &&  x < 6.766511 ) weight=2.223161;
4509                         if(x>= 6.766511 &&  x < 7.632396 ) weight=2.449445;
4510                         if(x>= 7.632396 &&  x < 8.609086 ) weight=2.661734;
4511                         if(x>= 8.609086 &&  x < 9.710759 ) weight=2.851935;
4512                         if(x>= 9.710759 &&  x < 10.953409 ) weight=2.974319;
4513                         if(x>= 10.953409 &&  x < 12.355077 ) weight=3.106314;
4514                         if(x>= 12.355077 &&  x < 13.936111 ) weight=3.126815;
4515                         if(x>= 13.936111 &&  x < 15.719464 ) weight=3.150053;
4516                         if(x>= 15.719464 &&  x < 17.731026 ) weight=3.218509;
4517                         if(x>= 17.731026 &&  x < 20.000000 ) weight=3.252141;
4518
4519                 }
4520             else if(pdg_particle==221){
4521                    weight=1;
4522                 }
4523             else weight=1;
4524
4525             return weight;
4526
4527 }
4528
4529
4530
4531