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