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