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