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