]> 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
Added histograms and the possibility to save or not the tree
[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),
56f370b0 48 fHistRCPhiEta(0),
49 fHistRCPt(0),
50 fHistRCPtExLJ(0),
51 fHistRCPtExPartialLJ(0),
52 fHistRCPtRand(0),
53 fHistRhoVSRCPt(0),
faf21244 54 fHistDeltaPtRCvsEP(0),
56f370b0 55 fHistDeltaPtRCExLJ(0),
56 fHistDeltaPtRCExPartialLJ(0),
57 fHistDeltaPtRCRand(0),
56f370b0 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),
faf21244 67 fHistDeltaPtEmbvsEP(0),
a487deae 68 fHistRCPtExLJVSDPhiLJ(0),
56f370b0 69 fHistRCPtExPartialLJVSDPhiLJ(0),
6a20534a 70 fHistEmbJetsPhiEta(0),
71 fHistLeadPartPhiEta(0)
a487deae 72{
73 // Default constructor.
74
ed062c95 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;
a487deae 94
95 SetMakeGeneralHistograms(kTRUE);
96}
97
98//________________________________________________________________________
99AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) :
9239b066 100 AliAnalysisTaskEmcalJet(name, kTRUE),
f660c2d6 101 fMCJetPtThreshold(1),
a487deae 102 fMinRC2LJ(-1),
a487deae 103 fRCperEvent(-1),
6421eeb0 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),
56f370b0 117 fHistRCPhiEta(0),
118 fHistRCPt(0),
119 fHistRCPtExLJ(0),
120 fHistRCPtExPartialLJ(0),
121 fHistRCPtRand(0),
122 fHistRhoVSRCPt(0),
faf21244 123 fHistDeltaPtRCvsEP(0),
56f370b0 124 fHistDeltaPtRCExLJ(0),
125 fHistDeltaPtRCExPartialLJ(0),
126 fHistDeltaPtRCRand(0),
56f370b0 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),
faf21244 136 fHistDeltaPtEmbvsEP(0),
a487deae 137 fHistRCPtExLJVSDPhiLJ(0),
56f370b0 138 fHistRCPtExPartialLJVSDPhiLJ(0),
6a20534a 139 fHistEmbJetsPhiEta(0),
140 fHistLeadPartPhiEta(0)
a487deae 141{
142 // Standard constructor.
143
ed062c95 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{
56f370b0 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];
faf21244 175 fHistDeltaPtRCvsEP = new TH2*[fNcentBins];
56f370b0 176 fHistDeltaPtRCExLJ = new TH1*[fNcentBins];
177 fHistDeltaPtRCExPartialLJ = new TH1*[fNcentBins];
178 fHistDeltaPtRCRand = new TH1*[fNcentBins];
56f370b0 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];
faf21244 188 fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
56f370b0 189
190 for (Int_t i = 0; i < fNcentBins; i++) {
a487deae 191 fHistRCPt[i] = 0;
192 fHistRCPtExLJ[i] = 0;
56f370b0 193 fHistRCPtExPartialLJ[i] = 0;
194 fHistRCPtRand[i] = 0;
59f16b27 195 fHistRhoVSRCPt[i] = 0;
faf21244 196 fHistDeltaPtRCvsEP[i] = 0;
a487deae 197 fHistDeltaPtRCExLJ[i] = 0;
56f370b0 198 fHistDeltaPtRCExPartialLJ[i] = 0;
a487deae 199 fHistDeltaPtRCRand[i] = 0;
1f9c287f 200 fHistEmbJetsPtArea[i] = 0;
201 fHistEmbJetsCorrPtArea[i] = 0;
6a20534a 202 fHistEmbPartPtvsJetPt[i] = 0;
ca5c29fa 203 fHistEmbPartPtvsJetCorrPt[i] = 0;
2103dc6a 204 fHistJetPtvsJetCorrPt[i] = 0;
6a20534a 205 fHistDistLeadPart2JetAxis[i] = 0;
1f9c287f 206 fHistEmbBkgArea[i] = 0;
59f16b27 207 fHistRhoVSEmbBkg[i] = 0;
1f9c287f 208 fHistDeltaPtEmbArea[i] = 0;
faf21244 209 fHistDeltaPtEmbvsEP[i] = 0;
a487deae 210 }
a487deae 211}
212
a487deae 213//________________________________________________________________________
214void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
215{
216 // Create user output.
217
9239b066 218 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
a487deae 219
ed062c95 220 AllocateHistogramArrays();
221
6421eeb0 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) {
5ce8ae64 232 fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
6a20534a 233 fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
234 fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
235 fOutput->Add(fHistRCPhiEta);
236
6421eeb0 237 if (fJetsCont) {
6a20534a 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);
56f370b0 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);
6a20534a 247 }
248 }
249
6421eeb0 250 if (fEmbJetsCont) {
5ce8ae64 251 fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
6a20534a 252 fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
253 fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
254 fOutput->Add(fHistEmbJetsPhiEta);
a487deae 255
5ce8ae64 256 fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
6a20534a 257 fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
258 fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
259 fOutput->Add(fHistLeadPartPhiEta);
a487deae 260 }
261
262 TString histname;
263
ca5c29fa 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);
9adcb46d 269 Float_t *binsArea = GenerateFixedBinArray(50, 0, 2);
ca5c29fa 270
6a20534a 271 for (Int_t i = 0; i < fNcentBins; i++) {
6421eeb0 272 if (fTracksCont || fCaloClustersCont) {
6a20534a 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
faf21244 287 histname = "fHistDeltaPtRCvsEP_";
6a20534a 288 histname += i;
9f6dc7fd 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}");
faf21244 291 fHistDeltaPtRCvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
292 fHistDeltaPtRCvsEP[i]->GetZaxis()->SetTitle("counts");
293 fOutput->Add(fHistDeltaPtRCvsEP[i]);
6a20534a 294
6421eeb0 295 if (fJetsCont) {
6a20534a 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]);
56f370b0 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]);
6a20534a 323 }
324 }
325
6421eeb0 326 if (fRandTracksCont || fRandCaloClustersCont) {
6a20534a 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
6421eeb0 342 if (fEmbJetsCont) {
1f9c287f 343 histname = "fHistEmbJetsPtArea_";
a487deae 344 histname += i;
9adcb46d 345 fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
1f9c287f 346 fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
6a20534a 347 fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
1f9c287f 348 fOutput->Add(fHistEmbJetsPtArea[i]);
a487deae 349
1f9c287f 350 histname = "fHistEmbJetsCorrPtArea_";
a487deae 351 histname += i;
9adcb46d 352 fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
1f9c287f 353 fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
6a20534a 354 fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
1f9c287f 355 fOutput->Add(fHistEmbJetsCorrPtArea[i]);
a487deae 356
6a20534a 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
ca5c29fa 365 histname = "fHistEmbPartPtvsJetCorrPt_";
a487deae 366 histname += i;
9adcb46d 367 fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
368 fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
ca5c29fa 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]);
a487deae 373
2103dc6a 374 histname = "fHistJetPtvsJetCorrPt_";
375 histname += i;
9adcb46d 376 fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
377 fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
2103dc6a 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
6a20534a 383 histname = "fHistDistLeadPart2JetAxis_";
a487deae 384 histname += i;
6a20534a 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]);
a487deae 389
1f9c287f 390 histname = "fHistEmbBkgArea_";
391 histname += i;
9adcb46d 392 fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
1f9c287f 393 fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
6a20534a 394 fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
1f9c287f 395 fOutput->Add(fHistEmbBkgArea[i]);
396
59f16b27 397 histname = "fHistRhoVSEmbBkg_";
398 histname += i;
1f9c287f 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})");
6a20534a 401 fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
59f16b27 402 fOutput->Add(fHistRhoVSEmbBkg[i]);
a487deae 403
1f9c287f 404 histname = "fHistDeltaPtEmbArea_";
a487deae 405 histname += i;
9adcb46d 406 fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
faf21244 407 50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
1f9c287f 408 fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
409 fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
faf21244 410 fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
1f9c287f 411 fOutput->Add(fHistDeltaPtEmbArea[i]);
faf21244 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]);
a487deae 421 }
422 }
423
cac0e10b 424 delete[] binsPt;
425 delete[] binsCorrPt;
426 delete[] binsArea;
427
a487deae 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
a487deae 436 // ************
437 // Random cones
438 // _________________________________
439
6421eeb0 440 const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
6a20534a 441 Float_t RCpt = 0;
442 Float_t RCeta = 0;
443 Float_t RCphi = 0;
a487deae 444
6421eeb0 445 if (fTracksCont || fCaloClustersCont) {
6a20534a 446
447 for (Int_t i = 0; i < fRCperEvent; i++) {
448 // Simple random cones
449 RCpt = 0;
450 RCeta = 0;
451 RCphi = 0;
6421eeb0 452 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
6a20534a 453 if (RCpt > 0) {
454 fHistRCPhiEta->Fill(RCeta, RCphi);
455 fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
456
457 fHistRCPt[fCentBin]->Fill(RCpt);
9f6dc7fd 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);
6a20534a 464 }
6421eeb0 465
466 if (fJetsCont) {
6a20534a 467
468 // Random cones far from leading jet
6421eeb0 469 AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
6a20534a 470
471 RCpt = 0;
472 RCeta = 0;
473 RCphi = 0;
6421eeb0 474 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
6a20534a 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 }
56f370b0 485
486 //partial exclusion
487 if(fBeamType == kpA) {
488
489 RCpt = 0;
490 RCeta = 0;
491 RCphi = 0;
6421eeb0 492 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet, kTRUE);
56f370b0 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 }
a487deae 505 }
a487deae 506 }
6a20534a 507 }
508
509 // Random cones with randomized particles
6421eeb0 510 if (fRandTracksCont || fRandCaloClustersCont) {
a487deae 511 RCpt = 0;
512 RCeta = 0;
513 RCphi = 0;
6421eeb0 514 GetRandomCone(RCpt, RCeta, RCphi, fRandTracksCont, fRandCaloClustersCont, 0);
a487deae 515 if (RCpt > 0) {
1f9c287f 516 fHistRCPtRand[fCentBin]->Fill(RCpt);
a487deae 517 fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
518 }
519 }
6a20534a 520
a487deae 521 // ************
522 // Embedding
523 // _________________________________
524
6421eeb0 525 if (fEmbJetsCont) {
6a20534a 526
6421eeb0 527 AliEmcalJet *embJet = NextEmbeddedJet(kTRUE);
6a20534a 528
529 while (embJet != 0) {
6421eeb0 530 TLorentzVector mom;
531 fEmbJetsCont->GetLeadingHadronMomentum(mom,embJet);
6a20534a 532
6421eeb0 533 Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - mom.Eta()) * (embJet->Eta() - mom.Eta()) + (embJet->Phi() - mom.Phi()) * (embJet->Phi() - mom.Phi()));
6a20534a 534
535 fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
ca5c29fa 536 fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
6421eeb0 537 fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi());
6a20534a 538 fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
539
6421eeb0 540 fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), mom.Pt());
541 fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), mom.Pt());
6a20534a 542 fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
2103dc6a 543 fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
6a20534a 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());
faf21244 548 fHistDeltaPtEmbvsEP[fCentBin]->Fill(embJet->Phi() - fEPV0, embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
f660c2d6 549
6a20534a 550 embJet = NextEmbeddedJet();
a487deae 551 }
a487deae 552 }
553
554 return kTRUE;
555}
556
557//________________________________________________________________________
6421eeb0 558AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Bool_t reset)
a487deae 559{
6421eeb0 560 // Get the next accepted embedded jet.
f660c2d6 561
6421eeb0 562 Int_t i = reset ? 0 : -1;
a487deae 563
6421eeb0 564 AliEmcalJet* jet = fEmbJetsCont->GetNextAcceptJet(i);
565 while (jet && jet->MCPt() < fMCJetPtThreshold) jet = fEmbJetsCont->GetNextAcceptJet();
f660c2d6 566
6421eeb0 567 return jet;
a487deae 568}
569
570//________________________________________________________________________
571void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
6421eeb0 572 AliParticleContainer* tracks, AliClusterContainer* clusters,
573 AliEmcalJet *jet, Bool_t bPartialExclusion) const
a487deae 574{
575 // Get rigid cone.
576
6421eeb0 577 eta = -999;
578 phi = -999;
a487deae 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
6421eeb0 592 Float_t maxEta = fConeMaxEta;
593 Float_t minEta = fConeMinEta;
594 Float_t maxPhi = fConeMaxPhi;
595 Float_t minPhi = fConeMinPhi;
a487deae 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;
56f370b0 602 Bool_t reject = kTRUE;
a487deae 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));
56f370b0 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
a487deae 623 repeats++;
56f370b0 624 } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
a487deae 625
626 if (repeats == 999) {
627 AliWarning(Form("%s: Could not get random cone!", GetName()));
628 return;
629 }
630
631 if (clusters) {
6421eeb0 632 AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
633 while (cluster) {
a487deae 634 TLorentzVector nPart;
635 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
a487deae 636
6421eeb0 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)
a487deae 647 pt += nPart.Pt();
6421eeb0 648
649 cluster = clusters->GetNextAcceptCluster();
a487deae 650 }
651 }
652
653 if (tracks) {
6421eeb0 654 AliVParticle* track = tracks->GetNextAcceptParticle(0);
655 while(track) {
a487deae 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));
6421eeb0 665 if (d <= fConeRadius)
a487deae 666 pt += track->Pt();
6421eeb0 667
668 track = tracks->GetNextAcceptParticle();
a487deae 669 }
670 }
671}
672
6421eeb0 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
a487deae 691//________________________________________________________________________
692void AliAnalysisTaskDeltaPt::ExecOnce()
693{
694 // Initialize the analysis.
695
9239b066 696 AliAnalysisTaskEmcalJet::ExecOnce();
a487deae 697
6421eeb0 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;
7cd832c7 704 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
705 if (fEmbJetsCont && fEmbJetsCont->GetArray() == 0) fEmbJetsCont = 0;
a487deae 706
a487deae 707 if (fRCperEvent < 0) {
6421eeb0 708 Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
709 Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
710 fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
a487deae 711 if (fRCperEvent == 0)
712 fRCperEvent = 1;
713 }
714
715 if (fMinRC2LJ < 0)
6421eeb0 716 fMinRC2LJ = fConeRadius * 1.5;
a487deae 717
6421eeb0 718 const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
a487deae 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
56f370b0 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 }
a487deae 752}