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