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