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