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