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