PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx

   1 // $Id$
   2
   3 #include "AliAnalysisTaskEMCALIsoPhoton.h"
   4
   5 #include <TCanvas.h>
   6 #include <TChain.h>
   7 #include <TFile.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <THnSparse.h>
  12 #include <TLorentzVector.h>
  13 #include <TList.h>
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliAnalysisTask.h"
  17 #include "AliEMCALGeometry.h"
  18 #include "AliEMCALRecoUtils.h"
  19 #include "AliESDCaloCells.h"
  20 #include "AliESDCaloCluster.h"
  21 #include "AliESDEvent.h"
  22 #include "AliESDHeader.h"
  23 #include "AliESDInputHandler.h"
  24 #include "AliESDUtils.h"
  25 #include "AliESDtrack.h"
  26 #include "AliESDtrackCuts.h"
  27 #include "AliAODEvent.h"
  28 #include "AliAODTrack.h"
  29 #include "AliMCEvent.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliStack.h"
  32 #include "AliAODMCParticle.h"
  33 #include "AliVEvent.h"
  34 #include "AliVTrack.h"
  35 #include "AliV0vertexer.h"
  36 #include "AliVCluster.h"
  37 #include "AliOADBContainer.h"
  38
  39 #include <iostream>
  40 using std::cout;
  41 using std::endl;
  42
  43 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
  44
  45 //________________________________________________________________________
  46 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
  47   AliAnalysisTaskSE(),
  48   fESDClusters(0),
  49   fAODClusters(0),
  50   fSelPrimTracks(0),
  51   fTracks(0),
  52   fAODMCParticles(0),
  53   fESDCells(0),
  54   fAODCells(0),
  55   fPrTrCuts(0),
  56   fGeom(0x0),
  57   fGeoName("EMCAL_COMPLETEV1"),
  58   fOADBContainer(0),
  59   fPeriod("LHC11c"),
  60   fTrigBit("kEMC7"),
  61   fIsTrain(0),
  62   fIsMc(0),
  63   fDebug(0),
  64   fPathStrOpt("/"),
  65   fExoticCut(0.97),
  66   fIsoConeR(0.4),
  67   fNDimensions(7),
  68   fECut(3.),
  69   fTrackMult(0),
  70   fMcIdFamily(""),
  71   fNClusForDirPho(0),
  72   fDirPhoPt(0),
  73   fHigherPtCone(0),
  74   fImportGeometryFromFile(0),
  75   fImportGeometryFilePath(""),
  76   fMaxPtTrack(0),
  77   fMaxEClus(0),
  78   fNCells50(0),
  79   fFilterBit(0),
  80   fSelHybrid(kFALSE),
  81   fFillQA(kFALSE),
  82   fClusIdFromTracks(""),
  83   fCpvFromTrack(kFALSE),
  84   fNBinsPt(200),
  85   fPtBinLowEdge(-0.25),
  86   fPtBinHighEdge(99.75),
  87   fRemMatchClus(kFALSE),
  88   fMinIsoClusE(0),
  89   fNCuts(5),
  90   fTrCoreRem(kTRUE),
  91   fClusTDiff(30e-9),
  92   fESD(0),
  93   fAOD(0),
  94   fVEvent(0),
  95   fMCEvent(0),
  96   fStack(0),
  97   fOutputList(0),
  98   fEvtSel(0),
  99   fNClusEt10(0),
 100   fRecoPV(0),
 101   fPVtxZ(0),
 102   fTrMultDist(0),
 103   fClusEtCPVSBGISO(0),
 104   fClusEtCPVBGISO(0),
 105   fMCDirPhotonPtEtaPhi(0),
 106   fMCIsoDirPhotonPtEtaPhi(0),
 107   fMCDirPhotonPtEtIso(0),
 108   fDecayPhotonPtMC(0),
 109   fCellAbsIdVsAmpl(0),
 110   fNClusHighClusE(0),
 111   fHigherPtConeM02(0),
 112   fClusEtMcPt(0),
 113   fClusMcDetaDphi(0),
 114   fNClusPerPho(0),
 115   fMcPtInConeBG(0),
 116   fMcPtInConeSBG(0),
 117   fMcPtInConeBGnoUE(0),
 118   fMcPtInConeSBGnoUE(0),
 119   fMcPtInConeTrBGnoUE(0),
 120   fMcPtInConeTrSBGnoUE(0),
 121   fMcPtInConeMcPhoPt(0),
 122   fAllIsoEtMcGamma(0),
 123   fAllIsoNoUeEtMcGamma(0),
 124   fMCDirPhotonPtEtaPhiNoClus(0),
 125   fHnOutput(0),
 126   fQAList(0),
 127   fNTracks(0),
 128   fEmcNCells(0),
 129   fEmcNClus(0),
 130   fEmcNClusCut(0),
 131   fNTracksECut(0),
 132   fEmcNCellsCut(0),
 133   fEmcClusETM1(0),
 134   fEmcClusETM2(0),
 135   fEmcClusNotExo(0),
 136   fEmcClusEClusCuts(0),
 137   fEmcClusEPhi(0),
 138   fEmcClusEPhiCut(0),
 139   fEmcClusEEta(0),
 140   fEmcClusEEtaCut(0),
 141   fTrackPtPhi(0),
 142   fTrackPtPhiCut(0),
 143   fTrackPtEta(0),
 144   fTrackPtEtaCut(0),
 145   fMaxCellEPhi(0)
 146 {
 147   // Default constructor.
 148   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 149 }
 150
 151 //________________________________________________________________________
 152 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
 153   AliAnalysisTaskSE(name),
 154   fESDClusters(0),
 155   fAODClusters(0),
 156   fSelPrimTracks(0),
 157   fTracks(0),
 158   fAODMCParticles(0),
 159   fESDCells(0),
 160   fAODCells(0),
 161   fPrTrCuts(0),
 162   fGeom(0x0),
 163   fGeoName("EMCAL_COMPLETEV1"),
 164   fOADBContainer(0),
 165   fPeriod("LHC11c"),
 166   fTrigBit("kEMC7"),
 167   fIsTrain(0),
 168   fIsMc(0),
 169   fDebug(0),
 170   fPathStrOpt("/"),
 171   fExoticCut(0.97),
 172   fIsoConeR(0.4),
 173   fNDimensions(7),
 174   fECut(3.),
 175   fTrackMult(0),
 176   fMcIdFamily(""),
 177   fNClusForDirPho(0),
 178   fDirPhoPt(0),
 179   fHigherPtCone(0),
 180   fImportGeometryFromFile(0),
 181   fImportGeometryFilePath(""),
 182   fMaxPtTrack(0),
 183   fMaxEClus(0),
 184   fNCells50(0),
 185   fFilterBit(0),
 186   fSelHybrid(kFALSE),
 187   fFillQA(kFALSE),
 188   fClusIdFromTracks(""),
 189   fCpvFromTrack(kFALSE),
 190   fNBinsPt(200),
 191   fPtBinLowEdge(-0.25),
 192   fPtBinHighEdge(99.75),
 193   fRemMatchClus(kFALSE),
 194   fMinIsoClusE(0),
 195   fNCuts(5),
 196   fTrCoreRem(kTRUE),
 197   fClusTDiff(30e-9),
 198   fESD(0),
 199   fAOD(0),
 200   fVEvent(0),
 201   fMCEvent(0),
 202   fStack(0),
 203   fOutputList(0),
 204   fEvtSel(0),
 205   fNClusEt10(0),
 206   fRecoPV(0),
 207   fPVtxZ(0),
 208   fTrMultDist(0),
 209   fClusEtCPVSBGISO(0),
 210   fClusEtCPVBGISO(0),
 211   fMCDirPhotonPtEtaPhi(0),
 212   fMCIsoDirPhotonPtEtaPhi(0),
 213   fMCDirPhotonPtEtIso(0),
 214   fDecayPhotonPtMC(0),
 215   fCellAbsIdVsAmpl(0),
 216   fNClusHighClusE(0),
 217   fHigherPtConeM02(0),
 218   fClusEtMcPt(0),
 219   fClusMcDetaDphi(0),
 220   fNClusPerPho(0),
 221   fMcPtInConeBG(0),
 222   fMcPtInConeSBG(0),
 223   fMcPtInConeBGnoUE(0),
 224   fMcPtInConeSBGnoUE(0),
 225   fMcPtInConeTrBGnoUE(0),
 226   fMcPtInConeTrSBGnoUE(0),
 227   fMcPtInConeMcPhoPt(0),
 228   fAllIsoEtMcGamma(0),
 229   fAllIsoNoUeEtMcGamma(0),
 230   fMCDirPhotonPtEtaPhiNoClus(0),
 231   fHnOutput(0),
 232   fQAList(0),
 233   fNTracks(0),
 234   fEmcNCells(0),
 235   fEmcNClus(0),
 236   fEmcNClusCut(0),
 237   fNTracksECut(0),
 238   fEmcNCellsCut(0),
 239   fEmcClusETM1(0),
 240   fEmcClusETM2(0),
 241   fEmcClusNotExo(0),
 242   fEmcClusEClusCuts(0),
 243   fEmcClusEPhi(0),
 244   fEmcClusEPhiCut(0),
 245   fEmcClusEEta(0),
 246   fEmcClusEEtaCut(0),
 247   fTrackPtPhi(0),
 248   fTrackPtPhiCut(0),
 249   fTrackPtEta(0),
 250   fTrackPtEtaCut(0),
 251   fMaxCellEPhi(0)
 252 {
 253   // Constructor
 254
 255   // Define input and output slots here
 256   // Input slot #0 works with a TChain
 257   DefineInput(0, TChain::Class());
 258   // Output slot #0 id reserved by the base class for AOD
 259   // Output slot #1 writes into a TH1 container
 260   DefineOutput(1, TList::Class());
 261   DefineOutput(2, TList::Class());
 262 }
 263
 264 //________________________________________________________________________
 265 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 266 {
 267   // Create histograms, called once.
 268
 269   fESDClusters = new TObjArray();
 270   fSelPrimTracks = new TObjArray();
 271
 272
 273   fOutputList = new TList();
 274   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 275
 276   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
 277   fOADBContainer = new AliOADBContainer("AliEMCALgeo");
 278   fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
 279
 280   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 281   fOutputList->Add(fEvtSel);
 282
 283   fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
 284   fOutputList->Add(fNClusEt10);
 285
 286   fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
 287   fOutputList->Add(fRecoPV);
 288
 289   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 290   fOutputList->Add(fPVtxZ);
 291
 292   fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5);
 293   fOutputList->Add(fTrMultDist);
 294
 295   fClusEtCPVSBGISO = new TH2F("hClusEtCPVSBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.1<#lambda_{0}^{2}<0.3;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100);
 296   fClusEtCPVSBGISO->Sumw2();
 297   fOutputList->Add(fClusEtCPVSBGISO);
 298
 299   fClusEtCPVBGISO = new TH2F("hClusEtCPVBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.5<#lambda_{0}^{2}<2.0;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100);
 300   fClusEtCPVBGISO->Sumw2();
 301   fOutputList->Add(fClusEtCPVBGISO);
 302
 303   fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",100,-0.5,99.5,154,-0.77,0.77,130,1.38,3.20);
 304   fMCDirPhotonPtEtaPhi->Sumw2();
 305   fOutputList->Add(fMCDirPhotonPtEtaPhi);
 306
 307   fMCIsoDirPhotonPtEtaPhi = new TH3F("hMCIsoDirPhotonPtEtaPhi","photon (gq->#gammaq, isolated@MC) p_{T}, #eta, #phi;GeV/c;#eta;#phi",100,-0.5,99.5,154,-0.77,0.77,130,1.38,3.20);
 308   fMCIsoDirPhotonPtEtaPhi->Sumw2();
 309   fOutputList->Add(fMCIsoDirPhotonPtEtaPhi);
 310
 311   fMCDirPhotonPtEtIso = new TH2F("hMCDirPhotonPtEtIso",Form("photon (gq->#gammaq @MC) p_{T}, E_{T}^{ISO} (R=%1.1f);GeV/c;E_{T}^{ISO} GeV/c",fIsoConeR),100,-0.5,99.5,20,-0.25,9.75);
 312   fMCDirPhotonPtEtIso->Sumw2();
 313   fOutputList->Add(fMCDirPhotonPtEtIso);
 314
 315
 316   fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 317   fDecayPhotonPtMC->Sumw2();
 318   fOutputList->Add(fDecayPhotonPtMC);
 319
 320   fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,24*48*10,-0.5,24*48*10-0.5);
 321   fOutputList->Add(fCellAbsIdVsAmpl);
 322
 323   fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5);
 324   fOutputList->Add(fNClusHighClusE);
 325
 326   fHigherPtConeM02 = new TH2F("hHigherPtConeM02","#lambda_{0}^{2} vs. in-cone-p_{T}^{max};p_{T}^{max} (GeV/c, in the cone);#lambda_{0}^{2}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,400,0,4);
 327   fOutputList->Add(fHigherPtConeM02);
 328
 329   fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100);
 330   fOutputList->Add(fClusEtMcPt);
 331
 332   fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7);
 333   fOutputList->Add(fClusMcDetaDphi);
 334
 335   fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5);
 336   fOutputList->Add(fNClusPerPho);
 337
 338   fMcPtInConeBG = new TH2F("hMcPtInConeBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 339   fOutputList->Add(fMcPtInConeBG);
 340
 341   fMcPtInConeSBG  = new TH2F("hMcPtInConeSBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 342   fOutputList->Add(fMcPtInConeSBG);
 343
 344   fMcPtInConeBGnoUE = new TH2F("hMcPtInConeBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 345   fOutputList->Add(fMcPtInConeBGnoUE);
 346
 347   fMcPtInConeSBGnoUE  = new TH2F("hMcPtInConeSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 348   fOutputList->Add(fMcPtInConeSBGnoUE);
 349
 350   fMcPtInConeTrBGnoUE = new TH2F("hMcPtInConeTrBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (BG template);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 351   fOutputList->Add(fMcPtInConeTrBGnoUE);
 352
 353   fMcPtInConeTrSBGnoUE  = new TH2F("hMcPtInConeTrSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (SBG range);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 354   fOutputList->Add(fMcPtInConeTrSBGnoUE);
 355
 356   fMcPtInConeMcPhoPt  = new TH2F("hMcPtInConeMcPhoPt","#sum_{in-cone}p_{T}^{mc-primaries} vs. prompt photon p_{T};p_{T}^{mc-#gamma} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,200,-0.25,99.75);
 357   fOutputList->Add(fMcPtInConeMcPhoPt);
 358
 359   fAllIsoEtMcGamma  = new TH2F("hAllIsoEtMcGammaE","ISO^{TRK+EMC} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50);
 360   fOutputList->Add(fAllIsoEtMcGamma);
 361
 362   fAllIsoNoUeEtMcGamma  = new TH2F("hAllIsoNoUeEtMcGammaE","ISO^{TRK+EMC}_{noue} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC}_{noue} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50);
 363   fOutputList->Add(fAllIsoNoUeEtMcGamma);
 364
 365
 366   fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and  #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20);
 367   fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus);
 368
 369   Int_t nEt=fNBinsPt*5, nM02=400, nCeIso=1000, nTrIso=1000,  nAllIso=1000,  nCeIsoNoUE=1000,  nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=100, nPhoMcPt=100;
 370   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt};
 371   fNDimensions = sizeof(bins)/sizeof(Int_t);
 372   const Int_t ndims =   fNDimensions;
 373   Double_t xmin[] = { fPtBinLowEdge,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-0.5};
 374   Double_t xmax[] = { fPtBinHighEdge, 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5};
 375   if(fPeriod.Contains("11h")){
 376     xmax[12]=3999.5;
 377   }
 378   fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult,mc-p_{T}^{#gamma}", ndims, bins, xmin, xmax);
 379   fOutputList->Add(fHnOutput);
 380
 381   //QA outputs
 382   fQAList = new TList();
 383   fQAList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 384
 385   fNTracks = new TH1F("hNTracks","# of selected tracks;n-tracks;counts",120,-0.5,119.5);
 386   fNTracks->Sumw2();
 387   fQAList->Add(fNTracks);
 388
 389   fEmcNCells = new TH1F("fEmcNCells",";n/event;count",120,-0.5,119.5);
 390   fEmcNCells->Sumw2();
 391   fQAList->Add(fEmcNCells);
 392   fEmcNClus = new TH1F("fEmcNClus",";n/event;count",120,-0.5,119.5);
 393   fEmcNClus->Sumw2();
 394   fQAList->Add(fEmcNClus);
 395   fEmcNClusCut = new TH1F("fEmcNClusCut",";n/event;count",120,-0.5,119.5);
 396   fEmcNClusCut->Sumw2();
 397   fQAList->Add(fEmcNClusCut);
 398   fNTracksECut = new TH1F("fNTracksECut",";n/event;count",120,-0.5,119.5);
 399   fNTracksECut->Sumw2();
 400   fQAList->Add(fNTracksECut);
 401   fEmcNCellsCut = new TH1F("fEmcNCellsCut",";n/event;count",120,-0.5,119.5);
 402   fEmcNCellsCut->Sumw2();
 403   fQAList->Add(fEmcNCellsCut);
 404   fEmcClusETM1 = new TH1F("fEmcClusETM1","(method clus->GetTrackDx,z);GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 405   fEmcClusETM1->Sumw2();
 406   fQAList->Add(fEmcClusETM1);
 407   fEmcClusETM2 = new TH1F("fEmcClusETM2","(method track->GetEMCALcluster());GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 408   fEmcClusETM2->Sumw2();
 409   fQAList->Add(fEmcClusETM2);
 410   fEmcClusNotExo  = new TH1F("fEmcClusNotExo","exotics removed;GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 411   fEmcClusNotExo->Sumw2();
 412   fQAList->Add(fEmcClusNotExo);
 413   fEmcClusEPhi = new TH2F("fEmcClusEPhi",";GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 414   fEmcClusEPhi->Sumw2();
 415   fQAList->Add(fEmcClusEPhi);
 416   fEmcClusEPhiCut = new TH2F("fEmcClusEPhiCut",";GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 417   fEmcClusEPhiCut->Sumw2();
 418   fQAList->Add(fEmcClusEPhiCut);
 419   fEmcClusEEta = new TH2F("fEmcClusEEta",";GeV;#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,19,-0.9,0.9);
 420   fEmcClusEEta->Sumw2();
 421   fQAList->Add(fEmcClusEEta);
 422   fEmcClusEEtaCut = new TH2F("fEmcClusEEtaCut",";GeV;#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
 423   fEmcClusEEtaCut->Sumw2();
 424   fQAList->Add(fEmcClusEEtaCut);
 425   fTrackPtPhi = new TH2F("fTrackPtPhi",";p_{T} [GeV/c];#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 426   fTrackPtPhi->Sumw2();
 427   fQAList->Add(fTrackPtPhi);
 428   fTrackPtPhiCut = new TH2F("fTrackPtPhiCut",";p_{T} [GeV/c];#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 429   fTrackPtPhiCut->Sumw2();
 430   fQAList->Add(fTrackPtPhiCut);
 431   fTrackPtEta = new TH2F("fTrackPtEta",";p_{T} [GeV/c];#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
 432   fTrackPtEta->Sumw2();
 433   fQAList->Add(fTrackPtEta);
 434   fTrackPtEtaCut = new TH2F("fTrackPtEtaCut",";p_{T} [GeV/c];#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
 435   fTrackPtEtaCut->Sumw2();
 436   fQAList->Add(fTrackPtEtaCut);
 437   fEmcClusEClusCuts = new TH2F("fEmcClusEClusCuts",";GeV;cut",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,fNCuts,-0.5,fNCuts-0.5);
 438   fEmcClusEClusCuts->Sumw2();
 439   fQAList->Add(fEmcClusEClusCuts);
 440
 441   fMaxCellEPhi = new TH2F("fMaxCellEPhi","Most energetic cell in event; GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 442   fMaxCellEPhi->Sumw2();
 443   fQAList->Add(fMaxCellEPhi);
 444
 445   PostData(1, fOutputList);
 446   PostData(2, fQAList);
 447 }
 448
 449 //________________________________________________________________________
 450 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 451 {
 452   // Main loop, called for each event.
 453   fESDClusters = 0;
 454   fESDCells = 0;
 455   fAODClusters = 0;
 456   fAODCells = 0;
 457   // event trigger selection
 458   Bool_t isSelected = 0;
 459   if(fPeriod.Contains("11a")){
 460     if(fTrigBit.Contains("kEMC"))
 461       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 462     else
 463       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 464   }
 465   else{
 466     if(fTrigBit.Contains("kEMC"))
 467       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 468     else
 469       if(fTrigBit.Contains("kMB"))
 470         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 471       else
 472         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
 473   }
 474   if(fPeriod.Contains("11h")){
 475     if(fTrigBit.Contains("kAny"))
 476       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
 477     if(fTrigBit.Contains("kAnyINT"))
 478       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT);
 479   }
 480   if(fIsMc)
 481     isSelected = kTRUE;
 482   if(fDebug)
 483     printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
 484   if(!isSelected )
 485         return;
 486   if(fIsMc){
 487     TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
 488     TFile *file = (TFile*)tree->GetCurrentFile();
 489     TString filename = file->GetName();
 490     if(!filename.Contains(fPathStrOpt.Data()))
 491       return;
 492   }
 493   fVEvent = (AliVEvent*)InputEvent();
 494   if (!fVEvent) {
 495     printf("ERROR: event not available\n");
 496     return;
 497   }
 498   Int_t   runnumber = InputEvent()->GetRunNumber() ;
 499   if(fDebug)
 500     printf("run number = %d\n",runnumber);
 501
 502   fESD = dynamic_cast<AliESDEvent*>(fVEvent);
 503   if(!fESD){
 504     fAOD = dynamic_cast<AliAODEvent*>(fVEvent);
 505     if(!fAOD){
 506       printf("ERROR: Invalid type of event!!!\n");
 507       return;
 508     }
 509     else if(fDebug)
 510       printf("AOD EVENT!\n");
 511   }
 512
 513   fEvtSel->Fill(0);
 514   if(fDebug)
 515     printf("event is ok,\n run number=%d\n",runnumber);
 516
 517
 518   AliVVertex *pv = (AliVVertex*)fVEvent->GetPrimaryVertex();
 519   Bool_t pvStatus = kTRUE;
 520   if(fESD){
 521     AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex();
 522     pvStatus = esdv->GetStatus();
 523   }
 524   /*if(fAOD){
 525     AliAODVertex *aodv = (AliAODVertex*)fAOD->GetPrimaryVertex();
 526     pvStatus = aodv->GetStatus();
 527     }*/
 528   if(!pv)
 529     return;
 530   if(!pvStatus)
 531     fRecoPV->Fill(0);
 532   else
 533     fRecoPV->Fill(1);
 534   fPVtxZ->Fill(pv->GetZ());
 535   if(TMath::Abs(pv->GetZ())>10)
 536     return;
 537   if(fDebug)
 538     printf("passed vertex cut\n");
 539
 540   fEvtSel->Fill(1);
 541   if(fESD)
 542     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 543   if(fAOD)
 544     fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
 545
 546   if(!fTracks){
 547     AliError("Track array in event is NULL!");
 548     if(fDebug)
 549       printf("returning due to not finding tracks!\n");
 550     return;
 551   }
 552   // Track loop to fill a pT spectrum
 553   const Int_t Ntracks = fTracks->GetEntriesFast();
 554   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 555     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 556     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 557     AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
 558     if (!track)
 559       continue;
 560     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
 561     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 562     if(esdTrack){
 563       if(esdTrack->GetEMCALcluster()>0)
 564         fClusIdFromTracks.Append(Form("%d ",esdTrack->GetEMCALcluster()));
 565       if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
 566         fSelPrimTracks->Add(track);
 567       }
 568     }
 569     else if(aodTrack){
 570       if (fSelHybrid && !aodTrack->IsHybridGlobalConstrainedGlobal())
 571         continue ;
 572       if(!fSelHybrid && !aodTrack->TestFilterBit(fFilterBit))
 573         continue;
 574       fSelPrimTracks->Add(track);
 575       /*if(fTrackMaxPt<track->Pt())
 576         fTrackMaxPt = track->Pt();*/
 577     }
 578   }
 579
 580   TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
 581   for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 582     if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 583       break;
 584     /*if(fESD)
 585       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 586       else*/
 587     // if(fVEvent->GetEMCALMatrix(mod))
 588     fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
 589     fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
 590   }
 591
 592   if(fESD){
 593     AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 594     fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 595     if(fDebug)
 596       printf("ESD Track mult= %d\n",fTrackMult);
 597   }
 598   else if(fAOD){
 599     fTrackMult = Ntracks;
 600     if(fDebug)
 601       printf("AOD Track mult= %d\n",fTrackMult);
 602   }
 603   fTrMultDist->Fill(fTrackMult);
 604
 605   if(fESD){
 606     TList *l = fESD->GetList();
 607     fESDClusters =  dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 608     fESDCells = fESD->GetEMCALCells();
 609     if(fDebug)
 610       printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast());
 611   }
 612   else if(fAOD){
 613     fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
 614     fAODCells = fAOD->GetEMCALCells();
 615     if(fDebug)
 616       printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast());
 617   }
 618
 619
 620   fMCEvent = MCEvent();
 621   if(fMCEvent){
 622     if(fDebug)
 623       std::cout<<"MCevent exists\n";
 624     fStack = (AliStack*)fMCEvent->Stack();
 625     if(!fStack)
 626       fAODMCParticles = (TClonesArray*)fVEvent->FindListObject("mcparticles");
 627   }
 628   else{
 629     if(fDebug)
 630       std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
 631   }
 632   LoopOnCells();
 633   FollowGamma();
 634   if(fDebug)
 635     printf("passed calling of FollowGamma\n");
 636   FillClusHists();
 637   if(fDebug)
 638     printf("passed calling of FillClusHists\n");
 639   FillMcHists();
 640   if(fDebug)
 641     printf("passed calling of FillMcHists\n");
 642   if(fFillQA)
 643     FillQA();
 644   if(fDebug)
 645     printf("passed calling of FillQA\n");
 646   /*if(fESD)
 647     fESDClusters->Clear();*/
 648   fSelPrimTracks->Clear();
 649   fNClusForDirPho = 0;
 650   fNCells50 = 0;
 651   fClusIdFromTracks = "";
 652
 653   PostData(1, fOutputList);
 654   PostData(2, fQAList);
 655 }
 656
 657 //________________________________________________________________________
 658 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 659 {
 660   if(fDebug)
 661     printf("Inside FillClusHists()....\n");
 662   // Fill cluster histograms.
 663   TObjArray *clusters = fESDClusters;
 664
 665   if (!clusters){
 666     clusters = fAODClusters;
 667     if(fDebug)
 668       printf("ESD clusters empty...");
 669   }
 670   if (!clusters){
 671     if(fDebug)
 672       printf("  and AOD clusters as well!!!\n");
 673     return;
 674   }
 675   if(fDebug)
 676     printf("\n");
 677
 678   const Int_t nclus = clusters->GetEntries();
 679   if(nclus==0)
 680     return;
 681   if(fDebug)
 682     printf("Inside FillClusHists and there are %d clusters\n",nclus);
 683   Double_t maxE = 0;
 684   Int_t nclus10 = 0;
 685   Double_t ptmc=-1;
 686   for(Int_t ic=0;ic<nclus;ic++){
 687     maxE=0;
 688     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 689     if(!c)
 690       continue;
 691     if(!c->IsEMCAL())
 692       continue;
 693     if(c->E()<fECut)
 694       continue;
 695     if(fCpvFromTrack && fClusIdFromTracks.Contains(Form("%d",ic)))
 696        continue;
 697     if(IsExotic(c))
 698       continue;
 699     Short_t id;
 700     Double_t Emax = GetMaxCellEnergy( c, id);
 701     if(fDebug)
 702       printf("cluster max cell E=%1.1f",Emax);
 703     Float_t clsPos[3] = {0,0,0};
 704     c->GetPosition(clsPos);
 705     TVector3 clsVec(clsPos);
 706     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 707     if(fDebug)
 708       printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
 709     if(Et>10)
 710       nclus10++;
 711     Float_t ceiso, cephiband, cecore;
 712     Float_t triso, trphiband, trcore;
 713     Float_t alliso, allphiband;//, allcore;
 714     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 715     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 716     GetCeIso(clsVec, id, ceiso, cephiband, cecore);
 717     GetTrIso(clsVec, triso, trphiband, trcore);
 718     Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
 719     if(Et>10 && Et<15 && dr>0.025){
 720       fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
 721       if(fDebug)
 722         printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
 723     }
 724     alliso = ceiso + triso;
 725     allphiband = cephiband + trphiband;
 726     //allcore = cecore + trcore;
 727     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 728     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 729     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 730     Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR);
 731     if(fDebug && Et>10)
 732       printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
 733     if(c->GetM02()>0.5 && c->GetM02()<2.0){
 734       fMcPtInConeBG->Fill(alliso-allisoue, mcptsum);
 735       fMcPtInConeBGnoUE->Fill(alliso, mcptsum);
 736       fMcPtInConeTrBGnoUE->Fill(triso, mcptsum);
 737     }
 738     if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){
 739       fMcPtInConeSBG->Fill(alliso-allisoue, mcptsum);
 740       fMcPtInConeSBGnoUE->Fill(alliso, mcptsum);
 741       fMcPtInConeTrSBGnoUE->Fill(triso, mcptsum);
 742       if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){
 743         fAllIsoEtMcGamma->Fill(Et, alliso-cecore-allisoue);
 744         fAllIsoNoUeEtMcGamma->Fill(Et, alliso-cecore);
 745       }
 746     }
 747     Bool_t isCPV = kFALSE;
 748     if(TMath::Abs(c->GetTrackDx())>0.03 || TMath::Abs(c->GetTrackDz())>0.02)
 749       isCPV = kTRUE;
 750     if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV)
 751       fClusEtCPVSBGISO->Fill(Et,alliso - trcore);
 752     if(c->GetM02()>0.5 && c->GetM02()<2.0 && isCPV)
 753       fClusEtCPVBGISO->Fill(Et,alliso - trcore);
 754     const Int_t ndims =   fNDimensions;
 755     Double_t outputValues[ndims];
 756     if(mcptsum<2)
 757       ptmc = GetClusSource(c);
 758     else
 759       ptmc = -0.1;
 760     outputValues[0] = Et;
 761     outputValues[1] = c->GetM02();
 762     outputValues[2] = ceiso/*cecore*/-ceisoue;
 763     outputValues[3] = triso-trisoue;
 764     outputValues[4] = alliso/*cecore*/-allisoue - trcore;
 765     outputValues[5] = ceiso;
 766     outputValues[6] = alliso - trcore;
 767     if(fDebug)
 768       printf("track-cluster dphi=%1.3f, deta=%1.3f\n",c->GetTrackDx(),c->GetTrackDz());
 769     if(TMath::Abs(c->GetTrackDx())<0.1)
 770       outputValues[7] = c->GetTrackDx();
 771     else
 772       outputValues[7] = 0.099*c->GetTrackDx()/TMath::Abs(c->GetTrackDx());
 773     if(TMath::Abs(c->GetTrackDz())<0.05)
 774       outputValues[8] = c->GetTrackDz();
 775     else
 776       outputValues[8] = 0.049*c->GetTrackDz()/TMath::Abs(c->GetTrackDz());
 777     outputValues[9] = clsVec.Eta();
 778     outputValues[10] = clsVec.Phi();
 779     if(fESDCells)
 780       outputValues[11] = fESDCells->GetCellTime(id);
 781     else if(fAODCells)
 782       outputValues[11] = fAODCells->GetCellTime(id);
 783     outputValues[12] = fTrackMult;
 784     outputValues[13] = ptmc;
 785     fHnOutput->Fill(outputValues);
 786     if(c->E()>maxE)
 787       maxE = c->E();
 788   }
 789   fNClusHighClusE->Fill(maxE,nclus);
 790   fMaxEClus = maxE;
 791   fNClusEt10->Fill(nclus10);
 792   fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
 793 }
 794
 795 //________________________________________________________________________
 796 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core)
 797 {
 798   // Get cell isolation.
 799   AliVCaloCells *cells = fESDCells;
 800   if (!cells){
 801     cells = fAODCells;
 802     if(fDebug)
 803       printf("ESD cells empty...");
 804   }
 805   if (!cells){
 806      if(fDebug)
 807       printf("  and AOD cells are empty  as well!!!\n");
 808     return;
 809   }
 810
 811   TObjArray *clusters = fESDClusters;
 812   if (!clusters)
 813     clusters = fAODClusters;
 814   if (!clusters)
 815     return;
 816
 817
 818   const Int_t nclus = clusters->GetEntries();
 819   //const Int_t ncells = cells->GetNumberOfCells();
 820   Float_t totiso=0;
 821   Float_t totband=0;
 822   Float_t totcore=0;
 823   Float_t etacl = vec.Eta();
 824   Float_t phicl = vec.Phi();
 825   if(phicl<0)
 826     phicl+=TMath::TwoPi();
 827   /*Int_t absid = -1;
 828   Float_t eta=-1, phi=-1;
 829   for(int icell=0;icell<ncells;icell++){
 830     absid = TMath::Abs(cells->GetCellNumber(icell));
 831     Float_t celltime = cells->GetCellTime(absid);
 832     //if(TMath::Abs(celltime)>2e-8 && (!fIsMc))
 833     if(TMath::Abs(celltime-maxtcl)>2e-8 )
 834       continue;
 835     if(!fGeom)
 836       return;
 837     fGeom->EtaPhiFromIndex(absid,eta,phi);
 838     Float_t dphi = TMath::Abs(phi-phicl);
 839     Float_t deta = TMath::Abs(eta-etacl);
 840     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/
 841   for(int ic=0;ic<nclus;ic++){
 842     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 843     if(!c)
 844       continue;
 845     if(!c->IsEMCAL())
 846       continue;
 847     if(c->E()<fMinIsoClusE)
 848       continue;
 849     Short_t id;
 850     GetMaxCellEnergy( c, id);
 851     Double_t maxct = cells->GetCellTime(id);
 852     if(TMath::Abs(maxct)>fClusTDiff/*2.5e-9*/ && (!fIsMc))
 853       continue;
 854     Float_t clsPos[3] = {0,0,0};
 855     c->GetPosition(clsPos);
 856     TVector3 cv(clsPos);
 857     Double_t Et = c->E()*TMath::Sin(cv.Theta());
 858     Float_t dphi = TMath::Abs(cv.Phi()-phicl);
 859     Float_t deta = TMath::Abs(cv.Eta()-etacl);
 860     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 861     if(R<0.007)
 862       continue;
 863     if(maxid==id)
 864       continue;
 865     Double_t matchedpt =  GetTrackMatchedPt(c->GetTrackMatchedIndex());
 866     if(matchedpt>0 && fRemMatchClus)
 867       continue;
 868     Double_t nEt = TMath::Max(Et-matchedpt, 0.0);
 869     if(nEt<0)
 870       printf("nEt=%1.1f\n",nEt);
 871     if(R<fIsoConeR){
 872       totiso += nEt;
 873       if(R<0.04)
 874         totcore += nEt;
 875     }
 876     else{
 877       if(dphi>fIsoConeR)
 878         continue;
 879       if(deta<fIsoConeR)
 880         continue;
 881       totband += nEt;
 882     }
 883   }
 884   iso = totiso;
 885   phiband = totband;
 886   core = totcore;
 887 }
 888 //________________________________________________________________________
 889 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 890 {
 891   // Get track isolation.
 892
 893   if(!fSelPrimTracks)
 894     return;
 895   fHigherPtCone = 0;
 896   const Int_t ntracks = fSelPrimTracks->GetEntries();
 897   Double_t totiso=0;
 898   Double_t totband=0;
 899   Double_t totcore=0;
 900   Double_t etacl = vec.Eta();
 901   Double_t phicl = vec.Phi();
 902   if(phicl<0)
 903     phicl+=TMath::TwoPi();
 904   for(int itrack=0;itrack<ntracks;itrack++){
 905     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
 906     if(!track)
 907       continue;
 908     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 909     Double_t deta = TMath::Abs(track->Eta()-etacl);
 910     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 911     Double_t pt = track->Pt();
 912     if(pt>fHigherPtCone)
 913       fHigherPtCone = pt;
 914     if(R<fIsoConeR){
 915       totiso += track->Pt();
 916       if(R<0.04 && this->fTrCoreRem)
 917         totcore += pt;
 918     }
 919     else{
 920       if(dphi>fIsoConeR)
 921         continue;
 922       if(deta<fIsoConeR)
 923         continue;
 924       totband += track->Pt();
 925     }
 926   }
 927   iso = totiso;
 928   phiband = totband;
 929   core = totcore;
 930 }
 931
 932 //________________________________________________________________________
 933 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 934 {
 935   // Calculate the energy of cross cells around the leading cell.
 936
 937   AliVCaloCells *cells = 0;
 938   cells = fESDCells;
 939   if (!cells)
 940     cells = fAODCells;
 941   if (!cells)
 942     return 0;
 943
 944   if (!fGeom)
 945     return 0;
 946
 947   Int_t iSupMod = -1;
 948   Int_t iTower  = -1;
 949   Int_t iIphi   = -1;
 950   Int_t iIeta   = -1;
 951   Int_t iphi    = -1;
 952   Int_t ieta    = -1;
 953   Int_t iphis   = -1;
 954   Int_t ietas   = -1;
 955
 956   Double_t crossEnergy = 0;
 957
 958   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 959   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 960
 961   Int_t ncells = cluster->GetNCells();
 962   for (Int_t i=0; i<ncells; i++) {
 963     Int_t cellAbsId = cluster->GetCellAbsId(i);
 964     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 965     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 966     Int_t aphidiff = TMath::Abs(iphi-iphis);
 967     if (aphidiff>1)
 968       continue;
 969     Int_t aetadiff = TMath::Abs(ieta-ietas);
 970     if (aetadiff>1)
 971       continue;
 972     if ( (aphidiff==1 && aetadiff==0) ||
 973         (aphidiff==0 && aetadiff==1) ) {
 974       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 975     }
 976   }
 977
 978   return crossEnergy;
 979 }
 980
 981 //________________________________________________________________________
 982 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 983 {
 984   // Get maximum energy of attached cell.
 985
 986   id = -1;
 987
 988   AliVCaloCells *cells = 0;
 989   cells = fESDCells;
 990   if (!cells)
 991     cells = fAODCells;
 992   if(!cells)
 993     return 0;
 994
 995   Double_t maxe = 0;
 996   Int_t ncells = cluster->GetNCells();
 997   for (Int_t i=0; i<ncells; i++) {
 998     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 999     if (e>maxe) {
1000       maxe = e;
1001       id   = cluster->GetCellAbsId(i);
1002     }
1003   }
1004   return maxe;
1005 }
1006
1007 //________________________________________________________________________
1008 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
1009 {
1010   if(!fStack && !fAODMCParticles)
1011     return;
1012   //ESD
1013   if(fStack){
1014     Int_t nTracks = fStack->GetNtrack();
1015     if(fDebug)
1016       printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1017     for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1018       TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1019       if(!mcp)
1020         continue;
1021       Int_t pdg = mcp->GetPdgCode();
1022       if(pdg!=22)
1023         continue;
1024       if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1025         continue;
1026       Int_t imom = mcp->GetMother(0);
1027       if(imom<0 || imom>nTracks)
1028         continue;
1029       TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
1030       if(!mcmom)
1031         continue;
1032       Int_t pdgMom = mcmom->GetPdgCode();
1033       Double_t mcphi = mcp->Phi();
1034       Double_t mceta = mcp->Eta();
1035       if((imom==6 || imom==7) && pdgMom==22) {
1036         fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1037         Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1038         fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1039         if(ptsum<2)
1040           fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1041         if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1042           fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1043         if(fNClusForDirPho==0)
1044           fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1045         if(fDebug){
1046           printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatusCode());
1047           printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatusCode());
1048         }
1049       }
1050       else{
1051         if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1052           fDecayPhotonPtMC->Fill(mcp->Pt());
1053       }
1054     }
1055   }
1056   //AOD
1057   else if(fAODMCParticles){
1058     Int_t nTracks = fAODMCParticles->GetEntriesFast();
1059     if(fDebug)
1060       printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1061     for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1062       AliAODMCParticle *mcp = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1063       if(!mcp)
1064         continue;
1065       Int_t pdg = mcp->GetPdgCode();
1066       if(pdg!=22)
1067         continue;
1068       if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1069         continue;
1070       Int_t imom = mcp->GetMother();
1071       if(imom<0 || imom>nTracks)
1072         continue;
1073       AliAODMCParticle *mcmom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(imom));
1074       if(!mcmom)
1075         continue;
1076       Int_t pdgMom = mcmom->GetPdgCode();
1077       Double_t mcphi = mcp->Phi();
1078       Double_t mceta = mcp->Eta();
1079       if((imom==6 || imom==7) && pdgMom==22) {
1080         fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1081         Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1082         fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1083         if(ptsum<2)
1084           fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1085         if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1086           fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1087         if(fNClusForDirPho==0)
1088           fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1089         if(fDebug){
1090           printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatus());
1091           printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatus());
1092         }
1093       }
1094       else{
1095         if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1096           fDecayPhotonPtMC->Fill(mcp->Pt());
1097       }
1098     }
1099   }
1100 }
1101 //________________________________________________________________________
1102 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
1103 {
1104   if(!c)
1105     return -0.1;
1106   if(!fStack && !fAODMCParticles)
1107     return -0.1;
1108   Int_t clabel = c->GetLabel();
1109   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
1110     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
1111   if(!fMcIdFamily.Contains(Form("%d",clabel)))
1112     return -0.1;
1113   fNClusForDirPho++;
1114   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
1115   Int_t partonpos = partonposstr.Atoi();
1116   if(fDebug)
1117     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
1118   Float_t clsPos[3] = {0,0,0};
1119   c->GetPosition(clsPos);
1120   TVector3 clsVec(clsPos);
1121   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
1122   //ESD
1123   if(fStack){
1124     Int_t nmcp = fStack->GetNtrack();
1125     if(clabel<0 || clabel>nmcp)
1126       return -0.1;
1127     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
1128     if(!mcp)
1129       return -0.1;
1130     if(fDebug){
1131       printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),clabel);
1132     }
1133     Int_t lab1 =  mcp->GetFirstDaughter();
1134     if(lab1<0 || lab1>nmcp)
1135       return -0.1;
1136     TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
1137     if(!mcd)
1138       return -0.1;
1139     if(fDebug)
1140       printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->Energy(),mcd->GetPdgCode(),lab1);
1141     if(mcd->GetPdgCode()==22){
1142       fClusEtMcPt->Fill(mcd->Pt(), Et);
1143       fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1144     }
1145     else{
1146       printf("Warning: daughter of photon parton is not a photon\n");
1147       return -0.1;
1148     }
1149   }
1150   //AOD
1151   else if(fAODMCParticles){
1152     Int_t nmcp = fAODMCParticles->GetEntriesFast();
1153     if(clabel<0 || clabel>nmcp)
1154       return -0.1;
1155     AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(partonpos));
1156     if(!mcp)
1157       return -0.1;
1158     if(fDebug){
1159       printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),clabel);
1160     }
1161     Int_t lab1 =  mcp->GetDaughter(0);
1162     if(lab1<0 || lab1>nmcp)
1163       return -0.1;
1164     AliAODMCParticle *mcd = static_cast<AliAODMCParticle*>(fAODMCParticles->At(lab1));
1165     if(!mcd)
1166       return -0.1;
1167     if(fDebug)
1168       printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->E(),mcd->GetPdgCode(),lab1);
1169     if(mcd->GetPdgCode()==22){
1170       fClusEtMcPt->Fill(mcd->Pt(), Et);
1171       fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1172     }
1173     else{
1174       printf("Warning: daughter of photon parton is not a photon\n");
1175       return -0.1;
1176     }
1177   }
1178   return fDirPhoPt;
1179 }
1180 //________________________________________________________________________
1181 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
1182 {
1183   if(!fStack && !fAODMCParticles)
1184     return;
1185   Int_t selfid = 6;
1186   //ESD
1187   if(fStack){
1188     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
1189     if(!mcp)
1190       return;
1191     if(mcp->GetPdgCode()!=22){
1192       mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
1193       if(!mcp)
1194         return;
1195     }
1196     Int_t daug0f =  mcp->GetFirstDaughter();
1197     Int_t daug0l =  mcp->GetLastDaughter();
1198     Int_t nd0 = daug0l - daug0f;
1199     if(fDebug)
1200       printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),nd0+1);
1201     fMcIdFamily = Form("%d,",selfid);
1202     GetDaughtersInfo(daug0f,daug0l, selfid,"");
1203     if(fDebug){
1204       printf("\t---- end of the prompt  gamma's family tree ----\n\n");
1205       printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1206     }
1207     TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
1208     if(!md)
1209       return;
1210     fDirPhoPt = md->Pt();
1211   }
1212   //AOD
1213   else   if(fAODMCParticles){
1214     AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(selfid));
1215     if(!mcp)
1216       return;
1217     if(mcp->GetPdgCode()!=22){
1218       mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(++selfid));
1219       if(!mcp)
1220         return;
1221     }
1222     Int_t daug0f =  mcp->GetDaughter(0);
1223     Int_t daug0l =  mcp->GetDaughter(mcp->GetNDaughters()-1);
1224     Int_t nd0 = daug0l - daug0f;
1225     if(fDebug)
1226       printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),nd0+1);
1227     fMcIdFamily = Form("%d,",selfid);
1228     GetDaughtersInfo(daug0f,daug0l, selfid,"");
1229     if(fDebug){
1230       printf("\t---- end of the prompt  gamma's family tree ----\n\n");
1231       printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1232     }
1233     AliAODMCParticle *md = static_cast<AliAODMCParticle*>(fAODMCParticles->At(daug0f));
1234     if(!md)
1235       return;
1236     fDirPhoPt = md->Pt();
1237   }
1238
1239 }
1240 //________________________________________________________________________
1241 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
1242 {
1243   if(fStack){
1244     int nmcp = fStack->GetNtrack();
1245     if(firstd<0 || firstd>nmcp)
1246       return;
1247     if(lastd<0 || lastd>nmcp)
1248       return;
1249     if(firstd>lastd){
1250       printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1251     return;
1252     }
1253     TString indenter = Form("\t%s",inputind);
1254     TParticle *dp = 0x0;
1255     if(fDebug)
1256       printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1257     for(int id=firstd; id<lastd+1; id++){
1258       dp =  static_cast<TParticle*>(fStack->Particle(id));
1259       if(!dp)
1260         continue;
1261       Int_t dfd = dp->GetFirstDaughter();
1262       Int_t dld = dp->GetLastDaughter();
1263       Int_t dnd =  dld - dfd + 1;
1264       Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
1265       if(fDebug)
1266         printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->Energy(), vr, dp->GetPdgCode(), dnd, dfd, dld);
1267       fMcIdFamily += Form("%d,",id);
1268       GetDaughtersInfo(dfd,dld,id,indenter.Data());
1269     }
1270   }
1271   if(fAODMCParticles){
1272     int nmcp = fAODMCParticles->GetEntriesFast();
1273     if(firstd<0 || firstd>nmcp)
1274       return;
1275     if(lastd<0 || lastd>nmcp)
1276       return;
1277     if(firstd>lastd){
1278       printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1279     return;
1280     }
1281     TString indenter = Form("\t%s",inputind);
1282     AliAODMCParticle *dp = 0x0;
1283     if(fDebug)
1284       printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1285     for(int id=firstd; id<lastd+1; id++){
1286       dp =  static_cast<AliAODMCParticle*>(fAODMCParticles->At(id));
1287       if(!dp)
1288         continue;
1289       Int_t dfd = dp->GetDaughter(0);
1290       Int_t dld = dp->GetDaughter(dp->GetNDaughters()-1);
1291       Int_t dnd =  dld - dfd + 1;
1292       Float_t vr = TMath::Sqrt(dp->Xv()*dp->Xv()+dp->Xv()*dp->Xv());
1293       if(fDebug)
1294         printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->E(), vr, dp->GetPdgCode(), dnd, dfd, dld);
1295       fMcIdFamily += Form("%d,",id);
1296       GetDaughtersInfo(dfd,dld,id,indenter.Data());
1297     }
1298   }
1299 }
1300
1301 //________________________________________________________________________
1302 Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
1303 {
1304   if(!fStack && !fAODMCParticles)
1305     return 0;
1306   if(fDebug)
1307     printf("Inside GetMcPtSumInCone!!\n");
1308   Float_t ptsum = 0;
1309   TString addpartlabels = "";
1310   //ESD
1311   if(fStack){
1312     Int_t nTracks = fStack->GetNtrack();
1313     for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1314       TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1315       if(!mcp)
1316         continue;
1317       Int_t status = mcp->GetStatusCode();
1318       if(status!=1)
1319         continue;
1320       Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz());
1321       if(mcvr>1)
1322         continue;
1323       /*else {
1324         if(fDebug)
1325         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1326         }*/
1327       Float_t dphi = mcp->Phi() - phiclus;
1328       Float_t deta = mcp->Eta() - etaclus;
1329       if(fDebug && TMath::Abs(dphi)<0.01)
1330         printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1331
1332       if(deta>R || dphi>R)
1333         continue;
1334       Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1335       if(dR>R)
1336         continue;
1337       addpartlabels += Form("%d,",iTrack);
1338       if(mcp->Pt()<0.2)
1339         continue;
1340       ptsum += mcp->Pt();
1341     }
1342   }
1343   //AOD
1344   if(fAODMCParticles){
1345     Int_t nTracks = fAODMCParticles->GetEntriesFast();
1346     for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1347       AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1348       if(!mcp)
1349         continue;
1350       Int_t status = mcp->GetStatus();
1351       if(status!=1)
1352         continue;
1353       Float_t mcvr = TMath::Sqrt(mcp->Xv()*mcp->Xv()+ mcp->Yv()* mcp->Yv() + mcp->Zv()*mcp->Zv());
1354       if(mcvr>1)
1355         continue;
1356       /*else {
1357         if(fDebug)
1358         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1359         }*/
1360       Float_t dphi = mcp->Phi() - phiclus;
1361       Float_t deta = mcp->Eta() - etaclus;
1362       if(fDebug && TMath::Abs(dphi)<0.01)
1363         printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1364
1365       if(deta>R || dphi>R)
1366         continue;
1367       Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1368       if(dR>R)
1369         continue;
1370       addpartlabels += Form("%d,",iTrack);
1371       if(mcp->Pt()<0.2)
1372         continue;
1373       ptsum += mcp->Pt();
1374     }
1375   }
1376   return ptsum;
1377 }
1378 //________________________________________________________________________
1379 void AliAnalysisTaskEMCALIsoPhoton::FillQA()
1380 {
1381
1382   TObjArray *clusters = fESDClusters;
1383   //"none", "exotic", "exo+cpv1", "exo+cpv1+time", "exo+cpv1+time+m02"),
1384   if (!clusters){
1385     clusters = fAODClusters;
1386     if(fDebug)
1387       printf("ESD clusters empty...");
1388   }
1389   if (!clusters){
1390     if(fDebug)
1391       printf("  and AOD clusters as well!!!\n");
1392     return;
1393   }
1394   if(!fSelPrimTracks)
1395     return;
1396   const int ntracks = fSelPrimTracks->GetEntriesFast();
1397   const int ncells = fNCells50;//fESDCells->GetNumberOfCells();
1398   const Int_t nclus = clusters->GetEntries();
1399   fNTracks->Fill(ntracks);
1400   fEmcNCells->Fill(ncells);
1401   fEmcNClus->Fill(nclus);
1402   if(fMaxEClus>fECut){
1403     fNTracksECut->Fill(ntracks);
1404     fEmcNCellsCut->Fill(ncells);
1405     fEmcNClusCut->Fill(nclus);
1406   }
1407   for(int it=0;it<ntracks;it++){
1408     AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it);
1409     if(!t)
1410       continue;
1411     fTrackPtPhi->Fill(t->Pt(),t->Phi());
1412     fTrackPtEta->Fill(t->Pt(),t->Eta());
1413     if(fMaxEClus>fECut){
1414       fTrackPtPhiCut->Fill(t->Pt(), t->Phi());
1415       fTrackPtEtaCut->Fill(t->Pt(), t->Eta());
1416     }
1417   }
1418   for(int ic=0;ic<nclus;ic++){
1419     AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(ic));
1420     //AliESDCaloCluster *c = (AliESDCaloCluster*)clusters->At(ic);
1421     if(!c)
1422       continue;
1423     if(!c->IsEMCAL())
1424       continue;
1425     Float_t clsPos[3] = {0,0,0};
1426     c->GetPosition(clsPos);
1427     TVector3 clsVec(clsPos);
1428     Double_t cphi = clsVec.Phi();
1429     Double_t ceta = clsVec.Eta();
1430     Short_t id;
1431     GetMaxCellEnergy( c, id);
1432     fEmcClusEClusCuts->Fill(c->E(),0);
1433     fEmcClusEPhi->Fill(c->E(), cphi);
1434     fEmcClusEEta->Fill(c->E(), ceta);
1435     if(fMaxEClus>fECut){
1436       fEmcClusEPhiCut->Fill(c->E(), cphi);
1437       fEmcClusEEtaCut->Fill(c->E(), ceta);
1438     }
1439     Double_t maxt=0;
1440     if(fESDCells)
1441       maxt = fESDCells->GetCellTime(id);
1442     else if(fAODCells)
1443       maxt = fAODCells->GetCellTime(id);
1444     if(IsExotic(c))
1445       continue;
1446     fEmcClusNotExo->Fill(c->E());
1447     fEmcClusEClusCuts->Fill(c->E(),1);
1448     if(fClusIdFromTracks.Contains(Form("%d",ic)))
1449       fEmcClusETM2->Fill(c->E());
1450     if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
1451       fEmcClusETM1->Fill(c->E());
1452       continue;
1453     }
1454     fEmcClusEClusCuts->Fill(c->E(),2);
1455     if(TMath::Abs(maxt)>30e-9)
1456       continue;
1457     fEmcClusEClusCuts->Fill(c->E(),3);
1458     if(c->GetM02()>0.1)
1459       fEmcClusEClusCuts->Fill(c->E(),4);
1460   }
1461 }
1462 //________________________________________________________________________
1463 Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex)
1464 {
1465   Double_t pt = 0;
1466   if(!fTracks)
1467     return pt;
1468   if(matchIndex<0 || matchIndex>fTracks->GetEntries()){
1469     if(fDebug)
1470       printf("track-matched index out of track array range!!!\n");
1471     return pt;
1472   }
1473   AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex));
1474   if(!track){
1475     if(fDebug)
1476       printf("track-matched pointer does not exist!!!\n");
1477     return pt;
1478   }
1479   if(fESD){
1480     if(!fPrTrCuts)
1481       return pt;
1482     if(!fPrTrCuts->IsSelected(track))
1483       return pt;
1484     pt = track->Pt();
1485   }
1486   return pt;
1487 }
1488 //________________________________________________________________________
1489 void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells()
1490 {
1491   AliVCaloCells *cells = 0;
1492   cells = fESDCells;
1493   if (!cells)
1494     cells = fAODCells;
1495   if(!cells)
1496     return;
1497   Double_t maxe = 0;
1498   Double_t maxphi = -10;
1499   Int_t ncells = cells->GetNumberOfCells();
1500   Double_t eta,phi;
1501   for (Int_t i=0; i<ncells; i++) {
1502     Short_t absid = TMath::Abs(cells->GetCellNumber(i));
1503     Double_t e = cells->GetCellAmplitude(absid);
1504     if(e>0.05)
1505       fNCells50++;
1506     else
1507       continue;
1508     fGeom->EtaPhiFromIndex(absid,eta,phi);
1509     if(maxe<e){
1510       maxe = e;
1511       maxphi = phi;
1512     }
1513   }
1514   fMaxCellEPhi->Fill(maxe,maxphi);
1515
1516 }
1517 //________________________________________________________________________
1518 bool AliAnalysisTaskEMCALIsoPhoton::IsExotic(AliVCluster *c)
1519 {
1520   bool isExo = 0;
1521   Short_t id;
1522   Double_t Emax = GetMaxCellEnergy( c, id);
1523   Double_t Ecross = GetCrossEnergy( c, id);
1524   if((1-Ecross/Emax)>fExoticCut)
1525     isExo = 1;
1526   return isExo;
1527 }
1528 //________________________________________________________________________
1529 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
1530 {
1531   // Called once at the end of the query.
1532 }