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