]> 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
updates 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>
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;
59 fHistEmbJetsPt[i] = 0;
60 fHistEmbJetsCorrPt[i] = 0;
61 fHistEmbJetsArea[i] = 0;
62 fHistEmbPartPt[i] = 0;
63 fHistDistEmbPartJetAxis[i] = 0;
59f16b27 64 fHistRhoVSEmbBkg[i] = 0;
a487deae 65 fHistDeltaPtEmb[i] = 0;
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;
105 fHistEmbJetsPt[i] = 0;
106 fHistEmbJetsCorrPt[i] = 0;
107 fHistEmbJetsArea[i] = 0;
108 fHistEmbPartPt[i] = 0;
109 fHistDistEmbPartJetAxis[i] = 0;
59f16b27 110 fHistRhoVSEmbBkg[i] = 0;
a487deae 111 fHistDeltaPtEmb[i] = 0;
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);
136 fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
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);
158 fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
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);
165 fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
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);
172 fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
173 fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
174 fOutput->Add(fHistRCPtRand[i]);
175
59f16b27 176 histname = "fHistRhoVSRCPt_";
177 histname += i;
178 fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
179 fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho [GeV/c]");
180 fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
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);
186 fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
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);
193 fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
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);
200 fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
201 fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
202 fOutput->Add(fHistDeltaPtRCRand[i]);
203
204 if (!fEmbJetsName.IsNull()) {
205 histname = "fHistEmbJetsPt_";
206 histname += i;
207 fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
208 fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
209 fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
210 fOutput->Add(fHistEmbJetsPt[i]);
211
212 histname = "fHistEmbJetsCorrPt_";
213 histname += i;
214 fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
215 fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
216 fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
217 fOutput->Add(fHistEmbJetsCorrPt[i]);
218
219 histname = "fHistEmbJetsArea_";
220 histname += i;
221 fHistEmbJetsArea[i] = new TH1F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
222 fHistEmbJetsArea[i]->GetXaxis()->SetTitle("area");
223 fHistEmbJetsArea[i]->GetYaxis()->SetTitle("counts");
224 fOutput->Add(fHistEmbJetsArea[i]);
225
226 histname = "fHistEmbPartPt_";
227 histname += i;
228 fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
229 fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
230 fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
231 fOutput->Add(fHistEmbPartPt[i]);
232
233 histname = "fHistDistEmbPartJetAxis_";
234 histname += i;
235 fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
236 fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
237 fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
238 fOutput->Add(fHistDistEmbPartJetAxis[i]);
239
240 histname = "fHistEmbNotFoundPhiEta_";
241 histname += i;
242 fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
243 fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
244 fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
245 fOutput->Add(fHistEmbNotFoundPhiEta[i]);
59f16b27 246
247 histname = "fHistRhoVSEmbBkg_";
248 histname += i;
249 fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
250 fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho [GeV/c]");
251 fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
252 fOutput->Add(fHistRhoVSEmbBkg[i]);
a487deae 253
254 histname = "fHistDeltaPtEmb_";
255 histname += i;
256 fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
257 fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
258 fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
259 fOutput->Add(fHistDeltaPtEmb[i]);
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
284 // Simple random cones
285 for (Int_t i = 0; i < fRCperEvent; i++) {
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);
59f16b27 292 fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal, RCpt);
a487deae 293
294 fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
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 }
310 fHistRCPtExLJ[fCentBin]->Fill(RCpt / rcArea);
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) {
320 fHistRCPtRand[fCentBin]->Fill(RCpt / rcArea);
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
370 fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
371 fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - fRhoVal * embJet->Area());
372 fHistEmbJetsArea[fCentBin]->Fill(embJet->Area());
373 fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
374
375 fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRhoVal - probePt);
59f16b27 376 fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal, embJet->Pt() - probePt);
a487deae 377 }
378 else {
379 if (fEmbTracks)
380 DoEmbTrackLoop();
381 if (fEmbCaloClusters)
382 DoEmbClusterLoop();
383 if (fEmbeddedTrackId >= 0) {
384 AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
385 fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
386 }
387 else if (fEmbeddedClusterId >= 0) {
388 AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
389 TLorentzVector nPart;
390 cluster2->GetMomentum(nPart, fVertex);
391 fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
392 }
393 else {
394 AliWarning(Form("%s - Embedded particle not found!", GetName()));
395 }
396 }
397
398 return kTRUE;
399}
400
401//________________________________________________________________________
402void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
403{
404 // Do track loop.
405
406 if (!fEmbTracks)
407 return;
408
409 fEmbeddedTrackId = -1;
410
411 Int_t ntracks = fEmbTracks->GetEntriesFast();
412
413 for (Int_t i = 0; i < ntracks; i++) {
414
415 AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track
416
417 if (!track) {
418 AliError(Form("Could not retrieve track %d",i));
419 continue;
420 }
421
422 AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
423
424 if (vtrack && !AcceptTrack(vtrack, kTRUE))
425 continue;
426
427 if (track->GetLabel() == 100) {
428 fEmbeddedTrackId = i;
429 continue;
430 }
431 }
432}
433
434//________________________________________________________________________
435void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
436{
437 // Do cluster loop.
438
439 if (!fEmbCaloClusters)
440 return;
441
442 fEmbeddedClusterId = -1;
443
444 Int_t nclusters = fEmbCaloClusters->GetEntriesFast();
445
446 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
447 AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
448 if (!cluster) {
449 AliError(Form("Could not receive cluster %d", iClusters));
450 continue;
451 }
452
453 if (!AcceptCluster(cluster, kTRUE))
454 continue;
455
456 if (cluster->Chi2() == 100) {
457 fEmbeddedClusterId = iClusters;
458 continue;
459 }
460 }
461}
462
463//________________________________________________________________________
464void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
465{
466 // Do the embedded jet loop.
467
468 if (!fEmbJets)
469 return;
470
471 TLorentzVector maxClusVect;
472
473 Int_t nembjets = fEmbJets->GetEntriesFast();
474
475 for (Int_t ij = 0; ij < nembjets; ij++) {
476
477 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
478
479 if (!jet) {
480 AliError(Form("Could not receive jet %d", ij));
481 continue;
482 }
483
484 if (!AcceptJet(jet))
485 continue;
486
487 if (!jet->IsMC())
488 continue;
489
490 AliVParticle *maxTrack = 0;
491
492 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
493 AliVParticle *track = jet->TrackAt(it, fEmbTracks);
494
495 if (!track)
496 continue;
497
498 if (track->GetLabel() != 100)
499 continue;
500
501 if (!maxTrack || track->Pt() > maxTrack->Pt())
502 maxTrack = track;
503 }
504
505 AliVCluster *maxClus = 0;
506
507 for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
508 AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
509
510 if (!cluster)
511 continue;
512
513 if (cluster->Chi2() != 100)
514 continue;
515
516 TLorentzVector nPart;
517 cluster->GetMomentum(nPart, fVertex);
518
519 if (!maxClus || nPart.Et() > maxClusVect.Et()) {
520 new (&maxClusVect) TLorentzVector(nPart);
521 maxClus = cluster;
522 }
523 }
524
525 if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
526 embPart = maxClus;
527 embJet = jet;
528 }
529 else if (maxTrack) {
530 embPart = maxTrack;
531 embJet = jet;
532 }
533
534 return; // MC jet found, exit
535 }
536}
537
538//________________________________________________________________________
539void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
540 AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
541{
542 // Get rigid cone.
543
544 if (!tracks)
545 tracks = fTracks;
546
547 if (!clusters)
548 clusters = fCaloClusters;
549
550 eta = 0;
551 phi = 0;
552 pt = 0;
553
554 if (!tracks && !clusters)
555 return;
556
557 Float_t LJeta = 999;
558 Float_t LJphi = 999;
559
560 if (jet) {
561 LJeta = jet->Eta();
562 LJphi = jet->Phi();
563 }
564
565 Float_t maxEta = fJetMaxEta;
566 Float_t minEta = fJetMinEta;
567 Float_t maxPhi = fJetMaxPhi;
568 Float_t minPhi = fJetMinPhi;
569
570 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
571 if (minPhi < 0) minPhi = 0;
572
573 Float_t dLJ = 0;
574 Int_t repeats = 0;
575 do {
576 eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
577 phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
578 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
579 repeats++;
580 } while (dLJ < fMinRC2LJ && repeats < 999);
581
582 if (repeats == 999) {
583 AliWarning(Form("%s: Could not get random cone!", GetName()));
584 return;
585 }
586
587 if (clusters) {
588 Int_t nclusters = clusters->GetEntriesFast();
589 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
590 AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
591 if (!cluster) {
592 AliError(Form("Could not receive cluster %d", iClusters));
593 continue;
594 }
595
596 if (!AcceptCluster(cluster, fMCAna))
597 continue;
598
599 TLorentzVector nPart;
600 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
601
602 Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
603
604 if (d <= fJetRadius)
605 pt += nPart.Pt();
606 }
607 }
608
609 if (tracks) {
610 Int_t ntracks = tracks->GetEntriesFast();
611 for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
612 AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
613 if(!track) {
614 AliError(Form("Could not retrieve track %d",iTracks));
615 continue;
616 }
617
618 if (!AcceptTrack(track, fMCAna))
619 continue;
620
621 Float_t tracketa = track->Eta();
622 Float_t trackphi = track->Phi();
623
624 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
625 trackphi += 2 * TMath::Pi();
626 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
627 trackphi -= 2 * TMath::Pi();
628
629 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
630 if (d <= fJetRadius)
631 pt += track->Pt();
632 }
633 }
634}
635
636//________________________________________________________________________
637void AliAnalysisTaskDeltaPt::ExecOnce()
638{
639 // Initialize the analysis.
640
641 if (!fEmbJetsName.IsNull() && !fEmbJets) {
642 fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
643 if (!fEmbJets) {
644 AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
645 return;
646 }
647 else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
648 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
649 fEmbJets = 0;
650 return;
651 }
652 }
653
654 if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
655 fEmbCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
656 if (!fEmbCaloClusters) {
657 AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
658 return;
659 }
660 else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
661 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
662 fEmbCaloClusters = 0;
663 return;
664 }
665 }
666
667 if (!fEmbTracksName.IsNull() && !fEmbTracks) {
668 fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
669 if (!fEmbTracks) {
670 AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
671 return;
672 }
673 else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
674 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
675 fEmbTracks = 0;
676 return;
677 }
678 }
679
680 if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
681 fRandCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
682 if (!fRandCaloClusters) {
683 AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
684 return;
685 }
686 else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
687 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
688 fRandCaloClusters = 0;
689 return;
690 }
691 }
692
693 if (!fRandTracksName.IsNull() && !fRandTracks) {
694 fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
695 if (!fRandTracks) {
696 AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
697 return;
698 }
699 else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
700 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
701 fRandTracks = 0;
702 return;
703 }
704 }
705
706 AliAnalysisTaskEmcalJet::ExecOnce();
707
708 if (fRCperEvent < 0) {
709 Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
710 Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
711 fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
712 if (fRCperEvent == 0)
713 fRCperEvent = 1;
714 }
715
716 if (fMinRC2LJ < 0)
717 fMinRC2LJ = fJetRadius * 1.5;
718
719 const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
720 if (fMinRC2LJ > maxDist) {
721 AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
722 "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
723 fMinRC2LJ = maxDist;
724 }
725}
726
727//________________________________________________________________________
728void AliAnalysisTaskDeltaPt::Terminate(Option_t *)
729{
730 // Called once at the end of the analysis.
731}