]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx
timing calibration for 2013 period added
[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() :
6421eeb0 31 AliAnalysisTaskEmcalJetDev("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) :
6421eeb0 122 AliAnalysisTaskEmcalJetDev(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
6421eeb0 216 AliAnalysisTaskEmcalJetDev::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;
faf21244 285 fHistDeltaPtRCvsEP[i] = new TH2F(histname.Data(), histname.Data(), 402, -TMath::Pi()*1.01, TMath::Pi()*3.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
286 fHistDeltaPtRCvsEP[i]->GetXaxis()->SetTitle("#phi_{RC} - #Psi_{EP}");
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);
faf21244 454 fHistDeltaPtRCvsEP[fCentBin]->Fill(RCphi - fEPV0, RCpt - rcArea * fRhoVal);
6a20534a 455 }
6421eeb0 456
457 if (fJetsCont) {
6a20534a 458
459 // Random cones far from leading jet
6421eeb0 460 AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
6a20534a 461
462 RCpt = 0;
463 RCeta = 0;
464 RCphi = 0;
6421eeb0 465 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
6a20534a 466 if (RCpt > 0) {
467 if (jet) {
468 Float_t dphi = RCphi - jet->Phi();
469 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
470 if (dphi < -1.6) dphi += TMath::Pi() * 2;
471 fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
472 }
473 fHistRCPtExLJ[fCentBin]->Fill(RCpt);
474 fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
475 }
56f370b0 476
477 //partial exclusion
478 if(fBeamType == kpA) {
479
480 RCpt = 0;
481 RCeta = 0;
482 RCphi = 0;
6421eeb0 483 GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet, kTRUE);
56f370b0 484
485 if (RCpt > 0) {
486 if (jet) {
487 Float_t dphi = RCphi - jet->Phi();
488 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
489 if (dphi < -1.6) dphi += TMath::Pi() * 2;
490 fHistRCPtExPartialLJVSDPhiLJ->Fill(RCpt, dphi);
491 }
492 fHistRCPtExPartialLJ[fCentBin]->Fill(RCpt);
493 fHistDeltaPtRCExPartialLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
494 }
495 }
a487deae 496 }
a487deae 497 }
6a20534a 498 }
499
500 // Random cones with randomized particles
6421eeb0 501 if (fRandTracksCont || fRandCaloClustersCont) {
a487deae 502 RCpt = 0;
503 RCeta = 0;
504 RCphi = 0;
6421eeb0 505 GetRandomCone(RCpt, RCeta, RCphi, fRandTracksCont, fRandCaloClustersCont, 0);
a487deae 506 if (RCpt > 0) {
1f9c287f 507 fHistRCPtRand[fCentBin]->Fill(RCpt);
a487deae 508 fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
509 }
510 }
6a20534a 511
a487deae 512 // ************
513 // Embedding
514 // _________________________________
515
6421eeb0 516 if (fEmbJetsCont) {
6a20534a 517
6421eeb0 518 AliEmcalJet *embJet = NextEmbeddedJet(kTRUE);
6a20534a 519
520 while (embJet != 0) {
6421eeb0 521 TLorentzVector mom;
522 fEmbJetsCont->GetLeadingHadronMomentum(mom,embJet);
6a20534a 523
6421eeb0 524 Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - mom.Eta()) * (embJet->Eta() - mom.Eta()) + (embJet->Phi() - mom.Phi()) * (embJet->Phi() - mom.Phi()));
6a20534a 525
526 fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
ca5c29fa 527 fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
6421eeb0 528 fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi());
6a20534a 529 fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
530
6421eeb0 531 fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), mom.Pt());
532 fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), mom.Pt());
6a20534a 533 fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
2103dc6a 534 fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
6a20534a 535
536 fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
537 fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
538 fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
faf21244 539 fHistDeltaPtEmbvsEP[fCentBin]->Fill(embJet->Phi() - fEPV0, embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
f660c2d6 540
6a20534a 541 embJet = NextEmbeddedJet();
a487deae 542 }
a487deae 543 }
544
545 return kTRUE;
546}
547
548//________________________________________________________________________
6421eeb0 549AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Bool_t reset)
a487deae 550{
6421eeb0 551 // Get the next accepted embedded jet.
f660c2d6 552
6421eeb0 553 Int_t i = reset ? 0 : -1;
a487deae 554
6421eeb0 555 AliEmcalJet* jet = fEmbJetsCont->GetNextAcceptJet(i);
556 while (jet && jet->MCPt() < fMCJetPtThreshold) jet = fEmbJetsCont->GetNextAcceptJet();
f660c2d6 557
6421eeb0 558 return jet;
a487deae 559}
560
561//________________________________________________________________________
562void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
6421eeb0 563 AliParticleContainer* tracks, AliClusterContainer* clusters,
564 AliEmcalJet *jet, Bool_t bPartialExclusion) const
a487deae 565{
566 // Get rigid cone.
567
6421eeb0 568 eta = -999;
569 phi = -999;
a487deae 570 pt = 0;
571
572 if (!tracks && !clusters)
573 return;
574
575 Float_t LJeta = 999;
576 Float_t LJphi = 999;
577
578 if (jet) {
579 LJeta = jet->Eta();
580 LJphi = jet->Phi();
581 }
582
6421eeb0 583 Float_t maxEta = fConeMaxEta;
584 Float_t minEta = fConeMinEta;
585 Float_t maxPhi = fConeMaxPhi;
586 Float_t minPhi = fConeMinPhi;
a487deae 587
588 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
589 if (minPhi < 0) minPhi = 0;
590
591 Float_t dLJ = 0;
592 Int_t repeats = 0;
56f370b0 593 Bool_t reject = kTRUE;
a487deae 594 do {
595 eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
596 phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
597 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
56f370b0 598
599 if(bPartialExclusion) {
600 reject = kFALSE;
601
602 TRandom3 rnd;
603 rnd.SetSeed(0);
604
605 Double_t ncoll = GetNColl();
606
607 Double_t prob = 0.;
608 if(ncoll>0)
609 prob = 1./ncoll;
610
611 if(rnd.Rndm()<=prob) reject = kTRUE; //reject cone
612 }
613
a487deae 614 repeats++;
56f370b0 615 } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
a487deae 616
617 if (repeats == 999) {
618 AliWarning(Form("%s: Could not get random cone!", GetName()));
619 return;
620 }
621
622 if (clusters) {
6421eeb0 623 AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
624 while (cluster) {
a487deae 625 TLorentzVector nPart;
626 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
a487deae 627
6421eeb0 628 Float_t cluseta = nPart.Eta();
629 Float_t clusphi = nPart.Phi();
630
631 if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
632 clusphi += 2 * TMath::Pi();
633 if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
634 clusphi -= 2 * TMath::Pi();
635
636 Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
637 if (d <= fConeRadius)
a487deae 638 pt += nPart.Pt();
6421eeb0 639
640 cluster = clusters->GetNextAcceptCluster();
a487deae 641 }
642 }
643
644 if (tracks) {
6421eeb0 645 AliVParticle* track = tracks->GetNextAcceptParticle(0);
646 while(track) {
a487deae 647 Float_t tracketa = track->Eta();
648 Float_t trackphi = track->Phi();
649
650 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
651 trackphi += 2 * TMath::Pi();
652 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
653 trackphi -= 2 * TMath::Pi();
654
655 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
6421eeb0 656 if (d <= fConeRadius)
a487deae 657 pt += track->Pt();
6421eeb0 658
659 track = tracks->GetNextAcceptParticle();
a487deae 660 }
661 }
662}
663
6421eeb0 664//________________________________________________________________________
665void AliAnalysisTaskDeltaPt::SetConeEtaPhiEMCAL()
666{
667 // Set default cuts for full cones
668
669 SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
670 SetConePhiLimits(1.4+fConeRadius,TMath::Pi()-fConeRadius);
671}
672
673//________________________________________________________________________
674void AliAnalysisTaskDeltaPt::SetConeEtaPhiTPC()
675{
676 // Set default cuts for charged cones
677
678 SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
679 SetConePhiLimits(-10, 10);
680}
681
a487deae 682//________________________________________________________________________
683void AliAnalysisTaskDeltaPt::ExecOnce()
684{
685 // Initialize the analysis.
686
6421eeb0 687 AliAnalysisTaskEmcalJetDev::ExecOnce();
a487deae 688
6421eeb0 689 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
690 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
691 if (fEmbTracksCont && fEmbTracksCont->GetArray() == 0) fEmbTracksCont = 0;
692 if (fEmbCaloClustersCont && fEmbCaloClustersCont->GetArray() == 0) fEmbCaloClustersCont = 0;
693 if (fRandTracksCont && fRandTracksCont->GetArray() == 0) fRandTracksCont = 0;
694 if (fRandCaloClustersCont && fRandCaloClustersCont->GetArray() == 0) fRandCaloClustersCont = 0;
7cd832c7 695 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
696 if (fEmbJetsCont && fEmbJetsCont->GetArray() == 0) fEmbJetsCont = 0;
a487deae 697
a487deae 698 if (fRCperEvent < 0) {
6421eeb0 699 Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
700 Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
701 fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
a487deae 702 if (fRCperEvent == 0)
703 fRCperEvent = 1;
704 }
705
706 if (fMinRC2LJ < 0)
6421eeb0 707 fMinRC2LJ = fConeRadius * 1.5;
a487deae 708
6421eeb0 709 const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
a487deae 710 if (fMinRC2LJ > maxDist) {
711 AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
712 "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
713 fMinRC2LJ = maxDist;
714 }
715}
716
56f370b0 717//________________________________________________________________________
718Double_t AliAnalysisTaskDeltaPt::GetNColl() const {
719 // Get NColl - returns value of corresponding bin
720 // only works for pA
721 // values taken from V0A slicing https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Tables_with_centrality_bins_from
722
723 if(fBeamType == kpA) {
724
725 const Int_t nNCollBins = 7;
726
727 Double_t centMin[nNCollBins] = {0.,5.,10.,20.,40.,60.,80.};
728 Double_t centMax[nNCollBins] = {5.,10.,20.,40.,60.,80.,100.};
729
730 Double_t nColl[nNCollBins] = {14.7,13.,11.7,9.38,6.49,3.96,1.52};
731
732 for(Int_t i = 0; i<nNCollBins; i++) {
733 if(fCent>=centMin[i] && fCent<centMax[i])
734 return nColl[i];
735 }
736
737 return -1.;
738 }
739 else {
740 AliWarning(Form("%s: Only works for pA analysis. Returning -1",GetName()));
741 return -1.;
742 }
a487deae 743}