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 <TLorentzVector.h>
  11 #include <TTree.h>
  12
  13 #include "AliAnalysisManager.h"
  14 #include "AliAnalysisTask.h"
  15 #include "AliEMCALGeometry.h"
  16 #include "AliEMCALRecoUtils.h"
  17 #include "AliESDCaloCells.h"
  18 #include "AliESDCaloCluster.h"
  19 #include "AliESDEvent.h"
  20 #include "AliESDHeader.h"
  21 #include "AliESDInputHandler.h"
  22 #include "AliESDUtils.h"
  23 #include "AliESDpid.h"
  24 #include "AliESDtrack.h"
  25 #include "AliESDtrackCuts.h"
  26 #include "AliESDv0.h"
  27 #include "AliKFParticle.h"
  28 #include "AliMCEvent.h"
  29 #include "AliMCEventHandler.h"
  30 #include "AliStack.h"
  31 #include "AliV0vertexer.h"
  32 #include "AliVCluster.h"
  33
  34 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
  35
  36 //________________________________________________________________________
  37 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
  38   AliAnalysisTaskSE(),
  39   fCaloClusters(0),
  40   fSelPrimTracks(0),
  41   fEMCalCells(0),
  42   fPrTrCuts(0),
  43   fGeom(0x0),
  44   fESD(0),
  45   fOutputList(0),
  46   fEvtSel(0),
  47   fPVtxZ(0),
  48   fCellAbsIdVsAmpl(0),
  49   fNClusHighClusE(0),
  50   fM02Et(0),
  51   fM02EtTM(0),
  52   fM02EtCeIso1(0),
  53   fM02EtCeIso2(0),
  54   fM02EtCeIso5(0),
  55   fM02EtTrIso1(0),
  56   fM02EtTrIso2(0),
  57   fM02EtTrIso5(0),
  58   fM02EtAllIso1(0),
  59   fM02EtAllIso2(0),
  60   fM02EtAllIso5(0),
  61   fCeIsoVsEtPho(0),
  62   fTrIsoVsEtPho(0),
  63   fAllIsoVsEtPho(0),
  64   fM02EtCeIso1TM(0),
  65   fM02EtCeIso2TM(0),
  66   fM02EtCeIso5TM(0),
  67   fM02EtTrIso1TM(0),
  68   fM02EtTrIso2TM(0),
  69   fM02EtTrIso5TM(0),
  70   fM02EtAllIso1TM(0),
  71   fM02EtAllIso2TM(0),
  72   fM02EtAllIso5TM(0),
  73   fCeIsoVsEtPhoTM(0),
  74   fTrIsoVsEtPhoTM(0),
  75   fAllIsoVsEtPhoTM(0)
  76 {
  77   // Default constructor.
  78 }
  79
  80 //________________________________________________________________________
  81 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
  82   AliAnalysisTaskSE(name),
  83   fCaloClusters(0),
  84   fSelPrimTracks(0),
  85   fEMCalCells(0),
  86   fPrTrCuts(0),
  87   fGeom(0x0),
  88   fGeoName("EMCAL_COMPLETEV1"),
  89   fPeriod("LHC11c"),
  90   fIsTrain(0),
  91   fExoticCut(0.97),
  92   fIsoConeR(0.4),
  93   fESD(0),
  94   fOutputList(0),
  95   fEvtSel(0),
  96   fPVtxZ(0),
  97   fCellAbsIdVsAmpl(0),
  98   fNClusHighClusE(0),
  99   fM02Et(0),
 100   fM02EtTM(0),
 101   fM02EtCeIso1(0),
 102   fM02EtCeIso2(0),
 103   fM02EtCeIso5(0),
 104   fM02EtTrIso1(0),
 105   fM02EtTrIso2(0),
 106   fM02EtTrIso5(0),
 107   fM02EtAllIso1(0),
 108   fM02EtAllIso2(0),
 109   fM02EtAllIso5(0),
 110   fCeIsoVsEtPho(0),
 111   fTrIsoVsEtPho(0),
 112   fAllIsoVsEtPho(0),
 113   fM02EtCeIso1TM(0),
 114   fM02EtCeIso2TM(0),
 115   fM02EtCeIso5TM(0),
 116   fM02EtTrIso1TM(0),
 117   fM02EtTrIso2TM(0),
 118   fM02EtTrIso5TM(0),
 119   fM02EtAllIso1TM(0),
 120   fM02EtAllIso2TM(0),
 121   fM02EtAllIso5TM(0),
 122   fCeIsoVsEtPhoTM(0),
 123   fTrIsoVsEtPhoTM(0),
 124   fAllIsoVsEtPhoTM(0)
 125 {
 126   // Constructor
 127
 128   // Define input and output slots here
 129   // Input slot #0 works with a TChain
 130   DefineInput(0, TChain::Class());
 131   // Output slot #0 id reserved by the base class for AOD
 132   // Output slot #1 writes into a TH1 container
 133   DefineOutput(1, TList::Class());
 134 }
 135
 136 //________________________________________________________________________
 137 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 138 {
 139   // Create histograms, called once.
 140
 141   fCaloClusters = new TRefArray();
 142   fSelPrimTracks = new TObjArray();
 143
 144
 145   fOutputList = new TList();
 146   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 147
 148   fGeom = AliEMCALGeometry::GetInstance(fGeoName);
 149
 150   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 151   fOutputList->Add(fEvtSel);
 152
 153
 154   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 155   fOutputList->Add(fPVtxZ);
 156
 157   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);
 158   fOutputList->Add(fPVtxZ);
 159
 160   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);
 161   fOutputList->Add(fNClusHighClusE);
 162
 163   fM02Et = new TH2F("fM02Et","M02 vs Et for all clusters;E_{T} (GeV);M02",1000,0,100,400,0,4);
 164   fM02Et->Sumw2();
 165   fOutputList->Add(fM02Et);
 166
 167   fM02EtTM = new TH2F("fM02EtTM","M02 vs Et for all track-matched clusters;E_{T} (GeV);M02",1000,0,100,400,0,4);
 168   fM02EtTM->Sumw2();
 169   fOutputList->Add(fM02EtTM);
 170
 171   fM02EtCeIso1 = new TH2F("fM02EtCeIso1","M02 vs Et for all clusters (ISO_{EMC}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 172   fM02EtCeIso1->Sumw2();
 173   fOutputList->Add(fM02EtCeIso1);
 174
 175   fM02EtCeIso2 = new TH2F("fM02EtCeIso2","M02 vs Et for all clusters (ISO_{EMC}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 176   fM02EtCeIso2->Sumw2();
 177   fOutputList->Add(fM02EtCeIso2);
 178
 179   fM02EtCeIso5 = new TH2F("fM02EtCeIso5","M02 vs Et for all clusters (ISO_{EMC}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 180   fM02EtCeIso5->Sumw2();
 181   fOutputList->Add(fM02EtCeIso5);
 182
 183   fM02EtTrIso1 = new TH2F("fM02EtTrIso1","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 184   fM02EtTrIso1->Sumw2();
 185   fOutputList->Add(fM02EtTrIso1);
 186
 187   fM02EtTrIso2 = new TH2F("fM02EtTrIso2","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 188   fM02EtTrIso2->Sumw2();
 189   fOutputList->Add(fM02EtTrIso2);
 190
 191   fM02EtTrIso5 = new TH2F("fM02EtTrIso5","M02 vs Et for all clusters (ISO_{TRK}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 192   fM02EtTrIso5->Sumw2();
 193   fOutputList->Add(fM02EtTrIso5);
 194
 195   fM02EtAllIso1 = new TH2F("fM02EtAllIso1","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 196   fM02EtAllIso1->Sumw2();
 197   fOutputList->Add(fM02EtAllIso1);
 198
 199   fM02EtAllIso2 = new TH2F("fM02EtAllIso2","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 200   fM02EtAllIso2->Sumw2();
 201   fOutputList->Add(fM02EtAllIso2);
 202
 203   fM02EtAllIso5 = new TH2F("fM02EtAllIso5","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 204   fM02EtAllIso5->Sumw2();
 205   fOutputList->Add(fM02EtAllIso5);
 206
 207   fCeIsoVsEtPho = new TH2F("fCeIsoVsEtPho","ISO_{EMC} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC}",1000,0,100,1000,-10,190);
 208   fCeIsoVsEtPho->Sumw2();
 209   fOutputList->Add(fCeIsoVsEtPho);
 210
 211   fTrIsoVsEtPho = new TH2F("fTrIsoVsEtPho","ISO_{TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{TRK}",1000,0,100,1000,-10,190);
 212   fTrIsoVsEtPho->Sumw2();
 213   fOutputList->Add(fTrIsoVsEtPho);
 214
 215   fAllIsoVsEtPho = new TH2F("fAllIsoVsEtPho","ISO_{EMC+TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC+TRK}",1000,0,100,1000,-10,190);
 216   fAllIsoVsEtPho->Sumw2();
 217   fOutputList->Add(fAllIsoVsEtPho);
 218
 219   fM02EtCeIso1TM = new TH2F("fM02EtCeIso1TM","M02 vs Et for all clusters (ISO_{EMC}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 220   fM02EtCeIso1TM->Sumw2();
 221   fOutputList->Add(fM02EtCeIso1TM);
 222
 223   fM02EtCeIso2TM = new TH2F("fM02EtCeIso2TM","M02 vs Et for all clusters (ISO_{EMC}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 224   fM02EtCeIso2TM->Sumw2();
 225   fOutputList->Add(fM02EtCeIso2TM);
 226
 227   fM02EtCeIso5TM = new TH2F("fM02EtCeIso5TM","M02 vs Et for all clusters (ISO_{EMC}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 228   fM02EtCeIso5TM->Sumw2();
 229   fOutputList->Add(fM02EtCeIso5TM);
 230
 231   fM02EtTrIso1TM = new TH2F("fM02EtTrIso1TM","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 232   fM02EtTrIso1TM->Sumw2();
 233   fOutputList->Add(fM02EtTrIso1TM);
 234
 235   fM02EtTrIso2TM = new TH2F("fM02EtTrIso2TM","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 236   fM02EtTrIso2TM->Sumw2();
 237   fOutputList->Add(fM02EtTrIso2TM);
 238
 239   fM02EtTrIso5TM = new TH2F("fM02EtTrIso5TM","M02 vs Et for all clusters (ISO_{TRK}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 240   fM02EtTrIso5TM->Sumw2();
 241   fOutputList->Add(fM02EtTrIso5TM);
 242
 243   fM02EtAllIso1TM = new TH2F("fM02EtAllIso1TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 244   fM02EtAllIso1TM->Sumw2();
 245   fOutputList->Add(fM02EtAllIso1TM);
 246
 247   fM02EtAllIso2TM = new TH2F("fM02EtAllIso2TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 248   fM02EtAllIso2TM->Sumw2();
 249   fOutputList->Add(fM02EtAllIso2TM);
 250
 251   fM02EtAllIso5TM = new TH2F("fM02EtAllIso5TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
 252   fM02EtAllIso5TM->Sumw2();
 253   fOutputList->Add(fM02EtAllIso5TM);
 254
 255   fCeIsoVsEtPhoTM = new TH2F("fCeIsoVsEtPhoTM","ISO_{EMC} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC}",1000,0,100,1000,-10,190);
 256   fCeIsoVsEtPhoTM->Sumw2();
 257   fOutputList->Add(fCeIsoVsEtPhoTM);
 258
 259   fTrIsoVsEtPhoTM = new TH2F("fTrIsoVsEtPhoTM","ISO_{TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{TRK}",1000,0,100,1000,-10,190);
 260   fTrIsoVsEtPhoTM->Sumw2();
 261   fOutputList->Add(fTrIsoVsEtPhoTM);
 262
 263   fAllIsoVsEtPhoTM = new TH2F("fAllIsoVsEtPhoTM","ISO_{EMC+TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC+TRK}",1000,0,100,1000,-10,190);
 264   fAllIsoVsEtPhoTM->Sumw2();
 265   fOutputList->Add(fAllIsoVsEtPhoTM);
 266
 267   PostData(1, fOutputList);
 268 }
 269
 270 //________________________________________________________________________
 271 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 272 {
 273   // Main loop, called for each event.
 274
 275   // event trigger selection
 276   Bool_t isSelected = 0;
 277   if(fPeriod.Contains("11a"))
 278     isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 279   else
 280     isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 281   if(!isSelected )
 282         return;
 283
 284   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 285   if (!fESD) {
 286     printf("ERROR: fESD not available\n");
 287     return;
 288   }
 289
 290   fEvtSel->Fill(0);
 291
 292   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
 293   if(!pv)
 294     return;
 295   fPVtxZ->Fill(pv->GetZ());
 296   if(TMath::Abs(pv->GetZ())>15)
 297     return;
 298
 299   fEvtSel->Fill(1);
 300
 301   // Track loop to fill a pT spectrum
 302   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 303     AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 304     if (!track)
 305       continue;
 306     if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
 307       fSelPrimTracks->Add(track);
 308       //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
 309     }
 310   }
 311
 312   if(!fIsTrain){
 313     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 314       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 315         break;
 316       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 317     }
 318   }
 319   fESD->GetEMCALClusters(fCaloClusters);
 320   fEMCalCells = fESD->GetEMCALCells();
 321
 322
 323   FillClusHists();
 324
 325   fCaloClusters->Clear();
 326   fSelPrimTracks->Clear();
 327   PostData(1, fOutputList);
 328 }
 329
 330 //________________________________________________________________________
 331 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 332 {
 333   // Fill cluster histograms.
 334
 335   if(!fCaloClusters)
 336     return;
 337   const Int_t nclus = fCaloClusters->GetEntries();
 338   if(nclus==0)
 339     return;
 340   Double_t maxE;
 341   Int_t nthresholds = 0;
 342   for(Int_t ic=0;ic<nclus;ic++){
 343     maxE=0;
 344     AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 345     if(!c)
 346       continue;
 347     if(!c->IsEMCAL())
 348       continue;
 349     Short_t id;
 350     Double_t Emax = GetMaxCellEnergy( c, id);
 351     Double_t Ecross = GetCrossEnergy( c, id);
 352     if((1-Ecross/Emax)>fExoticCut)
 353       continue;
 354     Float_t clsPos[3] = {0,0,0};
 355     c->GetPosition(clsPos);
 356     TVector3 clsVec(clsPos);
 357     Float_t ceiso, cephiband, cecore;
 358     Float_t triso, trphiband, trcore;
 359     Float_t alliso, allphiband, allcore;
 360     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 361     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 362     GetCeIso(clsVec, ceiso, cephiband, cecore);
 363     GetTrIso(clsVec, triso, trphiband, trcore);
 364     alliso = ceiso + triso;
 365     allphiband = cephiband + trphiband;
 366     allcore = cecore + trcore;
 367     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 368     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 369     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 370     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 371     fM02Et->Fill(Et, c->GetM02());
 372     if(ceiso>1)
 373       fM02EtCeIso1->Fill(Et, c->GetM02());
 374     if(ceiso>2)
 375       fM02EtCeIso2->Fill(Et, c->GetM02());
 376     if(ceiso>5)
 377       fM02EtCeIso5->Fill(Et, c->GetM02());
 378     if(triso>1)
 379       fM02EtTrIso1->Fill(Et, c->GetM02());
 380     if(triso>2)
 381       fM02EtTrIso2->Fill(Et, c->GetM02());
 382     if(triso>5)
 383       fM02EtTrIso5->Fill(Et, c->GetM02());
 384     if(alliso>1)
 385       fM02EtAllIso1->Fill(Et, c->GetM02());
 386     if(alliso>2)
 387       fM02EtAllIso2->Fill(Et, c->GetM02());
 388     if(alliso>5)
 389       fM02EtAllIso5->Fill(Et, c->GetM02());
 390     if(c->GetM02()>0.1 && c->GetM02()<0.3){
 391       fCeIsoVsEtPho->Fill(Et, ceiso - cecore - ceisoue);
 392       fTrIsoVsEtPho->Fill(Et, triso - trcore - trisoue);
 393       fAllIsoVsEtPho->Fill(Et, alliso - allcore - allisoue);
 394       }
 395     Double_t dR = TMath::Sqrt(pow(c->GetTrackDx(),2)+pow(c->GetTrackDz(),2));
 396     if(dR<0.014){
 397       fM02EtTM->Fill(Et, c->GetM02());
 398       if(ceiso>1)
 399         fM02EtCeIso1TM->Fill(Et, c->GetM02());
 400       if(ceiso>2)
 401         fM02EtCeIso2TM->Fill(Et, c->GetM02());
 402       if(ceiso>5)
 403         fM02EtCeIso5TM->Fill(Et, c->GetM02());
 404       if(triso>1)
 405         fM02EtTrIso1TM->Fill(Et, c->GetM02());
 406       if(triso>2)
 407         fM02EtTrIso2TM->Fill(Et, c->GetM02());
 408       if(triso>5)
 409         fM02EtTrIso5TM->Fill(Et, c->GetM02());
 410       if(alliso>1)
 411         fM02EtAllIso1TM->Fill(Et, c->GetM02());
 412       if(alliso>2)
 413         fM02EtAllIso2TM->Fill(Et, c->GetM02());
 414       if(alliso>5)
 415         fM02EtAllIso5TM->Fill(Et, c->GetM02());
 416       if(c->GetM02()>0.1 && c->GetM02()<0.3){
 417         fCeIsoVsEtPhoTM->Fill(Et, ceiso);
 418         fTrIsoVsEtPhoTM->Fill(Et, triso);
 419         fAllIsoVsEtPhoTM->Fill(Et, alliso);
 420       }
 421     }
 422     if(c->E()>maxE)
 423       maxE = c->E();
 424   }
 425   fNClusHighClusE->Fill(maxE,nclus);
 426 }
 427
 428 //________________________________________________________________________
 429 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 430 {
 431   // Get cell isolation.
 432
 433   if(!fEMCalCells)
 434     return;
 435   const Int_t ncells = fEMCalCells->GetNumberOfCells();
 436   Float_t totiso=0;
 437   Float_t totband=0;
 438   Float_t totcore=0;
 439   Float_t etacl = vec.Eta();
 440   Float_t phicl = vec.Phi();
 441   Float_t thetacl = vec.Theta();
 442   if(phicl<0)
 443     phicl+=TMath::TwoPi();
 444   Int_t absid = -1;
 445   Float_t eta=-1, phi=-1;
 446   for(int icell=0;icell<ncells;icell++){
 447     absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
 448     if(!fGeom)
 449       return;
 450     fGeom->EtaPhiFromIndex(absid,eta,phi);
 451     Float_t dphi = TMath::Abs(phi-phicl);
 452     Float_t deta = TMath::Abs(eta-etacl);
 453     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 454     Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 455     if(R<fIsoConeR){
 456       totiso += etcell;
 457       if(R<0.04)
 458         totcore += etcell;
 459     }
 460     else{
 461       if(dphi>fIsoConeR)
 462         continue;
 463       if(deta<fIsoConeR)
 464         continue;
 465       totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 466     }
 467   }
 468   iso = totiso;
 469   phiband = totband;
 470   core = totcore;
 471 }
 472
 473 //________________________________________________________________________
 474 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 475 {
 476   // Get track isolation.
 477
 478   if(!fSelPrimTracks)
 479     return;
 480   const Int_t ntracks = fSelPrimTracks->GetEntries();
 481   Double_t totiso=0;
 482   Double_t totband=0;
 483   Double_t totcore=0;
 484   Double_t etacl = vec.Eta();
 485   Double_t phicl = vec.Phi();
 486   if(phicl<0)
 487     phicl+=TMath::TwoPi();
 488   for(int itrack=0;itrack<ntracks;itrack++){
 489     AliESDtrack *track = static_cast<AliESDtrack*> (fSelPrimTracks->At(itrack));
 490     if(!track)
 491       continue;
 492     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 493     Double_t deta = TMath::Abs(track->Eta()-etacl);
 494     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 495     Double_t pt = track->Pt();
 496     if(R<fIsoConeR){
 497       totiso += track->Pt();
 498       if(R<0.04)
 499         totcore += pt;
 500     }
 501     else{
 502       if(dphi>fIsoConeR)
 503         continue;
 504       if(deta<fIsoConeR)
 505         continue;
 506       totband += track->Pt();
 507     }
 508   }
 509   iso = totiso;
 510   phiband = totband;
 511   core = totcore;
 512 }
 513
 514 //________________________________________________________________________
 515 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 516 {
 517   // Calculate the energy of cross cells around the leading cell.
 518
 519   AliVCaloCells *cells = 0;
 520   cells = fESD->GetEMCALCells();
 521   if (!cells)
 522     return 0;
 523
 524   if (!fGeom)
 525     return 0;
 526
 527   Int_t iSupMod = -1;
 528   Int_t iTower  = -1;
 529   Int_t iIphi   = -1;
 530   Int_t iIeta   = -1;
 531   Int_t iphi    = -1;
 532   Int_t ieta    = -1;
 533   Int_t iphis   = -1;
 534   Int_t ietas   = -1;
 535
 536   Double_t crossEnergy = 0;
 537
 538   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 539   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 540
 541   Int_t ncells = cluster->GetNCells();
 542   for (Int_t i=0; i<ncells; i++) {
 543     Int_t cellAbsId = cluster->GetCellAbsId(i);
 544     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 545     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 546     Int_t aphidiff = TMath::Abs(iphi-iphis);
 547     if (aphidiff>1)
 548       continue;
 549     Int_t aetadiff = TMath::Abs(ieta-ietas);
 550     if (aetadiff>1)
 551       continue;
 552     if ( (aphidiff==1 && aetadiff==0) ||
 553         (aphidiff==0 && aetadiff==1) ) {
 554       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 555     }
 556   }
 557
 558   return crossEnergy;
 559 }
 560
 561 //________________________________________________________________________
 562 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 563 {
 564   // Get maximum energy of attached cell.
 565
 566   id = -1;
 567
 568   AliVCaloCells *cells = 0;
 569   cells = fESD->GetEMCALCells();
 570   if (!cells)
 571     return 0;
 572
 573   Double_t maxe = 0;
 574   Int_t ncells = cluster->GetNCells();
 575   for (Int_t i=0; i<ncells; i++) {
 576     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 577     if (e>maxe) {
 578       maxe = e;
 579       id   = cluster->GetCellAbsId(i);
 580     }
 581   }
 582   return maxe;
 583 }
 584
 585 //________________________________________________________________________
 586 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
 587 {
 588   // Called once at the end of the query.
 589 }