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