]> 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
additional fixes from Salvatore
[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"
21
22#include "AliAnalysisTaskDeltaPt.h"
23
24ClassImp(AliAnalysisTaskDeltaPt)
25
26//________________________________________________________________________
27AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() :
28 AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
f660c2d6 29 fMCJetPtThreshold(1),
a487deae 30 fMinRC2LJ(-1),
31 fEmbJetsName(""),
32 fEmbTracksName(""),
33 fEmbCaloName(""),
34 fRandTracksName("TracksRandomized"),
35 fRandCaloName("CaloClustersRandomized"),
36 fRCperEvent(-1),
37 fEmbJets(0),
38 fEmbTracks(0),
39 fEmbCaloClusters(0),
40 fRandTracks(0),
41 fRandCaloClusters(0),
f660c2d6 42 fEmbeddedClusterNIds(0),
43 fEmbeddedTrackNIds(0),
a487deae 44 fHistRCPhiEta(0),
45 fHistRCPtExLJVSDPhiLJ(0),
6a20534a 46 fHistEmbJetsPhiEta(0),
47 fHistLeadPartPhiEta(0)
a487deae 48{
49 // Default constructor.
50
51 for (Int_t i = 0; i < 4; i++) {
52 fHistRCPt[i] = 0;
53 fHistRCPtExLJ[i] = 0;
54 fHistRCPtRand[i] = 0;
59f16b27 55 fHistRhoVSRCPt[i] = 0;
a487deae 56 fHistDeltaPtRC[i] = 0;
57 fHistDeltaPtRCExLJ[i] = 0;
58 fHistDeltaPtRCRand[i] = 0;
5ce8ae64 59 fHistEmbRejectedJetsPhiEta[i] = 0;
60 fHistEmbRejectedJetsPtArea[i] = 0;
6a20534a 61 fHistEmbNotFoundPt[i] = 0;
a487deae 62 fHistEmbNotFoundPhiEta[i] = 0;
1f9c287f 63 fHistEmbJetsPtArea[i] = 0;
64 fHistEmbJetsCorrPtArea[i] = 0;
6a20534a 65 fHistEmbPartPtvsJetPt[i] = 0;
ca5c29fa 66 fHistEmbPartPtvsJetCorrPt[i] = 0;
2103dc6a 67 fHistJetPtvsJetCorrPt[i] = 0;
6a20534a 68 fHistDistLeadPart2JetAxis[i] = 0;
1f9c287f 69 fHistEmbBkgArea[i] = 0;
59f16b27 70 fHistRhoVSEmbBkg[i] = 0;
1f9c287f 71 fHistDeltaPtEmbArea[i] = 0;
a487deae 72 }
73
74 SetMakeGeneralHistograms(kTRUE);
75}
76
77//________________________________________________________________________
78AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) :
79 AliAnalysisTaskEmcalJet(name, kTRUE),
f660c2d6 80 fMCJetPtThreshold(1),
a487deae 81 fMinRC2LJ(-1),
82 fEmbJetsName(""),
83 fEmbTracksName(""),
84 fEmbCaloName(""),
85 fRandTracksName("TracksRandomized"),
86 fRandCaloName("CaloClustersRandomized"),
87 fRCperEvent(-1),
88 fEmbJets(0),
89 fEmbTracks(0),
90 fEmbCaloClusters(0),
91 fRandTracks(0),
92 fRandCaloClusters(0),
f660c2d6 93 fEmbeddedClusterNIds(0),
94 fEmbeddedTrackNIds(0),
a487deae 95 fHistRCPhiEta(0),
96 fHistRCPtExLJVSDPhiLJ(0),
6a20534a 97 fHistEmbJetsPhiEta(0),
98 fHistLeadPartPhiEta(0)
a487deae 99{
100 // Standard constructor.
101
102 for (Int_t i = 0; i < 4; i++) {
103 fHistRCPt[i] = 0;
104 fHistRCPtExLJ[i] = 0;
105 fHistRCPtRand[i] = 0;
59f16b27 106 fHistRhoVSRCPt[i] = 0;
a487deae 107 fHistDeltaPtRC[i] = 0;
108 fHistDeltaPtRCExLJ[i] = 0;
109 fHistDeltaPtRCRand[i] = 0;
5ce8ae64 110 fHistEmbRejectedJetsPhiEta[i] = 0;
111 fHistEmbRejectedJetsPtArea[i] = 0;
6a20534a 112 fHistEmbNotFoundPt[i] = 0;
a487deae 113 fHistEmbNotFoundPhiEta[i] = 0;
1f9c287f 114 fHistEmbJetsPtArea[i] = 0;
115 fHistEmbJetsCorrPtArea[i] = 0;
6a20534a 116 fHistEmbPartPtvsJetPt[i] = 0;
ca5c29fa 117 fHistEmbPartPtvsJetCorrPt[i] = 0;
2103dc6a 118 fHistJetPtvsJetCorrPt[i] = 0;
6a20534a 119 fHistDistLeadPart2JetAxis[i] = 0;
1f9c287f 120 fHistEmbBkgArea[i] = 0;
59f16b27 121 fHistRhoVSEmbBkg[i] = 0;
1f9c287f 122 fHistDeltaPtEmbArea[i] = 0;
a487deae 123 }
124
125 SetMakeGeneralHistograms(kTRUE);
126}
127
128//________________________________________________________________________
129AliAnalysisTaskDeltaPt::~AliAnalysisTaskDeltaPt()
130{
131 // Destructor.
132}
133
134//________________________________________________________________________
135void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
136{
137 // Create user output.
138
139 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
140
6a20534a 141 if (!fTracksName.IsNull() || !fCaloName.IsNull()) {
5ce8ae64 142 fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
6a20534a 143 fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
144 fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
145 fOutput->Add(fHistRCPhiEta);
146
147 if (!fJetsName.IsNull()) {
148 fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
149 fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
150 fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
151 fOutput->Add(fHistRCPtExLJVSDPhiLJ);
152 }
153 }
154
155 if (!fEmbJetsName.IsNull()) {
5ce8ae64 156 fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
6a20534a 157 fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
158 fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
159 fOutput->Add(fHistEmbJetsPhiEta);
a487deae 160
5ce8ae64 161 fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
6a20534a 162 fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
163 fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
164 fOutput->Add(fHistLeadPartPhiEta);
a487deae 165 }
166
167 TString histname;
168
ca5c29fa 169 const Int_t nbinsZ = 12;
170 Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
171
172 Float_t *binsPt = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
173 Float_t *binsCorrPt = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
174 Float_t *binsArea = GenerateFixedBinArray(40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
175
6a20534a 176 for (Int_t i = 0; i < fNcentBins; i++) {
177 if (!fTracksName.IsNull() || !fCaloName.IsNull()) {
178 histname = "fHistRCPt_";
179 histname += i;
180 fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
181 fHistRCPt[i]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
182 fHistRCPt[i]->GetYaxis()->SetTitle("counts");
183 fOutput->Add(fHistRCPt[i]);
184
185 histname = "fHistRhoVSRCPt_";
186 histname += i;
187 fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
188 fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
189 fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
190 fOutput->Add(fHistRhoVSRCPt[i]);
191
192 histname = "fHistDeltaPtRC_";
193 histname += i;
194 fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
195 fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
196 fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
197 fOutput->Add(fHistDeltaPtRC[i]);
198
199 if (!fJetsName.IsNull()) {
200 histname = "fHistRCPtExLJ_";
201 histname += i;
202 fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
203 fHistRCPtExLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
204 fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
205 fOutput->Add(fHistRCPtExLJ[i]);
206
207 histname = "fHistDeltaPtRCExLJ_";
208 histname += i;
209 fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
210 fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
211 fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
212 fOutput->Add(fHistDeltaPtRCExLJ[i]);
213 }
214 }
215
216 if (!fRandTracksName.IsNull() || !fRandCaloName.IsNull()) {
217 histname = "fHistRCPtRand_";
218 histname += i;
219 fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
220 fHistRCPtRand[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
221 fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
222 fOutput->Add(fHistRCPtRand[i]);
223
224 histname = "fHistDeltaPtRCRand_";
225 histname += i;
226 fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
227 fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
228 fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
229 fOutput->Add(fHistDeltaPtRCRand[i]);
230 }
231
a487deae 232 if (!fEmbJetsName.IsNull()) {
1f9c287f 233 histname = "fHistEmbJetsPtArea_";
a487deae 234 histname += i;
ca5c29fa 235 fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 40, binsArea, fNbins, binsPt, nbinsZ, binsZ);
1f9c287f 236 fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
6a20534a 237 fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
1f9c287f 238 fOutput->Add(fHistEmbJetsPtArea[i]);
a487deae 239
1f9c287f 240 histname = "fHistEmbJetsCorrPtArea_";
a487deae 241 histname += i;
ca5c29fa 242 fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 40, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
1f9c287f 243 fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
6a20534a 244 fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
1f9c287f 245 fOutput->Add(fHistEmbJetsCorrPtArea[i]);
a487deae 246
6a20534a 247 histname = "fHistEmbPartPtvsJetPt_";
248 histname += i;
249 fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
250 fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
251 fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
252 fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
253 fOutput->Add(fHistEmbPartPtvsJetPt[i]);
254
ca5c29fa 255 histname = "fHistEmbPartPtvsJetCorrPt_";
a487deae 256 histname += i;
ca5c29fa 257 fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
258 fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
259 fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
260 fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
261 fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
a487deae 262
2103dc6a 263 histname = "fHistJetPtvsJetCorrPt_";
264 histname += i;
265 fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
266 fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
267 fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
268 fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
269 fOutput->Add(fHistJetPtvsJetCorrPt[i]);
270
6a20534a 271 histname = "fHistDistLeadPart2JetAxis_";
a487deae 272 histname += i;
6a20534a 273 fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
274 fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
275 fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
276 fOutput->Add(fHistDistLeadPart2JetAxis[i]);
a487deae 277
278 histname = "fHistEmbNotFoundPhiEta_";
279 histname += i;
5ce8ae64 280 fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
a487deae 281 fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
282 fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
283 fOutput->Add(fHistEmbNotFoundPhiEta[i]);
59f16b27 284
6a20534a 285 histname = "fHistEmbNotFoundPt_";
286 histname += i;
287 fHistEmbNotFoundPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
288 fHistEmbNotFoundPt[i]->GetXaxis()->SetTitle("#it{p}_{T,const}^{emb} (GeV/#it{c})");
289 fHistEmbNotFoundPt[i]->GetYaxis()->SetTitle("counts");
290 fOutput->Add(fHistEmbNotFoundPt[i]);
291
5ce8ae64 292 histname = "fHistEmbRejectedJetsPhiEta_";
293 histname += i;
294 fHistEmbRejectedJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
295 fHistEmbRejectedJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
296 fHistEmbRejectedJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
297 fOutput->Add(fHistEmbRejectedJetsPhiEta[i]);
298
299 histname = "fHistEmbRejectedJetsPtArea_";
300 histname += i;
301 fHistEmbRejectedJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
302 fHistEmbRejectedJetsPtArea[i]->GetXaxis()->SetTitle("area");
303 fHistEmbRejectedJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
304 fOutput->Add(fHistEmbRejectedJetsPtArea[i]);
305
1f9c287f 306 histname = "fHistEmbBkgArea_";
307 histname += i;
308 fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
309 fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
6a20534a 310 fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
1f9c287f 311 fOutput->Add(fHistEmbBkgArea[i]);
312
59f16b27 313 histname = "fHistRhoVSEmbBkg_";
314 histname += i;
1f9c287f 315 fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
316 fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
6a20534a 317 fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
59f16b27 318 fOutput->Add(fHistRhoVSEmbBkg[i]);
a487deae 319
1f9c287f 320 histname = "fHistDeltaPtEmbArea_";
a487deae 321 histname += i;
1f9c287f 322 fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins * 2, -fMaxBinPt, fMaxBinPt);
323 fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
324 fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
325 fOutput->Add(fHistDeltaPtEmbArea[i]);
a487deae 326 }
327 }
328
cac0e10b 329 delete[] binsPt;
330 delete[] binsCorrPt;
331 delete[] binsArea;
332
a487deae 333 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
334}
335
336//________________________________________________________________________
337Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
338{
339 // Fill histograms.
340
a487deae 341 // ************
342 // Random cones
343 // _________________________________
344
345 const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
6a20534a 346 Float_t RCpt = 0;
347 Float_t RCeta = 0;
348 Float_t RCphi = 0;
a487deae 349
6a20534a 350 if (fTracks || fCaloClusters) {
351
352 for (Int_t i = 0; i < fRCperEvent; i++) {
353 // Simple random cones
354 RCpt = 0;
355 RCeta = 0;
356 RCphi = 0;
357 GetRandomCone(RCpt, RCeta, RCphi, 0);
358 if (RCpt > 0) {
359 fHistRCPhiEta->Fill(RCeta, RCphi);
360 fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
361
362 fHistRCPt[fCentBin]->Fill(RCpt);
363 fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
364 }
365
366 if (fJets) {
367
368 // Random cones far from leading jet
ca5c29fa 369 static Int_t sortedJets[9999] = {-1};
370 GetSortedArray(sortedJets, fJets);
6a20534a 371
372 AliEmcalJet* jet = 0;
373
ca5c29fa 374 if (sortedJets[0] >= 0)
6a20534a 375 jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
376
377 RCpt = 0;
378 RCeta = 0;
379 RCphi = 0;
380 GetRandomCone(RCpt, RCeta, RCphi, jet);
381 if (RCpt > 0) {
382 if (jet) {
383 Float_t dphi = RCphi - jet->Phi();
384 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
385 if (dphi < -1.6) dphi += TMath::Pi() * 2;
386 fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
387 }
388 fHistRCPtExLJ[fCentBin]->Fill(RCpt);
389 fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
390 }
a487deae 391 }
a487deae 392 }
6a20534a 393 }
394
395 // Random cones with randomized particles
396 if (fRandTracks || fRandCaloClusters) {
a487deae 397 RCpt = 0;
398 RCeta = 0;
399 RCphi = 0;
400 GetRandomCone(RCpt, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
401 if (RCpt > 0) {
1f9c287f 402 fHistRCPtRand[fCentBin]->Fill(RCpt);
a487deae 403 fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
404 }
405 }
6a20534a 406
a487deae 407 // ************
408 // Embedding
409 // _________________________________
410
6a20534a 411 if (fEmbJets) {
412
413 AliEmcalJet *embJet = NextEmbeddedJet(0);
414
415 Int_t countEmbJets = 0;
416
417 while (embJet != 0) {
5be3857d 418 AliDebug(2,Form("Elaborating embedded jet n. %d", countEmbJets));
419 countEmbJets++;
420
5ce8ae64 421 if (!AcceptJet(embJet)) {
5be3857d 422 AliDebug(2,"Embedded jet not accepted, skipping...");
5ce8ae64 423 fHistEmbRejectedJetsPhiEta[fCentBin]->Fill(embJet->Eta(), embJet->Phi());
424 fHistEmbRejectedJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt());
425
426 embJet = NextEmbeddedJet();
427 continue;
428 }
6a20534a 429
430 Double_t maxClusterPt = 0;
431 Double_t maxClusterEta = 0;
432 Double_t maxClusterPhi = 0;
f660c2d6 433
6a20534a 434 Double_t maxTrackPt = 0;
435 Double_t maxTrackEta = 0;
436 Double_t maxTrackPhi = 0;
437
438 Double_t maxPartPt = 0;
439 Double_t maxPartEta = 0;
440 Double_t maxPartPhi = 0;
441
ca5c29fa 442 if (fLeadingHadronType == 1 || fLeadingHadronType == 2) {
443 AliVCluster *cluster = embJet->GetLeadingCluster(fEmbCaloClusters);
444 if (cluster) {
445 TLorentzVector nPart;
446 cluster->GetMomentum(nPart, fVertex);
447
448 maxClusterEta = nPart.Eta();
449 maxClusterPhi = nPart.Phi();
450 maxClusterPt = nPart.Pt();
451 }
6a20534a 452 }
453
ca5c29fa 454 if (fLeadingHadronType == 0 || fLeadingHadronType == 2) {
455 AliVParticle *track = embJet->GetLeadingTrack(fEmbTracks);
456 if (track) {
457 maxTrackEta = track->Eta();
458 maxTrackPhi = track->Phi();
459 maxTrackPt = track->Pt();
460 }
6a20534a 461 }
462
463 if (maxTrackPt > maxClusterPt) {
464 maxPartPt = maxTrackPt;
465 maxPartEta = maxTrackEta;
466 maxPartPhi = maxTrackPhi;
467 }
468 else {
469 maxPartPt = maxClusterPt;
470 maxPartEta = maxClusterEta;
471 maxPartPhi = maxClusterPhi;
472 }
473
474 Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - maxPartEta) * (embJet->Eta() - maxPartEta) + (embJet->Phi() - maxPartPhi) * (embJet->Phi() - maxPartPhi));
475
476 fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
ca5c29fa 477 fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
6a20534a 478 fHistLeadPartPhiEta->Fill(maxPartEta, maxPartPhi);
6a20534a 479 fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
480
ca5c29fa 481 fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), maxPartPt);
482 fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), maxPartPt);
6a20534a 483 fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
2103dc6a 484 fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
6a20534a 485
486 fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
487 fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
488 fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
f660c2d6 489
6a20534a 490 embJet = NextEmbeddedJet();
a487deae 491 }
f660c2d6 492
6a20534a 493 if (countEmbJets==0) {
ca5c29fa 494 AliDebug(1,"No embedded jets found!");
6a20534a 495 if (fEmbTracks) {
496 DoEmbTrackLoop();
497 for (Int_t i = 0; i < fEmbeddedTrackNIds; i++) {
ca5c29fa 498 AliDebug(2,Form("Embedded track %d found!",i));
6a20534a 499 AliVParticle *track2 = static_cast<AliVParticle*>(fEmbTracks->At(fEmbeddedTrackIds[i]));
500 if (!track2) continue;
501 fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
502 fHistEmbNotFoundPt[fCentBin]->Fill(track2->Pt());
503 }
f660c2d6 504 }
6a20534a 505
506 if (fEmbCaloClusters) {
507 DoEmbClusterLoop();
508 for (Int_t i = 0; i < fEmbeddedClusterNIds; i++) {
ca5c29fa 509 AliDebug(2,Form("Embedded cluster %d found!",i));
6a20534a 510 AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterIds[i]));
511 TLorentzVector nPart;
512 cluster2->GetMomentum(nPart, fVertex);
513 fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
514 fHistEmbNotFoundPt[fCentBin]->Fill(nPart.Pt());
515 }
f660c2d6 516 }
a487deae 517 }
518 }
519
520 return kTRUE;
521}
522
523//________________________________________________________________________
524void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
525{
526 // Do track loop.
527
f660c2d6 528 fEmbeddedTrackNIds = 0;
529
a487deae 530 if (!fEmbTracks)
531 return;
532
a487deae 533 Int_t ntracks = fEmbTracks->GetEntriesFast();
534
535 for (Int_t i = 0; i < ntracks; i++) {
536
537 AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track
538
539 if (!track) {
540 AliError(Form("Could not retrieve track %d",i));
541 continue;
542 }
543
544 AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
545
f660c2d6 546 if (vtrack && !AcceptTrack(vtrack))
a487deae 547 continue;
548
7f76e479 549 if (track->GetLabel() > 0) {
f660c2d6 550 fEmbeddedTrackIds[fEmbeddedTrackNIds] = i;
551 fEmbeddedTrackNIds++;
a487deae 552 }
553 }
554}
555
556//________________________________________________________________________
557void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
558{
559 // Do cluster loop.
560
f660c2d6 561 fEmbeddedClusterNIds = 0;
562
a487deae 563 if (!fEmbCaloClusters)
564 return;
565
a487deae 566 Int_t nclusters = fEmbCaloClusters->GetEntriesFast();
567
568 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
569 AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
570 if (!cluster) {
571 AliError(Form("Could not receive cluster %d", iClusters));
572 continue;
573 }
574
f660c2d6 575 if (!AcceptCluster(cluster))
a487deae 576 continue;
577
7f76e479 578 if (cluster->GetLabel() > 0) {
f660c2d6 579 fEmbeddedClusterIds[fEmbeddedClusterNIds] = iClusters;
580 fEmbeddedClusterNIds++;
a487deae 581 }
582 }
583}
584
585//________________________________________________________________________
f660c2d6 586AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Int_t i)
a487deae 587{
588 // Do the embedded jet loop.
589
f660c2d6 590 static Int_t iJet = 0;
591
592 if (i >= 0)
593 iJet = i;
7f76e479 594 else
595 iJet++;
f660c2d6 596
a487deae 597 if (!fEmbJets)
f660c2d6 598 return 0;
a487deae 599
600 TLorentzVector maxClusVect;
601
f660c2d6 602 const Int_t nembjets = fEmbJets->GetEntriesFast();
a487deae 603
f660c2d6 604 for (; iJet < nembjets; iJet++) {
a487deae 605
f660c2d6 606 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(iJet));
a487deae 607
608 if (!jet) {
f660c2d6 609 AliError(Form("Could not receive jet %d", iJet));
a487deae 610 continue;
611 }
ca5c29fa 612
613 if (jet->MCPt() < fMCJetPtThreshold)
614 continue;
5ce8ae64 615
ca5c29fa 616 return jet;
a487deae 617 }
f660c2d6 618
619 return 0;
a487deae 620}
621
622//________________________________________________________________________
623void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
624 AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
625{
626 // Get rigid cone.
627
628 if (!tracks)
629 tracks = fTracks;
630
631 if (!clusters)
632 clusters = fCaloClusters;
633
634 eta = 0;
635 phi = 0;
636 pt = 0;
637
638 if (!tracks && !clusters)
639 return;
640
641 Float_t LJeta = 999;
642 Float_t LJphi = 999;
643
644 if (jet) {
645 LJeta = jet->Eta();
646 LJphi = jet->Phi();
647 }
648
649 Float_t maxEta = fJetMaxEta;
650 Float_t minEta = fJetMinEta;
651 Float_t maxPhi = fJetMaxPhi;
652 Float_t minPhi = fJetMinPhi;
653
654 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
655 if (minPhi < 0) minPhi = 0;
656
657 Float_t dLJ = 0;
658 Int_t repeats = 0;
659 do {
660 eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
661 phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
662 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
663 repeats++;
664 } while (dLJ < fMinRC2LJ && repeats < 999);
665
666 if (repeats == 999) {
667 AliWarning(Form("%s: Could not get random cone!", GetName()));
668 return;
669 }
670
671 if (clusters) {
672 Int_t nclusters = clusters->GetEntriesFast();
673 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
674 AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
675 if (!cluster) {
676 AliError(Form("Could not receive cluster %d", iClusters));
677 continue;
678 }
679
f660c2d6 680 if (!AcceptCluster(cluster))
a487deae 681 continue;
682
683 TLorentzVector nPart;
684 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
685
686 Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
687
688 if (d <= fJetRadius)
689 pt += nPart.Pt();
690 }
691 }
692
693 if (tracks) {
694 Int_t ntracks = tracks->GetEntriesFast();
695 for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
696 AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
697 if(!track) {
698 AliError(Form("Could not retrieve track %d",iTracks));
699 continue;
700 }
7f76e479 701
f660c2d6 702 if (!AcceptTrack(track))
a487deae 703 continue;
704
705 Float_t tracketa = track->Eta();
706 Float_t trackphi = track->Phi();
707
708 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
709 trackphi += 2 * TMath::Pi();
710 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
711 trackphi -= 2 * TMath::Pi();
712
713 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
714 if (d <= fJetRadius)
715 pt += track->Pt();
716 }
717 }
718}
719
720//________________________________________________________________________
721void AliAnalysisTaskDeltaPt::ExecOnce()
722{
723 // Initialize the analysis.
724
725 if (!fEmbJetsName.IsNull() && !fEmbJets) {
726 fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
727 if (!fEmbJets) {
728 AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
729 return;
730 }
731 else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
732 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
733 fEmbJets = 0;
734 return;
735 }
736 }
737
738 if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
739 fEmbCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
740 if (!fEmbCaloClusters) {
741 AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
742 return;
743 }
744 else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
745 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
746 fEmbCaloClusters = 0;
747 return;
748 }
749 }
750
751 if (!fEmbTracksName.IsNull() && !fEmbTracks) {
752 fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
753 if (!fEmbTracks) {
754 AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
755 return;
756 }
757 else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
758 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
759 fEmbTracks = 0;
760 return;
761 }
762 }
763
764 if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
765 fRandCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
766 if (!fRandCaloClusters) {
767 AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
768 return;
769 }
770 else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
771 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
772 fRandCaloClusters = 0;
773 return;
774 }
775 }
776
777 if (!fRandTracksName.IsNull() && !fRandTracks) {
778 fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
779 if (!fRandTracks) {
780 AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
781 return;
782 }
783 else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
784 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
785 fRandTracks = 0;
786 return;
787 }
788 }
789
790 AliAnalysisTaskEmcalJet::ExecOnce();
791
792 if (fRCperEvent < 0) {
793 Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
794 Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
795 fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
796 if (fRCperEvent == 0)
797 fRCperEvent = 1;
798 }
799
800 if (fMinRC2LJ < 0)
801 fMinRC2LJ = fJetRadius * 1.5;
802
803 const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
804 if (fMinRC2LJ > maxDist) {
805 AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
806 "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
807 fMinRC2LJ = maxDist;
808 }
809}
810
811//________________________________________________________________________
812void AliAnalysisTaskDeltaPt::Terminate(Option_t *)
813{
814 // Called once at the end of the analysis.
815}