AliPhysics  9b6b435 (9b6b435)
AliAnalysisTaskEmcalJetSpectraQA.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <TClonesArray.h>
17 #include <TH2F.h>
18 #include <TH3F.h>
19 #include <THnSparse.h>
20 #include <THashList.h>
21 
22 #include <AliAnalysisManager.h>
23 #include <AliVEventHandler.h>
24 
25 #include "AliEmcalJet.h"
26 #include "AliVCluster.h"
27 #include "AliVParticle.h"
28 #include "AliLog.h"
29 #include "AliJetContainer.h"
30 #include "AliClusterContainer.h"
31 #include "AliParticleContainer.h"
32 
34 
35 // Definitions of class AliAnalysisTaskEmcalJetSpectraQA::AliEmcalJetInfo
36 
40  fArea(0),
41  fMCPt(0),
42  fNConstituents(0),
43  fNEF(0),
44  fCent(0),
45  fEP(0),
46  fCorrPt(0),
47  fZ(0),
48  fLeadingPt(0)
49 {
50 }
51 
57  fArea(jet.Area()),
58  fMCPt(jet.MCPt()),
59  fNConstituents(jet.GetNumberOfConstituents()),
60  fNEF(jet.NEF()),
61  fCent(0),
62  fEP(0),
63  fCorrPt(0),
64  fZ(0),
65  fLeadingPt(0)
66 {
67  jet.GetMomentum(*this);
68 }
69 
70 // Definitions of class AliAnalysisTaskEmcalJetSpectraQA
71 
75 
78  AliAnalysisTaskEmcalJetLight("AliAnalysisTaskEmcalJetSpectraQA", kTRUE),
80  fJetEPaxis(kFALSE),
81  fAreaAxis(kTRUE),
82  fPtBinWidth(0.5),
83  fMaxPt(250),
84  fIsEmbedded(kFALSE),
85  fHistManager()
86 
87 {
89 }
90 
95  AliAnalysisTaskEmcalJetLight(name, kTRUE),
97  fJetEPaxis(kFALSE),
98  fAreaAxis(kTRUE),
99  fPtBinWidth(0.5),
100  fMaxPt(250),
101  fIsEmbedded(kFALSE),
102  fHistManager(name)
103 {
105 }
106 
111 {
112  AliError("Tree output not implemented. Falling back to THnSparse output");
113  AllocateTHnSparse(jets);
114 }
115 
120 {
121  Double_t jetRadius = jets->GetJetRadius();
122 
123  TString title[30]= {""};
124  Int_t nbins[30] = {0};
125  Double_t min[30] = {0.};
126  Double_t max[30] = {0.};
127  Int_t dim = 0;
128 
129  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
130 
131  if (fForceBeamType != kpp) {
132  title[dim] = "Centrality (%)";
133  nbins[dim] = 20;
134  min[dim] = 0;
135  max[dim] = 100;
136  dim++;
137 
138  if (fJetEPaxis) {
139  title[dim] = "#phi_{jet} - #psi_{EP}";
140  nbins[dim] = nPtBins/5;
141  min[dim] = 0;
142  max[dim] = TMath::Pi();
143  dim++;
144  }
145  }
146 
147  title[dim] = "#eta_{jet}";
148  nbins[dim] = nPtBins/10;
149  min[dim] = -1;
150  max[dim] = 1;
151  dim++;
152 
153  title[dim] = "#phi_{jet} (rad)";
154  nbins[dim] = nPtBins/10*3;
155  min[dim] = 0;
156  max[dim] = 2*TMath::Pi();
157  dim++;
158 
159  title[dim] = "#it{p}_{T} (GeV/#it{c})";
160  nbins[dim] = nPtBins;
161  min[dim] = 0;
162  max[dim] = fMaxPt;
163  dim++;
164 
165  if (fIsEmbedded) {
166  title[dim] = "#it{p}_{T}^{MC} (GeV/#it{c})";
167  nbins[dim] = nPtBins;
168  min[dim] = 0;
169  max[dim] = fMaxPt;
170  dim++;
171  }
172 
173  if (fForceBeamType != kpp) {
174  title[dim] = "#it{p}_{T}^{corr} (GeV/#it{c})";
175  nbins[dim] = nPtBins;
176  min[dim] = -fMaxPt/2;
177  max[dim] = fMaxPt/2;
178  dim++;
179  }
180 
181  if (fAreaAxis) {
182  // area resolution is about 0.01 (w/ ghost area 0.005)
183  // for nPtBins = 250 use bin width 0.01
184  title[dim] = "#it{A}_{jet}";
185  nbins[dim] = TMath::CeilNint(2.0*jetRadius*jetRadius*TMath::Pi() / 0.01 * nPtBins / 250);
186  min[dim] = 0;
187  max[dim] = 2.0*jetRadius*jetRadius*TMath::Pi();
188  dim++;
189  }
190 
191  if (fClusterCollArray.size() > 0 && fParticleCollArray.size() > 0) {
192  title[dim] = "NEF";
193  nbins[dim] = nPtBins/5;
194  min[dim] = 0;
195  max[dim] = 1.0;
196  dim++;
197  }
198 
199  title[dim] = "#it{z}_{leading}";
200  nbins[dim] = nPtBins/5;
201  min[dim] = 0;
202  max[dim] = 1.0;
203  dim++;
204 
205  if (fForceBeamType != kpp) {
206  title[dim] = "No. of constituents";
207  nbins[dim] = 125;
208  min[dim] = 0;
209  max[dim] = 250;
210  dim++;
211  }
212  else {
213  title[dim] = "No. of constituents";
214  nbins[dim] = 50;
215  min[dim] = -0.5;
216  max[dim] = 49.5;
217  dim++;
218  }
219 
220  title[dim] = "#it{p}_{T,particle}^{leading} (GeV/#it{c})";
221  nbins[dim] = nPtBins/10*3;
222  min[dim] = 0;
223  max[dim] = 150;
224  dim++;
225 
226  TString histname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
227  THnSparse* hn = fHistManager.CreateTHnSparse(histname.Data(), histname.Data(), dim, nbins, min, max);
228  for (Int_t i = 0; i < dim; i++) {
229  hn->GetAxis(i)->SetTitle(title[i]);
230  }
231 }
232 
237 {
238  TString histname;
239  TString title;
240 
241  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
242 
243  for (Int_t i = 0; i < GetNCentBins(); i++) {
244  histname = TString::Format("%s/fHistJetPtEtaPhi_%d", jets->GetArrayName().Data(), i);
245  title = histname + ";#it{p}_{T} (GeV/#it{c});#eta;#phi (rad)";
246  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, -1, 1, 41, 0, 2*TMath::Pi()*41/40, nPtBins, 0, fMaxPt);
247 
248  histname = TString::Format("%s/fHistJetPtArea_%d", jets->GetArrayName().Data(), i);
249  title = histname + ";#it{p}_{T} (GeV/#it{c});#it{A}_{jet};counts";
250  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 150, 0, 1.5);
251 
252  histname = TString::Format("%s/fHistJetPtEP_%d", jets->GetArrayName().Data(), i);
253  title = histname + ";#it{p}_{T} (GeV/#it{c});#phi_{jet} - #psi_{EP};counts";
254  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 100, 0, TMath::Pi());;
255 
256  histname = TString::Format("%s/fHistJetPtNEF_%d", jets->GetArrayName().Data(), i);
257  title = histname + ";#it{p}_{T} (GeV/#it{c});NEF;counts";
258  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 102, 0, 1.02);
259 
260  histname = TString::Format("%s/fHistJetPtZ_%d", jets->GetArrayName().Data(), i);
261  title = histname + ";#it{p}_{T} (GeV/#it{c});#it{z}_{leading};counts";
262  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 102, 0, 1.02);
263 
264  histname = TString::Format("%s/fHistJetPtLeadingPartPt_%d", jets->GetArrayName().Data(), i);
265  title = histname + ";#it{p}_{T} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c});counts";
266  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 120, 0, 120);
267 
268  if (!jets->GetRhoName().IsNull()) {
269  histname = TString::Format("%s/fHistJetCorrPtEtaPhi_%d", jets->GetArrayName().Data(), i);
270  title = histname + ";#it{p}_{T,corr} (GeV/#it{c});#eta;#phi (rad)";
271  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, -1, 1, 41, 0, 2*TMath::Pi()*201/200, nPtBins*2, -fMaxPt, fMaxPt);
272 
273  histname = TString::Format("%s/fHistJetCorrPtArea_%d", jets->GetArrayName().Data(), i);
274  title = histname + ";#it{p}_{T,corr} (GeV/#it{c});#it{A}_{jet};counts";
275  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 150, 0, 1.5);
276 
277  histname = TString::Format("%s/fHistJetCorrPtEP_%d", jets->GetArrayName().Data(), i);
278  title = histname + ";#it{p}_{T,corr} (GeV/#it{c});#phi_{jet} - #psi_{EP};counts";
279  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 100, 0, TMath::Pi());;
280 
281  histname = TString::Format("%s/fHistJetCorrPtNEF_%d", jets->GetArrayName().Data(), i);
282  title = histname + ";#it{p}_{T,corr} (GeV/#it{c});NEF;counts";
283  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 102, 0, 1.02);
284 
285  histname = TString::Format("%s/fHistJetCorrPtZ_%d", jets->GetArrayName().Data(), i);
286  title = histname + ";#it{p}_{T,corr} (GeV/#it{c});#it{z}_{leading};counts";
287  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 102, 0, 1.02);
288 
289  histname = TString::Format("%s/fHistJetCorrPtLeadingPartPt_%d", jets->GetArrayName().Data(), i);
290  title = histname + ";#it{p}_{T,corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c});counts";
291  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 120, 0, 120);
292 
293  histname = TString::Format("%s/fHistJetPtCorrPt_%d", jets->GetArrayName().Data(), i);
294  title = histname + ";#it{p}_{T} (GeV/#it{c});#it{p}_{T,corr} (GeV/#it{c});counts";
295  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins*2, -fMaxPt, fMaxPt);
296 
297  if (fIsEmbedded) {
298  histname = TString::Format("%s/fHistJetMCPtCorrPt_%d", jets->GetArrayName().Data(), i);
299  title = histname + ";#it{p}_{T,MC} (GeV/#it{c});#it{p}_{T,corr} (GeV/#it{c});counts";
300  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins*2, -fMaxPt, fMaxPt);
301  }
302  }
303 
304  if (fIsEmbedded) {
305  histname = TString::Format("%s/fHistJetPtMCPt_%d", jets->GetArrayName().Data(), i);
306  title = histname + ";#it{p}_{T} (GeV/#it{c});#it{p}_{T,MC} (GeV/#it{c});counts";
307  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
308  }
309  }
310 }
311 
314 {
316 
317  Int_t maxTracks = 6000;
318  Int_t constituentsNbins = 250;
319  Double_t constituentsMax = 249.5;
320  Double_t maxRho = 500;
321 
322  if (fForceBeamType == kpp) {
323  constituentsNbins = 50;
324  constituentsMax = 49.5;
325  maxRho = 50;
326  maxTracks = 200;
327  }
328  else if (fForceBeamType == kpA) {
329  constituentsNbins = 100;
330  constituentsMax = 99.5;
331  maxRho = 200;
332  maxTracks = 500;
333  }
334 
335  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
336 
337  TString histname;
338  TString title;
339 
340  for (auto cont_it : fJetCollArray) {
341  AliJetContainer* jets = cont_it.second;
343 
344  switch (fHistoType) {
345  case kTH2:
346  AllocateTHX(jets);
347  break;
348  case kTHnSparse:
349  AllocateTHnSparse(jets);
350  break;
351  case kTTree:
352  AllocateTTree(jets);
353  break;
354  }
355 
356  TString histname;
357 
358  for (Int_t i = 0; i < GetNCentBins(); i++) {
359  if (jets->GetParticleContainer()) {
360  histname = TString::Format("%s/fHistTracksJetPt_%d", jets->GetArrayName().Data(), i);
361  title = histname + ";#it{p}_{T,track} (GeV/#it{c});#it{p}_{T,jet} (GeV/#it{c});counts";
362  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins / 2, 0, fMaxPt / 2, nPtBins, 0, fMaxPt);
363 
364  histname = TString::Format("%s/fHistTracksPtDist_%d", jets->GetArrayName().Data(), i);
365  title = histname + ";#it{p}_{T,track} (GeV/#it{c});#it{d};counts";
366  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins / 2, 0, fMaxPt / 2, 100, 0, 5);
367 
368  histname = TString::Format("%s/fHistTracksZJetPtJetConst_%d", jets->GetArrayName().Data(), i);
369  title = histname + ";#it{z}_{track} (GeV/#it{c});#it{d};No. of constituents";
370  fHistManager.CreateTH3(histname.Data(), title.Data(), 120, 0.0, 1.2, nPtBins, 0, fMaxPt, constituentsNbins, -0.5, constituentsMax);
371  }
372 
373  if (jets->GetClusterContainer()) {
374  histname = TString::Format("%s/fHistClustersJetPt_%d", jets->GetArrayName().Data(), i);
375  title = histname + ";#it{p}_{T,cluster} (GeV/#it{c});#it{p}_{T,jet} (GeV/#it{c});counts";
376  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins / 2, 0, fMaxPt / 2, nPtBins, 0, fMaxPt);
377 
378  histname = TString::Format("%s/fHistClustersPtDist_%d", jets->GetArrayName().Data(), i);
379  title = histname + ";#it{p}_{T,cluster} (GeV/#it{c});#it{d};counts";
380  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins / 2, 0, fMaxPt / 2, 100, 0, 5);
381 
382  histname = TString::Format("%s/fHistClustersZJetPtJetConst_%d", jets->GetArrayName().Data(), i);
383  title = histname + ";#it{z}_{cluster} (GeV/#it{c});#it{p}_{T,jet} (GeV/#it{c});No. of constituents";
384  fHistManager.CreateTH3(histname.Data(), title.Data(), 120, 0.0, 1.2, nPtBins, 0, fMaxPt, constituentsNbins, -0.5, constituentsMax);
385  }
386 
387  histname = TString::Format("%s/fHistRejectionReason_%d", jets->GetArrayName().Data(), i);
388  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
389  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 100, 0, 250);
390  SetRejectionReasonLabels(hist->GetXaxis());
391  }
392 
393  histname = TString::Format("%s/fHistLeadJetPtVsCent", jets->GetArrayName().Data());
394  title = histname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
395  fHistManager.CreateTH2(histname.Data(), title.Data(), 100, 0, 100, nPtBins, 0, fMaxPt);
396 
397  histname = TString::Format("%s/fHistLeadJetPtVsNTracks", jets->GetArrayName().Data());
398  title = histname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
399  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, 0, maxTracks, nPtBins, 0, fMaxPt);
400 
401  if (!jets->GetRhoName().IsNull()) {
402  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
403  title = histname + ";Centrality (%);#rho (GeV/#it{c} rad^{-1});counts";
404  fHistManager.CreateTH2(histname.Data(), title.Data(), 100, 0, 100, 1000, 0, maxRho);
405 
406  histname = TString::Format("%s/fHistRhoVsNTracks", jets->GetArrayName().Data());
407  title = histname + ";no. of tracks;#rho (GeV/#it{c} rad^{-1});counts";
408  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, 0, maxTracks, 1000, 0, maxRho);
409 
410  histname = TString::Format("%s/fHistRhoVsLeadJetPt", jets->GetArrayName().Data());
411  title = histname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} rad^{-1});counts";
412  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 1000, 0, maxRho);
413  }
414  }
415 
416  TIter nextElement(fHistManager.GetListOfHistograms());
417  TObject* obj = 0;
418  while ((obj = nextElement())) fOutput->Add(obj);
419  PostData(1, fOutput);
420 }
421 
426 {
427  TString histname;
428 
429  if (fCentBin < 0) {
430  AliError(Form("fCentBin is %d! fCent = %.3f. Fix the centrality bins to include all possible values of centrality.", fCentBin, fCent));
431  return kFALSE;
432  }
433 
434  for (auto cont_it : fJetCollArray) {
435  AliJetContainer* jets = cont_it.second;
436  Double_t rhoVal = 0;
437  if (jets->GetRhoParameter()) rhoVal = jets->GetRhoVal();
438 
439  Double_t leadJetPt = 0;
440  for (auto jet : jets->accepted()) {
441 
442  UInt_t rejectionReason = 0;
443  if (!jets->AcceptJet(jet, rejectionReason)) {
444  histname = TString::Format("%s/fHistRejectionReason_%d", jets->GetArrayName().Data(), fCentBin);
445  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
446  continue;
447  }
448 
449  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
450  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
451 
452  TLorentzVector leadPart;
453 
454  jets->GetLeadingHadronMomentum(leadPart, jet);
455 
456  // Fill THnSparse
457  Double_t ep = jet->Phi() - fEPV0;
458  while (ep < 0) ep += TMath::Pi();
459  while (ep >= TMath::Pi()) ep -= TMath::Pi();
460 
461  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
462  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
463 
464  AliEmcalJetInfo jetInfo(*jet);
465  jetInfo.fCent = fCent;
466  jetInfo.fEP = ep;
467  jetInfo.fCorrPt = corrPt;
468  jetInfo.fZ = z;
469  jetInfo.fLeadingPt = ptLeading;
470 
471  if (jet->Pt() > leadJetPt) leadJetPt = jet->Pt();
472 
473  FillJetHisto(jetInfo, jets);
474 
475  AliParticleContainer* tracks = jets->GetParticleContainer();
476  if (tracks) {
477  for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
478  AliVParticle *track = jet->TrackAt(it, tracks->GetArray());
479  if (track) {
480  histname = TString::Format("%s/fHistTracksJetPt_%d", jets->GetArrayName().Data(), fCentBin);
481  fHistManager.FillTH2(histname.Data(), track->Pt(), jet->Pt());
482 
483  Double_t dphi = TVector2::Phi_0_2pi(track->Phi() - jet->Phi());
484  Double_t deta = track->Eta() - jet->Eta();
485  Double_t dist = TMath::Sqrt(deta * deta + dphi * dphi);
486 
487  histname = TString::Format("%s/fHistTracksPtDist_%d", jets->GetArrayName().Data(), fCentBin);
488  fHistManager.FillTH2(histname.Data(), track->Pt(), dist);
489 
490  histname = TString::Format("%s/fHistTracksZJetPtJetConst_%d", jets->GetArrayName().Data(), fCentBin);
491  fHistManager.FillTH3(histname.Data(), GetParallelFraction(track, jet), jet->Pt(), jet->GetNumberOfConstituents());
492  }
493  }
494  }
495 
496  AliClusterContainer* clusters = jets->GetClusterContainer();
497  if (clusters) {
498  for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
499  AliVCluster *cluster = jet->ClusterAt(ic, clusters->GetArray());
500 
501  if (cluster) {
502  TLorentzVector nPart;
503  if (clusters->GetDefaultClusterEnergy() >=0 && clusters->GetDefaultClusterEnergy() < AliVCluster::kLastUserDefEnergy) {
504  cluster->GetMomentum(nPart, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
505  }
506  else {
507  cluster->GetMomentum(nPart, fVertex);
508  }
509 
510  histname = TString::Format("%s/fHistClustersJetPt_%d", jets->GetArrayName().Data(), fCentBin);
511  fHistManager.FillTH2(histname.Data(), nPart.Pt(), jet->Pt());
512 
513  Double_t dphi = TVector2::Phi_0_2pi(nPart.Phi() - jet->Phi());
514  Double_t deta = nPart.Eta() - jet->Eta();
515  Double_t dist = TMath::Sqrt(deta * deta + dphi * dphi);
516 
517  histname = TString::Format("%s/fHistClustersPtDist_%d", jets->GetArrayName().Data(), fCentBin);
518  fHistManager.FillTH2(histname.Data(), nPart.Pt(), dist);
519 
520  histname = TString::Format("%s/fHistClustersZJetPtJetConst_%d", jets->GetArrayName().Data(), fCentBin);
521  fHistManager.FillTH3(histname.Data(), GetParallelFraction(nPart.Vect(), jet), jet->Pt(), jet->GetNumberOfConstituents());
522  }
523  }
524  }
525  } //jet loop
526 
527  Int_t ntracks = 0;
528  for (auto cont : this->fParticleCollArray) ntracks += cont.second->GetNAcceptEntries();
529 
530  histname = TString::Format("%s/fHistLeadJetPtVsCent", jets->GetArrayName().Data());
531  fHistManager.FillTH2(histname.Data(), fCent, leadJetPt);
532 
533  histname = TString::Format("%s/fHistLeadJetPtVsNTracks", jets->GetArrayName().Data());
534  fHistManager.FillTH2(histname.Data(), ntracks, leadJetPt);
535 
536  if (jets->GetRhoParameter()) {
537  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
538  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
539 
540  histname = TString::Format("%s/fHistRhoVsNTracks", jets->GetArrayName().Data());
541  fHistManager.FillTH2(histname.Data(), ntracks, rhoVal);
542 
543  histname = TString::Format("%s/fHistRhoVsLeadJetPt", jets->GetArrayName().Data());
544  fHistManager.FillTH2(histname.Data(), leadJetPt, rhoVal);
545  }
546  }
547  return kTRUE;
548 }
549 
555 {
556  TString histname;
557 
558  histname = TString::Format("%s/fHistJetPtEtaPhi_%d", jets->GetArrayName().Data(), fCentBin);
559  fHistManager.FillTH3(histname.Data(), jet.Eta(), jet.Phi_0_2pi(), jet.Pt());
560 
561  histname = TString::Format("%s/fHistJetPtArea_%d", jets->GetArrayName().Data(), fCentBin);
562  fHistManager.FillTH2(histname.Data(), jet.Pt(), jet.fArea);
563 
564  histname = TString::Format("%s/fHistJetPtEP_%d", jets->GetArrayName().Data(), fCentBin);
565  fHistManager.FillTH2(histname.Data(), jet.Pt(), jet.fEP);
566 
567  histname = TString::Format("%s/fHistJetPtNEF_%d", jets->GetArrayName().Data(), fCentBin);
568  fHistManager.FillTH2(histname.Data(), jet.Pt(), jet.fNEF);
569 
570  histname = TString::Format("%s/fHistJetPtZ_%d", jets->GetArrayName().Data(), fCentBin);
571  fHistManager.FillTH2(histname.Data(), jet.Pt(), jet.fZ);
572 
573  histname = TString::Format("%s/fHistJetPtLeadingPartPt_%d", jets->GetArrayName().Data(), fCentBin);
574  fHistManager.FillTH2(histname.Data(), jet.Pt(), jet.fLeadingPt);
575 
576  if (fIsEmbedded) {
577  histname = TString::Format("%s/fHistJetPtMCPt_%d", jets->GetArrayName().Data(), fCentBin);
578  fHistManager.FillTH2(histname.Data(), jet.Pt(), jet.fMCPt);
579  }
580 
581  if (!jets->GetRhoName().IsNull()) {
582  histname = TString::Format("%s/fHistJetCorrPtEtaPhi_%d", jets->GetArrayName().Data(), fCentBin);
583  fHistManager.FillTH3(histname.Data(), jet.Eta(), jet.Phi_0_2pi(), jet.fCorrPt);
584 
585  histname = TString::Format("%s/fHistJetCorrPtArea_%d", jets->GetArrayName().Data(), fCentBin);
586  fHistManager.FillTH2(histname.Data(), jet.fCorrPt, jet.fArea);
587 
588  histname = TString::Format("%s/fHistJetCorrPtEP_%d", jets->GetArrayName().Data(), fCentBin);
589  fHistManager.FillTH2(histname.Data(), jet.fCorrPt, jet.fEP);
590 
591  histname = TString::Format("%s/fHistJetCorrPtNEF_%d", jets->GetArrayName().Data(), fCentBin);
592  fHistManager.FillTH2(histname.Data(), jet.fCorrPt, jet.fNEF);
593 
594  histname = TString::Format("%s/fHistJetCorrPtZ_%d", jets->GetArrayName().Data(), fCentBin);
595  fHistManager.FillTH2(histname.Data(), jet.fCorrPt, jet.fZ);
596 
597  histname = TString::Format("%s/fHistJetCorrPtLeadingPartPt_%d", jets->GetArrayName().Data(), fCentBin);
598  fHistManager.FillTH2(histname.Data(), jet.fCorrPt, jet.fLeadingPt);
599 
600  histname = TString::Format("%s/fHistJetPtCorrPt_%d", jets->GetArrayName().Data(), fCentBin);
601  fHistManager.FillTH2(histname.Data(), jet.Pt(), jet.fCorrPt);
602 
603  if (fIsEmbedded) {
604  histname = TString::Format("%s/fHistJetMCPtCorrPt_%d", jets->GetArrayName().Data(), fCentBin);
605  fHistManager.FillTH2(histname.Data(), jet.fMCPt, jet.fCorrPt);
606  }
607  }
608 }
609 
615 {
616  AliError("Tree output not implemented. Falling back to THnSparse output");
617  FillTHnSparse(jet, jets);
618 }
619 
625 {
626  TString histname;
627  Double_t contents[30]={0};
628 
629  histname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
630  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
631 
632  if (!histJetObservables) return;
633 
634  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
635  TString title(histJetObservables->GetAxis(i)->GetTitle());
636  if (title=="Centrality (%)")
637  contents[i] = jet.fCent;
638  else if (title=="#phi_{jet} - #psi_{EP}")
639  contents[i] = jet.fEP;
640  else if (title=="#eta_{jet}")
641  contents[i] = jet.Eta();
642  else if (title=="#phi_{jet} (rad)")
643  contents[i] = jet.Phi_0_2pi();
644  else if (title=="#it{p}_{T} (GeV/#it{c})")
645  contents[i] = jet.Pt();
646  else if (title=="#it{p}_{T}^{MC} (GeV/#it{c})")
647  contents[i] = jet.fMCPt;
648  else if (title=="#it{p}_{T}^{corr} (GeV/#it{c})")
649  contents[i] = jet.fCorrPt;
650  else if (title=="#it{A}_{jet}")
651  contents[i] = jet.fArea;
652  else if (title=="NEF")
653  contents[i] = jet.fNEF;
654  else if (title=="#it{z}_{leading}")
655  contents[i] = jet.fZ;
656  else if (title=="No. of constituents")
657  contents[i] = jet.fNConstituents;
658  else if (title=="#it{p}_{T,particle}^{leading} (GeV/#it{c})")
659  contents[i] = jet.fLeadingPt;
660  else
661  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
662  }
663 
664  histJetObservables->Fill(contents);
665 }
666 
672 {
673  switch (fHistoType) {
674  case kTH2:
675  FillTHX(jet, jets);
676  break;
677 
678  case kTHnSparse:
679  FillTHnSparse(jet, jets);
680  break;
681 
682  case kTTree:
683  FillTTree(jet, jets);
684  break;
685  }
686 }
687 
697  Double_t trackPtCut, Double_t clusECut, TString suffix)
698 {
699  // Get the pointer to the existing analysis manager via the static access method.
700  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
701  if (!mgr) {
702  ::Error("AddTaskEmcalJetSpectraQA", "No analysis manager to connect to.");
703  return nullptr;
704  }
705 
706  // Check the analysis type using the event handlers connected to the analysis manager.
707  AliVEventHandler* handler = mgr->GetInputEventHandler();
708  if (!handler) {
709  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
710  return nullptr;
711  }
712 
714 
715  if (handler->InheritsFrom("AliESDInputHandler")) {
716  dataType = kESD;
717  }
718  else if (handler->InheritsFrom("AliAODInputHandler")) {
719  dataType = kAOD;
720  }
721 
722  // Init the task and do settings
723 
724  if (trackName == "usedefault") {
725  if (dataType == kESD) {
726  trackName = "Tracks";
727  }
728  else if (dataType == kAOD) {
729  trackName = "tracks";
730  }
731  else {
732  trackName = "";
733  }
734  }
735 
736  if (clusName == "usedefault") {
737  if (dataType == kESD) {
738  clusName = "CaloClusters";
739  }
740  else if (dataType == kAOD) {
741  clusName = "caloClusters";
742  }
743  else {
744  clusName = "";
745  }
746  }
747 
748  TString name("AliAnalysisTaskEmcalJetSpectraQA");
749  if (strcmp(suffix,"")) {
750  name += "_";
751  name += suffix;
752  }
753 
755  jetTask->SetVzRange(-10,10);
756  jetTask->SetNeedEmcalGeom(kFALSE);
757  AliParticleContainer *partCont = jetTask->AddParticleContainer(trackName.Data());
758  if (partCont) partCont->SetParticlePtCut(trackPtCut);
759 
760  AliClusterContainer *clusterCont = jetTask->AddClusterContainer(clusName.Data());
761  if (clusterCont) {
762  clusterCont->SetClusECut(0.);
763  clusterCont->SetClusPtCut(0.);
764  clusterCont->SetClusHadCorrEnergyCut(clusECut);
765  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
766  }
767 
768  // Final settings, pass to manager and set the containers
769  mgr->AddTask(jetTask);
770 
771  // Create containers for input/output
772  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
773  TString contname(name);
774  contname += "_histos";
775  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
776  TList::Class(),AliAnalysisManager::kOutputContainer,
777  Form("%s", AliAnalysisManager::GetCommonFileName()));
778  mgr->ConnectInput (jetTask, 0, cinput1 );
779  mgr->ConnectOutput (jetTask, 1, coutput1 );
780 
781  return jetTask;
782 }
Declaration of class AliAnalysisTaskEmcalJetSpectraQA.
AliClusterContainer * AddClusterContainer(std::string branchName, std::string contName="")
THashList * CreateHistoGroup(const char *groupname)
Create a new group of histograms within a parent group.
void SetParticlePtCut(Double_t cut)
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
EDataType_t
Switch for the data type.
const char * title
Definition: MakeQAPdf.C:27
static AliAnalysisTaskEmcalJetSpectraQA * AddTaskEmcalJetSpectraQA(TString ntracks="usedefault", TString nclusters="usedefault", Double_t trackPtCut=0.15, Double_t clusECut=0.30, TString suffix="")
Double_t fEPV0
!event plane V0
Float_t fPtBinWidth
Histogram pt bin width.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
virtual void FillTHX(const AliEmcalJetInfo &jetInfo, const AliJetContainer *jets)
virtual void FillTHnSparse(const AliEmcalJetInfo &jetInfo, const AliJetContainer *jets)
AliClusterContainer * GetClusterContainer() const
void FillTH3(const char *hname, double x, double y, double z, double weight=1., Option_t *opt="")
Fill a 3D histogram within the container.
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
static UShort_t GetRejectionReasonBitPosition(UInt_t rejectionReason)
Returns the highest bit in the rejection map as reason why the object was rejected.
void UserCreateOutputObjects()
Overloads base class method. Creates output objects.
Container for particles within the EMCAL framework.
Int_t GetDefaultClusterEnergy() const
virtual void AllocateTTree(const AliJetContainer *jets)
void GetMomentum(TLorentzVector &vec) const
EBeamType_t fForceBeamType
forced beam type
const Int_t nPtBins
AliParticleContainer * GetParticleContainer() const
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
void FillJetHisto(const AliEmcalJetInfo &jetInfo, const AliJetContainer *jets)
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
THistManager fHistManager
Histogram manager.
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
float Float_t
Definition: External.C:68
Double_t Phi_0_2pi() const
Double_t fVertex[3]
!event vertex
TClonesArray * GetArray() const
virtual void FillTTree(const AliEmcalJetInfo &jetInfo, const AliJetContainer *jets)
Base task in the EMCAL jet framework (lighter version of AliAnalysisTaskEmcalJet) ...
Int_t fCentBin
!event centrality bin
Bool_t fJetEPaxis
whether a EP-jet axis should be included in the THnSparse
Bool_t fAreaAxis
whether the area axis should be included
AliRhoParameter * GetRhoParameter()
const TString & GetArrayName() const
AliParticleContainer * AddParticleContainer(std::string branchName, std::string contName="")
Float_t GetJetRadius() const
Definition: External.C:220
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
void SetVzRange(Double_t min, Double_t max)
Bool_t GetMomentum(TLorentzVector &mom, const AliVCluster *vc, Double_t mass) const
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
void SetClusPtCut(Double_t cut)
AliAnalysisTaskEmcalJetSpectraQA()
Default constructor for ROOT I/O purposes.
std::map< std::string, AliJetContainer * > fJetCollArray
jet collection array
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
void SetClusECut(Double_t cut)
void SetDefaultClusterEnergy(Int_t d)
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
Container structure for EMCAL clusters.
virtual void AllocateTHnSparse(const AliJetContainer *jets)
virtual void AllocateTHX(const AliJetContainer *jets)
Double_t fCent
!event centrality
Container for jet within the EMCAL jet framework.
Implementation of a task to perform QA on jet spectra.
TH3 * CreateTH3(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, int nbinsz, double zmin, double zmax, Option_t *opt="")
Create a new TH2 within the container.
void SetClusHadCorrEnergyCut(Double_t cut)
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)