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