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