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