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