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