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