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