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