cd95650eabd9d946526235f405763ec9e52682e3
[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   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 }