c90a58d110d36945b897103e24f43e2a230f56ca
[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
27 #include "AliAnalysisTaskSAQA.h"
28
29 ClassImp(AliAnalysisTaskSAQA)
30
31 //________________________________________________________________________
32 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() : 
33   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
34   fCellEnergyCut(0.1),
35   fDoTrigger(kFALSE),
36   fRepropagateTracks(kFALSE),
37   fTrgClusName("ClustersL1GAMMAFEE"),
38   fTrgClusters(0),
39   fHistCentrality(0),
40   fHistZVertex(0),
41   fHistTracksCent(0),
42   fHistClusCent(0),
43   fHistClusTracks(0),
44   fHistCellsCent(0),
45   fHistCellsTracks(0),
46   fHistMaxL1FastORCent(0),
47   fHistMaxL1ClusCent(0),
48   fHistMaxL1ThrCent(0),
49   fHistTracksPt(0),
50   fHistTrPhiEta(0),
51   fHistTrEmcPhiEta(0),
52   fHistTrPhiEtaNonProp(0),
53   fHistDeltaEtaPt(0),
54   fHistDeltaPhiPt(0),
55   fHistDeltaEtaNewProp(0),
56   fHistDeltaPhiNewProp(0),
57   fHistClusPhiEtaEnergy(0),
58   fHistNCellsEnergy(0),
59   fHistClusTimeEnergy(0),
60   fHistCellsEnergy(0),
61   fHistChVSneCells(0),
62   fHistChVSneClus(0),
63   fHistChVSneCorrCells(0)
64 {
65   // Default constructor.
66
67   for (Int_t i = 0; i < 5; i++) {
68     fHistTrackPhi[i] = 0;
69     fHistTrackEta[i] = 0;
70   }
71
72   for (Int_t i = 0; i < 6; i++) {
73     fHistTrackPhiPt[i] = 0;
74     fHistTrackEtaPt[i] = 0;
75   }
76
77   for (Int_t i = 0; i < 4; i++) {
78     fHistJetsPhiEta[i] = 0;
79     fHistJetsPtNonBias[i] = 0;
80     fHistJetsPtTrack[i] = 0;
81     fHistJetsPtClus[i] = 0;
82     fHistJetsPt[i] = 0;
83     fHistJetsPtArea[i] = 0;
84     fHistJetsPtAreaNonBias[i] = 0;
85   }
86 }
87
88 //________________________________________________________________________
89 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) : 
90   AliAnalysisTaskEmcalJet(name, kTRUE),
91   fCellEnergyCut(0.1),
92   fDoTrigger(kFALSE),
93   fRepropagateTracks(kFALSE),
94   fTrgClusName("ClustersL1GAMMAFEE"),
95   fTrgClusters(0),
96   fHistCentrality(0),
97   fHistZVertex(0),
98   fHistTracksCent(0),
99   fHistClusCent(0),
100   fHistClusTracks(0),
101   fHistCellsCent(0),
102   fHistCellsTracks(0),
103   fHistMaxL1FastORCent(0),
104   fHistMaxL1ClusCent(0),
105   fHistMaxL1ThrCent(0),
106   fHistTracksPt(0),
107   fHistTrPhiEta(0),
108   fHistTrEmcPhiEta(0),
109   fHistTrPhiEtaNonProp(0),
110   fHistDeltaEtaPt(0),
111   fHistDeltaPhiPt(0),
112   fHistDeltaEtaNewProp(0),
113   fHistDeltaPhiNewProp(0),
114   fHistClusPhiEtaEnergy(0),
115   fHistNCellsEnergy(0),
116   fHistClusTimeEnergy(0),
117   fHistCellsEnergy(0),
118   fHistChVSneCells(0),
119   fHistChVSneClus(0),
120   fHistChVSneCorrCells(0)
121 {
122   // Standard constructor.
123
124   for (Int_t i = 0; i < 5; i++) {
125     fHistTrackPhi[i] = 0;
126     fHistTrackEta[i] = 0;
127   }
128
129   for (Int_t i = 0; i < 6; i++) {
130     fHistTrackPhiPt[i] = 0;
131     fHistTrackEtaPt[i] = 0;
132   }
133
134   for (Int_t i = 0; i < 4; i++) {
135     fHistJetsPhiEta[i] = 0;
136     fHistJetsPtNonBias[i] = 0;
137     fHistJetsPtTrack[i] = 0;
138     fHistJetsPtClus[i] = 0;
139     fHistJetsPt[i] = 0;
140     fHistJetsPtArea[i] = 0;
141     fHistJetsPtAreaNonBias[i] = 0;
142   }
143 }
144
145 //________________________________________________________________________
146 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
147 {
148   // Destructor
149 }
150
151 //________________________________________________________________________
152 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
153 {
154   // Create histograms
155
156   OpenFile(1);
157   fOutput = new TList();
158   fOutput->SetOwner();
159
160   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
161   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
162   fHistCentrality->GetYaxis()->SetTitle("counts");
163   fOutput->Add(fHistCentrality);
164
165   fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
166   fHistZVertex->GetXaxis()->SetTitle("z");
167   fHistZVertex->GetYaxis()->SetTitle("counts");
168   fOutput->Add(fHistZVertex);
169
170   fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
171   fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
172   fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
173   fOutput->Add(fHistTracksCent);
174
175   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
176     fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
177     fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
178     fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
179     fOutput->Add(fHistClusCent);
180
181     fHistClusTracks = new TH2F("fHistClusTracks","Clusters vs. tracks", fNbins, 0, 4000, fNbins, 0, 2000);
182     fHistClusTracks->GetXaxis()->SetTitle("No. of tracks");
183     fHistClusTracks->GetYaxis()->SetTitle("No. of clusters");
184     fOutput->Add(fHistClusTracks);
185
186     fHistCellsCent = new TH2F("fHistCellsCent","Cells vs. centrality", 100, 0, 100, fNbins, 0, 6000);
187     fHistCellsCent->GetXaxis()->SetTitle("Centrality (%)");
188     fHistCellsCent->GetYaxis()->SetTitle("No. of EMCal cells");
189     fOutput->Add(fHistCellsCent);
190
191     fHistCellsTracks = new TH2F("fHistCellsTracks","Cells vs. tracks", fNbins, 0, 4000, fNbins, 0, 6000);
192     fHistCellsTracks->GetXaxis()->SetTitle("No. of tracks");
193     fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
194     fOutput->Add(fHistCellsTracks);
195
196     fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  -1e-6, 1e-6);
197     fHistClusTimeEnergy->GetXaxis()->SetTitle("Energy (GeV)");
198     fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
199     fOutput->Add(fHistClusTimeEnergy);
200
201     if (fDoTrigger) {
202       fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
203       fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
204       fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
205       fOutput->Add(fHistMaxL1FastORCent);
206       
207       fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
208       fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
209       fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
210       fOutput->Add(fHistMaxL1ClusCent);
211       
212       fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
213       fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
214       fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
215       fOutput->Add(fHistMaxL1ThrCent);
216     }
217   }
218
219   fHistTracksPt = new TH1F("fHistTracksPt","p_{T} spectrum of reconstructed tracks", fNbins, fMinBinPt, fMaxBinPt);
220   fHistTracksPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
221   fHistTracksPt->GetYaxis()->SetTitle("counts");
222   fOutput->Add(fHistTracksPt);
223        
224   fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 80, -2, 2, 128, 0, 6.4);
225   fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
226   fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
227   fOutput->Add(fHistTrPhiEta);
228
229   fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 80, -2, 2, 128, 0, 6.4);
230   fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
231   fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
232   fOutput->Add(fHistTrEmcPhiEta);
233
234   fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 80, -2, 2, 128, 0, 6.4);
235   fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
236   fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
237   fOutput->Add(fHistTrPhiEtaNonProp);
238
239   fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
240   fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
241   fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
242   fOutput->Add(fHistDeltaEtaPt);
243
244   fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
245   fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
246   fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
247   fOutput->Add(fHistDeltaPhiPt);
248
249   if (fRepropagateTracks) {
250     fHistDeltaEtaNewProp = new TH1F("fHistDeltaEtaNewProp","fHistDeltaEtaNewProp", 800, -0.4, 0.4);
251     fHistDeltaEtaNewProp->GetXaxis()->SetTitle("#delta#eta");
252     fHistDeltaEtaNewProp->GetYaxis()->SetTitle("counts");
253     fOutput->Add(fHistDeltaEtaNewProp);
254     
255     fHistDeltaPhiNewProp = new TH1F("fHistDeltaPhiNewProp","fHistDeltaPhiNewProp", 2560, -1.6, 3.2);
256     fHistDeltaPhiNewProp->GetXaxis()->SetTitle("#delta#phi");
257     fHistDeltaPhiNewProp->GetYaxis()->SetTitle("counts");
258     fOutput->Add(fHistDeltaPhiNewProp);
259   }
260
261   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
262     fHistClusPhiEtaEnergy = new TH3F("fHistClusPhiEtaEnergy","Phi-Eta-Energy distribution of clusters", fNbins, fMinBinPt, fMaxBinPt, 80, -2, 2, 128, 0, 6.4);
263     fHistClusPhiEtaEnergy->GetXaxis()->SetTitle("E [GeV]");
264     fHistClusPhiEtaEnergy->GetYaxis()->SetTitle("#eta");
265     fHistClusPhiEtaEnergy->GetZaxis()->SetTitle("#phi");
266     fOutput->Add(fHistClusPhiEtaEnergy);
267
268     fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
269     fHistNCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
270     fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
271     fOutput->Add(fHistNCellsEnergy);
272   }
273
274   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
275    
276     fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
277     fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
278     fHistCellsEnergy->GetYaxis()->SetTitle("counts");
279     fOutput->Add(fHistCellsEnergy);
280     
281     fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy", 
282                                 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
283     fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
284     fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
285     fOutput->Add(fHistChVSneCells);
286     
287     fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy", 
288                                (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
289     fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
290     fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
291     fOutput->Add(fHistChVSneClus);
292     
293     fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy", 
294                                     (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
295     fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
296     fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
297     fOutput->Add(fHistChVSneCorrCells);
298   }
299  
300   for (Int_t i = 0; i < 5; i++) {
301     TString histnamephi("fHistTrackPhi_");
302     histnamephi += i;
303     fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 128, 0, 6.4);
304     fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
305     fOutput->Add(fHistTrackPhi[i]);
306
307     TString histnameeta("fHistTrackEta_");
308     histnameeta += i;
309     fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
310     fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
311     fOutput->Add(fHistTrackEta[i]);
312   }
313
314   for (Int_t i = 0; i < 6; i++) {
315     TString histnamephipt("fHistTrackPhiPt_");
316     histnamephipt += i;
317     fHistTrackPhiPt[i] = new TH1F(histnamephipt.Data(),histnamephipt.Data(), 128, 0, 6.4);
318     fHistTrackPhiPt[i]->GetXaxis()->SetTitle("Phi");
319     fOutput->Add(fHistTrackPhiPt[i]);
320
321     TString histnameetapt("fHistTrackEtaPt_");
322     histnameetapt += i;
323     fHistTrackEtaPt[i] = new TH1F(histnameetapt.Data(),histnameetapt.Data(), 100, -2, 2);
324     fHistTrackEtaPt[i]->GetXaxis()->SetTitle("Eta");
325     fOutput->Add(fHistTrackEtaPt[i]);
326   }
327   
328   fHistTrackPhi[0]->SetLineColor(kRed);
329   fHistTrackEta[0]->SetLineColor(kRed);
330   fHistTrackPhi[1]->SetLineColor(kBlue);
331   fHistTrackEta[1]->SetLineColor(kBlue);
332   fHistTrackPhi[2]->SetLineColor(kGreen);
333   fHistTrackEta[2]->SetLineColor(kGreen);
334   fHistTrackPhi[3]->SetLineColor(kOrange);
335   fHistTrackEta[3]->SetLineColor(kOrange);
336   fHistTrackPhi[4]->SetLineColor(kBlack);
337   fHistTrackEta[4]->SetLineColor(kBlack);
338
339   if (!fJetsName.IsNull()) {
340
341     TString histname;
342
343     for (Int_t i = 0; i < 4; i++) {
344       histname = "fHistJetsPhiEta_";
345       histname += i;
346       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 80, -2, 2, 128, 0, 6.4);
347       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
348       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
349       fHistJetsPhiEta[i]->GetZaxis()->SetTitle("p_{T} [GeV/c]");
350       fOutput->Add(fHistJetsPhiEta[i]);
351
352       histname = "fHistJetsPtNonBias_";
353       histname += i;
354       fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
355       fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
356       fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
357       fOutput->Add(fHistJetsPtNonBias[i]);
358
359       histname = "fHistJetsPtTrack_";
360       histname += i;
361       fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
362       fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
363       fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
364       fOutput->Add(fHistJetsPtTrack[i]);
365
366       if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
367         histname = "fHistJetsPtClus_";
368         histname += i;
369         fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
370         fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
371         fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
372         fOutput->Add(fHistJetsPtClus[i]);
373       }
374
375       histname = "fHistJetsPt_";
376       histname += i;
377       fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
378       fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
379       fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
380       fOutput->Add(fHistJetsPt[i]);
381
382       histname = "fHistJetsPtAreaNonBias_";
383       histname += i;
384       fHistJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
385       fHistJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
386       fHistJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
387       fOutput->Add(fHistJetsPtAreaNonBias[i]);
388
389       histname = "fHistJetsPtArea_";
390       histname += i;
391       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
392       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
393       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
394       fOutput->Add(fHistJetsPtArea[i]);
395     }
396   }
397
398   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
399 }
400
401 //________________________________________________________________________
402 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
403 {
404   // Retrieve event objects.
405
406   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
407     return kFALSE;
408
409   if (!fTrgClusName.IsNull() && fDoTrigger && !fTrgClusters) {
410     fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
411     if (!fTrgClusters) {
412       AliError(Form("%s: Could not retrieve trigger clusters %s!", GetName(), fTrgClusName.Data())); 
413       return kFALSE;
414     }
415     else {
416       TClass *cl = fTrgClusters->GetClass();
417       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
418         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fTrgClusName.Data())); 
419         fTrgClusters = 0;
420         return kFALSE;
421       }
422     }
423   }
424
425   return kTRUE;
426 }
427
428 //________________________________________________________________________
429 Bool_t AliAnalysisTaskSAQA::FillHistograms()
430 {
431   // Fill histograms.
432
433   fHistCentrality->Fill(fCent);
434   if (fTracks)
435     fHistTracksCent->Fill(fCent, fTracks->GetEntriesFast());
436   if (fCaloClusters)
437     fHistClusCent->Fill(fCent, fCaloClusters->GetEntriesFast());
438
439   fHistZVertex->Fill(fVertex[2]);
440
441   Float_t trackSum = DoTrackLoop();
442
443   DoJetLoop();
444
445   if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
446
447     if (fTracks && fCaloClusters)
448       fHistClusTracks->Fill(fTracks->GetEntriesFast(), fCaloClusters->GetEntriesFast());
449
450     Float_t clusSum = DoClusterLoop();
451
452     Float_t cellSum = 0, cellCutSum = 0;
453     
454     Int_t ncells = DoCellLoop(cellSum, cellCutSum);
455
456     if (fTracks)
457       fHistCellsTracks->Fill(fTracks->GetEntriesFast(), ncells);
458
459     fHistCellsCent->Fill(fCent, ncells);
460     
461     fHistChVSneCells->Fill(cellSum, trackSum);
462     fHistChVSneClus->Fill(clusSum, trackSum);
463     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
464
465     if (fDoTrigger) {
466       Float_t maxTrgClus = DoTriggerClusLoop();
467       fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
468     
469       Int_t maxL1amp = -1;
470       Int_t maxL1thr = -1;
471     
472       DoTriggerPrimitives(maxL1amp, maxL1thr);
473       
474       if (maxL1amp > -1) 
475         fHistMaxL1FastORCent->Fill(fCent, maxL1amp);
476       
477       if (maxL1thr > -1) 
478         fHistMaxL1ThrCent->Fill(fCent, maxL1thr);
479     }
480   }
481
482   return kTRUE;
483 }
484
485 //________________________________________________________________________
486 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
487 {
488   // Do cell loop.
489
490   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
491
492   if (!cells)
493     return 0;
494
495   const Int_t ncells = cells->GetNumberOfCells();
496
497   for (Int_t pos = 0; pos < ncells; pos++) {
498     Float_t amp = cells->GetAmplitude(pos);
499     fHistCellsEnergy->Fill(amp);
500     sum += amp;
501     if (amp < fCellEnergyCut)
502       continue;
503     sum_cut += amp;
504   } 
505
506   return ncells;
507 }
508
509 //________________________________________________________________________
510 Float_t AliAnalysisTaskSAQA::DoClusterLoop()
511 {
512   // Do cluster loop.
513
514   if (!fCaloClusters)
515     return 0;
516
517   Float_t sum = 0;
518
519   // Cluster loop
520   Int_t nclusters =  fCaloClusters->GetEntriesFast();
521
522   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
523     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
524     if (!cluster) {
525       AliError(Form("Could not receive cluster %d", iClusters));
526       continue;
527     }  
528
529     if (!AcceptCluster(cluster, kTRUE))
530       continue;
531
532     sum += cluster->E();
533
534     TLorentzVector nPart;
535     cluster->GetMomentum(nPart, fVertex);
536
537     fHistClusPhiEtaEnergy->Fill(cluster->E(), nPart.Eta(), nPart.Phi());
538     fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
539
540     fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
541   }
542
543   return sum;
544 }
545
546 //________________________________________________________________________
547 Float_t AliAnalysisTaskSAQA::DoTrackLoop()
548 {
549   // Do track loop.
550
551   if (!fTracks)
552     return 0;
553
554   Float_t sum = 0;
555
556   Int_t ntracks = fTracks->GetEntriesFast();
557   Int_t nclusters = 0;
558   if (fCaloClusters)
559     nclusters = fCaloClusters->GetEntriesFast();
560
561   for (Int_t i = 0; i < ntracks; i++) {
562
563     AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track  
564
565     if (!track) {
566       AliError(Form("Could not retrieve track %d",i)); 
567       continue; 
568     }
569
570     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
571     
572     if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
573       continue;
574     
575     fHistTracksPt->Fill(track->Pt());
576
577     sum += track->P();
578
579     Int_t label = track->GetLabel();
580       
581     fHistTrPhiEta->Fill(track->Eta(), track->Phi());
582     
583     fHistTrackEta[4]->Fill(track->Eta());
584     fHistTrackPhi[4]->Fill(track->Phi());
585
586     if (label >= 0 && label < 4) {
587       fHistTrackEta[label]->Fill(track->Eta());
588       fHistTrackPhi[label]->Fill(track->Phi());
589     }
590
591     if (track->Pt() < 0.5) {
592       fHistTrackPhiPt[0]->Fill(track->Phi());
593       fHistTrackEtaPt[0]->Fill(track->Eta());
594     }
595     else if (track->Pt() < 1) {
596       fHistTrackPhiPt[1]->Fill(track->Phi());
597       fHistTrackEtaPt[1]->Fill(track->Eta());
598     }
599     else if (track->Pt() < 2) {
600       fHistTrackPhiPt[2]->Fill(track->Phi());
601       fHistTrackEtaPt[2]->Fill(track->Eta());
602     }
603     else if (track->Pt() < 3) {
604       fHistTrackPhiPt[3]->Fill(track->Phi());
605       fHistTrackEtaPt[3]->Fill(track->Eta());
606     }
607     else if (track->Pt() < 5) {
608       fHistTrackPhiPt[4]->Fill(track->Phi());
609       fHistTrackEtaPt[4]->Fill(track->Eta());
610     }
611     else {
612       fHistTrackPhiPt[5]->Fill(track->Phi());
613       fHistTrackEtaPt[5]->Fill(track->Eta());
614     }
615
616     if (!vtrack)
617       continue;
618
619     if (vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999)
620       fHistTrPhiEtaNonProp->Fill(vtrack->Eta(), vtrack->Phi());
621
622     fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
623     fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
624     fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
625
626     if (fRepropagateTracks && vtrack->GetTrackEtaOnEMCal() > -2) {    
627       Float_t propeta = -999, propphi = -999;
628       PropagateTrack(vtrack, propeta, propphi);
629       fHistDeltaEtaNewProp->Fill(propeta - vtrack->GetTrackEtaOnEMCal());
630       fHistDeltaPhiNewProp->Fill(propphi - vtrack->GetTrackPhiOnEMCal());
631     }
632   }
633   
634   return sum;
635 }
636
637 //____________________________________________________________________________
638 void AliAnalysisTaskSAQA::PropagateTrack(AliVTrack *track, Float_t &eta, Float_t &phi)
639 {
640   eta = -999;
641   phi = -999;
642
643   if (!track)
644     return;
645
646   if (track->Pt() == 0) 
647     return;
648
649   // init the magnetic field if not already on
650   if(!TGeoGlobalMagField::Instance()->GetField()) {
651     AliInfo("Init the magnetic field\n");
652     AliAODEvent* aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
653     if (aodevent) {
654       Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
655       Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
656       AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
657       TGeoGlobalMagField::Instance()->SetField(field);
658     }
659   }
660     
661   Double_t cv[21];
662   for (Int_t i = 0; i < 21; i++) cv[i] = 0;
663
664   Double_t pos[3], mom[3];
665   track->GetXYZ(pos);
666   track->GetPxPyPz(mom);
667   AliExternalTrackParam *trackParam = new AliExternalTrackParam(pos, mom, cv, track->Charge());
668  
669   if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, 430., 0.139, 20, kTRUE, 0.8, -1)) return;
670   Double_t trkPos[3] = {0., 0., 0.};
671   if(!trackParam->GetXYZ(trkPos)) return;
672   TVector3 trkPosVec(trkPos[0], trkPos[1], trkPos[2]);
673   eta = trkPosVec.Eta();
674   phi = trkPosVec.Phi();
675   if(phi < 0)
676     phi += 2 * TMath::Pi();
677 }
678
679 //________________________________________________________________________
680 void AliAnalysisTaskSAQA::DoJetLoop()
681 {
682   // Do jet loop.
683
684   if (!fJets)
685     return;
686
687   Int_t njets = fJets->GetEntriesFast();
688
689   for (Int_t ij = 0; ij < njets; ij++) {
690
691     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
692
693     if (!jet) {
694       AliError(Form("Could not receive jet %d", ij));
695       continue;
696     }  
697
698     if (!AcceptJet(jet, kFALSE))
699       continue;
700
701     fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
702     fHistJetsPtAreaNonBias[fCentBin]->Fill(jet->Pt(), jet->Area());
703
704     if (jet->MaxTrackPt() > fPtBiasJetTrack)
705       fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
706     
707     if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
708       fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
709     
710     if (!AcceptBiasJet(jet))
711       continue;
712
713     fHistJetsPt[fCentBin]->Fill(jet->Pt());
714     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
715
716     fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
717   }
718 }
719
720 //________________________________________________________________________
721 Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
722 {
723   // Do trigger cluster loop.
724
725   if (!fTrgClusters)
726     return 0;
727
728   Int_t ntrgclusters = fTrgClusters->GetEntriesFast();
729   Float_t maxTrgClus = 0;
730
731   for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
732     AliVCluster* cluster = static_cast<AliVCluster*>(fTrgClusters->At(iClusters));
733     if (!cluster) {
734       AliError(Form("Could not receive cluster %d", iClusters));
735       continue;
736     }  
737     
738     if (!(cluster->IsEMCAL())) continue;
739
740     if (cluster->E() > maxTrgClus)
741       maxTrgClus = cluster->E();
742
743   }
744   return maxTrgClus;
745 }
746
747 //________________________________________________________________________
748 void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
749 {
750   // Do trigger primitives loop.
751
752   AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
753
754   if (!triggers || triggers->GetEntries() == 0)
755     return;
756     
757   triggers->Reset();
758   Int_t L1amp = 0;
759   Int_t L1thr = 0;
760   maxL1amp = -1;
761   maxL1thr = -1;
762
763   while (triggers->Next()) {
764     triggers->GetL1TimeSum(L1amp);
765     if (maxL1amp < L1amp) 
766       maxL1amp = L1amp;
767
768     triggers->GetL1Threshold(L1thr);
769     if (maxL1thr < L1thr) 
770       maxL1thr = L1thr;
771   }
772 }
773
774 //________________________________________________________________________
775 void AliAnalysisTaskSAQA::Terminate(Option_t *) 
776 {
777   // Called once at the end of the analysis.
778 }