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