PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx

   1 // $Id$
   2
   3 #include "AliAnalysisTaskEMCALIsoPhoton.h"
   4
   5 #include <TCanvas.h>
   6 #include <TChain.h>
   7 #include <TFile.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <THnSparse.h>
  12 #include <TLorentzVector.h>
  13 #include <TList.h>
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliAnalysisTask.h"
  17 #include "AliEMCALGeometry.h"
  18 #include "AliEMCALRecoUtils.h"
  19 #include "AliESDCaloCells.h"
  20 #include "AliESDCaloCluster.h"
  21 #include "AliESDEvent.h"
  22 #include "AliESDHeader.h"
  23 #include "AliESDInputHandler.h"
  24 #include "AliESDUtils.h"
  25 #include "AliESDtrack.h"
  26 #include "AliESDtrackCuts.h"
  27 #include "AliAODEvent.h"
  28 #include "AliAODTrack.h"
  29 #include "AliMCEvent.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliStack.h"
  32 #include "AliAODMCParticle.h"
  33 #include "AliVEvent.h"
  34 #include "AliVTrack.h"
  35 #include "AliV0vertexer.h"
  36 #include "AliVCluster.h"
  37 #include "AliOADBContainer.h"
  38
  39 #include <iostream>
  40 using std::cout;
  41 using std::endl;
  42
  43 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
  44
  45 //________________________________________________________________________
  46 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
  47   AliAnalysisTaskSE(),
  48   fESDClusters(0),
  49   fAODClusters(0),
  50   fSelPrimTracks(0),
  51   fTracks(0),
  52   fAODMCParticles(0),
  53   fESDCells(0),
  54   fAODCells(0),
  55   fPrTrCuts(0),
  56   fCompTrCuts(0),
  57   fGeom(0x0),
  58   fGeoName("EMCAL_COMPLETEV1"),
  59   fOADBContainer(0),
  60   fVecPv(0.,0.,0.),
  61   fPeriod("LHC11c"),
  62   fTrigBit("kEMC7"),
  63   fIsTrain(0),
  64   fIsMc(0),
  65   fDebug(0),
  66   fPathStrOpt("/"),
  67   fExoticCut(0.97),
  68   fIsoConeR(0.4),
  69   fNDimensions(7),
  70   fECut(3.),
  71   fTrackMult(0),
  72   fMcIdFamily(""),
  73   fNClusForDirPho(0),
  74   fDirPhoPt(0),
  75   fHigherPtCone(0),
  76   fImportGeometryFromFile(0),
  77   fImportGeometryFilePath(""),
  78   fMaxPtTrack(0),
  79   fMaxEClus(0),
  80   fNCells50(0),
  81   fFilterBit(0),
  82   fSelHybrid(kFALSE),
  83   fFillQA(kFALSE),
  84   fClusIdFromTracks(""),
  85   fCpvFromTrack(kFALSE),
  86   fNBinsPt(200),
  87   fPtBinLowEdge(0),
  88   fPtBinHighEdge(100),
  89   fRemMatchClus(kFALSE),
  90   fMinIsoClusE(0),
  91   fNCuts(5),
  92   fTrCoreRem(kFALSE),
  93   fClusTDiff(30e-9),
  94   fPileUpRejSPD(kFALSE),
  95   fDistToBadChan(0),
  96   fInConeInvMass(""),
  97   fInConePairClEt(""),
  98   fESD(0),
  99   fAOD(0),
 100   fVEvent(0),
 101   fMCEvent(0),
 102   fStack(0),
 103   fOutputList(0),
 104   fEvtSel(0),
 105   fNClusEt10(0),
 106   fClusArrayNames(0),
 107   fRecoPV(0),
 108   fPVtxZ(0),
 109   fTrMultDist(0),
 110   fClusEtCPVSBGISO(0),
 111   fClusEtCPVBGISO(0),
 112   fMCDirPhotonPtEtaPhi(0),
 113   fMCIsoDirPhotonPtEtaPhi(0),
 114   fMCDirPhotonPtEtIso(0),
 115   fDecayPhotonPtMC(0),
 116   fCellAbsIdVsAmpl(0),
 117   fNClusHighClusE(0),
 118   fHigherPtConeM02(0),
 119   fClusEtMcPt(0),
 120   fClusMcDetaDphi(0),
 121   fNClusPerPho(0),
 122   fMcPtInConeBG(0),
 123   fMcPtInConeSBG(0),
 124   fMcPtInConeBGnoUE(0),
 125   fMcPtInConeSBGnoUE(0),
 126   fMcPtInConeTrBGnoUE(0),
 127   fMcPtInConeTrSBGnoUE(0),
 128   fMcPtInConeMcPhoPt(0),
 129   fAllIsoEtMcGamma(0),
 130   fAllIsoNoUeEtMcGamma(0),
 131   fMCDirPhotonPtEtaPhiNoClus(0),
 132   fEtCandIsoAndIsoWoPairEt(0),
 133   fInConePairedClusEtVsCandEt(0),
 134   fHnOutput(0),
 135   fQAList(0),
 136   fNTracks(0),
 137   fEmcNCells(0),
 138   fEmcNClus(0),
 139   fEmcNClusCut(0),
 140   fNTracksECut(0),
 141   fEmcNCellsCut(0),
 142   fEmcClusETM1(0),
 143   fEmcClusETM2(0),
 144   fEmcClusNotExo(0),
 145   fEmcClusEClusCuts(0),
 146   fEmcClusEPhi(0),
 147   fEmcClusEPhiCut(0),
 148   fEmcClusEEta(0),
 149   fEmcClusEEtaCut(0),
 150   fTrackPtPhi(0),
 151   fTrackPtPhiCut(0),
 152   fTrackPtEta(0),
 153   fTrackPtEtaCut(0),
 154   fMaxCellEPhi(0),
 155   fDetaDphiFromTM(0),
 156   fEoverPvsE(0)
 157 {
 158   // Default constructor.
 159   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 160 }
 161
 162 //________________________________________________________________________
 163 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
 164   AliAnalysisTaskSE(name),
 165   fESDClusters(0),
 166   fAODClusters(0),
 167   fSelPrimTracks(0),
 168   fTracks(0),
 169   fAODMCParticles(0),
 170   fESDCells(0),
 171   fAODCells(0),
 172   fPrTrCuts(0),
 173   fCompTrCuts(0),
 174   fGeom(0x0),
 175   fGeoName("EMCAL_COMPLETEV1"),
 176   fOADBContainer(0),
 177   fVecPv(0.,0.,0.),
 178   fPeriod("LHC11c"),
 179   fTrigBit("kEMC7"),
 180   fIsTrain(0),
 181   fIsMc(0),
 182   fDebug(0),
 183   fPathStrOpt("/"),
 184   fExoticCut(0.97),
 185   fIsoConeR(0.4),
 186   fNDimensions(7),
 187   fECut(3.),
 188   fTrackMult(0),
 189   fMcIdFamily(""),
 190   fNClusForDirPho(0),
 191   fDirPhoPt(0),
 192   fHigherPtCone(0),
 193   fImportGeometryFromFile(0),
 194   fImportGeometryFilePath(""),
 195   fMaxPtTrack(0),
 196   fMaxEClus(0),
 197   fNCells50(0),
 198   fFilterBit(0),
 199   fSelHybrid(kFALSE),
 200   fFillQA(kFALSE),
 201   fClusIdFromTracks(""),
 202   fCpvFromTrack(kFALSE),
 203   fNBinsPt(200),
 204   fPtBinLowEdge(0.),
 205   fPtBinHighEdge(100),
 206   fRemMatchClus(kFALSE),
 207   fMinIsoClusE(0),
 208   fNCuts(5),
 209   fTrCoreRem(kFALSE),
 210   fClusTDiff(30e-9),
 211   fPileUpRejSPD(kFALSE),
 212   fDistToBadChan(0),
 213   fInConeInvMass(""),
 214   fInConePairClEt(""),
 215   fESD(0),
 216   fAOD(0),
 217   fVEvent(0),
 218   fMCEvent(0),
 219   fStack(0),
 220   fOutputList(0),
 221   fEvtSel(0),
 222   fNClusEt10(0),
 223   fClusArrayNames(0),
 224   fRecoPV(0),
 225   fPVtxZ(0),
 226   fTrMultDist(0),
 227   fClusEtCPVSBGISO(0),
 228   fClusEtCPVBGISO(0),
 229   fMCDirPhotonPtEtaPhi(0),
 230   fMCIsoDirPhotonPtEtaPhi(0),
 231   fMCDirPhotonPtEtIso(0),
 232   fDecayPhotonPtMC(0),
 233   fCellAbsIdVsAmpl(0),
 234   fNClusHighClusE(0),
 235   fHigherPtConeM02(0),
 236   fClusEtMcPt(0),
 237   fClusMcDetaDphi(0),
 238   fNClusPerPho(0),
 239   fMcPtInConeBG(0),
 240   fMcPtInConeSBG(0),
 241   fMcPtInConeBGnoUE(0),
 242   fMcPtInConeSBGnoUE(0),
 243   fMcPtInConeTrBGnoUE(0),
 244   fMcPtInConeTrSBGnoUE(0),
 245   fMcPtInConeMcPhoPt(0),
 246   fAllIsoEtMcGamma(0),
 247   fAllIsoNoUeEtMcGamma(0),
 248   fMCDirPhotonPtEtaPhiNoClus(0),
 249   fEtCandIsoAndIsoWoPairEt(0),
 250   fInConePairedClusEtVsCandEt(0),
 251   fHnOutput(0),
 252   fQAList(0),
 253   fNTracks(0),
 254   fEmcNCells(0),
 255   fEmcNClus(0),
 256   fEmcNClusCut(0),
 257   fNTracksECut(0),
 258   fEmcNCellsCut(0),
 259   fEmcClusETM1(0),
 260   fEmcClusETM2(0),
 261   fEmcClusNotExo(0),
 262   fEmcClusEClusCuts(0),
 263   fEmcClusEPhi(0),
 264   fEmcClusEPhiCut(0),
 265   fEmcClusEEta(0),
 266   fEmcClusEEtaCut(0),
 267   fTrackPtPhi(0),
 268   fTrackPtPhiCut(0),
 269   fTrackPtEta(0),
 270   fTrackPtEtaCut(0),
 271   fMaxCellEPhi(0),
 272   fDetaDphiFromTM(0),
 273   fEoverPvsE(0)
 274 {
 275   // Constructor
 276
 277   // Define input and output slots here
 278   // Input slot #0 works with a TChain
 279   DefineInput(0, TChain::Class());
 280   // Output slot #0 id reserved by the base class for AOD
 281   // Output slot #1 writes into a TH1 container
 282   DefineOutput(1, TList::Class());
 283   DefineOutput(2, TList::Class());
 284 }
 285
 286 //________________________________________________________________________
 287 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 288 {
 289   // Create histograms, called once.
 290
 291   fESDClusters = new TObjArray();
 292   fAODClusters = new TObjArray();
 293   fSelPrimTracks = new TObjArray();
 294
 295
 296   fOutputList = new TList();
 297   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 298
 299   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
 300   fOADBContainer = new AliOADBContainer("AliEMCALgeo");
 301   fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
 302
 303   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 304   fOutputList->Add(fEvtSel);
 305
 306   fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
 307   fOutputList->Add(fNClusEt10);
 308
 309   fClusArrayNames = new TH1F("hClusArrayNames","cluster array names (0=CaloClusters,1=EmcCaloClusters,2=Others);option;#events",3,0,3);
 310   fOutputList->Add(fClusArrayNames);
 311
 312   fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
 313   fOutputList->Add(fRecoPV);
 314
 315   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 316   fOutputList->Add(fPVtxZ);
 317
 318   fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5);
 319   fOutputList->Add(fTrMultDist);
 320
 321   fClusEtCPVSBGISO = new TH2F("hClusEtCPVSBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.1<#lambda_{0}^{2}<0.3;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100);
 322   fClusEtCPVSBGISO->Sumw2();
 323   fOutputList->Add(fClusEtCPVSBGISO);
 324
 325   fClusEtCPVBGISO = new TH2F("hClusEtCPVBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.5<#lambda_{0}^{2}<2.0;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100);
 326   fClusEtCPVBGISO->Sumw2();
 327   fOutputList->Add(fClusEtCPVBGISO);
 328
 329   fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20);
 330   fMCDirPhotonPtEtaPhi->Sumw2();
 331   fOutputList->Add(fMCDirPhotonPtEtaPhi);
 332
 333   fMCIsoDirPhotonPtEtaPhi = new TH3F("hMCIsoDirPhotonPtEtaPhi","photon (gq->#gammaq, isolated@MC) p_{T}, #eta, #phi;GeV/c;#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20);
 334   fMCIsoDirPhotonPtEtaPhi->Sumw2();
 335   fOutputList->Add(fMCIsoDirPhotonPtEtaPhi);
 336
 337   fMCDirPhotonPtEtIso = new TH2F("hMCDirPhotonPtEtIso",Form("photon (gq->#gammaq @MC) p_{T}, E_{T}^{ISO} (R=%1.1f);GeV/c;E_{T}^{ISO} GeV/c",fIsoConeR),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,20,-0.25,9.75);
 338   fMCDirPhotonPtEtIso->Sumw2();
 339   fOutputList->Add(fMCDirPhotonPtEtIso);
 340
 341
 342   fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 343   fDecayPhotonPtMC->Sumw2();
 344   fOutputList->Add(fDecayPhotonPtMC);
 345
 346   fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,24*48*10,-0.5,24*48*10-0.5);
 347   fOutputList->Add(fCellAbsIdVsAmpl);
 348
 349   fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5);
 350   fOutputList->Add(fNClusHighClusE);
 351
 352   fHigherPtConeM02 = new TH2F("hHigherPtConeM02","#lambda_{0}^{2} vs. in-cone-p_{T}^{max};p_{T}^{max} (GeV/c, in the cone);#lambda_{0}^{2}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,400,0,4);
 353   fOutputList->Add(fHigherPtConeM02);
 354
 355   fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100);
 356   fOutputList->Add(fClusEtMcPt);
 357
 358   fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7);
 359   fOutputList->Add(fClusMcDetaDphi);
 360
 361   fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5);
 362   fOutputList->Add(fNClusPerPho);
 363
 364   fMcPtInConeBG = new TH2F("hMcPtInConeBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 365   fOutputList->Add(fMcPtInConeBG);
 366
 367   fMcPtInConeSBG  = new TH2F("hMcPtInConeSBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 368   fOutputList->Add(fMcPtInConeSBG);
 369
 370   fMcPtInConeBGnoUE = new TH2F("hMcPtInConeBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 371   fOutputList->Add(fMcPtInConeBGnoUE);
 372
 373   fMcPtInConeSBGnoUE  = new TH2F("hMcPtInConeSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 374   fOutputList->Add(fMcPtInConeSBGnoUE);
 375
 376   fMcPtInConeTrBGnoUE = new TH2F("hMcPtInConeTrBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (BG template);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 377   fOutputList->Add(fMcPtInConeTrBGnoUE);
 378
 379   fMcPtInConeTrSBGnoUE  = new TH2F("hMcPtInConeTrSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (SBG range);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 380   fOutputList->Add(fMcPtInConeTrSBGnoUE);
 381
 382   fMcPtInConeMcPhoPt  = new TH2F("hMcPtInConeMcPhoPt","#sum_{in-cone}p_{T}^{mc-primaries} vs. prompt photon p_{T};p_{T}^{mc-#gamma} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,200,-0.25,99.75);
 383   fOutputList->Add(fMcPtInConeMcPhoPt);
 384
 385   fAllIsoEtMcGamma  = new TH2F("hAllIsoEtMcGammaE","ISO^{TRK+EMC} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50);
 386   fOutputList->Add(fAllIsoEtMcGamma);
 387
 388   fAllIsoNoUeEtMcGamma  = new TH2F("hAllIsoNoUeEtMcGammaE","ISO^{TRK+EMC}_{noue} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC}_{noue} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50);
 389   fOutputList->Add(fAllIsoNoUeEtMcGamma);
 390
 391
 392   fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and  #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20);
 393   fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus);
 394
 395   fEtCandIsoAndIsoWoPairEt = new TH3F("hEtCandIsoAndIsoWoPairEt","E_{T}^{cand} vs. E_{T}^{ISO} (EMC+Trk) (0.1<M02<0.3, 0.110<m_{#gamma#gamma}<0.165 only);E_{T}^{cand}; E_{T}^{ISO}; E_{T}^{ISO} (w/o #pi^{0} pair E_{T})",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,200,1000,0,200);
 396   fOutputList->Add(fEtCandIsoAndIsoWoPairEt);
 397
 398   fInConePairedClusEtVsCandEt = new TH2F("hInConePairedClusEtVsCandEt","E_{T}^{partner} vs. E_{T}^{cand} (R<0.4, 0.110<m_{#gamma#gamma}<0.165);E_{T}^{cand};E_{T}^{partner}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,200,0,40);
 399   fOutputList->Add(fInConePairedClusEtVsCandEt);
 400
 401   Int_t nEt=fNBinsPt*5, nM02=400, nCeIso=1000, nTrIso=1000,  nAllIso=1000,  nCeIsoNoUE=1000,  nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=100, nPhoMcPt=fNBinsPt, nInConeMass=100;
 402   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt,nInConeMass};
 403   fNDimensions = sizeof(bins)/sizeof(Int_t);
 404   const Int_t ndims =   fNDimensions;
 405   Double_t xmin[] = { fPtBinLowEdge,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,fPtBinLowEdge,0.0};
 406   Double_t xmax[] = { fPtBinHighEdge, 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,fPtBinHighEdge, 1.0};
 407   if(fPeriod.Contains("11h")){
 408     xmax[12]=3999.5;
 409   }
 410   fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult,mc-p_{T}^{#gamma}", ndims, bins, xmin, xmax);
 411   fHnOutput->Sumw2();
 412   fOutputList->Add(fHnOutput);
 413
 414   //QA outputs
 415   fQAList = new TList();
 416   fQAList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 417
 418   fNTracks = new TH1F("hNTracks","# of selected tracks;n-tracks;counts",120,-0.5,119.5);
 419   fNTracks->Sumw2();
 420   fQAList->Add(fNTracks);
 421
 422   fEmcNCells = new TH1F("fEmcNCells",";n/event;count",120,-0.5,119.5);
 423   fEmcNCells->Sumw2();
 424   fQAList->Add(fEmcNCells);
 425   fEmcNClus = new TH1F("fEmcNClus",";n/event;count",120,-0.5,119.5);
 426   fEmcNClus->Sumw2();
 427   fQAList->Add(fEmcNClus);
 428   fEmcNClusCut = new TH1F("fEmcNClusCut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5);
 429   fEmcNClusCut->Sumw2();
 430   fQAList->Add(fEmcNClusCut);
 431   fNTracksECut = new TH1F("fNTracksECut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5);
 432   fNTracksECut->Sumw2();
 433   fQAList->Add(fNTracksECut);
 434   fEmcNCellsCut = new TH1F("fEmcNCellsCut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5);
 435   fEmcNCellsCut->Sumw2();
 436   fQAList->Add(fEmcNCellsCut);
 437   fEmcClusETM1 = new TH1F("fEmcClusETM1","(method clus->GetTrackDx,z);GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 438   fEmcClusETM1->Sumw2();
 439   fQAList->Add(fEmcClusETM1);
 440   fEmcClusETM2 = new TH1F("fEmcClusETM2","(method track->GetEMCALcluster());GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 441   fEmcClusETM2->Sumw2();
 442   fQAList->Add(fEmcClusETM2);
 443   fEmcClusNotExo  = new TH1F("fEmcClusNotExo","exotics removed;GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge);
 444   fEmcClusNotExo->Sumw2();
 445   fQAList->Add(fEmcClusNotExo);
 446   fEmcClusEPhi = new TH2F("fEmcClusEPhi",";GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 447   fEmcClusEPhi->Sumw2();
 448   fQAList->Add(fEmcClusEPhi);
 449   fEmcClusEPhiCut = new TH2F("fEmcClusEPhiCut",Form("(at least one E_{clus}>%1.1f);GeV;#phi",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 450   fEmcClusEPhiCut->Sumw2();
 451   fQAList->Add(fEmcClusEPhiCut);
 452   fEmcClusEEta = new TH2F("fEmcClusEEta",";GeV;#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,19,-0.9,0.9);
 453   fEmcClusEEta->Sumw2();
 454   fQAList->Add(fEmcClusEEta);
 455   fEmcClusEEtaCut = new TH2F("fEmcClusEEtaCut",Form("(at least one E_{clus}>%1.1f);GeV;#eta",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
 456   fEmcClusEEtaCut->Sumw2();
 457   fQAList->Add(fEmcClusEEtaCut);
 458   fTrackPtPhi = new TH2F("fTrackPtPhi",";p_{T} [GeV/c];#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 459   fTrackPtPhi->Sumw2();
 460   fQAList->Add(fTrackPtPhi);
 461   fTrackPtPhiCut = new TH2F("fTrackPtPhiCut",Form("(at least one E_{clus}>%1.1f);p_{T} [GeV/c];#phi",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 462   fTrackPtPhiCut->Sumw2();
 463   fQAList->Add(fTrackPtPhiCut);
 464   fTrackPtEta = new TH2F("fTrackPtEta",";p_{T} [GeV/c];#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
 465   fTrackPtEta->Sumw2();
 466   fQAList->Add(fTrackPtEta);
 467   fTrackPtEtaCut = new TH2F("fTrackPtEtaCut",Form("(at least one E_{clus}>%1.1f);p_{T} [GeV/c];#eta",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9);
 468   fTrackPtEtaCut->Sumw2();
 469   fQAList->Add(fTrackPtEtaCut);
 470   fEmcClusEClusCuts = new TH2F("fEmcClusEClusCuts",";GeV;cut",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,fNCuts,-0.5,fNCuts-0.5);
 471   fEmcClusEClusCuts->Sumw2();
 472   fQAList->Add(fEmcClusEClusCuts);
 473
 474   fMaxCellEPhi = new TH2F("fMaxCellEPhi","Most energetic cell in event; GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3);
 475   fMaxCellEPhi->Sumw2();
 476   fQAList->Add(fMaxCellEPhi);
 477
 478   fDetaDphiFromTM = new TH2F("fDetaDphiFromTM","d#phi vs. d#eta of clusters from track->GetEMCALcluster();d#eta;d#phi",100,-0.05,0.05,200,-0.1,0.1);
 479   fDetaDphiFromTM->Sumw2();
 480   fQAList->Add(fDetaDphiFromTM);
 481
 482   fEoverPvsE = new TH2F("fEoverPvsE","E^{clus}/p^{track} vs E^{clus} (80<TPCsignal<100);E^{clus} [GeV];E^{clus}/p^{track} [c^{-1}]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,100,0,2);
 483   fEoverPvsE->Sumw2();
 484   fQAList->Add(fEoverPvsE);
 485
 486   PostData(1, fOutputList);
 487   PostData(2, fQAList);
 488 }
 489
 490 //________________________________________________________________________
 491 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 492 {
 493   // Main loop, called for each event.
 494   fESDClusters = 0;
 495   fESDCells = 0;
 496   fAODClusters = 0;
 497   fAODCells = 0;
 498   // event trigger selection
 499   Bool_t isSelected = 0;
 500   if(fPeriod.Contains("11a")){
 501     if(fTrigBit.Contains("kEMC"))
 502       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 503     else
 504       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 505   }
 506   else{
 507     if(fTrigBit.Contains("kEMC"))
 508       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 509     else
 510       if(fTrigBit.Contains("kMB"))
 511         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 512       else
 513         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
 514   }
 515   if(fPeriod.Contains("11h")){
 516     if(fTrigBit.Contains("kAny"))
 517       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
 518     if(fTrigBit.Contains("kAnyINT"))
 519       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT);
 520   }
 521   if(fIsMc)
 522     isSelected = kTRUE;
 523   if(fDebug)
 524     printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
 525   if(!isSelected )
 526         return;
 527   if(fIsMc){
 528     TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
 529     TFile *file = (TFile*)tree->GetCurrentFile();
 530     TString filename = file->GetName();
 531     if(!filename.Contains(fPathStrOpt.Data()))
 532       return;
 533   }
 534   fVEvent = (AliVEvent*)InputEvent();
 535   if (!fVEvent) {
 536     printf("ERROR: event not available\n");
 537     return;
 538   }
 539   Int_t   runnumber = InputEvent()->GetRunNumber() ;
 540   if(fDebug)
 541     printf("run number = %d\n",runnumber);
 542
 543   fESD = dynamic_cast<AliESDEvent*>(fVEvent);
 544   if(!fESD){
 545     fAOD = dynamic_cast<AliAODEvent*>(fVEvent);
 546     if(!fAOD){
 547       printf("ERROR: Invalid type of event!!!\n");
 548       return;
 549     }
 550     else if(fDebug)
 551       printf("AOD EVENT!\n");
 552   }
 553
 554   fEvtSel->Fill(0);
 555   if(fDebug)
 556     printf("event is ok,\n run number=%d\n",runnumber);
 557
 558
 559   AliVVertex *pv = (AliVVertex*)fVEvent->GetPrimaryVertex();
 560   Bool_t pvStatus = kTRUE;
 561   if(fESD){
 562     AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex();
 563     pvStatus = esdv->GetStatus();
 564   }
 565   /*if(fAOD){
 566     AliAODVertex *aodv = (AliAODVertex*)fAOD->GetPrimaryVertex();
 567     pvStatus = aodv->GetStatus();
 568     }*/
 569   if(!pv)
 570     return;
 571   if(!pvStatus)
 572     fRecoPV->Fill(0);
 573   else
 574     fRecoPV->Fill(1);
 575   fPVtxZ->Fill(pv->GetZ());
 576   fVecPv.SetXYZ(pv->GetX(),pv->GetY(),pv->GetZ());
 577   if(TMath::Abs(pv->GetZ())>10)
 578     return;
 579   if(fDebug)
 580     printf("passed vertex cut\n");
 581
 582   fEvtSel->Fill(1);
 583   if(fVEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.) && fPileUpRejSPD){
 584     if(fDebug)
 585       printf("Event is SPD pile-up!***\n");
 586     return;
 587   }
 588   if(fESD)
 589     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 590   if(fAOD)
 591     fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
 592
 593   if(!fTracks){
 594     AliError("Track array in event is NULL!");
 595     if(fDebug)
 596       printf("returning due to not finding tracks!\n");
 597     return;
 598   }
 599   // Track loop to fill a pT spectrum
 600   const Int_t Ntracks = fTracks->GetEntriesFast();
 601   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 602     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 603     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 604     AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
 605     if (!track)
 606       continue;
 607     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
 608     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 609     if(esdTrack){
 610       if(esdTrack->GetEMCALcluster()>0)
 611         fClusIdFromTracks.Append(Form("%d ",esdTrack->GetEMCALcluster()));
 612       if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
 613         fSelPrimTracks->Add(track);
 614       } else if(fCompTrCuts && fCompTrCuts->IsSelected(track)) {
 615         fSelPrimTracks->Add(track);
 616       }
 617     }
 618     else if(aodTrack){
 619       if (fSelHybrid && !aodTrack->IsHybridGlobalConstrainedGlobal())
 620         continue ;
 621       if(!fSelHybrid && !aodTrack->TestFilterBit(fFilterBit))
 622         continue;
 623       fSelPrimTracks->Add(track);
 624       /*if(fTrackMaxPt<track->Pt())
 625         fTrackMaxPt = track->Pt();*/
 626     }
 627   }
 628
 629   TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
 630   for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 631     if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 632       break;
 633     /*if(fESD)
 634       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 635       else*/
 636     // if(fVEvent->GetEMCALMatrix(mod))
 637     fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
 638     fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
 639   }
 640
 641   if(fESD){
 642     AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 643     fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 644     if(fDebug)
 645       printf("ESD Track mult= %d\n",fTrackMult);
 646   }
 647   else if(fAOD){
 648     fTrackMult = Ntracks;
 649     if(fDebug)
 650       printf("AOD Track mult= %d\n",fTrackMult);
 651   }
 652   fTrMultDist->Fill(fTrackMult);
 653   TList *l = 0;
 654   TString clusArrayName = "";
 655   if(fESD){
 656     l = fESD->GetList();
 657     if(fDebug)
 658       l->Print();
 659     for(int nk=0;nk<l->GetEntries();nk++){
 660       TObject *obj = (TObject*)l->At(nk);
 661       TString oname = obj->GetName();
 662       if(oname.Contains("CaloClus"))
 663         clusArrayName = oname;
 664       else
 665         continue;
 666       if(clusArrayName=="CaloClusters")
 667         fClusArrayNames->Fill(0);
 668       else{
 669         if(clusArrayName=="EmcCaloClusters")
 670           fClusArrayNames->Fill(1);
 671         else
 672           fClusArrayNames->Fill(2);
 673       }
 674     }
 675     fESDClusters =  dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName));
 676     fESDCells = fESD->GetEMCALCells();
 677     if(fDebug)
 678       printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast());
 679   }
 680   else if(fAOD){
 681     l = fAOD->GetList();
 682     if(fDebug)
 683       l->Print();
 684     //fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
 685     for(int nk=0;nk<l->GetEntries();nk++){
 686       TObject *obj = (TObject*)l->At(nk);
 687       TString oname = obj->GetName();
 688       if(oname.Contains("aloClus"))
 689         clusArrayName = oname;
 690       else
 691         continue;
 692       if(clusArrayName=="caloClusters")
 693         fClusArrayNames->Fill(0);
 694       else{
 695         if(clusArrayName=="EmcCaloClusters")
 696           fClusArrayNames->Fill(1);
 697         else
 698           fClusArrayNames->Fill(2);
 699       }
 700     }
 701     fAODClusters = dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName));
 702     fAODCells = fAOD->GetEMCALCells();
 703     if(fDebug)
 704       printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast());
 705   }
 706   if(fDebug){
 707     printf("clus array is named %s +++++++++\n",clusArrayName.Data());
 708   }
 709
 710
 711   fMCEvent = MCEvent();
 712   if(fMCEvent){
 713     if(fDebug)
 714       std::cout<<"MCevent exists\n";
 715     fStack = (AliStack*)fMCEvent->Stack();
 716     if(!fStack)
 717       fAODMCParticles = (TClonesArray*)fVEvent->FindListObject("mcparticles");
 718   }
 719   else{
 720     if(fDebug)
 721       std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
 722   }
 723   LoopOnCells();
 724   FollowGamma();
 725   if(fDebug)
 726     printf("passed calling of FollowGamma\n");
 727   FillClusHists();
 728   if(fDebug)
 729     printf("passed calling of FillClusHists\n");
 730   FillMcHists();
 731   if(fDebug)
 732     printf("passed calling of FillMcHists\n");
 733   if(fFillQA)
 734     FillQA();
 735   if(fDebug)
 736     printf("passed calling of FillQA\n");
 737   if(fESD)
 738     fESDClusters->Clear();
 739   fSelPrimTracks->Clear();
 740   fNClusForDirPho = 0;
 741   fNCells50 = 0;
 742   fClusIdFromTracks = "";
 743   fVecPv.Clear();
 744
 745   PostData(1, fOutputList);
 746   PostData(2, fQAList);
 747 }
 748
 749 //________________________________________________________________________
 750 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 751 {
 752   if(fDebug)
 753     printf("Inside FillClusHists()....\n");
 754   // Fill cluster histograms.
 755   TObjArray *clusters = fESDClusters;
 756
 757   if (!clusters){
 758     clusters = fAODClusters;
 759     if(fDebug)
 760       printf("ESD clusters empty...");
 761   }
 762   if (!clusters){
 763     if(fDebug)
 764       printf("  and AOD clusters as well!!!\n");
 765     return;
 766   }
 767   if(fDebug)
 768     printf("\n");
 769
 770   const Int_t nclus = clusters->GetEntries();
 771   if(nclus==0)
 772     return;
 773   if(fDebug)
 774     printf("Inside FillClusHists and there are %d clusters\n",nclus);
 775   Double_t maxE = 0;
 776   Int_t nclus10 = 0;
 777   Double_t ptmc=-1;
 778   for(Int_t ic=0;ic<nclus;ic++){
 779     maxE=0;
 780     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 781     if(!c){
 782       if(fDebug)
 783         printf("cluster pointer does not exist! xxxx\n");
 784       continue;
 785     }
 786     if(!c->IsEMCAL()){
 787       if(fDebug)
 788         printf("cluster is not EMCAL! xxxx\n");
 789       continue;
 790     }
 791     if(c->E()<fECut){
 792       if(fDebug)
 793         printf("cluster has E<%1.1f! xxxx\n", fECut);
 794       continue;
 795     }
 796     if(fCpvFromTrack && fClusIdFromTracks.Contains(Form("%d",ic))){
 797       if(fDebug)
 798         printf("cluster does not pass CPV criterion! xxxx\n");
 799        continue;
 800     }
 801     if(IsExotic(c)){
 802       if(fDebug)
 803         printf("cluster is exotic! xxxx\n");
 804       continue;
 805     }
 806     if(c->GetDistanceToBadChannel()<fDistToBadChan){
 807       if(fDebug)
 808         printf("cluster distance to bad channel is %1.1f (<%1.1f) xxxx\n",c->GetDistanceToBadChannel(),fDistToBadChan);
 809       continue;
 810     }
 811     Short_t id;
 812     Double_t Emax = GetMaxCellEnergy( c, id);
 813     if(fDebug)
 814       printf("cluster max cell E=%1.1f",Emax);
 815     Float_t clsPos[3] = {0,0,0};
 816     c->GetPosition(clsPos);
 817     TVector3 clsVec(clsPos);
 818     clsVec -= fVecPv;
 819     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 820     if(fDebug)
 821       printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
 822     if(Et>10)
 823       nclus10++;
 824     Float_t ceiso=0, cephiband=0, cecore=0;
 825     Float_t triso=0, trphiband=0, trcore=0;
 826     Float_t alliso=0, allphiband=0;//, allcore;
 827     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 828     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 829     Bool_t isCPV = kFALSE;
 830     if(TMath::Abs(c->GetTrackDx())>0.03 || TMath::Abs(c->GetTrackDz())>0.02)
 831       isCPV = kTRUE;
 832     GetCeIso(clsVec, id, ceiso, cephiband, cecore, Et);
 833     GetTrIso(clsVec, triso, trphiband, trcore);
 834     Int_t nInConePairs = 0;
 835     Double_t onePairMass = 0;
 836     //---
 837     //if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV){
 838     TObjArray *inConeInvMassArr = (TObjArray*)fInConeInvMass.Tokenize(";");
 839     TObjArray *inConePairClEt =  (TObjArray*)fInConePairClEt.Tokenize(";");
 840     nInConePairs = inConeInvMassArr->GetEntriesFast();
 841     Int_t nInConePi0 = inConePairClEt->GetEntriesFast();
 842     if(nInConePairs != nInConePi0)
 843       printf("Inconsistent number of in cone pairs!!!\n");
 844     for(int ipair=0;ipair<nInConePairs;ipair++){
 845       TObjString *obs = (TObjString*)inConeInvMassArr->At(ipair);
 846       TObjString *obet = (TObjString*)inConePairClEt->At(ipair);
 847       TString smass = obs->GetString();
 848       TString spairEt = obet->GetString();
 849       Double_t pairmass = smass.Atof();
 850       Double_t pairEt = spairEt.Atof();//this must be zero when inv mass outside pi0 range
 851       if(0==ipair && nInConePairs==1)
 852         onePairMass = pairmass;
 853       if(fDebug)
 854         printf("=================+++++++++++++++Inv mass inside the cone for photon range: %1.1f,%1.1f,%1.1f+-++++-+-+-+-++-+-+-\n",Et,pairmass,ceiso+triso);
 855       fEtCandIsoAndIsoWoPairEt->Fill(Et,ceiso+triso,ceiso+triso-pairEt);
 856     }
 857     //}
 858     //---
 859     Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
 860     if(Et>10 && Et<15 && dr>0.025){
 861       fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
 862       if(fDebug)
 863         printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
 864     }
 865     alliso = ceiso + triso;
 866     allphiband = cephiband + trphiband;
 867     //allcore = cecore + trcore;
 868     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 869     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 870     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 871     Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR);
 872     if(fDebug && Et>10)
 873       printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
 874     if(c->GetM02()>0.5 && c->GetM02()<2.0){
 875       fMcPtInConeBG->Fill(alliso-allisoue, mcptsum);
 876       fMcPtInConeBGnoUE->Fill(alliso, mcptsum);
 877       fMcPtInConeTrBGnoUE->Fill(triso, mcptsum);
 878     }
 879     if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){
 880       fMcPtInConeSBG->Fill(alliso-allisoue, mcptsum);
 881       fMcPtInConeSBGnoUE->Fill(alliso, mcptsum);
 882       fMcPtInConeTrSBGnoUE->Fill(triso, mcptsum);
 883       if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){
 884         fAllIsoEtMcGamma->Fill(Et, alliso-cecore-allisoue);
 885         fAllIsoNoUeEtMcGamma->Fill(Et, alliso-cecore);
 886       }
 887     }
 888     if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV)
 889       fClusEtCPVSBGISO->Fill(Et,alliso - trcore);
 890     if(c->GetM02()>0.5 && c->GetM02()<2.0 && isCPV)
 891       fClusEtCPVBGISO->Fill(Et,alliso - trcore);
 892     const Int_t ndims =   fNDimensions;
 893     Double_t outputValues[ndims];
 894     if(mcptsum<2)
 895       ptmc = GetClusSource(c);
 896     else
 897       ptmc = -0.1;
 898     outputValues[0] = Et;
 899     outputValues[1] = c->GetM02();
 900     outputValues[2] = ceiso/*cecore*/-ceisoue;
 901     outputValues[3] = triso-trisoue;
 902     outputValues[4] = alliso/*cecore*/-allisoue - trcore;
 903     outputValues[5] = ceiso;
 904     outputValues[6] = alliso - trcore;
 905     if(fDebug)
 906       printf("track-cluster dphi=%1.3f, deta=%1.3f\n",c->GetTrackDx(),c->GetTrackDz());
 907     if(TMath::Abs(c->GetTrackDx())<0.1)
 908       outputValues[7] = c->GetTrackDx();
 909     else
 910       outputValues[7] = 0.099*c->GetTrackDx()/TMath::Abs(c->GetTrackDx());
 911     if(TMath::Abs(c->GetTrackDz())<0.05)
 912       outputValues[8] = c->GetTrackDz();
 913     else
 914       outputValues[8] = 0.049*c->GetTrackDz()/TMath::Abs(c->GetTrackDz());
 915     outputValues[9] = clsVec.Eta();
 916     outputValues[10] = clsVec.Phi();
 917     if(fESDCells)
 918       outputValues[11] = fESDCells->GetCellTime(id);
 919     else if(fAODCells)
 920       outputValues[11] = fAODCells->GetCellTime(id);
 921     outputValues[12] = fTrackMult;
 922     outputValues[13] = ptmc;
 923     if(nInConePairs == 1)
 924       outputValues[14] = onePairMass;
 925     else
 926       outputValues[14] = -1;
 927     fHnOutput->Fill(outputValues);
 928     if(c->E()>maxE)
 929       maxE = c->E();
 930   }
 931   fNClusHighClusE->Fill(maxE,nclus);
 932   fMaxEClus = maxE;
 933   fNClusEt10->Fill(nclus10);
 934   fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
 935 }
 936
 937 //________________________________________________________________________
 938 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core, Double_t EtCl)
 939 {
 940   if(fDebug)
 941     printf("....indside GetCeIso funtcion\n");
 942   // Get cell isolation.
 943   AliVCaloCells *cells = fESDCells;
 944   if (!cells){
 945     cells = fAODCells;
 946     if(fDebug)
 947       printf("ESD cells empty...\n");
 948   }
 949   if (!cells){
 950      if(fDebug)
 951       printf("  and AOD cells are empty  as well!!!\n");
 952     return;
 953   }
 954
 955   TObjArray *clusters = fESDClusters;
 956   if (!clusters)
 957     clusters = fAODClusters;
 958   if (!clusters)
 959     return;
 960
 961   fInConeInvMass = "";
 962   fInConePairClEt="";
 963   const Int_t nclus = clusters->GetEntries();
 964   //const Int_t ncells = cells->GetNumberOfCells();
 965   Float_t totiso=0;
 966   Float_t totband=0;
 967   Float_t totcore=0;
 968   Float_t etacl = vec.Eta();
 969   Float_t phicl = vec.Phi();
 970   if(phicl<0)
 971     phicl+=TMath::TwoPi();
 972   /*Int_t absid = -1;
 973   Float_t eta=-1, phi=-1;
 974   for(int icell=0;icell<ncells;icell++){
 975     absid = TMath::Abs(cells->GetCellNumber(icell));
 976     Float_t celltime = cells->GetCellTime(absid);
 977     //if(TMath::Abs(celltime)>2e-8 && (!fIsMc))
 978     if(TMath::Abs(celltime-maxtcl)>2e-8 )
 979       continue;
 980     if(!fGeom)
 981       return;
 982     fGeom->EtaPhiFromIndex(absid,eta,phi);
 983     Float_t dphi = TMath::Abs(phi-phicl);
 984     Float_t deta = TMath::Abs(eta-etacl);
 985     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/
 986   for(int ic=0;ic<nclus;ic++){
 987     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 988     if(!c)
 989       continue;
 990     if(!c->IsEMCAL())
 991       continue;
 992     if(c->E()<fMinIsoClusE)
 993       continue;
 994     Short_t id=-1;
 995     GetMaxCellEnergy( c, id);
 996     Double_t maxct = cells->GetCellTime(id);
 997     if(TMath::Abs(maxct)>fClusTDiff/*2.5e-9*/ && (!fIsMc))
 998       continue;
 999     Float_t clsPos[3] = {0,0,0};
1000     c->GetPosition(clsPos);
1001     TVector3 cv(clsPos);
1002     cv -= fVecPv;
1003     Double_t Et = c->E()*TMath::Sin(cv.Theta());
1004     Float_t dphi = (cv.Phi()-phicl);
1005     Float_t deta = (cv.Eta()-etacl);
1006     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
1007     if(R<0.007)
1008       continue;
1009     if(maxid==id)
1010       continue;
1011     Double_t matchedpt =  GetTrackMatchedPt(c->GetTrackMatchedIndex());
1012     if(fCpvFromTrack){
1013       if(matchedpt>0 && fRemMatchClus )
1014         continue;
1015     } else {
1016       if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
1017         if(fRemMatchClus){
1018           if(fDebug)
1019             printf("This isolation cluster is matched to a track!++++++++++++++++++++++++++++++++++++++++++++++++++\n");
1020           continue;
1021         }
1022       }
1023     }
1024     Double_t nEt = TMath::Max(Et-matchedpt, 0.0);
1025     if(nEt<0)
1026       printf("nEt=%1.1f\n",nEt);
1027     if(R<fIsoConeR){
1028       if(c->GetM02()>0.1 && c->GetM02()<0.3 && !(matchedpt>0)){
1029         TLorentzVector lv, lvec;
1030         lv.SetPtEtaPhiM(Et,cv.Eta(),cv.Phi(),0);
1031         lvec.SetPtEtaPhiM(EtCl,vec.Eta(),vec.Phi(),0);
1032         TLorentzVector lpair = lv + lvec;
1033         fInConeInvMass += Form("%f;",lpair.M());
1034         if(lpair.M()>0.11 && lpair.M()<0.165){
1035           fInConePairedClusEtVsCandEt->Fill(EtCl,Et);
1036           fInConePairClEt += Form("%f;",Et);
1037           //continue;
1038         }
1039         else
1040           fInConePairClEt += Form("%f;",0.0);
1041         /*if(lpair.M()>0.52 && lpair.M()<0.58)
1042           continue;*/
1043       }
1044       totiso += nEt;
1045       if(R<0.04)
1046         totcore += nEt;
1047     }
1048     else{
1049       if(dphi>fIsoConeR)
1050         continue;
1051       if(deta<fIsoConeR)
1052         continue;
1053       totband += nEt;
1054     }
1055   }
1056   iso = totiso;
1057   phiband = totband;
1058   core = totcore;
1059 }
1060 //________________________________________________________________________
1061 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
1062 {
1063   // Get track isolation.
1064
1065   if(!fSelPrimTracks)
1066     return;
1067   fHigherPtCone = 0;
1068   const Int_t ntracks = fSelPrimTracks->GetEntries();
1069   Double_t totiso=0;
1070   Double_t totband=0;
1071   Double_t totcore=0;
1072   Double_t etacl = vec.Eta();
1073   Double_t phicl = vec.Phi();
1074   if(phicl<0)
1075     phicl+=TMath::TwoPi();
1076   for(int itrack=0;itrack<ntracks;itrack++){
1077     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
1078     if(!track)
1079       continue;
1080     Double_t dphi = TMath::Abs(track->Phi()-phicl);
1081     Double_t deta = TMath::Abs(track->Eta()-etacl);
1082     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
1083     Double_t pt = track->Pt();
1084     if(pt>fHigherPtCone)
1085       fHigherPtCone = pt;
1086     if(R<fIsoConeR){
1087       totiso += track->Pt();
1088       if(R<0.04 && this->fTrCoreRem)
1089         totcore += pt;
1090     }
1091     else{
1092       if(dphi>fIsoConeR)
1093         continue;
1094       if(deta<fIsoConeR)
1095         continue;
1096       totband += track->Pt();
1097     }
1098   }
1099   iso = totiso;
1100   phiband = totband;
1101   core = totcore;
1102 }
1103
1104 //________________________________________________________________________
1105 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1106 {
1107   // Calculate the energy of cross cells around the leading cell.
1108
1109   AliVCaloCells *cells = 0;
1110   cells = fESDCells;
1111   if (!cells)
1112     cells = fAODCells;
1113   if (!cells)
1114     return 0;
1115
1116   if (!fGeom)
1117     return 0;
1118
1119   Int_t iSupMod = -1;
1120   Int_t iTower  = -1;
1121   Int_t iIphi   = -1;
1122   Int_t iIeta   = -1;
1123   Int_t iphi    = -1;
1124   Int_t ieta    = -1;
1125   Int_t iphis   = -1;
1126   Int_t ietas   = -1;
1127
1128   Double_t crossEnergy = 0;
1129
1130   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1131   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1132
1133   Int_t ncells = cluster->GetNCells();
1134   for (Int_t i=0; i<ncells; i++) {
1135     Int_t cellAbsId = cluster->GetCellAbsId(i);
1136     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1137     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1138     Int_t aphidiff = TMath::Abs(iphi-iphis);
1139     if (aphidiff>1)
1140       continue;
1141     Int_t aetadiff = TMath::Abs(ieta-ietas);
1142     if (aetadiff>1)
1143       continue;
1144     if ( (aphidiff==1 && aetadiff==0) ||
1145         (aphidiff==0 && aetadiff==1) ) {
1146       crossEnergy += cells->GetCellAmplitude(cellAbsId);
1147     }
1148   }
1149
1150   return crossEnergy;
1151 }
1152
1153 //________________________________________________________________________
1154 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1155 {
1156   // Get maximum energy of attached cell.
1157
1158   id = -1;
1159
1160   AliVCaloCells *cells = 0;
1161   cells = fESDCells;
1162   if (!cells)
1163     cells = fAODCells;
1164   if(!cells)
1165     return 0;
1166
1167   Double_t maxe = 0;
1168   Int_t ncells = cluster->GetNCells();
1169   for (Int_t i=0; i<ncells; i++) {
1170     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1171     if (e>maxe) {
1172       maxe = e;
1173       id   = cluster->GetCellAbsId(i);
1174     }
1175   }
1176   return maxe;
1177 }
1178
1179 //________________________________________________________________________
1180 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
1181 {
1182   if(!fStack && !fAODMCParticles)
1183     return;
1184   //ESD
1185   if(fStack){
1186     Int_t nTracks = fStack->GetNtrack();
1187     if(fDebug)
1188       printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1189     for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1190       TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1191       if(!mcp)
1192         continue;
1193       Int_t pdg = mcp->GetPdgCode();
1194       if(pdg!=22)
1195         continue;
1196       if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1197         continue;
1198       Int_t imom = mcp->GetMother(0);
1199       if(imom<0 || imom>nTracks)
1200         continue;
1201       TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
1202       if(!mcmom)
1203         continue;
1204       Int_t pdgMom = mcmom->GetPdgCode();
1205       Double_t mcphi = mcp->Phi();
1206       Double_t mceta = mcp->Eta();
1207       if((imom==6 || imom==7) && pdgMom==22) {
1208         fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1209         Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1210         fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1211         if(ptsum<2)
1212           fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1213         if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1214           fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1215         if(fNClusForDirPho==0)
1216           fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1217         if(fDebug){
1218           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());
1219           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());
1220         }
1221       }
1222       else{
1223         if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1224           fDecayPhotonPtMC->Fill(mcp->Pt());
1225       }
1226     }
1227   }
1228   //AOD
1229   else if(fAODMCParticles){
1230     Int_t nTracks = fAODMCParticles->GetEntriesFast();
1231     if(fDebug)
1232       printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
1233     for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
1234       AliAODMCParticle *mcp = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1235       if(!mcp)
1236         continue;
1237       Int_t pdg = mcp->GetPdgCode();
1238       if(pdg!=22)
1239         continue;
1240       if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
1241         continue;
1242       Int_t imom = mcp->GetMother();
1243       if(imom<0 || imom>nTracks)
1244         continue;
1245       AliAODMCParticle *mcmom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(imom));
1246       if(!mcmom)
1247         continue;
1248       Int_t pdgMom = mcmom->GetPdgCode();
1249       Double_t mcphi = mcp->Phi();
1250       Double_t mceta = mcp->Eta();
1251       if((imom==6 || imom==7) && pdgMom==22) {
1252         fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1253         Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR);
1254         fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum);
1255         if(ptsum<2)
1256           fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1257         if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR))
1258           fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum);
1259         if(fNClusForDirPho==0)
1260           fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
1261         if(fDebug){
1262           printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatus());
1263           printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatus());
1264         }
1265       }
1266       else{
1267         if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
1268           fDecayPhotonPtMC->Fill(mcp->Pt());
1269       }
1270     }
1271   }
1272 }
1273 //________________________________________________________________________
1274 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
1275 {
1276   if(!c)
1277     return -0.1;
1278   if(!fStack && !fAODMCParticles)
1279     return -0.1;
1280   Int_t clabel = c->GetLabel();
1281   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
1282     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
1283   if(!fMcIdFamily.Contains(Form("%d",clabel)))
1284     return -0.1;
1285   fNClusForDirPho++;
1286   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
1287   Int_t partonpos = partonposstr.Atoi();
1288   if(fDebug)
1289     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
1290   Float_t clsPos[3] = {0,0,0};
1291   c->GetPosition(clsPos);
1292   TVector3 clsVec(clsPos);
1293   clsVec -= fVecPv;
1294   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
1295   //ESD
1296   if(fStack){
1297     Int_t nmcp = fStack->GetNtrack();
1298     if(clabel<0 || clabel>nmcp)
1299       return -0.1;
1300     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
1301     if(!mcp)
1302       return -0.1;
1303     if(fDebug){
1304       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);
1305     }
1306     Int_t lab1 =  mcp->GetFirstDaughter();
1307     if(lab1<0 || lab1>nmcp)
1308       return -0.1;
1309     TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
1310     if(!mcd)
1311       return -0.1;
1312     if(fDebug)
1313       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);
1314     if(mcd->GetPdgCode()==22){
1315       fClusEtMcPt->Fill(mcd->Pt(), Et);
1316       fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1317     }
1318     else{
1319       if(fDebug)
1320         printf("Warning: daughter of photon parton is not a photon\n");
1321       return -0.1;
1322     }
1323   }
1324   //AOD
1325   else if(fAODMCParticles){
1326     Int_t nmcp = fAODMCParticles->GetEntriesFast();
1327     if(clabel<0 || clabel>nmcp)
1328       return -0.1;
1329     AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(partonpos));
1330     if(!mcp)
1331       return -0.1;
1332     if(fDebug){
1333       printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),clabel);
1334     }
1335     Int_t lab1 =  mcp->GetDaughter(0);
1336     if(lab1<0 || lab1>nmcp)
1337       return -0.1;
1338     AliAODMCParticle *mcd = static_cast<AliAODMCParticle*>(fAODMCParticles->At(lab1));
1339     if(!mcd)
1340       return -0.1;
1341     if(fDebug)
1342       printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->E(),mcd->GetPdgCode(),lab1);
1343     if(mcd->GetPdgCode()==22){
1344       fClusEtMcPt->Fill(mcd->Pt(), Et);
1345       fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
1346     }
1347     else{
1348       printf("Warning: daughter of photon parton is not a photon\n");
1349       return -0.1;
1350     }
1351   }
1352   return fDirPhoPt;
1353 }
1354 //________________________________________________________________________
1355 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
1356 {
1357   if(!fStack && !fAODMCParticles)
1358     return;
1359   Int_t selfid = 6;
1360   //ESD
1361   if(fStack){
1362     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
1363     if(!mcp)
1364       return;
1365     if(mcp->GetPdgCode()!=22){
1366       mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
1367       if(!mcp)
1368         return;
1369     }
1370     Int_t daug0f =  mcp->GetFirstDaughter();
1371     Int_t daug0l =  mcp->GetLastDaughter();
1372     Int_t nd0 = daug0l - daug0f;
1373     if(fDebug)
1374       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);
1375     fMcIdFamily = Form("%d,",selfid);
1376     GetDaughtersInfo(daug0f,daug0l, selfid,"");
1377     if(fDebug){
1378       printf("\t---- end of the prompt  gamma's family tree ----\n\n");
1379       printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1380     }
1381     TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
1382     if(!md)
1383       return;
1384     fDirPhoPt = md->Pt();
1385   }
1386   //AOD
1387   else   if(fAODMCParticles){
1388     AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(selfid));
1389     if(!mcp)
1390       return;
1391     if(mcp->GetPdgCode()!=22){
1392       mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(++selfid));
1393       if(!mcp)
1394         return;
1395     }
1396     Int_t daug0f =  mcp->GetDaughter(0);
1397     Int_t daug0l =  mcp->GetDaughter(mcp->GetNDaughters()-1);
1398     Int_t nd0 = daug0l - daug0f;
1399     if(fDebug)
1400       printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),nd0+1);
1401     fMcIdFamily = Form("%d,",selfid);
1402     GetDaughtersInfo(daug0f,daug0l, selfid,"");
1403     if(fDebug){
1404       printf("\t---- end of the prompt  gamma's family tree ----\n\n");
1405       printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1406     }
1407     AliAODMCParticle *md = static_cast<AliAODMCParticle*>(fAODMCParticles->At(daug0f));
1408     if(!md)
1409       return;
1410     fDirPhoPt = md->Pt();
1411   }
1412
1413 }
1414 //________________________________________________________________________
1415 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
1416 {
1417   if(fStack){
1418     int nmcp = fStack->GetNtrack();
1419     if(firstd<0 || firstd>nmcp)
1420       return;
1421     if(lastd<0 || lastd>nmcp)
1422       return;
1423     if(firstd>lastd){
1424       printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1425     return;
1426     }
1427     TString indenter = Form("\t%s",inputind);
1428     TParticle *dp = 0x0;
1429     if(fDebug)
1430       printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1431     for(int id=firstd; id<lastd+1; id++){
1432       dp =  static_cast<TParticle*>(fStack->Particle(id));
1433       if(!dp)
1434         continue;
1435       Int_t dfd = dp->GetFirstDaughter();
1436       Int_t dld = dp->GetLastDaughter();
1437       Int_t dnd =  dld - dfd + 1;
1438       Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
1439       if(fDebug)
1440         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);
1441       fMcIdFamily += Form("%d,",id);
1442       GetDaughtersInfo(dfd,dld,id,indenter.Data());
1443     }
1444   }
1445   if(fAODMCParticles){
1446     int nmcp = fAODMCParticles->GetEntriesFast();
1447     if(firstd<0 || firstd>nmcp)
1448       return;
1449     if(lastd<0 || lastd>nmcp)
1450       return;
1451     if(firstd>lastd){
1452       printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1453     return;
1454     }
1455     TString indenter = Form("\t%s",inputind);
1456     AliAODMCParticle *dp = 0x0;
1457     if(fDebug)
1458       printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1459     for(int id=firstd; id<lastd+1; id++){
1460       dp =  static_cast<AliAODMCParticle*>(fAODMCParticles->At(id));
1461       if(!dp)
1462         continue;
1463       Int_t dfd = dp->GetDaughter(0);
1464       Int_t dld = dp->GetDaughter(dp->GetNDaughters()-1);
1465       Int_t dnd =  dld - dfd + 1;
1466       Float_t vr = TMath::Sqrt(dp->Xv()*dp->Xv()+dp->Xv()*dp->Xv());
1467       if(fDebug)
1468         printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->E(), vr, dp->GetPdgCode(), dnd, dfd, dld);
1469       fMcIdFamily += Form("%d,",id);
1470       GetDaughtersInfo(dfd,dld,id,indenter.Data());
1471     }
1472   }
1473 }
1474
1475 //________________________________________________________________________
1476 Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
1477 {
1478   if(!fStack && !fAODMCParticles)
1479     return 0;
1480   if(fDebug)
1481     printf("Inside GetMcPtSumInCone!!\n");
1482   Float_t ptsum = 0;
1483   TString addpartlabels = "";
1484   //ESD
1485   if(fStack){
1486     Int_t nTracks = fStack->GetNtrack();
1487     for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1488       TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1489       if(!mcp)
1490         continue;
1491       Int_t status = mcp->GetStatusCode();
1492       if(status!=1)
1493         continue;
1494       Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz());
1495       if(mcvr>1)
1496         continue;
1497       /*else {
1498         if(fDebug)
1499         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1500         }*/
1501       Float_t dphi = mcp->Phi() - phiclus;
1502       Float_t deta = mcp->Eta() - etaclus;
1503       if(fDebug && TMath::Abs(dphi)<0.01)
1504         printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1505
1506       if(deta>R || dphi>R)
1507         continue;
1508       Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1509       if(dR>R)
1510         continue;
1511       addpartlabels += Form("%d,",iTrack);
1512       if(mcp->Pt()<0.2)
1513         continue;
1514       ptsum += mcp->Pt();
1515     }
1516   }
1517   //AOD
1518   if(fAODMCParticles){
1519     Int_t nTracks = fAODMCParticles->GetEntriesFast();
1520     for(Int_t iTrack=9;iTrack<nTracks;iTrack++){
1521       AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
1522       if(!mcp)
1523         continue;
1524       Int_t status = mcp->GetStatus();
1525       if(status!=1)
1526         continue;
1527       Float_t mcvr = TMath::Sqrt(mcp->Xv()*mcp->Xv()+ mcp->Yv()* mcp->Yv() + mcp->Zv()*mcp->Zv());
1528       if(mcvr>1)
1529         continue;
1530       /*else {
1531         if(fDebug)
1532         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1533         }*/
1534       Float_t dphi = mcp->Phi() - phiclus;
1535       Float_t deta = mcp->Eta() - etaclus;
1536       if(fDebug && TMath::Abs(dphi)<0.01)
1537         printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1538
1539       if(deta>R || dphi>R)
1540         continue;
1541       Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1542       if(dR>R)
1543         continue;
1544       addpartlabels += Form("%d,",iTrack);
1545       if(mcp->Pt()<0.2)
1546         continue;
1547       ptsum += mcp->Pt();
1548     }
1549   }
1550   return ptsum;
1551 }
1552 //________________________________________________________________________
1553 void AliAnalysisTaskEMCALIsoPhoton::FillQA()
1554 {
1555
1556   TObjArray *clusters = fESDClusters;
1557   //"none", "exotic", "exo+cpv1", "exo+cpv1+time", "exo+cpv1+time+m02"),
1558   if (!clusters){
1559     clusters = fAODClusters;
1560     if(fDebug)
1561       printf("ESD clusters empty...");
1562   }
1563   if (!clusters){
1564     if(fDebug)
1565       printf("  and AOD clusters as well!!!\n");
1566     return;
1567   }
1568   if(!fSelPrimTracks)
1569     return;
1570   const int ntracks = fSelPrimTracks->GetEntriesFast();
1571   const int ncells = fNCells50;//fESDCells->GetNumberOfCells();
1572   const Int_t nclus = clusters->GetEntries();
1573   fNTracks->Fill(ntracks);
1574   fEmcNCells->Fill(ncells);
1575   fEmcNClus->Fill(nclus);
1576   if(fMaxEClus>fECut){
1577     fNTracksECut->Fill(ntracks);
1578     fEmcNCellsCut->Fill(ncells);
1579     fEmcNClusCut->Fill(nclus);
1580   }
1581   for(int it=0;it<ntracks;it++){
1582     AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it);
1583     if(!t)
1584       continue;
1585     fTrackPtPhi->Fill(t->Pt(),t->Phi());
1586     fTrackPtEta->Fill(t->Pt(),t->Eta());
1587     if(fMaxEClus>fECut){
1588       fTrackPtPhiCut->Fill(t->Pt(), t->Phi());
1589       fTrackPtEtaCut->Fill(t->Pt(), t->Eta());
1590     }
1591     if(t->GetTPCsignal()<80 || t->GetTPCsignal()>100)
1592       continue;
1593     if(t->GetEMCALcluster()<=0 || t->GetEMCALcluster()>nclus)
1594       continue;
1595     AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(t->GetEMCALcluster()));
1596     if(!c)
1597       continue;
1598     fEoverPvsE->Fill(c->E(),c->E()/t->P());
1599   }
1600   for(int ic=0;ic<nclus;ic++){
1601     AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(ic));
1602     //AliESDCaloCluster *c = (AliESDCaloCluster*)clusters->At(ic);
1603     if(!c)
1604       continue;
1605     if(!c->IsEMCAL())
1606       continue;
1607     Float_t clsPos[3] = {0,0,0};
1608     c->GetPosition(clsPos);
1609     TVector3 clsVec(clsPos);
1610     clsVec -= fVecPv;
1611     Double_t cphi = clsVec.Phi();
1612     Double_t ceta = clsVec.Eta();
1613     Short_t id = -1;
1614     GetMaxCellEnergy( c, id);
1615     fEmcClusEClusCuts->Fill(c->E(),0);
1616     fEmcClusEPhi->Fill(c->E(), cphi);
1617     fEmcClusEEta->Fill(c->E(), ceta);
1618     if(fMaxEClus>fECut){
1619       fEmcClusEPhiCut->Fill(c->E(), cphi);
1620       fEmcClusEEtaCut->Fill(c->E(), ceta);
1621     }
1622     Double_t maxt=0;
1623     if(fESDCells)
1624       maxt = fESDCells->GetCellTime(id);
1625     else if(fAODCells)
1626       maxt = fAODCells->GetCellTime(id);
1627     if(IsExotic(c))
1628       continue;
1629     fEmcClusNotExo->Fill(c->E());
1630     fEmcClusEClusCuts->Fill(c->E(),1);
1631     if(fClusIdFromTracks.Contains(Form("%d",ic))){
1632       fEmcClusETM2->Fill(c->E());
1633       fDetaDphiFromTM->Fill(c->GetTrackDz(),c->GetTrackDx());
1634     }
1635     if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){
1636       fEmcClusETM1->Fill(c->E());
1637       continue;
1638     }
1639     fEmcClusEClusCuts->Fill(c->E(),2);
1640     if(TMath::Abs(maxt)>30e-9)
1641       continue;
1642     fEmcClusEClusCuts->Fill(c->E(),3);
1643     if(c->GetM02()>0.1)
1644       fEmcClusEClusCuts->Fill(c->E(),4);
1645   }
1646 }
1647 //________________________________________________________________________
1648 Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex)
1649 {
1650   Double_t pt = 0;
1651   if(!fTracks)
1652     return pt;
1653   if(matchIndex<0 || matchIndex>fTracks->GetEntries()){
1654     if(fDebug)
1655       printf("track-matched index out of track array range!!!\n");
1656     return pt;
1657   }
1658   AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex));
1659   if(!track){
1660     if(fDebug)
1661       printf("track-matched pointer does not exist!!!\n");
1662     return pt;
1663   }
1664   if(fESD){
1665     if(fPrTrCuts && fPrTrCuts->IsSelected(track))
1666       pt = track->Pt();
1667     else {
1668       if(fCompTrCuts && fCompTrCuts->IsSelected(track))
1669         pt = track->Pt();
1670       else
1671         return pt;
1672     }
1673   }
1674   return pt;
1675 }
1676 //________________________________________________________________________
1677 void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells()
1678 {
1679   AliVCaloCells *cells = 0;
1680   cells = fESDCells;
1681   if (!cells)
1682     cells = fAODCells;
1683   if(!cells)
1684     return;
1685   Double_t maxe = 0;
1686   Double_t maxphi = -10;
1687   Int_t ncells = cells->GetNumberOfCells();
1688   Double_t eta,phi;
1689   for (Int_t i=0; i<ncells; i++) {
1690     Short_t absid = TMath::Abs(cells->GetCellNumber(i));
1691     Double_t e = cells->GetCellAmplitude(absid);
1692     if(e>0.05)
1693       fNCells50++;
1694     else
1695       continue;
1696     fGeom->EtaPhiFromIndex(absid,eta,phi);
1697     if(maxe<e){
1698       maxe = e;
1699       maxphi = phi;
1700     }
1701   }
1702   fMaxCellEPhi->Fill(maxe,maxphi);
1703
1704 }
1705 //________________________________________________________________________
1706 bool AliAnalysisTaskEMCALIsoPhoton::IsExotic(AliVCluster *c)
1707 {
1708   bool isExo = 0;
1709   Short_t id = -1;
1710   Double_t Emax = GetMaxCellEnergy( c, id);
1711   Double_t Ecross = GetCrossEnergy( c, id);
1712   if((1-Ecross/Emax)>fExoticCut)
1713     isExo = 1;
1714   return isExo;
1715 }
1716 //________________________________________________________________________
1717 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
1718 {
1719   // Called once at the end of the query.
1720 }