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