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