PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx

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