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