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 = 0;
 341   for(Int_t ic=0;ic<nclus;ic++){
 342     maxE=0;
 343     AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 344     if(!c)
 345       continue;
 346     if(!c->IsEMCAL())
 347       continue;
 348     Short_t id;
 349     Double_t Emax = GetMaxCellEnergy( c, id);
 350     Double_t Ecross = GetCrossEnergy( c, id);
 351     if((1-Ecross/Emax)>fExoticCut)
 352       continue;
 353     Float_t clsPos[3] = {0,0,0};
 354     c->GetPosition(clsPos);
 355     TVector3 clsVec(clsPos);
 356     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 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     fM02Et->Fill(Et, c->GetM02());
 371     if(ceiso-cecore<1)
 372       fM02EtCeIso1->Fill(Et, c->GetM02());
 373     if(ceiso-cecore<2)
 374       fM02EtCeIso2->Fill(Et, c->GetM02());
 375     if(ceiso-cecore<5)
 376       fM02EtCeIso5->Fill(Et, c->GetM02());
 377     if(triso-trcore<1)
 378       fM02EtTrIso1->Fill(Et, c->GetM02());
 379     if(triso-trcore<2)
 380       fM02EtTrIso2->Fill(Et, c->GetM02());
 381     if(triso-trcore<5)
 382       fM02EtTrIso5->Fill(Et, c->GetM02());
 383     if(alliso-allcore<1)
 384       fM02EtAllIso1->Fill(Et, c->GetM02());
 385     if(alliso-allcore<2)
 386       fM02EtAllIso2->Fill(Et, c->GetM02());
 387     if(alliso-allcore<5)
 388       fM02EtAllIso5->Fill(Et, c->GetM02());
 389     if(c->GetM02()>0.1 && c->GetM02()<0.3){
 390       fCeIsoVsEtPho->Fill(Et, ceiso - cecore - ceisoue);
 391       fTrIsoVsEtPho->Fill(Et, triso - trcore - trisoue);
 392       fAllIsoVsEtPho->Fill(Et, alliso - allcore - allisoue);
 393       }
 394     Double_t dR = TMath::Sqrt(pow(c->GetTrackDx(),2)+pow(c->GetTrackDz(),2));
 395     if(dR<0.014){
 396       fM02EtTM->Fill(Et, c->GetM02());
 397       if(ceiso-cecore<1)
 398         fM02EtCeIso1TM->Fill(Et, c->GetM02());
 399       if(ceiso-cecore<2)
 400         fM02EtCeIso2TM->Fill(Et, c->GetM02());
 401       if(ceiso-cecore<5)
 402         fM02EtCeIso5TM->Fill(Et, c->GetM02());
 403       if(triso-trcore<1)
 404         fM02EtTrIso1TM->Fill(Et, c->GetM02());
 405       if(triso-trcore<2)
 406         fM02EtTrIso2TM->Fill(Et, c->GetM02());
 407       if(triso-trcore<5)
 408         fM02EtTrIso5TM->Fill(Et, c->GetM02());
 409       if(alliso-allcore<1)
 410         fM02EtAllIso1TM->Fill(Et, c->GetM02());
 411       if(alliso-allcore<2)
 412         fM02EtAllIso2TM->Fill(Et, c->GetM02());
 413       if(alliso-allcore<5)
 414         fM02EtAllIso5TM->Fill(Et, c->GetM02());
 415       if(c->GetM02()>0.1 && c->GetM02()<0.3){
 416         fCeIsoVsEtPhoTM->Fill(Et, ceiso);
 417         fTrIsoVsEtPhoTM->Fill(Et, triso);
 418         fAllIsoVsEtPhoTM->Fill(Et, alliso);
 419       }
 420     }
 421     if(c->E()>maxE)
 422       maxE = c->E();
 423   }
 424   fNClusHighClusE->Fill(maxE,nclus);
 425 }
 426
 427 //________________________________________________________________________
 428 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 429 {
 430   // Get cell isolation.
 431
 432   if(!fEMCalCells)
 433     return;
 434   const Int_t ncells = fEMCalCells->GetNumberOfCells();
 435   Float_t totiso=0;
 436   Float_t totband=0;
 437   Float_t totcore=0;
 438   Float_t etacl = vec.Eta();
 439   Float_t phicl = vec.Phi();
 440   Float_t thetacl = vec.Theta();
 441   if(phicl<0)
 442     phicl+=TMath::TwoPi();
 443   Int_t absid = -1;
 444   Float_t eta=-1, phi=-1;
 445   for(int icell=0;icell<ncells;icell++){
 446     absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
 447     if(!fGeom)
 448       return;
 449     fGeom->EtaPhiFromIndex(absid,eta,phi);
 450     Float_t dphi = TMath::Abs(phi-phicl);
 451     Float_t deta = TMath::Abs(eta-etacl);
 452     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 453     Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 454     if(R<fIsoConeR){
 455       totiso += etcell;
 456       if(R<0.04)
 457         totcore += etcell;
 458     }
 459     else{
 460       if(dphi>fIsoConeR)
 461         continue;
 462       if(deta<fIsoConeR)
 463         continue;
 464       totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 465     }
 466   }
 467   iso = totiso;
 468   phiband = totband;
 469   core = totcore;
 470 }
 471
 472 //________________________________________________________________________
 473 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 474 {
 475   // Get track isolation.
 476
 477   if(!fSelPrimTracks)
 478     return;
 479   const Int_t ntracks = fSelPrimTracks->GetEntries();
 480   Double_t totiso=0;
 481   Double_t totband=0;
 482   Double_t totcore=0;
 483   Double_t etacl = vec.Eta();
 484   Double_t phicl = vec.Phi();
 485   if(phicl<0)
 486     phicl+=TMath::TwoPi();
 487   for(int itrack=0;itrack<ntracks;itrack++){
 488     AliESDtrack *track = static_cast<AliESDtrack*> (fSelPrimTracks->At(itrack));
 489     if(!track)
 490       continue;
 491     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 492     Double_t deta = TMath::Abs(track->Eta()-etacl);
 493     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 494     Double_t pt = track->Pt();
 495     if(R<fIsoConeR){
 496       totiso += track->Pt();
 497       if(R<0.04)
 498         totcore += pt;
 499     }
 500     else{
 501       if(dphi>fIsoConeR)
 502         continue;
 503       if(deta<fIsoConeR)
 504         continue;
 505       totband += track->Pt();
 506     }
 507   }
 508   iso = totiso;
 509   phiband = totband;
 510   core = totcore;
 511 }
 512
 513 //________________________________________________________________________
 514 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 515 {
 516   // Calculate the energy of cross cells around the leading cell.
 517
 518   AliVCaloCells *cells = 0;
 519   cells = fESD->GetEMCALCells();
 520   if (!cells)
 521     return 0;
 522
 523   if (!fGeom)
 524     return 0;
 525
 526   Int_t iSupMod = -1;
 527   Int_t iTower  = -1;
 528   Int_t iIphi   = -1;
 529   Int_t iIeta   = -1;
 530   Int_t iphi    = -1;
 531   Int_t ieta    = -1;
 532   Int_t iphis   = -1;
 533   Int_t ietas   = -1;
 534
 535   Double_t crossEnergy = 0;
 536
 537   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 538   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 539
 540   Int_t ncells = cluster->GetNCells();
 541   for (Int_t i=0; i<ncells; i++) {
 542     Int_t cellAbsId = cluster->GetCellAbsId(i);
 543     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 544     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 545     Int_t aphidiff = TMath::Abs(iphi-iphis);
 546     if (aphidiff>1)
 547       continue;
 548     Int_t aetadiff = TMath::Abs(ieta-ietas);
 549     if (aetadiff>1)
 550       continue;
 551     if ( (aphidiff==1 && aetadiff==0) ||
 552         (aphidiff==0 && aetadiff==1) ) {
 553       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 554     }
 555   }
 556
 557   return crossEnergy;
 558 }
 559
 560 //________________________________________________________________________
 561 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 562 {
 563   // Get maximum energy of attached cell.
 564
 565   id = -1;
 566
 567   AliVCaloCells *cells = 0;
 568   cells = fESD->GetEMCALCells();
 569   if (!cells)
 570     return 0;
 571
 572   Double_t maxe = 0;
 573   Int_t ncells = cluster->GetNCells();
 574   for (Int_t i=0; i<ncells; i++) {
 575     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 576     if (e>maxe) {
 577       maxe = e;
 578       id   = cluster->GetCellAbsId(i);
 579     }
 580   }
 581   return maxe;
 582 }
 583
 584 //________________________________________________________________________
 585 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
 586 {
 587   // Called once at the end of the query.
 588 }