8dc65b5479bf93bf84738f788e59f328ce73224a
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskSAQA.cxx
1 // $Id$
2 //
3 // General QA task.
4 //
5 // Author: S.Aiola
6
7 #include <TChain.h>
8 #include <TClonesArray.h>
9 #include <TH1F.h>
10 #include <TH2F.h>
11 #include <TH3F.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14
15 #include "AliAnalysisManager.h"
16 #include "AliCentrality.h"
17 #include "AliVCluster.h"
18 #include "AliVParticle.h"
19 #include "AliVTrack.h"
20 #include "AliEmcalJet.h"
21 #include "AliVEventHandler.h"
22 #include "AliAODEvent.h"
23 #include "AliExternalTrackParam.h"
24 #include "AliTrackerBase.h"
25 #include "AliLog.h"
26 #include "AliEMCALGeometry.h"
27 #include "AliEMCALGeoParams.h"
28
29 #include "AliAnalysisTaskSAQA.h"
30
31 ClassImp(AliAnalysisTaskSAQA)
32
33 //________________________________________________________________________
34 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() : 
35   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
36   fCellEnergyCut(0.1),
37   fNclusters(0),
38   fNtracks(0),
39   fNjets(0),
40   fHistTracksCent(0),
41   fHistClusCent(0),
42   fHistJetsCent(0),
43   fHistClusTracks(0),
44   fHistJetsParts(0),
45   fHistCellsCent(0),
46   fHistCellsTracks(0),
47   fHistTrEmcPhiEta(0),
48   fHistTrPhiEtaNonProp(0),
49   fHistDeltaEtaPt(0),
50   fHistDeltaPhiPt(0),
51   fHistNCellsEnergy(0),
52   fHistFcrossEnergy(0),
53   fHistClusTimeEnergy(0),
54   fHistCellsAbsIdEnergy(0),
55   fHistChVSneCells(0),
56   fHistChVSneClus(0),
57   fHistChVSneCorrCells(0)
58 {
59   // Default constructor.
60
61   for (Int_t i = 0; i < 4; i++) {
62     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
63     fHistClusPhiEtaEnergy[i] = 0;
64     fHistJetsPhiEtaPt[i] = 0;
65     fHistJetsPtArea[i] = 0;
66   }
67
68   SetMakeGeneralHistograms(kTRUE);
69 }
70
71 //________________________________________________________________________
72 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) : 
73   AliAnalysisTaskEmcalJet(name, kTRUE),
74   fCellEnergyCut(0.1),
75   fNclusters(0),
76   fNtracks(0),
77   fNjets(0),
78   fHistTracksCent(0),
79   fHistClusCent(0),
80   fHistJetsCent(0),
81   fHistClusTracks(0),
82   fHistJetsParts(0),
83   fHistCellsCent(0),
84   fHistCellsTracks(0),
85   fHistTrEmcPhiEta(0),
86   fHistTrPhiEtaNonProp(0),
87   fHistDeltaEtaPt(0),
88   fHistDeltaPhiPt(0),
89   fHistNCellsEnergy(0),
90   fHistFcrossEnergy(0),
91   fHistClusTimeEnergy(0),
92   fHistCellsAbsIdEnergy(0),
93   fHistChVSneCells(0),
94   fHistChVSneClus(0),
95   fHistChVSneCorrCells(0)
96 {
97   // Standard constructor.
98
99   for (Int_t i = 0; i < 4; i++) {
100     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
101     fHistClusPhiEtaEnergy[i] = 0;
102     fHistJetsPhiEtaPt[i] = 0;
103     fHistJetsPtArea[i] = 0;
104   }
105
106   SetMakeGeneralHistograms(kTRUE);
107 }
108
109 //________________________________________________________________________
110 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
111 {
112   // Destructor
113 }
114
115 //________________________________________________________________________
116 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
117 {
118   // Create histograms
119
120   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
121
122   if (!fTracksName.IsNull()) {
123     fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
124     fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
125     fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
126     fOutput->Add(fHistTracksCent);
127
128     TString histname;
129
130     Int_t nlabels = 4;
131     if (fParticleLevel)
132       nlabels = 1;
133
134     for (Int_t i = 0; i < fNcentBins; i++) {
135       for (Int_t j = 0; j < nlabels; j++) {
136         histname = Form("fHistTrPhiEtaPt_%d_%d",i,j);
137         fHistTrPhiEtaPt[i][j] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
138         fHistTrPhiEtaPt[i][j]->GetXaxis()->SetTitle("#eta");
139         fHistTrPhiEtaPt[i][j]->GetYaxis()->SetTitle("#phi");
140         fHistTrPhiEtaPt[i][j]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
141         fOutput->Add(fHistTrPhiEtaPt[i][j]);
142       }
143     }
144
145     fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
146     fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
147     fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
148     fOutput->Add(fHistTrEmcPhiEta);
149
150     fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
151     fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
152     fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
153     fOutput->Add(fHistTrPhiEtaNonProp);
154     
155     fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
156     fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
157     fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
158     fOutput->Add(fHistDeltaEtaPt);
159     
160     fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
161     fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
162     fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
163     fOutput->Add(fHistDeltaPhiPt);
164   }
165
166   if (!fCaloName.IsNull()) {
167     fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
168     fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
169     fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
170     fOutput->Add(fHistClusCent);
171
172     fHistClusTracks = new TH2F("fHistClusTracks","Clusters vs. tracks", fNbins, 0, 4000, fNbins, 0, 2000);
173     fHistClusTracks->GetXaxis()->SetTitle("No. of tracks");
174     fHistClusTracks->GetYaxis()->SetTitle("No. of clusters");
175     fOutput->Add(fHistClusTracks);
176
177     fHistCellsCent = new TH2F("fHistCellsCent","Cells vs. centrality", 100, 0, 100, fNbins, 0, 6000);
178     fHistCellsCent->GetXaxis()->SetTitle("Centrality (%)");
179     fHistCellsCent->GetYaxis()->SetTitle("No. of EMCal cells");
180     fOutput->Add(fHistCellsCent);
181
182     fHistCellsTracks = new TH2F("fHistCellsTracks","Cells vs. tracks", fNbins, 0, 4000, fNbins, 0, 6000);
183     fHistCellsTracks->GetXaxis()->SetTitle("No. of tracks");
184     fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
185     fOutput->Add(fHistCellsTracks);
186
187     TString histname;
188
189     for (Int_t i = 0; i < fNcentBins; i++) {
190       histname = "fHistClusPhiEtaEnergy_";
191       histname += i;
192       fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
193       fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
194       fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
195       fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("E_{cluster} (GeV)");
196       fOutput->Add(fHistClusPhiEtaEnergy[i]);
197     }
198
199     fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  -1e-6, 1e-6);
200     fHistClusTimeEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
201     fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
202     fOutput->Add(fHistClusTimeEnergy);
203
204     fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
205     fHistNCellsEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
206     fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
207     fOutput->Add(fHistNCellsEnergy); 
208
209     fHistFcrossEnergy = new TH2F("fHistFcrossEnergy","fHistFcrossEnergy", fNbins, fMinBinPt, fMaxBinPt, 200, -3.5, 1.5);
210     fHistFcrossEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
211     fHistFcrossEnergy->GetYaxis()->SetTitle("F_{cross}");
212     fOutput->Add(fHistFcrossEnergy); 
213      
214     fHistCellsAbsIdEnergy = new TH2F("fHistCellsAbsIdEnergy","fHistCellsAbsIdEnergy", 11600,0,11599,(Int_t)(fNbins / 2), fMinBinPt, fMaxBinPt / 2);
215     fHistCellsAbsIdEnergy->GetXaxis()->SetTitle("cell abs. Id");
216     fHistCellsAbsIdEnergy->GetYaxis()->SetTitle("E_{cluster} (GeV)");
217     fHistCellsAbsIdEnergy->GetZaxis()->SetTitle("counts");    
218     fOutput->Add(fHistCellsAbsIdEnergy);
219     
220     fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy", 
221                                 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
222     fHistChVSneCells->GetXaxis()->SetTitle("Energy (GeV)");
223     fHistChVSneCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
224     fOutput->Add(fHistChVSneCells);
225     
226     fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy", 
227                                (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
228     fHistChVSneClus->GetXaxis()->SetTitle("Energy (GeV)");
229     fHistChVSneClus->GetYaxis()->SetTitle("Momentum (GeV/c)");
230     fOutput->Add(fHistChVSneClus);
231     
232     fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy", 
233                                     (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
234     fHistChVSneCorrCells->GetXaxis()->SetTitle("Energy (GeV)");
235     fHistChVSneCorrCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
236     fOutput->Add(fHistChVSneCorrCells);
237   }
238        
239   if (!fJetsName.IsNull()) {
240     fHistJetsCent = new TH2F("fHistJetsCent","Jets vs. centrality", 100, 0, 100, 150, -0.5, 149.5);
241     fHistJetsCent->GetXaxis()->SetTitle("Centrality (%)");
242     fHistJetsCent->GetYaxis()->SetTitle("No. of jets");
243     fOutput->Add(fHistJetsCent);
244     
245     fHistJetsParts = new TH2F("fHistJetsParts","Jets vs. centrality", fNbins, 0, 6000, 150, -0.5, 149.5);
246     fHistJetsParts->GetXaxis()->SetTitle("No. of particles");
247     fHistJetsParts->GetYaxis()->SetTitle("No. of jets");
248     fOutput->Add(fHistJetsParts);
249
250     TString histname;
251
252     for (Int_t i = 0; i < fNcentBins; i++) {
253       histname = "fHistJetsPhiEtaPt_";
254       histname += i;
255       fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01,  (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
256       fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
257       fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
258       fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
259       fOutput->Add(fHistJetsPhiEtaPt[i]);
260
261       histname = "fHistJetsPtArea_";
262       histname += i;
263       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
264       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
265       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
266       fOutput->Add(fHistJetsPtArea[i]);
267     }
268   }
269
270   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
271 }
272
273 //________________________________________________________________________
274 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
275 {
276   // Retrieve event objects.
277
278   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
279     return kFALSE;
280
281   fNclusters = 0;
282   fNtracks = 0;
283   fNjets = 0;
284
285   return kTRUE;
286 }
287
288
289 //________________________________________________________________________
290 Bool_t AliAnalysisTaskSAQA::FillHistograms()
291 {
292   // Fill histograms.
293
294   Float_t clusSum = 0;
295   Float_t trackSum = 0;
296
297   if (fTracks) {
298     trackSum = DoTrackLoop();
299
300     fHistTracksCent->Fill(fCent, fNtracks);
301   } 
302
303   if (fCaloClusters) {
304     clusSum = DoClusterLoop();
305
306     fHistClusCent->Fill(fCent, fNclusters);
307     fHistClusTracks->Fill(fNtracks, fNclusters);
308
309     Float_t cellSum = 0, cellCutSum = 0;
310     
311     Int_t ncells = DoCellLoop(cellSum, cellCutSum);
312
313     if (fTracks)
314       fHistCellsTracks->Fill(fNtracks, ncells);
315
316     fHistCellsCent->Fill(fCent, ncells);
317     
318     fHistChVSneCells->Fill(cellSum, trackSum);
319     fHistChVSneClus->Fill(clusSum, trackSum);
320     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
321   }
322
323   if (fJets) {
324     DoJetLoop();
325     
326     fHistJetsCent->Fill(fCent, fNjets);
327     fHistJetsParts->Fill(fNtracks + fNclusters, fNjets);
328   }
329
330   return kTRUE;
331 }
332
333 //________________________________________________________________________
334 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
335 {
336   // Do cell loop.
337
338   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
339
340   if (!cells)
341     return 0;
342
343   const Int_t ncells = cells->GetNumberOfCells();
344
345   for (Int_t pos = 0; pos < ncells; pos++) {
346     Float_t amp   = cells->GetAmplitude(pos);
347     Int_t   absId = cells->GetCellNumber(pos);
348     fHistCellsAbsIdEnergy->Fill(absId,amp);
349     sum += amp;
350     if (amp < fCellEnergyCut)
351       continue;
352     sum_cut += amp;
353   } 
354
355   return ncells;
356 }
357
358 //________________________________________________________________________
359 Double_t AliAnalysisTaskSAQA::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
360 {
361   Int_t    AbsIdseed  = -1;
362   Double_t Eseed      = 0;
363   for (Int_t i = 0; i < cluster->GetNCells(); i++) {
364     if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > AbsIdseed) {
365       Eseed     = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
366       AbsIdseed = cluster->GetCellAbsId(i);
367     }
368   }
369
370   if (Eseed < 1e-9)
371     return 100;
372
373   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
374   fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta); 
375   fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);  
376   
377   //Get close cells index and energy, not in corners
378   
379   Int_t absID1 = -1;
380   Int_t absID2 = -1;
381   
382   if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
383   if (iphi > 0)                                 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
384   
385   // In case of cell in eta = 0 border, depending on SM shift the cross cell index
386   
387   Int_t absID3 = -1;
388   Int_t absID4 = -1;
389   
390   if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
391     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
392     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1); 
393   }
394   else if (ieta == 0 && imod%2) {
395     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
396     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1); 
397   }
398   else {
399     if (ieta < AliEMCALGeoParams::fgkEMCALCols-1) 
400       absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
401     if (ieta > 0)                                 
402       absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1); 
403   }
404   
405   Double_t  ecell1 = cells->GetCellAmplitude(absID1);
406   Double_t  ecell2 = cells->GetCellAmplitude(absID2);
407   Double_t  ecell3 = cells->GetCellAmplitude(absID3);
408   Double_t  ecell4 = cells->GetCellAmplitude(absID4);
409
410   Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
411   
412   Double_t Fcross = 1 - Ecross/Eseed;
413
414   return Fcross;
415 }
416
417 //________________________________________________________________________
418 Float_t AliAnalysisTaskSAQA::DoClusterLoop()
419 {
420   // Do cluster loop.
421
422   if (!fCaloClusters)
423     return 0;
424
425   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
426
427   Float_t sum = 0;
428
429   // Cluster loop
430   Int_t nclusters =  fCaloClusters->GetEntriesFast();
431
432   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
433     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
434     if (!cluster) {
435       AliError(Form("Could not receive cluster %d", iClusters));
436       continue;
437     }  
438
439     if (!AcceptCluster(cluster, kTRUE))
440       continue;
441
442     sum += cluster->E();
443
444     TLorentzVector nPart;
445     cluster->GetMomentum(nPart, fVertex);
446
447     fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
448     fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
449
450     fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
451
452     if (cells)
453       fHistFcrossEnergy->Fill(cluster->E(), GetFcross(cluster, cells));
454
455     fNclusters++;
456   }
457
458   return sum;
459 }
460
461 //________________________________________________________________________
462 Float_t AliAnalysisTaskSAQA::DoTrackLoop()
463 {
464   // Do track loop.
465
466   if (!fTracks)
467     return 0;
468
469   Float_t sum = 0;
470
471   Int_t ntracks = fTracks->GetEntriesFast();
472   Int_t nclusters = 0;
473   if (fCaloClusters)
474     nclusters = fCaloClusters->GetEntriesFast();
475
476   for (Int_t i = 0; i < ntracks; i++) {
477
478     AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track  
479
480     if (!track) {
481       AliError(Form("Could not retrieve track %d",i)); 
482       continue; 
483     }
484
485     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
486     
487     if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
488       continue;
489
490     fNtracks++;
491
492     sum += track->P();
493
494     if (fParticleLevel) {
495       fHistTrPhiEtaPt[fCentBin][0]->Fill(track->Eta(), track->Phi(), track->Pt());
496     }
497     else {
498       fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
499
500       Int_t label = track->GetLabel();
501       if (label >= 0 && label < 3)
502         fHistTrPhiEtaPt[fCentBin][label]->Fill(track->Eta(), track->Phi(), track->Pt());
503       else
504         AliWarning(Form("%s: track label %d not recognized!", GetName(), label));
505     }
506
507     if (!vtrack)
508       continue;
509
510     if (vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999)
511       fHistTrPhiEtaNonProp->Fill(vtrack->Eta(), vtrack->Phi());
512
513     fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
514     fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
515     fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
516   }
517   
518   return sum;
519 }
520
521 //________________________________________________________________________
522 void AliAnalysisTaskSAQA::DoJetLoop()
523 {
524   // Do jet loop.
525
526   if (!fJets)
527     return;
528
529   Int_t njets = fJets->GetEntriesFast();
530
531   for (Int_t ij = 0; ij < njets; ij++) {
532
533     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
534
535     if (!jet) {
536       AliError(Form("Could not receive jet %d", ij));
537       continue;
538     }  
539
540     if (!AcceptJet(jet))
541       continue;
542
543     fNjets++;
544
545     fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
546     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
547   }
548 }
549
550
551 //________________________________________________________________________
552 void AliAnalysisTaskSAQA::Terminate(Option_t *) 
553 {
554   // Called once at the end of the analysis.
555 }