]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskDeltaPt.cxx
... / ...
CommitLineData
1// $Id$
2//
3// Jet deltaPt task.
4//
5// Author: S.Aiola
6
7#include <TClonesArray.h>
8#include <TH1F.h>
9#include <TH2F.h>
10#include <TH3F.h>
11#include <TList.h>
12#include <TLorentzVector.h>
13#include <TRandom3.h>
14
15#include "AliVCluster.h"
16#include "AliVParticle.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
25#include "AliAnalysisTaskDeltaPt.h"
26
27ClassImp(AliAnalysisTaskDeltaPt)
28
29//________________________________________________________________________
30AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() :
31 AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
32 fMCJetPtThreshold(1),
33 fMinRC2LJ(-1),
34 fRCperEvent(-1),
35 fConeRadius(0.2),
36 fConeMinEta(-0.9),
37 fConeMaxEta(0.9),
38 fConeMinPhi(0),
39 fConeMaxPhi(TMath::Pi()*2),
40 fJetsCont(0),
41 fTracksCont(0),
42 fCaloClustersCont(0),
43 fEmbJetsCont(0),
44 fEmbTracksCont(0),
45 fEmbCaloClustersCont(0),
46 fRandTracksCont(0),
47 fRandCaloClustersCont(0),
48 fHistRCPhiEta(0),
49 fHistRCPt(0),
50 fHistRCPtExLJ(0),
51 fHistRCPtExPartialLJ(0),
52 fHistRCPtRand(0),
53 fHistRhoVSRCPt(0),
54 fHistDeltaPtRCvsEP(0),
55 fHistDeltaPtRCExLJ(0),
56 fHistDeltaPtRCExPartialLJ(0),
57 fHistDeltaPtRCRand(0),
58 fHistEmbJetsPtArea(0),
59 fHistEmbJetsCorrPtArea(0),
60 fHistEmbPartPtvsJetPt(0),
61 fHistEmbPartPtvsJetCorrPt(0),
62 fHistJetPtvsJetCorrPt(0),
63 fHistDistLeadPart2JetAxis(0),
64 fHistEmbBkgArea(0),
65 fHistRhoVSEmbBkg(0),
66 fHistDeltaPtEmbArea(0),
67 fHistDeltaPtEmbvsEP(0),
68 fHistRCPtExLJVSDPhiLJ(0),
69 fHistRCPtExPartialLJVSDPhiLJ(0),
70 fHistEmbJetsPhiEta(0),
71 fHistLeadPartPhiEta(0)
72{
73 // Default constructor.
74
75 fHistRCPt = 0;
76 fHistRCPtExLJ = 0;
77 fHistRCPtExPartialLJ = 0;
78 fHistRCPtRand = 0;
79 fHistRhoVSRCPt = 0;
80 fHistDeltaPtRCvsEP = 0;
81 fHistDeltaPtRCExLJ = 0;
82 fHistDeltaPtRCExPartialLJ = 0;
83 fHistDeltaPtRCRand = 0;
84 fHistEmbJetsPtArea = 0;
85 fHistEmbJetsCorrPtArea = 0;
86 fHistEmbPartPtvsJetPt = 0;
87 fHistEmbPartPtvsJetCorrPt = 0;
88 fHistJetPtvsJetCorrPt = 0;
89 fHistDistLeadPart2JetAxis = 0;
90 fHistEmbBkgArea = 0;
91 fHistRhoVSEmbBkg = 0;
92 fHistDeltaPtEmbArea = 0;
93 fHistDeltaPtEmbvsEP = 0;
94
95 SetMakeGeneralHistograms(kTRUE);
96}
97
98//________________________________________________________________________
99AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) :
100 AliAnalysisTaskEmcalJet(name, kTRUE),
101 fMCJetPtThreshold(1),
102 fMinRC2LJ(-1),
103 fRCperEvent(-1),
104 fConeRadius(0.2),
105 fConeMinEta(-0.9),
106 fConeMaxEta(0.9),
107 fConeMinPhi(0),
108 fConeMaxPhi(TMath::Pi()*2),
109 fJetsCont(0),
110 fTracksCont(0),
111 fCaloClustersCont(0),
112 fEmbJetsCont(0),
113 fEmbTracksCont(0),
114 fEmbCaloClustersCont(0),
115 fRandTracksCont(0),
116 fRandCaloClustersCont(0),
117 fHistRCPhiEta(0),
118 fHistRCPt(0),
119 fHistRCPtExLJ(0),
120 fHistRCPtExPartialLJ(0),
121 fHistRCPtRand(0),
122 fHistRhoVSRCPt(0),
123 fHistDeltaPtRCvsEP(0),
124 fHistDeltaPtRCExLJ(0),
125 fHistDeltaPtRCExPartialLJ(0),
126 fHistDeltaPtRCRand(0),
127 fHistEmbJetsPtArea(0),
128 fHistEmbJetsCorrPtArea(0),
129 fHistEmbPartPtvsJetPt(0),
130 fHistEmbPartPtvsJetCorrPt(0),
131 fHistJetPtvsJetCorrPt(0),
132 fHistDistLeadPart2JetAxis(0),
133 fHistEmbBkgArea(0),
134 fHistRhoVSEmbBkg(0),
135 fHistDeltaPtEmbArea(0),
136 fHistDeltaPtEmbvsEP(0),
137 fHistRCPtExLJVSDPhiLJ(0),
138 fHistRCPtExPartialLJVSDPhiLJ(0),
139 fHistEmbJetsPhiEta(0),
140 fHistLeadPartPhiEta(0)
141{
142 // Standard constructor.
143
144 fHistRCPt = 0;
145 fHistRCPtExLJ = 0;
146 fHistRCPtExPartialLJ = 0;
147 fHistRCPtRand = 0;
148 fHistRhoVSRCPt = 0;
149 fHistDeltaPtRCvsEP = 0;
150 fHistDeltaPtRCExLJ = 0;
151 fHistDeltaPtRCExPartialLJ = 0;
152 fHistDeltaPtRCRand = 0;
153 fHistEmbJetsPtArea = 0;
154 fHistEmbJetsCorrPtArea = 0;
155 fHistEmbPartPtvsJetPt = 0;
156 fHistEmbPartPtvsJetCorrPt = 0;
157 fHistJetPtvsJetCorrPt = 0;
158 fHistDistLeadPart2JetAxis = 0;
159 fHistEmbBkgArea = 0;
160 fHistRhoVSEmbBkg = 0;
161 fHistDeltaPtEmbArea = 0;
162 fHistDeltaPtEmbvsEP = 0;
163
164 SetMakeGeneralHistograms(kTRUE);
165}
166
167//________________________________________________________________________
168void AliAnalysisTaskDeltaPt::AllocateHistogramArrays()
169{
170 fHistRCPt = new TH1*[fNcentBins];
171 fHistRCPtExLJ = new TH1*[fNcentBins];
172 fHistRCPtExPartialLJ = new TH1*[fNcentBins];
173 fHistRCPtRand = new TH1*[fNcentBins];
174 fHistRhoVSRCPt = new TH2*[fNcentBins];
175 fHistDeltaPtRCvsEP = new TH2*[fNcentBins];
176 fHistDeltaPtRCExLJ = new TH1*[fNcentBins];
177 fHistDeltaPtRCExPartialLJ = new TH1*[fNcentBins];
178 fHistDeltaPtRCRand = new TH1*[fNcentBins];
179 fHistEmbJetsPtArea = new TH3*[fNcentBins];
180 fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
181 fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
182 fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
183 fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
184 fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
185 fHistEmbBkgArea = new TH2*[fNcentBins];
186 fHistRhoVSEmbBkg = new TH2*[fNcentBins];
187 fHistDeltaPtEmbArea = new TH2*[fNcentBins];
188 fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
189
190 for (Int_t i = 0; i < fNcentBins; i++) {
191 fHistRCPt[i] = 0;
192 fHistRCPtExLJ[i] = 0;
193 fHistRCPtExPartialLJ[i] = 0;
194 fHistRCPtRand[i] = 0;
195 fHistRhoVSRCPt[i] = 0;
196 fHistDeltaPtRCvsEP[i] = 0;
197 fHistDeltaPtRCExLJ[i] = 0;
198 fHistDeltaPtRCExPartialLJ[i] = 0;
199 fHistDeltaPtRCRand[i] = 0;
200 fHistEmbJetsPtArea[i] = 0;
201 fHistEmbJetsCorrPtArea[i] = 0;
202 fHistEmbPartPtvsJetPt[i] = 0;
203 fHistEmbPartPtvsJetCorrPt[i] = 0;
204 fHistJetPtvsJetCorrPt[i] = 0;
205 fHistDistLeadPart2JetAxis[i] = 0;
206 fHistEmbBkgArea[i] = 0;
207 fHistRhoVSEmbBkg[i] = 0;
208 fHistDeltaPtEmbArea[i] = 0;
209 fHistDeltaPtEmbvsEP[i] = 0;
210 }
211}
212
213//________________________________________________________________________
214void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
215{
216 // Create user output.
217
218 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
219
220 AllocateHistogramArrays();
221
222 fJetsCont = GetJetContainer("Jets");
223 fTracksCont = GetParticleContainer("Tracks");
224 fCaloClustersCont = GetClusterContainer("CaloClusters");
225 fEmbJetsCont = GetJetContainer("EmbJets");
226 fEmbTracksCont = GetParticleContainer("EmbTracks");
227 fEmbCaloClustersCont = GetClusterContainer("EmbCaloClusters");
228 fRandTracksCont = GetParticleContainer("RandTracks");
229 fRandCaloClustersCont = GetClusterContainer("RandCaloClusters");
230
231 if (fTracksCont || fCaloClustersCont) {
232 fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
233 fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
234 fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
235 fOutput->Add(fHistRCPhiEta);
236
237 if (fJetsCont) {
238 fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
239 fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
240 fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
241 fOutput->Add(fHistRCPtExLJVSDPhiLJ);
242
243 fHistRCPtExPartialLJVSDPhiLJ = new TH2F("fHistRCPtExPartialLJVSDPhiLJ","fHistRCPtExPartialLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
244 fHistRCPtExPartialLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
245 fHistRCPtExPartialLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
246 fOutput->Add(fHistRCPtExPartialLJVSDPhiLJ);
247 }
248 }
249
250 if (fEmbJetsCont) {
251 fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
252 fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
253 fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
254 fOutput->Add(fHistEmbJetsPhiEta);
255
256 fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
257 fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
258 fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
259 fOutput->Add(fHistLeadPartPhiEta);
260 }
261
262 TString histname;
263
264 const Int_t nbinsZ = 12;
265 Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
266
267 Float_t *binsPt = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
268 Float_t *binsCorrPt = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
269 Float_t *binsArea = GenerateFixedBinArray(50, 0, 2);
270
271 for (Int_t i = 0; i < fNcentBins; i++) {
272 if (fTracksCont || fCaloClustersCont) {
273 histname = "fHistRCPt_";
274 histname += i;
275 fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
276 fHistRCPt[i]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
277 fHistRCPt[i]->GetYaxis()->SetTitle("counts");
278 fOutput->Add(fHistRCPt[i]);
279
280 histname = "fHistRhoVSRCPt_";
281 histname += i;
282 fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
283 fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
284 fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
285 fOutput->Add(fHistRhoVSRCPt[i]);
286
287 histname = "fHistDeltaPtRCvsEP_";
288 histname += i;
289 fHistDeltaPtRCvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
290 fHistDeltaPtRCvsEP[i]->GetXaxis()->SetTitle("#phi_{RC} - #psi_{RP}");
291 fHistDeltaPtRCvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
292 fHistDeltaPtRCvsEP[i]->GetZaxis()->SetTitle("counts");
293 fOutput->Add(fHistDeltaPtRCvsEP[i]);
294
295 if (fJetsCont) {
296 histname = "fHistRCPtExLJ_";
297 histname += i;
298 fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
299 fHistRCPtExLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
300 fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
301 fOutput->Add(fHistRCPtExLJ[i]);
302
303 histname = "fHistDeltaPtRCExLJ_";
304 histname += i;
305 fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
306 fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
307 fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
308 fOutput->Add(fHistDeltaPtRCExLJ[i]);
309
310 histname = "fHistRCPtExPartialLJ_";
311 histname += i;
312 fHistRCPtExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
313 fHistRCPtExPartialLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
314 fHistRCPtExPartialLJ[i]->GetYaxis()->SetTitle("counts");
315 fOutput->Add(fHistRCPtExPartialLJ[i]);
316
317 histname = "fHistDeltaPtRCExPartialLJ_";
318 histname += i;
319 fHistDeltaPtRCExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
320 fHistDeltaPtRCExPartialLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
321 fHistDeltaPtRCExPartialLJ[i]->GetYaxis()->SetTitle("counts");
322 fOutput->Add(fHistDeltaPtRCExPartialLJ[i]);
323 }
324 }
325
326 if (fRandTracksCont || fRandCaloClustersCont) {
327 histname = "fHistRCPtRand_";
328 histname += i;
329 fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
330 fHistRCPtRand[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
331 fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
332 fOutput->Add(fHistRCPtRand[i]);
333
334 histname = "fHistDeltaPtRCRand_";
335 histname += i;
336 fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
337 fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
338 fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
339 fOutput->Add(fHistDeltaPtRCRand[i]);
340 }
341
342 if (fEmbJetsCont) {
343 histname = "fHistEmbJetsPtArea_";
344 histname += i;
345 fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
346 fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
347 fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
348 fOutput->Add(fHistEmbJetsPtArea[i]);
349
350 histname = "fHistEmbJetsCorrPtArea_";
351 histname += i;
352 fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
353 fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
354 fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
355 fOutput->Add(fHistEmbJetsCorrPtArea[i]);
356
357 histname = "fHistEmbPartPtvsJetPt_";
358 histname += i;
359 fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
360 fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
361 fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
362 fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
363 fOutput->Add(fHistEmbPartPtvsJetPt[i]);
364
365 histname = "fHistEmbPartPtvsJetCorrPt_";
366 histname += i;
367 fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
368 fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
369 fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
370 fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
371 fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
372 fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
373
374 histname = "fHistJetPtvsJetCorrPt_";
375 histname += i;
376 fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
377 fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
378 fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
379 fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
380 fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
381 fOutput->Add(fHistJetPtvsJetCorrPt[i]);
382
383 histname = "fHistDistLeadPart2JetAxis_";
384 histname += i;
385 fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
386 fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
387 fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
388 fOutput->Add(fHistDistLeadPart2JetAxis[i]);
389
390 histname = "fHistEmbBkgArea_";
391 histname += i;
392 fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
393 fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
394 fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
395 fOutput->Add(fHistEmbBkgArea[i]);
396
397 histname = "fHistRhoVSEmbBkg_";
398 histname += i;
399 fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
400 fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
401 fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
402 fOutput->Add(fHistRhoVSEmbBkg[i]);
403
404 histname = "fHistDeltaPtEmbArea_";
405 histname += i;
406 fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
407 50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
408 fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
409 fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
410 fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
411 fOutput->Add(fHistDeltaPtEmbArea[i]);
412
413 histname = "fHistDeltaPtEmbvsEP_";
414 histname += i;
415 fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(),
416 402, -TMath::Pi()*1.01, TMath::Pi()*3.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
417 fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
418 fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
419 fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
420 fOutput->Add(fHistDeltaPtEmbvsEP[i]);
421 }
422 }
423
424 delete[] binsPt;
425 delete[] binsCorrPt;
426 delete[] binsArea;
427
428 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
429}
430
431//________________________________________________________________________
432Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
433{
434 // Fill histograms.
435
436 // ************
437 // Random cones
438 // _________________________________
439
440 const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
441 Float_t RCpt = 0;
442 Float_t RCeta = 0;
443 Float_t RCphi = 0;
444
445 if (fTracksCont || fCaloClustersCont) {
446
447 for (Int_t i = 0; i < fRCperEvent; i++) {
448 // Simple random cones
449 RCpt = 0;
450 RCeta = 0;
451 RCphi = 0;
452 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
453 if (RCpt > 0) {
454 fHistRCPhiEta->Fill(RCeta, RCphi);
455 fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
456
457 fHistRCPt[fCentBin]->Fill(RCpt);
458
459 Double_t ep = RCphi - fEPV0;
460 while (ep < 0) ep += TMath::Pi();
461 while (ep >= TMath::Pi()) ep -= TMath::Pi();
462
463 fHistDeltaPtRCvsEP[fCentBin]->Fill(ep, RCpt - rcArea * fRhoVal);
464 }
465
466 if (fJetsCont) {
467
468 // Random cones far from leading jet
469 AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
470
471 RCpt = 0;
472 RCeta = 0;
473 RCphi = 0;
474 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
475 if (RCpt > 0) {
476 if (jet) {
477 Float_t dphi = RCphi - jet->Phi();
478 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
479 if (dphi < -1.6) dphi += TMath::Pi() * 2;
480 fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
481 }
482 fHistRCPtExLJ[fCentBin]->Fill(RCpt);
483 fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
484 }
485
486 //partial exclusion
487 if(fBeamType == kpA) {
488
489 RCpt = 0;
490 RCeta = 0;
491 RCphi = 0;
492 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet, kTRUE);
493
494 if (RCpt > 0) {
495 if (jet) {
496 Float_t dphi = RCphi - jet->Phi();
497 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
498 if (dphi < -1.6) dphi += TMath::Pi() * 2;
499 fHistRCPtExPartialLJVSDPhiLJ->Fill(RCpt, dphi);
500 }
501 fHistRCPtExPartialLJ[fCentBin]->Fill(RCpt);
502 fHistDeltaPtRCExPartialLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
503 }
504 }
505 }
506 }
507 }
508
509 // Random cones with randomized particles
510 if (fRandTracksCont || fRandCaloClustersCont) {
511 RCpt = 0;
512 RCeta = 0;
513 RCphi = 0;
514 GetRandomCone(RCpt, RCeta, RCphi, fRandTracksCont, fRandCaloClustersCont, 0);
515 if (RCpt > 0) {
516 fHistRCPtRand[fCentBin]->Fill(RCpt);
517 fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
518 }
519 }
520
521 // ************
522 // Embedding
523 // _________________________________
524
525 if (fEmbJetsCont) {
526
527 AliEmcalJet *embJet = NextEmbeddedJet(kTRUE);
528
529 while (embJet != 0) {
530 TLorentzVector mom;
531 fEmbJetsCont->GetLeadingHadronMomentum(mom,embJet);
532
533 Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - mom.Eta()) * (embJet->Eta() - mom.Eta()) + (embJet->Phi() - mom.Phi()) * (embJet->Phi() - mom.Phi()));
534
535 fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
536 fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
537 fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi());
538 fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
539
540 fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), mom.Pt());
541 fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), mom.Pt());
542 fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
543 fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
544
545 fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
546 fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
547 fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
548 fHistDeltaPtEmbvsEP[fCentBin]->Fill(embJet->Phi() - fEPV0, embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
549
550 embJet = NextEmbeddedJet();
551 }
552 }
553
554 return kTRUE;
555}
556
557//________________________________________________________________________
558AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Bool_t reset)
559{
560 // Get the next accepted embedded jet.
561
562 Int_t i = reset ? 0 : -1;
563
564 AliEmcalJet* jet = fEmbJetsCont->GetNextAcceptJet(i);
565 while (jet && jet->MCPt() < fMCJetPtThreshold) jet = fEmbJetsCont->GetNextAcceptJet();
566
567 return jet;
568}
569
570//________________________________________________________________________
571void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
572 AliParticleContainer* tracks, AliClusterContainer* clusters,
573 AliEmcalJet *jet, Bool_t bPartialExclusion) const
574{
575 // Get rigid cone.
576
577 eta = -999;
578 phi = -999;
579 pt = 0;
580
581 if (!tracks && !clusters)
582 return;
583
584 Float_t LJeta = 999;
585 Float_t LJphi = 999;
586
587 if (jet) {
588 LJeta = jet->Eta();
589 LJphi = jet->Phi();
590 }
591
592 Float_t maxEta = fConeMaxEta;
593 Float_t minEta = fConeMinEta;
594 Float_t maxPhi = fConeMaxPhi;
595 Float_t minPhi = fConeMinPhi;
596
597 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
598 if (minPhi < 0) minPhi = 0;
599
600 Float_t dLJ = 0;
601 Int_t repeats = 0;
602 Bool_t reject = kTRUE;
603 do {
604 eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
605 phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
606 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
607
608 if(bPartialExclusion) {
609 reject = kFALSE;
610
611 TRandom3 rnd;
612 rnd.SetSeed(0);
613
614 Double_t ncoll = GetNColl();
615
616 Double_t prob = 0.;
617 if(ncoll>0)
618 prob = 1./ncoll;
619
620 if(rnd.Rndm()<=prob) reject = kTRUE; //reject cone
621 }
622
623 repeats++;
624 } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
625
626 if (repeats == 999) {
627 AliWarning(Form("%s: Could not get random cone!", GetName()));
628 return;
629 }
630
631 if (clusters) {
632 AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
633 while (cluster) {
634 TLorentzVector nPart;
635 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
636
637 Float_t cluseta = nPart.Eta();
638 Float_t clusphi = nPart.Phi();
639
640 if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
641 clusphi += 2 * TMath::Pi();
642 if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
643 clusphi -= 2 * TMath::Pi();
644
645 Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
646 if (d <= fConeRadius)
647 pt += nPart.Pt();
648
649 cluster = clusters->GetNextAcceptCluster();
650 }
651 }
652
653 if (tracks) {
654 AliVParticle* track = tracks->GetNextAcceptParticle(0);
655 while(track) {
656 Float_t tracketa = track->Eta();
657 Float_t trackphi = track->Phi();
658
659 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
660 trackphi += 2 * TMath::Pi();
661 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
662 trackphi -= 2 * TMath::Pi();
663
664 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
665 if (d <= fConeRadius)
666 pt += track->Pt();
667
668 track = tracks->GetNextAcceptParticle();
669 }
670 }
671}
672
673//________________________________________________________________________
674void AliAnalysisTaskDeltaPt::SetConeEtaPhiEMCAL()
675{
676 // Set default cuts for full cones
677
678 SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
679 SetConePhiLimits(1.4+fConeRadius,TMath::Pi()-fConeRadius);
680}
681
682//________________________________________________________________________
683void AliAnalysisTaskDeltaPt::SetConeEtaPhiTPC()
684{
685 // Set default cuts for charged cones
686
687 SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
688 SetConePhiLimits(-10, 10);
689}
690
691//________________________________________________________________________
692void AliAnalysisTaskDeltaPt::ExecOnce()
693{
694 // Initialize the analysis.
695
696 AliAnalysisTaskEmcalJet::ExecOnce();
697
698 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
699 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
700 if (fEmbTracksCont && fEmbTracksCont->GetArray() == 0) fEmbTracksCont = 0;
701 if (fEmbCaloClustersCont && fEmbCaloClustersCont->GetArray() == 0) fEmbCaloClustersCont = 0;
702 if (fRandTracksCont && fRandTracksCont->GetArray() == 0) fRandTracksCont = 0;
703 if (fRandCaloClustersCont && fRandCaloClustersCont->GetArray() == 0) fRandCaloClustersCont = 0;
704 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
705 if (fEmbJetsCont && fEmbJetsCont->GetArray() == 0) fEmbJetsCont = 0;
706
707 if (fRCperEvent < 0) {
708 Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
709 Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
710 fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
711 if (fRCperEvent == 0)
712 fRCperEvent = 1;
713 }
714
715 if (fMinRC2LJ < 0)
716 fMinRC2LJ = fConeRadius * 1.5;
717
718 const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
719 if (fMinRC2LJ > maxDist) {
720 AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
721 "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
722 fMinRC2LJ = maxDist;
723 }
724}
725
726//________________________________________________________________________
727Double_t AliAnalysisTaskDeltaPt::GetNColl() const {
728 // Get NColl - returns value of corresponding bin
729 // only works for pA
730 // values taken from V0A slicing https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Tables_with_centrality_bins_from
731
732 if(fBeamType == kpA) {
733
734 const Int_t nNCollBins = 7;
735
736 Double_t centMin[nNCollBins] = {0.,5.,10.,20.,40.,60.,80.};
737 Double_t centMax[nNCollBins] = {5.,10.,20.,40.,60.,80.,100.};
738
739 Double_t nColl[nNCollBins] = {14.7,13.,11.7,9.38,6.49,3.96,1.52};
740
741 for(Int_t i = 0; i<nNCollBins; i++) {
742 if(fCent>=centMin[i] && fCent<centMax[i])
743 return nColl[i];
744 }
745
746 return -1.;
747 }
748 else {
749 AliWarning(Form("%s: Only works for pA analysis. Returning -1",GetName()));
750 return -1.;
751 }
752}