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