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