]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx
Fix by Raphaelle for JIRA ticket ALIROOT-5590
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetv2QA.cxx
CommitLineData
51a8c8fd 1// Jet v2 task using QA method, based on jet sample task (S.Aiola).
2//
3// Authors: Jason Mueller (CERN summer student 2014) & Alice Ohlson
4
5
6#include <TClonesArray.h>
7#include <TH1F.h>
8#include <TH2F.h>
9#include <TH3F.h>
10#include <TList.h>
11#include <TLorentzVector.h>
12#include <TF1.h>
13
14#include "AliVCluster.h"
15#include "AliAODCaloCluster.h"
16#include "AliESDCaloCluster.h"
17#include "AliVTrack.h"
18#include "AliEmcalJet.h"
19#include "AliRhoParameter.h"
20#include "AliLog.h"
21#include "AliJetContainer.h"
22#include "AliParticleContainer.h"
23#include "AliClusterContainer.h"
24#include "AliPicoTrack.h"
25
26#include "AliAnalysisTaskEmcalJetv2QA.h"
27
28ClassImp(AliAnalysisTaskEmcalJetv2QA)
29
30//________________________________________________________________________
31AliAnalysisTaskEmcalJetv2QA::AliAnalysisTaskEmcalJetv2QA() :
32AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetv2QA", kTRUE),
33 nCentBins(0),
34 nJetPtBins(0),
35 centBins(0x0),
36 jetPtBins(0x0),
37 fJetv2(0x0),
38 doPtWeight(kFALSE),
39 fHistTracksPt(0),
40 fHistClustersPt(0),
41 fHistLeadingJetPt(0),
42 fHistLeadingJetPtCorr(0),
43 fHistJetsPhiEta(0),
44 fHistJetsPtArea(0),
45 fHistJetsPtLeadHad(0),
46 fHistJetsCorrPtArea(0),
47 fHistPtDEtaDPhiTrackClus(0),
48 fHistPtDEtaDPhiClusTrack(0),
49 fDPhiJet(0),
50 fDPhiJetPythia(0),
51 fDPhiEP(0),
52 hGx(0),
53 hGy2(0),
54 hGxGy2(0),
55 hGy4(0),
56 hGx2(0),
57 hGx2Gy2(0),
58 hGxGy4(0),
59 hGy6(0),
60 hGx2Gy4(0),
61 hGxGy6(0),
62 hGy8(0),
63 hGy(0),
64 hN(0),
65 htv2std(0),
66 htjv2std(0),
67 htj2v2std(0),
68 hV0jv2std(0),
69 htdPsi(0),
70 htjdPsi(0),
71 htj2dPsi(0),
72 hV0jdPsi(0),
73 hAx(0),
74 hAxDijet(0),
75 hQx(0),
76 hQy(0),
77 hEventData(0),
78 hNTracks(0),
79 hNTracksCent(0),
80 hGxTracks(0),
81 hGyTracks(0),
82 hGy2Tracks(0),
83 hGxGy2Tracks(0),
84 hGy4Tracks(0),
85 htEPRes(0),
86 htjEPRes(0),
87 htj2EPRes(0),
88 fJetsCont(0),
89 fTracksCont(0),
90 fCaloClustersCont(0)
91{
92 // Default constructor.
93
94 cout << endl;
95 cout << "*******************************************************" << endl;
96 cout << "*** AliAnalysisTaskEmcalJetv2QA constructor! ***" << endl;
97 cout << "*******************************************************" << endl;
98 cout << endl;
99
100
101
102 SetMakeGeneralHistograms(kTRUE);
103}
104
105//________________________________________________________________________
106AliAnalysisTaskEmcalJetv2QA::AliAnalysisTaskEmcalJetv2QA(const char *name) :
107 AliAnalysisTaskEmcalJet(name, kTRUE),
108 nCentBins(0),
109 nJetPtBins(0),
110 centBins(0x0),
111 jetPtBins(0x0),
112 fJetv2(0x0),
113 doPtWeight(kFALSE),
114 fHistTracksPt(0),
115 fHistClustersPt(0),
116 fHistLeadingJetPt(0),
117 fHistLeadingJetPtCorr(0),
118 fHistJetsPhiEta(0),
119 fHistJetsPtArea(0),
120 fHistJetsPtLeadHad(0),
121 fHistJetsCorrPtArea(0),
122 fHistPtDEtaDPhiTrackClus(0),
123 fHistPtDEtaDPhiClusTrack(0),
124 fDPhiJet(0),
125 fDPhiJetPythia(0),
126 fDPhiEP(0),
127 hGx(0),
128 hGy2(0),
129 hGxGy2(0),
130 hGy4(0),
131 hGx2(0),
132 hGx2Gy2(0),
133 hGxGy4(0),
134 hGy6(0),
135 hGx2Gy4(0),
136 hGxGy6(0),
137 hGy8(0),
138 hGy(0),
139 hN(0),
140 htv2std(0),
141 htjv2std(0),
142 htj2v2std(0),
143 hV0jv2std(0),
144 htdPsi(0),
145 htjdPsi(0),
146 htj2dPsi(0),
147 hV0jdPsi(0),
148 hAx(0),
149 hAxDijet(0),
150 hQx(0),
151 hQy(0),
152 hEventData(0),
153 hNTracks(0),
154 hNTracksCent(0),
155 hGxTracks(0),
156 hGyTracks(0),
157 hGy2Tracks(0),
158 hGxGy2Tracks(0),
159 hGy4Tracks(0),
160 htEPRes(0),
161 htjEPRes(0),
162 htj2EPRes(0),
163 fJetsCont(0),
164 fTracksCont(0),
165 fCaloClustersCont(0)
166{
167 // Standard constructor.
168
169
170 SetMakeGeneralHistograms(kTRUE);
171}
172
173//________________________________________________________________________
174AliAnalysisTaskEmcalJetv2QA::~AliAnalysisTaskEmcalJetv2QA()
175{
176 // Destructor.
177}
178
179//________________________________________________________________________
180void AliAnalysisTaskEmcalJetv2QA::UserCreateOutputObjects()
181{
182 // Create user output.
183
184 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
185
186 fJetsCont = GetJetContainer(0);
187 if(fJetsCont) { //get particles and clusters connected to jets
188 fTracksCont = fJetsCont->GetParticleContainer();
189 fCaloClustersCont = fJetsCont->GetClusterContainer();
190 } else { //no jets, just analysis tracks and clusters
191 fTracksCont = GetParticleContainer(0);
192 fCaloClustersCont = GetClusterContainer(0);
193 }
194
195 nCentBins = 6; // set number of centrality bins
196 nJetPtBins = 4; // set number of jetPtBins
197 centBins = new Double_t[nCentBins+1];
198 jetPtBins = new Double_t[nJetPtBins+1];
199 centBins[0] = 0.; centBins[1] = 5.; centBins[2] = 10.; centBins[3] = 20.; centBins[4] = 30.; centBins[5] = 40.; centBins[6] = 50.; // set edges of bins
200 jetPtBins[0] = 50.; jetPtBins[1] = 70.; jetPtBins[2] = 90.; jetPtBins[3] = 120.; jetPtBins[4] = 200.; // set edges of bins
201
202 fTracksCont->SetClassName("AliVTrack");
203 fCaloClustersCont->SetClassName("AliVCluster");
204
205 fHistTracksPt = new TH1F("fHistTracksPt","fHistTracksPt", fNbins / 2, fMinBinPt, fMaxBinPt / 2);
206 fHistTracksPt->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
207 fHistTracksPt->GetYaxis()->SetTitle("counts");
208 fOutput->Add(fHistTracksPt);
209
210 fHistClustersPt = new TH1F("fHistClustersPt", "fHistClustersPt", fNbins / 2, fMinBinPt, fMaxBinPt / 2);
211 fHistClustersPt->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
212 fHistClustersPt->GetYaxis()->SetTitle("counts");
213 fOutput->Add(fHistClustersPt);
214
215 fHistLeadingJetPt = new TH1F("fHistLeadingJetPt", "fHistLeadingJetPt", fNbins, fMinBinPt, fMaxBinPt);
216 fHistLeadingJetPt->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
217 fHistLeadingJetPt->GetYaxis()->SetTitle("counts");
218 fOutput->Add(fHistLeadingJetPt);
219
220 fHistLeadingJetPtCorr = new TH1F("fHistLeadingJetPtCorr", "fHistLeadingJetPtCorr", fNbins, fMinBinPt, fMaxBinPt);
221 fHistLeadingJetPtCorr->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
222 fHistLeadingJetPtCorr->GetYaxis()->SetTitle("counts");
223 fOutput->Add(fHistLeadingJetPtCorr);
224
225 fHistJetsPhiEta = new TH2F("fHistJetsPhiEta", "fHistJetsPhiEta", 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
226 fHistJetsPhiEta->GetXaxis()->SetTitle("#eta");
227 fHistJetsPhiEta->GetYaxis()->SetTitle("#phi");
228 fOutput->Add(fHistJetsPhiEta);
229
230 fHistJetsPtArea = new TH2F("fHistJetsPtArea", "fHistJetsPtArea", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
231 fHistJetsPtArea->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
232 fHistJetsPtArea->GetYaxis()->SetTitle("area");
233 fOutput->Add(fHistJetsPtArea);
234
235 fHistJetsPtLeadHad = new TH2F("fHistJetsPtLeadHad", "fHistJetsPtLeadHad", fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
236 fHistJetsPtLeadHad->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
237 fHistJetsPtLeadHad->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
238 fHistJetsPtLeadHad->GetZaxis()->SetTitle("counts");
239 fOutput->Add(fHistJetsPtLeadHad);
240
241 fHistJetsCorrPtArea = new TH2F("fHistJetsCorrPtArea", "fHistJetsCorrPtArea", fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
242 fHistJetsCorrPtArea->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
243 fHistJetsCorrPtArea->GetYaxis()->SetTitle("area");
244 fOutput->Add(fHistJetsCorrPtArea);
245
246 fHistPtDEtaDPhiTrackClus = new TH3F("fHistPtDEtaDPhiTrackClus","fHistPtDEtaDPhiTrackClus;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
247 fOutput->Add(fHistPtDEtaDPhiTrackClus);
248
249 fHistPtDEtaDPhiClusTrack = new TH3F("fHistPtDEtaDPhiClusTrack","fHistPtDEtaDPhiClusTrack;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
250 fOutput->Add(fHistPtDEtaDPhiClusTrack);
251
252 fDPhiJet = new TH1F("fDPhiJet","fDPhiJet",90, -TMath::Pi()/3, 5*TMath::Pi()/3);
253 fDPhiJet->GetXaxis()->SetTitle("#Delta#varphi");
254 fDPhiJet->GetYaxis()->SetTitle("counts");
255 fOutput->Add(fDPhiJet);
256
257 fDPhiJetPythia = new TH1F("fDPhiJetPythia","fDPhiJetPythia",90, -TMath::Pi()/3, 5*TMath::Pi()/3);
258 fDPhiJetPythia->GetXaxis()->SetTitle("#Delta#varphi");
259 fDPhiJetPythia->GetYaxis()->SetTitle("counts");
260 fOutput->Add(fDPhiJetPythia);
261
262 fDPhiEP = new TH1F("fDPhiEP","fDPhiEP",90, 0, 2*TMath::Pi());
263 fDPhiEP->GetXaxis()->SetTitle("#Delta#varphi");
264 fDPhiEP->GetYaxis()->SetTitle("counts");
265 fOutput->Add(fDPhiEP);
266
267 hGx = new TH2D("hGx", "Gx v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
268 hGx->GetXaxis()->SetTitle("Centrality (%)");
269 hGx->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
270 fOutput->Add(hGx);
271
272 hGy2 = new TH2D("hGy2", "Gy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
273 hGy2->GetXaxis()->SetTitle("Centrality (%)");
274 hGy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
275 fOutput->Add(hGy2);
276
277 hGxGy2 = new TH2D("hGxGy2", "GxGy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
278 hGxGy2->GetXaxis()->SetTitle("Centrality (%)");
279 hGxGy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
280 fOutput->Add(hGxGy2);
281
282 hGy4 = new TH2D("hGy4", "Gy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
283 hGy4->GetXaxis()->SetTitle("Centrality (%)");
284 hGy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
285 fOutput->Add(hGy4);
286
287 hGx2 = new TH2D("hGx2", "Gx2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
288 hGx2->GetXaxis()->SetTitle("Centrality (%)");
289 hGx2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
290 fOutput->Add(hGx2);
291
292 hGx2Gy2 = new TH2D("hGx2Gy2", "Gx2Gy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
293 hGx2Gy2->GetXaxis()->SetTitle("Centrality (%)");
294 hGx2Gy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
295 fOutput->Add(hGx2Gy2);
296
297 hGxGy4 = new TH2D("hGxGy4", "GxGy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
298 hGxGy4->GetXaxis()->SetTitle("Centrality (%)");
299 hGxGy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
300 fOutput->Add(hGxGy4);
301
302 hGy6 = new TH2D("hGy6", "Gy6 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
303 hGy6->GetXaxis()->SetTitle("Centrality (%)");
304 hGy6->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
305 fOutput->Add(hGy6);
306
307 hGx2Gy4 = new TH2D("hGx2Gy4", "Gx2Gy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
308 hGx2Gy4->GetXaxis()->SetTitle("Centrality (%)");
309 hGx2Gy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
310 fOutput->Add(hGx2Gy4);
311
312 hGxGy6 = new TH2D("hGxGy6", "GxGy6 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
313 hGxGy6->GetXaxis()->SetTitle("Centrality (%)");
314 hGxGy6->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
315 fOutput->Add(hGxGy6);
316
317 hGy8 = new TH2D("hGy8", "Gy8 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
318 hGy8->GetXaxis()->SetTitle("Centrality (%)");
319 hGy8->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
320 fOutput->Add(hGy8);
321
322 hGy = new TH2D("hGy", "Gy v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
323 hGy->GetXaxis()->SetTitle("Centrality (%)");
324 hGy->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
325 fOutput->Add(hGy);
326
327 hN = new TH2D("hN", "N v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
328 hN->GetXaxis()->SetTitle("Centrality (%)");
329 hN->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
330 fOutput->Add(hN);
331
332 htv2std = new TH2D("htv2std", "v2std v Centrality v JetPt w/o jets", nCentBins, centBins, nJetPtBins, jetPtBins);
333 htv2std->GetXaxis()->SetTitle("Centrality (%)");
334 htv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
335 fOutput->Add(htv2std);
336
337 htjv2std = new TH2D("htjv2std", "v2std v Centrality v JetPt w/ jets", nCentBins, centBins, nJetPtBins, jetPtBins);
338 htjv2std->GetXaxis()->SetTitle("Centrality (%)");
339 htjv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
340 fOutput->Add(htjv2std);
341
342 htj2v2std = new TH2D("htj2v2std", "v2std v Centrality v JetPt w/ trackPt < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins);
343 htj2v2std->GetXaxis()->SetTitle("Centrality (%)");
344 htj2v2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
345 fOutput->Add(htj2v2std);
346
347 hV0jv2std = new TH2D("hV0jv2std", "v2std v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
348 hV0jv2std->GetXaxis()->SetTitle("Centrality (%)");
349 hV0jv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
350 fOutput->Add(hV0jv2std);
351
352 Int_t ndpsibins = 100;
353 Double_t dpsibins[101];
354 for(Int_t t = 0; t < 101; t++) dpsibins[t] = TMath::Pi()*t/50.;
355
356 htdPsi = new TH3D("htdPsi", "JetAxis - EventPlane w/o jets", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
357 htdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
358 fOutput->Add(htdPsi);
359
360 htjdPsi = new TH3D("htjdPsi", "JetAxis - EventPlane w/ jets", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
361 htjdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
362 fOutput->Add(htjdPsi);
363
364 htj2dPsi = new TH3D("htj2dPsi", "JetAxis - EventPlane w/ trackPt < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
365 htj2dPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
366 fOutput->Add(htj2dPsi);
367
368 hV0jdPsi = new TH3D("hV0jdPsi", "JetAxis - EventPlane", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
369 hV0jdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
370 fOutput->Add(hV0jdPsi);
371
372 hQx = new TH2D("hQx", "Qx Distribution in EP frame", 100, -0.3, 0.3, nCentBins, centBins);
373 hQx->GetXaxis()->SetTitle("Qx");
374 hQx->GetYaxis()->SetTitle("Centrality (%)");
375 fOutput->Add(hQx);
376
377 hQy = new TH2D("hQy", "Qy Distribution in EP frame", 100, -0.3, 0.3, nCentBins, centBins);
378 hQy->GetXaxis()->SetTitle("Qy");
379 hQy->GetYaxis()->SetTitle("Centrality (%)");
380 fOutput->Add(hQy);
381
382 hAx = new TH2D("hAx", "Ax Distribution in EP frame w/o Dijets", 100, -35, 70, nJetPtBins, jetPtBins);
383 hAx->GetXaxis()->SetTitle("Ax");
384 hAx->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
385 fOutput->Add(hAx);
386
387 hAxDijet = new TH2D("hAxDijet", "Ax Distribution in EP frame w/ Dijets", 100, -35, 70, nJetPtBins, jetPtBins);
388 hAxDijet->GetXaxis()->SetTitle("Ax");
389 hAxDijet->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
390 fOutput->Add(hAxDijet);
391
392 hEventData = new TH1F("hEventData","Events Kept and Discarded", 9, 0, 9);
393 hEventData->GetYaxis()->SetTitle("counts");
394 fOutput->Add(hEventData);
395
396 hNTracks = new TH2F("hNTracks","Number of Tracks Per Event", 100, 0, 3000, 3, 0, 3);
397 hNTracks->GetXaxis()->SetTitle("# tracks");
398 fOutput->Add(hNTracks);
399
400 hNTracksCent = new TH2D("hNTracksCent", "NTracks by centrality", 100, 0, 3000, (Int_t)centBins[nCentBins], centBins[0], centBins[nCentBins]);
401 hNTracksCent->GetXaxis()->SetTitle("# tracks");
402 hNTracksCent->GetYaxis()->SetTitle("Centrality (%)");
403 fOutput->Add(hNTracksCent);
404
405 hGxTracks = new TH2D("hGxTracks", "Gx by NTracks", 200, -200, 200, 100, 0, 3000);
406 hGxTracks->GetXaxis()->SetTitle("Gx");
407 hGxTracks->GetYaxis()->SetTitle("# tracks");
408 fOutput->Add(hGxTracks);
409
410 hGyTracks = new TH2D("hGyTracks", "Gy by NTracks", 200, -200, 200, 100, 0, 3000);
411 hGyTracks->GetXaxis()->SetTitle("Gy");
412 hGyTracks->GetYaxis()->SetTitle("# tracks");
413 fOutput->Add(hGyTracks);
414
415 hGy2Tracks = new TH2D("hGy2Tracks", "Gy2 by NTracks", 100, 0, 20000, 100, 0, 3000);
416 hGy2Tracks->GetXaxis()->SetTitle("Gy2");
417 hGy2Tracks->GetYaxis()->SetTitle("# tracks");
418 fOutput->Add(hGy2Tracks);
419
420 hGxGy2Tracks = new TH2D("hGxGy2Tracks", "GxGy2 by NTracks", 100, -100000, 100000, 100, 0, 3000);
421 hGxGy2Tracks->GetXaxis()->SetTitle("GxGy2");
422 hGxGy2Tracks->GetYaxis()->SetTitle("# tracks");
423 fOutput->Add(hGxGy2Tracks);
424
425 hGy4Tracks = new TH2D("hGy4Tracks", "Gy4 by NTracks", 100, 0, 100000000, 100, 0, 3000);
426 hGy4Tracks->GetXaxis()->SetTitle("Gy4");
427 hGy4Tracks->GetYaxis()->SetTitle("# tracks");
428 fOutput->Add(hGy4Tracks);
429
430 htEPRes = new TH2D("htEPRes", "EP Resolution w/o Jets", nCentBins, centBins, nJetPtBins, jetPtBins);
431 htEPRes->GetXaxis()->SetTitle("Centrality (%)");
432 htEPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
433 fOutput->Add(htEPRes);
434
435 htjEPRes = new TH2D("htjEPRes", "EP Resolution w/ Jets", nCentBins, centBins, nJetPtBins, jetPtBins);
436 htjEPRes->GetXaxis()->SetTitle("Centrality (%)");
437 htjEPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
438 fOutput->Add(htjEPRes);
439
440 htj2EPRes = new TH2D("htj2EPRes", "EP Resolution w/ trackPT < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins);
441 htj2EPRes->GetXaxis()->SetTitle("Centrality (%)");
442 htj2EPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
443 fOutput->Add(htj2EPRes);
444
445 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
446}
447
448//________________________________________________________________________
449Bool_t AliAnalysisTaskEmcalJetv2QA::FillHistograms()
450{
451 // Fill histograms.
452
453 return kTRUE;
454}
455
456
457//________________________________________________________________________
458void AliAnalysisTaskEmcalJetv2QA::ExecOnce() {
459
460 AliAnalysisTaskEmcalJet::ExecOnce();
461
462 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
463 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
464 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
465
466}
467
468//________________________________________________________________________
469Bool_t AliAnalysisTaskEmcalJetv2QA::Run() // this part loops over each event
470{
471 // Run analysis code here, if needed. It will be executed before FillHistograms().
472
473 Double_t jetPhi = -999; // angle of leading jet
474 Double_t jetPt = -999; // pt of leading jet
475 Double_t jetArea = -999;
476 Double_t trackPt = -999;
477 Double_t phi = -999; // track phi
478 Double_t dPhi = -999; // track phi - jet axis
479 Double_t dPhiQA = -999; // track phi - EP
480 Double_t tSin = 0; // used for std EP calc
481 Double_t tCos = 0; // used for std EP calc
482 Double_t jSin = 0; // used for std EP calc
483 Double_t jCos = 0; // used for std EP calc
484 Double_t tSin2 = 0; // used for std EP calc with trackPt < 2 GeV
485 Double_t tCos2 = 0; // used for std EP calc with trackPt < 2 GeV
486 Double_t qx = 0; // used for Qx distribution
487 Double_t qy = 0; // used for Qy distribution
488 Double_t ax = 0; // used for Ax distribution
489 Double_t tEP = -999; // EP w/o jets
490 Double_t tjEP = -999; // EP w/ jets
491 Double_t tjEP2 = -999; // EP w/ jets
492 Double_t dPsi = -999; // jet axis - EP
493 Double_t gx = 0; // used for G moment calc
494 Double_t gy = 0; // used for G moment calc
495 Int_t isDijet = 0; // if 0, no dijet. if 1, dijet.
496 Int_t nTracksBkgnd = 0; // used to keep track of number of tracks in background
497 Int_t nTracksJet =0; // used to keep track of number of tracks in jets
498
499 if(fJetsCont)
500 {
501
502 AliEmcalJet *jettest = fJetsCont->GetNextAcceptJet(0);
503 while(jettest) {
504
505 fHistJetsPtArea->Fill(jettest->Pt(), jettest->Area());
506 fHistJetsPhiEta->Fill(jettest->Eta(), jettest->Phi());
507
508 Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jettest);
509 fHistJetsPtLeadHad->Fill(jettest->Pt(), ptLeading);
510
511 if (fHistJetsCorrPtArea) {
512 Float_t corrPt = jettest->Pt() - fJetsCont->GetRhoVal() * jettest->Area();
513 fHistJetsCorrPtArea->Fill(corrPt, jettest->Area());
514 }
515 jettest = fJetsCont->GetNextAcceptJet();
516 }
517
518
519 AliEmcalJet *jet = fJetsCont->GetLeadingJet();
520 if(jet)
521 {
522 jetPhi = jet->Phi(); // get leading jet phi (jet axis)
523 jetPt = jet->Pt(); // get leading jet pT to filter out low pT jet events
524 jetArea = jet->Area();
525 }
526 }
527
528 // event cuts
529 if(!fTracksCont)
530 {
531 hEventData->Fill("!fTracksCont",1);
532 return kFALSE;
533 }
534 if(!fJetsCont)
535 {
536 hEventData->Fill("!fJetsCont",1);
537 return kFALSE;
538 }
539 if(jetPt == -999)
540 {
541 hEventData->Fill("jetPt=-999",1);
542 return kFALSE;
543 }
544 if(jetPt < jetPtBins[0])
545 {
546 hEventData->Fill("leadingJetPt<jetPtMin",1);
547 return kFALSE;
548 }
549 if(jetPt > jetPtBins[nJetPtBins])
550 {
551 hEventData->Fill("leadingJetPt>jetPtMax",1);
552 return kFALSE;
553 }
554 if(fCent < centBins[0])
555 {
556 hEventData->Fill("cent<centMin",1);
557 return kFALSE;
558 }
559 if(fCent > centBins[nCentBins])
560 {
561 hEventData->Fill("cent>centMax",1);
562 return kFALSE;
563 }
564 hEventData->Fill("good event",1);
565
566 AliEmcalJet *dijet = fJetsCont->GetNextAcceptJet(0); // check for dijet events
567 while(dijet)
568 {
569 if(dijet->Pt() > 50 && fabs(jetPhi-dijet->Phi()-TMath::Pi()) < 0.4) // loop over jets with pT>50 and exclude leading jet and check that angular separation is < 0.4
570 isDijet = 1;
571 dijet = fJetsCont->GetNextAcceptJet();
572 }
573
574 if (fCaloClustersCont)
575 {
576 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
577 while(cluster) {
578 TLorentzVector nPart;
579 cluster->GetMomentum(nPart, fVertex);
580 fHistClustersPt->Fill(nPart.Pt());
581 cluster = fCaloClustersCont->GetNextAcceptCluster();
582 }
583 }
584
585 fHistLeadingJetPt->Fill(jetPt);
586 fHistLeadingJetPtCorr->Fill(jetPt-fJetsCont->GetRhoVal()*jetArea);
587
588 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
589 while(track)
590 { // loop over all particles (including jet tracks)
591 trackPt = track->Pt();
592 fHistTracksPt->Fill(trackPt);
593 phi = track->Phi(); // get track phi
594
595 dPhi = phi-jetPhi; // get track phi - jet axis
596 if(dPhi < 0) dPhi += TMath::TwoPi();
597 if(dPhi > TMath::TwoPi()) dPhi -= TMath::TwoPi();
598
599 dPhiQA = phi-fEPV0; // get track phi - EP
600 if(dPhiQA < 0) dPhiQA += TMath::TwoPi();
601 if(dPhiQA > TMath::TwoPi()) dPhiQA -= TMath::TwoPi();
602 fDPhiEP->Fill(dPhiQA);
603
604 // fill jet-hadron correlation just to check if track labels make sense...
605 if(track->Pt()>1.)
606 {
607 Double_t dphiJet = dPhi;
608 if(dphiJet > 5*TMath::Pi()/3) dphiJet -= 2*TMath::Pi();
609 if(dphiJet < -TMath::Pi()/3) dphiJet += 2*TMath::Pi();
610 if(track->GetLabel() == 0)
611 fDPhiJet->Fill(dphiJet);
612 else
613 fDPhiJetPythia->Fill(dphiJet);
614 }
615
616 Double_t weight = 1.;
617 if(doPtWeight) weight = trackPt;
618
619 gx += weight*cos(2*dPhi);
620 gy += weight*sin(2*dPhi);
621
622 if(track->GetLabel() == 0)
623 { // sum for std EP method
624 tSin += weight*sin(2*phi); // bkgnd has label = 0
625 tCos += weight*cos(2*phi);
626 qx += weight*cos(2*dPhiQA);
627 qy += weight*sin(2*dPhiQA);
628 nTracksBkgnd += 1;
629 }
630 else
631 {
632 jSin += weight*sin(2*phi); // jets have label =/= 0
633 jCos += weight*cos(2*phi);
634 ax += weight*cos(2*dPhi);
635 nTracksJet += 1;
636 }
637
638 if(track->Pt() < 2.)
639 { // sum for std EP method w/ trackPt < 2 GeV
640 tSin2 += weight*sin(2*phi);
641 tCos2 += weight*cos(2*phi);
642 }
643
644
645 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()); // increment to next track
646 } // close loop over particles
647
648 hNTracks->Fill(nTracksBkgnd,"Bkgnd Tracks",1);
649 hNTracks->Fill(nTracksJet,"Jet Tracks",1);
650 hNTracks->Fill(nTracksBkgnd+nTracksJet,"Total Tracks",1);
651 hNTracksCent->Fill(nTracksBkgnd+nTracksJet,fCent,1);
652
653 if(nTracksBkgnd == 0)
654 {
655 hEventData->Fill("no tracks",1);
656 hEventData->Fill("good event",-1);
657 return kFALSE;
658 }
659
660 Double_t v2weight = 1+2*fJetv2*cos(2*(jetPhi-fEPV0)); // set v2 weight for event
661
662 hGx->Fill(fCent, jetPt, v2weight*gx); // fill histograms for QA method
663 hGy2->Fill(fCent, jetPt, v2weight*gy*gy);
664 hGxGy2->Fill(fCent, jetPt, v2weight*gx*gy*gy);
665 hGy4->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy);
666 hGx2->Fill(fCent, jetPt, v2weight*gx*gx);
667 hGx2Gy2->Fill(fCent, jetPt, v2weight*gx*gx*gy*gy);
668 hGxGy4->Fill(fCent, jetPt, v2weight*gx*gy*gy*gy*gy);
669 hGy6->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy*gy*gy);
670 hGx2Gy4->Fill(fCent, jetPt, v2weight*gx*gx*gy*gy*gy*gy);
671 hGxGy6->Fill(fCent, jetPt, v2weight*gx*gy*gy*gy*gy*gy*gy);
672 hGy8->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy*gy*gy*gy*gy);
673 hGy->Fill(fCent, jetPt, v2weight*gy);
674 hN->Fill(fCent, jetPt, v2weight);
675
676 hGxTracks->Fill(gx,nTracksBkgnd+nTracksJet);
677 hGyTracks->Fill(gy,nTracksBkgnd+nTracksJet);
678 hGy2Tracks->Fill(gy*gy,nTracksBkgnd+nTracksJet);
679 hGxGy2Tracks->Fill(gx*gy*gy,nTracksBkgnd+nTracksJet);
680 hGy4Tracks->Fill(gy*gy*gy*gy,nTracksBkgnd+nTracksJet);
681 hQx->Fill(qx/nTracksBkgnd,fCent);
682 hQy->Fill(qy/nTracksBkgnd,fCent);
683
684 if(isDijet == 0)
685 hAx->Fill(ax,jetPt);
686 if(isDijet == 1)
687 hAxDijet->Fill(ax,jetPt);
688
689 tEP = 0.5*atan2(tSin,tCos); // calculate EP w/o jets
690 tjEP = 0.5*atan2((tSin+jSin),(tCos+jCos)); // calculate EP w/ jets
691 tjEP2 = 0.5*atan2(tSin2,tCos2); // calculate EP w/ trackPt < 2 GeV
692
693 htEPRes->Fill(fCent, jetPt, v2weight*cos(2*(tEP-fEPV0)));
694 htjEPRes->Fill(fCent, jetPt, v2weight*cos(2*(tjEP-fEPV0)));
695 htj2EPRes->Fill(fCent, jetPt, v2weight*cos(2*(tjEP2-fEPV0)));
696
697
698 dPsi = jetPhi-tEP;
699 if(dPsi < 0) dPsi += TMath::TwoPi();
700 if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
701
702 htv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/o jets
703 htdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/o jets
704
705 dPsi = jetPhi-tjEP;
706 if(dPsi < 0) dPsi += TMath::TwoPi();
707 if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
708
709 htjv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/ jets
710 htjdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/ jets
711
712 dPsi = jetPhi-tjEP2;
713 if(dPsi < 0) dPsi += TMath::TwoPi();
714 if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
715
716 htj2v2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/ trackPt < 2 GeV
717 htj2dPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/ trackPt < 2 GeV
718
719 dPsi = jetPhi-fEPV0;
720 if(dPsi < 0) dPsi += TMath::TwoPi();
721 if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
722
723 hV0jv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data
724 hV0jdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EPV0
725
726
727 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
728}
729
730//________________________________________________________________________
731void AliAnalysisTaskEmcalJetv2QA::Terminate(Option_t *)
732{
733 // Called once at the end of the analysis.
734 if(centBins) delete [] centBins;
735 if(jetPtBins) delete [] jetPtBins;
736}