AliPhysics  5a28df1 (5a28df1)
AliAnalysisTaskEmcalJetPerformance.cxx
Go to the documentation of this file.
1 /**********************************************************************************
2 * Copyright (C) 2016, Copyright Holders of the ALICE Collaboration *
3 * All rights reserved. *
4 * *
5 * Redistribution and use in source and binary forms, with or without *
6 * modification, are permitted provided that the following conditions are met: *
7 * * Redistributions of source code must retain the above copyright *
8 * notice, this list of conditions and the following disclaimer. *
9 * * Redistributions in binary form must reproduce the above copyright *
10 * notice, this list of conditions and the following disclaimer in the *
11 * documentation and/or other materials provided with the distribution. *
12 * * Neither the name of the <organization> nor the *
13 * names of its contributors may be used to endorse or promote products *
14 * derived from this software without specific prior written permission. *
15 * *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19 * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26 * *********************************************************************************/
27 
28 #include <vector>
29 
30 #include <TClonesArray.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3.h>
34 #include <TList.h>
35 #include <THnSparse.h>
36 #include <TRandom3.h>
37 #include <TGrid.h>
38 #include <TFile.h>
39 
40 #include <AliVCluster.h>
41 #include <AliVParticle.h>
42 #include <AliLog.h>
43 
44 #include "AliAnalysisManager.h"
45 #include <AliVEventHandler.h>
46 #include "AliTLorentzVector.h"
47 #include "AliEmcalJet.h"
48 #include "AliRhoParameter.h"
49 #include "AliJetContainer.h"
50 #include "AliParticleContainer.h"
51 #include "AliClusterContainer.h"
52 #include "AliEMCALGeometry.h"
53 #include "AliOADBContainer.h"
55 #include "AliEMCALTriggerPatchInfo.h"
57 #include "AliMCEvent.h"
58 
60 
64 
70  fPlotJetHistograms(kFALSE),
71  fPlotClusterHistograms(kFALSE),
72  fPlotParticleCompositionHistograms(kFALSE),
73  fComputeBackground(kFALSE),
74  fDoTriggerSimulation(kFALSE),
75  fPlotMatchedJetHistograms(kFALSE),
76  fComputeMBDownscaling(kFALSE),
77  fPlotDCal(kFALSE),
78  fMaxPt(250),
79  fNEtaBins(40),
80  fNPhiBins(200),
81  fNCentHistBins(0),
82  fCentHistBins(0),
83  fNPtHistBins(0),
84  fPtHistBins(0),
85  fNM02HistBins(0),
86  fM02HistBins(0),
87  fNEoverPBins(0),
88  fEoverPBins(0),
89  fTrackMatchingDeltaEtaMax(0.015),
90  fTrackMatchingDeltaPhiMax(0.030),
91  fMBUpscaleFactor(1.),
92  fMedianEMCal(0.),
93  fMedianDCal(0.),
94  fkEMCEJE(kFALSE),
95  fEmbeddingQA(),
96  fMinSharedMomentumFraction(0.5),
97  fMaxMatchedJetDistance(0.3),
98  fUseResponseMaker(kFALSE),
99  fMCJetContainer(nullptr),
100  fUseAliEventCuts(kTRUE),
101  fEventCuts(0),
102  fEventCutList(0),
103  fUseManualEventCuts(kFALSE),
104  fGeneratorLevel(0),
105  fHistManager()
106 {
108 }
109 
116  AliAnalysisTaskEmcalJet(name, kTRUE),
117  fPlotJetHistograms(kFALSE),
118  fPlotClusterHistograms(kFALSE),
120  fComputeBackground(kFALSE),
121  fDoTriggerSimulation(kFALSE),
123  fComputeMBDownscaling(kFALSE),
124  fPlotDCal(kFALSE),
125  fMaxPt(250),
126  fNEtaBins(40),
127  fNPhiBins(200),
128  fNCentHistBins(0),
129  fCentHistBins(0),
130  fNPtHistBins(0),
131  fPtHistBins(0),
132  fNM02HistBins(0),
133  fM02HistBins(0),
134  fNEoverPBins(0),
135  fEoverPBins(0),
138  fMBUpscaleFactor(1.),
139  fMedianEMCal(0.),
140  fMedianDCal(0.),
141  fkEMCEJE(kFALSE),
142  fEmbeddingQA(),
145  fUseResponseMaker(kFALSE),
147  fUseAliEventCuts(kTRUE),
148  fEventCuts(0),
149  fEventCutList(0),
150  fUseManualEventCuts(kFALSE),
151  fGeneratorLevel(0),
152  fHistManager(name)
153 {
155 }
156 
161 {
162 }
163 
168 {
169  fNCentHistBins = 4;
171  fCentHistBins[0] = 0;
172  fCentHistBins[1] = 10;
173  fCentHistBins[2] = 30;
174  fCentHistBins[3] = 50;
175  fCentHistBins[4] = 90;
176 
177  fNPtHistBins = 82;
180  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
181  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
182  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
183  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
184  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
185  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
186 
187  fNM02HistBins = 81;
189  GenerateFixedBinArray(35, 0, 0.7, fM02HistBins);
190  GenerateFixedBinArray(6, 0.7, 1., fM02HistBins+35);
191  GenerateFixedBinArray(20, 1., 3., fM02HistBins+41);
192  GenerateFixedBinArray(10, 3., 5., fM02HistBins+61);
193  GenerateFixedBinArray(10, 5., 10., fM02HistBins+71);
194 
195  fNEoverPBins = 47;
197  GenerateFixedBinArray(30, 0, 1.5, fEoverPBins);
198  GenerateFixedBinArray(10, 1.5, 3.5, fEoverPBins+30);
199  GenerateFixedBinArray(7, 3.5, 10.5, fEoverPBins+40);
200 }
201 
207 {
209 
210  // Intialize AliEventCuts
211  if (fUseAliEventCuts) {
212  fEventCutList = new TList();
213  fEventCutList ->SetOwner();
214  fEventCutList ->SetName("EventCutOutput");
215 
216  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
217  if(fUseManualEventCuts==1)
218  {
219  fEventCuts.SetManualMode();
220  fEventCuts.fMC = false;
221  fEventCuts.SetupLHC15o();
222  fEventCuts.fUseVariablesCorrelationCuts = true;
223  }
224  fEventCuts.AddQAplotsToList(fEventCutList);
225  fOutput->Add(fEventCutList);
226  }
227 
228  // Get the MC particle branch, in case it exists
229  fGeneratorLevel = GetMCParticleContainer("mcparticles");
230 
231  // Get MC jet container, in order to check the jet acceptance criteria
233  for (Int_t i=0; i<2; i++) {
234  auto jetCont = GetJetContainer(i);
235  TString jetContName = jetCont->GetName();
236  if (jetContName.Contains("mcparticles")) {
237  fMCJetContainer = jetCont;
238  }
239  }
240  if (!fMCJetContainer) {
241  Printf("No MC jet container found!");
242  }
243  Printf("mcJetContainer: %s", fMCJetContainer->GetName());
244  }
245 
246  // Allocate histograms
247  if (fPlotJetHistograms) {
249  }
252  }
255  }
256  if (fComputeBackground) {
258  }
259  if (fDoTriggerSimulation) {
261  }
264  }
265 
266  // Initialize embedding QA
268  if (embeddingHelper) {
269  bool res = fEmbeddingQA.Initialize();
270  if (res) {
272  }
273  }
274 
275  TIter next(fHistManager.GetListOfHistograms());
276  TObject* obj = 0;
277  while ((obj = next())) {
278  fOutput->Add(obj);
279  }
280 
281  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
282 }
283 
284 /*
285  * This function allocates the histograms for single jets.
286  * A set of histograms is allocated per each jet container.
287  */
289 {
290  TString histname;
291  TString title;
292 
293  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
294 
295  AliJetContainer* jets = 0;
296  TIter nextJetColl(&fJetCollArray);
297  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
298 
299  // Jet rejection reason
300  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
301  if (fForceBeamType != kpp) {
302  title = histname + ";Centrality (%);Rejection reason;#it{p}_{T,corr} (GeV/#it{c});counts";
303  TH3* hist = fHistManager.CreateTH3(histname.Data(), title.Data(), 10, 0, 100, 32, 0, 32, 50, 0, fMaxPt);
304  SetRejectionReasonLabels(hist->GetYaxis());
305  } else {
306  title = histname + ";Rejection reason;#it{p}_{T,corr} (GeV/#it{c});counts";
307  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
308  SetRejectionReasonLabels(hist->GetXaxis());
309  }
310 
311  // Rho vs. Centrality
312  if (!jets->GetRhoName().IsNull()) {
313  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
314  if (fForceBeamType != kpp) {
315  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
316  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
317  }
318  }
319 
320  // (Centrality, pT, NEF)
321  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
322  Int_t nbinsy = fMaxPt; Int_t miny = 0; Int_t maxy = fMaxPt;
323  Int_t nbinsz = 50; Int_t minz = 0; Int_t maxz = 1.;
324 
325  histname = TString::Format("%s/JetHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
326  if (fForceBeamType != kpp) {
327  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
328  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
329  }
330  else {
331  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
332  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
333  }
334 
335  if (fPlotDCal) {
336  histname = TString::Format("%s/JetHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
337  if (fForceBeamType != kpp) {
338  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
339  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
340  }
341  else {
342  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
343  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
344  }
345  }
346 
347  // (Centrality, pT upscaled, calo type)
348  if (fComputeMBDownscaling) {
349  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
350  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
351  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, 2, -0.5, 1.5, "s");
352  }
353 
354  // pT-leading vs. pT
355  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
356  if (fForceBeamType != kpp) {
357  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
358  fHistManager.CreateTH3(histname.Data(), title.Data(), 10, 0, 100, nPtBins, 0, fMaxPt, fMaxPt, 0, fMaxPt);
359  }
360  else {
361  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
362  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt, 0, fMaxPt);
363  }
364 
365  // A vs. pT
366  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
367  if (fForceBeamType != kpp) {
368  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
369  fHistManager.CreateTH3(histname.Data(), title.Data(), 10, 0, 100, nPtBins, 0, fMaxPt, 50, 0, 0.5);
370  }
371  else {
372  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
373  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 50, 0, 0.5);
374  }
375 
376  // (Centrality, pT, z-leading (charged))
377  nbinsx = 20; minx = 0; maxx = 100;
378  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
379  nbinsz = 50; minz = 0; maxz = 1.;
380 
381  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
382  if (fForceBeamType != kpp) {
383  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
384  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
385  }
386  else {
387  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
388  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
389  }
390 
391  if (fPlotDCal) {
392  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
393  if (fForceBeamType != kpp) {
394  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
395  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
396  }
397  else {
398  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
399  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
400  }
401  }
402 
403  // (Centrality, pT, z (charged))
404  nbinsx = 20; minx = 0; maxx = 100;
405  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
406  nbinsz = 50; minz = 0; maxz = 1.;
407 
408  histname = TString::Format("%s/JetHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
409  if (fForceBeamType != kpp) {
410  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
411  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
412  }
413  else {
414  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
415  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
416  }
417 
418  if (fPlotDCal) {
419  histname = TString::Format("%s/JetHistograms/hZVsPtDCal", jets->GetArrayName().Data());
420  if (fForceBeamType != kpp) {
421  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
422  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
423  }
424  else {
425  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
426  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
427  }
428  }
429 
430  // (Centrality, pT, Nconst)
431  nbinsx = 20; minx = 0; maxx = 100;
432  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
433  nbinsz = 50; minz = 0; maxz = fMaxPt;
434 
435  histname = TString::Format("%s/JetHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
436  if (fForceBeamType != kpp) {
437  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
438  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
439  }
440  else {
441  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
442  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
443  }
444 
445  if (fPlotDCal) {
446  histname = TString::Format("%s/JetHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
447  if (fForceBeamType != kpp) {
448  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
449  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
450  }
451  else {
452  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
453  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
454  }
455  }
456 
457  // (Centrality, pT) for eta<0 and eta>0
458  if (fForceBeamType != kpp) {
459  nbinsx = 20; minx = 0; maxx = 100;
460  nbinsy = fMaxPt; miny = 0; maxy = fMaxPt;
461 
462  histname = TString::Format("%s/JetHistograms/hEtaPosVsPtEMCal", jets->GetArrayName().Data());
463  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c})";
464  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy);
465 
466  histname = TString::Format("%s/JetHistograms/hEtaNegVsPtEMCal", jets->GetArrayName().Data());
467  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c})";
468  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy);
469  }
470 
471  // (Centrality, jet pT, Enonlincorr - Ehadcorr)
472  if (fForceBeamType != kpp) {
473  nbinsx = 20; minx = 0; maxx = 100;
474  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
475  nbinsz = nPtBins; minz = 0; maxz = fMaxPt;
476 
477  histname = TString::Format("%s/JetHistograms/hDeltaEHadCorr", jets->GetArrayName().Data());
478  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#sum#it{E}_{nonlincorr} - #it{E}_{hadcorr}";
479  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
480  }
481 
482  // (Median patch energy, calo type, jet pT, centrality)
483  if (fDoTriggerSimulation) {
484  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
485  title = histname + ";#it{E}_{patch,med};type;#it{p}_{T}^{corr} (GeV/#it{c});Centrality (%)";
486  Int_t nbins5[4] = {100, 2, nPtBins, 50};
487  Double_t min5[4] = {0,-0.5, 0, 0};
488  Double_t max5[4] = {50,1.5, fMaxPt, 100};
489  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins5, min5, max5);
490  }
491 
492  }
493 
494  // MB downscale factor histogram
495  if (fComputeMBDownscaling) {
496  histname = "Trigger/hMBDownscaleFactor";
497  title = histname + ";Downscale factor;counts";
498  fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 200);
499  }
500 
501 }
502 
503 /*
504  * This function allocates the histograms for the calorimeter performance study.
505  */
507 {
508  TString histname;
509  TString htitle;
510 
511  const Int_t nRcorrBins = 50;
512  Double_t *RcorrBins = GenerateFixedBinArray(nRcorrBins, 0., 1.);
513  const Int_t nCellBins = 30;
514  Double_t *cellBins = GenerateFixedBinArray(nCellBins, -0.5, 29.5);
515  const Int_t nMatchedTrackBins = 5;
516  Double_t *matchedTrackBins = GenerateFixedBinArray(nMatchedTrackBins, -0.5, 4.5);
517  const Int_t nDeltaEtaBins = 60;
518  Double_t *deltaEtaBins = GenerateFixedBinArray(nDeltaEtaBins, -0.015, 0.015);
519 
522 
523  // Plot M02 distribution (centrality, Eclus nonlincorr, M02)
524  histname = "ClusterHistograms/hM02";
525  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); M02";
527 
528  // Plot Ncell distribution for M02 > 0.4 and 0.1 < M02 < 0.4 (centrality, Eclus nonlincorr, Ncells)
529  histname = "ClusterHistograms/hNcellsM02G04";
530  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); Ncells";
531  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nCellBins, cellBins);
532 
533  histname = "ClusterHistograms/hNcellsM02L04";
534  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); Ncells";
535  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nCellBins, cellBins);
536 
539 
540  // Plot matched track pT for all clusters, M02 > 0.4 clusters, and 0.1 < M02 < 0.4 clusters (centrality, Eclus nonlincorr, trackPsum)
541  histname = "ClusterHistograms/hMatchedTrackPt";
542  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); #Sigma#it{p}_{track} (GeV/c)";
544 
545  histname = "ClusterHistograms/hMatchedTrackPtM02G04";
546  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); #Sigma#it{p}_{track} (GeV/c)";
548 
549  histname = "ClusterHistograms/hMatchedTrackPtM02L04";
550  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); #Sigma#it{p}_{track} (GeV/c)";
552 
553  // Plot number of matched tracks for all clusters, M02 > 0.4 clusters, and 0.1 < M02 < 0.4 clusters (centrality, Eclus nonlincorr, N matches)
554  histname = "ClusterHistograms/hMatchedTrackN";
555  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); N_{tracks}";
556  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nMatchedTrackBins, matchedTrackBins);
557 
558  histname = "ClusterHistograms/hMatchedTrackNM02G04";
559  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); N_{tracks}";
560  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nMatchedTrackBins, matchedTrackBins);
561 
562  histname = "ClusterHistograms/hMatchedTrackNM02L04";
563  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); N_{tracks}";
564  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nMatchedTrackBins, matchedTrackBins);
565 
566  // Plot M02 distribution for clusters with matched tracks (centrality, Eclus nonlincorr, M02)
567  histname = "ClusterHistograms/hM02Matched";
568  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); M02";
570 
571  // Plot M02 distribution for clusters without matched tracks (centrality, Eclus nonlincorr, M02)
572  histname = "ClusterHistograms/hM02Unmatched";
573  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); M02";
575 
576  // Plot clus-track deltaEta of matched tracks (deltaEta, Eclus, M02)
577  histname = "ClusterHistograms/hDeltaEtaCentral";
578  htitle = histname + ";#eta_{track} - #eta_{clus};#it{E}_{clus} (GeV); M02";
579  fHistManager.CreateTH3(histname.Data(), htitle.Data(), nDeltaEtaBins, deltaEtaBins, fNPtHistBins, fPtHistBins, fNM02HistBins, fM02HistBins);
580 
581  histname = "ClusterHistograms/hDeltaEtaPeripheral";
582  htitle = histname + ";#eta_{track} - #eta_{clus};#it{E}_{clus} (GeV); M02";
583  fHistManager.CreateTH3(histname.Data(), htitle.Data(), nDeltaEtaBins, deltaEtaBins, fNPtHistBins, fPtHistBins, fNM02HistBins, fM02HistBins);
584 
587 
588  // Plot E/p vs. M02 for 0-10% and 50-90% (Eclus nonlincorr, Eclus nonlincorr / trackPsum, M02)
589  histname = "ClusterHistograms/hEoverPM02Central";
590  htitle = histname + ";#it{E}_{clus} (GeV); #it{E}_{clus} / #Sigma#it{p}_{track} (GeV); M02";
592 
593  histname = "ClusterHistograms/hEoverPM02Peripheral";
594  htitle = histname + ";#it{E}_{clus} (GeV); #it{E}_{clus} / #Sigma#it{p}_{track} (GeV); M02";
596 
599 
600  // Plot Rcorr distribution (centrality, trackPSum, Rcorr = (Enonlincorr - Ehadcorr) / trackPSum)
601  histname = "ClusterHistograms/hRcorrVsCent";
602  htitle = histname + ";Centrality (%);#Sigma#it{p}_{track} (GeV); R_{corr} = #frac{#DeltaE_{clus}}{#Sigmap_{track}}";
603  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
604 
605  // Plot Rcorr distribution for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr)
606  histname = "ClusterHistograms/hRcorr0-10";
607  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); R_{corr} = #frac{#DeltaE_{clus}}{#Sigmap_{track}}";
608  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
609 
610  // Plot Rcorr distribution for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr)
611  histname = "ClusterHistograms/hRcorr50-90";
612  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); R_{corr} = #frac{#DeltaE_{clus}}{#Sigmap_{track}}";
613  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
614 
615  // Plot also Rcorr-clus (centrality, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
616  histname = "ClusterHistograms/hRcorrClusVsCent";
617  htitle = histname + ";Centrality (%);#Sigma#it{p}_{track} (GeV); #frac{#DeltaE_{clus}}{E_{clus}}";
618  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
619 
620  // Rcorr-clus for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
621  histname = "ClusterHistograms/hRcorrClus0-10";
622  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); #frac{#DeltaE_{clus}}{E_{clus}}";
623  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
624 
625  // Rcorr-clus for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
626  histname = "ClusterHistograms/hRcorrClus50-90";
627  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); #frac{#DeltaE_{clus}}{E_{clus}}";
628  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
629 
630  // Plot total track multiplicity
631  histname = "ClusterHistograms/hTrackMultiplicity";
632  htitle = histname + ";N_{tracks};Centrality (%)";
633  fHistManager.CreateTH2(histname.Data(), htitle.Data(), 1000, 0, 10000, 20, 0, 100);
634 }
635 
636 /*
637  * This function allocates the histograms for the jet composition study.
638  */
640 {
641  TString histname;
642  TString htitle;
643  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
644 
645  const Int_t nRejBins = 32;
646  Double_t* rejReasonBins = new Double_t[nRejBins+1];
647  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
648  const Int_t nContributorTypes = 11;
649  Double_t *contributorTypeBins = GenerateFixedBinArray(nContributorTypes, -0.5, 10.5);
650  const Int_t nParticleTypes = 17;
651  Double_t *particleTypeBins = GenerateFixedBinArray(nParticleTypes, -0.5, 16.5);
652 
653  AliEmcalContainer* cont = 0;
654  TIter nextClusColl(&fClusterCollArray);
655  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
656 
657  histname = "ClusterHistogramsMC/hClusterRejectionReasonMC";
658  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
659  TH2* histMC2 = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
660  SetRejectionReasonLabels(histMC2->GetXaxis());
661  }
662 
663  // M02 vs. Energy vs. Particle type
664  histname = "ClusterHistogramsMC/hM02VsParticleTypeCentral";
665  htitle = histname + ";M02;#it{E}_{clus} (GeV); Particle type";
666  TH3* hM02VsParticleTypeCentral = fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins, nParticleTypes, particleTypeBins);
667  SetParticleTypeLabels(hM02VsParticleTypeCentral->GetZaxis());
668 
669  histname = "ClusterHistogramsMC/hM02VsParticleTypePeripheral";
670  htitle = histname + ";M02;#it{E}_{clus} (GeV); Particle type";
671  TH3* hM02VsParticleTypePeripheral = fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins, nParticleTypes, particleTypeBins);
672  SetParticleTypeLabels(hM02VsParticleTypePeripheral->GetZaxis());
673 
674  // Plot photon energy in photon-hadron overlap clusters (Centrality, Photon energy, M02)
675  histname = "ClusterHistogramsMC/hPhotonHadronPhotonEnergy";
676  htitle = histname + ";Centrality (%);M02;#it{E}_{photon} (GeV)";
678 
679  // Plot hadron energy in hadron-photon overlap clusters (Centrality, Photon energy, M02)
680  histname = "ClusterHistogramsMC/hHadronPhotonHadronEnergy";
681  htitle = histname + ";Centrality (%);M02;#it{E}_{hadron} (GeV)";
683 
684  if (fPlotJetHistograms) {
685 
686  // M02 vs. Energy vs. Particle type vs. Jet pT, for particles inside jets
687  Int_t dim = 0;
688  TString title[20];
689  Int_t nbins[20] = {0};
690  Double_t min[30] = {0.};
691  Double_t max[30] = {0.};
692  Double_t *binEdges[20] = {0};
693 
694  title[dim] = "M02";
695  nbins[dim] = fNM02HistBins;
696  binEdges[dim] = fM02HistBins;
697  min[dim] = fM02HistBins[0];
698  max[dim] = fM02HistBins[fNM02HistBins];
699  dim++;
700 
701  title[dim] = "#it{E}_{clus} (GeV)";
702  nbins[dim] = fNPtHistBins;
703  binEdges[dim] = fPtHistBins;
704  min[dim] = fPtHistBins[0];
705  max[dim] = fPtHistBins[fNPtHistBins];
706  dim++;
707 
708  title[dim] = "Contributor type";
709  nbins[dim] = nContributorTypes;
710  min[dim] = -0.5;
711  max[dim] = 7.5;
712  binEdges[dim] = contributorTypeBins;
713  dim++;
714 
715  title[dim] = "#it{p}_{T,jet}^{corr}";
716  nbins[dim] = nPtBins;
717  min[dim] = 0;
718  max[dim] = fMaxPt;
719  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
720  dim++;
721 
722  TString thnname = "JetPerformanceMC/hM02VsContributorTypeJets";
723  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
724  for (Int_t i = 0; i < dim; i++) {
725  hn->GetAxis(i)->SetTitle(title[i]);
726  hn->SetBinEdges(i, binEdges[i]);
727  }
728 
729  // Particle composition inside each jet -- jet pT vs. particle type vs. particle number vs. particle pT sum
730  // (One entry per jet for each particle type)
731  dim = 0;
732 
733  title[dim] = "#it{p}_{T,jet}^{corr}";
734  nbins[dim] = nPtBins;
735  min[dim] = 0;
736  max[dim] = fMaxPt;
737  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
738  dim++;
739 
740  title[dim] = "Contributor type";
741  nbins[dim] = nContributorTypes;
742  min[dim] = -0.5;
743  max[dim] = 7.5;
744  binEdges[dim] = contributorTypeBins;
745  dim++;
746 
747  title[dim] = "N";
748  nbins[dim] = 30;
749  min[dim] = -0.5;
750  max[dim] = 29.5;
751  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
752  dim++;
753 
754  title[dim] = "#it{p}_{T,sum} (GeV)";
755  nbins[dim] = fNPtHistBins;
756  binEdges[dim] = fPtHistBins;
757  min[dim] = fPtHistBins[0];
758  max[dim] = fPtHistBins[fNPtHistBins];
759  dim++;
760 
761  thnname = "JetPerformanceMC/hJetComposition";
762  THnSparse* thn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
763  for (Int_t i = 0; i < dim; i++) {
764  thn->GetAxis(i)->SetTitle(title[i]);
765  thn->SetBinEdges(i, binEdges[i]);
766  }
767 
768  // Hadronic calo energy in each jet
769 
770  // Jet pT vs. Summed energy of hadronic clusters without a matched track
771  histname = "JetPerformance/hHadCaloEnergyUnmatched";
772  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
773  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
774 
775  // Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction)
776  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
777  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
778  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
779 
780  // Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction)
781  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
782  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
783  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
784 
785  }
786 
787 }
788 
789 /*
790  * This function sets axis labels for particle type histograms.
791  */
793 {
794  axis->SetBinLabel(1, "SinglePhoton");
795  axis->SetBinLabel(2, "SingleElectron");
796  axis->SetBinLabel(3, "SingleChargedPion");
797  axis->SetBinLabel(4, "SingleProton");
798  axis->SetBinLabel(5, "SingleAntiProton");
799  axis->SetBinLabel(6, "SingleChargedKaon");
800  axis->SetBinLabel(7, "SingleK0L");
801  axis->SetBinLabel(8, "SingleNeutron");
802  axis->SetBinLabel(9, "SingleAntiNeutron");
803  axis->SetBinLabel(10, "SingleOther");
804  axis->SetBinLabel(11, "PhotonHadron");
805  axis->SetBinLabel(12, "HadronPhoton");
806  axis->SetBinLabel(13, "MergedPi0");
807  axis->SetBinLabel(14, "PhotonPhotonOther");
808  axis->SetBinLabel(15, "HadronHadron");
809  axis->SetBinLabel(16, "TwoContributorsOther");
810  axis->SetBinLabel(17, "MoreThanTwoContributors");
811 }
812 
813 /*
814  * This function allocates background subtraction histograms, if enabled.
815  * A set of histograms is allocated per each jet container.
816  */
818 {
819  TString histname;
820  TString title;
821 
822  AliJetContainer* jets = 0;
823  TIter nextJetColl(&fJetCollArray);
824  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
825 
826  // EMCal
827  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
828  title = histname + ";Centrality;Scale factor;counts";
829  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
830 
831  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
832  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
833  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
834 
835  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCalFid", jets->GetArrayName().Data());
836  title = histname + ";Centrality;Scale factor;counts";
837  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
838 
839  // DCal
840  if (fPlotDCal) {
841  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCal", jets->GetArrayName().Data());
842  title = histname + ";Centrality;Scale factor;counts";
843  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
844 
845  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtDCal", jets->GetArrayName().Data());
846  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
847  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
848 
849  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCalFid", jets->GetArrayName().Data());
850  title = histname + ";Centrality;Scale factor;counts";
851  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
852  }
853 
854  }
855 }
856 
857 /*
858  * This function allocates the histograms for single jets, when the "simulated" trigger has been fired.
859  * A set of histograms is allocated per each jet container.
860  */
862 {
863  TString histname;
864  TString title;
865  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
866 
867  //----------------------------------------------
868  // Trigger patch histograms
869 
870  // patch eta vs. phi
871  histname = "TriggerSimHistograms/hEtaVsPhi";
872  title = histname + ";#eta_{patch} (rad);#phi_{patch} (rad)";
873  fHistManager.CreateTH2(histname.Data(), title.Data(), 140, -0.7, 0.7, 500, 1., 6.);
874 
875  // N patches
876  histname = "TriggerSimHistograms/hNPatches";
877  title = histname + ";#it{N}_{patches};type";
878  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, 0, 200, 2, -0.5, 1.5);
879 
880  // patch E vs. centrality
881  histname = "TriggerSimHistograms/hPatchE";
882  title = histname + ";Centrality (%);#it{E}_{patch} (GeV)";
883  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, nPtBins, 0, fMaxPt);
884 
885  // patch median vs. Centrality
886  histname = "TriggerSimHistograms/hPatchMedianE";
887  title = histname + ";Centrality (%);#it{E}_{patch,med} (GeV);type";
888  fHistManager.CreateTH3(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 50, 2, -0.5, 1.5);
889 
890  //----------------------------------------------
891  // Jet histograms for "triggered" events
892  AliJetContainer* jets = 0;
893  TIter nextJetColl(&fJetCollArray);
894  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
895 
896  // Jet rejection reason
897  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
898  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
899  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
900  SetRejectionReasonLabels(hist->GetXaxis());
901 
902  // Rho vs. Centrality
903  if (!jets->GetRhoName().IsNull()) {
904  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
905  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
906  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
907  }
908 
909  // (Centrality, pT, NEF)
910  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
911  Int_t nbinsy = nPtBins; Int_t miny = 0; Int_t maxy = fMaxPt;
912  Int_t nbinsz = 50; Int_t minz = 0; Int_t maxz = 1.;
913 
914  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
915  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
916  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
917 
918  if (fPlotDCal) {
919  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
920  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
921  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
922  }
923 
924  // pT-leading vs. pT
925  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
926  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
927  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
928 
929  // A vs. pT
930  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
931  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
932  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 50, 0, 0.5);
933 
934  // (Centrality, pT, z-leading (charged))
935  nbinsx = 20; minx = 0; maxx = 100;
936  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
937  nbinsz = 50; minz = 0; maxz = 1.;
938 
939  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
940  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
941  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
942 
943  if (fPlotDCal) {
944  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
945  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
946  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
947  }
948 
949  // z (charged) vs. pT
950  nbinsx = 20; minx = 0; maxx = 100;
951  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
952  nbinsz = 50; minz = 0; maxz = 1.;
953 
954  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
955  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
956  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
957 
958  if (fPlotDCal) {
959  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtDCal", jets->GetArrayName().Data());
960  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
961  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
962  }
963 
964  // (Centrality, pT, Nconst)
965  nbinsx = 20; minx = 0; maxx = 100;
966  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
967  nbinsz = 50; minz = 0; maxz = fMaxPt;
968 
969  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
970  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
971  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
972 
973  if (fPlotDCal) {
974  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
975  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
976  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
977  }
978 
979  }
980 }
981 
982 /*
983  * This function allocates histograms for matched truth-det jets in the case of embedding.
984  * The jet matching information must be previously filled by another task, such as AliJetResponseMaker.
985  */
987 {
988  TString histname;
989  TString title;
990  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
991 
992  // Response matrix, (centrality, pT-truth, pT-det)
993  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
994  Int_t nbinsy = fMaxPt; Int_t miny = 0; Int_t maxy = fMaxPt;
995  Int_t nbinsz = fMaxPt; Int_t minz = 0; Int_t maxz = fMaxPt;
996 
997  histname = "MatchedJetHistograms/hResponseMatrixEMCal";
998  if (fForceBeamType != kpp) {
999  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
1000  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1001  }
1002  else {
1003  title = histname + ";#it{p}_{T,corr}^{det} (GeV/#it{c});#it{p}_{T}^{truth} (GeV/#it{c})";
1004  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1005  }
1006 
1007  if (fPlotDCal) {
1008  histname = "MatchedJetHistograms/hResponseMatrixDCal";
1009  if (fForceBeamType != kpp) {
1010  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
1011  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1012  }
1013  else {
1014  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
1015  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1016  }
1017  }
1018 
1019  // JES shift, (centrality, pT-truth, (pT-det - pT-truth) / pT-truth)
1020  nbinsx = 20; minx = 0; maxx = 100;
1021  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
1022  nbinsz = 250; minz = -5.; maxz = 5.;
1023 
1024  histname = "MatchedJetHistograms/hJESshiftEMCal";
1025  if (fForceBeamType != kpp) {
1026  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1027  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1028  }
1029  else {
1030  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1031  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1032  }
1033 
1034  if (fPlotDCal) {
1035  histname = "MatchedJetHistograms/hJESshiftDCal";
1036  if (fForceBeamType != kpp) {
1037  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1038  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1039  }
1040  else {
1041  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1042  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1043  }
1044  }
1045 
1046  // NEF of det-level matched jets, (centrality, pT-truth, NEF)
1047  nbinsx = 20; minx = 0; maxx = 100;
1048  nbinsy = fMaxPt; miny = 0; maxy = fMaxPt;
1049  nbinsz = 50; minz = 0; maxz = 1.;
1050 
1051  histname = "MatchedJetHistograms/hNEFVsPt";
1052  if (fForceBeamType != kpp) {
1053  title = histname + ";Centrality (%);#it{p}_{T,corr}^{det} (GeV/#it{c});Calo energy fraction";
1054  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1055  }
1056  else {
1057  title = histname + ";#it{p}_{T,corr}^{det} (GeV/#it{c});Calo energy fraction";
1058  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1059  }
1060 
1061  // z-leading (charged) of det-level matched jets, (centrality, pT-truth, z-leading)
1062  nbinsx = 20; minx = 0; maxx = 100;
1063  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
1064  nbinsz = 50; minz = 0; maxz = 1.;
1065 
1066  histname = "MatchedJetHistograms/hZLeadingVsPt";
1067  if (fForceBeamType != kpp) {
1068  title = histname + ";Centrality (%);#it{p}_{T,corr}^{det} (GeV/#it{c});#it{z}_{leading}";
1069  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1070  }
1071  else {
1072  title = histname + ";#it{p}_{T,corr}^{det} (GeV/#it{c});#it{z}_{leading}";
1073  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1074  }
1075 
1076  // Matching distance, (centrality, pT-truth, R)
1077  nbinsx = 20; minx = 0; maxx = 100;
1078  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
1079  nbinsz = 50; minz = 0; maxz = 1.;
1080 
1081  histname = "MatchedJetHistograms/hMatchingDistance";
1082  if (fForceBeamType != kpp) {
1083  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});R";
1084  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1085 
1086  histname = "MatchedJetHistograms/hMatchingDistancepp";
1087  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});R";
1088  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1089  }
1090  else {
1091  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});R";
1092  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1093  }
1094 
1095  // Jet matching QA (copied from AliAnalysisTaskEmcalJetHCorrelations.cxx)
1096  if (fForceBeamType != kpp) {
1097  histname = "MatchedJetHistograms/fHistJetMatchingQA";
1098  title = histname;
1099  std::vector<std::string> binLabels = {"noMatch", "matchedJet", "sharedMomentumFraction", "partLevelMatchedJet", "jetDistance", "passedAllCuts"};
1100  auto histMatchedJetCuts = fHistManager.CreateTH1(histname.Data(), title.Data(), binLabels.size(), 0, binLabels.size());
1101  // Set label names
1102  for (unsigned int i = 1; i <= binLabels.size(); i++) {
1103  histMatchedJetCuts->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
1104  }
1105  histMatchedJetCuts->GetYaxis()->SetTitle("Number of jets");
1106  }
1107 
1108 }
1109 
1115 {
1116  // Configure base class to set fTriggerPatchInfo to array of trigger patches, each event
1117  // (Need to call this before base class ExecOnce)
1118  if (fDoTriggerSimulation) {
1119  this->SetCaloTriggerPatchInfoName("EmcalTriggers");
1120  }
1121 
1123 
1124  fNeedEmcalGeom = kTRUE;
1125 
1126  // Check if trigger patches are loaded
1127  if (fDoTriggerSimulation) {
1128  if (fTriggerPatchInfo) {
1129  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1130  TClass cls(objname);
1131  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1132  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1133  GetName(), cls.GetName(), "EmcalTriggers"));
1134  fTriggerPatchInfo = 0;
1135  }
1136  }
1137  if (!fTriggerPatchInfo) {
1138  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1139  return;
1140  }
1141  }
1142 }
1143 
1148 
1149  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
1150 
1152 
1153  // Get instance of the downscale factor helper class
1155  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
1156 
1157  // There are two possible min bias triggers for LHC15o
1158  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
1159  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
1160  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
1161 
1162  // Get the downscale factor for whichever MB trigger exists in the given run
1163  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
1164  Double_t downscalefactor;
1165  for (auto i : runtriggers) {
1166  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
1167  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
1168  break;
1169  }
1170  }
1171 
1172  // Store the inverse of the downscale factor, used later to weight the pT spectrum
1173  fMBUpscaleFactor = 1/downscalefactor;
1174 
1175  TString histname = "Trigger/hMBDownscaleFactor";
1176  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
1177 
1178  }
1179 
1180 }
1181 
1186 {
1187  if (fUseAliEventCuts) {
1188  if (!fEventCuts.AcceptEvent(InputEvent()))
1189  {
1190  PostData(1, fOutput);
1191  return kFALSE;
1192  }
1193  }
1194  else {
1196  }
1197  return kTRUE;
1198 }
1199 
1208 {
1209  TString histname;
1210  AliJetContainer* jetCont = 0;
1211  TIter next(&fJetCollArray);
1212  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1213  TString jetContName = jetCont->GetName();
1214 
1215  // Do a simple trigger simulation (if requested)
1216  if (fDoTriggerSimulation) {
1218  }
1219 
1220  }
1221 
1222  // Compute the full jet background scale factor and delta-pt
1223  if (fComputeBackground) {
1225  }
1226 
1227  // Only fill the embedding qa plots if:
1228  // - We are using the embedding helper
1229  // - The class has been initialized
1230  // - Both jet collections are available
1231  if (fEmbeddingQA.IsInitialized()) {
1233  }
1234 
1235  return kTRUE;
1236 }
1237 
1242 {
1243  TString histname;
1244 
1245  // Check if trigger patches are loaded
1246  if (fTriggerPatchInfo) {
1247  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1248  TClass cls(objname);
1249  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1250  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1251  GetName(), cls.GetName(), "EmcalTriggers"));
1252  fTriggerPatchInfo = 0;
1253  }
1254  }
1255  if (!fTriggerPatchInfo) {
1256  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1257  return;
1258  }
1259 
1260  // Compute patches in EMCal, DCal (I want offline simple trigger patch, i.e. patch calculated using FEE energy)
1261  std::vector<Double_t> vecEMCal;
1262  std::vector<Double_t> vecDCal;
1263  for(auto p : *fTriggerPatchInfo){
1264  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1265  if (recpatch) {
1266 
1267  if(!recpatch->IsJetHighSimple()) continue;
1268 
1269  histname = "TriggerSimHistograms/hEtaVsPhi";
1270  fHistManager.FillTH2(histname.Data(), recpatch->GetEtaGeo(), recpatch->GetPhiGeo());
1271 
1272  histname = "TriggerSimHistograms/hPatchE";
1273  fHistManager.FillTH2(histname.Data(), fCent, recpatch->GetPatchE());
1274 
1275  if (recpatch->IsEMCal()) {
1276  vecEMCal.push_back(recpatch->GetPatchE());
1277  } else {
1278  vecDCal.push_back(recpatch->GetPatchE());
1279  }
1280 
1281  }
1282  }
1283 
1284  // Compute the median in each calorimeter
1285  const Int_t nBkgPatchesEMCal = vecEMCal.size(); // 6*8;
1286  const Int_t nBkgPatchesDCal = vecDCal.size(); // 4*5;
1287  fMedianEMCal = TMath::Median(nBkgPatchesEMCal, &vecEMCal[0]); // point to array used internally by vector
1288  fMedianDCal = TMath::Median(nBkgPatchesDCal, &vecDCal[0]);
1289 
1290  histname = "TriggerSimHistograms/hPatchMedianE";
1291  fHistManager.FillTH3(histname.Data(), fCent, fMedianEMCal, kEMCal);
1292  fHistManager.FillTH3(histname.Data(), fCent, fMedianDCal, kDCal);
1293 
1294  histname = "TriggerSimHistograms/hNPatches";
1295  fHistManager.FillTH2(histname.Data(), nBkgPatchesEMCal, kEMCal);
1296  fHistManager.FillTH2(histname.Data(), nBkgPatchesDCal, kDCal);
1297 
1298  // Then compute background subtracted patches, by subtracting from each patch the median patch E from the opposite hemisphere
1299  // If a patch is above threshold, the event is "triggered"
1300  Bool_t fkEMCEJE = kFALSE;
1301  Double_t threshold = 20;
1302  for(auto p : *fTriggerPatchInfo){
1303  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1304  if (recpatch) {
1305 
1306  if(!recpatch->IsJetHighSimple()) continue;
1307 
1308  if (recpatch->IsEMCal()) {
1309  if ((recpatch->GetPatchE() - fMedianDCal) > threshold) {
1310  fkEMCEJE = kTRUE;
1311  break;
1312  }
1313  } else {
1314  if ((recpatch->GetPatchE() - fMedianEMCal) > threshold) {
1315  fkEMCEJE = kTRUE;
1316  break;
1317  }
1318  }
1319  }
1320  }
1321 
1322  if (fkEMCEJE) {
1324  }
1325 
1326 }
1327 
1335 {
1336 
1337  if (fPlotJetHistograms) {
1339  }
1340  if (fPlotClusterHistograms) {
1342  }
1345  }
1348  }
1349 
1350  return kTRUE;
1351 }
1352 
1358 {
1359  TString histname;
1360  AliJetContainer* jets = 0;
1361  TIter nextJetColl(&fJetCollArray);
1362  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1363  TString jetContName = jets->GetName();
1364 
1365  Double_t rhoVal = 0;
1366  if (jets->GetRhoParameter()) {
1367  rhoVal = jets->GetRhoVal();
1368  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
1369  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1370  }
1371 
1372  for (auto jet : jets->all()) {
1373 
1374  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1375  Float_t corrPt = GetJetPt(jet, rhoVal);
1376 
1377  // compute jet acceptance type
1378  Double_t type = GetJetType(jet);
1379  if ( type != kEMCal ) {
1380  if ( type != kDCal || !fPlotDCal ) {
1381  continue;
1382  }
1383  }
1384 
1385  // (Centrality, Area, pT) (fill before area cut)
1386  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
1387  if (fForceBeamType != kpp) {
1388  fHistManager.FillTH3(histname.Data(), fCent, corrPt, jet->Area());
1389  }
1390  else {
1391  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
1392  }
1393 
1394  // (Centrality, pT-leading, pT) (before leading hadron cuts)
1395  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1396  if (fForceBeamType != kpp) {
1397  fHistManager.FillTH3(histname.Data(), fCent, corrPt, ptLeading);
1398  }
1399  else {
1400  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1401  }
1402 
1403  // (Centrality, pT, z-leading (charged)) (before leading hadron cuts)
1404  if (type == kEMCal) {
1405  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
1406  }
1407  else if (type == kDCal) {
1408  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
1409  }
1410  TLorentzVector leadPart;
1411  jets->GetLeadingHadronMomentum(leadPart, jet);
1412  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1413  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
1414  if (fForceBeamType != kpp) {
1415  fHistManager.FillTH3(histname, fCent, corrPt, z);
1416  }
1417  else {
1418  fHistManager.FillTH2(histname, corrPt, z);
1419  }
1420 
1421  // Rejection reason
1422  UInt_t rejectionReason = 0;
1423  if (!jets->AcceptJet(jet, rejectionReason)) {
1424  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1425  if (fForceBeamType != kpp) {
1426  fHistManager.FillTH3(histname.Data(), fCent, jets->GetRejectionReasonBitPosition(rejectionReason), corrPt);
1427  }
1428  else {
1429  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), corrPt);
1430  }
1431  continue;
1432  }
1433 
1434  // (Centrality, pT, NEF)
1435  if (type == kEMCal) {
1436  histname = TString::Format("%s/JetHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
1437  }
1438  else if (type == kDCal) {
1439  histname = TString::Format("%s/JetHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
1440  }
1441  if (fForceBeamType != kpp) {
1442  fHistManager.FillTH3(histname, fCent, corrPt, jet->NEF());
1443  }
1444  else {
1445  fHistManager.FillTH2(histname, corrPt, jet->NEF());
1446  }
1447 
1448  // (Centrality, pT upscaled, calo type)
1449  if (fComputeMBDownscaling) {
1450  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
1451  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type, fMBUpscaleFactor);
1452  }
1453 
1454  // (Centrality, pT, z (charged))
1455  if (type == kEMCal) {
1456  histname = TString::Format("%s/JetHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
1457  }
1458  else if (type == kDCal) {
1459  histname = TString::Format("%s/JetHistograms/hZVsPtDCal", jets->GetArrayName().Data());
1460  }
1461  const AliVTrack* track;
1462  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1463  track = static_cast<AliVTrack*>(jet->Track(i));
1464  z = track->Pt() / TMath::Abs(corrPt);
1465  if (fForceBeamType != kpp) {
1466  fHistManager.FillTH3(histname, fCent, corrPt, z);
1467  }
1468  else {
1469  fHistManager.FillTH2(histname, corrPt, z);
1470  }
1471  }
1472 
1473  // (Centrality, pT, Nconst)
1474  if (type == kEMCal) {
1475  histname = TString::Format("%s/JetHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
1476  }
1477  else if (type == kDCal) {
1478  histname = TString::Format("%s/JetHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
1479  }
1480  if (fForceBeamType != kpp) {
1481  fHistManager.FillTH3(histname, fCent, corrPt, 1.*jet->GetNumberOfConstituents());
1482  }
1483  else {
1484  fHistManager.FillTH2(histname, corrPt, 1.*jet->GetNumberOfConstituents());
1485  }
1486 
1487  // (Centrality, pT) for eta<0 and eta>0
1488  if (fForceBeamType != kpp) {
1489  if (type == kEMCal) {
1490  if (jet->Eta() > 0) {
1491  histname = TString::Format("%s/JetHistograms/hEtaPosVsPtEMCal", jets->GetArrayName().Data());
1492  fHistManager.FillTH2(histname, fCent, corrPt);
1493  }
1494  else if (jet->Eta() < 0) {
1495  histname = TString::Format("%s/JetHistograms/hEtaNegVsPtEMCal", jets->GetArrayName().Data());
1496  fHistManager.FillTH2(histname, fCent, corrPt);
1497  }
1498  }
1499  }
1500 
1501  // (Centrality, jet pT, Enonlincorr - Ehadcorr)
1502  if (fForceBeamType != kpp) {
1503  Double_t deltaEhadcorr = 0;
1504  const AliVCluster* clus = nullptr;
1505  Int_t nClusters = jet->GetNumberOfClusters();
1506  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1507  clus = jet->Cluster(iClus);
1508  deltaEhadcorr += (clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy());
1509  }
1510 
1511  histname = TString::Format("%s/JetHistograms/hDeltaEHadCorr", jets->GetArrayName().Data());
1512  fHistManager.FillTH3(histname, fCent, corrPt, deltaEhadcorr);
1513  }
1514 
1515 
1516  // (Median patch energy, calo type, jet pT, centrality)
1517  if (fDoTriggerSimulation) {
1518  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
1519  Double_t x[4] = {fMedianEMCal, kEMCal, corrPt, fCent};
1520  fHistManager.FillTHnSparse(histname, x);
1521  Double_t y[4] = {fMedianDCal, kDCal, corrPt, fCent};
1522  fHistManager.FillTHnSparse(histname, y);
1523  }
1524 
1525  } //jet loop
1526 
1527  }
1528 }
1529 
1530 /*
1531  * This function fills the histograms for the calorimeter performance study.
1532  */
1534 {
1535  TString histname;
1536 
1537  // Loop through clusters
1539  const AliVCluster* clus;
1540  for (auto it : clusters->accepted_momentum()) {
1541 
1542  clus = it.second;
1543  Double_t clusPhi = it.first.Phi_0_2pi();
1544  Double_t clusEta = it.first.Eta();
1545 
1546  // Include only EMCal/DCal clusters
1547  if (!clus->IsEMCAL()) {
1548  continue;
1549  }
1550 
1551  // Compute the sum of matched track momentum, and various track matching / hadronic corretion quantities
1552  Double_t trackPSum = 0;
1553  Int_t nTracksMatched = 0;
1554  const AliVTrack* track = nullptr;
1555  for (Int_t itrack=0; itrack < clus->GetNTracksMatched(); itrack++) {
1556  track = dynamic_cast<AliVTrack*>(clus->GetTrackMatched(itrack));
1557  if (!track) {
1558  continue;
1559  }
1560 
1561  Double_t trackPhi = TVector2::Phi_0_2pi(track->GetTrackPhiOnEMCal());
1562  Double_t trackEta = track->GetTrackEtaOnEMCal();
1563  Double_t deta = TMath::Abs(clusEta - trackEta);
1564  Double_t dphi = TMath::Abs(clusPhi - trackPhi);
1565 
1567  trackPSum += track->P();
1568  nTracksMatched++;
1569  }
1570  }
1571 
1572  Double_t EoverP = 0;
1573  if (trackPSum > 1e-3) {
1574  EoverP = clus->GetNonLinCorrEnergy() / trackPSum;
1575  }
1576 
1577  Double_t deltaE = clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy();
1578  Double_t Rcorr = 0;
1579  if (trackPSum > 1e-3) {
1580  Rcorr = deltaE / trackPSum;
1581  }
1582  Double_t RcorrClus = deltaE / clus->GetNonLinCorrEnergy();
1583 
1586 
1587  // Fill M02 distribution (centrality, Eclus nonlincorr, M02)
1588  histname = "ClusterHistograms/hM02";
1589  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1590 
1591  // Plot Ncell distribution for M02 > 0.4 or M02 < 0.4 (centrality, Eclus nonlincorr, Ncells)
1592  if (clus->GetM02() > 0.4) {
1593  histname = "ClusterHistograms/hNcellsM02G04";
1594  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetNCells());
1595  }
1596  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1597  histname = "ClusterHistograms/hNcellsM02L04";
1598  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetNCells());
1599  }
1600 
1603 
1604  // Plot matched track pT (centrality, Eclus nonlincorr, trackPsum)
1605  histname = "ClusterHistograms/hMatchedTrackPt";
1606  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1607 
1608  if (clus->GetM02() > 0.4) {
1609  histname = "ClusterHistograms/hMatchedTrackPtM02G04";
1610  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1611  }
1612  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1613  histname = "ClusterHistograms/hMatchedTrackPtM02L04";
1614  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1615  }
1616 
1617  // Plot number of matched tracks (centrality, Eclus nonlincorr, N matches)
1618  histname = "ClusterHistograms/hMatchedTrackN";
1619  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1620 
1621  if (clus->GetM02() > 0.4) {
1622  histname = "ClusterHistograms/hMatchedTrackNM02G04";
1623  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1624  }
1625  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1626  histname = "ClusterHistograms/hMatchedTrackNM02L04";
1627  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1628  }
1629 
1630  // Plot M02 distribution for clusters with matched tracks (centrality, Eclus nonlincorr, M02)
1631  histname = "ClusterHistograms/hM02Matched";
1632  if (nTracksMatched > 0) {
1633  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1634  }
1635 
1636  // Plot M02 distribution for clusters without matched tracks (centrality, Eclus nonlincorr, M02)
1637  histname = "ClusterHistograms/hM02Unmatched";
1638  if (nTracksMatched == 0) {
1639  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1640  }
1641 
1642  // Plot clus-track deltaEta if there is one matched track (deltaEta, Eclus, M02)
1643  if (nTracksMatched == 1) {
1644 
1645  const AliVTrack* track = dynamic_cast<AliVTrack*>(clus->GetTrackMatched(0));
1646  if (track) {
1647  Double_t trackEta = track->GetTrackEtaOnEMCal();
1648  Double_t deta = trackEta - clusEta;
1649 
1650  if (fCent > 0 && fCent < 10) {
1651  histname = "ClusterHistograms/hDeltaEtaCentral";
1652  fHistManager.FillTH3(histname.Data(), deta, clus->GetNonLinCorrEnergy(), clus->GetM02());
1653  }
1654 
1655  if (fCent > 50 && fCent < 90) {
1656  histname = "ClusterHistograms/hDeltaEtaPeripheral";
1657  fHistManager.FillTH3(histname.Data(), deta, clus->GetNonLinCorrEnergy(), clus->GetM02());
1658  }
1659  }
1660  }
1661 
1664 
1665  // Plot E/p vs. M02 for 0-10% and 50-90% (Eclus nonlincorr, Eclus nonlincorr / trackPsum, M02)
1666  if (fCent > 0 && fCent < 10) {
1667  histname = "ClusterHistograms/hEoverPM02Central";
1668  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), EoverP, clus->GetM02());
1669  }
1670  if (fCent > 50 && fCent < 90) {
1671  histname = "ClusterHistograms/hEoverPM02Peripheral";
1672  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), EoverP, clus->GetM02());
1673  }
1674 
1677 
1678  // Fill Rcorr distribution (centrality, trackPSum, Rcorr = (Enonlincorr - Ehadcorr) / trackPSum)
1679  histname = "ClusterHistograms/hRcorrVsCent";
1680  fHistManager.FillTH3(histname, fCent, trackPSum, Rcorr);
1681 
1682  // Fill Rcorr distribution for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr)
1683  if (fCent > 0 && fCent < 10) {
1684  histname = "ClusterHistograms/hRcorr0-10";
1685  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, Rcorr);
1686  }
1687 
1688  // Fill Rcorr distribution for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr)
1689  if (fCent > 50 && fCent < 90) {
1690  histname = "ClusterHistograms/hRcorr50-90";
1691  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, Rcorr);
1692  }
1693 
1694  // Fill also Rcorr-clus (centrality, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1695  histname = "ClusterHistograms/hRcorrClusVsCent";
1696  fHistManager.FillTH3(histname, fCent, trackPSum, RcorrClus);
1697 
1698  // Fill Rcorr-clus for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1699  if (fCent > 0 && fCent < 10) {
1700  histname = "ClusterHistograms/hRcorrClus0-10";
1701  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, RcorrClus);
1702  }
1703 
1704  // Fill Rcorr-clus for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1705  if (fCent > 50 && fCent < 90) {
1706  histname = "ClusterHistograms/hRcorrClus50-90";
1707  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, RcorrClus);
1708  }
1709 
1710  }
1711 
1712  // Fill total track multiplicity
1713  histname = "ClusterHistograms/hTrackMultiplicity";
1714  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1715  Int_t nTracks = trackCont->GetNAcceptedTracks();
1716  fHistManager.FillTH2(histname.Data(), nTracks, fCent);
1717 
1718 }
1719 
1720 /*
1721  * This function fills particle composition histograms for the calorimeter performance study in MC.
1722  */
1724 {
1725  // If MC, get the MC event
1726  const AliMCEvent* mcevent = nullptr;
1727  if (fGeneratorLevel) {
1728  mcevent = MCEvent();
1729  }
1730  else {
1731  return;
1732  }
1733 
1734  // Loop through clusters, and plot M02 for each particle type
1736 
1737  // Loop through jets, to fill various histograms
1738  if (fPlotJetHistograms) {
1740  }
1741 
1742 }
1743 
1744 /*
1745  * Loop through clusters, and plot M02 for each particle type
1746  */
1748 {
1749  TString histname;
1751  const AliVCluster* clus;
1752  std::vector<ContributorType> vecContributorTypes;
1753  std::vector<Int_t> vecContributorLabels;
1754  for (auto it : clusters->accepted_momentum()) {
1755 
1756  clus = it.second;
1757 
1758  // Include only EMCal/DCal clusters (reject PHOS clusters)
1759  if (!clus->IsEMCAL()) {
1760  continue;
1761  }
1762 
1763  // Loop through the cluster's contributors in order to classify its type
1764  ParticleType particleType = kNotDefined;
1765  ContributorType contributorType = kUndefined;
1766  const Int_t nLabels = clus->GetNLabels();
1767 
1768  // Create a vector to store the contributor types for PhysicalPrimary particles
1769  vecContributorTypes.clear();
1770  vecContributorLabels.clear();
1771  for (Int_t iLabel=0; iLabel<nLabels; iLabel++) {
1772 
1773  Int_t label = clus->GetLabels()[iLabel];
1774  if (TMath::Abs(label) > 0) { // if the particle has a truth-level match, the label is nonzero
1775  contributorType = GetContributorType(clus, mcevent, label);
1776  if (contributorType != kUndefined) {
1777  vecContributorTypes.push_back(contributorType);
1778  vecContributorLabels.push_back(label);
1779  }
1780  }
1781  }
1782 
1783  Int_t nLabelsPhysPrim = vecContributorTypes.size();
1784 
1785  if (nLabelsPhysPrim == 1) {
1786 
1787  contributorType = vecContributorTypes[0];
1788 
1789  if (contributorType == kPhoton) {
1790  particleType = kSinglePhoton;
1791  }
1792  else if (contributorType == kElectron) {
1793  particleType = kSingleElectron;
1794  }
1795  else if (contributorType == kChargedPion) {
1796  particleType = kSingleChargedPion;
1797  }
1798  else if (contributorType == kProton) {
1799  particleType = kSingleProton;
1800  }
1801  else if (contributorType == kAntiProton) {
1802  particleType = kSingleAntiProton;
1803  }
1804  else if (contributorType == kChargedKaon) {
1805  particleType = kSingleChargedKaon;
1806  }
1807  else if (contributorType == kK0L) {
1808  particleType = kSingleK0L;
1809  }
1810  else if (contributorType == kNeutron) {
1811  particleType = kSingleNeutron;
1812  }
1813  else if (contributorType == kAntiNeutron) {
1814  particleType = kSingleAntiNeutron;
1815  }
1816  else {
1817  particleType = kSingleOther;
1818  }
1819 
1820  }
1821  else if (nLabelsPhysPrim == 2) {
1822 
1823  // Get the contributor particle types
1824  ContributorType contributorType1 = vecContributorTypes[0];
1825  ContributorType contributorType2 = vecContributorTypes[1];
1826 
1827  // Get the fraction of cluster energy from each contributor
1828  //Double_t frac0 = clus->GetClusterMCEdepFraction(0);
1829  Double_t frac1 = clus->GetClusterMCEdepFraction(1);
1830 
1831  // Check whether the leading/subleading contributors are photons/hadrons
1832  Bool_t isHadron1 = IsHadron(contributorType1);
1833  Bool_t isHadron2 = IsHadron(contributorType2);
1834  Bool_t isPhoton1 = contributorType1 == kPhoton;
1835  Bool_t isPhoton2 = contributorType2 == kPhoton;
1836 
1837  if (isHadron1 && isHadron2) {
1838  particleType = kHadronHadron;
1839  }
1840  else if (isPhoton1 && isHadron2) {
1841  particleType = kPhotonHadron;
1842 
1843  // Plot cluster energy when subleading hadron is subtracted
1844  Double_t photonEnergy = clus->GetNonLinCorrEnergy() * (1 - frac1);
1845  histname = "ClusterHistogramsMC/hPhotonHadronPhotonEnergy";
1846  fHistManager.FillTH3(histname.Data(), fCent, clus->GetM02(), photonEnergy);
1847  }
1848  else if (isHadron1 && isPhoton2) {
1849  particleType = kHadronPhoton;
1850 
1851  // Plot cluster energy when subleading hadron is subtracted
1852  Double_t hadronEnergy = clus->GetNonLinCorrEnergy() * (1 - frac1);
1853  histname = "ClusterHistogramsMC/hHadronPhotonHadronEnergy";
1854  fHistManager.FillTH3(histname.Data(), fCent, clus->GetM02(), hadronEnergy);
1855  }
1856  else if (isPhoton1 && isPhoton2) {
1857 
1858  // By default, assume the two photons are not a merged pi0
1859  particleType = kPhotonPhotonOther;
1860 
1861  // Using the vector of accepted contributor labels, check whether the two photons are from the same pi0
1862  AliAODMCParticle *part1 = fGeneratorLevel->GetMCParticleWithLabel(vecContributorLabels[0]);
1863  AliAODMCParticle *part2 = fGeneratorLevel->GetMCParticleWithLabel(vecContributorLabels[1]);
1864  if (part1 && part2) {
1865  Int_t iMother1 = part1->GetMother();
1866  Int_t iMother2 = part2->GetMother();
1867  AliVParticle *mother1 = mcevent->GetTrack(iMother1);
1868  AliVParticle *mother2 = mcevent->GetTrack(iMother2);
1869 
1870  if (mother1 && mother2) {
1871  if ( (mother1->PdgCode() == 111) && (mother2->PdgCode() == 111) ) {
1872  if (iMother1 == iMother2) {
1873  particleType = kMergedPi0;
1874  }
1875  }
1876  }
1877  }
1878  }
1879  else {
1880  particleType = kTwoContributorsOther; // this includes partially contained conversion overlaps
1881  }
1882 
1883  }
1884  else if (nLabelsPhysPrim > 2) {
1885  particleType = kMoreThanTwoContributors;
1886  }
1887 
1888  // (M02, Eclus, part type)
1889  if (fCent > 0 && fCent < 10) {
1890  histname = "ClusterHistogramsMC/hM02VsParticleTypeCentral";
1891  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType);
1892  }
1893  if (fCent > 50 && fCent < 90) {
1894  histname = "ClusterHistogramsMC/hM02VsParticleTypePeripheral";
1895  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType);
1896  }
1897 
1898  }
1899 }
1900 
1901 /*
1902  * Loop through jets, to fill particle composition histograms.
1903  */
1905 {
1906  TString histname;
1907 
1908  const AliVCluster* clus;
1909  AliJetContainer* jets = GetJetContainer(0); // there is only a single, det-level jet finder here
1910  for (const auto jet : jets->accepted()) {
1911 
1912  Double_t jetPt = GetJetPt(jet, 0);
1913 
1914  // Keep track of hadronic calo energy in each jet
1915  Double_t hadCaloEnergyUnmatched = 0;
1916  Double_t hadCaloEnergyMatchedNonlincorr = 0;
1917  Double_t hadCaloEnergyMatchedHadCorr = 0;
1918 
1919  // Loop through clusters in each jet, to plot several histograms
1920  Int_t nClusters = jet->GetNumberOfClusters();
1921  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1922 
1923  clus = jet->Cluster(iClus);
1924 
1925  // Get the particle type of the cluster
1926  ContributorType contributorType = kUndefined;
1927  Int_t label = TMath::Abs(clus->GetLabel());
1928  if (label > 0) {
1929  contributorType = GetContributorType(clus, mcevent, label);
1930  }
1931 
1932  // Plot M02 for each particle type
1933  histname = "JetPerformanceMC/hM02VsContributorTypeJets";
1934  Double_t x[4] = {clus->GetM02(), clus->GetNonLinCorrEnergy(), contributorType, jetPt};
1935  fHistManager.FillTHnSparse(histname, x);
1936 
1937  // If the cluster is a hadron, sum its energy to compute the jet's hadronic calo energy
1938  Bool_t isHadron = IsHadron(contributorType);
1939  if (isHadron) {
1940  Bool_t hasMatchedTrack = (clus->GetNTracksMatched() > 0);
1941  //Bool_t hasMatchedTrack = ((clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy()) > 1e-3);
1942  if (hasMatchedTrack) {
1943  hadCaloEnergyMatchedNonlincorr += clus->GetNonLinCorrEnergy();
1944  hadCaloEnergyMatchedHadCorr += clus->GetHadCorrEnergy();
1945  }
1946  else {
1947  hadCaloEnergyUnmatched += clus->GetNonLinCorrEnergy();
1948  }
1949  }
1950 
1951  }
1952 
1953  // Fill hadronic calo energy in each jet
1954 
1955  // (Jet pT, Summed energy of hadronic clusters without a matched track)
1956  histname = "JetPerformance/hHadCaloEnergyUnmatched";
1957  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyUnmatched);
1958 
1959  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction))
1960  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
1961  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedNonlincorr);
1962 
1963  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction))
1964  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
1965  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedHadCorr);
1966 
1967  // Loop through particle types, and plot jet composition for each particle type
1968  histname = "JetPerformanceMC/hJetComposition";
1969  for (Int_t type = 0; type < 8; type++) {
1970 
1971  ContributorType contributorType = kUndefined;
1972  Double_t nSum = 0;
1973  Double_t pTsum = 0;
1974 
1975  // Loop through clusters in jet, and add to sum if cluster matches particle type
1976  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1977 
1978  clus = jet->Cluster(iClus);
1979 
1980  Int_t label = TMath::Abs(clus->GetLabel());
1981  if (label > 0) {
1982  contributorType = GetContributorType(clus, mcevent, label);
1983  }
1984 
1985  if (type == contributorType) {
1986  nSum++;
1987  pTsum += clus->GetNonLinCorrEnergy();
1988  }
1989  }
1990 
1991  // Fill jet composition histogram
1992  Double_t x[4] = {jetPt, 1.*type, nSum, pTsum};
1993  fHistManager.FillTHnSparse(histname, x);
1994 
1995  }
1996  }
1997 }
1998 
2003 {
2004  // Loop over tracks and clusters in order to:
2005  // (1) Compute scale factor for full jets
2006  // (2) Compute delta-pT for full jets, with the random cone method
2007 
2008  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
2009  Double_t etaTPC = 0.9;
2010  Double_t etaEMCal = 0.7;
2011  Double_t etaMinDCal = 0.22;
2012  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
2013  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
2014  Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
2015  Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
2016 
2017  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
2018  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
2019  Double_t accDCal = 2 * (etaEMCal - etaMinDCal) * (phiMaxDCal - phiMinDCal);
2020 
2021  // Loop over jet containers
2022  AliJetContainer* jetCont = 0;
2023  TIter nextJetColl(&fJetCollArray);
2024  while ((jetCont = static_cast<AliJetContainer*>(nextJetColl()))) {
2025 
2026  // Define fiducial acceptances, to be used to generate random cones, and for scale factor studies
2027  TRandom3* r = new TRandom3(0);
2028  Double_t jetR = jetCont->GetJetRadius();
2029  Double_t etaEMCalfid = etaEMCal - jetR;
2030  Double_t etaMinDCalfid = etaMinDCal + jetR;
2031  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
2032  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
2033  Double_t phiMinDCalfid = phiMinDCal + jetR;
2034  Double_t phiMaxDCalfid = phiMaxDCal - jetR;
2035  Double_t accEMCalfid = 2 * etaEMCalfid * (phiMaxEMCalfid - phiMinEMCalfid);
2036  Double_t accDCalfid = 2 * (etaEMCalfid - etaMinDCalfid) * (phiMaxDCalfid - phiMinDCalfid);
2037  if ( (etaEMCalfid - etaMinDCalfid) < 0) {
2038  accDCalfid = 0;
2039  }
2040 
2041  // Generate EMCal random cone eta-phi
2042  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
2043  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
2044 
2045  // Generate DCal random cone eta-phi
2046  Double_t etaDCalRC = r->Uniform(etaMinDCalfid, etaEMCalfid);
2047  Double_t sign = r->Uniform(-1., 1.);
2048  if (sign < 0) {
2049  etaDCalRC = -1*etaDCalRC;
2050  }
2051  Double_t phiDCalRC = r->Uniform(phiMinDCalfid, phiMaxDCalfid);
2052 
2053  // Initialize the various sums to 0
2054  Double_t trackPtSumTPC = 0;
2055  Double_t trackPtSumEMCal = 0;
2056  Double_t trackPtSumEMCalfid = 0;
2057  Double_t trackPtSumDCal = 0;
2058  Double_t trackPtSumDCalfid = 0;
2059  Double_t trackPtSumEMCalRC = 0;
2060  Double_t trackPtSumDCalRC = 0;
2061  Double_t clusPtSumEMCal = 0;
2062  Double_t clusPtSumEMCalfid = 0;
2063  Double_t clusPtSumDCal = 0;
2064  Double_t clusPtSumDCalfid = 0;
2065  Double_t clusPtSumEMCalRC = 0;
2066  Double_t clusPtSumDCalRC = 0;
2067 
2068  // Loop over tracks. Sum the track pT:
2069  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal fiducial volume, (4) in the DCal, (5) in the DCal fiducial volume, (6) in the EMCal random cone, (7) in the DCal random cone
2070  // Note: Loops over all det-level track containers. For data there should be only one. For embedding, there should be signal+background tracks.
2071  AliParticleContainer * partCont = 0;
2072  AliTLorentzVector track;
2073  Double_t trackEta;
2074  Double_t trackPhi;
2075  Double_t trackPt;
2076  Double_t deltaR;
2077  TIter nextPartCont(&fParticleCollArray);
2078  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
2079 
2080  TString partContName = partCont->GetName();
2081  if (!partContName.CompareTo("tracks")) {
2082 
2083  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(partCont);
2084  for (auto trackIterator : trackCont->accepted_momentum() ) {
2085 
2086  track.Clear();
2087  track = trackIterator.first;
2088  trackEta = track.Eta();
2089  trackPhi = track.Phi_0_2pi();
2090  trackPt = track.Pt();
2091 
2092  // (1)
2093  if (TMath::Abs(trackEta) < etaTPC) {
2094  trackPtSumTPC += trackPt;
2095  }
2096 
2097  // (2)
2098  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
2099  trackPtSumEMCal += trackPt;
2100  }
2101 
2102  // (3)
2103  if (TMath::Abs(trackEta) < etaEMCalfid && trackPhi > phiMinEMCalfid && trackPhi < phiMaxEMCalfid) {
2104  trackPtSumEMCalfid += trackPt;
2105  }
2106 
2107  // (4)
2108  if (fPlotDCal) {
2109  if (TMath::Abs(trackEta) > etaMinDCal && TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinDCal && trackPhi < phiMaxDCal) {
2110  trackPtSumDCal += trackPt;
2111  }
2112  }
2113 
2114  // (5)
2115  if (fPlotDCal) {
2116  if (TMath::Abs(trackEta) > etaMinDCalfid && TMath::Abs(trackEta) < etaEMCalfid && trackPhi > phiMinDCalfid && trackPhi < phiMaxDCalfid) {
2117  trackPtSumDCalfid += trackPt;
2118  }
2119  }
2120 
2121  // (6)
2122  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
2123  if (deltaR < jetR) {
2124  trackPtSumEMCalRC += trackPt;
2125  }
2126 
2127  // (7)
2128  if (fPlotDCal) {
2129  deltaR = GetDeltaR(&track, etaDCalRC, phiDCalRC);
2130  if (deltaR < jetR) {
2131  trackPtSumDCalRC += trackPt;
2132  }
2133  }
2134  }
2135  }
2136  }
2137 
2138  // Loop over clusters, if the jet container is for full jets. Sum the cluster ET:
2139  // (1) in the EMCal, (2) in the EMCal fiducial volume, (3) in the DCal, (4), in the DCal fiducial volume, (5) in the EMCal random cone, (6) in the DCal random cone
2140  if (jetCont->GetClusterContainer()) {
2142 
2143  AliTLorentzVector clus;
2144  Double_t clusEta;
2145  Double_t clusPhi;
2146  Double_t clusPt;
2147  for (auto clusIterator : clusCont->accepted_momentum() ) {
2148 
2149  clus.Clear();
2150  clus = clusIterator.first;
2151  clusEta = clus.Eta();
2152  clusPhi = clus.Phi_0_2pi();
2153  clusPt = clus.Pt();
2154 
2155  // (1)
2156  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
2157  clusPtSumEMCal += clusPt;
2158  }
2159 
2160  // (2)
2161  if (TMath::Abs(clusEta) < etaEMCalfid && clusPhi > phiMinEMCalfid && clusPhi < phiMaxEMCalfid) {
2162  clusPtSumEMCalfid += clusPt;
2163  }
2164 
2165  // (3)
2166  if (fPlotDCal) {
2167  if (TMath::Abs(clusEta) > etaMinDCal && TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinDCal && clusPhi < phiMaxDCal) {
2168  clusPtSumDCal += clusPt;
2169  }
2170  }
2171 
2172  // (4)
2173  if (fPlotDCal) {
2174  if (TMath::Abs(clusEta) > etaMinDCalfid && TMath::Abs(clusEta) < etaEMCalfid && clusPhi > phiMinDCalfid && clusPhi < phiMaxDCalfid) {
2175  clusPtSumDCalfid += clusPt;
2176  }
2177  }
2178 
2179  // (5)
2180  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
2181  if (deltaR < jetR) {
2182  clusPtSumEMCalRC += clusPt;
2183  }
2184 
2185  // (6)
2186  if (fPlotDCal) {
2187  deltaR = GetDeltaR(&clus, etaDCalRC, phiDCalRC);
2188  if (deltaR < jetR) {
2189  clusPtSumDCalRC += clusPt;
2190  }
2191  }
2192 
2193  }
2194  }
2195 
2196  // Compute the scale factor, as a function of centrality, for (1) EMCal, (2) EMCalfid, (3) DCal, (4) DCalfid
2197  // (1)
2198  Double_t numerator = (trackPtSumEMCal + clusPtSumEMCal) / accEMCal;
2199  Double_t denominator = trackPtSumTPC / accTPC;
2200  Double_t scaleFactor = numerator / denominator;
2201  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
2202  fHistManager.FillTH2(histname, fCent, scaleFactor);
2203 
2204  // (2)
2205  if (accEMCalfid > 1e-3) {
2206  numerator = (trackPtSumEMCalfid + clusPtSumEMCalfid) / accEMCalfid;
2207  scaleFactor = numerator / denominator;
2208  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCalFid", jetCont->GetArrayName().Data());
2209  fHistManager.FillTH2(histname, fCent, scaleFactor);
2210  }
2211 
2212  // (3)
2213  if (fPlotDCal) {
2214  numerator = (trackPtSumDCal + clusPtSumDCal) / accDCal;
2215  scaleFactor = numerator / denominator;
2216  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCal", jetCont->GetArrayName().Data());
2217  fHistManager.FillTH2(histname, fCent, scaleFactor);
2218  }
2219 
2220  // (4)
2221  if (fPlotDCal) {
2222  if (accDCalfid > 1e-3) {
2223  numerator = (trackPtSumDCalfid + clusPtSumDCalfid) / accDCalfid;
2224  scaleFactor = numerator / denominator;
2225  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCalFid", jetCont->GetArrayName().Data());
2226  fHistManager.FillTH2(histname, fCent, scaleFactor);
2227  }
2228  }
2229 
2230  // Compute delta pT, as a function of centrality
2231 
2232  // EMCal
2233  Double_t rho = jetCont->GetRhoVal();
2234  Double_t deltaPt = trackPtSumEMCalRC + clusPtSumEMCalRC - rho * TMath::Pi() * jetR * jetR;
2235  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
2236  fHistManager.FillTH2(histname, fCent, deltaPt);
2237 
2238  // DCal
2239  if (fPlotDCal) {
2240  if (accDCalfid > 1e-3) {
2241  deltaPt = trackPtSumDCalRC + clusPtSumDCalRC - rho * TMath::Pi() * jetR * jetR;
2242  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtDCal", jetCont->GetArrayName().Data());
2243  fHistManager.FillTH2(histname, fCent, deltaPt);
2244  }
2245  }
2246 
2247  delete r;
2248 
2249  }
2250 
2251 }
2252 
2258 {
2259  TString histname;
2260  AliJetContainer* jets = 0;
2261  TIter nextJetColl(&fJetCollArray);
2262  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
2263  TString jetContName = jets->GetName();
2264 
2265  Double_t rhoVal = 0;
2266  if (jets->GetRhoParameter()) {
2267  rhoVal = jets->GetRhoVal();
2268  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
2269  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
2270  }
2271 
2272  for (auto jet : jets->all()) {
2273 
2274  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
2275  Float_t corrPt = GetJetPt(jet, rhoVal);
2276 
2277  // A vs. pT (fill before area cut)
2278  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
2279  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
2280 
2281  // Rejection reason
2282  UInt_t rejectionReason = 0;
2283  if (!jets->AcceptJet(jet, rejectionReason)) {
2284  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
2285  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
2286  continue;
2287  }
2288 
2289  // compute jet acceptance type
2290  Double_t type = GetJetType(jet);
2291  if ( type != kEMCal ) {
2292  if ( type != kDCal || !fPlotDCal ) {
2293  continue;
2294  }
2295  }
2296 
2297  // (Centrality, pT, NEF, calo type)
2298  if (type == kEMCal) {
2299  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
2300  }
2301  else if (type == kDCal) {
2302  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
2303  }
2304  fHistManager.FillTH3(histname, fCent, corrPt, jet->NEF());
2305 
2306  // pT-leading vs. pT
2307  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
2308  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
2309 
2310  // (Centrality, pT, z-leading (charged), calo type)
2311  if (type == kEMCal) {
2312  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
2313  }
2314  else if (type == kDCal) {
2315  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
2316  }
2317  TLorentzVector leadPart;
2318  jets->GetLeadingHadronMomentum(leadPart, jet);
2319  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
2320  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
2321  fHistManager.FillTH3(histname, fCent, corrPt, z);
2322 
2323  // (Centrality, pT, z (charged), calo type)
2324  if (type == kEMCal) {
2325  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
2326  }
2327  else if (type == kDCal) {
2328  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtDCal", jets->GetArrayName().Data());
2329  }
2330  const AliVTrack* track;
2331  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
2332  track = static_cast<AliVTrack*>(jet->Track(i));
2333  z = track->Pt() / TMath::Abs(corrPt);
2334  fHistManager.FillTH3(histname, fCent, corrPt, z);
2335  }
2336 
2337  // (Centrality, pT, Nconst)
2338  if (type == kEMCal) {
2339  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
2340  }
2341  else if (type == kDCal) {
2342  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
2343  }
2344  fHistManager.FillTH3(histname, fCent, corrPt, 1.*jet->GetNumberOfConstituents());
2345 
2346  } //jet loop
2347 
2348  }
2349 }
2350 
2356 {
2357  TString histname;
2358  AliJetContainer* jets = 0;
2359  const AliEmcalJet* matchedPartLevelJet = nullptr;
2360  TIter nextJetColl(&fJetCollArray);
2361  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
2362  TString jetContName = jets->GetName();
2363 
2364  // Only loop over jets in the detector-level jet container
2365  if (jetContName.Contains("mcparticles")) {
2366  continue;
2367  }
2368 
2369  Double_t rhoVal = 0;
2370  if (jets->GetRhoParameter()) {
2371  rhoVal = jets->GetRhoVal();
2372  }
2373 
2374  for (auto jet : jets->accepted()) {
2375 
2376  if (fUseResponseMaker) {
2377  // Get the matched part-level jet, based on JetResponseMaker's geometrical criteria
2378  matchedPartLevelJet = jet->MatchedJet();
2379  }
2380  else {
2381  // Get the matched part-level jet, if one exists, subject to fMinSharedMomentumFraction, fMaxMatchedJetDistance criteria
2382  matchedPartLevelJet = GetMatchedPartLevelJet(jets, jet, "MatchedJetHistograms/fHistJetMatchingQA");
2383  }
2384 
2385  // Check that the matched jet exists, and is accepted
2386  if (!matchedPartLevelJet) {
2387  continue;
2388  }
2389  UInt_t rejectionReason = 0;
2390  if (!fMCJetContainer->AcceptJet(matchedPartLevelJet, rejectionReason)) {
2391  continue;
2392  }
2393 
2394  // compute jet acceptance type
2395  Double_t type = GetJetType(jet);
2396  if ( type != kEMCal ) {
2397  if ( type != kDCal || !fPlotDCal ) {
2398  continue;
2399  }
2400  }
2401 
2402  Float_t detPt = GetJetPt(jet, rhoVal);
2403  Float_t truthPt = matchedPartLevelJet->Pt();
2404 
2405  // Fill response matrix (centrality, pT-truth, pT-det)
2406  if (type == kEMCal) {
2407  histname = "MatchedJetHistograms/hResponseMatrixEMCal";
2408  }
2409  else if (type == kDCal) {
2410  histname = "MatchedJetHistograms/hResponseMatrixDCal";
2411  }
2412  if (fForceBeamType != kpp) {
2413  fHistManager.FillTH3(histname, fCent, truthPt, detPt);
2414  }
2415  else {
2416  fHistManager.FillTH2(histname, detPt, truthPt);
2417  }
2418 
2419  // Fill JES shift (centrality, pT-truth, (pT-det - pT-truth) / pT-truth)
2420  if (type == kEMCal) {
2421  histname = "MatchedJetHistograms/hJESshiftEMCal";
2422  }
2423  else if (type == kDCal) {
2424  histname = "MatchedJetHistograms/hJESshiftDCal";
2425  }
2426  if (fForceBeamType != kpp) {
2427  fHistManager.FillTH3(histname, fCent, truthPt, (detPt-truthPt)/truthPt );
2428  }
2429  else {
2430  fHistManager.FillTH2(histname, truthPt, (detPt-truthPt)/truthPt );
2431  }
2432 
2433  // Fill NEF of det-level matched jets (centrality, pT-truth, NEF)
2434  histname = "MatchedJetHistograms/hNEFVsPt";
2435  if (fForceBeamType != kpp) {
2436  fHistManager.FillTH3(histname, fCent, truthPt, jet->NEF());
2437  }
2438  else {
2439  fHistManager.FillTH2(histname, truthPt, jet->NEF());
2440  }
2441 
2442  // Fill z-leading (charged) of det-level matched jets (centrality, pT-truth, z-leading)
2443  histname = "MatchedJetHistograms/hZLeadingVsPt";
2444  TLorentzVector leadPart;
2445  jets->GetLeadingHadronMomentum(leadPart, jet);
2446  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
2447  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
2448  if (fForceBeamType != kpp) {
2449  fHistManager.FillTH3(histname, fCent, truthPt, z);
2450  }
2451  else {
2452  fHistManager.FillTH2(histname, truthPt, z);
2453  }
2454 
2455  // Fill matching distance between combined jet and pp det-level jet (centrality, pT-truth, R)
2456  histname = "MatchedJetHistograms/hMatchingDistance";
2457  if (fForceBeamType != kpp) {
2458  fHistManager.FillTH3(histname, fCent, truthPt, jet->ClosestJetDistance());
2459  }
2460  else {
2461  fHistManager.FillTH2(histname, truthPt, jet->ClosestJetDistance());
2462  }
2463 
2464  // Fill matching distance between pp det-level jet and pp truth-level jet (centrality, pT-truth, R)
2465  if (fForceBeamType != kpp) {
2466  histname = "MatchedJetHistograms/hMatchingDistancepp";
2467  fHistManager.FillTH3(histname, fCent, truthPt, matchedPartLevelJet->ClosestJetDistance());
2468  }
2469 
2470  } //jet loop
2471  }
2472 }
2473 
2474 /*
2475  * Compute the MC particle type of a given cluster contributor, using the MC particle container
2476  */
2478 {
2479  ContributorType contributorType = kUndefined;
2480 
2481  AliAODMCParticle *part = fGeneratorLevel->GetMCParticleWithLabel(label);
2482  if (part) {
2483 
2484  TString histname = "ClusterHistogramsMC/hClusterRejectionReasonMC";
2485  UInt_t rejectionReason = 0;
2486  if (!fGeneratorLevel->AcceptMCParticle(part, rejectionReason)) {
2487  fHistManager.FillTH2(histname, fGeneratorLevel->GetRejectionReasonBitPosition(rejectionReason), clus->GetNonLinCorrEnergy());
2488  return contributorType;
2489  }
2490 
2491  if (part->GetGeneratorIndex() == 0) { // generator index in cocktail
2492 
2493  // select charged pions, protons, kaons, electrons, muons
2494  Int_t pdg = part->PdgCode();
2495 
2496  if (pdg == 22) { // gamma 22
2497  contributorType = kPhoton;
2498  }
2499  else if (TMath::Abs(pdg) == 211) { // pi+ 211 (abs value ensures both particles and antiparticles are included)
2500  contributorType = kChargedPion;
2501  }
2502  else if (pdg == 2212) { // proton 2212
2503  contributorType = kProton;
2504  }
2505  else if (pdg == -2212) {
2506  contributorType = kAntiProton;
2507  }
2508  else if (TMath::Abs(pdg) == 321) { // K+ 321
2509  contributorType = kChargedKaon;
2510  }
2511  else if (pdg == 130) { // K0L 130
2512  contributorType = kK0L;
2513  }
2514  else if (pdg == 2112) { // neutron 2112
2515  contributorType = kNeutron;
2516  }
2517  else if (pdg == -2112) {
2518  contributorType = kAntiNeutron;
2519  }
2520  else if (TMath::Abs(pdg) == 11) { // e- 11
2521  contributorType = kElectron;
2522  }
2523  else if (TMath::Abs(pdg) == 13) { // mu- 13
2524  contributorType = kMuon;
2525  }
2526  else {
2527  contributorType = kOther;
2528  }
2529  }
2530  }
2531  return contributorType;
2532 }
2533 
2538 {
2539  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
2540  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
2541  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
2542  return deltaR;
2543 }
2544 
2549 {
2550  UInt_t jetType = jet->GetJetAcceptanceType();
2551  Double_t type = -1;
2552  if (jetType & AliEmcalJet::kEMCALfid) {
2553  type = kEMCal;
2554  }
2555  else if (jetType & AliEmcalJet::kDCALonlyfid) {
2556  type = kDCal;
2557  }
2558 
2559  return type;
2560 }
2561 
2566 {
2567  Double_t pT = jet->Pt() - rho * jet->Area();
2568  return pT;
2569 }
2570 
2575 {
2576  return (contributor == kChargedPion) || (contributor == kProton) || (contributor == kAntiProton) || (contributor == kChargedKaon) || (contributor == kK0L) || (contributor == kNeutron) || (contributor == kAntiNeutron);
2577 }
2578 
2598 const AliEmcalJet* AliAnalysisTaskEmcalJetPerformance::GetMatchedPartLevelJet(const AliJetContainer * jets, const AliEmcalJet * jet, const std::string & histName)
2599 {
2600  bool returnValue = false;
2601  const AliEmcalJet* partLevelJet = nullptr;
2602 
2603  // First, check if combined jet has a pp det-level match assigned
2604  if (jet->ClosestJet()) {
2605  fHistManager.FillTH1(histName.c_str(), "matchedJet");
2606  returnValue = true;
2607  AliDebugStream(4) << "Jet is matched!\nJet: " << jet->toString() << "\n";
2608 
2609  // Check shared momentum fraction
2610  // We explicitly want to use indices instead of geometric matching
2611  double sharedFraction = jets->GetFractionSharedPt(jet, nullptr);
2612  if (sharedFraction < fMinSharedMomentumFraction) {
2613  AliDebugStream(4) << "Jet rejected due to shared momentum fraction of " << sharedFraction << ", which is smaller than the min momentum fraction of " << fMinSharedMomentumFraction << "\n";
2614  returnValue = false;
2615  }
2616  else {
2617  AliDebugStream(4) << "Passed shared momentum fraction with value of " << sharedFraction << "\n";
2618  fHistManager.FillTH1(histName.c_str(), "sharedMomentumFraction");
2619  }
2620 
2621  // Check that the combined jet has a particle-level match
2622  AliEmcalJet * detLevelJet = jet->ClosestJet();
2623  partLevelJet = detLevelJet->ClosestJet();
2624  if (!partLevelJet) {
2625  AliDebugStream(4) << "Jet rejected due to no matching part level jet.\n";
2626  returnValue = false;
2627  }
2628  else {
2629  AliDebugStream(4) << "Det level jet has a required match to a part level jet.\n" << "Part level jet: " << partLevelJet->toString() << "\n";
2630  fHistManager.FillTH1(histName.c_str(), "partLevelMatchedJet");
2631  }
2632 
2633  // Check the matching distance between the combined and pp det-level jets, if a value has been set
2634  if (fMaxMatchedJetDistance > 0) {
2635  double matchedJetDistance = jet->ClosestJetDistance();
2636  if (matchedJetDistance > fMaxMatchedJetDistance) {
2637  AliDebugStream(4) << "Jet rejected due to matching distance of " << matchedJetDistance << ", which is larger than the max distance of " << fMaxMatchedJetDistance << "\n";
2638  returnValue = false;
2639  }
2640  else {
2641  AliDebugStream(4) << "Jet passed distance cut with distance of " << matchedJetDistance << "\n";
2642  fHistManager.FillTH1(histName.c_str(), "jetDistance");
2643  }
2644  }
2645 
2646  // Record all cuts passed
2647  if (returnValue == true) {
2648  fHistManager.FillTH1(histName.c_str(), "passedAllCuts");
2649  }
2650  }
2651  else {
2652  AliDebugStream(5) << "Rejected jet because it was not matched to a external event jet.\n";
2653  fHistManager.FillTH1(histName.c_str(), "noMatch");
2654  returnValue = false;
2655  }
2656 
2657  if (returnValue) {
2658  return partLevelJet;
2659  }
2660  return 0;
2661 }
2662 
2667  const char *ntracks,
2668  const char *nclusters,
2669  const char *nGenLev,
2670  const Double_t minTrPt,
2671  const Double_t minClPt,
2672  const char *suffix)
2673 {
2674  // Get the pointer to the existing analysis manager via the static access method.
2675  //==============================================================================
2676  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2677  if (!mgr)
2678  {
2679  ::Error("AddTaskEmcalJetPerformance", "No analysis manager to connect to.");
2680  return 0;
2681  }
2682 
2683  // Check the analysis type using the event handlers connected to the analysis manager.
2684  //==============================================================================
2685  AliVEventHandler* handler = mgr->GetInputEventHandler();
2686  if (!handler)
2687  {
2688  ::Error("AddTaskEmcalJetPerformance", "This task requires an input event handler");
2689  return 0;
2690  }
2691 
2692  enum EDataType_t {
2693  kUnknown,
2694  kESD,
2695  kAOD
2696  };
2697 
2698  EDataType_t dataType = kUnknown;
2699 
2700  if (handler->InheritsFrom("AliESDInputHandler")) {
2701  dataType = kESD;
2702  }
2703  else if (handler->InheritsFrom("AliAODInputHandler")) {
2704  dataType = kAOD;
2705  }
2706 
2707  //-------------------------------------------------------
2708  // Init the task and do settings
2709  //-------------------------------------------------------
2710 
2711  TString trackName(ntracks);
2712  TString clusName(nclusters);
2713 
2714  if (trackName == "usedefault") {
2715  if (dataType == kESD) {
2716  trackName = "Tracks";
2717  }
2718  else if (dataType == kAOD) {
2719  trackName = "tracks";
2720  }
2721  else {
2722  trackName = "";
2723  }
2724  }
2725 
2726  if (clusName == "usedefault") {
2727  if (dataType == kESD) {
2728  clusName = "CaloClusters";
2729  }
2730  else if (dataType == kAOD) {
2731  clusName = "caloClusters";
2732  }
2733  else {
2734  clusName = "";
2735  }
2736  }
2737 
2738  TString name("AliAnalysisTaskEmcalJetPerformance");
2739  if (!trackName.IsNull()) {
2740  name += "_";
2741  name += trackName;
2742  }
2743  if (!clusName.IsNull()) {
2744  name += "_";
2745  name += clusName;
2746  }
2747  if (strcmp(suffix,"") != 0) {
2748  name += "_";
2749  name += suffix;
2750  }
2751 
2753  // Configure jet performance task
2755 
2757  // Create track and cluster containers with the standard cuts
2758 
2759  AliParticleContainer* partCont = 0;
2760  if (trackName == "mcparticles") {
2761  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
2762  partCont = mcpartCont;
2763  }
2764  else if (trackName == "tracks" || trackName == "Tracks") {
2765  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
2766  partCont = trackCont;
2767  }
2768  if (partCont) partCont->SetParticlePtCut(minTrPt);
2769  if (partCont) task->AdoptParticleContainer(partCont);
2770 
2771  // Add the generator-level container, if specified
2772  if (nGenLev && strcmp(nGenLev,"")!=0) {
2773  AliMCParticleContainer* mcpartCont = task->AddMCParticleContainer(nGenLev);
2774  mcpartCont->SelectPhysicalPrimaries(kTRUE);
2775  mcpartCont->SetParticlePtCut(0);
2776  }
2777 
2778  // Add the cluster container
2779  AliClusterContainer* clusCont = 0;
2780  if (!clusName.IsNull()) {
2781  clusCont = new AliClusterContainer(clusName);
2782  clusCont->SetClusECut(0.);
2783  clusCont->SetClusPtCut(0.);
2784  }
2785  if (clusCont) task->AdoptClusterContainer(clusCont);
2786 
2787  //-------------------------------------------------------
2788  // Final settings, pass to manager and set the containers
2789  //-------------------------------------------------------
2790 
2791  mgr->AddTask(task);
2792 
2793  // Create containers for input/output
2794  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
2795  TString contname(name);
2796  contname += "_histos";
2797  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
2798  TList::Class(),AliAnalysisManager::kOutputContainer,
2799  Form("%s", AliAnalysisManager::GetCommonFileName()));
2800  mgr->ConnectInput (task, 0, cinput1 );
2801  mgr->ConnectOutput (task, 1, coutput1 );
2802 
2803  return task;
2804 }
2805 
Int_t pdg
static AliAnalysisTaskEmcalJetPerformance * AddTaskEmcalJetPerformance(const char *ntracks="usedefault", const char *nclusters="usedefault", const char *nGenLev="mcparticles", const Double_t minTrPt=0.15, const Double_t minClPt=0.30, const char *suffix="")
Double_t Area() const
Definition: AliEmcalJet.h:130
TObjArray fClusterCollArray
cluster collection array
void SetParticlePtCut(Double_t cut)
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
TList * fEventCutList
! Output list for event cut histograms
Double_t GetJetPt(const AliEmcalJet *jet, Double_t rho)
const char * title
Definition: MakeQAPdf.C:27
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
Bool_t fUseAliEventCuts
Flag to use AliEventCuts (otherwise AliAnalysisTaskEmcal will be used)
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:328
UInt_t fOffTrigger
offline trigger for event selection
void AdoptParticleContainer(AliParticleContainer *cont)
const AliClusterIterableMomentumContainer accepted_momentum() const
Int_t nbinsy
void FillParticleCompositionJetHistograms(const AliMCEvent *mcevent)
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
AliEmcalJet * MatchedJet() const
Definition: AliEmcalJet.h:331
Declaration of class AliTLorentzVector.
static AliEmcalDownscaleFactorsOCDB * Instance()
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.
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
AliMCParticleContainer * fGeneratorLevel
! generator level container
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
Container for particles within the EMCAL framework.
void SetCaloTriggerPatchInfoName(const char *n)
Bool_t fPlotClusterHistograms
Set whether to plot cluster histograms.
Double_t GetDeltaR(const AliTLorentzVector *part, Double_t etaRef, Double_t phiRef)
void AdoptClusterContainer(AliClusterContainer *cont)
TObjArray fParticleCollArray
particle/track collection array
virtual Bool_t AcceptMCParticle(const AliAODMCParticle *vp, UInt_t &rejectionReason) const
bool AddQAPlotsToList(TList *list)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const Int_t nPtBins
Bool_t fComputeMBDownscaling
Set whether to compute and plot MB downscaling factors.
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
AliEmcalEmbeddingQA fEmbeddingQA
! QA hists for embedding (will only be added if embedding)
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
bool IsInitialized() const
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.
int Int_t
Definition: External.C:63
AliJetContainer * fMCJetContainer
!Pointer to jet container of truth-level jets
Double_t fTrackMatchingDeltaEtaMax
Maximum delta-eta to consider a track to be matched to a cluster.
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
AliEMCALGeometry * fGeom
!emcal geometry
UInt_t GetJetAcceptanceType() const
Definition: AliEmcalJet.h:367
Double_t Phi_0_2pi() const
Implementation of task to embed external events.
TString toString() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Double_t fTrackMatchingDeltaPhiMax
Maximum delta-phi to consider a track to be matched to a cluster.
Photon+Photon (not from same pi0) are the only contributors.
BeamType fForceBeamType
forced beam type
Bool_t fUseManualEventCuts
Flag to use manual event cuts.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Double_t fCent
!event centrality
Photon (direct or decay) is the only contributor.
Two particles from merged pi0 are the only contributors.
Bool_t fPlotJetHistograms
Set whether to enable inclusive jet histograms.
Int_t fNEoverPBins
! number of variable E/p bins
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
Calculates the fraction of momentum z of part 1 w.r.t. part 2 in the direction of part 2...
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
AliEventCuts fEventCuts
event selection utility
Bool_t fPlotMatchedJetHistograms
Set whether to plot matched jet histograms (must run ResponseMaker first)
TObjArray fJetCollArray
jet collection array
Bool_t fkEMCEJE
! flag telling whether the event is "triggered" or not in "simulation"
AliRhoParameter * GetRhoParameter()
const AliEmcalJet * GetMatchedPartLevelJet(const AliJetContainer *jets, const AliEmcalJet *jet, const std::string &histName)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Photon+Hadron are the only contributors, with Photon leading.
Int_t nbinsx
virtual Bool_t IsEventSelected()
Performing event selection.
Bool_t fUseResponseMaker
Flag to use Response Maker rather than JetTagger.
Bool_t IsHadron(const ContributorType contributor)
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
Bool_t fComputeBackground
Set whether to enable study of background.
Float_t GetJetRadius() const
Double_t fMinSharedMomentumFraction
Minimum shared momentum (pp det-level track pT in combined jet) / (pp det-level track pT) ...
AliEmcalList * fOutput
!output list
ContributorType GetContributorType(const AliVCluster *clus, const AliMCEvent *mcevent, Int_t label)
Handler for downscale factors for various triggers obtained from the OCDB.
AliMCParticleContainer * GetMCParticleContainer(Int_t i=0) const
Int_t fNPtHistBins
! number of variable pt bins
Definition: External.C:220
Bool_t fPlotParticleCompositionHistograms
Set whether to plot jet composition histograms.
void SelectPhysicalPrimaries(Bool_t s)
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Int_t fNM02HistBins
! number of variable M02 bins
Double_t fMBUpscaleFactor
! inverse of downscale factor, for MB trigger
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
const AliTrackIterableMomentumContainer accepted_momentum() const
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
void SetClusPtCut(Double_t cut)
void SetRejectionReasonLabels(TAxis *axis)
void UserCreateOutputObjects()
Main initialization function on the worker.
Bool_t fDoTriggerSimulation
Set whether to perform a simple trigger simulation.
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
EDataType_t
Switch for the data type.
void SetClusECut(Double_t cut)
Double_t fMedianDCal
! median patch energy in DCal, per event
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.
Double_t fMedianEMCal
! median patch energy in EMCal, per event
Container structure for EMCAL clusters.
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:74
Container for MC-true particles within the EMCAL framework.
virtual AliAODMCParticle * GetMCParticleWithLabel(Int_t lab) const
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
Container for jet within the EMCAL jet framework.
Hadron+Photon are the only contributors, with Hadron leading.
void FillParticleCompositionClusterHistograms(const AliMCEvent *mcevent)
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.
Bool_t fPlotDCal
Set whether to enable several DCal-specific histograms.
Double_t fMaxMatchedJetDistance
Maximum distance between two matched jets.
const AliJetIterableContainer all() const
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()