AliPhysics  b555aef (b555aef)
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 = nPtBins; 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  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
709  title = histname + ";Centrality;Scale factor;counts";
710  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
711 
712  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
713  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
714  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
715 
716  }
717 }
718 
719 /*
720  * This function allocates the histograms for single jets, when the "simulated" trigger has been fired.
721  * A set of histograms is allocated per each jet container.
722  */
724 {
725  TString histname;
726  TString title;
727  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
728 
729  //----------------------------------------------
730  // Trigger patch histograms
731 
732  // patch eta vs. phi
733  histname = "TriggerSimHistograms/hEtaVsPhi";
734  title = histname + ";#eta_{patch} (rad);#phi_{patch} (rad)";
735  fHistManager.CreateTH2(histname.Data(), title.Data(), 140, -0.7, 0.7, 500, 1., 6.);
736 
737  // N patches
738  histname = "TriggerSimHistograms/hNPatches";
739  title = histname + ";#it{N}_{patches};type";
740  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, 0, 200, 2, -0.5, 1.5);
741 
742  // patch E vs. centrality
743  histname = "TriggerSimHistograms/hPatchE";
744  title = histname + ";Centrality (%);#it{E}_{patch} (GeV)";
745  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, nPtBins, 0, fMaxPt);
746 
747  // patch median vs. Centrality
748  histname = "TriggerSimHistograms/hPatchMedianE";
749  title = histname + ";Centrality (%);#it{E}_{patch,med} (GeV);type";
750  fHistManager.CreateTH3(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 50, 2, -0.5, 1.5);
751 
752  //----------------------------------------------
753  // Jet histograms for "triggered" events
754  AliJetContainer* jets = 0;
755  TIter nextJetColl(&fJetCollArray);
756  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
757 
758  // Jet rejection reason
759  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
760  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
761  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
762  SetRejectionReasonLabels(hist->GetXaxis());
763 
764  // Rho vs. Centrality
765  if (!jets->GetRhoName().IsNull()) {
766  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
767  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
768  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
769  }
770 
771  // (Centrality, pT, NEF)
772  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
773  Int_t nbinsy = nPtBins; Int_t miny = 0; Int_t maxy = fMaxPt;
774  Int_t nbinsz = 50; Int_t minz = 0; Int_t maxz = 1.;
775 
776  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
777  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
778  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
779 
780  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
781  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
782  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
783 
784  // pT-leading vs. pT
785  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
786  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
787  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
788 
789  // A vs. pT
790  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
791  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
792  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/3, 0, 0.5);
793 
794  // (Centrality, pT, z-leading (charged))
795  nbinsx = 20; minx = 0; maxx = 100;
796  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
797  nbinsz = 50; minz = 0; maxz = 1.;
798 
799  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
800  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
801  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
802 
803  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
804  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
805  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
806 
807  // z (charged) vs. pT
808  nbinsx = 20; minx = 0; maxx = 100;
809  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
810  nbinsz = 50; minz = 0; maxz = 1.;
811 
812  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
813  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
814  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
815 
816  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtDCal", jets->GetArrayName().Data());
817  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
818  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
819 
820  // (Centrality, pT, Nconst)
821  nbinsx = 20; minx = 0; maxx = 100;
822  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
823  nbinsz = 50; minz = 0; maxz = fMaxPt;
824 
825  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
826  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
827  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
828 
829  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
830  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
831  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
832 
833  }
834 }
835 
836 /*
837  * This function allocates histograms for matched truth-det jets in the case of embedding.
838  * The jet matching information must be previously filled by another task, such as AliJetResponseMaker.
839  */
841 {
842  TString histname;
843  TString title;
844  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
845 
846  // Response matrix, (centrality, pT-truth, pT-det)
847  Int_t nbinsx = 20; Int_t minx = 0; Int_t maxx = 100;
848  Int_t nbinsy = nPtBins; Int_t miny = 0; Int_t maxy = fMaxPt;
849  Int_t nbinsz = nPtBins; Int_t minz = 0; Int_t maxz = fMaxPt;
850 
851  histname = "MatchedJetHistograms/hResponseMatrixEMCal";
852  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
853  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
854 
855  histname = "MatchedJetHistograms/hResponseMatrixDCal";
856  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#it{p}_{T,corr}^{det} (GeV/#it{c})";
857  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
858 
859  // JES shift, (centrality, pT-truth, (pT-det - pT-truth) / pT-truth)
860  nbinsx = 20; minx = 0; maxx = 100;
861  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
862  nbinsz = 250; minz = -5.; maxz = 5.;
863 
864  histname = "MatchedJetHistograms/hJESshiftEMCal";
865  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
866  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
867 
868  histname = "MatchedJetHistograms/hJESshiftDCal";
869  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});#frac{#it{p}_{T,corr}^{det} - #it{p}_{T}^{truth}}{#it{p}_{T}^{truth}}";
870  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
871 
872  // NEF of det-level matched jets, (centrality, pT-truth, NEF)
873  nbinsx = 20; minx = 0; maxx = 100;
874  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
875  nbinsz = 50; minz = 0; maxz = 1.;
876 
877  histname = "MatchedJetHistograms/hNEFVsPt";
878  title = histname + ";Centrality (%);#it{p}_{T,corr}^{det} (GeV/#it{c});Calo energy fraction";
879  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
880 
881  // z-leading (charged) of det-level matched jets, (centrality, pT-truth, z-leading)
882  nbinsx = 20; minx = 0; maxx = 100;
883  nbinsy = nPtBins; miny = 0; maxy = fMaxPt;
884  nbinsz = 50; minz = 0; maxz = 1.;
885 
886  histname = "MatchedJetHistograms/hZLeadingVsPt";
887  title = histname + ";Centrality (%);#it{p}_{T,corr}^{det} (GeV/#it{c});#it{z}_{leading}";
888  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
889 
890  // Matching distance, (centrality, pT-truth, R)
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/hMatchingDistance";
896  title = histname + ";Centrality (%);#it{p}_{T}^{truth} (GeV/#it{c});R";
897  fHistManager.CreateTH3(histname.Data(), title.Data(), nbinsx, minx, maxx, nbinsy, miny, maxy, nbinsz, minz, maxz);
898 
899 }
900 
906 {
907  // Configure base class to set fTriggerPatchInfo to array of trigger patches, each event
908  // (Need to call this before base class ExecOnce)
909  if (fDoTriggerSimulation) {
910  this->SetCaloTriggerPatchInfoName("EmcalTriggers");
911  }
912 
914 
915  fNeedEmcalGeom = kTRUE;
916 
917  // Check if trigger patches are loaded
918  if (fDoTriggerSimulation) {
919  if (fTriggerPatchInfo) {
920  TString objname(fTriggerPatchInfo->GetClass()->GetName());
921  TClass cls(objname);
922  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
923  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
924  GetName(), cls.GetName(), "EmcalTriggers"));
925  fTriggerPatchInfo = 0;
926  }
927  }
928  if (!fTriggerPatchInfo) {
929  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
930  return;
931  }
932  }
933 }
934 
939 
940  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
941 
943 
944  // Get instance of the downscale factor helper class
946  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
947 
948  // There are two possible min bias triggers for LHC15o
949  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
950  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
951  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
952 
953  // Get the downscale factor for whichever MB trigger exists in the given run
954  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
955  Double_t downscalefactor;
956  for (auto i : runtriggers) {
957  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
958  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
959  break;
960  }
961  }
962 
963  // Store the inverse of the downscale factor, used later to weight the pT spectrum
964  fMBUpscaleFactor = 1/downscalefactor;
965 
966  TString histname = "Trigger/hMBDownscaleFactor";
967  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
968 
969  }
970 
971 }
972 
977 {
978  if (fUseAliEventCuts) {
979  if (!fEventCuts.AcceptEvent(InputEvent()))
980  {
981  PostData(1, fOutput);
982  return kFALSE;
983  }
984  }
985  else {
987  }
988  return kTRUE;
989 }
990 
999 {
1000  TString histname;
1001  AliJetContainer* jetCont = 0;
1002  TIter next(&fJetCollArray);
1003  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1004  TString jetContName = jetCont->GetName();
1005 
1006  // Do a simple trigger simulation (if requested)
1007  if (fDoTriggerSimulation) {
1009  }
1010 
1011  }
1012 
1013  // Compute the full jet background scale factor and delta-pt
1014  if (fComputeBackground) {
1016  }
1017 
1018  // Only fill the embedding qa plots if:
1019  // - We are using the embedding helper
1020  // - The class has been initialized
1021  // - Both jet collections are available
1022  if (fEmbeddingQA.IsInitialized()) {
1024  }
1025 
1026  return kTRUE;
1027 }
1028 
1033 {
1034  TString histname;
1035 
1036  // Check if trigger patches are loaded
1037  if (fTriggerPatchInfo) {
1038  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1039  TClass cls(objname);
1040  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1041  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1042  GetName(), cls.GetName(), "EmcalTriggers"));
1043  fTriggerPatchInfo = 0;
1044  }
1045  }
1046  if (!fTriggerPatchInfo) {
1047  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1048  return;
1049  }
1050 
1051  // Compute patches in EMCal, DCal (I want offline simple trigger patch, i.e. patch calculated using FEE energy)
1052  std::vector<Double_t> vecEMCal;
1053  std::vector<Double_t> vecDCal;
1054  for(auto p : *fTriggerPatchInfo){
1055  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1056  if (recpatch) {
1057 
1058  if(!recpatch->IsJetHighSimple()) continue;
1059 
1060  histname = "TriggerSimHistograms/hEtaVsPhi";
1061  fHistManager.FillTH2(histname.Data(), recpatch->GetEtaGeo(), recpatch->GetPhiGeo());
1062 
1063  histname = "TriggerSimHistograms/hPatchE";
1064  fHistManager.FillTH2(histname.Data(), fCent, recpatch->GetPatchE());
1065 
1066  if (recpatch->IsEMCal()) {
1067  vecEMCal.push_back(recpatch->GetPatchE());
1068  } else {
1069  vecDCal.push_back(recpatch->GetPatchE());
1070  }
1071 
1072  }
1073  }
1074 
1075  // Compute the median in each calorimeter
1076  const Int_t nBkgPatchesEMCal = vecEMCal.size(); // 6*8;
1077  const Int_t nBkgPatchesDCal = vecDCal.size(); // 4*5;
1078  fMedianEMCal = TMath::Median(nBkgPatchesEMCal, &vecEMCal[0]); // point to array used internally by vector
1079  fMedianDCal = TMath::Median(nBkgPatchesDCal, &vecDCal[0]);
1080 
1081  histname = "TriggerSimHistograms/hPatchMedianE";
1082  fHistManager.FillTH3(histname.Data(), fCent, fMedianEMCal, kEMCal);
1083  fHistManager.FillTH3(histname.Data(), fCent, fMedianDCal, kDCal);
1084 
1085  histname = "TriggerSimHistograms/hNPatches";
1086  fHistManager.FillTH2(histname.Data(), nBkgPatchesEMCal, kEMCal);
1087  fHistManager.FillTH2(histname.Data(), nBkgPatchesDCal, kDCal);
1088 
1089  // Then compute background subtracted patches, by subtracting from each patch the median patch E from the opposite hemisphere
1090  // If a patch is above threshold, the event is "triggered"
1091  Bool_t fkEMCEJE = kFALSE;
1092  Double_t threshold = 20;
1093  for(auto p : *fTriggerPatchInfo){
1094  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1095  if (recpatch) {
1096 
1097  if(!recpatch->IsJetHighSimple()) continue;
1098 
1099  if (recpatch->IsEMCal()) {
1100  if ((recpatch->GetPatchE() - fMedianDCal) > threshold) {
1101  fkEMCEJE = kTRUE;
1102  break;
1103  }
1104  } else {
1105  if ((recpatch->GetPatchE() - fMedianEMCal) > threshold) {
1106  fkEMCEJE = kTRUE;
1107  break;
1108  }
1109  }
1110  }
1111  }
1112 
1113  if (fkEMCEJE) {
1115  }
1116 
1117 }
1118 
1126 {
1127 
1128  if (fPlotJetHistograms) {
1130  }
1131  if (fPlotClusterHistograms) {
1133  }
1136  }
1139  }
1140 
1141  return kTRUE;
1142 }
1143 
1149 {
1150  TString histname;
1151  AliJetContainer* jets = 0;
1152  TIter nextJetColl(&fJetCollArray);
1153  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1154  TString jetContName = jets->GetName();
1155 
1156  Double_t rhoVal = 0;
1157  if (jets->GetRhoParameter()) {
1158  rhoVal = jets->GetRhoVal();
1159  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
1160  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1161  }
1162 
1163  for (auto jet : jets->all()) {
1164 
1165  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1166  Float_t corrPt = GetJetPt(jet, rhoVal);
1167 
1168  // A vs. pT (fill before area cut)
1169  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
1170  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
1171 
1172 
1173  // Rejection reason
1174  UInt_t rejectionReason = 0;
1175  if (!jets->AcceptJet(jet, rejectionReason)) {
1176  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1177  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1178  continue;
1179  }
1180 
1181  // compute jet acceptance type
1182  Double_t type = GetJetType(jet);
1183  if ( (type != kEMCal) && (type != kDCal) ) {
1184  continue;
1185  }
1186 
1187  // (Centrality, pT, NEF)
1188  if (type == kEMCal) {
1189  histname = TString::Format("%s/JetHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
1190  }
1191  else if (type == kDCal) {
1192  histname = TString::Format("%s/JetHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
1193  }
1194  fHistManager.FillTH3(histname, fCent, corrPt, jet->NEF());
1195 
1196  // (Centrality, pT upscaled, calo type)
1197  if (fComputeMBDownscaling) {
1198  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
1199  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type, fMBUpscaleFactor);
1200  }
1201 
1202  // pT-leading vs. pT
1203  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1204  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1205 
1206  // (Centrality, pT, z-leading (charged))
1207  if (type == kEMCal) {
1208  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
1209  }
1210  else if (type == kDCal) {
1211  histname = TString::Format("%s/JetHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
1212  }
1213  TLorentzVector leadPart;
1214  jets->GetLeadingHadronMomentum(leadPart, jet);
1215  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1216  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
1217  fHistManager.FillTH3(histname, fCent, corrPt, z);
1218 
1219  // (Centrality, pT, z (charged))
1220  if (type == kEMCal) {
1221  histname = TString::Format("%s/JetHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
1222  }
1223  else if (type == kDCal) {
1224  histname = TString::Format("%s/JetHistograms/hZVsPtDCal", jets->GetArrayName().Data());
1225  }
1226  const AliVTrack* track;
1227  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1228  track = static_cast<AliVTrack*>(jet->Track(i));
1229  z = track->Pt() / TMath::Abs(corrPt);
1230  fHistManager.FillTH3(histname, fCent, corrPt, z);
1231  }
1232 
1233  // (Centrality, pT, Nconst)
1234  if (type == kEMCal) {
1235  histname = TString::Format("%s/JetHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
1236  }
1237  else if (type == kDCal) {
1238  histname = TString::Format("%s/JetHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
1239  }
1240  fHistManager.FillTH3(histname, fCent, corrPt, 1.*jet->GetNumberOfConstituents());
1241 
1242  // (Centrality, jet pT, Enonlincorr - Ehadcorr)
1243  Double_t deltaEhadcorr = 0;
1244  const AliVCluster* clus = nullptr;
1245  Int_t nClusters = jet->GetNumberOfClusters();
1246  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1247  clus = jet->Cluster(iClus);
1248  deltaEhadcorr += (clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy());
1249  }
1250 
1251  histname = TString::Format("%s/JetHistograms/hDeltaEHadCorr", jets->GetArrayName().Data());
1252  fHistManager.FillTH3(histname, fCent, corrPt, deltaEhadcorr);
1253 
1254 
1255  // (Median patch energy, calo type, jet pT, centrality)
1256  if (fDoTriggerSimulation) {
1257  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
1258  Double_t x[4] = {fMedianEMCal, kEMCal, corrPt, fCent};
1259  fHistManager.FillTHnSparse(histname, x);
1260  Double_t y[4] = {fMedianDCal, kDCal, corrPt, fCent};
1261  fHistManager.FillTHnSparse(histname, y);
1262  }
1263 
1264  } //jet loop
1265 
1266  }
1267 }
1268 
1269 /*
1270  * This function fills the histograms for the calorimeter performance study.
1271  */
1273 {
1274  TString histname;
1275 
1276  // Loop through clusters
1278  const AliVCluster* clus;
1279  for (auto it : clusters->accepted_momentum()) {
1280 
1281  clus = it.second;
1282  Double_t clusPhi = it.first.Phi_0_2pi();
1283  Double_t clusEta = it.first.Eta();
1284 
1285  // Include only EMCal/DCal clusters
1286  if (!clus->IsEMCAL()) {
1287  continue;
1288  }
1289 
1290  // Compute the sum of matched track momentum, and various track matching / hadronic corretion quantities
1291  Double_t trackPSum = 0;
1292  Int_t nTracksMatched = 0;
1293  const AliVTrack* track = nullptr;
1294  for (Int_t itrack=0; itrack < clus->GetNTracksMatched(); itrack++) {
1295  track = dynamic_cast<AliVTrack*>(clus->GetTrackMatched(itrack));
1296  if (!track) {
1297  continue;
1298  }
1299 
1300  Double_t trackPhi = TVector2::Phi_0_2pi(track->GetTrackPhiOnEMCal());
1301  Double_t trackEta = track->GetTrackEtaOnEMCal();
1302  Double_t deta = TMath::Abs(clusEta - trackEta);
1303  Double_t dphi = TMath::Abs(clusPhi - trackPhi);
1304 
1306  trackPSum += track->P();
1307  nTracksMatched++;
1308  }
1309  }
1310 
1311  Double_t EoverP = 0;
1312  if (trackPSum > 1e-3) {
1313  EoverP = clus->GetNonLinCorrEnergy() / trackPSum;
1314  }
1315 
1316  Double_t deltaE = clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy();
1317  Double_t Rcorr = 0;
1318  if (trackPSum > 1e-3) {
1319  Rcorr = deltaE / trackPSum;
1320  }
1321  Double_t RcorrClus = deltaE / clus->GetNonLinCorrEnergy();
1322 
1325 
1326  // Fill M02 distribution (centrality, Eclus nonlincorr, M02)
1327  histname = "ClusterHistograms/hM02";
1328  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1329 
1330  // Plot Ncell distribution for M02 > 0.4 or M02 < 0.4 (centrality, Eclus nonlincorr, Ncells)
1331  if (clus->GetM02() > 0.4) {
1332  histname = "ClusterHistograms/hNcellsM02G04";
1333  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetNCells());
1334  }
1335  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1336  histname = "ClusterHistograms/hNcellsM02L04";
1337  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetNCells());
1338  }
1339 
1342 
1343  // Plot matched track pT (centrality, Eclus nonlincorr, trackPsum)
1344  histname = "ClusterHistograms/hMatchedTrackPt";
1345  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1346 
1347  if (clus->GetM02() > 0.4) {
1348  histname = "ClusterHistograms/hMatchedTrackPtM02G04";
1349  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1350  }
1351  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1352  histname = "ClusterHistograms/hMatchedTrackPtM02L04";
1353  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), trackPSum);
1354  }
1355 
1356  // Plot number of matched tracks (centrality, Eclus nonlincorr, N matches)
1357  histname = "ClusterHistograms/hMatchedTrackN";
1358  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1359 
1360  if (clus->GetM02() > 0.4) {
1361  histname = "ClusterHistograms/hMatchedTrackNM02G04";
1362  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1363  }
1364  if (clus->GetM02() > 0.1 && clus->GetM02() < 0.4) {
1365  histname = "ClusterHistograms/hMatchedTrackNM02L04";
1366  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), nTracksMatched);
1367  }
1368 
1369  // Plot M02 distribution for clusters with matched tracks (centrality, Eclus nonlincorr, M02)
1370  histname = "ClusterHistograms/hM02Matched";
1371  if (nTracksMatched > 0) {
1372  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1373  }
1374 
1375  // Plot M02 distribution for clusters without matched tracks (centrality, Eclus nonlincorr, M02)
1376  histname = "ClusterHistograms/hM02Unmatched";
1377  if (nTracksMatched == 0) {
1378  fHistManager.FillTH3(histname, fCent, clus->GetNonLinCorrEnergy(), clus->GetM02());
1379  }
1380 
1381  // Plot clus-track deltaEta if there is one matched track (deltaEta, Eclus, M02)
1382  if (nTracksMatched == 1) {
1383 
1384  const AliVTrack* track = dynamic_cast<AliVTrack*>(clus->GetTrackMatched(0));
1385  if (track) {
1386  Double_t trackEta = track->GetTrackEtaOnEMCal();
1387  Double_t deta = trackEta - clusEta;
1388 
1389  if (fCent > 0 && fCent < 10) {
1390  histname = "ClusterHistograms/hDeltaEtaCentral";
1391  fHistManager.FillTH3(histname.Data(), deta, clus->GetNonLinCorrEnergy(), clus->GetM02());
1392  }
1393 
1394  if (fCent > 50 && fCent < 90) {
1395  histname = "ClusterHistograms/hDeltaEtaPeripheral";
1396  fHistManager.FillTH3(histname.Data(), deta, clus->GetNonLinCorrEnergy(), clus->GetM02());
1397  }
1398  }
1399  }
1400 
1403 
1404  // Plot E/p vs. M02 for 0-10% and 50-90% (Eclus nonlincorr, Eclus nonlincorr / trackPsum, M02)
1405  if (fCent > 0 && fCent < 10) {
1406  histname = "ClusterHistograms/hEoverPM02Central";
1407  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), EoverP, clus->GetM02());
1408  }
1409  if (fCent > 50 && fCent < 90) {
1410  histname = "ClusterHistograms/hEoverPM02Peripheral";
1411  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), EoverP, clus->GetM02());
1412  }
1413 
1416 
1417  // Fill Rcorr distribution (centrality, trackPSum, Rcorr = (Enonlincorr - Ehadcorr) / trackPSum)
1418  histname = "ClusterHistograms/hRcorrVsCent";
1419  fHistManager.FillTH3(histname, fCent, trackPSum, Rcorr);
1420 
1421  // Fill Rcorr distribution for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr)
1422  if (fCent > 0 && fCent < 10) {
1423  histname = "ClusterHistograms/hRcorr0-10";
1424  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, Rcorr);
1425  }
1426 
1427  // Fill Rcorr distribution for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr)
1428  if (fCent > 50 && fCent < 90) {
1429  histname = "ClusterHistograms/hRcorr50-90";
1430  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, Rcorr);
1431  }
1432 
1433  // Fill also Rcorr-clus (centrality, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1434  histname = "ClusterHistograms/hRcorrClusVsCent";
1435  fHistManager.FillTH3(histname, fCent, trackPSum, RcorrClus);
1436 
1437  // Fill Rcorr-clus for 0-10% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1438  if (fCent > 0 && fCent < 10) {
1439  histname = "ClusterHistograms/hRcorrClus0-10";
1440  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, RcorrClus);
1441  }
1442 
1443  // Fill Rcorr-clus for 50-90% centrality (Eclus nonlincorr, trackPSum, Rcorr-clus = (Enonlincorr - Ehadcorr) / Enonlincorr )
1444  if (fCent > 50 && fCent < 90) {
1445  histname = "ClusterHistograms/hRcorrClus50-90";
1446  fHistManager.FillTH3(histname, clus->GetNonLinCorrEnergy(), trackPSum, RcorrClus);
1447  }
1448 
1449  }
1450 
1451  // Fill total track multiplicity
1452  histname = "ClusterHistograms/hTrackMultiplicity";
1453  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1454  Int_t nTracks = trackCont->GetNAcceptedTracks();
1455  fHistManager.FillTH2(histname.Data(), nTracks, fCent);
1456 
1457 }
1458 
1459 /*
1460  * This function fills particle composition histograms for the calorimeter performance study in MC.
1461  */
1463 {
1464  // If MC, get the MC event
1465  const AliMCEvent* mcevent = nullptr;
1466  if (fGeneratorLevel) {
1467  mcevent = MCEvent();
1468  }
1469  else {
1470  return;
1471  }
1472 
1473  // Loop through clusters, and plot M02 for each particle type
1475 
1476  // Loop through jets, to fill various histograms
1477  if (fPlotJetHistograms) {
1479  }
1480 
1481 }
1482 
1483 /*
1484  * Loop through clusters, and plot M02 for each particle type
1485  */
1487 {
1488  TString histname;
1490  const AliVCluster* clus;
1491  std::vector<ContributorType> vecContributorTypes;
1492  std::vector<Int_t> vecContributorLabels;
1493  for (auto it : clusters->accepted_momentum()) {
1494 
1495  clus = it.second;
1496 
1497  // Include only EMCal/DCal clusters (reject PHOS clusters)
1498  if (!clus->IsEMCAL()) {
1499  continue;
1500  }
1501 
1502  // Loop through the cluster's contributors in order to classify its type
1503  ParticleType particleType = kNotDefined;
1504  ContributorType contributorType = kUndefined;
1505  const Int_t nLabels = clus->GetNLabels();
1506 
1507  // Create a vector to store the contributor types for PhysicalPrimary particles
1508  vecContributorTypes.clear();
1509  vecContributorLabels.clear();
1510  for (Int_t iLabel=0; iLabel<nLabels; iLabel++) {
1511 
1512  Int_t label = clus->GetLabels()[iLabel];
1513  if (TMath::Abs(label) > 0) { // if the particle has a truth-level match, the label is nonzero
1514  contributorType = GetContributorType(clus, mcevent, label);
1515  if (contributorType != kUndefined) {
1516  vecContributorTypes.push_back(contributorType);
1517  vecContributorLabels.push_back(label);
1518  }
1519  }
1520  }
1521 
1522  Int_t nLabelsPhysPrim = vecContributorTypes.size();
1523 
1524  if (nLabelsPhysPrim == 1) {
1525 
1526  contributorType = vecContributorTypes[0];
1527 
1528  if (contributorType == kPhoton) {
1529  particleType = kSinglePhoton;
1530  }
1531  else if (contributorType == kElectron) {
1532  particleType = kSingleElectron;
1533  }
1534  else if (contributorType == kChargedPion) {
1535  particleType = kSingleChargedPion;
1536  }
1537  else if (contributorType == kProton) {
1538  particleType = kSingleProton;
1539  }
1540  else if (contributorType == kAntiProton) {
1541  particleType = kSingleAntiProton;
1542  }
1543  else if (contributorType == kChargedKaon) {
1544  particleType = kSingleChargedKaon;
1545  }
1546  else if (contributorType == kK0L) {
1547  particleType = kSingleK0L;
1548  }
1549  else if (contributorType == kNeutron) {
1550  particleType = kSingleNeutron;
1551  }
1552  else if (contributorType == kAntiNeutron) {
1553  particleType = kSingleAntiNeutron;
1554  }
1555  else {
1556  particleType = kSingleOther;
1557  }
1558 
1559  }
1560  else if (nLabelsPhysPrim == 2) {
1561 
1562  // Get the contributor particle types
1563  ContributorType contributorType1 = vecContributorTypes[0];
1564  ContributorType contributorType2 = vecContributorTypes[1];
1565 
1566  // Get the fraction of cluster energy from each contributor
1567  //Double_t frac0 = clus->GetClusterMCEdepFraction(0);
1568  Double_t frac1 = clus->GetClusterMCEdepFraction(1);
1569 
1570  // Check whether the leading/subleading contributors are photons/hadrons
1571  Bool_t isHadron1 = IsHadron(contributorType1);
1572  Bool_t isHadron2 = IsHadron(contributorType2);
1573  Bool_t isPhoton1 = contributorType1 == kPhoton;
1574  Bool_t isPhoton2 = contributorType2 == kPhoton;
1575 
1576  if (isHadron1 && isHadron2) {
1577  particleType = kHadronHadron;
1578  }
1579  else if (isPhoton1 && isHadron2) {
1580  particleType = kPhotonHadron;
1581 
1582  // Plot cluster energy when subleading hadron is subtracted
1583  Double_t photonEnergy = clus->GetNonLinCorrEnergy() * (1 - frac1);
1584  histname = "ClusterHistogramsMC/hPhotonHadronPhotonEnergy";
1585  fHistManager.FillTH3(histname.Data(), fCent, clus->GetM02(), photonEnergy);
1586  }
1587  else if (isHadron1 && isPhoton2) {
1588  particleType = kHadronPhoton;
1589 
1590  // Plot cluster energy when subleading hadron is subtracted
1591  Double_t hadronEnergy = clus->GetNonLinCorrEnergy() * (1 - frac1);
1592  histname = "ClusterHistogramsMC/hHadronPhotonHadronEnergy";
1593  fHistManager.FillTH3(histname.Data(), fCent, clus->GetM02(), hadronEnergy);
1594  }
1595  else if (isPhoton1 && isPhoton2) {
1596 
1597  // By default, assume the two photons are not a merged pi0
1598  particleType = kPhotonPhotonOther;
1599 
1600  // Using the vector of accepted contributor labels, check whether the two photons are from the same pi0
1601  AliAODMCParticle *part1 = fGeneratorLevel->GetMCParticleWithLabel(vecContributorLabels[0]);
1602  AliAODMCParticle *part2 = fGeneratorLevel->GetMCParticleWithLabel(vecContributorLabels[1]);
1603  if (part1 && part2) {
1604  Int_t iMother1 = part1->GetMother();
1605  Int_t iMother2 = part2->GetMother();
1606  AliVParticle *mother1 = mcevent->GetTrack(iMother1);
1607  AliVParticle *mother2 = mcevent->GetTrack(iMother2);
1608 
1609  if (mother1 && mother2) {
1610  if ( (mother1->PdgCode() == 111) && (mother2->PdgCode() == 111) ) {
1611  if (iMother1 == iMother2) {
1612  particleType = kMergedPi0;
1613  }
1614  }
1615  }
1616  }
1617  }
1618  else {
1619  particleType = kTwoContributorsOther; // this includes partially contained conversion overlaps
1620  }
1621 
1622  }
1623  else if (nLabelsPhysPrim > 2) {
1624  particleType = kMoreThanTwoContributors;
1625  }
1626 
1627  // (M02, Eclus, part type)
1628  if (fCent > 0 && fCent < 10) {
1629  histname = "ClusterHistogramsMC/hM02VsParticleTypeCentral";
1630  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType);
1631  }
1632  if (fCent > 50 && fCent < 90) {
1633  histname = "ClusterHistogramsMC/hM02VsParticleTypePeripheral";
1634  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType);
1635  }
1636 
1637  }
1638 }
1639 
1640 /*
1641  * Loop through jets, to fill particle composition histograms.
1642  */
1644 {
1645  TString histname;
1646 
1647  const AliVCluster* clus;
1648  AliJetContainer* jets = GetJetContainer(0); // there is only a single, det-level jet finder here
1649  for (const auto jet : jets->accepted()) {
1650 
1651  Double_t jetPt = GetJetPt(jet, 0);
1652 
1653  // Keep track of hadronic calo energy in each jet
1654  Double_t hadCaloEnergyUnmatched = 0;
1655  Double_t hadCaloEnergyMatchedNonlincorr = 0;
1656  Double_t hadCaloEnergyMatchedHadCorr = 0;
1657 
1658  // Loop through clusters in each jet, to plot several histograms
1659  Int_t nClusters = jet->GetNumberOfClusters();
1660  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1661 
1662  clus = jet->Cluster(iClus);
1663 
1664  // Get the particle type of the cluster
1665  ContributorType contributorType = kUndefined;
1666  Int_t label = TMath::Abs(clus->GetLabel());
1667  if (label > 0) {
1668  contributorType = GetContributorType(clus, mcevent, label);
1669  }
1670 
1671  // Plot M02 for each particle type
1672  histname = "JetPerformanceMC/hM02VsContributorTypeJets";
1673  Double_t x[4] = {clus->GetM02(), clus->GetNonLinCorrEnergy(), contributorType, jetPt};
1674  fHistManager.FillTHnSparse(histname, x);
1675 
1676  // If the cluster is a hadron, sum its energy to compute the jet's hadronic calo energy
1677  Bool_t isHadron = IsHadron(contributorType);
1678  if (isHadron) {
1679  Bool_t hasMatchedTrack = (clus->GetNTracksMatched() > 0);
1680  //Bool_t hasMatchedTrack = ((clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy()) > 1e-3);
1681  if (hasMatchedTrack) {
1682  hadCaloEnergyMatchedNonlincorr += clus->GetNonLinCorrEnergy();
1683  hadCaloEnergyMatchedHadCorr += clus->GetHadCorrEnergy();
1684  }
1685  else {
1686  hadCaloEnergyUnmatched += clus->GetNonLinCorrEnergy();
1687  }
1688  }
1689 
1690  }
1691 
1692  // Fill hadronic calo energy in each jet
1693 
1694  // (Jet pT, Summed energy of hadronic clusters without a matched track)
1695  histname = "JetPerformance/hHadCaloEnergyUnmatched";
1696  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyUnmatched);
1697 
1698  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction))
1699  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
1700  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedNonlincorr);
1701 
1702  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction))
1703  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
1704  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedHadCorr);
1705 
1706  // Loop through particle types, and plot jet composition for each particle type
1707  histname = "JetPerformanceMC/hJetComposition";
1708  for (Int_t type = 0; type < 8; type++) {
1709 
1710  ContributorType contributorType = kUndefined;
1711  Double_t nSum = 0;
1712  Double_t pTsum = 0;
1713 
1714  // Loop through clusters in jet, and add to sum if cluster matches particle type
1715  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1716 
1717  clus = jet->Cluster(iClus);
1718 
1719  Int_t label = TMath::Abs(clus->GetLabel());
1720  if (label > 0) {
1721  contributorType = GetContributorType(clus, mcevent, label);
1722  }
1723 
1724  if (type == contributorType) {
1725  nSum++;
1726  pTsum += clus->GetNonLinCorrEnergy();
1727  }
1728  }
1729 
1730  // Fill jet composition histogram
1731  Double_t x[4] = {jetPt, 1.*type, nSum, pTsum};
1732  fHistManager.FillTHnSparse(histname, x);
1733 
1734  }
1735  }
1736 }
1737 
1742 {
1743  // Loop over tracks and clusters in order to:
1744  // (1) Compute scale factor for full jets
1745  // (2) Compute delta-pT for full jets, with the random cone method
1746 
1747  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
1748  Double_t etaTPC = 0.9;
1749  Double_t etaEMCal = 0.7;
1750  //Double_t etaMinDCal = 0.22;
1751  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
1752  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
1753  //Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
1754  //Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
1755 
1756  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
1757  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
1758  //Double_t accDCalRegion = 2 * etaEMCal * (phiMaxDCal - phiMinDCal);
1759 
1760  // Loop over jet containers
1761  AliJetContainer* jetCont = 0;
1762  TIter nextJetColl(&fJetCollArray);
1763  while ((jetCont = static_cast<AliJetContainer*>(nextJetColl()))) {
1764 
1765  // Define fiducial acceptances, to be used to generate random cones
1766  TRandom3* r = new TRandom3(0);
1767  Double_t jetR = jetCont->GetJetRadius();
1768  Double_t etaEMCalfid = etaEMCal - jetR;
1769  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
1770  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
1771 
1772  // Generate EMCal random cone eta-phi
1773  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
1774  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
1775 
1776  // Initialize the various sums to 0
1777  Double_t trackPtSumTPC = 0;
1778  Double_t trackPtSumEMCal = 0;
1779  Double_t trackPtSumEMCalRC = 0;
1780  Double_t clusESumEMCal = 0;
1781  Double_t clusESumEMCalRC = 0;
1782 
1783  // Define a 2D vector (initialized to 0) to store the sum of track pT, and another for cluster ET
1784  std::vector<std::vector<Double_t>> trackPtSumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1785  std::vector<std::vector<Double_t>> clusESumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1786 
1787  // Loop over tracks. Sum the track pT:
1788  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal random cone,
1789  // Note: Loops over all det-level track containers. For data there should be only one. For embedding, there should be signal+background tracks.
1790  AliParticleContainer * partCont = 0;
1791  AliTLorentzVector track;
1792  Double_t trackEta;
1793  Double_t trackPhi;
1794  Double_t trackPt;
1795  Double_t deltaR;
1796  TIter nextPartCont(&fParticleCollArray);
1797  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
1798 
1799  TString partContName = partCont->GetName();
1800  if (!partContName.CompareTo("tracks")) {
1801 
1802  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(partCont);
1803  for (auto trackIterator : trackCont->accepted_momentum() ) {
1804 
1805  track.Clear();
1806  track = trackIterator.first;
1807  trackEta = track.Eta();
1808  trackPhi = track.Phi_0_2pi();
1809  trackPt = track.Pt();
1810 
1811  // (1)
1812  if (TMath::Abs(trackEta) < etaTPC) {
1813  trackPtSumTPC += trackPt;
1814  }
1815 
1816  // (2)
1817  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
1818  trackPtSumEMCal += trackPt;
1819  }
1820 
1821  // (3)
1822  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
1823  if (deltaR < jetR) {
1824  trackPtSumEMCalRC += trackPt;
1825  }
1826  }
1827  }
1828  }
1829 
1830  // Loop over clusters. Sum the cluster ET:
1831  // (1) in the EMCal, (2) in the EMCal random cone
1833  AliTLorentzVector clus;
1834  Double_t clusEta;
1835  Double_t clusPhi;
1836  Double_t clusE;
1837  for (auto clusIterator : clusCont->accepted_momentum() ) {
1838 
1839  clus.Clear();
1840  clus = clusIterator.first;
1841  clusEta = clus.Eta();
1842  clusPhi = clus.Phi_0_2pi();
1843  clusE = clus.E();
1844 
1845  // (1)
1846  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
1847  clusESumEMCal += clusE;
1848  }
1849 
1850  // (2)
1851  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
1852  if (deltaR < jetR) {
1853  clusESumEMCalRC += clusE;
1854  }
1855 
1856  }
1857 
1858  // Compute the scale factor for EMCal, as a function of centrality
1859  Double_t numerator = (trackPtSumEMCal + clusESumEMCal) / accEMCal;
1860  Double_t denominator = trackPtSumTPC / accTPC;
1861  Double_t scaleFactor = numerator / denominator;
1862  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
1863  fHistManager.FillTH2(histname, fCent, scaleFactor);
1864 
1865  // Compute delta pT for EMCal, as a function of centrality
1866  Double_t rho = jetCont->GetRhoVal();
1867  Double_t deltaPt = trackPtSumEMCalRC + clusESumEMCalRC - rho * TMath::Pi() * jetR * jetR;
1868  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
1869  fHistManager.FillTH2(histname, fCent, deltaPt);
1870 
1871  delete r;
1872 
1873  }
1874 
1875 }
1876 
1882 {
1883  TString histname;
1884  AliJetContainer* jets = 0;
1885  TIter nextJetColl(&fJetCollArray);
1886  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1887  TString jetContName = jets->GetName();
1888 
1889  Double_t rhoVal = 0;
1890  if (jets->GetRhoParameter()) {
1891  rhoVal = jets->GetRhoVal();
1892  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
1893  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1894  }
1895 
1896  for (auto jet : jets->all()) {
1897 
1898  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1899  Float_t corrPt = GetJetPt(jet, rhoVal);
1900 
1901  // A vs. pT (fill before area cut)
1902  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
1903  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
1904 
1905  // Rejection reason
1906  UInt_t rejectionReason = 0;
1907  if (!jets->AcceptJet(jet, rejectionReason)) {
1908  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1909  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1910  continue;
1911  }
1912 
1913  // compute jet acceptance type
1914  Double_t type = GetJetType(jet);
1915  if ( (type != kEMCal) && (type != kDCal) ) {
1916  continue;
1917  }
1918 
1919  // (Centrality, pT, NEF, calo type)
1920  if (type == kEMCal) {
1921  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtEMCal", jets->GetArrayName().Data());
1922  }
1923  else if (type == kDCal) {
1924  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtDCal", jets->GetArrayName().Data());
1925  }
1926  fHistManager.FillTH3(histname, fCent, corrPt, jet->NEF());
1927 
1928  // pT-leading vs. pT
1929  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1930  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1931 
1932  // (Centrality, pT, z-leading (charged), calo type)
1933  if (type == kEMCal) {
1934  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtEMCal", jets->GetArrayName().Data());
1935  }
1936  else if (type == kDCal) {
1937  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPtDCal", jets->GetArrayName().Data());
1938  }
1939  TLorentzVector leadPart;
1940  jets->GetLeadingHadronMomentum(leadPart, jet);
1941  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1942  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
1943  fHistManager.FillTH3(histname, fCent, corrPt, z);
1944 
1945  // (Centrality, pT, z (charged), calo type)
1946  if (type == kEMCal) {
1947  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtEMCal", jets->GetArrayName().Data());
1948  }
1949  else if (type == kDCal) {
1950  histname = TString::Format("%s/TriggerSimHistograms/hZVsPtDCal", jets->GetArrayName().Data());
1951  }
1952  const AliVTrack* track;
1953  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1954  track = static_cast<AliVTrack*>(jet->Track(i));
1955  z = track->Pt() / TMath::Abs(corrPt);
1956  fHistManager.FillTH3(histname, fCent, corrPt, z);
1957  }
1958 
1959  // (Centrality, pT, Nconst)
1960  if (type == kEMCal) {
1961  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtEMCal", jets->GetArrayName().Data());
1962  }
1963  else if (type == kDCal) {
1964  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPtDCal", jets->GetArrayName().Data());
1965  }
1966  fHistManager.FillTH3(histname, fCent, corrPt, 1.*jet->GetNumberOfConstituents());
1967 
1968  } //jet loop
1969 
1970  }
1971 }
1972 
1978 {
1979  TString histname;
1980  AliJetContainer* jets = 0;
1981  const AliEmcalJet* matchedJet = nullptr;
1982  TIter nextJetColl(&fJetCollArray);
1983  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1984  TString jetContName = jets->GetName();
1985 
1986  // Only loop over jets in the detector-level jet container
1987  if (jetContName.Contains("mcparticles")) {
1988  continue;
1989  }
1990 
1991  Double_t rhoVal = 0;
1992  if (jets->GetRhoParameter()) {
1993  rhoVal = jets->GetRhoVal();
1994  }
1995 
1996  for (auto jet : jets->accepted()) {
1997 
1998  // Get the matched jet, if it exists
1999  matchedJet = jet->MatchedJet();
2000  if (!matchedJet) {
2001  continue;
2002  }
2003 
2004  // compute jet acceptance type
2005  Double_t type = GetJetType(jet);
2006  if ( (type != kEMCal) && (type != kDCal) ) {
2007  continue;
2008  }
2009 
2010  Float_t detPt = GetJetPt(jet, rhoVal);
2011  Float_t truthPt = matchedJet->Pt();
2012 
2013  // Fill response matrix (centrality, pT-truth, pT-det)
2014  if (type == kEMCal) {
2015  histname = "MatchedJetHistograms/hResponseMatrixEMCal";
2016  }
2017  else if (type == kDCal) {
2018  histname = "MatchedJetHistograms/hResponseMatrixDCal";
2019  }
2020  fHistManager.FillTH3(histname, fCent, truthPt, detPt);
2021 
2022  // Fill JES shift (centrality, pT-truth, (pT-det - pT-truth) / pT-truth)
2023  if (type == kEMCal) {
2024  histname = "MatchedJetHistograms/hJESshiftEMCal";
2025  }
2026  else if (type == kDCal) {
2027  histname = "MatchedJetHistograms/hJESshiftDCal";
2028  }
2029  fHistManager.FillTH3(histname, fCent, truthPt, (detPt-truthPt)/truthPt );
2030 
2031  // Fill NEF of det-level matched jets (centrality, pT-truth, NEF)
2032  histname = "MatchedJetHistograms/hNEFVsPt";
2033  fHistManager.FillTH3(histname, fCent, truthPt, jet->NEF());
2034 
2035  // Fill z-leading (charged) of det-level matched jets (centrality, pT-truth, z-leading)
2036  histname = "MatchedJetHistograms/hZLeadingVsPt";
2037  TLorentzVector leadPart;
2038  jets->GetLeadingHadronMomentum(leadPart, jet);
2039  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
2040  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin <1
2041  fHistManager.FillTH3(histname, fCent, truthPt, z);
2042 
2043  // Fill matching distance (centrality, pT-truth, R)
2044  histname = "MatchedJetHistograms/hMatchingDistance";
2045  fHistManager.FillTH3(histname, fCent, truthPt, jet->ClosestJetDistance());
2046 
2047  } //jet loop
2048  }
2049 }
2050 
2051 /*
2052  * Compute the MC particle type of a given cluster contributor, using the MC particle container
2053  */
2055 {
2056  ContributorType contributorType = kUndefined;
2057 
2058  AliAODMCParticle *part = fGeneratorLevel->GetMCParticleWithLabel(label);
2059  if (part) {
2060 
2061  TString histname = "ClusterHistogramsMC/hClusterRejectionReasonMC";
2062  UInt_t rejectionReason = 0;
2063  if (!fGeneratorLevel->AcceptMCParticle(part, rejectionReason)) {
2064  fHistManager.FillTH2(histname, fGeneratorLevel->GetRejectionReasonBitPosition(rejectionReason), clus->GetNonLinCorrEnergy());
2065  return contributorType;
2066  }
2067 
2068  if (part->GetGeneratorIndex() == 0) { // generator index in cocktail
2069 
2070  // select charged pions, protons, kaons, electrons, muons
2071  Int_t pdg = part->PdgCode();
2072 
2073  if (pdg == 22) { // gamma 22
2074  contributorType = kPhoton;
2075  }
2076  else if (TMath::Abs(pdg) == 211) { // pi+ 211 (abs value ensures both particles and antiparticles are included)
2077  contributorType = kChargedPion;
2078  }
2079  else if (pdg == 2212) { // proton 2212
2080  contributorType = kProton;
2081  }
2082  else if (pdg == -2212) {
2083  contributorType = kAntiProton;
2084  }
2085  else if (TMath::Abs(pdg) == 321) { // K+ 321
2086  contributorType = kChargedKaon;
2087  }
2088  else if (pdg == 130) { // K0L 130
2089  contributorType = kK0L;
2090  }
2091  else if (pdg == 2112) { // neutron 2112
2092  contributorType = kNeutron;
2093  }
2094  else if (pdg == -2112) {
2095  contributorType = kAntiNeutron;
2096  }
2097  else if (TMath::Abs(pdg) == 11) { // e- 11
2098  contributorType = kElectron;
2099  }
2100  else if (TMath::Abs(pdg) == 13) { // mu- 13
2101  contributorType = kMuon;
2102  }
2103  else {
2104  contributorType = kOther;
2105  }
2106  }
2107  }
2108  return contributorType;
2109 }
2110 
2115 {
2116  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
2117  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
2118  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
2119  return deltaR;
2120 }
2121 
2126 {
2127  UInt_t jetType = jet->GetJetAcceptanceType();
2128  Double_t type = -1;
2129  if (jetType & AliEmcalJet::kEMCAL) {
2130  type = kEMCal;
2131  }
2132  else if (jetType & AliEmcalJet::kDCALonly) {
2133  type = kDCal;
2134  }
2135 
2136  return type;
2137 }
2138 
2143 {
2144  Double_t pT = jet->Pt() - rho * jet->Area();
2145  return pT;
2146 }
2147 
2152 {
2153  return (contributor == kChargedPion) || (contributor == kProton) || (contributor == kAntiProton) || (contributor == kChargedKaon) || (contributor == kK0L) || (contributor == kNeutron) || (contributor == kAntiNeutron);
2154 }
2155 
2160  const char *ntracks,
2161  const char *nclusters,
2162  const char *nGenLev,
2163  const Double_t minTrPt,
2164  const Double_t minClPt,
2165  const char *suffix)
2166 {
2167  // Get the pointer to the existing analysis manager via the static access method.
2168  //==============================================================================
2169  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2170  if (!mgr)
2171  {
2172  ::Error("AddTaskEmcalJetPerformance", "No analysis manager to connect to.");
2173  return 0;
2174  }
2175 
2176  // Check the analysis type using the event handlers connected to the analysis manager.
2177  //==============================================================================
2178  AliVEventHandler* handler = mgr->GetInputEventHandler();
2179  if (!handler)
2180  {
2181  ::Error("AddTaskEmcalJetPerformance", "This task requires an input event handler");
2182  return 0;
2183  }
2184 
2185  enum EDataType_t {
2186  kUnknown,
2187  kESD,
2188  kAOD
2189  };
2190 
2191  EDataType_t dataType = kUnknown;
2192 
2193  if (handler->InheritsFrom("AliESDInputHandler")) {
2194  dataType = kESD;
2195  }
2196  else if (handler->InheritsFrom("AliAODInputHandler")) {
2197  dataType = kAOD;
2198  }
2199 
2200  //-------------------------------------------------------
2201  // Init the task and do settings
2202  //-------------------------------------------------------
2203 
2204  TString trackName(ntracks);
2205  TString clusName(nclusters);
2206 
2207  if (trackName == "usedefault") {
2208  if (dataType == kESD) {
2209  trackName = "Tracks";
2210  }
2211  else if (dataType == kAOD) {
2212  trackName = "tracks";
2213  }
2214  else {
2215  trackName = "";
2216  }
2217  }
2218 
2219  if (clusName == "usedefault") {
2220  if (dataType == kESD) {
2221  clusName = "CaloClusters";
2222  }
2223  else if (dataType == kAOD) {
2224  clusName = "caloClusters";
2225  }
2226  else {
2227  clusName = "";
2228  }
2229  }
2230 
2231  TString name("AliAnalysisTaskEmcalJetPerformance");
2232  if (!trackName.IsNull()) {
2233  name += "_";
2234  name += trackName;
2235  }
2236  if (!clusName.IsNull()) {
2237  name += "_";
2238  name += clusName;
2239  }
2240  if (strcmp(suffix,"") != 0) {
2241  name += "_";
2242  name += suffix;
2243  }
2244 
2246  // Configure jet performance task
2248 
2250  // Create track and cluster containers with the standard cuts
2251 
2252  AliParticleContainer* partCont = 0;
2253  if (trackName == "mcparticles") {
2254  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
2255  partCont = mcpartCont;
2256  }
2257  else if (trackName == "tracks" || trackName == "Tracks") {
2258  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
2259  partCont = trackCont;
2260  }
2261  if (partCont) partCont->SetParticlePtCut(minTrPt);
2262  if (partCont) task->AdoptParticleContainer(partCont);
2263 
2264  // Add the generator-level container, if specified
2265  if (nGenLev && strcmp(nGenLev,"")!=0) {
2266  AliMCParticleContainer* mcpartCont = task->AddMCParticleContainer(nGenLev);
2267  mcpartCont->SelectPhysicalPrimaries(kTRUE);
2268  mcpartCont->SetParticlePtCut(0);
2269  }
2270 
2271  // Add the cluster container
2272  AliClusterContainer* clusCont = 0;
2273  if (!clusName.IsNull()) {
2274  clusCont = new AliClusterContainer(clusName);
2275  clusCont->SetClusECut(0.);
2276  clusCont->SetClusPtCut(0.);
2277  }
2278  if (clusCont) task->AdoptClusterContainer(clusCont);
2279 
2280  //-------------------------------------------------------
2281  // Final settings, pass to manager and set the containers
2282  //-------------------------------------------------------
2283 
2284  mgr->AddTask(task);
2285 
2286  // Create containers for input/output
2287  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
2288  TString contname(name);
2289  contname += "_histos";
2290  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
2291  TList::Class(),AliAnalysisManager::kOutputContainer,
2292  Form("%s", AliAnalysisManager::GetCommonFileName()));
2293  mgr->ConnectInput (task, 0, cinput1 );
2294  mgr->ConnectOutput (task, 1, coutput1 );
2295 
2296  return task;
2297 }
2298 
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)
static AliEmcalDownscaleFactorsOCDB * Instance()
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.
void FillTH3(const char *hname, double x, double y, double z, double weight=1., Option_t *opt="")
Fill a 3D histogram within the container.
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
AliMCParticleContainer * fGeneratorLevel
! generator level container
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
Container for particles within the EMCAL framework.
void SetCaloTriggerPatchInfoName(const char *n)
Bool_t fPlotClusterHistograms
Set whether to plot cluster histograms.
Double_t GetDeltaR(const AliTLorentzVector *part, Double_t etaRef, Double_t phiRef)
void AdoptClusterContainer(AliClusterContainer *cont)
TObjArray fParticleCollArray
particle/track collection array
virtual Bool_t AcceptMCParticle(const AliAODMCParticle *vp, UInt_t &rejectionReason) const
bool AddQAPlotsToList(TList *list)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const Int_t nPtBins
Bool_t fComputeMBDownscaling
Set whether to compute and plot MB downscaling factors.
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.
std::vector< TString > GetTriggerClasses() const
Double_t fCent
!event centrality
Photon (direct or decay) is the only contributor.
Two particles from merged pi0 are the only contributors.
Bool_t 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.
AliEmcalList * fOutput
!output list
ContributorType GetContributorType(const AliVCluster *clus, const AliMCEvent *mcevent, Int_t label)
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.
Handler for downscale factors for various triggers obtained from the OCDB.
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()