AliPhysics  b555aef (b555aef)
AliAnalysisTaskDeltaPt.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Jet deltaPt task.
4 //
5 // Author: S.Aiola
6 
7 #include <TClonesArray.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <TList.h>
12 #include <TLorentzVector.h>
13 #include <TRandom3.h>
14 
15 #include "AliVCluster.h"
16 #include "AliVParticle.h"
17 #include "AliVTrack.h"
18 #include "AliEmcalJet.h"
19 #include "AliRhoParameter.h"
20 #include "AliLog.h"
21 #include "AliJetContainer.h"
22 #include "AliParticleContainer.h"
23 #include "AliClusterContainer.h"
24 
25 #include "AliAnalysisTaskDeltaPt.h"
26 
27 ClassImp(AliAnalysisTaskDeltaPt)
28 
29 //________________________________________________________________________
32  fMCJetPtThreshold(1),
33  fMinRC2LJ(-1),
34  fRCperEvent(-1),
35  fConeRadius(0.2),
36  fConeMinEta(-0.9),
37  fConeMaxEta(0.9),
38  fConeMinPhi(0),
39  fConeMaxPhi(TMath::Pi()*2),
40  fJetsCont(0),
41  fTracksCont(0),
42  fCaloClustersCont(0),
43  fEmbJetsCont(0),
44  fEmbTracksCont(0),
45  fEmbCaloClustersCont(0),
46  fRandTracksCont(0),
47  fRandCaloClustersCont(0),
48  fHistRhovsCent(0),
49  fHistRCPhiEta(0),
50  fHistRCPt(0),
51  fHistRCPtExLJ(0),
52  fHistRCPtExPartialLJ(0),
53  fHistRCPtRand(0),
54  fHistRhoVSRCPt(0),
55  fHistDeltaPtRCvsEP(0),
56  fHistDeltaPtRCExLJ(0),
57  fHistDeltaPtRCExPartialLJ(0),
58  fHistDeltaPtRCRand(0),
59  fHistEmbJetsPtArea(0),
60  fHistEmbJetsCorrPtArea(0),
61  fHistEmbPartPtvsJetPt(0),
62  fHistEmbPartPtvsJetCorrPt(0),
63  fHistJetPtvsJetCorrPt(0),
64  fHistDistLeadPart2JetAxis(0),
65  fHistEmbBkgArea(0),
66  fHistRhoVSEmbBkg(0),
67  fHistDeltaPtEmbArea(0),
68  fHistDeltaPtEmbvsEP(0),
69  fHistRCPtExLJVSDPhiLJ(0),
70  fHistRCPtExPartialLJVSDPhiLJ(0),
71  fHistEmbJetsPhiEta(0),
72  fHistLeadPartPhiEta(0)
73 {
74  // Default constructor.
75 
76  fHistRCPt = 0;
77  fHistRCPtExLJ = 0;
78  fHistRCPtExPartialLJ = 0;
79  fHistRCPtRand = 0;
80  fHistRhoVSRCPt = 0;
81  fHistDeltaPtRCvsEP = 0;
82  fHistDeltaPtRCExLJ = 0;
83  fHistDeltaPtRCExPartialLJ = 0;
84  fHistDeltaPtRCRand = 0;
85  fHistEmbJetsPtArea = 0;
86  fHistEmbJetsCorrPtArea = 0;
87  fHistEmbPartPtvsJetPt = 0;
88  fHistEmbPartPtvsJetCorrPt = 0;
89  fHistJetPtvsJetCorrPt = 0;
90  fHistDistLeadPart2JetAxis = 0;
91  fHistEmbBkgArea = 0;
92  fHistRhoVSEmbBkg = 0;
93  fHistDeltaPtEmbArea = 0;
94  fHistDeltaPtEmbvsEP = 0;
95 
96  SetMakeGeneralHistograms(kTRUE);
97 }
98 
99 //________________________________________________________________________
101  AliAnalysisTaskEmcalJet(name, kTRUE),
102  fMCJetPtThreshold(1),
103  fMinRC2LJ(-1),
104  fRCperEvent(-1),
105  fConeRadius(0.2),
106  fConeMinEta(-0.9),
107  fConeMaxEta(0.9),
108  fConeMinPhi(0),
109  fConeMaxPhi(TMath::Pi()*2),
110  fJetsCont(0),
111  fTracksCont(0),
112  fCaloClustersCont(0),
113  fEmbJetsCont(0),
114  fEmbTracksCont(0),
115  fEmbCaloClustersCont(0),
116  fRandTracksCont(0),
117  fRandCaloClustersCont(0),
118  fHistRhovsCent(0),
119  fHistRCPhiEta(0),
120  fHistRCPt(0),
121  fHistRCPtExLJ(0),
122  fHistRCPtExPartialLJ(0),
123  fHistRCPtRand(0),
124  fHistRhoVSRCPt(0),
125  fHistDeltaPtRCvsEP(0),
126  fHistDeltaPtRCExLJ(0),
127  fHistDeltaPtRCExPartialLJ(0),
128  fHistDeltaPtRCRand(0),
129  fHistEmbJetsPtArea(0),
130  fHistEmbJetsCorrPtArea(0),
131  fHistEmbPartPtvsJetPt(0),
132  fHistEmbPartPtvsJetCorrPt(0),
133  fHistJetPtvsJetCorrPt(0),
134  fHistDistLeadPart2JetAxis(0),
135  fHistEmbBkgArea(0),
136  fHistRhoVSEmbBkg(0),
137  fHistDeltaPtEmbArea(0),
138  fHistDeltaPtEmbvsEP(0),
139  fHistRCPtExLJVSDPhiLJ(0),
140  fHistRCPtExPartialLJVSDPhiLJ(0),
141  fHistEmbJetsPhiEta(0),
142  fHistLeadPartPhiEta(0)
143 {
144  // Standard constructor.
145 
146  fHistRCPt = 0;
147  fHistRCPtExLJ = 0;
149  fHistRCPtRand = 0;
150  fHistRhoVSRCPt = 0;
151  fHistDeltaPtRCvsEP = 0;
152  fHistDeltaPtRCExLJ = 0;
154  fHistDeltaPtRCRand = 0;
155  fHistEmbJetsPtArea = 0;
161  fHistEmbBkgArea = 0;
162  fHistRhoVSEmbBkg = 0;
165 
167 }
168 
169 //________________________________________________________________________
171 {
172  fHistRCPt = new TH1*[fNcentBins];
173  fHistRCPtExLJ = new TH1*[fNcentBins];
175  fHistRCPtRand = new TH1*[fNcentBins];
176  fHistRhoVSRCPt = new TH2*[fNcentBins];
191 
192  for (Int_t i = 0; i < fNcentBins; i++) {
193  fHistRCPt[i] = 0;
194  fHistRCPtExLJ[i] = 0;
195  fHistRCPtExPartialLJ[i] = 0;
196  fHistRCPtRand[i] = 0;
197  fHistRhoVSRCPt[i] = 0;
198  fHistDeltaPtRCvsEP[i] = 0;
199  fHistDeltaPtRCExLJ[i] = 0;
201  fHistDeltaPtRCRand[i] = 0;
202  fHistEmbJetsPtArea[i] = 0;
203  fHistEmbJetsCorrPtArea[i] = 0;
204  fHistEmbPartPtvsJetPt[i] = 0;
206  fHistJetPtvsJetCorrPt[i] = 0;
208  fHistEmbBkgArea[i] = 0;
209  fHistRhoVSEmbBkg[i] = 0;
210  fHistDeltaPtEmbArea[i] = 0;
211  fHistDeltaPtEmbvsEP[i] = 0;
212  }
213 }
214 
215 //________________________________________________________________________
217 {
218  // Create user output.
219 
221 
223 
224  fHistRhovsCent = new TH2F("fHistRhovsCent", "fHistRhovsCent", 101, -1, 100, fNbins, 0, fMaxBinPt*2);
225  fHistRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
226  fHistRhovsCent->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
227  fOutput->Add(fHistRhovsCent);
228 
229  fJetsCont = GetJetContainer("Jets");
230  fTracksCont = GetParticleContainer("Tracks");
231  fCaloClustersCont = GetClusterContainer("CaloClusters");
232  fEmbJetsCont = GetJetContainer("EmbJets");
233  fEmbTracksCont = GetParticleContainer("EmbTracks");
234  fEmbCaloClustersCont = GetClusterContainer("EmbCaloClusters");
235  fRandTracksCont = GetParticleContainer("RandTracks");
236  fRandCaloClustersCont = GetClusterContainer("RandCaloClusters");
237 
239  fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
240  fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
241  fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
242  fOutput->Add(fHistRCPhiEta);
243 
244  if (fJetsCont) {
245  fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
246  fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
247  fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
249 
250  fHistRCPtExPartialLJVSDPhiLJ = new TH2F("fHistRCPtExPartialLJVSDPhiLJ","fHistRCPtExPartialLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
251  fHistRCPtExPartialLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
252  fHistRCPtExPartialLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
254  }
255  }
256 
257  if (fEmbJetsCont) {
258  fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
259  fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
260  fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
262 
263  fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
264  fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
265  fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
267  }
268 
269  TString histname;
270 
271  const Int_t nbinsZ = 12;
272  Double_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
273 
276  Double_t *binsArea = GenerateFixedBinArray(50, 0, 2);
277 
278  for (Int_t i = 0; i < fNcentBins; i++) {
280  histname = "fHistRCPt_";
281  histname += i;
282  fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
283  fHistRCPt[i]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
284  fHistRCPt[i]->GetYaxis()->SetTitle("counts");
285  fOutput->Add(fHistRCPt[i]);
286 
287  histname = "fHistRhoVSRCPt_";
288  histname += i;
289  fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
290  fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
291  fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
292  fOutput->Add(fHistRhoVSRCPt[i]);
293 
294  histname = "fHistDeltaPtRCvsEP_";
295  histname += i;
296  fHistDeltaPtRCvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
297  fHistDeltaPtRCvsEP[i]->GetXaxis()->SetTitle("#phi_{RC} - #psi_{RP}");
298  fHistDeltaPtRCvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
299  fHistDeltaPtRCvsEP[i]->GetZaxis()->SetTitle("counts");
300  fOutput->Add(fHistDeltaPtRCvsEP[i]);
301 
302  if (fJetsCont) {
303  histname = "fHistRCPtExLJ_";
304  histname += i;
305  fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
306  fHistRCPtExLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
307  fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
308  fOutput->Add(fHistRCPtExLJ[i]);
309 
310  histname = "fHistDeltaPtRCExLJ_";
311  histname += i;
312  fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
313  fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
314  fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
315  fOutput->Add(fHistDeltaPtRCExLJ[i]);
316 
317  histname = "fHistRCPtExPartialLJ_";
318  histname += i;
319  fHistRCPtExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
320  fHistRCPtExPartialLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
321  fHistRCPtExPartialLJ[i]->GetYaxis()->SetTitle("counts");
322  fOutput->Add(fHistRCPtExPartialLJ[i]);
323 
324  histname = "fHistDeltaPtRCExPartialLJ_";
325  histname += i;
326  fHistDeltaPtRCExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
327  fHistDeltaPtRCExPartialLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
328  fHistDeltaPtRCExPartialLJ[i]->GetYaxis()->SetTitle("counts");
330  }
331  }
332 
334  histname = "fHistRCPtRand_";
335  histname += i;
336  fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
337  fHistRCPtRand[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
338  fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
339  fOutput->Add(fHistRCPtRand[i]);
340 
341  histname = "fHistDeltaPtRCRand_";
342  histname += i;
343  fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
344  fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
345  fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
346  fOutput->Add(fHistDeltaPtRCRand[i]);
347  }
348 
349  if (fEmbJetsCont) {
350  histname = "fHistEmbJetsPtArea_";
351  histname += i;
352  fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
353  fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
354  fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
355  fOutput->Add(fHistEmbJetsPtArea[i]);
356 
357  histname = "fHistEmbJetsCorrPtArea_";
358  histname += i;
359  fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
360  fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
361  fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
363 
364  histname = "fHistEmbPartPtvsJetPt_";
365  histname += i;
366  fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
367  fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
368  fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
369  fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
371 
372  histname = "fHistEmbPartPtvsJetCorrPt_";
373  histname += i;
374  fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
375  fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
376  fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
377  fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
378  fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
380 
381  histname = "fHistJetPtvsJetCorrPt_";
382  histname += i;
383  fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
384  fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
385  fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
386  fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
387  fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
389 
390  histname = "fHistDistLeadPart2JetAxis_";
391  histname += i;
392  fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
393  fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
394  fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
396 
397  histname = "fHistEmbBkgArea_";
398  histname += i;
399  fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
400  fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
401  fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
402  fOutput->Add(fHistEmbBkgArea[i]);
403 
404  histname = "fHistRhoVSEmbBkg_";
405  histname += i;
406  fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
407  fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
408  fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
409  fOutput->Add(fHistRhoVSEmbBkg[i]);
410 
411  histname = "fHistDeltaPtEmbArea_";
412  histname += i;
413  fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
414  50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
415  fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
416  fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
417  fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
418  fOutput->Add(fHistDeltaPtEmbArea[i]);
419 
420  histname = "fHistDeltaPtEmbvsEP_";
421  histname += i;
422  fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
423  fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
424  fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
425  fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
426  fOutput->Add(fHistDeltaPtEmbvsEP[i]);
427  }
428  }
429 
430  delete[] binsPt;
431  delete[] binsCorrPt;
432  delete[] binsArea;
433 
434  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
435 }
436 
437 //________________________________________________________________________
439 {
440  // Fill histograms.
441 
443 
444  // ************
445  // Random cones
446  // _________________________________
447 
448  const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
449  Float_t RCpt = 0;
450  Float_t RCeta = 0;
451  Float_t RCphi = 0;
452 
454 
455  for (Int_t i = 0; i < fRCperEvent; i++) {
456  // Simple random cones
457  RCpt = 0;
458  RCeta = 0;
459  RCphi = 0;
460  GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
461  if (RCpt > 0) {
462  fHistRCPhiEta->Fill(RCeta, RCphi);
463  fHistRhoVSRCPt[fCentBin]->Fill(fJetsCont->GetRhoVal() * rcArea, RCpt);
464 
465  fHistRCPt[fCentBin]->Fill(RCpt);
466 
467  Double_t ep = RCphi - fEPV0;
468  while (ep < 0) ep += TMath::Pi();
469  while (ep >= TMath::Pi()) ep -= TMath::Pi();
470 
471  fHistDeltaPtRCvsEP[fCentBin]->Fill(ep, RCpt - rcArea * fJetsCont->GetRhoVal());
472  }
473 
474  if (fJetsCont) {
475 
476  // Random cones far from leading jet
477  AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
478 
479  RCpt = 0;
480  RCeta = 0;
481  RCphi = 0;
482  GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
483  if (RCpt > 0) {
484  if (jet) {
485  Float_t dphi = RCphi - jet->Phi();
486  if (dphi > 4.8) dphi -= TMath::Pi() * 2;
487  if (dphi < -1.6) dphi += TMath::Pi() * 2;
488  fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
489  }
490  fHistRCPtExLJ[fCentBin]->Fill(RCpt);
491  fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fJetsCont->GetRhoVal());
492  }
493 
494  //partial exclusion
495  if(fBeamType == kpA) {
496 
497  RCpt = 0;
498  RCeta = 0;
499  RCphi = 0;
500  GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet, kTRUE);
501 
502  if (RCpt > 0) {
503  if (jet) {
504  Float_t dphi = RCphi - jet->Phi();
505  if (dphi > 4.8) dphi -= TMath::Pi() * 2;
506  if (dphi < -1.6) dphi += TMath::Pi() * 2;
507  fHistRCPtExPartialLJVSDPhiLJ->Fill(RCpt, dphi);
508  }
509  fHistRCPtExPartialLJ[fCentBin]->Fill(RCpt);
510  fHistDeltaPtRCExPartialLJ[fCentBin]->Fill(RCpt - rcArea * fJetsCont->GetRhoVal());
511  }
512  }
513  }
514  }
515  }
516 
517  // Random cones with randomized particles
519  RCpt = 0;
520  RCeta = 0;
521  RCphi = 0;
522  GetRandomCone(RCpt, RCeta, RCphi, fRandTracksCont, fRandCaloClustersCont, 0);
523  if (RCpt > 0) {
524  fHistRCPtRand[fCentBin]->Fill(RCpt);
525  fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fJetsCont->GetRhoVal());
526  }
527  }
528 
529  // ************
530  // Embedding
531  // _________________________________
532 
533  if (fEmbJetsCont) {
534 
535  AliEmcalJet *embJet = NextEmbeddedJet(kTRUE);
536 
537  while (embJet != 0) {
538  TLorentzVector mom;
540 
541  Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - mom.Eta()) * (embJet->Eta() - mom.Eta()) + (embJet->Phi() - mom.Phi()) * (embJet->Phi() - mom.Phi()));
542 
543  fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
544  fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fJetsCont->GetRhoVal());
545  fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi());
546  fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
547 
548  fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), mom.Pt());
549  fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fJetsCont->GetRhoVal() * embJet->Area(), mom.Pt());
550  fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
551  fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fJetsCont->GetRhoVal() * embJet->Area());
552 
553  fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
554  fHistRhoVSEmbBkg[fCentBin]->Fill(fJetsCont->GetRhoVal() * embJet->Area(), embJet->Pt() - embJet->MCPt());
555  fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fJetsCont->GetRhoVal() - embJet->MCPt());
556 
557  Double_t ep = embJet->Phi() - fEPV0;
558  while (ep < 0) ep += TMath::Pi();
559  while (ep >= TMath::Pi()) ep -= TMath::Pi();
560 
561  fHistDeltaPtEmbvsEP[fCentBin]->Fill(ep, embJet->Pt() - embJet->Area() * fJetsCont->GetRhoVal() - embJet->MCPt());
562 
563  embJet = NextEmbeddedJet();
564  }
565  }
566 
567  return kTRUE;
568 }
569 
570 //________________________________________________________________________
572 {
573  // Get the next accepted embedded jet.
574 
575  if(reset)
576  fEmbJetsCont->ResetCurrentID();
577 
579  while (jet && jet->MCPt() < fMCJetPtThreshold) jet = fEmbJetsCont->GetNextAcceptJet();
580 
581  return jet;
582 }
583 
584 //________________________________________________________________________
586  AliParticleContainer* tracks, AliClusterContainer* clusters,
587  AliEmcalJet *jet, Bool_t bPartialExclusion) const
588 {
589  // Get rigid cone.
590 
591  eta = -999;
592  phi = -999;
593  pt = 0;
594 
595  if (!tracks && !clusters)
596  return;
597 
598  Float_t LJeta = 999;
599  Float_t LJphi = 999;
600 
601  if (jet) {
602  LJeta = jet->Eta();
603  LJphi = jet->Phi();
604  }
605 
606  Float_t maxEta = fConeMaxEta;
607  Float_t minEta = fConeMinEta;
608  Float_t maxPhi = fConeMaxPhi;
609  Float_t minPhi = fConeMinPhi;
610 
611  if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
612  if (minPhi < 0) minPhi = 0;
613 
614  Float_t dLJ = 0;
615  Int_t repeats = 0;
616  Bool_t reject = kTRUE;
617  do {
618  eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
619  phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
620  dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
621 
622  if(bPartialExclusion) {
623  reject = kFALSE;
624 
625  TRandom3 rnd;
626  rnd.SetSeed(0);
627 
628  Double_t ncoll = GetNColl();
629 
630  Double_t prob = 0.;
631  if(ncoll>0)
632  prob = 1./ncoll;
633 
634  if(rnd.Rndm()<=prob) reject = kTRUE; //reject cone
635  }
636 
637  repeats++;
638  } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
639 
640  if (repeats == 999) {
641  AliWarning(Form("%s: Could not get random cone!", GetName()));
642  return;
643  }
644 
645  if (clusters) {
646  clusters->ResetCurrentID();
647  AliVCluster* cluster = clusters->GetNextAcceptCluster();
648  while (cluster) {
649  TLorentzVector nPart;
650  cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
651 
652  Float_t cluseta = nPart.Eta();
653  Float_t clusphi = nPart.Phi();
654 
655  if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
656  clusphi += 2 * TMath::Pi();
657  if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
658  clusphi -= 2 * TMath::Pi();
659 
660  Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
661  if (d <= fConeRadius)
662  pt += nPart.Pt();
663 
664  cluster = clusters->GetNextAcceptCluster();
665  }
666  }
667 
668  if (tracks) {
669  tracks->ResetCurrentID();
670  AliVParticle* track = tracks->GetNextAcceptParticle();
671  while(track) {
672  Float_t tracketa = track->Eta();
673  Float_t trackphi = track->Phi();
674 
675  if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
676  trackphi += 2 * TMath::Pi();
677  if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
678  trackphi -= 2 * TMath::Pi();
679 
680  Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
681  if (d <= fConeRadius)
682  pt += track->Pt();
683 
684  track = tracks->GetNextAcceptParticle();
685  }
686  }
687 }
688 
689 //________________________________________________________________________
691 {
692  // Set default cuts for full cones
693 
695  SetConePhiLimits(1.4+fConeRadius,TMath::Pi()-fConeRadius);
696 }
697 
698 //________________________________________________________________________
700 {
701  // Set default cuts for charged cones
702 
704  SetConePhiLimits(-10, 10);
705 }
706 
707 //________________________________________________________________________
709 {
710  // Initialize the analysis.
711 
713 
714  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
715  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
716  if (fEmbTracksCont && fEmbTracksCont->GetArray() == 0) fEmbTracksCont = 0;
718  if (fRandTracksCont && fRandTracksCont->GetArray() == 0) fRandTracksCont = 0;
720  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
721  if (fEmbJetsCont && fEmbJetsCont->GetArray() == 0) fEmbJetsCont = 0;
722 
723  if (fRCperEvent < 0) {
725  Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
726  fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
727  if (fRCperEvent == 0)
728  fRCperEvent = 1;
729  }
730 
731  if (fMinRC2LJ < 0)
732  fMinRC2LJ = fConeRadius * 1.5;
733 
734  const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
735  if (fMinRC2LJ > maxDist) {
736  AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
737  "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
738  fMinRC2LJ = maxDist;
739  }
740 }
741 
742 //________________________________________________________________________
744  // Get NColl - returns value of corresponding bin
745  // only works for pA
746  // values taken from V0A slicing https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Tables_with_centrality_bins_from
747 
748  if(fBeamType == kpA) {
749 
750  const Int_t nNCollBins = 7;
751 
752  Double_t centMin[nNCollBins] = {0.,5.,10.,20.,40.,60.,80.};
753  Double_t centMax[nNCollBins] = {5.,10.,20.,40.,60.,80.,100.};
754 
755  Double_t nColl[nNCollBins] = {14.7,13.,11.7,9.38,6.49,3.96,1.52};
756 
757  for(Int_t i = 0; i<nNCollBins; i++) {
758  if(fCent>=centMin[i] && fCent<centMax[i])
759  return nColl[i];
760  }
761 
762  return -1.;
763  }
764  else {
765  AliWarning(Form("%s: Only works for pA analysis. Returning -1",GetName()));
766  return -1.;
767  }
768 }
AliParticleContainer * fRandTracksCont
Embedded clusters.
AliParticleContainer * fEmbTracksCont
Embedded jets.
TH2 ** fHistEmbPartPtvsJetPt
Pt-rho*A vs. area of embedded jets.
Double_t Area() const
Definition: AliEmcalJet.h:130
virtual AliVParticle * GetNextAcceptParticle()
TH1 ** fHistDistLeadPart2JetAxis
Pt vs jet pt - rho*A.
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
Definition: External.C:260
AliClusterContainer * fCaloClustersCont
Tracks.
TH1 ** fHistDeltaPtRCRand
deltaPt = Pt(RC) - A * rho, imposing min distance from leading jet with 1/ncoll probability ...
Definition: External.C:236
Double_t MCPt() const
Definition: AliEmcalJet.h:153
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
TH2 ** fHistEmbBkgArea
Distance between leading particle and jet axis.
TH2 * fHistRCPtExPartialLJVSDPhiLJ
Random cone pt, imposing min distance from leading jet, vs. deltaPhi leading jet. ...
Double_t Eta() const
Definition: AliEmcalJet.h:121
TH1 ** fHistRCPt
Phi-Eta distribution of random cones.
TH2 * fHistRCPtExLJVSDPhiLJ
deltaPt = Pt(embjet) - Area(embjet) * rho - Pt(embtrack) vs. event plane
TH1 ** fHistDeltaPtRCExPartialLJ
deltaPt = Pt(RC) - A * rho, imposing min distance from leading jet
Double_t Phi() const
Definition: AliEmcalJet.h:117
TH2 * fHistEmbJetsPhiEta
Random cone pt, imposing min distance from leading jet, vs. deltaPhi leading jet with 1/ncoll probabi...
Double_t fMinBinPt
min pt in histograms
Double_t fEPV0
!event plane V0
Int_t fCentBin
!event centrality bin
AliClusterContainer * fRandCaloClustersCont
Randomized tracks.
TH2 ** fHistDeltaPtEmbvsEP
deltaPt = Pt(embjet) - Area(embjet) * rho - Pt(embtrack) vs. Area(embjet)
TRandom * gRandom
Container for particles within the EMCAL framework.
TH1 ** fHistRCPtRand
Random cone pt, imposing min distance from leading jet with 1/ncoll probability.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TH3 ** fHistEmbJetsPtArea
deltaPt = Pt(RC) - A * rho, randomzied particles
AliEmcalJet * GetLeadingJet(const char *opt="")
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
TH2 ** fHistDeltaPtRCvsEP
Area(RC) * rho vs. Pt(RC)
int Int_t
Definition: External.C:63
TH1 ** fHistRCPtExLJ
Random cone pt.
float Float_t
Definition: External.C:68
AliEmcalJet * NextEmbeddedJet(Bool_t reset=kFALSE)
TH3 ** fHistEmbJetsCorrPtArea
Pt vs. area of embedded jets.
TH2 ** fHistRhoVSEmbBkg
Pt(embjet) - Pt(embtrack) vs. area of embedded jets.
TH1 ** fHistDeltaPtRCExLJ
deltaPt = Pt(RC) - A * rho vs. event plane
TH1 ** fHistRCPtExPartialLJ
Random cone pt, imposing min distance from leading jet.
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
AliClusterContainer * fEmbCaloClustersCont
Embedded tracks.
BeamType fBeamType
!event beam type
Double_t fCent
!event centrality
TH2 ** fHistDeltaPtEmbArea
Area(embjet) * rho vs. Pt(embjet) - Pt(embtrack)
TH2 * fHistRCPhiEta
Rho vs. centrality.
TH2 ** fHistJetPtvsJetCorrPt
MC jet pt total jet pt - rho*A.
AliEmcalJet * GetNextAcceptJet()
void SetConePhiLimits(Float_t min, Float_t max)
AliParticleContainer * fTracksCont
Jets.
Double_t Pt() const
Definition: AliEmcalJet.h:109
void GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer *tracks, AliClusterContainer *clusters, AliEmcalJet *jet=0, Bool_t bPartialExclusion=0) const
Double_t centMax
Bool_t FillHistograms()
Function filling histograms.
TH2 * fHistRhovsCent
Randomized clusters.
void SetConeEtaLimits(Float_t min, Float_t max)
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
Definition: External.C:220
Double_t fVertex[3]
!event vertex
void SetMakeGeneralHistograms(Bool_t g)
AliJetContainer * fEmbJetsCont
Clusters.
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
void ExecOnce()
Perform steps needed to initialize the analysis.
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
TH2 ** fHistEmbPartPtvsJetCorrPt
MC jet pt total jet pt.
Container structure for EMCAL clusters.
AliVCluster * GetNextAcceptCluster()
TH2 * fHistLeadPartPhiEta
Phi-Eta distribution of embedded jets.
Definition: External.C:196
Int_t fNbins
no. of pt bins
TH2 ** fHistRhoVSRCPt
Random cone pt, randomized particles.
Bool_t reject
Double_t centMin