AliPhysics  a56b849 (a56b849)
 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 "AliTLorentzVector.h"
45 #include "AliEmcalJet.h"
46 #include "AliRhoParameter.h"
47 #include "AliJetContainer.h"
48 #include "AliParticleContainer.h"
49 #include "AliClusterContainer.h"
50 #include "AliEMCALGeometry.h"
51 #include "AliOADBContainer.h"
53 #include "AliEMCALTriggerPatchInfo.h"
55 #include "AliMCAnalysisUtils.h"
56 
58 
62 
68  fPlotJetHistograms(kFALSE),
69  fPlotClusterHistograms(kFALSE),
70  fComputeBackground(kFALSE),
71  fDoTriggerSimulation(kFALSE),
72  fComputeMBDownscaling(kFALSE),
73  fMaxPt(200),
74  fNEtaBins(40),
75  fNPhiBins(200),
76  fNCentHistBins(0),
77  fCentHistBins(0),
78  fNPtHistBins(0),
79  fPtHistBins(0),
80  fNM02HistBins(0),
81  fM02HistBins(0),
82  fMBUpscaleFactor(1.),
83  fMedianEMCal(0.),
84  fMedianDCal(0.),
85  fkEMCEJE(kFALSE),
86  fEmbeddingQA(),
87  fUseAliEventCuts(kTRUE),
88  fEventCuts(0),
89  fEventCutList(0),
90  fUseManualEventCuts(kFALSE),
91  fGeneratorLevel(0),
92  fHistManager()
93 {
95 }
96 
103  AliAnalysisTaskEmcalJet(name, kTRUE),
104  fPlotJetHistograms(kFALSE),
105  fPlotClusterHistograms(kFALSE),
106  fComputeBackground(kFALSE),
107  fDoTriggerSimulation(kFALSE),
108  fComputeMBDownscaling(kFALSE),
109  fMaxPt(200),
110  fNEtaBins(40),
111  fNPhiBins(200),
112  fNCentHistBins(0),
113  fCentHistBins(0),
114  fNPtHistBins(0),
115  fPtHistBins(0),
116  fNM02HistBins(0),
117  fM02HistBins(0),
118  fMBUpscaleFactor(1.),
119  fMedianEMCal(0.),
120  fMedianDCal(0.),
121  fkEMCEJE(kFALSE),
122  fEmbeddingQA(),
123  fUseAliEventCuts(kTRUE),
124  fEventCuts(0),
125  fEventCutList(0),
126  fUseManualEventCuts(kFALSE),
127  fGeneratorLevel(0),
128  fHistManager(name)
129 {
131 }
132 
137 {
138 }
139 
144 {
145  fNCentHistBins = 4;
147  fCentHistBins[0] = 0;
148  fCentHistBins[1] = 10;
149  fCentHistBins[2] = 30;
150  fCentHistBins[3] = 50;
151  fCentHistBins[4] = 90;
152 
153  fNPtHistBins = 82;
156  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
157  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
158  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
159  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
160  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
161  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
162 
163  fNM02HistBins = 81;
165  GenerateFixedBinArray(35, 0, 0.7, fM02HistBins);
166  GenerateFixedBinArray(6, 0.7, 1., fM02HistBins+35);
167  GenerateFixedBinArray(20, 1., 3., fM02HistBins+41);
168  GenerateFixedBinArray(10, 3., 5., fM02HistBins+61);
169  GenerateFixedBinArray(10, 5., 10., fM02HistBins+71);
170 }
171 
177 {
179 
180  // Intialize AliEventCuts
181  if (fUseAliEventCuts) {
182  fEventCutList = new TList();
183  fEventCutList ->SetOwner();
184  fEventCutList ->SetName("EventCutOutput");
185 
186  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
187  if(fUseManualEventCuts==1)
188  {
189  fEventCuts.SetManualMode();
190  // Configure manual settings here
191  // ...
192  }
193  fEventCuts.AddQAplotsToList(fEventCutList);
194  fOutput->Add(fEventCutList);
195  }
196 
197  // Get the MC particle branch, in case it exists
198  fGeneratorLevel = GetMCParticleContainer("mcparticles");
199 
200  if (fPlotJetHistograms) {
202  }
205  }
206  if (fComputeBackground) {
208  }
209  if (fDoTriggerSimulation) {
211  }
212 
213  // Initialize embedding QA
215  if (embeddingHelper) {
216  bool res = fEmbeddingQA.Initialize();
217  if (res) {
219  }
220  }
221 
222  TIter next(fHistManager.GetListOfHistograms());
223  TObject* obj = 0;
224  while ((obj = next())) {
225  fOutput->Add(obj);
226  }
227 
228  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
229 }
230 
231 /*
232  * This function allocates the histograms for single jets.
233  * A set of histograms is allocated per each jet container.
234  */
236 {
237  TString histname;
238  TString title;
239 
240  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
241 
242  AliJetContainer* jets = 0;
243  TIter nextJetColl(&fJetCollArray);
244  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
245 
246  // Jet rejection reason
247  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
248  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
249  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
250  SetRejectionReasonLabels(hist->GetXaxis());
251 
252  // Rho vs. Centrality
253  if (!jets->GetRhoName().IsNull()) {
254  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
255  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
256  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
257  }
258 
259  // (Centrality, pT, calo type)
260  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
261  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
262  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, 3, -0.5, 2.5);
263 
264  // (Centrality, eta, phi, pT, NEF, calo type)
265  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
266  title = histname + ";Centrality (%);#eta_{jet} (rad);#phi_{jet} (rad);#it{p}_{T}^{corr} (GeV/#it{c});NEF;type";
267  Int_t nbins1[6] = {20, fNEtaBins, fNPhiBins, nPtBins, 50, 3};
268  Double_t min1[6] = {0, -0.5,1., 0, 0, -0.5};
269  Double_t max1[6] = {100, 0.5,6., fMaxPt, 1., 2.5};
270  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 6, nbins1, min1, max1);
271 
272  // (Centrality, pT upscaled, calo type)
273  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
274  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
275  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, 3, -0.5, 2.5, "s");
276 
277  // pT-leading vs. pT
278  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
279  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
280  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
281 
282  // A vs. pT
283  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
284  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
285  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/3, 0, 0.5);
286 
287  // (Centrality, pT, z-leading (charged), calo type)
288  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
289  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading};type";
290  Int_t nbins2[4] = {20, nPtBins, 50, 3};
291  Double_t min2[4] = {0, 0, 0, -0.5};
292  Double_t max2[4] = {100, fMaxPt, 1., 2.5};
293  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins2, min2, max2);
294 
295  // (Centrality, pT, z (charged), calo type)
296  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
297  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z};type";
298  Int_t nbins3[4] = {20, nPtBins, 50, 3};
299  Double_t min3[4] = {0, 0, 0, -0.5};
300  Double_t max3[4] = {100, fMaxPt, 1., 2.5};
301  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins3, min3, max3);
302 
303  // (Centrality, pT, Nconst, calo type)
304  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
305  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents;type";
306  Int_t nbins4[4] = {20, nPtBins, 50, 3};
307  Double_t min4[4] = {0, 0, 0, -0.5};
308  Double_t max4[4] = {100, fMaxPt, fMaxPt, 2.5};
309  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins4, min4, max4);
310 
311  // (Median patch energy, calo type, jet pT, centrality)
312  if (fDoTriggerSimulation) {
313  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
314  title = histname + ";#it{E}_{patch,med};type;#it{p}_{T}^{corr} (GeV/#it{c});Centrality (%)";
315  Int_t nbins5[4] = {100, 2, nPtBins, 50};
316  Double_t min5[4] = {0,-0.5, 0, 0};
317  Double_t max5[4] = {50,1.5, fMaxPt, 100};
318  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins5, min5, max5);
319  }
320 
321  }
322 
323  // MB downscale factor histogram
324  if (fComputeMBDownscaling) {
325  histname = "Trigger/hMBDownscaleFactor";
326  title = histname + ";Downscale factor;counts";
327  fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 200);
328  }
329 
330 }
331 
332 /*
333  * This function allocates the histograms for the calorimeter performance study.
334  */
336 {
337  TString histname;
338  TString htitle;
339  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
340 
341  const Int_t nRejBins = 32;
342  Double_t* rejReasonBins = new Double_t[nRejBins+1];
343  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
344 
345  AliEmcalContainer* cont = 0;
346  TIter nextClusColl(&fClusterCollArray);
347  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
348 
349  // rejection reason plot, to make efficiency correction
350  histname = TString::Format("%s/hClusterRejectionReasonEMCal", cont->GetArrayName().Data());
351  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
352  TH2* hist = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
353  SetRejectionReasonLabels(hist->GetXaxis());
354 
355  histname = TString::Format("%s/hClusterRejectionReasonMC", cont->GetArrayName().Data());
356  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
357  TH2* histMC = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
358  SetRejectionReasonLabels(histMC->GetXaxis());
359  }
360 
361  const Int_t nParticleTypes = 8;
362  Double_t *particleTypeBins = GenerateFixedBinArray(nParticleTypes, -0.5, 7.5);
363 
364  // M02 vs. Energy vs. Particle type
365  histname = "JetPerformance/hM02VsParticleType";
366  htitle = histname + ";M02;#it{E}_{clus} (GeV); Particle type";
367  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNM02HistBins, fM02HistBins, fNPtHistBins, fPtHistBins, nParticleTypes, particleTypeBins);
368 
369  // M02 vs. Energy vs. Particle type vs. Jet pT, for particles inside jets
370  Int_t dim = 0;
371  TString title[20];
372  Int_t nbins[20] = {0};
373  Double_t min[30] = {0.};
374  Double_t max[30] = {0.};
375  Double_t *binEdges[20] = {0};
376 
377  title[dim] = "M02";
378  nbins[dim] = fNM02HistBins;
379  binEdges[dim] = fM02HistBins;
380  min[dim] = fM02HistBins[0];
381  max[dim] = fM02HistBins[fNM02HistBins];
382  dim++;
383 
384  title[dim] = "#it{E}_{clus} (GeV)";
385  nbins[dim] = fNPtHistBins;
386  binEdges[dim] = fPtHistBins;
387  min[dim] = fPtHistBins[0];
388  max[dim] = fPtHistBins[fNPtHistBins];
389  dim++;
390 
391  title[dim] = "Particle type";
392  nbins[dim] = nParticleTypes;
393  min[dim] = -0.5;
394  max[dim] = 7.5;
395  binEdges[dim] = particleTypeBins;
396  dim++;
397 
398  title[dim] = "#it{p}_{T,jet}^{corr}";
399  nbins[dim] = nPtBins;
400  min[dim] = 0;
401  max[dim] = fMaxPt;
402  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
403  dim++;
404 
405  TString thnname = "JetPerformance/hM02VsParticleTypeJets";
406  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
407  for (Int_t i = 0; i < dim; i++) {
408  hn->GetAxis(i)->SetTitle(title[i]);
409  hn->SetBinEdges(i, binEdges[i]);
410  }
411 
412  // Particle composition inside each jet -- jet pT vs. particle type vs. particle number vs. particle pT sum
413  // (One entry per jet for each particle type)
414  dim = 0;
415 
416  title[dim] = "#it{p}_{T,jet}^{corr}";
417  nbins[dim] = nPtBins;
418  min[dim] = 0;
419  max[dim] = fMaxPt;
420  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
421  dim++;
422 
423  title[dim] = "Particle type";
424  nbins[dim] = nParticleTypes;
425  min[dim] = -0.5;
426  max[dim] = 7.5;
427  binEdges[dim] = particleTypeBins;
428  dim++;
429 
430  title[dim] = "N";
431  nbins[dim] = 30;
432  min[dim] = -0.5;
433  max[dim] = 29.5;
434  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
435  dim++;
436 
437  title[dim] = "#it{p}_{T,sum} (GeV)";
438  nbins[dim] = fNPtHistBins;
439  binEdges[dim] = fPtHistBins;
440  min[dim] = fPtHistBins[0];
441  max[dim] = fPtHistBins[fNPtHistBins];
442  dim++;
443 
444  thnname = "JetPerformance/hJetComposition";
445  THnSparse* thn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
446  for (Int_t i = 0; i < dim; i++) {
447  thn->GetAxis(i)->SetTitle(title[i]);
448  thn->SetBinEdges(i, binEdges[i]);
449  }
450 
451  // Hadronic calo energy in each jet
452 
453  // Jet pT vs. Summed energy of hadronic clusters without a matched track
454  histname = "JetPerformance/hHadCaloEnergyUnmatched";
455  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
456  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
457 
458  // Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction)
459  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
460  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
461  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
462 
463  // Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction)
464  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
465  htitle = histname + ";#it{p}_{T,jet} (GeV);#it{p}_{T,had} (GeV)";
466  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNPtHistBins, fPtHistBins);
467 
468 }
469 
470 /*
471  * This function allocates background subtraction histograms, if enabled.
472  * A set of histograms is allocated per each jet container.
473  */
475 {
476  TString histname;
477  TString title;
478 
479  AliJetContainer* jets = 0;
480  TIter nextJetColl(&fJetCollArray);
481  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
482 
483  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
484  title = histname + ";Centrality;Scale factor;counts";
485  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
486 
487  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
488  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
489  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
490 
491  }
492 }
493 
494 /*
495  * This function allocates the histograms for single jets, when the "simulated" trigger has been fired.
496  * A set of histograms is allocated per each jet container.
497  */
499 {
500  TString histname;
501  TString title;
502  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
503 
504  //----------------------------------------------
505  // Trigger patch histograms
506 
507  // patch eta vs. phi
508  histname = "TriggerSimHistograms/hEtaVsPhi";
509  title = histname + ";#eta_{patch} (rad);#phi_{patch} (rad)";
510  fHistManager.CreateTH2(histname.Data(), title.Data(), 140, -0.7, 0.7, 500, 1., 6.);
511 
512  // N patches
513  histname = "TriggerSimHistograms/hNPatches";
514  title = histname + ";#it{N}_{patches};type";
515  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, 0, 200, 2, -0.5, 1.5);
516 
517  // patch E vs. centrality
518  histname = "TriggerSimHistograms/hPatchE";
519  title = histname + ";Centrality (%);#it{E}_{patch} (GeV)";
520  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, nPtBins, 0, fMaxPt);
521 
522  // patch median vs. Centrality
523  histname = "TriggerSimHistograms/hPatchMedianE";
524  title = histname + ";Centrality (%);#it{E}_{patch,med} (GeV);type";
525  fHistManager.CreateTH3(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 50, 2, -0.5, 1.5);
526 
527  //----------------------------------------------
528  // Jet histograms for "triggered" events
529  AliJetContainer* jets = 0;
530  TIter nextJetColl(&fJetCollArray);
531  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
532 
533  // Jet rejection reason
534  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
535  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
536  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
537  SetRejectionReasonLabels(hist->GetXaxis());
538 
539  // Rho vs. Centrality
540  if (!jets->GetRhoName().IsNull()) {
541  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
542  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
543  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
544  }
545 
546  // (Centrality, pT, calo type)
547  histname = TString::Format("%s/TriggerSimHistograms/hCentVsPt", jets->GetArrayName().Data());
548  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
549  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, 3, -0.5, 2.5);
550 
551  // (Centrality, eta, phi, pT, NEF, calo type)
552  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
553  title = histname + ";Centrality (%);#eta_{jet} (rad);#phi_{jet} (rad);#it{p}_{T}^{corr} (GeV/#it{c});NEF;type";
554  Int_t nbins1[6] = {20, fNEtaBins, fNPhiBins, nPtBins, 50, 3};
555  Double_t min1[6] = {0, -0.5,1., 0, 0, -0.5};
556  Double_t max1[6] = {100, 0.5,6., fMaxPt, 1., 2.5};
557  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 6, nbins1, min1, max1);
558 
559  // pT-leading vs. pT
560  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
561  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
562  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
563 
564  // A vs. pT
565  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
566  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
567  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/3, 0, 0.5);
568 
569  // (Centrality, pT, z-leading (charged), calo type)
570  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
571  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading};type";
572  Int_t nbins2[4] = {20, nPtBins, 50, 3};
573  Double_t min2[4] = {0, 0, 0, -0.5};
574  Double_t max2[4] = {100, fMaxPt, 1., 2.5};
575  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins2, min2, max2);
576 
577  // z (charged) vs. pT
578  histname = TString::Format("%s/TriggerSimHistograms/hZVsPt", jets->GetArrayName().Data());
579  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z};type";
580  Int_t nbins3[4] = {20, nPtBins, 50, 3};
581  Double_t min3[4] = {0, 0, 0, -0.5};
582  Double_t max3[4] = {100, fMaxPt, 1., 2.5};
583  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins3, min3, max3);
584 
585  // (Centrality, pT, Nconst, calo type)
586  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPt", jets->GetArrayName().Data());
587  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents;type";
588  Int_t nbins4[4] = {20, nPtBins, 50, 3};
589  Double_t min4[4] = {0, 0, 0, -0.5};
590  Double_t max4[4] = {100, fMaxPt, fMaxPt, 2.5};
591  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins4, min4, max4);
592 
593  }
594 }
595 
601 {
602  // Configure base class to set fTriggerPatchInfo to array of trigger patches, each event
603  // (Need to call this before base class ExecOnce)
604  if (fDoTriggerSimulation) {
605  this->SetCaloTriggerPatchInfoName("EmcalTriggers");
606  }
607 
609 
610  fNeedEmcalGeom = kTRUE;
611 
612  // Check if trigger patches are loaded
613  if (fDoTriggerSimulation) {
614  if (fTriggerPatchInfo) {
615  TString objname(fTriggerPatchInfo->GetClass()->GetName());
616  TClass cls(objname);
617  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
618  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
619  GetName(), cls.GetName(), "EmcalTriggers"));
620  fTriggerPatchInfo = 0;
621  }
622  }
623  if (!fTriggerPatchInfo) {
624  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
625  return;
626  }
627  }
628 }
629 
634 
635  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
636 
638 
639  // Get instance of the downscale factor helper class
641  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
642 
643  // There are two possible min bias triggers for LHC15o
644  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
645  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
646  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
647 
648  // Get the downscale factor for whichever MB trigger exists in the given run
649  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
650  Double_t downscalefactor;
651  for (auto i : runtriggers) {
652  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
653  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
654  break;
655  }
656  }
657 
658  // Store the inverse of the downscale factor, used later to weight the pT spectrum
659  fMBUpscaleFactor = 1/downscalefactor;
660 
661  TString histname = "Trigger/hMBDownscaleFactor";
662  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
663 
664  }
665 
666 }
667 
672 {
673  if (fUseAliEventCuts) {
674  if (!fEventCuts.AcceptEvent(InputEvent()))
675  {
676  PostData(1, fOutput);
677  return kFALSE;
678  }
679  }
680  else {
682  }
683  return kTRUE;
684 }
685 
694 {
695  TString histname;
696  AliJetContainer* jetCont = 0;
697  TIter next(&fJetCollArray);
698  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
699  TString jetContName = jetCont->GetName();
700 
701  // Compute the full jet background scale factor and delta-pt
702  if (fComputeBackground) {
704  }
705 
706  // Do a simple trigger simulation (if requested)
707  if (fDoTriggerSimulation) {
709  }
710 
711  }
712 
713  // Only fill the embedding qa plots if:
714  // - We are using the embedding helper
715  // - The class has been initialized
716  // - Both jet collections are available
717  if (fEmbeddingQA.IsInitialized()) {
719  }
720 
721  return kTRUE;
722 }
723 
728 {
729  TString histname;
730 
731  // Check if trigger patches are loaded
732  if (fTriggerPatchInfo) {
733  TString objname(fTriggerPatchInfo->GetClass()->GetName());
734  TClass cls(objname);
735  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
736  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
737  GetName(), cls.GetName(), "EmcalTriggers"));
738  fTriggerPatchInfo = 0;
739  }
740  }
741  if (!fTriggerPatchInfo) {
742  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
743  return;
744  }
745 
746  // Compute patches in EMCal, DCal (I want offline simple trigger patch, i.e. patch calculated using FEE energy)
747  std::vector<Double_t> vecEMCal;
748  std::vector<Double_t> vecDCal;
749  for(auto p : *fTriggerPatchInfo){
750  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
751  if (recpatch) {
752 
753  if(!recpatch->IsJetHighSimple()) continue;
754 
755  histname = "TriggerSimHistograms/hEtaVsPhi";
756  fHistManager.FillTH2(histname.Data(), recpatch->GetEtaGeo(), recpatch->GetPhiGeo());
757 
758  histname = "TriggerSimHistograms/hPatchE";
759  fHistManager.FillTH2(histname.Data(), fCent, recpatch->GetPatchE());
760 
761  if (recpatch->IsEMCal()) {
762  vecEMCal.push_back(recpatch->GetPatchE());
763  } else {
764  vecDCal.push_back(recpatch->GetPatchE());
765  }
766 
767  }
768  }
769 
770  // Compute the median in each calorimeter
771  const Int_t nBkgPatchesEMCal = vecEMCal.size(); // 6*8;
772  const Int_t nBkgPatchesDCal = vecDCal.size(); // 4*5;
773  fMedianEMCal = TMath::Median(nBkgPatchesEMCal, &vecEMCal[0]); // point to array used internally by vector
774  fMedianDCal = TMath::Median(nBkgPatchesDCal, &vecDCal[0]);
775 
776  histname = "TriggerSimHistograms/hPatchMedianE";
777  fHistManager.FillTH3(histname.Data(), fCent, fMedianEMCal, kEMCal);
778  fHistManager.FillTH3(histname.Data(), fCent, fMedianDCal, kDCal);
779 
780  histname = "TriggerSimHistograms/hNPatches";
781  fHistManager.FillTH2(histname.Data(), nBkgPatchesEMCal, kEMCal);
782  fHistManager.FillTH2(histname.Data(), nBkgPatchesDCal, kDCal);
783 
784  // Then compute background subtracted patches, by subtracting from each patch the median patch E from the opposite hemisphere
785  // If a patch is above threshold, the event is "triggered"
786  Bool_t fkEMCEJE = kFALSE;
787  Double_t threshold = 20;
788  for(auto p : *fTriggerPatchInfo){
789  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
790  if (recpatch) {
791 
792  if(!recpatch->IsJetHighSimple()) continue;
793 
794  if (recpatch->IsEMCal()) {
795  if ((recpatch->GetPatchE() - fMedianDCal) > threshold) {
796  fkEMCEJE = kTRUE;
797  break;
798  }
799  } else {
800  if ((recpatch->GetPatchE() - fMedianEMCal) > threshold) {
801  fkEMCEJE = kTRUE;
802  break;
803  }
804  }
805  }
806  }
807 
808  if (fkEMCEJE) {
810  }
811 
812 }
813 
821 {
822 
823  if (fPlotJetHistograms) {
825  }
828  }
829 
830  return kTRUE;
831 }
832 
838 {
839  TString histname;
840  AliJetContainer* jets = 0;
841  TIter nextJetColl(&fJetCollArray);
842  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
843  TString jetContName = jets->GetName();
844 
845  Double_t rhoVal = 0;
846  if (jets->GetRhoParameter()) {
847  rhoVal = jets->GetRhoVal();
848  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
849  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
850  }
851 
852  for (auto jet : jets->all()) {
853 
854  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
855  Float_t corrPt = GetJetPt(jet, rhoVal);
856 
857  // A vs. pT (fill before area cut)
858  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
859  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
860 
861 
862  // Rejection reason
863  UInt_t rejectionReason = 0;
864  if (!jets->AcceptJet(jet, rejectionReason)) {
865  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
866  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
867  continue;
868  }
869 
870  // compute jet acceptance type
871  Double_t type = GetJetType(jet);
872 
873  // (Centrality, pT, calo type)
874  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
875  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type);
876 
877  // (Centrality, eta, phi, pT, NEF, calo type)
878  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
879  Double_t x[6] = {fCent, jet->Eta(), jet->Phi_0_2pi(), corrPt, jet->NEF(), type};
880  fHistManager.FillTHnSparse(histname, x);
881 
882  // (Centrality, pT upscaled, calo type)
883  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
884  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type, fMBUpscaleFactor);
885 
886  // pT-leading vs. pT
887  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
888  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
889 
890  // (Centrality, pT, z-leading (charged), calo type)
891  TLorentzVector leadPart;
892  jets->GetLeadingHadronMomentum(leadPart, jet);
893  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
894  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
895  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
896  Double_t y[4] = {fCent, corrPt, z, type};
897  fHistManager.FillTHnSparse(histname, y);
898 
899  // (Centrality, pT, z (charged), calo type)
900  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
901  AliVTrack* track;
902  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
903  track = static_cast<AliVTrack*>(jet->Track(i));
904  z = track->Pt() / TMath::Abs(corrPt);
905  Double_t y2[4] = {fCent, corrPt, z, type};
906  fHistManager.FillTHnSparse(histname, y2);
907  }
908 
909  // (Centrality, pT, Nconst, calo type)
910  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
911  Double_t a[4] = {fCent, corrPt, 1.*jet->GetNumberOfConstituents(), type};
912  fHistManager.FillTHnSparse(histname, a);
913 
914  // (Median patch energy, calo type, jet pT, centrality)
915  if (fDoTriggerSimulation) {
916  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
917  Double_t x[4] = {fMedianEMCal, kEMCal, corrPt, fCent};
918  fHistManager.FillTHnSparse(histname, x);
919  Double_t y[4] = {fMedianDCal, kDCal, corrPt, fCent};
920  fHistManager.FillTHnSparse(histname, y);
921  }
922 
923  } //jet loop
924 
925  }
926 }
927 
928 /*
929  * This function fills the histograms for the calorimeter performance study.
930  */
932 {
933  TString histname;
934 
935  // If MC, get the MC event
936  const AliMCEvent* mcevent = nullptr;
937  if (fGeneratorLevel) {
938  mcevent = MCEvent();
939  }
940  else {
941  return;
942  }
943 
944  // Loop through clusters, and plot M02 for each particle type
946  const AliVCluster* clus;
947  for (auto it : clusters->accepted_momentum()) {
948 
949  clus = it.second;
950 
951  // Include only EMCal clusters (reject DCal and PHOS clusters)
952  if (!clus->IsEMCAL()) {
953  continue;
954  }
955  if (it.first.Phi_0_2pi() > fgkEMCalDCalPhiDivide) {
956  continue;
957  }
958 
959  // If MC, determine the particle type
960  // Each detector-level cluster contains an array of labels of each truth-level particle contributing to the cluster
961  ParticleType particleType1 = kUndefined;
962 
963  Int_t label = TMath::Abs(clus->GetLabel()); // returns mc-label of particle that deposited the most energy in the cluster
964  if (label > 0) { // if the particle has a truth-level match, the label is > 0
965 
966  // Method 1: Use AliMCAnalysisUtils to identify all particles
967  particleType1 = GetParticleType1(clus, mcevent, clusters->GetArray());
968 
969  }
970 
971  // (M02, Eclus, part type)
972  histname = "JetPerformance/hM02VsParticleType";
973  fHistManager.FillTH3(histname, clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType1);
974 
975  }
976 
977  // Loop through jets, to fill various histograms
978  AliJetContainer* jets = GetJetContainer(0); // there is only a single, det-level jet finder here
979  for (const auto jet : jets->accepted()) {
980 
981  Double_t jetPt = GetJetPt(jet, 0);
982 
983  // Keep track of hadronic calo energy in each jet
984  Double_t hadCaloEnergyUnmatched = 0;
985  Double_t hadCaloEnergyMatchedNonlincorr = 0;
986  Double_t hadCaloEnergyMatchedHadCorr = 0;
987 
988  // Loop through clusters in each jet, to plot several histograms
989  Int_t nClusters = jet->GetNumberOfClusters();
990  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
991 
992  clus = jet->Cluster(iClus);
993 
994  // Get the particle type of the cluster
995  ParticleType particleType1 = kUndefined;
996  Int_t label = TMath::Abs(clus->GetLabel());
997  if (label > 0) {
998  particleType1 = GetParticleType1(clus, mcevent, clusters->GetArray());
999  }
1000 
1001  // Plot M02 for each particle type
1002  histname = "JetPerformance/hM02VsParticleTypeJets";
1003  Double_t x[4] = {clus->GetM02(), clus->GetNonLinCorrEnergy(), particleType1, jetPt};
1004  fHistManager.FillTHnSparse(histname, x);
1005 
1006  // If the cluster is a hadron, sum its energy to compute the jet's hadronic calo energy
1007  if (particleType1 == kHadron) {
1008  Bool_t hasMatchedTrack = (clus->GetNTracksMatched() > 0);
1009  //Bool_t hasMatchedTrack = ((clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy()) > 1e-3);
1010  if (hasMatchedTrack) {
1011  hadCaloEnergyMatchedNonlincorr += clus->GetNonLinCorrEnergy();
1012  hadCaloEnergyMatchedHadCorr += clus->GetHadCorrEnergy();
1013  }
1014  else {
1015  hadCaloEnergyUnmatched += clus->GetNonLinCorrEnergy();
1016  }
1017  }
1018 
1019  }
1020 
1021  // Fill hadronic calo energy in each jet
1022 
1023  // (Jet pT, Summed energy of hadronic clusters without a matched track)
1024  histname = "JetPerformance/hHadCaloEnergyUnmatched";
1025  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyUnmatched);
1026 
1027  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (before hadronic correction))
1028  histname = "JetPerformance/hHadCaloEnergyMatchedNonlincorr";
1029  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedNonlincorr);
1030 
1031  // (Jet pT vs. Summed energy of hadronic clusters with a matched track (after hadronic correction))
1032  histname = "JetPerformance/hHadCaloEnergyMatchedHadCorr";
1033  fHistManager.FillTH2(histname, jetPt, hadCaloEnergyMatchedHadCorr);
1034 
1035  // Loop through particle types, and plot jet composition for each particle type
1036  histname = "JetPerformance/hJetComposition";
1037  for (Int_t type = 0; type < 8; type++) {
1038 
1039  ParticleType particleType1 = kUndefined;
1040  Double_t nSum = 0;
1041  Double_t pTsum = 0;
1042 
1043  // Loop through clusters in jet, and add to sum if cluster matches particle type
1044  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1045 
1046  clus = jet->Cluster(iClus);
1047 
1048  Int_t label = TMath::Abs(clus->GetLabel());
1049  if (label > 0) {
1050  particleType1 = GetParticleType1(clus, mcevent, clusters->GetArray());
1051  }
1052 
1053  if (type == particleType1) {
1054  nSum++;
1055  pTsum += clus->GetNonLinCorrEnergy();
1056  }
1057  }
1058 
1059  // Fill jet composition histogram
1060  Double_t x[4] = {jetPt, 1.*type, nSum, pTsum};
1061  fHistManager.FillTHnSparse(histname, x);
1062 
1063  }
1064  }
1065 }
1066 
1071 {
1072  // Loop over tracks and clusters in order to:
1073  // (1) Compute scale factor for full jets
1074  // (2) Compute delta-pT for full jets, with the random cone method
1075 
1076  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
1077  Double_t etaTPC = 0.9;
1078  Double_t etaEMCal = 0.7;
1079  //Double_t etaMinDCal = 0.22;
1080  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
1081  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
1082  //Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
1083  //Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
1084 
1085  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
1086  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
1087  //Double_t accDCalRegion = 2 * etaEMCal * (phiMaxDCal - phiMinDCal);
1088 
1089  // Loop over jet containers
1090  AliJetContainer* jetCont = 0;
1091  TIter nextJetColl(&fJetCollArray);
1092  while ((jetCont = static_cast<AliJetContainer*>(nextJetColl()))) {
1093 
1094  // Define fiducial acceptances, to be used to generate random cones
1095  TRandom3* r = new TRandom3(0);
1096  Double_t jetR = jetCont->GetJetRadius();
1097  Double_t etaEMCalfid = etaEMCal - jetR;
1098  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
1099  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
1100 
1101  // Generate EMCal random cone eta-phi
1102  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
1103  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
1104 
1105  // Initialize the various sums to 0
1106  Double_t trackPtSumTPC = 0;
1107  Double_t trackPtSumEMCal = 0;
1108  Double_t trackPtSumEMCalRC = 0;
1109  Double_t clusESumEMCal = 0;
1110  Double_t clusESumEMCalRC = 0;
1111 
1112  // Define a 2D vector (initialized to 0) to store the sum of track pT, and another for cluster ET
1113  std::vector<std::vector<Double_t>> trackPtSumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1114  std::vector<std::vector<Double_t>> clusESumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1115 
1116  // Loop over tracks. Sum the track pT:
1117  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal random cone,
1118  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1119  AliTLorentzVector track;
1120  Double_t trackEta;
1121  Double_t trackPhi;
1122  Double_t trackPt;
1123  Double_t deltaR;
1124  for (auto trackIterator : trackCont->accepted_momentum() ) {
1125 
1126  track.Clear();
1127  track = trackIterator.first;
1128  trackEta = track.Eta();
1129  trackPhi = track.Phi_0_2pi();
1130  trackPt = track.Pt();
1131 
1132  // (1)
1133  if (TMath::Abs(trackEta) < etaTPC) {
1134  trackPtSumTPC += trackPt;
1135  }
1136 
1137  // (2)
1138  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
1139  trackPtSumEMCal += trackPt;
1140  }
1141 
1142  // (3)
1143  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
1144  if (deltaR < jetR) {
1145  trackPtSumEMCalRC += trackPt;
1146  }
1147 
1148  }
1149 
1150  // Loop over clusters. Sum the cluster ET:
1151  // (1) in the EMCal, (2) in the EMCal random cone
1153  AliTLorentzVector clus;
1154  Double_t clusEta;
1155  Double_t clusPhi;
1156  Double_t clusE;
1157  for (auto clusIterator : clusCont->accepted_momentum() ) {
1158 
1159  clus.Clear();
1160  clus = clusIterator.first;
1161  clusEta = clus.Eta();
1162  clusPhi = clus.Phi_0_2pi();
1163  clusE = clus.E();
1164 
1165  // (1)
1166  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
1167  clusESumEMCal += clusE;
1168  }
1169 
1170  // (2)
1171  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
1172  if (deltaR < jetR) {
1173  clusESumEMCalRC += clusE;
1174  }
1175 
1176  }
1177 
1178  // Compute the scale factor for EMCal, as a function of centrality
1179  Double_t numerator = (trackPtSumEMCal + clusESumEMCal) / accEMCal;
1180  Double_t denominator = trackPtSumTPC / accTPC;
1181  Double_t scaleFactor = numerator / denominator;
1182  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
1183  fHistManager.FillTH2(histname, fCent, scaleFactor);
1184 
1185  // Compute delta pT for EMCal, as a function of centrality
1186  Double_t rho = jetCont->GetRhoVal();
1187  Double_t deltaPt = trackPtSumEMCalRC + clusESumEMCalRC - rho * TMath::Pi() * jetR * jetR;
1188  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
1189  fHistManager.FillTH2(histname, fCent, deltaPt);
1190 
1191  delete r;
1192 
1193  }
1194 
1195 }
1196 
1202 {
1203  TString histname;
1204  AliJetContainer* jets = 0;
1205  TIter nextJetColl(&fJetCollArray);
1206  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1207  TString jetContName = jets->GetName();
1208 
1209  Double_t rhoVal = 0;
1210  if (jets->GetRhoParameter()) {
1211  rhoVal = jets->GetRhoVal();
1212  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
1213  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1214  }
1215 
1216  for (auto jet : jets->all()) {
1217 
1218  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1219  Float_t corrPt = GetJetPt(jet, rhoVal);
1220 
1221  // A vs. pT (fill before area cut)
1222  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
1223  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
1224 
1225  // Rejection reason
1226  UInt_t rejectionReason = 0;
1227  if (!jets->AcceptJet(jet, rejectionReason)) {
1228  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1229  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1230  continue;
1231  }
1232 
1233  // compute jet acceptance type
1234  Double_t type = GetJetType(jet);
1235 
1236  // Centrality vs. pT
1237  histname = TString::Format("%s/TriggerSimHistograms/hCentVsPt", jets->GetArrayName().Data());
1238  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type);
1239 
1240  // (Centrality, eta, phi, pT, NEF, calo type)
1241  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
1242  Double_t x[6] = {fCent, jet->Eta(), jet->Phi_0_2pi(), corrPt, jet->NEF(), type};
1243  fHistManager.FillTHnSparse(histname, x);
1244 
1245  // pT-leading vs. pT
1246  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1247  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1248 
1249  // (Centrality, pT, z-leading (charged), calo type)
1250  TLorentzVector leadPart;
1251  jets->GetLeadingHadronMomentum(leadPart, jet);
1252  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1253  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
1254  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
1255  Double_t y[4] = {fCent, corrPt, z, type};
1256  fHistManager.FillTHnSparse(histname, y);
1257 
1258  // (Centrality, pT, z (charged), calo type)
1259  histname = TString::Format("%s/TriggerSimHistograms/hZVsPt", jets->GetArrayName().Data());
1260  AliVTrack* track;
1261  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1262  track = static_cast<AliVTrack*>(jet->Track(i));
1263  z = track->Pt() / TMath::Abs(corrPt);
1264  Double_t y2[4] = {fCent, corrPt, z, type};
1265  fHistManager.FillTHnSparse(histname, y2);
1266  }
1267 
1268  // (Centrality, pT, Nconst, calo type)
1269  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPt", jets->GetArrayName().Data());
1270  Double_t a[4] = {fCent, corrPt, 1.*jet->GetNumberOfConstituents(), type};
1271  fHistManager.FillTHnSparse(histname, a);
1272 
1273  } //jet loop
1274 
1275  }
1276 }
1277 
1278 /*
1279  * Compute the MC particle type using AliMCAnalysisUtils
1280  */
1281 AliAnalysisTaskEmcalJetPerformance::ParticleType AliAnalysisTaskEmcalJetPerformance::GetParticleType1(const AliVCluster* clus, const AliMCEvent* mcevent, const TClonesArray* clusArray)
1282 {
1283  ParticleType particleType = kUndefined;
1284 
1285  AliMCAnalysisUtils mcUtils;
1286  Int_t tag = mcUtils.CheckOrigin(clus->GetLabels(), clus->GetNLabels(), mcevent, clusArray);
1287 
1288  Bool_t isPhoton = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton);
1289  Bool_t isPi0 = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0);
1290  Bool_t isConversion = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion);
1291  Bool_t isEta = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCEta);
1292  Bool_t isPion = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPion);
1293  Bool_t isKaon = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCKaon);
1294  Bool_t isProton = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCProton);
1295  Bool_t isAntiProton = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCAntiProton);
1296  Bool_t isNeutron = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCNeutron);
1297  Bool_t isAntiNeutron = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCAntiNeutron);
1298  Bool_t isElectron = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron);
1299  Bool_t isMuon = mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCMuon);
1300 
1301  if (isPi0) {
1302  if (isConversion) {
1303  particleType = kPi0Conversion;
1304  }
1305  else {
1306  particleType = kPi0;
1307  }
1308  }
1309  else if (isEta) {
1310  particleType = kEta;
1311  }
1312  else if (isPhoton) {
1313  particleType = kPhoton;
1314  }
1315  else if (isPion || isKaon || isProton || isAntiProton || isNeutron || isAntiNeutron) {
1316  particleType = kHadron;
1317  }
1318  else if (isElectron) {
1319  particleType = kElectron;
1320  }
1321  else if (isMuon) {
1322  particleType = kMuon;
1323  }
1324  else {
1325  particleType = kOther;
1326  }
1327  return particleType;
1328 }
1329 
1330 /*
1331  * Compute the MC particle type using the MC particle container (and only AliMCAnalysisUtils to find merged pi0)
1332  */
1334 {
1335  ParticleType particleType = kUndefined;
1336 
1337  AliAODMCParticle *part = fGeneratorLevel->GetMCParticleWithLabel(label);
1338  if (part) {
1339 
1340  TString histname = TString::Format("%s/hClusterRejectionReasonMC", clusters->GetArrayName().Data());
1341  UInt_t rejectionReason = 0;
1342  if (!fGeneratorLevel->AcceptMCParticle(part, rejectionReason)) {
1343  fHistManager.FillTH2(histname, fGeneratorLevel->GetRejectionReasonBitPosition(rejectionReason), clus->GetNonLinCorrEnergy());
1344  return particleType;
1345  }
1346 
1347  if (part->GetGeneratorIndex() == 0) { // generator index in cocktail
1348 
1349  // select charged pions, protons, kaons, electrons, muons
1350  Int_t pdg = TMath::Abs(part->PdgCode()); // abs value ensures both particles and antiparticles are included
1351 
1352  if (pdg == 22) { // gamma 22
1353 
1354  AliMCAnalysisUtils mcUtils;
1355  Int_t tag;
1356  mcUtils.CheckOverlapped2GammaDecay(clus->GetLabels(), clus->GetNLabels(), part->GetMother(), mcevent, tag);
1357 
1358  if (mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)) {
1359  if (mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)) {
1360  particleType = kPi0Conversion;
1361  }
1362  else {
1363  particleType = kPi0;
1364  }
1365  }
1366  else if (mcUtils.CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)) {
1367  particleType = kEta;
1368  }
1369  else { // direct photon
1370  particleType = kPhoton;
1371  }
1372 
1373  }
1374  else if (pdg == 211 || 2212 || 321 || 2112) { // pi+ 211, proton 2212, K+ 321, neutron 2112
1375  particleType = kHadron;
1376  }
1377  else if (pdg == 11) { // e- 11
1378  particleType = kElectron;
1379  }
1380  else if (pdg == 13) { // mu- 13
1381  particleType = kMuon;
1382  }
1383  else {
1384  particleType = kOther;
1385  }
1386  }
1387  }
1388  return particleType;
1389 }
1390 
1395 {
1396  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
1397  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
1398  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1399  return deltaR;
1400 }
1401 
1406 {
1407  UInt_t jetType = jet->GetJetAcceptanceType();
1408  Double_t type = -1;
1409  if (jetType & AliEmcalJet::kEMCAL) {
1410  type = kEMCal;
1411  }
1412  else if (jetType & AliEmcalJet::kDCALonly) {
1413  type = kDCal;
1414  }
1415 
1416  return type;
1417 }
1418 
1423 {
1424  Double_t pT = jet->Pt() - rho * jet->Area();
1425  return pT;
1426 }
Int_t pdg
Double_t Area() const
Definition: AliEmcalJet.h:123
ParticleType GetParticleType2(const AliVCluster *clus, const AliMCEvent *mcevent, Int_t label, const AliClusterContainer *clusters)
TObjArray fClusterCollArray
cluster collection array
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
UInt_t fOffTrigger
offline trigger for event selection
static AliEmcalDownscaleFactorsOCDB * Instance()
const AliClusterIterableMomentumContainer accepted_momentum() const
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Declaration of class AliTLorentzVector.
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
void SetCaloTriggerPatchInfoName(const char *n)
Bool_t fPlotClusterHistograms
Set whether to plot cluster histograms.
Int_t CheckOrigin(Int_t label, const AliMCEvent *mcevent)
Double_t GetDeltaR(const AliTLorentzVector *part, Double_t etaRef, Double_t phiRef)
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
ParticleType GetParticleType1(const AliVCluster *clus, const AliMCEvent *mcevent, const TClonesArray *clusArray)
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:262
Double_t Phi_0_2pi() const
EMCal acceptance.
Definition: AliEmcalJet.h:62
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.
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
Bool_t fPlotJetHistograms
Set whether to enable inclusive jet histograms.
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
TObjArray fJetCollArray
jet collection array
Bool_t fkEMCEJE
! flag telling whether the event is "triggered" or not in "simulation"
AliRhoParameter * GetRhoParameter()
void CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex, const AliMCEvent *mcevent, Int_t &tag)
Double_t Pt() const
Definition: AliEmcalJet.h:102
virtual Bool_t IsEventSelected()
Performing event selection.
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
Bool_t fComputeBackground
Set whether to enable study of background.
Pi0 (merged pi0) with conversion of one photon (may be only partially contained in cluster) ...
AliEmcalList * fOutput
!output list
AliMCParticleContainer * GetMCParticleContainer(Int_t i=0) const
Int_t fNPtHistBins
! number of variable pt bins
Definition: External.C:220
Handler for downscale factors for various triggers obtained from the OCDB.
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:44
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 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
DCal acceptance – spans ONLY DCal (no PHOS or gap)
Definition: AliEmcalJet.h:66
Double_t fMedianDCal
! median patch energy in DCal, per event
Class with analysis utils for simulations.
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.
virtual AliAODMCParticle * GetMCParticleWithLabel(Int_t lab) const
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Bool_t CheckTagBit(Int_t tag, UInt_t test) const
Container for jet within the EMCAL jet framework.
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()
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal