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