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