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