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