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     for(int nk=0;nk<l->GetEntries();nk++){
 658       TObject *obj = (TObject*)l->At(nk);
 659       TString oname = obj->GetName();
 660       if(oname.Contains("CaloClus"))
 661         clusArrayName = oname;
 662     }
 663     fESDClusters =  dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName));
 664     fESDCells = fESD->GetEMCALCells();
 665     if(fDebug)
 666       printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast());
 667   }
 668   else if(fAOD){
 669     l = fAOD->GetList();
 670     //l->Print();
 671     //fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
 672     for(int nk=0;nk<l->GetEntries();nk++){
 673       TObject *obj = (TObject*)l->At(nk);
 674       TString oname = obj->GetName();
 675       if(oname.Contains("aloClus"))
 676         clusArrayName = oname;
 677     }
 678     fAODClusters = dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName));
 679     fAODCells = fAOD->GetEMCALCells();
 680     if(fDebug)
 681       printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast());
 682   }
 683   if(clusArrayName=="CaloClusters")
 684     fClusArrayNames->Fill(0);
 685   else{
 686     if(clusArrayName=="EmcCaloClusters")
 687       fClusArrayNames->Fill(1);
 688     else
 689       fClusArrayNames->Fill(2);
 690   }
 691   if(fDebug){
 692         printf("clus array is named %s +++++++++\n",clusArrayName.Data());
 693   }
 694
 695
 696   fMCEvent = MCEvent();
 697   if(fMCEvent){
 698     if(fDebug)
 699       std::cout<<"MCevent exists\n";
 700     fStack = (AliStack*)fMCEvent->Stack();
 701     if(!fStack)
 702       fAODMCParticles = (TClonesArray*)fVEvent->FindListObject("mcparticles");
 703   }
 704   else{
 705     if(fDebug)
 706       std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
 707   }
 708   LoopOnCells();
 709   FollowGamma();
 710   if(fDebug)
 711     printf("passed calling of FollowGamma\n");
 712   FillClusHists();
 713   if(fDebug)
 714     printf("passed calling of FillClusHists\n");
 715   FillMcHists();
 716   if(fDebug)
 717     printf("passed calling of FillMcHists\n");
 718   if(fFillQA)
 719     FillQA();
 720   if(fDebug)
 721     printf("passed calling of FillQA\n");
 722   if(fESD)
 723     fESDClusters->Clear();
 724   fSelPrimTracks->Clear();
 725   fNClusForDirPho = 0;
 726   fNCells50 = 0;
 727   fClusIdFromTracks = "";
 728   fVecPv.Clear();
 729
 730   PostData(1, fOutputList);
 731   PostData(2, fQAList);
 732 }
 733
 734 //________________________________________________________________________
 735 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 736 {
 737   if(fDebug)
 738     printf("Inside FillClusHists()....\n");
 739   // Fill cluster histograms.
 740   TObjArray *clusters = fESDClusters;
 741
 742   if (!clusters){
 743     clusters = fAODClusters;
 744     if(fDebug)
 745       printf("ESD clusters empty...");
 746   }
 747   if (!clusters){
 748     if(fDebug)
 749       printf("  and AOD clusters as well!!!\n");
 750     return;
 751   }
 752   if(fDebug)
 753     printf("\n");
 754
 755   const Int_t nclus = clusters->GetEntries();
 756   if(nclus==0)
 757     return;
 758   if(fDebug)
 759     printf("Inside FillClusHists and there are %d clusters\n",nclus);
 760   Double_t maxE = 0;
 761   Int_t nclus10 = 0;
 762   Double_t ptmc=-1;
 763   for(Int_t ic=0;ic<nclus;ic++){
 764     maxE=0;
 765     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 766     if(!c){
 767       if(fDebug)
 768         printf("cluster pointer does not exist! xxxx\n");
 769       continue;
 770     }
 771     if(!c->IsEMCAL()){
 772       if(fDebug)
 773         printf("cluster is not EMCAL! xxxx\n");
 774       continue;
 775     }
 776     if(c->E()<fECut){
 777       if(fDebug)
 778         printf("cluster has E<%1.1f! xxxx\n", fECut);
 779       continue;
 780     }
 781     if(fCpvFromTrack && fClusIdFromTracks.Contains(Form("%d",ic))){
 782       if(fDebug)
 783         printf("cluster does not pass CPV criterion! xxxx\n");
 784        continue;
 785     }
 786     if(IsExotic(c)){
 787       if(fDebug)
 788         printf("cluster is exotic! xxxx\n");
 789       continue;
 790     }
 791     if(c->GetDistanceToBadChannel()<fDistToBadChan){
 792       if(fDebug)
 793         printf("cluster distance to bad channel is %1.1f (<%1.1f) xxxx\n",c->GetDistanceToBadChannel(),fDistToBadChan);
 794       continue;
 795     }
 796     Short_t id;
 797     Double_t Emax = GetMaxCellEnergy( c, id);
 798     if(fDebug)
 799       printf("cluster max cell E=%1.1f",Emax);
 800     Float_t clsPos[3] = {0,0,0};
 801     c->GetPosition(clsPos);
 802     TVector3 clsVec(clsPos);
 803     clsVec -= fVecPv;
 804     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 805     if(fDebug)
 806       printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
 807     if(Et>10)
 808       nclus10++;
 809     Float_t ceiso=0, cephiband=0, cecore=0;
 810     Float_t triso=0, trphiband=0, trcore=0;
 811     Float_t alliso=0, allphiband=0;//, allcore;
 812     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 813     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 814     Bool_t isCPV = kFALSE;
 815     if(TMath::Abs(c->GetTrackDx())>0.03 || TMath::Abs(c->GetTrackDz())>0.02)
 816       isCPV = kTRUE;
 817     GetCeIso(clsVec, id, ceiso, cephiband, cecore, Et);
 818     GetTrIso(clsVec, triso, trphiband, trcore);
 819     Int_t nInConePairs = 0;
 820     Double_t onePairMass = 0;
 821     if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV){
 822       TObjArray *inConeInvMassArr = (TObjArray*)fInConeInvMass.Tokenize(";");
 823       TObjArray *inConePairClEt =  (TObjArray*)fInConePairClEt.Tokenize(";");
 824       nInConePairs = inConeInvMassArr->GetEntriesFast();
 825       Int_t nInConePi0 = inConePairClEt->GetEntriesFast();
 826       if(nInConePairs != nInConePi0)
 827         printf("Inconsistent number of in cone pairs!!!\n");
 828       for(int ipair=0;ipair<nInConePairs;ipair++){
 829         TObjString *obs = (TObjString*)inConeInvMassArr->At(ipair);
 830         TObjString *obet = (TObjString*)inConePairClEt->At(ipair);
 831         TString smass = obs->GetString();
 832         TString spairEt = obet->GetString();
 833         Double_t pairmass = smass.Atof();
 834         Double_t pairEt = spairEt.Atof();//this must be zero when inv mass outside pi0 range
 835         if(0==ipair && nInConePairs==1)
 836           onePairMass = pairmass;
 837         if(fDebug)
 838           printf("=================+++++++++++++++Inv mass inside the cone for photon range: %1.1f,%1.1f,%1.1f+-++++-+-+-+-++-+-+-\n",Et,pairmass,ceiso+triso);
 839         fEtCandIsoAndIsoWoPairEt->Fill(Et,ceiso+triso,ceiso+triso-pairEt);
 840       }
 841     }
 842     Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
 843     if(Et>10 && Et<15 && dr>0.025){
 844       fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
 845       if(fDebug)
 846         printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
 847     }
 848     alliso = ceiso + triso;
 849     allphiband = cephiband + trphiband;
 850     //allcore = cecore + trcore;
 851     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 852     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 853     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 854     Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR);
 855     if(fDebug && Et>10)
 856       printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
 857     if(c->GetM02()>0.5 && c->GetM02()<2.0){
 858       fMcPtInConeBG->Fill(alliso-allisoue, mcptsum);
 859       fMcPtInConeBGnoUE->Fill(alliso, mcptsum);
 860       fMcPtInConeTrBGnoUE->Fill(triso, mcptsum);
 861     }
 862     if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){
 863       fMcPtInConeSBG->Fill(alliso-allisoue, mcptsum);
 864       fMcPtInConeSBGnoUE->Fill(alliso, mcptsum);
 865       fMcPtInConeTrSBGnoUE->Fill(triso, mcptsum);
 866       if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){
 867         fAllIsoEtMcGamma->Fill(Et, alliso-cecore-allisoue);
 868         fAllIsoNoUeEtMcGamma->Fill(Et, alliso-cecore);
 869       }
 870     }
 871     if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV)
 872       fClusEtCPVSBGISO->Fill(Et,alliso - trcore);
 873     if(c->GetM02()>0.5 && c->GetM02()<2.0 && isCPV)
 874       fClusEtCPVBGISO->Fill(Et,alliso - trcore);
 875     const Int_t ndims =   fNDimensions;
 876     Double_t outputValues[ndims];
 877     if(mcptsum<2)
 878       ptmc = GetClusSource(c);
 879     else
 880       ptmc = -0.1;
 881     outputValues[0] = Et;
 882     outputValues[1] = c->GetM02();
 883     outputValues[2] = ceiso/*cecore*/-ceisoue;
 884     outputValues[3] = triso-trisoue;
 885     outputValues[4] = alliso/*cecore*/-allisoue - trcore;
 886     outputValues[5] = ceiso;
 887     outputValues[6] = alliso - trcore;
 888     if(fDebug)
 889       printf("track-cluster dphi=%1.3f, deta=%1.3f\n",c->GetTrackDx(),c->GetTrackDz());
 890     if(TMath::Abs(c->GetTrackDx())<0.1)
 891       outputValues[7] = c->GetTrackDx();
 892     else
 893       outputValues[7] = 0.099*c->GetTrackDx()/TMath::Abs(c->GetTrackDx());
 894     if(TMath::Abs(c->GetTrackDz())<0.05)
 895       outputValues[8] = c->GetTrackDz();
 896     else
 897       outputValues[8] = 0.049*c->GetTrackDz()/TMath::Abs(c->GetTrackDz());
 898     outputValues[9] = clsVec.Eta();
 899     outputValues[10] = clsVec.Phi();
 900     if(fESDCells)
 901       outputValues[11] = fESDCells->GetCellTime(id);
 902     else if(fAODCells)
 903       outputValues[11] = fAODCells->GetCellTime(id);
 904     outputValues[12] = fTrackMult;
 905     outputValues[13] = ptmc;
 906     if(nInConePairs == 1)
 907       outputValues[14] = onePairMass;
 908     else
 909       outputValues[14] = -1;
 910     fHnOutput->Fill(outputValues);
 911     if(c->E()>maxE)
 912       maxE = c->E();
 913   }
 914   fNClusHighClusE->Fill(maxE,nclus);
 915   fMaxEClus = maxE;
 916   fNClusEt10->Fill(nclus10);
 917   fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
 918 }
 919
 920 //________________________________________________________________________
 921 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core, Double_t EtCl)
 922 {
 923   if(fDebug)
 924     printf("....indside GetCeIso funtcion\n");
 925   // Get cell isolation.
 926   AliVCaloCells *cells = fESDCells;
 927   if (!cells){
 928     cells = fAODCells;
 929     if(fDebug)
 930       printf("ESD cells empty...\n");
 931   }
 932   if (!cells){
 933      if(fDebug)
 934       printf("  and AOD cells are empty  as well!!!\n");
 935     return;
 936   }
 937
 938   TObjArray *clusters = fESDClusters;
 939   if (!clusters)
 940     clusters = fAODClusters;
 941   if (!clusters)
 942     return;
 943
 944   fInConeInvMass = "";
 945   fInConePairClEt="";
 946   const Int_t nclus = clusters->GetEntries();
 947   //const Int_t ncells = cells->GetNumberOfCells();
 948   Float_t totiso=0;
 949   Float_t totband=0;
 950   Float_t totcore=0;
 951   Float_t etacl = vec.Eta();
 952   Float_t phicl = vec.Phi();
 953   if(phicl<0)
 954     phicl+=TMath::TwoPi();
 955   /*Int_t absid = -1;
 956   Float_t eta=-1, phi=-1;
 957   for(int icell=0;icell<ncells;icell++){
 958     absid = TMath::Abs(cells->GetCellNumber(icell));
 959     Float_t celltime = cells->GetCellTime(absid);
 960     //if(TMath::Abs(celltime)>2e-8 && (!fIsMc))
 961     if(TMath::Abs(celltime-maxtcl)>2e-8 )
 962       continue;
 963     if(!fGeom)
 964       return;
 965     fGeom->EtaPhiFromIndex(absid,eta,phi);
 966     Float_t dphi = TMath::Abs(phi-phicl);
 967     Float_t deta = TMath::Abs(eta-etacl);
 968     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/
 969   for(int ic=0;ic<nclus;ic++){
 970     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 971     if(!c)
 972       continue;
 973     if(!c->IsEMCAL())
 974       continue;
 975     if(c->E()<fMinIsoClusE)
 976       continue;
 977     Short_t id=-1;
 978     GetMaxCellEnergy( c, id);
 979     Double_t maxct = cells->GetCellTime(id);
 980     if(TMath::Abs(maxct)>fClusTDiff/*2.5e-9*/ && (!fIsMc))
 981       continue;
 982     Float_t clsPos[3] = {0,0,0};
 983     c->GetPosition(clsPos);
 984     TVector3 cv(clsPos);
 985     cv -= fVecPv;
 986     Double_t Et = c->E()*TMath::Sin(cv.Theta());
 987     Float_t dphi = (cv.Phi()-phicl);
 988     Float_t deta = (cv.Eta()-etacl);
 989     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 990     if(R<0.007)
 991       continue;
 992     if(maxid==id)
 993       continue;
 994     Double_t matchedpt =  GetTrackMatchedPt(c->GetTrackMatchedIndex());
 995     if(fCpvFromTrack){
 996       if(matchedpt>0 && fRemMatchClus )
 997         continue;
 998     } else {
 999       if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
1000         if(fRemMatchClus){
1001           if(fDebug)
1002             printf("This isolation cluster is matched to a track!++++++++++++++++++++++++++++++++++++++++++++++++++\n");
1003           continue;
1004         }
1005       }
1006     }
1007     Double_t nEt = TMath::Max(Et-matchedpt, 0.0);
1008     if(nEt<0)
1009       printf("nEt=%1.1f\n",nEt);
1010     if(R<fIsoConeR){
1011       if(c->GetM02()>0.1 && c->GetM02()<0.3 && !(matchedpt>0)){
1012         TLorentzVector lv, lvec;
1013         lv.SetPtEtaPhiM(Et,cv.Eta(),cv.Phi(),0);
1014         lvec.SetPtEtaPhiM(EtCl,vec.Eta(),vec.Phi(),0);
1015         TLorentzVector lpair = lv + lvec;
1016         fInConeInvMass += Form("%f;",lpair.M());
1017         if(lpair.M()>0.11 && lpair.M()<0.165){
1018           fInConePairedClusEtVsCandEt->Fill(EtCl,Et);
1019           fInConePairClEt += Form("%f;",Et);
1020           continue;
1021         }
1022         else
1023           fInConePairClEt += Form("%f;",0.0);
1024         if(lpair.M()>0.52 && lpair.M()<0.58)
1025           continue;
1026       }
1027       totiso += nEt;
1028       if(R<0.04)
1029         totcore += nEt;
1030     }
1031     else{
1032       if(dphi>fIsoConeR)
1033         continue;
1034       if(deta<fIsoConeR)
1035         continue;
1036       totband += nEt;
1037     }
1038   }
1039   iso = totiso;
1040   phiband = totband;
1041   core = totcore;
1042 }
1043 //________________________________________________________________________
1044 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
1045 {
1046   // Get track isolation.
1047
1048   if(!fSelPrimTracks)
1049     return;
1050   fHigherPtCone = 0;
1051   const Int_t ntracks = fSelPrimTracks->GetEntries();
1052   Double_t totiso=0;
1053   Double_t totband=0;
1054   Double_t totcore=0;
1055   Double_t etacl = vec.Eta();
1056   Double_t phicl = vec.Phi();
1057   if(phicl<0)
1058     phicl+=TMath::TwoPi();
1059   for(int itrack=0;itrack<ntracks;itrack++){
1060     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
1061     if(!track)
1062       continue;
1063     Double_t dphi = TMath::Abs(track->Phi()-phicl);
1064     Double_t deta = TMath::Abs(track->Eta()-etacl);
1065     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
1066     Double_t pt = track->Pt();
1067     if(pt>fHigherPtCone)
1068       fHigherPtCone = pt;
1069     if(R<fIsoConeR){
1070       totiso += track->Pt();
1071       if(R<0.04 && this->fTrCoreRem)
1072         totcore += pt;
1073     }
1074     else{
1075       if(dphi>fIsoConeR)
1076         continue;
1077       if(deta<fIsoConeR)
1078         continue;
1079       totband += track->Pt();
1080     }
1081   }
1082   iso = totiso;
1083   phiband = totband;
1084   core = totcore;
1085 }
1086
1087 //________________________________________________________________________
1088 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1089 {
1090   // Calculate the energy of cross cells around the leading cell.
1091
1092   AliVCaloCells *cells = 0;
1093   cells = fESDCells;
1094   if (!cells)
1095     cells = fAODCells;
1096   if (!cells)
1097     return 0;
1098
1099   if (!fGeom)
1100     return 0;
1101
1102   Int_t iSupMod = -1;
1103   Int_t iTower  = -1;
1104   Int_t iIphi   = -1;
1105   Int_t iIeta   = -1;
1106   Int_t iphi    = -1;
1107   Int_t ieta    = -1;
1108   Int_t iphis   = -1;
1109   Int_t ietas   = -1;
1110
1111   Double_t crossEnergy = 0;
1112
1113   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1114   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1115
1116   Int_t ncells = cluster->GetNCells();
1117   for (Int_t i=0; i<ncells; i++) {
1118     Int_t cellAbsId = cluster->GetCellAbsId(i);
1119     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1120     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1121     Int_t aphidiff = TMath::Abs(iphi-iphis);
1122     if (aphidiff>1)
1123       continue;
1124     Int_t aetadiff = TMath::Abs(ieta-ietas);
1125     if (aetadiff>1)
1126       continue;
1127     if ( (aphidiff==1 && aetadiff==0) ||
1128         (aphidiff==0 && aetadiff==1) ) {
1129       crossEnergy += cells->GetCellAmplitude(cellAbsId);
1130     }
1131   }
1132
1133   return crossEnergy;
1134 }
1135
1136 //________________________________________________________________________
1137 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1138 {
1139   // Get maximum energy of attached cell.
1140
1141   id = -1;
1142
1143   AliVCaloCells *cells = 0;
1144   cells = fESDCells;
1145   if (!cells)
1146     cells = fAODCells;
1147   if(!cells)
1148     return 0;
1149
1150   Double_t maxe = 0;
1151   Int_t ncells = cluster->GetNCells();
1152   for (Int_t i=0; i<ncells; i++) {
1153     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1154     if (e>maxe) {
1155       maxe = e;
1156       id   = cluster->GetCellAbsId(i);
1157     }
1158   }
1159   return maxe;
1160 }
1161
1162 //________________________________________________________________________
1163 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
1164 {
1165   if(!fStack && !fAODMCParticles)
1166     return;
1167   //ESD
1168   if(fStack){
1169     Int_t nTracks = fStack->GetNtrack();
1170     if(fDebug)
1171       printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1172     for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1173       TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1174       if(!mcp)
1175         continue;
1176       Int_t pdg = mcp->GetPdgCode();
1177       if(pdg!=22)
1178         continue;
1179       if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1180         continue;
1181       Int_t imom = mcp->GetMother(0);
1182       if(imom<0 || imom>nTracks)
1183         continue;
1184       TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
1185       if(!mcmom)
1186         continue;
1187       Int_t pdgMom = mcmom->GetPdgCode();
1188       Double_t mcphi = mcp->Phi();
1189       Double_t mceta = mcp->Eta();
1190       if((imom==6 || imom==7) && pdgMom==22) {
1191         fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1192         Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1193         fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1194         if(ptsum<2)
1195           fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1196         if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1197           fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1198         if(fNClusForDirPho==0)
1199           fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1200         if(fDebug){
1201           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());
1202           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());
1203         }
1204       }
1205       else{
1206         if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1207           fDecayPhotonPtMC->Fill(mcp->Pt());
1208       }
1209     }
1210   }
1211   //AOD
1212   else if(fAODMCParticles){
1213     Int_t nTracks = fAODMCParticles->GetEntriesFast();
1214     if(fDebug)
1215       printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1216     for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1217       AliAODMCParticle *mcp = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1218       if(!mcp)
1219         continue;
1220       Int_t pdg = mcp->GetPdgCode();
1221       if(pdg!=22)
1222         continue;
1223       if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1224         continue;
1225       Int_t imom = mcp->GetMother();
1226       if(imom<0 || imom>nTracks)
1227         continue;
1228       AliAODMCParticle *mcmom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(imom));
1229       if(!mcmom)
1230         continue;
1231       Int_t pdgMom = mcmom->GetPdgCode();
1232       Double_t mcphi = mcp->Phi();
1233       Double_t mceta = mcp->Eta();
1234       if((imom==6 || imom==7) && pdgMom==22) {
1235         fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1236         Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1237         fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1238         if(ptsum<2)
1239           fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1240         if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1241           fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1242         if(fNClusForDirPho==0)
1243           fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1244         if(fDebug){
1245           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());
1246           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());
1247         }
1248       }
1249       else{
1250         if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1251           fDecayPhotonPtMC->Fill(mcp->Pt());
1252       }
1253     }
1254   }
1255 }
1256 //________________________________________________________________________
1257 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
1258 {
1259   if(!c)
1260     return -0.1;
1261   if(!fStack && !fAODMCParticles)
1262     return -0.1;
1263   Int_t clabel = c->GetLabel();
1264   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
1265     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
1266   if(!fMcIdFamily.Contains(Form("%d",clabel)))
1267     return -0.1;
1268   fNClusForDirPho++;
1269   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
1270   Int_t partonpos = partonposstr.Atoi();
1271   if(fDebug)
1272     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
1273   Float_t clsPos[3] = {0,0,0};
1274   c->GetPosition(clsPos);
1275   TVector3 clsVec(clsPos);
1276   clsVec -= fVecPv;
1277   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
1278   //ESD
1279   if(fStack){
1280     Int_t nmcp = fStack->GetNtrack();
1281     if(clabel<0 || clabel>nmcp)
1282       return -0.1;
1283     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
1284     if(!mcp)
1285       return -0.1;
1286     if(fDebug){
1287       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);
1288     }
1289     Int_t lab1 =  mcp->GetFirstDaughter();
1290     if(lab1<0 || lab1>nmcp)
1291       return -0.1;
1292     TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
1293     if(!mcd)
1294       return -0.1;
1295     if(fDebug)
1296       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);
1297     if(mcd->GetPdgCode()==22){
1298       fClusEtMcPt->Fill(mcd->Pt(), Et);
1299       fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1300     }
1301     else{
1302       if(fDebug)
1303         printf("Warning: daughter of photon parton is not a photon\n");
1304       return -0.1;
1305     }
1306   }
1307   //AOD
1308   else if(fAODMCParticles){
1309     Int_t nmcp = fAODMCParticles->GetEntriesFast();
1310     if(clabel<0 || clabel>nmcp)
1311       return -0.1;
1312     AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(partonpos));
1313     if(!mcp)
1314       return -0.1;
1315     if(fDebug){
1316       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);
1317     }
1318     Int_t lab1 =  mcp->GetDaughter(0);
1319     if(lab1<0 || lab1>nmcp)
1320       return -0.1;
1321     AliAODMCParticle *mcd = static_cast<AliAODMCParticle*>(fAODMCParticles->At(lab1));
1322     if(!mcd)
1323       return -0.1;
1324     if(fDebug)
1325       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);
1326     if(mcd->GetPdgCode()==22){
1327       fClusEtMcPt->Fill(mcd->Pt(), Et);
1328       fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1329     }
1330     else{
1331       printf("Warning: daughter of photon parton is not a photon\n");
1332       return -0.1;
1333     }
1334   }
1335   return fDirPhoPt;
1336 }
1337 //________________________________________________________________________
1338 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
1339 {
1340   if(!fStack && !fAODMCParticles)
1341     return;
1342   Int_t selfid = 6;
1343   //ESD
1344   if(fStack){
1345     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
1346     if(!mcp)
1347       return;
1348     if(mcp->GetPdgCode()!=22){
1349       mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
1350       if(!mcp)
1351         return;
1352     }
1353     Int_t daug0f =  mcp->GetFirstDaughter();
1354     Int_t daug0l =  mcp->GetLastDaughter();
1355     Int_t nd0 = daug0l - daug0f;
1356     if(fDebug)
1357       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);
1358     fMcIdFamily = Form("%d,",selfid);
1359     GetDaughtersInfo(daug0f,daug0l, selfid,"");
1360     if(fDebug){
1361       printf("\t---- end of the prompt  gamma's family tree ----\n\n");
1362       printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1363     }
1364     TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
1365     if(!md)
1366       return;
1367     fDirPhoPt = md->Pt();
1368   }
1369   //AOD
1370   else   if(fAODMCParticles){
1371     AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(selfid));
1372     if(!mcp)
1373       return;
1374     if(mcp->GetPdgCode()!=22){
1375       mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(++selfid));
1376       if(!mcp)
1377         return;
1378     }
1379     Int_t daug0f =  mcp->GetDaughter(0);
1380     Int_t daug0l =  mcp->GetDaughter(mcp->GetNDaughters()-1);
1381     Int_t nd0 = daug0l - daug0f;
1382     if(fDebug)
1383       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);
1384     fMcIdFamily = Form("%d,",selfid);
1385     GetDaughtersInfo(daug0f,daug0l, selfid,"");
1386     if(fDebug){
1387       printf("\t---- end of the prompt  gamma's family tree ----\n\n");
1388       printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1389     }
1390     AliAODMCParticle *md = static_cast<AliAODMCParticle*>(fAODMCParticles->At(daug0f));
1391     if(!md)
1392       return;
1393     fDirPhoPt = md->Pt();
1394   }
1395
1396 }
1397 //________________________________________________________________________
1398 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
1399 {
1400   if(fStack){
1401     int nmcp = fStack->GetNtrack();
1402     if(firstd<0 || firstd>nmcp)
1403       return;
1404     if(lastd<0 || lastd>nmcp)
1405       return;
1406     if(firstd>lastd){
1407       printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1408     return;
1409     }
1410     TString indenter = Form("\t%s",inputind);
1411     TParticle *dp = 0x0;
1412     if(fDebug)
1413       printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1414     for(int id=firstd; id<lastd+1; id++){
1415       dp =  static_cast<TParticle*>(fStack->Particle(id));
1416       if(!dp)
1417         continue;
1418       Int_t dfd = dp->GetFirstDaughter();
1419       Int_t dld = dp->GetLastDaughter();
1420       Int_t dnd =  dld - dfd + 1;
1421       Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
1422       if(fDebug)
1423         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);
1424       fMcIdFamily += Form("%d,",id);
1425       GetDaughtersInfo(dfd,dld,id,indenter.Data());
1426     }
1427   }
1428   if(fAODMCParticles){
1429     int nmcp = fAODMCParticles->GetEntriesFast();
1430     if(firstd<0 || firstd>nmcp)
1431       return;
1432     if(lastd<0 || lastd>nmcp)
1433       return;
1434     if(firstd>lastd){
1435       printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1436     return;
1437     }
1438     TString indenter = Form("\t%s",inputind);
1439     AliAODMCParticle *dp = 0x0;
1440     if(fDebug)
1441       printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1442     for(int id=firstd; id<lastd+1; id++){
1443       dp =  static_cast<AliAODMCParticle*>(fAODMCParticles->At(id));
1444       if(!dp)
1445         continue;
1446       Int_t dfd = dp->GetDaughter(0);
1447       Int_t dld = dp->GetDaughter(dp->GetNDaughters()-1);
1448       Int_t dnd =  dld - dfd + 1;
1449       Float_t vr = TMath::Sqrt(dp->Xv()*dp->Xv()+dp->Xv()*dp->Xv());
1450       if(fDebug)
1451         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);
1452       fMcIdFamily += Form("%d,",id);
1453       GetDaughtersInfo(dfd,dld,id,indenter.Data());
1454     }
1455   }
1456 }
1457
1458 //________________________________________________________________________
1459 Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
1460 {
1461   if(!fStack && !fAODMCParticles)
1462     return 0;
1463   if(fDebug)
1464     printf("Inside GetMcPtSumInCone!!\n");
1465   Float_t ptsum = 0;
1466   TString addpartlabels = "";
1467   //ESD
1468   if(fStack){
1469     Int_t nTracks = fStack->GetNtrack();
1470     for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1471       TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1472       if(!mcp)
1473         continue;
1474       Int_t status = mcp->GetStatusCode();
1475       if(status!=1)
1476         continue;
1477       Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz());
1478       if(mcvr>1)
1479         continue;
1480       /*else {
1481         if(fDebug)
1482         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1483         }*/
1484       Float_t dphi = mcp->Phi() - phiclus;
1485       Float_t deta = mcp->Eta() - etaclus;
1486       if(fDebug && TMath::Abs(dphi)<0.01)
1487         printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1488
1489       if(deta>R || dphi>R)
1490         continue;
1491       Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1492       if(dR>R)
1493         continue;
1494       addpartlabels += Form("%d,",iTrack);
1495       if(mcp->Pt()<0.2)
1496         continue;
1497       ptsum += mcp->Pt();
1498     }
1499   }
1500   //AOD
1501   if(fAODMCParticles){
1502     Int_t nTracks = fAODMCParticles->GetEntriesFast();
1503     for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1504       AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1505       if(!mcp)
1506         continue;
1507       Int_t status = mcp->GetStatus();
1508       if(status!=1)
1509         continue;
1510       Float_t mcvr = TMath::Sqrt(mcp->Xv()*mcp->Xv()+ mcp->Yv()* mcp->Yv() + mcp->Zv()*mcp->Zv());
1511       if(mcvr>1)
1512         continue;
1513       /*else {
1514         if(fDebug)
1515         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1516         }*/
1517       Float_t dphi = mcp->Phi() - phiclus;
1518       Float_t deta = mcp->Eta() - etaclus;
1519       if(fDebug && TMath::Abs(dphi)<0.01)
1520         printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1521
1522       if(deta>R || dphi>R)
1523         continue;
1524       Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1525       if(dR>R)
1526         continue;
1527       addpartlabels += Form("%d,",iTrack);
1528       if(mcp->Pt()<0.2)
1529         continue;
1530       ptsum += mcp->Pt();
1531     }
1532   }
1533   return ptsum;
1534 }
1535 //________________________________________________________________________
1536 void AliAnalysisTaskEMCALIsoPhoton::FillQA()
1537 {
1538
1539   TObjArray *clusters = fESDClusters;
1540   //"none", "exotic", "exo+cpv1", "exo+cpv1+time", "exo+cpv1+time+m02"),
1541   if (!clusters){
1542     clusters = fAODClusters;
1543     if(fDebug)
1544       printf("ESD clusters empty...");
1545   }
1546   if (!clusters){
1547     if(fDebug)
1548       printf("  and AOD clusters as well!!!\n");
1549     return;
1550   }
1551   if(!fSelPrimTracks)
1552     return;
1553   const int ntracks = fSelPrimTracks->GetEntriesFast();
1554   const int ncells = fNCells50;//fESDCells->GetNumberOfCells();
1555   const Int_t nclus = clusters->GetEntries();
1556   fNTracks->Fill(ntracks);
1557   fEmcNCells->Fill(ncells);
1558   fEmcNClus->Fill(nclus);
1559   if(fMaxEClus>fECut){
1560     fNTracksECut->Fill(ntracks);
1561     fEmcNCellsCut->Fill(ncells);
1562     fEmcNClusCut->Fill(nclus);
1563   }
1564   for(int it=0;it<ntracks;it++){
1565     AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it);
1566     if(!t)
1567       continue;
1568     fTrackPtPhi->Fill(t->Pt(),t->Phi());
1569     fTrackPtEta->Fill(t->Pt(),t->Eta());
1570     if(fMaxEClus>fECut){
1571       fTrackPtPhiCut->Fill(t->Pt(), t->Phi());
1572       fTrackPtEtaCut->Fill(t->Pt(), t->Eta());
1573     }
1574     if(t->GetTPCsignal()<80 || t->GetTPCsignal()>100)
1575       continue;
1576     if(t->GetEMCALcluster()<=0 || t->GetEMCALcluster()>nclus)
1577       continue;
1578     AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(t->GetEMCALcluster()));
1579     if(!c)
1580       continue;
1581     fEoverPvsE->Fill(c->E(),c->E()/t->P());
1582   }
1583   for(int ic=0;ic<nclus;ic++){
1584     AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(ic));
1585     //AliESDCaloCluster *c = (AliESDCaloCluster*)clusters->At(ic);
1586     if(!c)
1587       continue;
1588     if(!c->IsEMCAL())
1589       continue;
1590     Float_t clsPos[3] = {0,0,0};
1591     c->GetPosition(clsPos);
1592     TVector3 clsVec(clsPos);
1593     clsVec -= fVecPv;
1594     Double_t cphi = clsVec.Phi();
1595     Double_t ceta = clsVec.Eta();
1596     Short_t id = -1;
1597     GetMaxCellEnergy( c, id);
1598     fEmcClusEClusCuts->Fill(c->E(),0);
1599     fEmcClusEPhi->Fill(c->E(), cphi);
1600     fEmcClusEEta->Fill(c->E(), ceta);
1601     if(fMaxEClus>fECut){
1602       fEmcClusEPhiCut->Fill(c->E(), cphi);
1603       fEmcClusEEtaCut->Fill(c->E(), ceta);
1604     }
1605     Double_t maxt=0;
1606     if(fESDCells)
1607       maxt = fESDCells->GetCellTime(id);
1608     else if(fAODCells)
1609       maxt = fAODCells->GetCellTime(id);
1610     if(IsExotic(c))
1611       continue;
1612     fEmcClusNotExo->Fill(c->E());
1613     fEmcClusEClusCuts->Fill(c->E(),1);
1614     if(fClusIdFromTracks.Contains(Form("%d",ic))){
1615       fEmcClusETM2->Fill(c->E());
1616       fDetaDphiFromTM->Fill(c->GetTrackDz(),c->GetTrackDx());
1617     }
1618     if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
1619       fEmcClusETM1->Fill(c->E());
1620       continue;
1621     }
1622     fEmcClusEClusCuts->Fill(c->E(),2);
1623     if(TMath::Abs(maxt)>30e-9)
1624       continue;
1625     fEmcClusEClusCuts->Fill(c->E(),3);
1626     if(c->GetM02()>0.1)
1627       fEmcClusEClusCuts->Fill(c->E(),4);
1628   }
1629 }
1630 //________________________________________________________________________
1631 Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex)
1632 {
1633   Double_t pt = 0;
1634   if(!fTracks)
1635     return pt;
1636   if(matchIndex<0 || matchIndex>fTracks->GetEntries()){
1637     if(fDebug)
1638       printf("track-matched index out of track array range!!!\n");
1639     return pt;
1640   }
1641   AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex));
1642   if(!track){
1643     if(fDebug)
1644       printf("track-matched pointer does not exist!!!\n");
1645     return pt;
1646   }
1647   if(fESD){
1648     if(!fPrTrCuts && !fCompTrCuts)
1649       return pt;
1650     if(!fPrTrCuts->IsSelected(track) && !fCompTrCuts->IsSelected(track))
1651       return pt;
1652     pt = track->Pt();
1653   }
1654   return pt;
1655 }
1656 //________________________________________________________________________
1657 void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells()
1658 {
1659   AliVCaloCells *cells = 0;
1660   cells = fESDCells;
1661   if (!cells)
1662     cells = fAODCells;
1663   if(!cells)
1664     return;
1665   Double_t maxe = 0;
1666   Double_t maxphi = -10;
1667   Int_t ncells = cells->GetNumberOfCells();
1668   Double_t eta,phi;
1669   for (Int_t i=0; i<ncells; i++) {
1670     Short_t absid = TMath::Abs(cells->GetCellNumber(i));
1671     Double_t e = cells->GetCellAmplitude(absid);
1672     if(e>0.05)
1673       fNCells50++;
1674     else
1675       continue;
1676     fGeom->EtaPhiFromIndex(absid,eta,phi);
1677     if(maxe<e){
1678       maxe = e;
1679       maxphi = phi;
1680     }
1681   }
1682   fMaxCellEPhi->Fill(maxe,maxphi);
1683
1684 }
1685 //________________________________________________________________________
1686 bool AliAnalysisTaskEMCALIsoPhoton::IsExotic(AliVCluster *c)
1687 {
1688   bool isExo = 0;
1689   Short_t id = -1;
1690   Double_t Emax = GetMaxCellEnergy( c, id);
1691   Double_t Ecross = GetCrossEnergy( c, id);
1692   if((1-Ecross/Emax)>fExoticCut)
1693     isExo = 1;
1694   return isExo;
1695 }
1696 //________________________________________________________________________
1697 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
1698 {
1699   // Called once at the end of the query.
1700 }