AliPhysics  a60a912 (a60a912)
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 #include "AliAnalysisManager.h"
44 #include <AliVEventHandler.h>
45 #include "AliTLorentzVector.h"
46 #include "AliEmcalJet.h"
47 #include "AliRhoParameter.h"
48 #include "AliJetContainer.h"
49 #include "AliParticleContainer.h"
50 #include "AliClusterContainer.h"
51 #include "AliEMCALGeometry.h"
52 #include "AliOADBContainer.h"
53 #include "AliEmcalDownscaleFactorsOCDB.h"
54 #include "AliEMCALTriggerPatchInfo.h"
55 #include "AliAnalysisTaskEmcalEmbeddingHelper.h"
56 #include "AliMCEvent.h"
57 
58 #include "AliAnalysisTaskEmcalJetPerformance.h"
59 
61 ClassImp(AliAnalysisTaskEmcalJetPerformance);
63 
67 AliAnalysisTaskEmcalJetPerformance::AliAnalysisTaskEmcalJetPerformance() :
68  AliAnalysisTaskEmcalJet(),
69  fPlotJetHistograms(kFALSE),
70  fPlotClusterHistograms(kFALSE),
71  fPlotCellNonlinearityHistograms(kFALSE),
72  fPlotParticleCompositionHistograms(kFALSE),
73  fComputeBackground(kFALSE),
74  fDoTriggerSimulation(kFALSE),
75  fPlotMatchedJetHistograms(kFALSE),
76  fComputeMBDownscaling(kFALSE),
77  fPlotDCal(kFALSE),
78  fDoClosureTest(kFALSE),
79  fFillChargedFluctuations(kFALSE),
80  fMinPt(-100),
81  fMaxPt(250),
82  fNEtaBins(40),
83  fNPhiBins(200),
84  fNCentHistBins(0),
85  fCentHistBins(0),
86  fNPtHistBins(0),
87  fPtHistBins(0),
88  fNM02HistBins(0),
89  fM02HistBins(0),
90  fNEoverPBins(0),
91  fEoverPBins(0),
92  fTrackMatchingDeltaEtaMax(0.015),
93  fTrackMatchingDeltaPhiMax(0.030),
94  fMBUpscaleFactor(1.),
95  fMedianEMCal(0.),
96  fMedianDCal(0.),
97  fkEMCEJE(kFALSE),
98  fDoTriggerResponse(kFALSE),
99  fDoJetMatchingGeometrical(kFALSE),
100  fDoJetMatchingMCFraction(kFALSE),
101  fEmbeddingQA(),
102  fMCJetContainer(nullptr),
103  fDetJetContainer(nullptr),
104  fDetJetContainerPPIntermediate(nullptr),
105  fRequireMatchedJetAccepted(kFALSE),
106  fJetMatchingR(0.),
107  fMinSharedMomentumFraction(0.5),
108  fMCJetMinMatchingPt(0.15),
109  fPlotJetMatchCandThresh(1.),
110  fUseAliEventCuts(kTRUE),
111  fEventCuts(0),
112  fEventCutList(0),
113  fUseManualEventCuts(kFALSE),
114  fGeneratorLevel(0),
115  fHistManager()
116 {
117  GenerateHistoBins();
118 }
119 
125 AliAnalysisTaskEmcalJetPerformance::AliAnalysisTaskEmcalJetPerformance(const char *name) :
126  AliAnalysisTaskEmcalJet(name, kTRUE),
127  fPlotJetHistograms(kFALSE),
128  fPlotClusterHistograms(kFALSE),
129  fPlotCellNonlinearityHistograms(kFALSE),
130  fPlotParticleCompositionHistograms(kFALSE),
131  fComputeBackground(kFALSE),
132  fDoTriggerSimulation(kFALSE),
133  fPlotMatchedJetHistograms(kFALSE),
134  fComputeMBDownscaling(kFALSE),
135  fPlotDCal(kFALSE),
136  fDoClosureTest(kFALSE),
137  fFillChargedFluctuations(kFALSE),
138  fMinPt(-100),
139  fMaxPt(250),
140  fNEtaBins(40),
141  fNPhiBins(200),
142  fNCentHistBins(0),
143  fCentHistBins(0),
144  fNPtHistBins(0),
145  fPtHistBins(0),
146  fNM02HistBins(0),
147  fM02HistBins(0),
148  fNEoverPBins(0),
149  fEoverPBins(0),
150  fTrackMatchingDeltaEtaMax(0.015),
151  fTrackMatchingDeltaPhiMax(0.030),
152  fMBUpscaleFactor(1.),
153  fMedianEMCal(0.),
154  fMedianDCal(0.),
155  fkEMCEJE(kFALSE),
156  fDoTriggerResponse(kFALSE),
157  fDoJetMatchingGeometrical(kFALSE),
158  fDoJetMatchingMCFraction(kFALSE),
159  fEmbeddingQA(),
160  fMCJetContainer(nullptr),
161  fDetJetContainer(nullptr),
162  fDetJetContainerPPIntermediate(nullptr),
163  fRequireMatchedJetAccepted(kFALSE),
164  fJetMatchingR(0.),
165  fMinSharedMomentumFraction(0.5),
166  fMCJetMinMatchingPt(0.15),
167  fPlotJetMatchCandThresh(1.),
168  fUseAliEventCuts(kTRUE),
169  fEventCuts(0),
170  fEventCutList(0),
171  fUseManualEventCuts(kFALSE),
172  fGeneratorLevel(0),
173  fHistManager(name)
174 {
175  GenerateHistoBins();
176 }
177 
181 AliAnalysisTaskEmcalJetPerformance::~AliAnalysisTaskEmcalJetPerformance()
182 {
183 }
184 
188 void AliAnalysisTaskEmcalJetPerformance::GenerateHistoBins()
189 {
190  fNCentHistBins = 4;
191  fCentHistBins = new Double_t[fNCentHistBins+1];
192  fCentHistBins[0] = 0;
193  fCentHistBins[1] = 10;
194  fCentHistBins[2] = 30;
195  fCentHistBins[3] = 50;
196  fCentHistBins[4] = 90;
197 
198  fNPtHistBins = 82;
199  fPtHistBins = new Double_t[fNPtHistBins+1];
200  GenerateFixedBinArray(6, 0, 0.3, fPtHistBins);
201  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
202  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
203  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
204  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
205  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
206  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
207 
208  fNM02HistBins = 81;
209  fM02HistBins = new Double_t[fNM02HistBins+1];
210  GenerateFixedBinArray(35, 0, 0.7, fM02HistBins);
211  GenerateFixedBinArray(6, 0.7, 1., fM02HistBins+35);
212  GenerateFixedBinArray(20, 1., 3., fM02HistBins+41);
213  GenerateFixedBinArray(10, 3., 5., fM02HistBins+61);
214  GenerateFixedBinArray(10, 5., 10., fM02HistBins+71);
215 
216  fNEoverPBins = 47;
217  fEoverPBins = new Double_t[fNEoverPBins+1];
218  GenerateFixedBinArray(30, 0, 1.5, fEoverPBins);
219  GenerateFixedBinArray(10, 1.5, 3.5, fEoverPBins+30);
220  GenerateFixedBinArray(7, 3.5, 10.5, fEoverPBins+40);
221 }
222 
227 void AliAnalysisTaskEmcalJetPerformance::UserCreateOutputObjects()
228 {
229  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
230 
231  // Intialize AliEventCuts
232  if (fUseAliEventCuts) {
233  fEventCutList = new TList();
234  fEventCutList ->SetOwner();
235  fEventCutList ->SetName("EventCutOutput");
236 
237  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
238  if(fUseManualEventCuts==1)
239  {
240  fEventCuts.SetManualMode();
241  fEventCuts.fMC = false;
242  fEventCuts.SetupLHC15o();
243  fEventCuts.fUseVariablesCorrelationCuts = true;
244  }
245  fEventCuts.AddQAplotsToList(fEventCutList);
246  fOutput->Add(fEventCutList);
247  }
248 
249  // Get the MC particle branch, in case it exists
250  fGeneratorLevel = GetMCParticleContainer("mcparticles");
251 
252  // Get jet containers, in order to check the jet acceptance criteria
253  // For geometrical matching, it is expected that there is one det-level and one truth-level container
254  // For MC fraction matching, it is expected that there is one Pb-Pb det-level, one pp det-level, and one pp truth-level container
255  if (fPlotMatchedJetHistograms) {
256 
257  if (fDoJetMatchingGeometrical) {
258  Printf("Geometrical jet matching enabled.");
259 
260  for (Int_t i=0; i<2; i++) {
261  auto jetCont = GetJetContainer(i);
262  TString jetContName = jetCont->GetName();
263  if (jetContName.Contains("mcparticles")) {
264  fMCJetContainer = jetCont;
265  }
266  else {
267  fDetJetContainer = jetCont;
268  }
269  }
270 
271  }
272  else if (fDoJetMatchingMCFraction) {
273  Printf("MC-fraction jet matching enabled.");
274 
275  for (Int_t i=0; i<3; i++) {
276  auto jetCont = GetJetContainer(i);
277  TString jetContName = jetCont->GetName();
278  if (jetContName.Contains("mcparticles")) {
279  fMCJetContainer = jetCont;
280  }
281  else if (jetContName.Contains("Combined")) {
282  fDetJetContainer = jetCont;
283  }
284  else {
285  fDetJetContainerPPIntermediate = jetCont;
286  }
287  }
288 
289  }
290 
291  if (!fMCJetContainer) {
292  Printf("No MC jet container found!");
293  }
294  Printf("mcJetContainer: %s", fMCJetContainer->GetName());
295 
296  if (!fDetJetContainer) {
297  Printf("No det-level jet container found!");
298  }
299  Printf("det-level JetContainer: %s", fDetJetContainer->GetName());
300 
301  if (fDoJetMatchingMCFraction) {
302  if (!fDetJetContainerPPIntermediate) {
303  Printf("No intermediate pp det-level jet container found, despite MC-fraction matching enabled!");
304  }
305  Printf("Intermediate pp det-level JetContainer: %s", fDetJetContainerPPIntermediate->GetName());
306  }
307 
308  }
309 
310  // Allocate histograms
311  if (fPlotJetHistograms) {
312  AllocateJetHistograms();
313  }
314  if (fPlotClusterHistograms) {
315  AllocateClusterHistograms();
316  }
317  if (fPlotCellNonlinearityHistograms) {
318  AllocateCellNonlinearityHistograms();
319  }
320  if (fPlotParticleCompositionHistograms) {
321  AllocateParticleCompositionHistograms();
322  }
323  if (fComputeBackground) {
324  AllocateBackgroundHistograms();
325  }
326  if (fDoTriggerSimulation) {
327  AllocateTriggerSimHistograms();
328  }
329  if (fPlotMatchedJetHistograms) {
330  AllocateMatchedJetHistograms();
331  }
332 
333  // Initialize embedding QA
334  const AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
335  if (embeddingHelper) {
336  bool res = fEmbeddingQA.Initialize();
337  if (res) {
338  fEmbeddingQA.AddQAPlotsToList(fOutput);
339  }
340  }
341 
342  TIter next(fHistManager.GetListOfHistograms());
343  TObject* obj = 0;
344  while ((obj = next())) {
345  fOutput->Add(obj);
346  }
347 
348  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
349 }
350 
351 /*
352  * This function allocates the histograms for single jets.
353  * A set of histograms is allocated per each jet container.
354  */
355 void AliAnalysisTaskEmcalJetPerformance::AllocateJetHistograms()
356 {
357  TString histname;
358  TString title;
359 
360  Int_t nPtBins1 = TMath::CeilNint(fMaxPt-fMinPt);
361  Int_t nPtBins2 = TMath::CeilNint((fMaxPt-fMinPt)/2);
362  Int_t nPtBins5 = TMath::CeilNint((fMaxPt-fMinPt)/5);
363 
364  AliJetContainer* jets = 0;
365  TIter nextJetColl(&fJetCollArray);
366  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
367 
368  if (fPlotMatchedJetHistograms && fDoJetMatchingMCFraction && jets != fMCJetContainer) {
369  if (!fDoClosureTest) {
370  continue; // don't plot det-level histograms if embedding
371  }
372  }
373 
374  // Jet rejection reason
375  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
376  if (fForceBeamType == kAA) {
377  title = histname + ";Centrality (%);Rejection reason;#it{p}_{T,corr} (GeV/#it{c});counts";
378  TH3* hist = fHistManager.CreateTH3(histname.Data(), title.Data(), 10, 0, 100, 32, 0, 32, nPtBins5, fMinPt, fMaxPt);
379  SetRejectionReasonLabels(hist->GetYaxis());
380  } else {
381  title = histname + ";Rejection reason;#it{p}_{T,corr} (GeV/#it{c});counts";
382  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, nPtBins5, fMinPt, fMaxPt);
383  SetRejectionReasonLabels(hist->GetXaxis());
384  }
385 
386  // Rho vs. Centrality
387  if (!jets->GetRhoName().IsNull()) {
388  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
389  if (fForceBeamType == kAA) {
390  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
391  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
392  }
393  else{
394  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
395  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 250, 0, 50);
396  }
397  }
398 
399  // (Centrality, pT, NEF)
400  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
401  Int_t nbinsy = nPtBins1; Int_t miny = fMinPt; Int_t maxy = fMaxPt;
402  Int_t nbinsz = 50; Int_t minz = 0; Int_t maxz = 1.;
403 
404  histname = TString::Format("%s/JetHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
405  if (fForceBeamType == kAA) {
406  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
407  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
408  }
409  else {
410  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
411  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
412  }
413 
414  if (fPlotDCal) {
415  histname = TString::Format("%s/JetHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
416  if (fForceBeamType == kAA) {
417  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
418  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
419  }
420  else {
421  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
422  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
423  }
424  }
425 
426  // (Centrality, pT upscaled, calo type)
427  if (fComputeMBDownscaling) {
428  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
429  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
430  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins2, fMinPt, fMaxPt, 2, -0.5, 1.5, "s");
431  }
432 
433  // pT-leading vs. pT
434  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
435  if (fForceBeamType == kAA) {
436  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
437  fHistManager.CreateTH3(histname.Data(), title.Data(), 10, 0, 100, nPtBins2, fMinPt, fMaxPt, fMaxPt, 0, fMaxPt);
438  }
439  else {
440  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
441  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins2, fMinPt, fMaxPt, nPtBins1, fMinPt, fMaxPt);
442  }
443 
444  // A vs. pT
445  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
446  if (fForceBeamType == kAA) {
447  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}/#piR^{2}";
448  fHistManager.CreateTH3(histname.Data(), title.Data(), 10, 0, 100, nPtBins2, fMinPt, fMaxPt, 75, 0, 3);
449  }
450  else {
451  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}/#piR^{2}";
452  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins2, fMinPt, fMaxPt, 100, 0, 3);
453  }
454 
455  // (Centrality, pT, z-leading (charged))
456  nbinsx = 20; minx = 0; maxx = 100;
457  nbinsy = nPtBins2; miny = fMinPt; maxy = fMaxPt;
458  nbinsz = 50; minz = 0; maxz = 1.;
459 
460  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
461  if (fForceBeamType == kAA) {
462  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
463  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
464  }
465  else {
466  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
467  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
468  }
469 
470  if (fPlotDCal) {
471  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
472  if (fForceBeamType == kAA) {
473  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
474  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
475  }
476  else {
477  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
478  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
479  }
480  }
481 
482  // (Centrality, pT, z (charged))
483  nbinsx = 20; minx = 0; maxx = 100;
484  nbinsy = nPtBins2; miny = fMinPt; maxy = fMaxPt;
485  nbinsz = 50; minz = 0; maxz = 1.;
486 
487  histname = TString::Format("%s/JetHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
488  if (fForceBeamType == kAA) {
489  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
490  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
491  }
492  else {
493  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
494  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
495  }
496 
497  if (fPlotDCal) {
498  histname = TString::Format("%s/JetHistograms/hZVsPtDCal", jets->GetArrayName().Data());
499  if (fForceBeamType == kAA) {
500  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
501  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
502  }
503  else {
504  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
505  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
506  }
507  }
508 
509  // (Centrality, pT, Nconst)
510  nbinsx = 20; minx = 0; maxx = 100;
511  nbinsy = nPtBins2; miny = fMinPt; maxy = fMaxPt;
512  nbinsz = 50; minz = 0; maxz = fMaxPt;
513 
514  histname = TString::Format("%s/JetHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
515  if (fForceBeamType == kAA) {
516  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
517  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
518  }
519  else {
520  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
521  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
522  }
523 
524  if (fPlotDCal) {
525  histname = TString::Format("%s/JetHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
526  if (fForceBeamType == kAA) {
527  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
528  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
529  }
530  else {
531  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
532  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
533  }
534  }
535 
536  // (Centrality, pT) for eta<0 and eta>0
537  if (fForceBeamType == kAA) {
538  nbinsx = 20; minx = 0; maxx = 100;
539  nbinsy = nPtBins1; miny = fMinPt; maxy = fMaxPt;
540 
541  histname = TString::Format("%s/JetHistograms/hEtaPosVsPtEMCal", jets->GetArrayName().Data());
542  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c})";
543  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy);
544 
545  histname = TString::Format("%s/JetHistograms/hEtaNegVsPtEMCal", jets->GetArrayName().Data());
546  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c})";
547  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy);
548  }
549 
550  // (Centrality, jet pT, Enonlincorr - Ehadcorr)
551  if (fForceBeamType == kAA) {
552  nbinsx = 20; minx = 0; maxx = 100;
553  nbinsy = nPtBins2; miny = fMinPt; maxy = fMaxPt;
554  nbinsz = nPtBins2; minz = fMinPt; maxz = fMaxPt;
555 
556  histname = TString::Format("%s/JetHistograms/hDeltaEHadCorr", jets->GetArrayName().Data());
557  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#sum#it{E}_{nonlincorr} - #it{E}_{hadcorr}";
558  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
559  }
560 
561  // (Median patch energy, calo type, jet pT, centrality)
562  if (fDoTriggerSimulation) {
563  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
564  title = histname + ";#it{E}_{patch,med};type;#it{p}_{T}^{corr} (GeV/#it{c});Centrality (%)";
565  Int_t nbins5[4] = {100, 2, nPtBins2, 50};
566  Double_t min5[4] = {0,-0.5, fMinPt, 0};
567  Double_t max5[4] = {50,1.5, fMaxPt, 100};
568  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins5, min5, max5);
569  }
570 
571  }
572 
573  // MB downscale factor histogram
574  if (fComputeMBDownscaling) {
575  histname = "Trigger/hMBDownscaleFactor";
576  title = histname + ";Downscale factor;counts";
577  fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 200);
578  }
579 
580 }
581 
582 /*
583  * This function allocates the histograms for the calorimeter performance study.
584  */
585 void AliAnalysisTaskEmcalJetPerformance::AllocateClusterHistograms()
586 {
587  TString histname;
588  TString htitle;
589 
590  const Int_t nRcorrBins = 50;
591  Double_t *RcorrBins = GenerateFixedBinArray(nRcorrBins, 0., 1.);
592  const Int_t nCellBins = 30;
593  Double_t *cellBins = GenerateFixedBinArray(nCellBins, -0.5, 29.5);
594  const Int_t nMatchedTrackBins = 5;
595  Double_t *matchedTrackBins = GenerateFixedBinArray(nMatchedTrackBins, -0.5, 4.5);
596  const Int_t nDeltaEtaBins = 60;
597  Double_t *deltaEtaBins = GenerateFixedBinArray(nDeltaEtaBins, -0.015, 0.015);
598 
601 
602  // Plot M02 distribution (centrality, Eclus nonlincorr, M02)
603  histname = "ClusterHistograms/hM02";
604  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); M02";
605  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, fNM02HistBins, fM02HistBins);
606 
607  // Plot Ncell distribution for M02 > 0.4 and 0.1 < M02 < 0.4 (centrality, Eclus nonlincorr, Ncells)
608  histname = "ClusterHistograms/hNcellsM02G04";
609  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); Ncells";
610  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nCellBins, cellBins);
611 
612  histname = "ClusterHistograms/hNcellsM02L04";
613  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); Ncells";
614  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nCellBins, cellBins);
615 
618 
619  // Plot matched track pT for all clusters, M02 > 0.4 clusters, and 0.1 < M02 < 0.4 clusters (centrality, Eclus nonlincorr, trackPsum)
620  histname = "ClusterHistograms/hMatchedTrackPt";
621  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); #Sigma#it{p}_{track} (GeV/c)";
622  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
623 
624  histname = "ClusterHistograms/hMatchedTrackPtM02G04";
625  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); #Sigma#it{p}_{track} (GeV/c)";
626  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
627 
628  histname = "ClusterHistograms/hMatchedTrackPtM02L04";
629  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); #Sigma#it{p}_{track} (GeV/c)";
630  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
631 
632  // Plot number of matched tracks for all clusters, M02 > 0.4 clusters, and 0.1 < M02 < 0.4 clusters (centrality, Eclus nonlincorr, N matches)
633  histname = "ClusterHistograms/hMatchedTrackN";
634  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); N_{tracks}";
635  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nMatchedTrackBins, matchedTrackBins);
636 
637  histname = "ClusterHistograms/hMatchedTrackNM02G04";
638  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); N_{tracks}";
639  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nMatchedTrackBins, matchedTrackBins);
640 
641  histname = "ClusterHistograms/hMatchedTrackNM02L04";
642  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); N_{tracks}";
643  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nMatchedTrackBins, matchedTrackBins);
644 
645  // Plot M02 distribution for clusters with matched tracks (centrality, Eclus nonlincorr, M02)
646  histname = "ClusterHistograms/hM02Matched";
647  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); M02";
648  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, fNM02HistBins, fM02HistBins);
649 
650  // Plot M02 distribution for clusters without matched tracks (centrality, Eclus nonlincorr, M02)
651  histname = "ClusterHistograms/hM02Unmatched";
652  htitle = histname + ";Centrality (%);#it{E}_{clus} (GeV); M02";
653  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, fNM02HistBins, fM02HistBins);
654 
655  // Plot clus-track deltaEta of matched tracks (deltaEta, Eclus, M02)
656  histname = "ClusterHistograms/hDeltaEtaCentral";
657  htitle = histname + ";#eta_{track} - #eta_{clus};#it{E}_{clus} (GeV); M02";
658  fHistManager.CreateTH3(histname.Data(), htitle.Data(), nDeltaEtaBins, deltaEtaBins, fNPtHistBins, fPtHistBins, fNM02HistBins, fM02HistBins);
659 
660  histname = "ClusterHistograms/hDeltaEtaPeripheral";
661  htitle = histname + ";#eta_{track} - #eta_{clus};#it{E}_{clus} (GeV); M02";
662  fHistManager.CreateTH3(histname.Data(), htitle.Data(), nDeltaEtaBins, deltaEtaBins, fNPtHistBins, fPtHistBins, fNM02HistBins, fM02HistBins);
663 
666 
667  // Plot E/p vs. M02 for 0-10% and 50-90% (Eclus nonlincorr, Eclus nonlincorr / trackPsum, M02)
668  histname = "ClusterHistograms/hEoverPM02Central";
669  htitle = histname + ";#it{E}_{clus} (GeV); #it{E}_{clus} / #Sigma#it{p}_{track} (GeV); M02";
670  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNEoverPBins, fEoverPBins, fNM02HistBins, fM02HistBins);
671 
672  histname = "ClusterHistograms/hEoverPM02Peripheral";
673  htitle = histname + ";#it{E}_{clus} (GeV); #it{E}_{clus} / #Sigma#it{p}_{track} (GeV); M02";
674  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNEoverPBins, fEoverPBins, fNM02HistBins, fM02HistBins);
675 
678 
679  // Plot Rcorr distribution (centrality, trackPSum, Rcorr = (Enonlincorr - Ehadcorr) / trackPSum)
680  histname = "ClusterHistograms/hRcorrVsCent";
681  htitle = histname + ";Centrality (%);#Sigma#it{p}_{track} (GeV); R_{corr} = #frac{#DeltaE_{clus}}{#Sigmap_{track}}";
682  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
683 
684  // Plot Rcorr distribution for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr)
685  histname = "ClusterHistograms/hRcorr0-10";
686  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); R_{corr} = #frac{#DeltaE_{clus}}{#Sigmap_{track}}";
687  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
688 
689  // Plot Rcorr distribution for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr)
690  histname = "ClusterHistograms/hRcorr50-90";
691  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); R_{corr} = #frac{#DeltaE_{clus}}{#Sigmap_{track}}";
692  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
693 
694  // Plot also Rcorr-clus (centrality, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
695  histname = "ClusterHistograms/hRcorrClusVsCent";
696  htitle = histname + ";Centrality (%);#Sigma#it{p}_{track} (GeV); #frac{#DeltaE_{clus}}{E_{clus}}";
697  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
698 
699  // Rcorr-clus for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
700  histname = "ClusterHistograms/hRcorrClus0-10";
701  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); #frac{#DeltaE_{clus}}{E_{clus}}";
702  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
703 
704  // Rcorr-clus for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
705  histname = "ClusterHistograms/hRcorrClus50-90";
706  htitle = histname + ";#it{E}_{clus} (GeV);#Sigma#it{p}_{track} (GeV); #frac{#DeltaE_{clus}}{E_{clus}}";
707  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins, nRcorrBins, RcorrBins);
708 
709  // Plot total track multiplicity
710  histname = "ClusterHistograms/hTrackMultiplicity";
711  htitle = histname + ";N_{tracks};Centrality (%)";
712  fHistManager.CreateTH2(histname.Data(), htitle.Data(), 1000, 0, 10000, 20, 0, 100);
713 }
714 
715 /*
716  * This function allocates the histograms for the cell-level non-linearity study for embedding.
717  * This should be run over MC.
718  */
719 void AliAnalysisTaskEmcalJetPerformance::AllocateCellNonlinearityHistograms()
720 {
721  TString histname;
722  TString htitle;
723 
724  // Plot cluster non-linearity scale factor
725  histname = "CellNonlinearityHistograms/hClusterNonlinearity";
726  htitle = histname + ";E_{cluster}; E_{NonLinCorr}/E";
727  fHistManager.CreateTH2(histname.Data(), htitle.Data(), 1500, 0, 150, 2000, 0.5, 1.5);
728 
729  // Plot cell non-linearity scale factor
730  histname = "CellNonlinearityHistograms/hCellNonlinearity";
731  htitle = histname + ";E_{cell}; E_{NonLinCorr}/E";
732  fHistManager.CreateTH2(histname.Data(), htitle.Data(), 1500, 0, 150, 2000, 0.5, 1.5);
733 
734  // Plot cell non-linearity scale factor for leading cell
735  histname = "CellNonlinearityHistograms/hLeadingCellNonlinearity";
736  htitle = histname + ";E_{leading cell}; E_{NonLinCorr}/E";
737  fHistManager.CreateTH2(histname.Data(), htitle.Data(), 1500, 0, 150, 2000, 0.5, 1.5);
738 
739 }
740 
741 /*
742  * This function allocates the histograms for the jet composition study.
743  */
744 void AliAnalysisTaskEmcalJetPerformance::AllocateParticleCompositionHistograms()
745 {
746  TString histname;
747  TString htitle;
748  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
749 
750  const Int_t nRejBins = 32;
751  Double_t* rejReasonBins = new Double_t[nRejBins+1];
752  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
753  const Int_t nContributorTypes = 11;
754  Double_t *contributorTypeBins = GenerateFixedBinArray(nContributorTypes, -0.5, 10.5);
755  const Int_t nParticleTypes = 17;
756  Double_t *particleTypeBins = GenerateFixedBinArray(nParticleTypes, -0.5, 16.5);
757 
758  AliEmcalContainer* cont = 0;
759  TIter nextClusColl(&fClusterCollArray);
760  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
761 
762  histname = "ClusterHistogramsMC/hClusterRejectionReasonMC";
763  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
764  TH2* histMC2 = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
765  SetRejectionReasonLabels(histMC2->GetXaxis());
766  }
767 
768  // M02 vs. Energy vs. Particle type
769  histname = "ClusterHistogramsMC/hM02VsParticleTypeCentral";
770  htitle = histname + ";M02;#it{E}_{clus} (GeV); Particle type";
771  TH3* hM02VsParticleTypeCentral = fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins, nParticleTypes, particleTypeBins);
772  SetParticleTypeLabels(hM02VsParticleTypeCentral->GetZaxis());
773 
774  histname = "ClusterHistogramsMC/hM02VsParticleTypePeripheral";
775  htitle = histname + ";M02;#it{E}_{clus} (GeV); Particle type";
776  TH3* hM02VsParticleTypePeripheral = fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins, nParticleTypes, particleTypeBins);
777  SetParticleTypeLabels(hM02VsParticleTypePeripheral->GetZaxis());
778 
779  // Plot photon energy in photon-hadron overlap clusters (Centrality, Photon energy, M02)
780  histname = "ClusterHistogramsMC/hPhotonHadronPhotonEnergy";
781  htitle = histname + ";Centrality (%);M02;#it{E}_{photon} (GeV)";
782  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins);
783 
784  // Plot hadron energy in hadron-photon overlap clusters (Centrality, Photon energy, M02)
785  histname = "ClusterHistogramsMC/hHadronPhotonHadronEnergy";
786  htitle = histname + ";Centrality (%);M02;#it{E}_{hadron} (GeV)";
787  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins);
788 
789  if (fPlotJetHistograms) {
790 
791  // M02 vs. Energy vs. Particle type vs. Jet pT, for particles inside jets
792  Int_t dim = 0;
793  TString title[20];
794  Int_t nbins[20] = {0};
795  Double_t min[30] = {0.};
796  Double_t max[30] = {0.};
797  Double_t *binEdges[20] = {0};
798 
799  title[dim] = "M02";
800  nbins[dim] = fNM02HistBins;
801  binEdges[dim] = fM02HistBins;
802  min[dim] = fM02HistBins[0];
803  max[dim] = fM02HistBins[fNM02HistBins];
804  dim++;
805 
806  title[dim] = "#it{E}_{clus} (GeV)";
807  nbins[dim] = fNPtHistBins;
808  binEdges[dim] = fPtHistBins;
809  min[dim] = fPtHistBins[0];
810  max[dim] = fPtHistBins[fNPtHistBins];
811  dim++;
812 
813  title[dim] = "Contributor type";
814  nbins[dim] = nContributorTypes;
815  min[dim] = -0.5;
816  max[dim] = 7.5;
817  binEdges[dim] = contributorTypeBins;
818  dim++;
819 
820  title[dim] = "#it{p}_{T,jet}^{corr}";
821  nbins[dim] = nPtBins;
822  min[dim] = 0;
823  max[dim] = fMaxPt;
824  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
825  dim++;
826 
827  TString thnname = "JetPerformanceMC/hM02VsContributorTypeJets";
828  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
829  for (Int_t i = 0; i < dim; i++) {
830  hn->GetAxis(i)->SetTitle(title[i]);
831  hn->SetBinEdges(i, binEdges[i]);
832  }
833 
834  // Particle composition inside each jet -- jet pT vs. particle type vs. particle number vs. particle pT sum
835  // (One entry per jet for each particle type)
836  dim = 0;
837 
838  title[dim] = "#it{p}_{T,jet}^{corr}";
839  nbins[dim] = nPtBins;
840  min[dim] = 0;
841  max[dim] = fMaxPt;
842  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
843  dim++;
844 
845  title[dim] = "Contributor type";
846  nbins[dim] = nContributorTypes;
847  min[dim] = -0.5;
848  max[dim] = 7.5;
849  binEdges[dim] = contributorTypeBins;
850  dim++;
851 
852  title[dim] = "N";
853  nbins[dim] = 30;
854  min[dim] = -0.5;
855  max[dim] = 29.5;
856  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
857  dim++;
858 
859  title[dim] = "#it{p}_{T,sum} (GeV)";
860  nbins[dim] = fNPtHistBins;
861  binEdges[dim] = fPtHistBins;
862  min[dim] = fPtHistBins[0];
863  max[dim] = fPtHistBins[fNPtHistBins];
864  dim++;
865 
866  thnname = "JetPerformanceMC/hJetComposition";
867  THnSparse* thn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
868  for (Int_t i = 0; i < dim; i++) {
869  thn->GetAxis(i)->SetTitle(title[i]);
870  thn->SetBinEdges(i, binEdges[i]);
871  }
872 
873  // Hadronic calo energy in each jet
874 
875  // Jet pT vs. Summed energy of hadronic clusters without a matched track
876  histname = "JetPerformance/hHadCaloEnergyUnmatched";
877  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
878  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
879 
880  // Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction)
881  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
882  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
883  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
884 
885  // Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction)
886  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
887  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
888  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
889 
890  }
891 
892 }
893 
894 /*
895  * This function sets axis labels for particle type histograms.
896  */
897 void AliAnalysisTaskEmcalJetPerformance::SetParticleTypeLabels(TAxis* axis)
898 {
899  axis->SetBinLabel(1, "SinglePhoton");
900  axis->SetBinLabel(2, "SingleElectron");
901  axis->SetBinLabel(3, "SingleChargedPion");
902  axis->SetBinLabel(4, "SingleProton");
903  axis->SetBinLabel(5, "SingleAntiProton");
904  axis->SetBinLabel(6, "SingleChargedKaon");
905  axis->SetBinLabel(7, "SingleK0L");
906  axis->SetBinLabel(8, "SingleNeutron");
907  axis->SetBinLabel(9, "SingleAntiNeutron");
908  axis->SetBinLabel(10, "SingleOther");
909  axis->SetBinLabel(11, "PhotonHadron");
910  axis->SetBinLabel(12, "HadronPhoton");
911  axis->SetBinLabel(13, "MergedPi0");
912  axis->SetBinLabel(14, "PhotonPhotonOther");
913  axis->SetBinLabel(15, "HadronHadron");
914  axis->SetBinLabel(16, "TwoContributorsOther");
915  axis->SetBinLabel(17, "MoreThanTwoContributors");
916 }
917 
918 /*
919  * This function allocates background subtraction histograms, if enabled.
920  * A set of histograms is allocated per each jet container.
921  */
922 void AliAnalysisTaskEmcalJetPerformance::AllocateBackgroundHistograms()
923 {
924  TString histname;
925  TString title;
926 
927  AliJetContainer* jets = 0;
928  TIter nextJetColl(&fJetCollArray);
929  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
930 
931  // EMCal
932  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
933  title = histname + ";Centrality;Scale factor;counts";
934  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
935  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
936  if (fForceBeamType == kAA) {
937  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
938  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
939  }
940  else{
941  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
942  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 1300, -25, 40);
943  }
944  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCalExcl", jets->GetArrayName().Data());
945  if (fForceBeamType == kAA) {
946  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
947  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
948  }
949  else{
950  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
951  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 1300, -25, 40);
952  }
953  if(fFillChargedFluctuations){
954  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCalCharged", jets->GetArrayName().Data());
955  if (fForceBeamType == kAA) {
956  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
957  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
958  }
959  else{
960  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
961  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 1300, -25, 40);
962  }
963  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCalChargedExcl", jets->GetArrayName().Data());
964  if (fForceBeamType == kAA) {
965  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
966  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
967  }
968  else{
969  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
970  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 1300, -25, 40);
971  }
972  }
973  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCalFid", jets->GetArrayName().Data());
974  title = histname + ";Centrality;Scale factor;counts";
975  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
976 
977  // DCal
978  if (fPlotDCal) {
979  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCal", jets->GetArrayName().Data());
980  title = histname + ";Centrality;Scale factor;counts";
981  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
982 
983  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtDCal", jets->GetArrayName().Data());
984  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
985  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
986 
987  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCalFid", jets->GetArrayName().Data());
988  title = histname + ";Centrality;Scale factor;counts";
989  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
990  }
991 
992  }
993 }
994 
995 /*
996  * This function allocates the histograms for single jets, when the "simulated" trigger has been fired.
997  * A set of histograms is allocated per each jet container.
998  */
999 void AliAnalysisTaskEmcalJetPerformance::AllocateTriggerSimHistograms()
1000 {
1001  TString histname;
1002  TString title;
1003  Int_t nPtBins1 = TMath::CeilNint(fMaxPt-fMinPt);
1004  Int_t nPtBins2 = TMath::CeilNint((fMaxPt-fMinPt)/2);
1005  Int_t nPtBins5 = TMath::CeilNint((fMaxPt-fMinPt)/5);
1006 
1007  //----------------------------------------------
1008  // Trigger patch histograms
1009 
1010  // patch eta vs. phi
1011  histname = "TriggerSimHistograms/hEtaVsPhi";
1012  title = histname + ";#eta_{patch} (rad);#phi_{patch} (rad)";
1013  fHistManager.CreateTH2(histname.Data(), title.Data(), 140, -0.7, 0.7, 500, 1., 6.);
1014 
1015  // N patches
1016  histname = "TriggerSimHistograms/hNPatches";
1017  title = histname + ";#it{N}_{patches};type";
1018  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, 0, 200, 2, -0.5, 1.5);
1019 
1020  // patch E vs. centrality
1021  histname = "TriggerSimHistograms/hPatchE";
1022  title = histname + ";Centrality (%);#it{E}_{patch} (GeV)";
1023  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, nPtBins2, 0, fMaxPt);
1024 
1025  // patch median vs. Centrality
1026  histname = "TriggerSimHistograms/hPatchMedianE";
1027  title = histname + ";Centrality (%);#it{E}_{patch,med} (GeV);type";
1028  fHistManager.CreateTH3(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 50, 2, -0.5, 1.5);
1029 
1030  if (fDoTriggerResponse) {
1031  histname = "TriggerSimHistograms/hMaxPatchResponseMatrix";
1032  title = histname + ";#it{E}_{maxpatch,det} (GeV);#it{E}_{maxpatch,truth} (GeV)";
1033  fHistManager.CreateTH2(histname.Data(), title.Data(), 400, 0, 200, 400, 0, 200);
1034 
1035  histname = "TriggerSimHistograms/hEMCalEnergyResponseMatrix";
1036  title = histname + ";#Sigma#it{E}_{clus,EMCal} (GeV);#it{E}_{EMCal,truth} (GeV)";
1037  fHistManager.CreateTH2(histname.Data(), title.Data(), 500, 0, 500, 500, 0, 500);
1038  }
1039 
1040  //----------------------------------------------
1041  // Jet histograms for "triggered" events
1042  AliJetContainer* jets = 0;
1043  TIter nextJetColl(&fJetCollArray);
1044  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1045 
1046  // Jet rejection reason
1047  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1048  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
1049  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, nPtBins5, fMinPt, fMaxPt);
1050  SetRejectionReasonLabels(hist->GetXaxis());
1051 
1052  // Rho vs. Centrality
1053  if (!jets->GetRhoName().IsNull()) {
1054  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
1055  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
1056  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
1057  }
1058 
1059  // (Centrality, pT, NEF)
1060  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
1061  Int_t nbinsy = nPtBins1; Int_t miny = fMinPt; Int_t maxy = fMaxPt;
1062  Int_t nbinsz = 50; Int_t minz = 0; Int_t maxz = 1.;
1063 
1064  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
1065  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
1066  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1067 
1068  if (fPlotDCal) {
1069  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
1070  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
1071  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1072  }
1073 
1074  // pT-leading vs. pT
1075  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1076  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
1077  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins2, fMinPt, fMaxPt, nPtBins2, fMinPt, fMaxPt);
1078 
1079  // A vs. pT
1080  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
1081  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}/#piR^{2}";
1082  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins2, fMinPt, fMaxPt, 75, 0, 3);
1083 
1084  // (Centrality, pT, z-leading (charged))
1085  nbinsx = 20; minx = 0; maxx = 100;
1086  nbinsy = nPtBins2; miny = fMinPt; maxy = fMaxPt;
1087  nbinsz = 50; minz = 0; maxz = 1.;
1088 
1089  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
1090  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
1091  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1092 
1093  if (fPlotDCal) {
1094  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
1095  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
1096  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1097  }
1098 
1099  // z (charged) vs. pT
1100  nbinsx = 20; minx = 0; maxx = 100;
1101  nbinsy = nPtBins2; miny = fMinPt; maxy = fMaxPt;
1102  nbinsz = 50; minz = 0; maxz = 1.;
1103 
1104  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
1105  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
1106  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1107 
1108  if (fPlotDCal) {
1109  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtDCal", jets->GetArrayName().Data());
1110  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
1111  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1112  }
1113 
1114  // (Centrality, pT, Nconst)
1115  nbinsx = 20; minx = 0; maxx = 100;
1116  nbinsy = nPtBins2; miny = fMinPt; maxy = fMaxPt;
1117  nbinsz = 50; minz = 0; maxz = fMaxPt;
1118 
1119  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
1120  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
1121  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1122 
1123  if (fPlotDCal) {
1124  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
1125  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
1126  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1127  }
1128 
1129  }
1130 }
1131 
1132 /*
1133  * This function allocates histograms for matched truth-det jets in the case of embedding.
1134  */
1135 void AliAnalysisTaskEmcalJetPerformance::AllocateMatchedJetHistograms()
1136 {
1137  TString histname;
1138  TString title;
1139  Int_t nPtBins1 = TMath::CeilNint(fMaxPt-fMinPt);
1140  Int_t nPtBinsTruth2 = TMath::CeilNint(fMaxPt/2);
1141 
1142  // Response matrix, (centrality, pT-truth, pT-det)
1143  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
1144  Int_t nbinsy = fMaxPt; Int_t miny = 0; Int_t maxy = fMaxPt;
1145  Int_t nbinsz = nPtBins1; Int_t minz = fMinPt; Int_t maxz = fMaxPt;
1146 
1147  histname = "MatchedJetHistograms/hResponseMatrixEMCal";
1148  if (fForceBeamType == kAA) {
1149  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
1150  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1151  }
1152  else {
1153  title = histname + ";#it{p}_{T,corr}^{det} (GeV/#it{c});#it{p}_{T}^{truth} (GeV/#it{c})";
1154  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsy, miny, maxy);
1155  }
1156 
1157  if (fPlotDCal) {
1158  histname = "MatchedJetHistograms/hResponseMatrixDCal";
1159  if (fForceBeamType == kAA) {
1160  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
1161  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1162  }
1163  else {
1164  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
1165  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1166  }
1167  }
1168 
1169  // JES shift, (centrality, pT-truth, (pT-det - pT-truth) / pT-truth)
1170  nbinsx = 20; minx = 0; maxx = 100;
1171  nbinsy = nPtBinsTruth2; miny = 0; maxy = fMaxPt;
1172  nbinsz = 250; minz = -5.; maxz = 5.;
1173 
1174  histname = "MatchedJetHistograms/hJESshiftEMCal";
1175  if (fForceBeamType == kAA) {
1176  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1177  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1178  }
1179  else {
1180  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1181  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1182  }
1183 
1184  if (fPlotDCal) {
1185  histname = "MatchedJetHistograms/hJESshiftDCal";
1186  if (fForceBeamType == kAA) {
1187  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1188  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1189  }
1190  else {
1191  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
1192  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1193  }
1194  }
1195 
1196  // NEF of det-level matched jets, (centrality, pT-truth, NEF)
1197  nbinsx = 20; minx = 0; maxx = 100;
1198  nbinsy = fMaxPt; miny = 0; maxy = fMaxPt;
1199  nbinsz = 50; minz = 0; maxz = 1.;
1200 
1201  histname = "MatchedJetHistograms/hNEFVsPt";
1202  if (fForceBeamType == kAA) {
1203  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});Calo energy fraction";
1204  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1205  }
1206  else {
1207  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});Calo energy fraction";
1208  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1209  }
1210 
1211  // z-leading (charged) of det-level matched jets, (centrality, pT-truth, z-leading)
1212  nbinsx = 20; minx = 0; maxx = 100;
1213  nbinsy = nPtBinsTruth2; miny = 0; maxy = fMaxPt;
1214  nbinsz = 50; minz = 0; maxz = 1.;
1215 
1216  histname = "MatchedJetHistograms/hZLeadingVsPt";
1217  if (fForceBeamType == kAA) {
1218  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#it{z}_{leading}";
1219  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1220  }
1221  else {
1222  title = histname + ";#it{p}_{T}^{truth} (GeV/#it{c});#it{z}_{leading}";
1223  fHistManager.CreateTH2(histname.Data(), title.Data(), nbinsy, miny, maxy, nbinsz, minz, maxz);
1224  }
1225 
1226  if (fDoJetMatchingGeometrical) {
1227 
1228  // Matching distance, (pT-det, pT-truth, deltaR)
1229  nbinsx = nPtBins1; minx = fMinPt; maxx = fMaxPt;
1230  nbinsy = nPtBinsTruth2; miny = 0; maxy = fMaxPt;
1231  nbinsz = 100; minz = 0; maxz = 4.;
1232  histname = "MatchedJetHistograms/hMatchingDistance";
1233  title = histname + ";#it{p}_{T}^{det} (GeV/#it{c});#it{p}_{T}^{truth} (GeV/#it{c});R";
1234  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1235 
1236  // Jet matching QA (copied from AliAnalysisTaskEmcalJetHCorrelations.cxx)
1237  histname = "MatchedJetHistograms/fHistJetMatchingQA";
1238  title = histname;
1239  std::vector<std::string> binLabels = {"noMatch", "matchedJet", "uniqueMatch", "jetDistance", "passedAllCuts"};
1240  auto histMatchedJetCuts = fHistManager.CreateTH1(histname.Data(), title.Data(), binLabels.size(), 0, binLabels.size());
1241  // Set label names
1242  for (unsigned int i = 1; i <= binLabels.size(); i++) {
1243  histMatchedJetCuts->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
1244  }
1245  histMatchedJetCuts->GetYaxis()->SetTitle("Number of jets");
1246  }
1247 
1248  if (fDoJetMatchingMCFraction) {
1249 
1250  // (det pT, shared MC fraction, deltaR) of closest jets
1251  nbinsx = nPtBins1; minx = fMinPt; maxx = fMaxPt;
1252  nbinsy = 20; miny = 0; maxy = 1.;
1253  nbinsz = 100; minz = 0; maxz = 4.;
1254  histname = "MatchedJetHistograms/hMatchingDistanceVsMCFraction";
1255  title = histname + ";#it{p}_{T}^{det} (GeV/#it{c});MC fraction;#DeltaR";
1256  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
1257 
1258  // Jet matching QA (copied from AliAnalysisTaskEmcalJetHCorrelations.cxx)
1259  histname = "MatchedJetHistograms/fHistJetMatchingQA";
1260  title = histname;
1261  std::vector<std::string> binLabels = {"noMatch", "matchedJet", "sharedMomentumFraction", "partLevelMatchedJet", "jetDistancePPdet", "jetDistancePPtruth", "passedAllCuts"};
1262  auto histMatchedJetCuts = fHistManager.CreateTH1(histname.Data(), title.Data(), binLabels.size(), 0, binLabels.size());
1263  // Set label names
1264  for (unsigned int i = 1; i <= binLabels.size(); i++) {
1265  histMatchedJetCuts->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
1266  }
1267  histMatchedJetCuts->GetYaxis()->SetTitle("Number of jets");
1268  }
1269 
1270 }
1271 
1276 void AliAnalysisTaskEmcalJetPerformance::ExecOnce()
1277 {
1278  // Configure base class to set fTriggerPatchInfo to array of trigger patches, each event
1279  // (Need to call this before base class ExecOnce)
1280  if (fDoTriggerSimulation) {
1281  this->SetCaloTriggerPatchInfoName("EmcalTriggers");
1282  }
1283 
1284  AliAnalysisTaskEmcalJet::ExecOnce();
1285 
1286  fNeedEmcalGeom = kTRUE;
1287 
1288  // Check if trigger patches are loaded
1289  if (fDoTriggerSimulation) {
1290  if (fTriggerPatchInfo) {
1291  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1292  TClass cls(objname);
1293  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1294  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1295  GetName(), cls.GetName(), "EmcalTriggers"));
1296  fTriggerPatchInfo = 0;
1297  }
1298  }
1299  if (!fTriggerPatchInfo) {
1300  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1301  return;
1302  }
1303  }
1304 }
1305 
1309 void AliAnalysisTaskEmcalJetPerformance::RunChanged(Int_t run){
1310 
1311  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
1312 
1313  if (fPlotJetHistograms && fComputeMBDownscaling) {
1314 
1315  // Get instance of the downscale factor helper class
1316  PWG::EMCAL::AliEmcalDownscaleFactorsOCDB *downscaleOCDB = PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::Instance();
1317  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
1318 
1319  // There are two possible min bias triggers for LHC15o
1320  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
1321  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
1322  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
1323 
1324  // Get the downscale factor for whichever MB trigger exists in the given run
1325  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
1326  Double_t downscalefactor;
1327  for (auto i : runtriggers) {
1328  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
1329  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
1330  break;
1331  }
1332  }
1333 
1334  // Store the inverse of the downscale factor, used later to weight the pT spectrum
1335  fMBUpscaleFactor = 1/downscalefactor;
1336 
1337  TString histname = "Trigger/hMBDownscaleFactor";
1338  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
1339 
1340  }
1341 
1342 }
1343 
1347 Bool_t AliAnalysisTaskEmcalJetPerformance::IsEventSelected()
1348 {
1349  if (fUseAliEventCuts) {
1350  if (!fEventCuts.AcceptEvent(InputEvent()))
1351  {
1352  PostData(1, fOutput);
1353  return kFALSE;
1354  }
1355  }
1356  else {
1357  Bool_t answer = AliAnalysisTaskEmcal::IsEventSelected();
1358  return answer;
1359  }
1360  return kTRUE;
1361 }
1362 
1370 Bool_t AliAnalysisTaskEmcalJetPerformance::Run()
1371 {
1372  TString histname;
1373  AliJetContainer* jetCont = 0;
1374  TIter next(&fJetCollArray);
1375  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1376  TString jetContName = jetCont->GetName();
1377 
1378  // Do a simple trigger simulation (if requested)
1379  if (fDoTriggerSimulation) {
1380  DoTriggerSimulation();
1381  }
1382 
1383  }
1384 
1385  // Compute the full jet background scale factor and delta-pt
1386  if (fComputeBackground) {
1387  ComputeBackground();
1388  }
1389 
1390  // Only fill the embedding qa plots if:
1391  // - We are using the embedding helper
1392  // - The class has been initialized
1393  // - Both jet collections are available
1394  if (fEmbeddingQA.IsInitialized()) {
1395  fEmbeddingQA.RecordEmbeddedEventProperties();
1396  }
1397  return kTRUE;
1398 }
1399 
1403 void AliAnalysisTaskEmcalJetPerformance::DoTriggerSimulation()
1404 {
1405  TString histname;
1406 
1407  // Check if trigger patches are loaded
1408  if (fTriggerPatchInfo) {
1409  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1410  TClass cls(objname);
1411  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1412  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1413  GetName(), cls.GetName(), "EmcalTriggers"));
1414  fTriggerPatchInfo = 0;
1415  }
1416  }
1417  if (!fTriggerPatchInfo) {
1418  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1419  return;
1420  }
1421 
1422  // Compute patches in EMCal, DCal (I want offline simple trigger patch, i.e. patch calculated using FEE energy)
1423  std::vector<Double_t> vecEMCal;
1424  std::vector<Double_t> vecDCal;
1425  Double_t maxPatchEnergy = 0.;
1426  Double_t patchE;
1427  AliEMCALTriggerPatchInfo *maxPatch = nullptr;
1428  for(auto p : *fTriggerPatchInfo){
1429  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1430  if (recpatch) {
1431 
1432  if(!recpatch->IsJetHighSimple()) continue;
1433 
1434  patchE = recpatch->GetPatchE();
1435 
1436  histname = "TriggerSimHistograms/hEtaVsPhi";
1437  fHistManager.FillTH2(histname.Data(), recpatch->GetEtaGeo(), recpatch->GetPhiGeo());
1438 
1439  histname = "TriggerSimHistograms/hPatchE";
1440  fHistManager.FillTH2(histname.Data(), fCent, patchE);
1441 
1442  if (recpatch->IsEMCal()) {
1443  vecEMCal.push_back(patchE);
1444  if (patchE > maxPatchEnergy) {
1445  maxPatchEnergy = patchE;
1446  maxPatch = recpatch;
1447  }
1448  } else {
1449  vecDCal.push_back(patchE);
1450  }
1451 
1452  }
1453  }
1454 
1455  // Compute the median in each calorimeter
1456  const Int_t nBkgPatchesEMCal = vecEMCal.size(); // 6*8;
1457  const Int_t nBkgPatchesDCal = vecDCal.size(); // 4*5;
1458  fMedianEMCal = TMath::Median(nBkgPatchesEMCal, &vecEMCal[0]); // point to array used internally by vector
1459  fMedianDCal = TMath::Median(nBkgPatchesDCal, &vecDCal[0]);
1460 
1461  histname = "TriggerSimHistograms/hPatchMedianE";
1462  fHistManager.FillTH3(histname.Data(), fCent, fMedianEMCal, kEMCal);
1463  fHistManager.FillTH3(histname.Data(), fCent, fMedianDCal, kDCal);
1464 
1465  histname = "TriggerSimHistograms/hNPatches";
1466  fHistManager.FillTH2(histname.Data(), nBkgPatchesEMCal, kEMCal);
1467  fHistManager.FillTH2(histname.Data(), nBkgPatchesDCal, kDCal);
1468 
1469  // Fill max patch response and total EMCal energy response, if requested
1470  if (fGeneratorLevel && fDoTriggerResponse) {
1471 
1472  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
1473  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
1474  Double_t detEnergyEMCal = 0;
1475  Double_t truthEnergyEMCal = 0;
1476  Double_t truthEnergyPatch = 0;
1477 
1478  AliTLorentzVector part;
1479  Double_t partEta;
1480  Double_t partPhi;
1481  Double_t partE;
1482  for (auto mcpart:fGeneratorLevel->accepted_momentum()) {
1483  part.Clear();
1484  part = mcpart.first;
1485  partEta = part.Eta();
1486  partPhi = part.Phi_0_2pi();
1487  partE = part.E();
1488  if (maxPatch) {
1489  if (partPhi < maxPatch->GetPhiMax() && partPhi > maxPatch->GetPhiMin()) {
1490  if (partEta < maxPatch->GetEtaMax() && partEta > maxPatch->GetEtaMin()) {
1491  truthEnergyPatch += partE;
1492  }
1493  }
1494  }
1495  if (partPhi < phiMaxEMCal && partPhi > phiMinEMCal) {
1496  if (partEta < 0.7 && partEta > -0.7) {
1497  truthEnergyEMCal += partE;
1498  }
1499  }
1500  }
1501 
1502  AliTLorentzVector clusFourVec;
1503  const AliVCluster* clus;
1504  Double_t clusEta;
1505  Double_t clusPhi;
1506  AliClusterContainer* clusters = GetClusterContainer(0);
1507  for (auto clusIterator : clusters->accepted_momentum() ) {
1508  clusFourVec.Clear();
1509  clusFourVec = clusIterator.first;
1510  clusEta = clusFourVec.Eta();
1511  clusPhi = clusFourVec.Phi_0_2pi();
1512  clus = clusIterator.second;
1513  if (clus->IsEMCAL()) {
1514  if (clusPhi < phiMaxEMCal && clusPhi > phiMinEMCal) {
1515  if (clusEta < 0.7 && clusEta > -0.7) {
1516  detEnergyEMCal += clusFourVec.E();
1517  }
1518  }
1519  }
1520  }
1521 
1522  if (maxPatch) {
1523  histname = "TriggerSimHistograms/hMaxPatchResponseMatrix";
1524  fHistManager.FillTH2(histname.Data(), maxPatchEnergy, truthEnergyPatch);
1525  }
1526 
1527  histname = "TriggerSimHistograms/hEMCalEnergyResponseMatrix";
1528  fHistManager.FillTH2(histname.Data(), detEnergyEMCal, truthEnergyEMCal);
1529 
1530  }
1531 
1532  // Then compute background subtracted patches, by subtracting from each patch the median patch E from the opposite hemisphere
1533  // If a patch is above threshold, the event is "triggered"
1534  Bool_t fkEMCEJE = kFALSE;
1535  Double_t threshold = 20;
1536  for(auto p : *fTriggerPatchInfo){
1537  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1538  if (recpatch) {
1539 
1540  if(!recpatch->IsJetHighSimple()) continue;
1541 
1542  if (recpatch->IsEMCal()) {
1543  if ((recpatch->GetPatchE() - fMedianDCal) > threshold) {
1544  fkEMCEJE = kTRUE;
1545  break;
1546  }
1547  } else {
1548  if ((recpatch->GetPatchE() - fMedianEMCal) > threshold) {
1549  fkEMCEJE = kTRUE;
1550  break;
1551  }
1552  }
1553  }
1554  }
1555 
1556  if (fkEMCEJE) {
1557  FillTriggerSimHistograms();
1558  }
1559 
1560 }
1561 
1568 Bool_t AliAnalysisTaskEmcalJetPerformance::FillHistograms()
1569 {
1570 
1571  if (fPlotJetHistograms) {
1572  FillJetHistograms();
1573  }
1574  if (fPlotClusterHistograms) {
1575  FillClusterHistograms();
1576  }
1577  if (fPlotCellNonlinearityHistograms) {
1578  FillCellNonlinearityHistograms();
1579  }
1580  if (fPlotParticleCompositionHistograms) {
1581  FillParticleCompositionHistograms();
1582  }
1583  if (fPlotMatchedJetHistograms) {
1584  FillMatchedJetHistograms();
1585  }
1586 
1587  return kTRUE;
1588 }
1589 
1594 void AliAnalysisTaskEmcalJetPerformance::FillJetHistograms()
1595 {
1596  TString histname;
1597  AliJetContainer* jets = 0;
1598  TIter nextJetColl(&fJetCollArray);
1599  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1600  TString jetContName = jets->GetName();
1601 
1602  if (fPlotMatchedJetHistograms && fDoJetMatchingMCFraction && jets != fMCJetContainer) {
1603  if (!fDoClosureTest) {
1604  continue; // don't plot det-level histograms if embedding
1605  }
1606  }
1607 
1608  Double_t rhoVal = 0;
1609  if (jets->GetRhoParameter()) {
1610  rhoVal = jets->GetRhoVal();
1611  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
1612  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1613  }
1614 
1615  for (auto jet : jets->all()) {
1616 
1617  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1618  Float_t corrPt = GetJetPt(jet, rhoVal);
1619  Double_t jetR = jets->GetJetRadius();
1620 
1621  // compute jet acceptance type
1622  Double_t type = GetJetType(jet);
1623  if ( type != kEMCal ) {
1624  if ( type != kDCal || !fPlotDCal ) {
1625  continue;
1626  }
1627  }
1628 
1629  // (Centrality, Area, pT) (fill before area cut)
1630  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
1631  if (fForceBeamType == kAA) {
1632  fHistManager.FillTH3(histname.Data(), fCent, corrPt, jet->Area()/(TMath::Pi()*pow(jetR,2)));
1633  }
1634  else {
1635  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area()/(TMath::Pi()*pow(jetR,2)));
1636  }
1637 
1638  // (Centrality, pT-leading, pT) (before leading hadron cuts)
1639  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1640  if (fForceBeamType == kAA) {
1641  fHistManager.FillTH3(histname.Data(), fCent, corrPt, ptLeading);
1642  }
1643  else {
1644  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1645  }
1646 
1647  // (Centrality, pT, z-leading (charged)) (before leading hadron cuts)
1648  if (type == kEMCal) {
1649  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
1650  }
1651  else if (type == kDCal) {
1652  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
1653  }
1654  TLorentzVector leadPart;
1655  jets->GetLeadingHadronMomentum(leadPart, jet);
1656  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1657  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
1658  if (fForceBeamType == kAA) {
1659  fHistManager.FillTH3(histname, fCent, corrPt, z);
1660  }
1661  else {
1662  fHistManager.FillTH2(histname, corrPt, z);
1663  }
1664 
1665  // Rejection reason
1666  UInt_t rejectionReason = 0;
1667  if (!jets->AcceptJet(jet, rejectionReason)) {
1668  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1669  if (fForceBeamType == kAA) {
1670  fHistManager.FillTH3(histname.Data(), fCent, jets->GetRejectionReasonBitPosition(rejectionReason), corrPt);
1671  }
1672  else {
1673  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), corrPt);
1674  }
1675  continue;
1676  }
1677 
1678  // (Centrality, pT, NEF)
1679  if (type == kEMCal) {
1680  histname = TString::Format("%s/JetHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
1681  }
1682  else if (type == kDCal) {
1683  histname = TString::Format("%s/JetHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
1684  }
1685  if (fForceBeamType == kAA) {
1686  fHistManager.FillTH3(histname, fCent, corrPt, jet->NEF());
1687  }
1688  else {
1689  fHistManager.FillTH2(histname, corrPt, jet->NEF());
1690  }
1691 
1692  // (Centrality, pT upscaled, calo type)
1693  if (fComputeMBDownscaling) {
1694  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
1695  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type, fMBUpscaleFactor);
1696  }
1697 
1698  // (Centrality, pT, z (charged))
1699  if (type == kEMCal) {
1700  histname = TString::Format("%s/JetHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
1701  }
1702  else if (type == kDCal) {
1703  histname = TString::Format("%s/JetHistograms/hZVsPtDCal", jets->GetArrayName().Data());
1704  }
1705  const AliVTrack* track;
1706  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1707  track = static_cast<AliVTrack*>(jet->Track(i));
1708  z = track->Pt() / TMath::Abs(corrPt);
1709  if (fForceBeamType == kAA) {
1710  fHistManager.FillTH3(histname, fCent, corrPt, z);
1711  }
1712  else {
1713  fHistManager.FillTH2(histname, corrPt, z);
1714  }
1715  }
1716 
1717  // (Centrality, pT, Nconst)
1718  if (type == kEMCal) {
1719  histname = TString::Format("%s/JetHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
1720  }
1721  else if (type == kDCal) {
1722  histname = TString::Format("%s/JetHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
1723  }
1724  if (fForceBeamType == kAA) {
1725  fHistManager.FillTH3(histname, fCent, corrPt, 1.*jet->GetNumberOfConstituents());
1726  }
1727  else {
1728  fHistManager.FillTH2(histname, corrPt, 1.*jet->GetNumberOfConstituents());
1729  }
1730 
1731  // (Centrality, pT) for eta<0 and eta>0
1732  if (fForceBeamType == kAA) {
1733  if (type == kEMCal) {
1734  if (jet->Eta() > 0) {
1735  histname = TString::Format("%s/JetHistograms/hEtaPosVsPtEMCal", jets->GetArrayName().Data());
1736  fHistManager.FillTH2(histname, fCent, corrPt);
1737  }
1738  else if (jet->Eta() < 0) {
1739  histname = TString::Format("%s/JetHistograms/hEtaNegVsPtEMCal", jets->GetArrayName().Data());
1740  fHistManager.FillTH2(histname, fCent, corrPt);
1741  }
1742  }
1743  }
1744 
1745  // (Centrality, jet pT, Enonlincorr - Ehadcorr)
1746  if (fForceBeamType == kAA) {
1747  Double_t deltaEhadcorr = 0;
1748  const AliVCluster* clus = nullptr;
1749  Int_t nClusters = jet->GetNumberOfClusters();
1750  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1751  clus = jet->Cluster(iClus);
1752  deltaEhadcorr += (clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy());
1753  }
1754 
1755  histname = TString::Format("%s/JetHistograms/hDeltaEHadCorr", jets->GetArrayName().Data());
1756  fHistManager.FillTH3(histname, fCent, corrPt, deltaEhadcorr);
1757  }
1758 
1759 
1760  // (Median patch energy, calo type, jet pT, centrality)
1761  if (fDoTriggerSimulation) {
1762  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
1763  Double_t x[4] = {fMedianEMCal, kEMCal, corrPt, fCent};
1764  fHistManager.FillTHnSparse(histname, x);
1765  Double_t y[4] = {fMedianDCal, kDCal, corrPt, fCent};
1766  fHistManager.FillTHnSparse(histname, y);
1767  }
1768 
1769  } //jet loop
1770 
1771  }
1772 }
1773 
1774 /*
1775  * This function fills the histograms for the calorimeter performance study.
1776  */
1777 void AliAnalysisTaskEmcalJetPerformance::FillClusterHistograms()
1778 {
1779  TString histname;
1780 
1781  // Loop through clusters
1782  AliClusterContainer* clusters = GetClusterContainer(0);
1783  const AliVCluster* clus;
1784  for (auto it : clusters->accepted_momentum()) {
1785 
1786  clus = it.second;
1787  Double_t clusPhi = it.first.Phi_0_2pi();
1788  Double_t clusEta = it.first.Eta();
1789 
1790  // Include only EMCal/DCal clusters
1791  if (!clus->IsEMCAL()) {
1792  continue;
1793  }
1794 
1795  // Compute the sum of matched track momentum, and various track matching / hadronic corretion quantities
1796  Double_t trackPSum = 0;
1797  Int_t nTracksMatched = 0;
1798  const AliVTrack* track = nullptr;
1799  for (Int_t itrack=0; itrack < clus->GetNTracksMatched(); itrack++) {
1800  track = dynamic_cast<AliVTrack*>(clus->GetTrackMatched(itrack));
1801  if (!track) {
1802  continue;
1803  }
1804 
1805  Double_t trackPhi = TVector2::Phi_0_2pi(track->GetTrackPhiOnEMCal());
1806  Double_t trackEta = track->GetTrackEtaOnEMCal();
1807  Double_t deta = TMath::Abs(clusEta - trackEta);
1808  Double_t dphi = TMath::Abs(clusPhi - trackPhi);
1809 
1810  if (deta < fTrackMatchingDeltaEtaMax && dphi < fTrackMatchingDeltaPhiMax) {
1811  trackPSum += track->P();
1812  nTracksMatched++;
1813  }
1814  }
1815 
1816  Double_t EoverP = 0;
1817  if (trackPSum > 1e-3) {
1818  EoverP = clus->GetNonLinCorrEnergy() / trackPSum;
1819  }
1820 
1821  Double_t deltaE = clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy();
1822  Double_t Rcorr = 0;
1823  if (trackPSum > 1e-3) {
1824  Rcorr = deltaE / trackPSum;
1825  }
1826  Double_t RcorrClus = deltaE / clus->GetNonLinCorrEnergy();
1827 
1830 
1831  // Fill M02 distribution (centrality, Eclus nonlincorr, M02)
1832  histname = "ClusterHistograms/hM02";
1833  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1834 
1835  // Plot Ncell distribution for M02 > 0.4 or M02 < 0.4 (centrality, Eclus nonlincorr, Ncells)
1836  if (clus->GetM02() > 0.4) {
1837  histname = "ClusterHistograms/hNcellsM02G04";
1838  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetNCells());
1839  }
1840  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1841  histname = "ClusterHistograms/hNcellsM02L04";
1842  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetNCells());
1843  }
1844 
1847 
1848  // Plot matched track pT (centrality, Eclus nonlincorr, trackPsum)
1849  histname = "ClusterHistograms/hMatchedTrackPt";
1850  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1851 
1852  if (clus->GetM02() > 0.4) {
1853  histname = "ClusterHistograms/hMatchedTrackPtM02G04";
1854  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1855  }
1856  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1857  histname = "ClusterHistograms/hMatchedTrackPtM02L04";
1858  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1859  }
1860 
1861  // Plot number of matched tracks (centrality, Eclus nonlincorr, N matches)
1862  histname = "ClusterHistograms/hMatchedTrackN";
1863  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1864 
1865  if (clus->GetM02() > 0.4) {
1866  histname = "ClusterHistograms/hMatchedTrackNM02G04";
1867  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1868  }
1869  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1870  histname = "ClusterHistograms/hMatchedTrackNM02L04";
1871  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1872  }
1873 
1874  // Plot M02 distribution for clusters with matched tracks (centrality, Eclus nonlincorr, M02)
1875  histname = "ClusterHistograms/hM02Matched";
1876  if (nTracksMatched > 0) {
1877  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1878  }
1879 
1880  // Plot M02 distribution for clusters without matched tracks (centrality, Eclus nonlincorr, M02)
1881  histname = "ClusterHistograms/hM02Unmatched";
1882  if (nTracksMatched == 0) {
1883  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1884  }
1885 
1886  // Plot clus-track deltaEta if there is one matched track (deltaEta, Eclus, M02)
1887  if (nTracksMatched == 1) {
1888 
1889  const AliVTrack* track = dynamic_cast<AliVTrack*>(clus->GetTrackMatched(0));
1890  if (track) {
1891  Double_t trackEta = track->GetTrackEtaOnEMCal();
1892  Double_t deta = trackEta - clusEta;
1893 
1894  if (fCent > 0 && fCent < 10) {
1895  histname = "ClusterHistograms/hDeltaEtaCentral";
1896  fHistManager.FillTH3(histname.Data(), deta, clus->GetNonLinCorrEnergy(), clus->GetM02());
1897  }
1898 
1899  if (fCent > 50 && fCent < 90) {
1900  histname = "ClusterHistograms/hDeltaEtaPeripheral";
1901  fHistManager.FillTH3(histname.Data(), deta, clus->GetNonLinCorrEnergy(), clus->GetM02());
1902  }
1903  }
1904  }
1905 
1908 
1909  // Plot E/p vs. M02 for 0-10% and 50-90% (Eclus nonlincorr, Eclus nonlincorr / trackPsum, M02)
1910  if (fCent > 0 && fCent < 10) {
1911  histname = "ClusterHistograms/hEoverPM02Central";
1912  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), EoverP, clus->GetM02());
1913  }
1914  if (fCent > 50 && fCent < 90) {
1915  histname = "ClusterHistograms/hEoverPM02Peripheral";
1916  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), EoverP, clus->GetM02());
1917  }
1918 
1921 
1922  // Fill Rcorr distribution (centrality, trackPSum, Rcorr = (Enonlincorr - Ehadcorr) / trackPSum)
1923  histname = "ClusterHistograms/hRcorrVsCent";
1924  fHistManager.FillTH3(histname, fCent, trackPSum, Rcorr);
1925 
1926  // Fill Rcorr distribution for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr)
1927  if (fCent > 0 && fCent < 10) {
1928  histname = "ClusterHistograms/hRcorr0-10";
1929  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, Rcorr);
1930  }
1931 
1932  // Fill Rcorr distribution for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr)
1933  if (fCent > 50 && fCent < 90) {
1934  histname = "ClusterHistograms/hRcorr50-90";
1935  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, Rcorr);
1936  }
1937 
1938  // Fill also Rcorr-clus (centrality, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1939  histname = "ClusterHistograms/hRcorrClusVsCent";
1940  fHistManager.FillTH3(histname, fCent, trackPSum, RcorrClus);
1941 
1942  // Fill Rcorr-clus for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1943  if (fCent > 0 && fCent < 10) {
1944  histname = "ClusterHistograms/hRcorrClus0-10";
1945  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, RcorrClus);
1946  }
1947 
1948  // Fill Rcorr-clus for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1949  if (fCent > 50 && fCent < 90) {
1950  histname = "ClusterHistograms/hRcorrClus50-90";
1951  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, RcorrClus);
1952  }
1953 
1954  }
1955 
1956  // Fill total track multiplicity
1957  histname = "ClusterHistograms/hTrackMultiplicity";
1958  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1959  Int_t nTracks = trackCont->GetNAcceptedTracks();
1960  fHistManager.FillTH2(histname.Data(), nTracks, fCent);
1961 
1962 }
1963 
1964 void AliAnalysisTaskEmcalJetPerformance::FillCellNonlinearityHistograms()
1965 {
1966  // Get cells from event
1967  fCaloCells = InputEvent()->GetEMCALCells();
1968 
1969  // Loop through clusters
1970  AliClusterContainer* clusters = GetClusterContainer(0);
1971  AliTLorentzVector clusVec;
1972  const AliVCluster* clus;
1973  for (auto it : clusters->accepted_momentum()) {
1974 
1975  clusVec = it.first;
1976  clus = it.second;
1977 
1978  // Include only EMCal clusters
1979  if (!clus->IsEMCAL()) {
1980  continue;
1981  }
1982  if (clusVec.Phi_0_2pi() > 4.) {
1983  continue;
1984  }
1985 
1986  Double_t clusEnergy = clus->E();
1987  Double_t clusEnergyNonLinCorr = clus->GetNonLinCorrEnergy();
1988  Double_t correctionFactor = clusEnergyNonLinCorr/clusEnergy;
1989 
1990  // Plot cluster non-linearity scale factor
1991  TString histname = "CellNonlinearityHistograms/hClusterNonlinearity";
1992  fHistManager.FillTH2(histname.Data(), clusEnergy, correctionFactor);
1993 
1994  // Plot cell non-linearity scale factor for all cells
1995  histname = "CellNonlinearityHistograms/hCellNonlinearity";
1996  Double_t leadEcell = 0;
1997  for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1998 
1999  Int_t absId = clus->GetCellAbsId(iCell);
2000  Double_t eCell = fCaloCells->GetCellAmplitude(absId);
2001  if (eCell > leadEcell) {
2002  leadEcell = eCell;
2003  }
2004 
2005  fHistManager.FillTH2(histname.Data(), eCell, correctionFactor);
2006  }
2007 
2008  // Plot cell non-linearity scale factor for leading cell
2009  histname = "CellNonlinearityHistograms/hLeadingCellNonlinearity";
2010  fHistManager.FillTH2(histname.Data(), leadEcell, correctionFactor);
2011 
2012  }
2013 
2014 }
2015 
2016 /*
2017  * This function fills particle composition histograms for the calorimeter performance study in MC.
2018  */
2019 void AliAnalysisTaskEmcalJetPerformance::FillParticleCompositionHistograms()
2020 {
2021  // If MC, get the MC event
2022  const AliMCEvent* mcevent = nullptr;
2023  if (fGeneratorLevel) {
2024  mcevent = MCEvent();
2025  }
2026  else {
2027  return;
2028  }
2029 
2030  // Loop through clusters, and plot M02 for each particle type
2031  FillParticleCompositionClusterHistograms(mcevent);
2032 
2033  // Loop through jets, to fill various histograms
2034  if (fPlotJetHistograms) {
2035  FillParticleCompositionJetHistograms(mcevent);
2036  }
2037 
2038 }
2039 
2040 /*
2041  * Loop through clusters, and plot M02 for each particle type
2042  */
2043 void AliAnalysisTaskEmcalJetPerformance::FillParticleCompositionClusterHistograms(const AliMCEvent* mcevent)
2044 {
2045  TString histname;
2046  AliClusterContainer* clusters = GetClusterContainer(0);
2047  const AliVCluster* clus;
2048  std::vector<ContributorType> vecContributorTypes;
2049  std::vector<Int_t> vecContributorLabels;
2050  for (auto it : clusters->accepted_momentum()) {
2051 
2052  clus = it.second;
2053 
2054  // Include only EMCal/DCal clusters (reject PHOS clusters)
2055  if (!clus->IsEMCAL()) {
2056  continue;
2057  }
2058 
2059  // Loop through the cluster's contributors in order to classify its type
2060  ParticleType particleType = kNotDefined;
2061  ContributorType contributorType = kUndefined;
2062  const Int_t nLabels = clus->GetNLabels();
2063 
2064  // Create a vector to store the contributor types for PhysicalPrimary particles
2065  vecContributorTypes.clear();
2066  vecContributorLabels.clear();
2067  for (Int_t iLabel=0; iLabel<nLabels; iLabel++) {
2068 
2069  Int_t label = clus->GetLabels()[iLabel];
2070  if (TMath::Abs(label) > 0) { // if the particle has a truth-level match, the label is nonzero
2071  contributorType = GetContributorType(clus, mcevent, label);
2072  if (contributorType != kUndefined) {
2073  vecContributorTypes.push_back(contributorType);
2074  vecContributorLabels.push_back(label);
2075  }
2076  }
2077  }
2078 
2079  Int_t nLabelsPhysPrim = vecContributorTypes.size();
2080 
2081  if (nLabelsPhysPrim == 1) {
2082 
2083  contributorType = vecContributorTypes[0];
2084 
2085  if (contributorType == kPhoton) {
2086  particleType = kSinglePhoton;
2087  }
2088  else if (contributorType == kElectron) {
2089  particleType = kSingleElectron;
2090  }
2091  else if (contributorType == kChargedPion) {
2092  particleType = kSingleChargedPion;
2093  }
2094  else if (contributorType == kProton) {
2095  particleType = kSingleProton;
2096  }
2097  else if (contributorType == kAntiProton) {
2098  particleType = kSingleAntiProton;
2099  }
2100  else if (contributorType == kChargedKaon) {
2101  particleType = kSingleChargedKaon;
2102  }
2103  else if (contributorType == kK0L) {
2104  particleType = kSingleK0L;
2105  }
2106  else if (contributorType == kNeutron) {
2107  particleType = kSingleNeutron;
2108  }
2109  else if (contributorType == kAntiNeutron) {
2110  particleType = kSingleAntiNeutron;
2111  }
2112  else {
2113  particleType = kSingleOther;
2114  }
2115 
2116  }
2117  else if (nLabelsPhysPrim == 2) {
2118 
2119  // Get the contributor particle types
2120  ContributorType contributorType1 = vecContributorTypes[0];
2121  ContributorType contributorType2 = vecContributorTypes[1];
2122 
2123  // Get the fraction of cluster energy from each contributor
2124  //Double_t frac0 = clus->GetClusterMCEdepFraction(0);
2125  Double_t frac1 = clus->GetClusterMCEdepFraction(1);
2126 
2127  // Check whether the leading/subleading contributors are photons/hadrons
2128  Bool_t isHadron1 = IsHadron(contributorType1);
2129  Bool_t isHadron2 = IsHadron(contributorType2);
2130  Bool_t isPhoton1 = contributorType1 == kPhoton;
2131  Bool_t isPhoton2 = contributorType2 == kPhoton;
2132 
2133  if (isHadron1 && isHadron2) {
2134  particleType = kHadronHadron;
2135  }
2136  else if (isPhoton1 && isHadron2) {
2137  particleType = kPhotonHadron;
2138 
2139  // Plot cluster energy when subleading hadron is subtracted
2140  Double_t photonEnergy = clus->GetNonLinCorrEnergy() * (1 - frac1);
2141  histname = "ClusterHistogramsMC/hPhotonHadronPhotonEnergy";
2142  fHistManager.FillTH3(histname.Data(), fCent, clus->GetM02(), photonEnergy);
2143  }
2144  else if (isHadron1 && isPhoton2) {
2145  particleType = kHadronPhoton;
2146 
2147  // Plot cluster energy when subleading hadron is subtracted
2148  Double_t hadronEnergy = clus->GetNonLinCorrEnergy() * (1 - frac1);
2149  histname = "ClusterHistogramsMC/hHadronPhotonHadronEnergy";
2150  fHistManager.FillTH3(histname.Data(), fCent, clus->GetM02(), hadronEnergy);
2151  }
2152  else if (isPhoton1 && isPhoton2) {
2153 
2154  // By default, assume the two photons are not a merged pi0
2155  particleType = kPhotonPhotonOther;
2156 
2157  // Using the vector of accepted contributor labels, check whether the two photons are from the same pi0
2158  AliAODMCParticle *part1 = fGeneratorLevel->GetMCParticleWithLabel(vecContributorLabels[0]);
2159  AliAODMCParticle *part2 = fGeneratorLevel->GetMCParticleWithLabel(vecContributorLabels[1]);
2160  if (part1 && part2) {
2161  Int_t iMother1 = part1->GetMother();
2162  Int_t iMother2 = part2->GetMother();
2163  AliVParticle *mother1 = mcevent->GetTrack(iMother1);
2164  AliVParticle *mother2 = mcevent->GetTrack(iMother2);
2165 
2166  if (mother1 && mother2) {
2167  if ( (mother1->PdgCode() == 111) && (mother2->PdgCode() == 111) ) {
2168  if (iMother1 == iMother2) {
2169  particleType = kMergedPi0;
2170  }
2171  }
2172  }
2173  }
2174  }
2175  else {
2176  particleType = kTwoContributorsOther; // this includes partially contained conversion overlaps
2177  }
2178 
2179  }
2180  else if (nLabelsPhysPrim > 2) {
2181  particleType = kMoreThanTwoContributors;
2182  }
2183 
2184  // (M02, Eclus, part type)
2185  if (fCent > 0 && fCent < 10) {
2186  histname = "ClusterHistogramsMC/hM02VsParticleTypeCentral";
2187  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType);
2188  }
2189  if (fCent > 50 && fCent < 90) {
2190  histname = "ClusterHistogramsMC/hM02VsParticleTypePeripheral";
2191  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType);
2192  }
2193 
2194  }
2195 }
2196 
2197 /*
2198  * Loop through jets, to fill particle composition histograms.
2199  */
2200 void AliAnalysisTaskEmcalJetPerformance::FillParticleCompositionJetHistograms(const AliMCEvent* mcevent)
2201 {
2202  TString histname;
2203 
2204  const AliVCluster* clus;
2205  AliJetContainer* jets = GetJetContainer(0); // there is only a single, det-level jet finder here
2206  for (const auto jet : jets->accepted()) {
2207 
2208  Double_t jetPt = GetJetPt(jet, 0);
2209 
2210  // Keep track of hadronic calo energy in each jet
2211  Double_t hadCaloEnergyUnmatched = 0;
2212  Double_t hadCaloEnergyMatchedNonlincorr = 0;
2213  Double_t hadCaloEnergyMatchedHadCorr = 0;
2214 
2215  // Loop through clusters in each jet, to plot several histograms
2216  Int_t nClusters = jet->GetNumberOfClusters();
2217  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
2218 
2219  clus = jet->Cluster(iClus);
2220 
2221  // Get the particle type of the cluster
2222  ContributorType contributorType = kUndefined;
2223  Int_t label = TMath::Abs(clus->GetLabel());
2224  if (label > 0) {
2225  contributorType = GetContributorType(clus, mcevent, label);
2226  }
2227 
2228  // Plot M02 for each particle type
2229  histname = "JetPerformanceMC/hM02VsContributorTypeJets";
2230  Double_t x[4] = {clus->GetM02(), clus->GetNonLinCorrEnergy(), contributorType, jetPt};
2231  fHistManager.FillTHnSparse(histname, x);
2232 
2233  // If the cluster is a hadron, sum its energy to compute the jet's hadronic calo energy
2234  Bool_t isHadron = IsHadron(contributorType);
2235  if (isHadron) {
2236  Bool_t hasMatchedTrack = (clus->GetNTracksMatched() > 0);
2237  //Bool_t hasMatchedTrack = ((clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy()) > 1e-3);
2238  if (hasMatchedTrack) {
2239  hadCaloEnergyMatchedNonlincorr += clus->GetNonLinCorrEnergy();
2240  hadCaloEnergyMatchedHadCorr += clus->GetHadCorrEnergy();
2241  }
2242  else {
2243  hadCaloEnergyUnmatched += clus->GetNonLinCorrEnergy();
2244  }
2245  }
2246 
2247  }
2248 
2249  // Fill hadronic calo energy in each jet
2250 
2251  // (Jet pT, Summed energy of hadronic clusters without a matched track)
2252  histname = "JetPerformance/hHadCaloEnergyUnmatched";
2253  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyUnmatched);
2254 
2255  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction))
2256  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
2257  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedNonlincorr);
2258 
2259  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction))
2260  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
2261  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedHadCorr);
2262 
2263  // Loop through particle types, and plot jet composition for each particle type
2264  histname = "JetPerformanceMC/hJetComposition";
2265  for (Int_t type = 0; type < 8; type++) {
2266 
2267  ContributorType contributorType = kUndefined;
2268  Double_t nSum = 0;
2269  Double_t pTsum = 0;
2270 
2271  // Loop through clusters in jet, and add to sum if cluster matches particle type
2272  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
2273 
2274  clus = jet->Cluster(iClus);
2275 
2276  Int_t label = TMath::Abs(clus->GetLabel());
2277  if (label > 0) {
2278  contributorType = GetContributorType(clus, mcevent, label);
2279  }
2280 
2281  if (type == contributorType) {
2282  nSum++;
2283  pTsum += clus->GetNonLinCorrEnergy();
2284  }
2285  }
2286 
2287  // Fill jet composition histogram
2288  Double_t x[4] = {jetPt, 1.*type, nSum, pTsum};
2289  fHistManager.FillTHnSparse(histname, x);
2290 
2291  }
2292  }
2293 }
2294 
2298 void AliAnalysisTaskEmcalJetPerformance::ComputeBackground()
2299 {
2300  // Loop over tracks and clusters in order to:
2301  // (1) Compute scale factor for full jets
2302  // (2) Compute delta-pT for full jets, with the random cone method
2303 
2304  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
2305  Double_t etaTPC = 0.9;
2306  Double_t etaEMCal = 0.7;
2307  Double_t etaMinDCal = 0.22;
2308  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
2309  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
2310  Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
2311  Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
2312 
2313  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
2314  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
2315  Double_t accDCal = 2 * (etaEMCal - etaMinDCal) * (phiMaxDCal - phiMinDCal);
2316 
2317  // Loop over jet containers
2318  AliJetContainer* jetCont = 0;
2319  TIter nextJetColl(&fJetCollArray);
2320  while ((jetCont = static_cast<AliJetContainer*>(nextJetColl()))) {
2321  if(jetCont->GetNJets()<1){
2322  continue;//Don't enter if there are no jets in the container
2323  }
2324  Int_t maxJetIds[] = {-1, -1};
2325  Float_t maxJetPts[] = { 0, 0};
2326  //loop over the jets in the jet container to find the leading and subleading jet
2327  //alternative: AliEmcalJet* leadingJet =jetCont->GetLeadingJet(); but there is no GetSubLeadingJet() function
2328  for (Int_t ij = 0; ij < jetCont->GetNJets(); ++ij) {
2329  //AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
2330  AliEmcalJet *jet = static_cast<AliEmcalJet*>(jetCont->GetJet(ij));
2331  //if (!AcceptJet(jet)) continue;
2332  if (jet->Pt() > maxJetPts[0]) {
2333  maxJetPts[1] = maxJetPts[0];
2334  maxJetIds[1] = maxJetIds[0];
2335  maxJetPts[0] = jet->Pt();
2336  maxJetIds[0] = ij;
2337  }
2338  else if (jet->Pt() > maxJetPts[1]) {
2339  maxJetPts[1] = jet->Pt();
2340  maxJetIds[1] = ij;
2341  }
2342  }
2343 
2344  // Define fiducial acceptances, to be used to generate random cones, and for scale factor studies
2345  TRandom3* r = new TRandom3(0);
2346  Double_t jetR = jetCont->GetJetRadius();
2347  Double_t etaEMCalfid = etaEMCal - jetR;
2348  Double_t etaMinDCalfid = etaMinDCal + jetR;
2349  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
2350  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
2351  Double_t phiMinDCalfid = phiMinDCal + jetR;
2352  Double_t phiMaxDCalfid = phiMaxDCal - jetR;
2353  Double_t accEMCalfid = 2 * etaEMCalfid * (phiMaxEMCalfid - phiMinEMCalfid);
2354  Double_t accDCalfid = 2 * (etaEMCalfid - etaMinDCalfid) * (phiMaxDCalfid - phiMinDCalfid);
2355  if ( (etaEMCalfid - etaMinDCalfid) < 0) {
2356  accDCalfid = 0;
2357  }
2358 
2359  // Generate EMCal random cone eta-phi
2360  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
2361  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
2362 
2363  // Generate DCal random cone eta-phi
2364  Double_t etaDCalRC = r->Uniform(etaMinDCalfid, etaEMCalfid);
2365  Double_t sign = r->Uniform(-1., 1.);
2366  if (sign < 0) {
2367  etaDCalRC = -1*etaDCalRC;
2368  }
2369  Double_t phiDCalRC = r->Uniform(phiMinDCalfid, phiMaxDCalfid);
2370 
2371  // Initialize the various sums to 0
2372  Double_t trackPtSumTPC = 0;
2373  Double_t trackPtSumEMCal = 0;
2374  Double_t trackPtSumEMCalfid = 0;
2375  Double_t trackPtSumDCal = 0;
2376  Double_t trackPtSumDCalfid = 0;
2377  Double_t trackPtSumEMCalRC = 0;
2378  Double_t trackPtSumDCalRC = 0;
2379  Double_t clusPtSumEMCal = 0;
2380  Double_t clusPtSumEMCalfid = 0;
2381  Double_t clusPtSumDCal = 0;
2382  Double_t clusPtSumDCalfid = 0;
2383  Double_t clusPtSumEMCalRC = 0;
2384  Double_t clusPtSumDCalRC = 0;
2385  Bool_t overlap=0;
2386  // Loop over tracks. Sum the track pT:
2387  // (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
2388  // Note: Loops over all det-level track containers. For data there should be only one. For embedding, there should be signal+background tracks.
2389  AliParticleContainer * partCont = 0;
2390  AliTLorentzVector trackVec;
2391  AliVParticle* track;
2392  Double_t trackEta;
2393  Double_t trackPhi;
2394  Double_t trackPt;
2395  Double_t deltaR;
2396  Double_t trackID=-1;
2397  TIter nextPartCont(&fParticleCollArray);
2398  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
2399 
2400  TString partContName = partCont->GetName();
2401  if (!partContName.CompareTo("tracks") || !partContName.CompareTo("thermalparticles")) {
2402  for (auto trackIterator : partCont->accepted_momentum() ) {
2403 
2404  trackVec.Clear();
2405  trackVec = trackIterator.first;
2406  trackEta = trackVec.Eta();
2407  trackPhi = trackVec.Phi_0_2pi();
2408  trackPt = trackVec.Pt();
2409  //To get track ID and particle pointer
2410  track = trackIterator.second;
2411  TClonesArray* fTracksContArray = partCont->GetArray();
2412  trackID = fTracksContArray->IndexOf(track);
2413 
2414  // (1)
2415  if (TMath::Abs(trackEta) < etaTPC) {
2416  trackPtSumTPC += trackPt;
2417  }
2418 
2419  // (2)
2420  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
2421  trackPtSumEMCal += trackPt;
2422  }
2423 
2424  // (3)
2425  if (TMath::Abs(trackEta) < etaEMCalfid && trackPhi > phiMinEMCalfid && trackPhi < phiMaxEMCalfid) {
2426  trackPtSumEMCalfid += trackPt;
2427  }
2428 
2429  // (4)
2430  if (fPlotDCal) {
2431  if (TMath::Abs(trackEta) > etaMinDCal && TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinDCal && trackPhi < phiMaxDCal) {
2432  trackPtSumDCal += trackPt;
2433  }
2434  }
2435 
2436  // (5)
2437  if (fPlotDCal) {
2438  if (TMath::Abs(trackEta) > etaMinDCalfid && TMath::Abs(trackEta) < etaEMCalfid && trackPhi > phiMinDCalfid && trackPhi < phiMaxDCalfid) {
2439  trackPtSumDCalfid += trackPt;
2440  }
2441  }
2442 
2443  // (6)
2444  deltaR = GetDeltaR(&trackVec, etaEMCalRC, phiEMCalRC);
2445  if (deltaR < jetR) {
2446  trackPtSumEMCalRC += trackPt;
2447  //Check if there is an overlap with a leading/subleading signal jet
2448  //if so set the flag to 1 and not take this random cone into accound
2449  if(IsSignalJetOverlap(1,trackID,jetCont,maxJetIds))
2450  {
2451  overlap=1;
2452  }
2453  }
2454 
2455  // (7)
2456  if (fPlotDCal) {
2457  deltaR = GetDeltaR(&trackVec, etaDCalRC, phiDCalRC);
2458  if (deltaR < jetR) {
2459  trackPtSumDCalRC += trackPt;
2460  }
2461  }
2462  }
2463  }
2464  }
2465 
2466  // Loop over clusters, if the jet container is for full jets. Sum the cluster ET:
2467  // (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
2468  TString jetContName = jetCont->GetName();
2469  if (!jetContName.Contains("Charged")) {
2470  AliClusterContainer* clusCont = GetClusterContainer(0);
2471  AliTLorentzVector clusVec;
2472  AliVCluster* clus;
2473  Double_t clusEta;
2474  Double_t clusPhi;
2475  Double_t clusPt;
2476  Double_t clusID=-1;
2477  for (auto clusIterator : clusCont->accepted_momentum() ) {
2478 
2479  clusVec.Clear();
2480  clusVec = clusIterator.first;
2481  clusEta = clusVec.Eta();
2482  clusPhi = clusVec.Phi_0_2pi();
2483  clusPt = clusVec.Pt();
2484  clus = clusIterator.second;
2485  TClonesArray* fClusContArray = clusCont->GetArray();
2486  clusID = fClusContArray->IndexOf(clus);
2487 
2488  // (1)
2489  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
2490  clusPtSumEMCal += clusPt;
2491  }
2492 
2493  // (2)
2494  if (TMath::Abs(clusEta) < etaEMCalfid && clusPhi > phiMinEMCalfid && clusPhi < phiMaxEMCalfid) {
2495  clusPtSumEMCalfid += clusPt;
2496  }
2497 
2498  // (3)
2499  if (fPlotDCal) {
2500  if (TMath::Abs(clusEta) > etaMinDCal && TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinDCal && clusPhi < phiMaxDCal) {
2501  clusPtSumDCal += clusPt;
2502  }
2503  }
2504 
2505  // (4)
2506  if (fPlotDCal) {
2507  if (TMath::Abs(clusEta) > etaMinDCalfid && TMath::Abs(clusEta) < etaEMCalfid && clusPhi > phiMinDCalfid && clusPhi < phiMaxDCalfid) {
2508  clusPtSumDCalfid += clusPt;
2509  }
2510  }
2511 
2512  // (5)
2513  deltaR = GetDeltaR(&clusVec, etaEMCalRC, phiEMCalRC);
2514  if (deltaR < jetR) {
2515  clusPtSumEMCalRC += clusPt;
2516  //Check if there is an overlap with a leading/subleading signal jet
2517  //if so set the flag to 1 and not take this random cone into accound
2518  if(IsSignalJetOverlap(0,clusID,jetCont,maxJetIds))
2519  {
2520  overlap=1;
2521  }
2522  }
2523  // (6)
2524  if (fPlotDCal) {
2525  deltaR = GetDeltaR(&clusVec, etaDCalRC, phiDCalRC);
2526  if (deltaR < jetR) {
2527  clusPtSumDCalRC += clusPt;
2528  }
2529  }
2530  }
2531  }
2532 
2533  // Compute the scale factor, as a function of centrality, for (1) EMCal, (2) EMCalfid, (3) DCal, (4) DCalfid
2534  // (1)
2535  Double_t numerator = (trackPtSumEMCal + clusPtSumEMCal) / accEMCal;
2536  Double_t denominator = trackPtSumTPC / accTPC;
2537  Double_t scaleFactor = numerator / denominator;
2538  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
2539  fHistManager.FillTH2(histname, fCent, scaleFactor);
2540 
2541  // (2)
2542  if (accEMCalfid > 1e-3) {
2543  numerator = (trackPtSumEMCalfid + clusPtSumEMCalfid) / accEMCalfid;
2544  scaleFactor = numerator / denominator;
2545  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCalFid", jetCont->GetArrayName().Data());
2546  fHistManager.FillTH2(histname, fCent, scaleFactor);
2547  }
2548 
2549  // (3)
2550  if (fPlotDCal) {
2551  numerator = (trackPtSumDCal + clusPtSumDCal) / accDCal;
2552  scaleFactor = numerator / denominator;
2553  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCal", jetCont->GetArrayName().Data());
2554  fHistManager.FillTH2(histname, fCent, scaleFactor);
2555  }
2556 
2557  // (4)
2558  if (fPlotDCal) {
2559  if (accDCalfid > 1e-3) {
2560  numerator = (trackPtSumDCalfid + clusPtSumDCalfid) / accDCalfid;
2561  scaleFactor = numerator / denominator;
2562  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCalFid", jetCont->GetArrayName().Data());
2563  fHistManager.FillTH2(histname, fCent, scaleFactor);
2564  }
2565  }
2566 
2567  // Compute delta pT, as a function of centrality
2568  // EMCal acceptance only charged component
2569  Double_t rho = jetCont->GetRhoVal();
2570  Double_t deltaPt = trackPtSumEMCalRC - rho * TMath::Pi() * jetR * jetR;
2571  if(fFillChargedFluctuations){
2572  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCalCharged", jetCont->GetArrayName().Data());
2573  fHistManager.FillTH2(histname, fCent, deltaPt);
2574 
2575  // EMCal acceptance only charged component
2576  // excluding random cones that overlap with signal jets
2577  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCalChargedExcl", jetCont->GetArrayName().Data());
2578  if(overlap==0)fHistManager.FillTH2(histname, fCent, deltaPt);
2579  }
2580  // EMCal acceptance charged and neutral component
2581  deltaPt += clusPtSumEMCalRC;
2582  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
2583  fHistManager.FillTH2(histname, fCent, deltaPt);
2584 
2585  // EMCal acceptance charged and neutral component
2586  // excluding random cones that overlap with signal jets
2587  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCalExcl", jetCont->GetArrayName().Data());
2588  if(overlap==0)fHistManager.FillTH2(histname, fCent, deltaPt);
2589 
2590  // DCal
2591  if (fPlotDCal) {
2592  if (accDCalfid > 1e-3) {
2593  deltaPt = trackPtSumDCalRC + clusPtSumDCalRC - rho * TMath::Pi() * jetR * jetR;
2594  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtDCal", jetCont->GetArrayName().Data());
2595  fHistManager.FillTH2(histname, fCent, deltaPt);
2596  }
2597  }
2598 
2599  delete r;
2600 
2601  }
2602 
2603 }
2604 
2609 void AliAnalysisTaskEmcalJetPerformance::FillTriggerSimHistograms()
2610 {
2611  TString histname;
2612  AliJetContainer* jets = 0;
2613  TIter nextJetColl(&fJetCollArray);
2614  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
2615  TString jetContName = jets->GetName();
2616  Double_t jetR = jets->GetJetRadius();
2617 
2618  Double_t rhoVal = 0;
2619  if (jets->GetRhoParameter()) {
2620  rhoVal = jets->GetRhoVal();
2621  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
2622  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
2623  }
2624 
2625  for (auto jet : jets->all()) {
2626 
2627  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
2628  Float_t corrPt = GetJetPt(jet, rhoVal);
2629 
2630  // A vs. pT (fill before area cut)
2631  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
2632  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area()/(TMath::Pi()*pow(jetR,2)));
2633 
2634  // Rejection reason
2635  UInt_t rejectionReason = 0;
2636  if (!jets->AcceptJet(jet, rejectionReason)) {
2637  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
2638  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
2639  continue;
2640  }
2641 
2642  // compute jet acceptance type
2643  Double_t type = GetJetType(jet);
2644  if ( type != kEMCal ) {
2645  if ( type != kDCal || !fPlotDCal ) {
2646  continue;
2647  }
2648  }
2649 
2650  // (Centrality, pT, NEF, calo type)
2651  if (type == kEMCal) {
2652  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
2653  }
2654  else if (type == kDCal) {
2655  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
2656  }
2657  fHistManager.FillTH3(histname, fCent, corrPt, jet->NEF());
2658 
2659  // pT-leading vs. pT
2660  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
2661  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
2662 
2663  // (Centrality, pT, z-leading (charged), calo type)
2664  if (type == kEMCal) {
2665  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
2666  }
2667  else if (type == kDCal) {
2668  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
2669  }
2670  TLorentzVector leadPart;
2671  jets->GetLeadingHadronMomentum(leadPart, jet);
2672  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
2673  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
2674  fHistManager.FillTH3(histname, fCent, corrPt, z);
2675 
2676  // (Centrality, pT, z (charged), calo type)
2677  if (type == kEMCal) {
2678  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
2679  }
2680  else if (type == kDCal) {
2681  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtDCal", jets->GetArrayName().Data());
2682  }
2683  const AliVTrack* track;
2684  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
2685  track = static_cast<AliVTrack*>(jet->Track(i));
2686  z = track->Pt() / TMath::Abs(corrPt);
2687  fHistManager.FillTH3(histname, fCent, corrPt, z);
2688  }
2689 
2690  // (Centrality, pT, Nconst)
2691  if (type == kEMCal) {
2692  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
2693  }
2694  else if (type == kDCal) {
2695  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
2696  }
2697  fHistManager.FillTH3(histname, fCent, corrPt, 1.*jet->GetNumberOfConstituents());
2698 
2699  } //jet loop
2700 
2701  }
2702 }
2703 
2712 void AliAnalysisTaskEmcalJetPerformance::FillMatchedJetHistograms()
2713 {
2714  // Loop over all jets and fill the ClosestJet(), i.e. the matching candidate.
2715  // Note: Allow truth jets to be outside of EMCALfid or fail 5 GeV requirement, since these can still contribute accepted det-jets
2716  // (but for the jet reconstruction efficiency denominator, the criteria should be enforced).
2717  if (fDoJetMatchingGeometrical) {
2718  if (fRequireMatchedJetAccepted) {
2719  fMCJetContainer->SetJetAcceptanceType(AliEmcalJet::kTPC);
2720  ComputeJetMatches(fDetJetContainer, fMCJetContainer, kTRUE);
2721  fMCJetContainer->SetJetAcceptanceType(AliEmcalJet::kEMCALfid);
2722  }
2723  else {
2724  ComputeJetMatches(fDetJetContainer, fMCJetContainer, kFALSE);
2725  }
2726  }
2727  else if (fDoJetMatchingMCFraction) {
2728  // First match PbPb-det to pp-det
2729  ComputeJetMatches(fDetJetContainer, fDetJetContainerPPIntermediate, kTRUE);
2730 
2731  // Then match pp-det to pp-truth
2732  if (fRequireMatchedJetAccepted) { // Require pp-truth be accepted (i.e. leading track req), but still allow geometrical acceptance migration
2733  fMCJetContainer->SetJetAcceptanceType(AliEmcalJet::kTPC);
2734  ComputeJetMatches(fDetJetContainerPPIntermediate, fMCJetContainer, kTRUE);
2735  fMCJetContainer->SetJetAcceptanceType(AliEmcalJet::kEMCALfid);
2736  }
2737  else{ // Don't require pp-truth jet to be accepted
2738  ComputeJetMatches(fDetJetContainerPPIntermediate, fMCJetContainer, kFALSE);
2739  }
2740 
2741  }
2742 
2743  // Loop through accepted det-level jets, and retrieve matching candidate.
2744  // It match passes criteria (i.e. matching distance, uniqueness, MC fraction), fill matching histos.
2745  Double_t rhoVal = 0;
2746  if (fDetJetContainer->GetRhoParameter()) {
2747  rhoVal = fDetJetContainer->GetRhoVal();
2748  }
2749  for (auto jet : fDetJetContainer->accepted()) {
2750 
2751  Float_t detPt = GetJetPt(jet, rhoVal);
2752 
2753  // Get the matched part-level jet
2754  const AliEmcalJet* matchedPartLevelJet = GetMatchedPartLevelJet(jet, detPt);
2755  if (!matchedPartLevelJet) {
2756  continue;
2757  }
2758  Float_t truthPt = matchedPartLevelJet->Pt();
2759 
2760  // compute jet acceptance type
2761  Double_t type = GetJetType(jet);
2762  if ( type != kEMCal ) {
2763  if ( type != kDCal || !fPlotDCal ) {
2764  continue;
2765  }
2766  }
2767 
2768  // Fill response matrix (centrality, pT-truth, pT-det)
2769  TString histname;
2770  if (type == kEMCal) {
2771  histname = "MatchedJetHistograms/hResponseMatrixEMCal";
2772  }
2773  else if (type == kDCal) {
2774  histname = "MatchedJetHistograms/hResponseMatrixDCal";
2775  }
2776  if (fForceBeamType == kAA) {
2777  fHistManager.FillTH3(histname, fCent, truthPt, detPt);
2778  }
2779  else {
2780  fHistManager.FillTH2(histname, detPt, truthPt);
2781  }
2782 
2783  // Fill JES shift (centrality, pT-truth, (pT-det - pT-truth) / pT-truth)
2784  if (type == kEMCal) {
2785  histname = "MatchedJetHistograms/hJESshiftEMCal";
2786  }
2787  else if (type == kDCal) {
2788  histname = "MatchedJetHistograms/hJESshiftDCal";
2789  }
2790  if (fForceBeamType == kAA) {
2791  fHistManager.FillTH3(histname, fCent, truthPt, (detPt-truthPt)/truthPt );
2792  }
2793  else {
2794  fHistManager.FillTH2(histname, truthPt, (detPt-truthPt)/truthPt );
2795  }
2796 
2797  // Fill NEF of det-level matched jets (centrality, pT-truth, NEF)
2798  histname = "MatchedJetHistograms/hNEFVsPt";
2799  if (fForceBeamType == kAA) {
2800  fHistManager.FillTH3(histname, fCent, truthPt, jet->NEF());
2801  }
2802  else {
2803  fHistManager.FillTH2(histname, truthPt, jet->NEF());
2804  }
2805 
2806  // Fill z-leading (charged) of det-level matched jets (centrality, pT-truth, z-leading)
2807  histname = "MatchedJetHistograms/hZLeadingVsPt";
2808  TLorentzVector leadPart;
2809  fDetJetContainer->GetLeadingHadronMomentum(leadPart, jet);
2810  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
2811  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
2812  if (fForceBeamType == kAA) {
2813  fHistManager.FillTH3(histname, fCent, truthPt, z);
2814  }
2815  else {
2816  fHistManager.FillTH2(histname, truthPt, z);
2817  }
2818 
2819  } //jet loop
2820 }
2821 
2822 /*
2823  * Loop over jets of two specified jet collections, and fill the ClosestJet(), i.e. the matching candidate.
2824  * The first collection always uses the container acceptance criteria.
2825  * The second collection can be configured to use the container acceptance criteria or not.
2826  */
2827 void AliAnalysisTaskEmcalJetPerformance::ComputeJetMatches(AliJetContainer* jetCont1, AliJetContainer* jetCont2, Bool_t bUseJetCont2Acceptance) {
2828 
2829  for (auto jet1 : jetCont1->all()) {
2830  jet1->ResetMatching();
2831  }
2832  for (auto jet2 : jetCont2->all()) {
2833  jet2->ResetMatching();
2834  }
2835 
2836  for (auto jet1 : jetCont1->accepted()) {
2837 
2838  if (bUseJetCont2Acceptance) {
2839  for (auto jet2 : jetCont2->accepted()) {
2840  SetJetClosestCandidate(jet1, jet2);
2841  }
2842  }
2843  else {
2844  for (auto jet2 : jetCont2->all()) {
2845  if (jet2->Pt() < fMCJetMinMatchingPt) {
2846  continue;
2847  }
2848  SetJetClosestCandidate(jet1, jet2);
2849  }
2850  }
2851  }
2852 }
2853 
2854 /*
2855  * Given two jets, set them as closest if they are closer than the current closest jets.
2856  */
2857 void AliAnalysisTaskEmcalJetPerformance::SetJetClosestCandidate(AliEmcalJet* jet1, AliEmcalJet* jet2) {
2858 
2859  Double_t deltaR = jet1->DeltaR(jet2);
2860  if (deltaR > 0.) {
2861 
2862  if (deltaR < jet1->ClosestJetDistance()) {
2863  jet1->SetClosestJet(jet2, deltaR);
2864  }
2865 
2866  if (deltaR < jet2->ClosestJetDistance()) {
2867  jet2->SetClosestJet(jet1, deltaR);
2868  }
2869  }
2870 }
2871 
2872 /*
2873  * Return a pointer to the matched truth-level jet, if it passes the matching criteria.
2874  * For fDoJetMatchingGeometrical, this means (1) within R = fJetMatchingR, (2) unique match.
2875  * For fDoJetMatchingMCFraction, this means also shared MC fraction requirement. That is, if the jet (combined jet)
2876  * is matched to a pp det-level jet, which is matched to a pp truth-level jet -- it must satisfy
2877  * (1) The shared momentum fraction being larger than some minimum value fMinSharedMomentumFraction, and
2878  * (2) Their matched distance being below the max matching distance fMaxMatchedJetDistance
2879  *
2880  * @param[in] detJet det-level jet to be checked for a successful truth-level match.
2881  * @return Pointer to truth-level matched jet, if it exists. False otherwise.
2882 */
2883 const AliEmcalJet* AliAnalysisTaskEmcalJetPerformance::GetMatchedPartLevelJet(const AliEmcalJet* detJet, Double_t detJetPt) {
2884 
2885  // Track in histogram how many matches pass each distinct matching criteria
2886  TString histNameQA = "MatchedJetHistograms/fHistJetMatchingQA";
2887 
2888  // Geometrical matching case (pp, p-Pb)
2889  if (fDoJetMatchingGeometrical) {
2890 
2891  bool returnValue = false;
2892  const AliEmcalJet* partLevelJet = detJet->ClosestJet();
2893  if (partLevelJet) {
2894  fHistManager.FillTH1(histNameQA, "matchedJet");
2895  returnValue = true;
2896 
2897  // Check if match is unique
2898  if (partLevelJet->ClosestJet() != detJet) {
2899  returnValue = false;
2900  }
2901  else {
2902  fHistManager.FillTH1(histNameQA, "uniqueMatch");
2903 
2904  // Fill matching distance between unique matches, without imposing deltaR cut (centrality, pT-truth, R)
2905  TString histname = "MatchedJetHistograms/hMatchingDistance";
2906  fHistManager.FillTH3(histname, detJetPt, partLevelJet->Pt(), detJet->ClosestJetDistance());
2907  }
2908 
2909  // Check if the matching distance cut is passed
2910  double matchedJetDistance = detJet->ClosestJetDistance();
2911  if (matchedJetDistance > fJetMatchingR) {
2912  returnValue = false;
2913  }
2914  else {
2915  fHistManager.FillTH1(histNameQA, "jetDistance");
2916  }
2917 
2918  // Record all cuts passed
2919  if (returnValue == true) {
2920  fHistManager.FillTH1(histNameQA, "passedAllCuts");
2921  }
2922  }
2923  else {
2924  fHistManager.FillTH1(histNameQA, "noMatch");
2925  returnValue = false;
2926  }
2927 
2928  if (returnValue) {
2929  return partLevelJet;
2930  }
2931  }
2932 
2933  // Shared MC fraction matching case (Pb-Pb)
2934  else if (fDoJetMatchingMCFraction) { // This function is essentially copied from AliAnalysisTaskEmcalJetHCorrelations::CheckForMatchedJet
2935 
2936  bool returnValue = false;
2937  const AliEmcalJet* partLevelJet = nullptr;
2938 
2939  // First, check if combined jet has a pp det-level match assigned
2940  if (detJet->ClosestJet()) {
2941  fHistManager.FillTH1(histNameQA, "matchedJet");
2942  returnValue = true;
2943 
2944  // Check shared momentum fraction.
2945  double sharedFraction = fDetJetContainer->GetFractionSharedPt(detJet, nullptr);
2946  if (sharedFraction < fMinSharedMomentumFraction) {
2947  returnValue = false;
2948  }
2949  else {
2950  fHistManager.FillTH1(histNameQA, "sharedMomentumFraction");
2951  }
2952 
2953  // Check that the combined jet has a particle-level match
2954  AliEmcalJet * detLevelJetPP = detJet->ClosestJet();
2955  partLevelJet = detLevelJetPP->ClosestJet();
2956  if (!partLevelJet) {
2957  returnValue = false;
2958  }
2959  else {
2960  fHistManager.FillTH1(histNameQA, "partLevelMatchedJet");
2961  }
2962 
2963  // Check the matching distance between the combined and pp det-level jets
2964  double matchedJetDistance = detJet->ClosestJetDistance();
2965  if (matchedJetDistance > fJetMatchingR) {
2966  returnValue = false;
2967  }
2968  else {
2969  fHistManager.FillTH1(histNameQA, "jetDistancePPdet");
2970  }
2971 
2972  // Check the matching distance between the combined and pp truth-level jets
2973  if (partLevelJet) {
2974  Double_t deltaR = detJet->DeltaR(partLevelJet);
2975  if (deltaR > fJetMatchingR) {
2976  returnValue = false;
2977  }
2978  else {
2979  fHistManager.FillTH1(histNameQA, "jetDistancePPtruth");
2980  }
2981  }
2982 
2983  // Record all cuts passed
2984  if (returnValue == true) {
2985  fHistManager.FillTH1(histNameQA, "passedAllCuts");
2986  }
2987 
2988  // Fill (det pT, shared MC fraction, deltaR) of closest jets
2989  TString histname = "MatchedJetHistograms/hMatchingDistanceVsMCFraction";
2990  fHistManager.FillTH3(histname, detJetPt, sharedFraction, matchedJetDistance);
2991 
2992  }
2993  else {
2994  fHistManager.FillTH1(histNameQA, "noMatch");
2995  returnValue = false;
2996  }
2997 
2998  if (returnValue) {
2999  return partLevelJet;
3000  }
3001  }
3002  return 0;
3003 }
3004 
3005 /*
3006  * Compute the MC particle type of a given cluster contributor, using the MC particle container
3007  */
3008 AliAnalysisTaskEmcalJetPerformance::ContributorType AliAnalysisTaskEmcalJetPerformance::GetContributorType(const AliVCluster* clus, const AliMCEvent* mcevent, Int_t label)
3009 {
3010  ContributorType contributorType = kUndefined;
3011 
3012  AliAODMCParticle *part = fGeneratorLevel->GetMCParticleWithLabel(label);
3013  if (part) {
3014 
3015  TString histname = "ClusterHistogramsMC/hClusterRejectionReasonMC";
3016  UInt_t rejectionReason = 0;
3017  if (!fGeneratorLevel->AcceptMCParticle(part, rejectionReason)) {
3018  fHistManager.FillTH2(histname, fGeneratorLevel->GetRejectionReasonBitPosition(rejectionReason), clus->GetNonLinCorrEnergy());
3019  return contributorType;
3020  }
3021 
3022  if (part->GetGeneratorIndex() == 0) { // generator index in cocktail
3023 
3024  // select charged pions, protons, kaons, electrons, muons
3025  Int_t pdg = part->PdgCode();
3026 
3027  if (pdg == 22) { // gamma 22
3028  contributorType = kPhoton;
3029  }
3030  else if (TMath::Abs(pdg) == 211) { // pi+ 211 (abs value ensures both particles and antiparticles are included)
3031  contributorType = kChargedPion;
3032  }
3033  else if (pdg == 2212) { // proton 2212
3034  contributorType = kProton;
3035  }
3036  else if (pdg == -2212) {
3037  contributorType = kAntiProton;
3038  }
3039  else if (TMath::Abs(pdg) == 321) { // K+ 321
3040  contributorType = kChargedKaon;
3041  }
3042  else if (pdg == 130) { // K0L 130
3043  contributorType = kK0L;
3044  }
3045  else if (pdg == 2112) { // neutron 2112
3046  contributorType = kNeutron;
3047  }
3048  else if (pdg == -2112) {
3049  contributorType = kAntiNeutron;
3050  }
3051  else if (TMath::Abs(pdg) == 11) { // e- 11
3052  contributorType = kElectron;
3053  }
3054  else if (TMath::Abs(pdg) == 13) { // mu- 13
3055  contributorType = kMuon;
3056  }
3057  else {
3058  contributorType = kOther;
3059  }
3060  }
3061  }
3062  return contributorType;
3063 }
3064 
3068 Double_t AliAnalysisTaskEmcalJetPerformance::GetDeltaR(const AliTLorentzVector* part, Double_t etaRef, Double_t phiRef)
3069 {
3070  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
3071  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
3072  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
3073  return deltaR;
3074 }
3075 
3079 Double_t AliAnalysisTaskEmcalJetPerformance::GetJetType(const AliEmcalJet* jet)
3080 {
3081  UInt_t jetType = jet->GetJetAcceptanceType();
3082  Double_t type = -1;
3083  if (jetType & AliEmcalJet::kEMCALfid) {
3084  type = kEMCal;
3085  }
3086  else if (jetType & AliEmcalJet::kDCALonlyfid) {
3087  type = kDCal;
3088  }
3089 
3090  return type;
3091 }
3092 
3096 Double_t AliAnalysisTaskEmcalJetPerformance::GetJetPt(const AliEmcalJet* jet, Double_t rho)
3097 {
3098  Double_t pT = jet->Pt() - rho * jet->Area();
3099  return pT;
3100 }
3101 
3105 Bool_t AliAnalysisTaskEmcalJetPerformance::IsHadron(const ContributorType contributor)
3106 {
3107  return (contributor == kChargedPion) || (contributor == kProton) || (contributor == kAntiProton) || (contributor == kChargedKaon) || (contributor == kK0L) || (contributor == kNeutron) || (contributor == kAntiNeutron);
3108 }
3109 
3115 Bool_t AliAnalysisTaskEmcalJetPerformance::IsSignalJetOverlap(Bool_t isTrack, Int_t particleID, const AliJetContainer* jetCont, Int_t maxJetIds[])
3116 {
3117  //loop over the two leading anti-kT jets in the event
3118  for (Int_t i=0; i<2; i++)
3119  {
3120  if(maxJetIds[i]==-1)continue;
3121  AliEmcalJet* jet= jetCont->GetJet(maxJetIds[i]);
3122  Int_t ClusterConstituentID=-1;
3123  Int_t TrackConstituentID =-1;
3124  //check if the particle in the random cone is a constituent of a signal-Anti-KT jet
3125  if(isTrack)
3126  {
3127  TrackConstituentID = jet->ContainsTrack(particleID);
3128  }
3129  else
3130  {
3131  ClusterConstituentID = jet->ContainsCluster(particleID);
3132  }
3133  if(ClusterConstituentID>-1 || TrackConstituentID>-1)
3134  {
3135  return kTRUE;
3136  }
3137  }
3138  return kFALSE;
3139 }
3143 AliAnalysisTaskEmcalJetPerformance* AliAnalysisTaskEmcalJetPerformance::AddTaskEmcalJetPerformance(
3144  const char *ntracks,
3145  const char *nclusters,
3146  const char *nGenLev,
3147  const Double_t minTrPt,
3148  const Double_t minClPt,
3149  const char *suffix)
3150 {
3151  // Get the pointer to the existing analysis manager via the static access method.
3152  //==============================================================================
3153  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3154  if (!mgr)
3155  {
3156  ::Error("AddTaskEmcalJetPerformance", "No analysis manager to connect to.");
3157  return 0;
3158  }
3159 
3160  // Check the analysis type using the event handlers connected to the analysis manager.
3161  //==============================================================================
3162  AliVEventHandler* handler = mgr->GetInputEventHandler();
3163  if (!handler)
3164  {
3165  ::Error("AddTaskEmcalJetPerformance", "This task requires an input event handler");
3166  return 0;
3167  }
3168 
3169  enum EDataType_t {
3170  kUnknown,
3171  kESD,
3172  kAOD
3173  };
3174 
3175  EDataType_t dataType = kUnknown;
3176 
3177  if (handler->InheritsFrom("AliESDInputHandler")) {
3178  dataType = kESD;
3179  }
3180  else if (handler->InheritsFrom("AliAODInputHandler")) {
3181  dataType = kAOD;
3182  }
3183 
3184  //-------------------------------------------------------
3185  // Init the task and do settings
3186  //-------------------------------------------------------
3187 
3188  TString trackName(ntracks);
3189  TString clusName(nclusters);
3190 
3191  if (trackName == "usedefault") {
3192  if (dataType == kESD) {
3193  trackName = "Tracks";
3194  }
3195  else if (dataType == kAOD) {
3196  trackName = "tracks";
3197  }
3198  else {
3199  trackName = "";
3200  }
3201  }
3202 
3203  if (clusName == "usedefault") {
3204  if (dataType == kESD) {
3205  clusName = "CaloClusters";
3206  }
3207  else if (dataType == kAOD) {
3208  clusName = "caloClusters";
3209  }
3210  else {
3211  clusName = "";
3212  }
3213  }
3214 
3215  TString name("AliAnalysisTaskEmcalJetPerformance");
3216  if (!trackName.IsNull()) {
3217  name += "_";
3218  name += trackName;
3219  }
3220  if (!clusName.IsNull()) {
3221  name += "_";
3222  name += clusName;
3223  }
3224  if (strcmp(suffix,"") != 0) {
3225  name += "_";
3226  name += suffix;
3227  }
3228 
3230  // Configure jet performance task
3231  AliAnalysisTaskEmcalJetPerformance* task = new AliAnalysisTaskEmcalJetPerformance(name);
3232 
3234  // Create track and cluster containers with the standard cuts
3235 
3236  AliParticleContainer* partCont = 0;
3237  if (trackName == "mcparticles") {
3238  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
3239  partCont = mcpartCont;
3240  }
3241  else if (trackName == "tracks" || trackName == "Tracks") {
3242  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
3243  partCont = trackCont;
3244  }
3245  if (partCont) partCont->SetParticlePtCut(minTrPt);
3246  if (partCont) task->AdoptParticleContainer(partCont);
3247 
3248  // Add the generator-level container, if specified
3249  if (nGenLev && strcmp(nGenLev,"")!=0) {
3250  AliMCParticleContainer* mcpartCont = task->AddMCParticleContainer(nGenLev);
3251  mcpartCont->SelectPhysicalPrimaries(kTRUE);
3252  mcpartCont->SetParticlePtCut(0);
3253  }
3254 
3255  // Add the cluster container
3256  AliClusterContainer* clusCont = 0;
3257  if (!clusName.IsNull()) {
3258  clusCont = new AliClusterContainer(clusName);
3259  clusCont->SetClusECut(0.);
3260  clusCont->SetClusPtCut(0.);
3261  }
3262  if (clusCont) task->AdoptClusterContainer(clusCont);
3263 
3264  //-------------------------------------------------------
3265  // Final settings, pass to manager and set the containers
3266  //-------------------------------------------------------
3267 
3268  mgr->AddTask(task);
3269 
3270  // Create containers for input/output
3271  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
3272  TString contname(name);
3273  contname += "_histos";
3274  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
3275  TList::Class(),AliAnalysisManager::kOutputContainer,
3276  Form("%s", AliAnalysisManager::GetCommonFileName()));
3277  mgr->ConnectInput (task, 0, cinput1 );
3278  mgr->ConnectOutput (task, 1, coutput1 );
3279 
3280  return task;
3281 }
3282 
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
AliAnalysisTaskEmcalJetPerformance::AddTaskEmcalJetPerformance
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="")
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:3143
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:130
AliAnalysisTaskEmcal::fClusterCollArray
TObjArray fClusterCollArray
cluster collection array
Definition: AliAnalysisTaskEmcal.h:942
AliParticleContainer::SetParticlePtCut
void SetParticlePtCut(Double_t cut)
Definition: AliParticleContainer.h:58
AliJetContainer::GetRhoVal
Double_t GetRhoVal() const
Definition: AliJetContainer.h:165
AliJetContainer::GetRhoName
const TString & GetRhoName() const
Definition: AliJetContainer.h:166
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskEmcalJetPerformance::fEventCutList
TList * fEventCutList
! Output list for event cut histograms
Definition: AliAnalysisTaskEmcalJetPerformance.h:228
AliAnalysisTaskEmcalJetPerformance::GetJetPt
Double_t GetJetPt(const AliEmcalJet *jet, Double_t rho)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:3096
title
const char * title
Definition: MakeQAPdf.C:27
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliAnalysisTaskEmcalJetPerformance::fUseAliEventCuts
Bool_t fUseAliEventCuts
Flag to use AliEventCuts (otherwise AliAnalysisTaskEmcal will be used)
Definition: AliAnalysisTaskEmcalJetPerformance.h:226
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:377
TH3
Definition: External.C:244
AliEmcalJet::ClosestJetDistance
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:328
AliAnalysisTaskEmcal::fOffTrigger
UInt_t fOffTrigger
offline trigger for event selection
Definition: AliAnalysisTaskEmcal.h:919
AliAnalysisTaskEmcal::AdoptParticleContainer
void AdoptParticleContainer(AliParticleContainer *cont)
Definition: AliAnalysisTaskEmcal.h:214
AliClusterContainer::accepted_momentum
const AliClusterIterableMomentumContainer accepted_momentum() const
Definition: AliClusterContainer.cxx:511
AliMCParticleContainer::accepted_momentum
const AliMCParticleIterableMomentumContainer accepted_momentum() const
Definition: AliMCParticleContainer.cxx:242
nbinsy
Int_t nbinsy
Definition: CalculateAveragePt.C:62
AliAnalysisTaskEmcalJetPerformance::FillParticleCompositionJetHistograms
void FillParticleCompositionJetHistograms(const AliMCEvent *mcevent)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:2200
AliTrackContainer
Container with name, TClonesArray and cuts for particles.
Definition: AliTrackContainer.h:56
THistManager::FillTH2
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Definition: THistManager.cxx:458
AliTLorentzVector.h
Declaration of class AliTLorentzVector.
PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::Instance
static AliEmcalDownscaleFactorsOCDB * Instance()
Definition: AliEmcalDownscaleFactorsOCDB.cxx:53
AliAnalysisTaskEmcalJetPerformance::IsSignalJetOverlap
Bool_t IsSignalJetOverlap(Bool_t isTrack, Int_t particleID, const AliJetContainer *jet, Int_t maxJetIds[])
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:3115
THistManager::FillTH3
void FillTH3(const char *hname, double x, double y, double z, double weight=1., Option_t *opt="")
Fill a 3D histogram within the container.
Definition: THistManager.cxx:533
AliAnalysisTaskEmcalEmbeddingHelper.h
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
AliAnalysisTaskEmcalJetPerformance::fGeneratorLevel
AliMCParticleContainer * fGeneratorLevel
! generator level container
Definition: AliAnalysisTaskEmcalJetPerformance.h:232
AliJetContainer::AcceptJet
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Definition: AliJetContainer.cxx:492
PWG::EMCAL::AliEmcalDownscaleFactorsOCDB::GetDownscaleFactorForTriggerClass
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
Definition: AliEmcalDownscaleFactorsOCDB.cxx:78
AliParticleContainer
Container for particles within the EMCAL framework.
Definition: AliParticleContainer.h:31
AliAnalysisTaskEmcal::SetCaloTriggerPatchInfoName
void SetCaloTriggerPatchInfoName(const char *n)
Definition: AliAnalysisTaskEmcal.h:254
AliAnalysisTaskEmcalJetPerformance::fPlotClusterHistograms
Bool_t fPlotClusterHistograms
Set whether to plot cluster histograms.
Definition: AliAnalysisTaskEmcalJetPerformance.h:177
AliAnalysisTaskEmcalJetPerformance::GetDeltaR
Double_t GetDeltaR(const AliTLorentzVector *part, Double_t etaRef, Double_t phiRef)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:3068
AliEmcalJet::ContainsCluster
Int_t ContainsCluster(AliVCluster *cluster, TClonesArray *clusters) const
Definition: AliEmcalJet.cxx:761
AliAnalysisTaskEmcal::AdoptClusterContainer
void AdoptClusterContainer(AliClusterContainer *cont)
Definition: AliAnalysisTaskEmcal.h:217
AliAnalysisTaskEmcal::fParticleCollArray
TObjArray fParticleCollArray
particle/track collection array
Definition: AliAnalysisTaskEmcal.h:941
AliMCParticleContainer::AcceptMCParticle
virtual Bool_t AcceptMCParticle(const AliAODMCParticle *vp, UInt_t &rejectionReason) const
Definition: AliMCParticleContainer.cxx:155
AliEmcalEmbeddingQA::AddQAPlotsToList
bool AddQAPlotsToList(TList *list)
Definition: AliEmcalEmbeddingQA.cxx:75
AliAnalysisTaskEmcalJetPerformance::fDoJetMatchingGeometrical
Bool_t fDoJetMatchingGeometrical
Do geometrical matching between det-level and truth-level jet container.
Definition: AliAnalysisTaskEmcalJetPerformance.h:213
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1671
nPtBins
const Int_t nPtBins
Definition: ComputeDmesonYield.C:30
AliEmcalDownscaleFactorsOCDB.h
AliAnalysisTaskEmcalJetPerformance::fComputeMBDownscaling
Bool_t fComputeMBDownscaling
Set whether to compute and plot MB downscaling factors.
Definition: AliAnalysisTaskEmcalJetPerformance.h:183
AliJetContainer::GetLeadingHadronMomentum
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
Definition: AliJetContainer.cxx:653
AliAnalysisTaskEmcalJetPerformance::fEmbeddingQA
AliEmcalEmbeddingQA fEmbeddingQA
! QA hists for embedding (will only be added if embedding)
Definition: AliAnalysisTaskEmcalJetPerformance.h:215
THistManager::FillTHnSparse
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
Definition: THistManager.cxx:591
AliEmcalEmbeddingQA::IsInitialized
bool IsInitialized() const
Definition: AliEmcalEmbeddingQA.h:30
THistManager::CreateTH2
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.
Definition: THistManager.cxx:141
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskEmcalJetPerformance::fMCJetContainer
AliJetContainer * fMCJetContainer
!Pointer to jet container of truth-level jets
Definition: AliAnalysisTaskEmcalJetPerformance.h:216
AliAnalysisTaskEmcalJetPerformance::fTrackMatchingDeltaEtaMax
Double_t fTrackMatchingDeltaEtaMax
Maximum delta-eta to consider a track to be matched to a cluster.
Definition: AliAnalysisTaskEmcalJetPerformance.h:202
AliEmcalJet::kTPC
TPC acceptance.
Definition: AliEmcalJet.h:67
UInt_t
unsigned int UInt_t
Definition: External.C:33
THistManager::GetListOfHistograms
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
AliJetContainer::GetLeadingHadronPt
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
Definition: AliJetContainer.cxx:631
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskEmcalJetPerformance.h
Studies of full jets.
AliAnalysisTaskEmcal::fGeom
AliEMCALGeometry * fGeom
!emcal geometry
Definition: AliAnalysisTaskEmcal.h:961
AliEmcalJet::GetJetAcceptanceType
UInt_t GetJetAcceptanceType() const
Definition: AliEmcalJet.h:367
AliTLorentzVector::Phi_0_2pi
Double_t Phi_0_2pi() const
Definition: AliTLorentzVector.cxx:67
AliAnalysisTaskEmcalJetPerformance::fJetMatchingR
Double_t fJetMatchingR
Jet matching R threshold.
Definition: AliAnalysisTaskEmcalJetPerformance.h:220
AliParticleContainer.h
AliAnalysisTaskEmcalJetPerformance::ComputeJetMatches
void ComputeJetMatches(AliJetContainer *jetCont1, AliJetContainer *jetCont2, Bool_t bUseJetCont2Acceptance)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:2827
AliAnalysisTaskEmcalEmbeddingHelper
Implementation of task to embed external events.
Definition: AliAnalysisTaskEmcalEmbeddingHelper.h:71
AliJetContainer::GetNJets
Int_t GetNJets() const
Definition: AliJetContainer.h:156
THistManager::CreateTH1
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Definition: THistManager.cxx:80
AliAnalysisTaskEmcalJetPerformance::fTrackMatchingDeltaPhiMax
Double_t fTrackMatchingDeltaPhiMax
Maximum delta-phi to consider a track to be matched to a cluster.
Definition: AliAnalysisTaskEmcalJetPerformance.h:203
AliAnalysisTaskEmcalJetPerformance::kPhotonPhotonOther
Photon+Photon (not from same pi0) are the only contributors.
Definition: AliAnalysisTaskEmcalJetPerformance.h:95
AliAnalysisTaskEmcal::fForceBeamType
BeamType fForceBeamType
forced beam type
Definition: AliAnalysisTaskEmcal.h:900
AliAnalysisTaskEmcalJetPerformance::fUseManualEventCuts
Bool_t fUseManualEventCuts
Flag to use manual event cuts.
Definition: AliAnalysisTaskEmcalJetPerformance.h:229
AliAnalysisTaskEmcal::GetClusterContainer
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1678
AliJetContainer.h
AliEmcalJet::ContainsTrack
Int_t ContainsTrack(AliVParticle *track, TClonesArray *tracks) const
Definition: AliEmcalJet.cxx:718
AliAnalysisTaskEmcal::fCent
Double_t fCent
!event centrality
Definition: AliAnalysisTaskEmcal.h:967
AliAnalysisTaskEmcalJetPerformance::kSinglePhoton
Photon (direct or decay) is the only contributor.
Definition: AliAnalysisTaskEmcalJetPerformance.h:82
AliAnalysisTaskEmcalJetPerformance::kMergedPi0
Two particles from merged pi0 are the only contributors.
Definition: AliAnalysisTaskEmcalJetPerformance.h:94
AliRhoParameter.h
nullptr
#define nullptr
Definition: AliAnalysisTaskEmcalOfflinePatchesRef.cxx:35
AliAnalysisTaskEmcalJetPerformance::fRequireMatchedJetAccepted
Bool_t fRequireMatchedJetAccepted
Flag to require matched truth jet to be accepted (aside from geometrical acceptance) ...
Definition: AliAnalysisTaskEmcalJetPerformance.h:219
AliAnalysisTaskEmcalJetPerformance::fPlotJetHistograms
Bool_t fPlotJetHistograms
Set whether to enable inclusive jet histograms.
Definition: AliAnalysisTaskEmcalJetPerformance.h:176
AliAnalysisTaskEmcalJetPerformance::fNEoverPBins
Int_t fNEoverPBins
! number of variable E/p bins
Definition: AliAnalysisTaskEmcalJetPerformance.h:198
AliAnalysisTaskEmcal::AddMCParticleContainer
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
Definition: AliAnalysisTaskEmcal.cxx:1627
AliAnalysisTaskEmcal::GetParallelFraction
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...
Definition: AliAnalysisTaskEmcal.cxx:1911
THistManager::FillTH1
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Definition: THistManager.cxx:411
AliAnalysisTaskEmcalJetPerformance::fEventCuts
AliEventCuts fEventCuts
event selection utility
Definition: AliAnalysisTaskEmcalJetPerformance.h:227
AliAnalysisTaskEmcalJetPerformance::fDoJetMatchingMCFraction
Bool_t fDoJetMatchingMCFraction
Do MC-fraction based matching using PbPb det-level, pp det-level, and pp truth-level jet containers...
Definition: AliAnalysisTaskEmcalJetPerformance.h:214
AliAnalysisTaskEmcalJetPerformance::fPlotMatchedJetHistograms
Bool_t fPlotMatchedJetHistograms
Set whether to plot matched jet histograms.
Definition: AliAnalysisTaskEmcalJetPerformance.h:182
AliAnalysisTaskEmcalJet::fJetCollArray
TObjArray fJetCollArray
jet collection array
Definition: AliAnalysisTaskEmcalJet.h:97
AliEmcalJet::DeltaR
Double_t DeltaR(const AliVParticle *part) const
Definition: AliEmcalJet.cxx:401
AliAnalysisTaskEmcalJetPerformance::SetJetClosestCandidate
void SetJetClosestCandidate(AliEmcalJet *jet1, AliEmcalJet *jet2)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:2857
AliAnalysisTaskEmcalJetPerformance::fkEMCEJE
Bool_t fkEMCEJE
! flag telling whether the event is "triggered" or not in "simulation"
Definition: AliAnalysisTaskEmcalJetPerformance.h:209
AliAnalysisTaskEmcal::fCaloCells
AliVCaloCells * fCaloCells
!cells
Definition: AliAnalysisTaskEmcal.h:964
AliJetContainer::GetRhoParameter
AliRhoParameter * GetRhoParameter()
Definition: AliJetContainer.h:164
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliAnalysisTaskEmcalJetPerformance::GetMatchedPartLevelJet
const AliEmcalJet * GetMatchedPartLevelJet(const AliEmcalJet *jet, Double_t detJetPt)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:2883
AliAnalysisTaskEmcalJetPerformance::kPhotonHadron
Photon+Hadron are the only contributors, with Photon leading.
Definition: AliAnalysisTaskEmcalJetPerformance.h:92
nbinsx
Int_t nbinsx
Definition: CalculateAveragePt.C:61
AliAnalysisTaskEmcalJetPerformance::fMCJetMinMatchingPt
Double_t fMCJetMinMatchingPt
Min jet pT for MC jets being matched, for when container criteria is not applied. ...
Definition: AliAnalysisTaskEmcalJetPerformance.h:222
AliAnalysisTaskEmcal::IsEventSelected
virtual Bool_t IsEventSelected()
Performing event selection.
Definition: AliAnalysisTaskEmcal.cxx:1113
AliAnalysisTaskEmcalJetPerformance::IsHadron
Bool_t IsHadron(const ContributorType contributor)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:3105
AliAnalysisTaskEmcal::GenerateFixedBinArray
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
Definition: AliAnalysisTaskEmcal.h:1069
AliAnalysisTaskEmcalJetPerformance::fComputeBackground
Bool_t fComputeBackground
Set whether to enable study of background.
Definition: AliAnalysisTaskEmcalJetPerformance.h:180
AliJetContainer::GetJetRadius
Float_t GetJetRadius() const
Definition: AliJetContainer.h:174
AliAnalysisTaskEmcalJetPerformance::fMinSharedMomentumFraction
Double_t fMinSharedMomentumFraction
Minimum shared momentum (pp det-level track pT in combined jet) / (pp det-level track pT) ...
Definition: AliAnalysisTaskEmcalJetPerformance.h:221
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:990
AliAnalysisTaskEmcalJetPerformance::fDetJetContainerPPIntermediate
AliJetContainer * fDetJetContainerPPIntermediate
!Pointer to jet container of intermediate pp det-level jets, if MC-fraction matching ...
Definition: AliAnalysisTaskEmcalJetPerformance.h:218
AliAnalysisTaskEmcalJetPerformance::GetContributorType
ContributorType GetContributorType(const AliVCluster *clus, const AliMCEvent *mcevent, Int_t label)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:3008
PWG::EMCAL::AliEmcalDownscaleFactorsOCDB
Handler for downscale factors for various triggers obtained from the OCDB.
Definition: AliEmcalDownscaleFactorsOCDB.h:65
AliEmcalJet::SetClosestJet
void SetClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:323
AliAnalysisTaskEmcal::GetMCParticleContainer
AliMCParticleContainer * GetMCParticleContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcal.h:246
AliAnalysisTaskEmcalJetPerformance::fNPtHistBins
Int_t fNPtHistBins
! number of variable pt bins
Definition: AliAnalysisTaskEmcalJetPerformance.h:194
TH2
Definition: External.C:220
AliAnalysisTaskEmcalJetPerformance::fFillChargedFluctuations
Bool_t fFillChargedFluctuations
Set whether to fill also charged component of background fluctuations.
Definition: AliAnalysisTaskEmcalJetPerformance.h:186
AliAnalysisTaskEmcalJetPerformance::fPlotParticleCompositionHistograms
Bool_t fPlotParticleCompositionHistograms
Set whether to plot jet composition histograms.
Definition: AliAnalysisTaskEmcalJetPerformance.h:179
AliParticleContainer::accepted_momentum
const AliParticleIterableMomentumContainer accepted_momentum() const
Definition: AliParticleContainer.cxx:489
AliMCParticleContainer::SelectPhysicalPrimaries
void SelectPhysicalPrimaries(Bool_t s)
Definition: AliMCParticleContainer.h:51
AliAnalysisTaskEmcalJetPerformance::fDoClosureTest
Bool_t fDoClosureTest
Set whether to do thermal model closure test.
Definition: AliAnalysisTaskEmcalJetPerformance.h:185
AliAnalysisTaskEmcal::fTriggerPatchInfo
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Definition: AliAnalysisTaskEmcal.h:966
AliAnalysisTaskEmcalJetPerformance::fDoTriggerResponse
Bool_t fDoTriggerResponse
flag whether to compute max patch response, in case of MC
Definition: AliAnalysisTaskEmcalJetPerformance.h:210
AliAnalysisTaskEmcalJet
Base task in the EMCAL jet framework.
Definition: AliAnalysisTaskEmcalJet.h:30
AliEmcalJet
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
AliAnalysisTaskEmcalJetPerformance::fNM02HistBins
Int_t fNM02HistBins
! number of variable M02 bins
Definition: AliAnalysisTaskEmcalJetPerformance.h:196
dataType
TString dataType[]
Definition: CompareShowerShapeLongPerAnaPerSM.C:68
AliAnalysisTaskEmcalJetPerformance::fMBUpscaleFactor
Double_t fMBUpscaleFactor
! inverse of downscale factor, for MB trigger
Definition: AliAnalysisTaskEmcalJetPerformance.h:206
AliAnalysisTaskEmcalJetPerformance::fDetJetContainer
AliJetContainer * fDetJetContainer
!Pointer to jet container of det-level jets
Definition: AliAnalysisTaskEmcalJetPerformance.h:217
AliJetContainer::GetFractionSharedPt
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Definition: AliJetContainer.cxx:821
AliAnalysisTaskEmcalJetPerformance::fPlotCellNonlinearityHistograms
Bool_t fPlotCellNonlinearityHistograms
Set whether to plot cell non-linearity histograms for embedding studies.
Definition: AliAnalysisTaskEmcalJetPerformance.h:178
GetRunNumber
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
AliClusterContainer::SetClusPtCut
void SetClusPtCut(Double_t cut)
Definition: AliClusterContainer.h:54
AliAnalysisTaskEmcal::SetRejectionReasonLabels
void SetRejectionReasonLabels(TAxis *axis)
Definition: AliAnalysisTaskEmcal.cxx:1875
AliAnalysisTaskEmcal::UserCreateOutputObjects
void UserCreateOutputObjects()
Main initialization function on the worker.
Definition: AliAnalysisTaskEmcal.cxx:362
AliAnalysisTaskEmcalJetPerformance::fDoTriggerSimulation
Bool_t fDoTriggerSimulation
Set whether to perform a simple trigger simulation.
Definition: AliAnalysisTaskEmcalJetPerformance.h:181
AliJetContainer::SetJetAcceptanceType
void SetJetAcceptanceType(UInt_t type)
Definition: AliJetContainer.h:98
nbins
const Int_t nbins
Definition: AverageDmesonRaa.C:63
AliJetContainer::accepted
const AliJetIterableContainer accepted() const
Definition: AliJetContainer.cxx:966
Bool_t
bool Bool_t
Definition: External.C:53
suffix
TString suffix
Definition: ComputeEfficiencyFromCombinHF.C:24
AliAnalysisTaskEmcal::EDataType_t
EDataType_t
Switch for the data type.
Definition: AliAnalysisTaskEmcal.h:101
AliClusterContainer::SetClusECut
void SetClusECut(Double_t cut)
Definition: AliClusterContainer.h:53
AliAnalysisTaskEmcalJetPerformance::fMedianDCal
Double_t fMedianDCal
! median patch energy in DCal, per event
Definition: AliAnalysisTaskEmcalJetPerformance.h:208
THistManager::CreateTHnSparse
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.
Definition: THistManager.cxx:280
AliAnalysisTaskEmcalJetPerformance::fMedianEMCal
Double_t fMedianEMCal
! median patch energy in EMCal, per event
Definition: AliAnalysisTaskEmcalJetPerformance.h:207
AliClusterContainer
Container structure for EMCAL clusters.
Definition: AliClusterContainer.h:33
AliEmcalJet::kDCALonlyfid
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:74
AliMCParticleContainer
Container for MC-true particles within the EMCAL framework.
Definition: AliMCParticleContainer.h:24
AliMCParticleContainer::GetMCParticleWithLabel
virtual AliAODMCParticle * GetMCParticleWithLabel(Int_t lab) const
Definition: AliMCParticleContainer.cxx:48
AliAnalysisTaskEmcal::fNeedEmcalGeom
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Definition: AliAnalysisTaskEmcal.h:940
AliEmcalJet::kEMCALfid
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
AliEmcalJet.h
AliJetContainer
Container for jet within the EMCAL jet framework.
Definition: AliJetContainer.h:37
AliAnalysisTaskEmcalJetPerformance::kHadronPhoton
Hadron+Photon are the only contributors, with Hadron leading.
Definition: AliAnalysisTaskEmcalJetPerformance.h:93
AliAnalysisTaskEmcalJetPerformance::FillParticleCompositionClusterHistograms
void FillParticleCompositionClusterHistograms(const AliMCEvent *mcevent)
Definition: AliAnalysisTaskEmcalJetPerformance.cxx:2043
AliClusterContainer.h
THistManager::CreateTH3
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.
Definition: THistManager.cxx:208
AliAnalysisTaskEmcalJetPerformance::fPlotDCal
Bool_t fPlotDCal
Set whether to enable several DCal-specific histograms.
Definition: AliAnalysisTaskEmcalJetPerformance.h:184
AliAnalysisTaskEmcalJetPerformance::fPlotJetMatchCandThresh
Double_t fPlotJetMatchCandThresh
Threshold for jet R to count candidates, affects plotting only.
Definition: AliAnalysisTaskEmcalJetPerformance.h:223
AliJetContainer::GetJet
AliEmcalJet * GetJet(Int_t i) const
Definition: AliJetContainer.cxx:279
AliJetContainer::all
const AliJetIterableContainer all() const
Definition: AliJetContainer.cxx:957
AliAnalysisTaskEmcalEmbeddingHelper::GetInstance
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()
Definition: AliAnalysisTaskEmcalEmbeddingHelper.h:87