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