AliPhysics  7c67529 (7c67529)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalDijetImbalance.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <vector>
17 
18 #include <TClonesArray.h>
19 #include <TH1F.h>
20 #include <TH2F.h>
21 #include <TH3.h>
22 #include <TList.h>
23 #include <THnSparse.h>
24 #include <TRandom3.h>
25 #include <TGrid.h>
26 #include <TFile.h>
27 
28 #include <AliVCluster.h>
29 #include <AliVParticle.h>
30 #include <AliLog.h>
31 
32 #include "AliTLorentzVector.h"
33 #include "AliEmcalJet.h"
34 #include "AliRhoParameter.h"
35 #include "AliJetContainer.h"
36 #include "AliParticleContainer.h"
37 #include "AliClusterContainer.h"
38 #include "AliEMCALGeometry.h"
40 #include "AliPHOSGeometry.h"
41 #include "AliOADBContainer.h"
42 #include "AliEMCALTriggerPatchInfo.h"
44 
46 
50 
56  fHistManager(),
57  fEventCuts(0),
58  fEventCutList(0),
59  fUseManualEventCuts(kFALSE),
60  fUseAliEventCuts(kTRUE),
61  fDeltaPhiMin(0),
62  fMinTrigJetPt(0),
63  fMinAssJetPt(0),
64  fDijetLeadingHadronPt(0),
65  fMaxPt(200),
66  fNCentHistBins(0),
67  fCentHistBins(0),
68  fNPtHistBins(0),
69  fPtHistBins(0),
70  fPlotJetHistograms(kFALSE),
71  fPlotDijetCandHistograms(kFALSE),
72  fPlotDijetImbalanceHistograms(kFALSE),
73  fComputeBackground(kFALSE),
74  fDoMomentumBalance(kFALSE),
75  fDoGeometricalMatching(kFALSE),
76  fMatchingJetR(0.2),
77  fTrackConstituentThreshold(0),
78  fClusterConstituentThreshold(0),
79  fMBUpscaleFactor(1.),
80  fNEtaBins(40),
81  fNPhiBins(200),
82  fLoadBackgroundScalingWeights(kTRUE),
83  fBackgroundScalingWeights(0),
84  fGapJetScalingWeights(0),
85  fComputeMBDownscaling(kFALSE),
86  fDoCaloStudy(kFALSE),
87  fDoTriggerSimulation(kFALSE),
88  fMedianEMCal(0),
89  fMedianDCal(0),
90  fkEMCEJE(kFALSE),
91  fPlotNeutralJets(kFALSE),
92  fPlotClustersInJets(kFALSE),
93  fPlotClusterTHnSparse(kTRUE),
94  fPlotClusWithoutNonLinCorr(kFALSE),
95  fPlotExotics(kFALSE),
96  fEmbeddingQA(),
97  fPHOSGeo(nullptr)
98 {
100  Dijet_t fDijet;
102 }
103 
110  AliAnalysisTaskEmcalJet(name, kTRUE),
111  fHistManager(name),
112  fEventCuts(0),
113  fEventCutList(0),
114  fUseManualEventCuts(kFALSE),
115  fUseAliEventCuts(kTRUE),
116  fDeltaPhiMin(0),
117  fMinTrigJetPt(0),
118  fMinAssJetPt(0),
119  fDijetLeadingHadronPt(0),
120  fMaxPt(200),
121  fNCentHistBins(0),
122  fCentHistBins(0),
123  fNPtHistBins(0),
124  fPtHistBins(0),
125  fPlotJetHistograms(kFALSE),
126  fPlotDijetCandHistograms(kFALSE),
127  fPlotDijetImbalanceHistograms(kFALSE),
128  fComputeBackground(kFALSE),
129  fDoMomentumBalance(kFALSE),
130  fDoGeometricalMatching(kFALSE),
131  fMatchingJetR(0.2),
132  fTrackConstituentThreshold(0),
133  fClusterConstituentThreshold(0),
134  fMBUpscaleFactor(1.),
135  fNEtaBins(40),
136  fNPhiBins(200),
137  fLoadBackgroundScalingWeights(kTRUE),
138  fBackgroundScalingWeights(0),
139  fGapJetScalingWeights(0),
140  fComputeMBDownscaling(kFALSE),
141  fDoCaloStudy(kFALSE),
142  fDoTriggerSimulation(kFALSE),
143  fMedianEMCal(0),
144  fMedianDCal(0),
145  fkEMCEJE(kFALSE),
146  fPlotNeutralJets(kFALSE),
147  fPlotClustersInJets(kFALSE),
148  fPlotClusterTHnSparse(kTRUE),
149  fPlotClusWithoutNonLinCorr(kFALSE),
150  fPlotExotics(kFALSE),
151  fEmbeddingQA(),
152  fPHOSGeo(nullptr)
153 {
155  Dijet_t fDijet;
157 }
158 
163 {
164 }
165 
170 {
171  fNCentHistBins = 4;
173  fCentHistBins[0] = 0;
174  fCentHistBins[1] = 10;
175  fCentHistBins[2] = 30;
176  fCentHistBins[3] = 50;
177  fCentHistBins[4] = 90;
178 
179  fNPtHistBins = 82;
182  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
183  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
184  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
185  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
186  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
187  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
188 }
189 
195 {
197 
205 
206  TIter next(fHistManager.GetListOfHistograms());
207  TObject* obj = 0;
208  while ((obj = next())) {
209  fOutput->Add(obj);
210  }
211 
212  // Intialize AliEventCuts
213  if (fUseAliEventCuts) {
214  fEventCutList = new TList();
215  fEventCutList ->SetOwner();
216  fEventCutList ->SetName("EventCutOutput");
217 
218  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
219  if(fUseManualEventCuts==1)
220  {
221  fEventCuts.SetManualMode();
222  // Configure manual settings here
223  // ...
224  }
225  fEventCuts.AddQAplotsToList(fEventCutList);
226  fOutput->Add(fEventCutList);
227  }
228 
229  // Load eta-phi background scale factors from histogram on AliEn
232  }
233 
234  // Initialize embedding QA
236  if (embeddingHelper) {
237  bool res = fEmbeddingQA.Initialize();
238  if (res) {
240  }
241  }
242 
243  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
244 }
245 
246 /*
247  * This function allocates the histograms for single jets.
248  * A set of histograms is allocated per each jet container.
249  */
251 {
252  TString histname;
253  TString title;
254 
255  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
256 
257  AliJetContainer* jets = 0;
258  TIter nextJetColl(&fJetCollArray);
259  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
260 
261  // Jet rejection reason
262  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
263  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
264  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
265  SetRejectionReasonLabels(hist->GetXaxis());
266 
267  // Rho vs. Centrality
268  if (!jets->GetRhoName().IsNull()) {
269  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
270  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
271  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
272  }
273 
274  // (Centrality, pT, calo type)
275  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
276  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
277  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, 3, -0.5, 2.5);
278 
279  // (Centrality, eta, phi, pT, NEF, calo type)
280  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
281  title = histname + ";Centrality (%);#eta_{jet} (rad);#phi_{jet} (rad);#it{p}_{T}^{corr} (GeV/#it{c});NEF;type";
282  Int_t nbins1[6] = {20, fNEtaBins, fNPhiBins, nPtBins, 50, 3};
283  Double_t min1[6] = {0, -0.5,1., 0, 0, -0.5};
284  Double_t max1[6] = {100, 0.5,6., fMaxPt, 1., 2.5};
285  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 6, nbins1, min1, max1);
286 
287  // (Centrality, pT upscaled, calo type)
288  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
289  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
290  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, 3, -0.5, 2.5, "s");
291 
292  // pT-leading vs. pT
293  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
294  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
295  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
296 
297  // A vs. pT
298  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
299  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
300  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/3, 0, 0.5);
301 
302  // (Centrality, pT, z-leading (charged), calo type)
303  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
304  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading};type";
305  Int_t nbins2[4] = {20, nPtBins, 50, 3};
306  Double_t min2[4] = {0, 0, 0, -0.5};
307  Double_t max2[4] = {100, fMaxPt, 1., 2.5};
308  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins2, min2, max2);
309 
310  // (Centrality, pT, z (charged), calo type)
311  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
312  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z};type";
313  Int_t nbins3[4] = {20, nPtBins, 50, 3};
314  Double_t min3[4] = {0, 0, 0, -0.5};
315  Double_t max3[4] = {100, fMaxPt, 1., 2.5};
316  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins3, min3, max3);
317 
318  // (Centrality, pT, Nconst, calo type)
319  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
320  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents;type";
321  Int_t nbins4[4] = {20, nPtBins, 50, 3};
322  Double_t min4[4] = {0, 0, 0, -0.5};
323  Double_t max4[4] = {100, fMaxPt, fMaxPt, 2.5};
324  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins4, min4, max4);
325 
326  // (Median patch energy, calo type, jet pT, centrality)
327  if (fDoTriggerSimulation) {
328  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
329  title = histname + ";#it{E}_{patch,med};type;#it{p}_{T}^{corr} (GeV/#it{c});Centrality (%)";
330  Int_t nbins5[4] = {100, 2, nPtBins, 50};
331  Double_t min5[4] = {0,-0.5, 0, 0};
332  Double_t max5[4] = {50,1.5, fMaxPt, 100};
333  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins5, min5, max5);
334  }
335 
336  // Allocate background subtraction histograms, if enabled
337  if (fComputeBackground) {
338 
339  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
340  title = histname + ";Centrality;Scale factor;counts";
341  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
342 
343  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
344  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
345  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
346 
347  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jets->GetArrayName().Data());
348  title = histname + ";#eta;#phi;Centrality;Scale factor;";
349  Int_t nbins[4] = {fNEtaBins, fNPhiBins, 50, 400};
350  Double_t min[4] = {-0.5,1., 0, 0};
351  Double_t max[4] = {0.5,6., 100, 20};
352  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins, min, max);
353 
354  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jets->GetArrayName().Data());
355  title = histname + ";#eta;#phi;Centrality;#delta#it{p}_{T} (GeV/#it{c})";
356  Int_t nbinsDpT[4] = {fNEtaBins, fNPhiBins, 10, 400};
357  Double_t minDpT[4] = {-0.5,1., 0, -50};
358  Double_t maxDpT[4] = {0.5,6., 100, 150};
359  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbinsDpT, minDpT, maxDpT);
360 
361  }
362 
363  }
364 
365  // MB downscale factor histogram
366  if (fComputeMBDownscaling) {
367  histname = "Trigger/hMBDownscaleFactor";
368  title = histname + ";Downscale factor;counts";
369  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 200);
370  }
371 
372 }
373 
374 /*
375  * This function allocates the histograms for dijet candidates, i.e. dijet pairs with the leading jet
376  * passing the trigger condition, but no condition on the subleading jet. In particular this histogram is
377  * designed to study the kinematic selection of dijets.
378  */
380 {
381  // Allocate dijet THnSparse
382  AliJetContainer* jets = 0;
383  TIter nextJetColl(&fJetCollArray);
384  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
385  TString axisTitle[30]= {""};
386  Int_t nbins[30] = {0};
387  Double_t min[30] = {0.};
388  Double_t max[30] = {0.};
389  Double_t *binEdges[20] = {0};
390  Int_t dim = 0;
391 
392  if (fForceBeamType != kpp) {
393  axisTitle[dim] = "Centrality (%)";
394  nbins[dim] = fNCentHistBins;
395  binEdges[dim] = fCentHistBins;
396  min[dim] = fCentHistBins[0];
397  max[dim] = fCentHistBins[fNCentHistBins];
398  dim++;
399  }
400 
401  axisTitle[dim] = "LeadingHadronRequired";
402  nbins[dim] = 2;
403  min[dim] = -0.5;
404  max[dim] = 1.5;
405  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
406  dim++;
407 
408  axisTitle[dim] = "#it{p}_{T,trig jet} (GeV/#it{c})";
409  nbins[dim] = TMath::CeilNint(fMaxPt/2);
410  min[dim] = 0;
411  max[dim] = fMaxPt;
412  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
413  dim++;
414 
415  axisTitle[dim] = "#it{p}_{T,ass jet} (GeV/#it{c})";
416  nbins[dim] = TMath::CeilNint(fMaxPt/2);
417  min[dim] = 0;
418  max[dim] = fMaxPt;
419  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
420  dim++;
421 
422  axisTitle[dim] = "#phi_{trig jet}";
423  nbins[dim] = TMath::CeilNint(fMaxPt/2);
424  min[dim] = 0;
425  max[dim] = TMath::TwoPi();
426  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
427  dim++;
428 
429  axisTitle[dim] = "#phi_{ass jet}";
430  nbins[dim] = TMath::CeilNint(fMaxPt/2);
431  min[dim] = 0;
432  max[dim] = TMath::TwoPi();
433  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
434  dim++;
435 
436  axisTitle[dim] = "#eta_{trig jet}";
437  nbins[dim] = TMath::CeilNint(fMaxPt/2);
438  min[dim] = -1;
439  max[dim] = 1;
440  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
441  dim++;
442 
443  axisTitle[dim] = "#eta_{ass jet}";
444  nbins[dim] = TMath::CeilNint(fMaxPt/2);
445  min[dim] = -1;
446  max[dim] = 1;
447  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
448  dim++;
449 
450  TString thnname = TString::Format("%s/DijetCandObservables", jets->GetArrayName().Data());
451  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
452  for (Int_t i = 0; i < dim; i++) {
453  hn->GetAxis(i)->SetTitle(axisTitle[i]);
454  hn->SetBinEdges(i, binEdges[i]);
455  }
456  }
457 }
458 
459 /*
460  * This function allocates the histograms for accepted dijets.
461  * The purpose is to study in detail the imbalance properties of the accepted dijets.
462  */
464 {
465  AliJetContainer* jets = 0;
466  TIter nextJetColl(&fJetCollArray);
467  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
468 
469  // Allocate dijet imbalance THnSparse, unless doing trigger simulation
470  if (!fDoTriggerSimulation) {
471  TString axisTitle[30]= {""};
472  Int_t nbins[30] = {0};
473  Double_t min[30] = {0.};
474  Double_t max[30] = {0.};
475  Double_t *binEdges[20] = {0};
476  Int_t dim = 0;
477 
478  if (fForceBeamType != kpp) {
479  axisTitle[dim] = "Centrality (%)";
480  nbins[dim] = fNCentHistBins;
481  binEdges[dim] = fCentHistBins;
482  min[dim] = fCentHistBins[0];
483  max[dim] = fCentHistBins[fNCentHistBins];
484  dim++;
485  }
486 
487  axisTitle[dim] = "#Delta#phi";
488  nbins[dim] = 100;
489  min[dim] = 0;
490  max[dim] = 4;
491  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
492  dim++;
493 
494  axisTitle[dim] = "#Delta#eta";
495  nbins[dim] = 100;
496  min[dim] = -2;
497  max[dim] = 2;
498  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
499  dim++;
500 
501  axisTitle[dim] = "A_{J}";
502  nbins[dim] = 100;
503  min[dim] = 0;
504  max[dim] = 1;
505  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
506  dim++;
507 
508  axisTitle[dim] = "x_{J}";
509  nbins[dim] = 100;
510  min[dim] = 0;
511  max[dim] = 1;
512  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
513  dim++;
514 
515  axisTitle[dim] = "k_{Ty} (GeV)";
516  nbins[dim] = 100;
517  min[dim] = 0;
518  max[dim] = 100;
519  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
520  dim++;
521 
522  axisTitle[dim] = "N_{tracks, trig jet}";
523  nbins[dim] = fMaxPt/5;
524  min[dim] = 0;
525  max[dim] = 100;
526  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
527  dim++;
528 
529  axisTitle[dim] = "N_{tracks, ass jet}";
530  nbins[dim] = fMaxPt/5;
531  min[dim] = 0;
532  max[dim] = 100;
533  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
534  dim++;
535 
536  TString thnname = TString::Format("%s/DijetImbalanceObservables", jets->GetArrayName().Data());
537  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
538  for (Int_t i = 0; i < dim; i++) {
539  hn->GetAxis(i)->SetTitle(axisTitle[i]);
540  hn->SetBinEdges(i, binEdges[i]);
541  }
542  }
543 
544  // Now, allocate 2d pt spectrum, upscaled according to MB downscaling factors (for trigger efficiency)
545  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
546 
547  // (Centrality, pT1, pT2) (upscaled)
548  TString histname = TString::Format("%s/DijetJetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
549  TString title = histname + ";Centrality (%);#it{p}_{T,1}^{corr} (GeV/#it{c});#it{p}_{T,2}^{corr} (GeV/#it{c})";
550  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt, "s");
551 
552  // Now, plot jet histograms for the leading jet within the dijet (for comparison to single jets, and triggered to MB)
553 
554  // (Centrality, pT, NEF, calo type)
555  histname = TString::Format("%s/DijetJetHistograms/hNEFVsPt", jets->GetArrayName().Data());
556  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF;type";
557  Int_t nbins1[4] = {20, nPtBins, 50, 3};
558  Double_t min1[4] = {0, 0, 0, -0.5};
559  Double_t max1[4] = {100, fMaxPt, 1., 2.5};
560  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins1, min1, max1);
561 
562  // (Centrality, pT, z-leading (charged), calo type)
563  histname = TString::Format("%s/DijetJetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
564  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading};type";
565  Int_t nbins2[4] = {20, nPtBins, 50, 3};
566  Double_t min2[4] = {0, 0, 0, -0.5};
567  Double_t max2[4] = {100, fMaxPt, 1., 2.5};
568  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins2, min2, max2);
569 
570  // (Centrality, pT, z (charged), calo type)
571  histname = TString::Format("%s/DijetJetHistograms/hZVsPt", jets->GetArrayName().Data());
572  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z};type";
573  Int_t nbins3[4] = {20, nPtBins, 50, 3};
574  Double_t min3[4] = {0, 0, 0, -0.5};
575  Double_t max3[4] = {100, fMaxPt, 1., 2.5};
576  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins3, min3, max3);
577 
578  // (Centrality, pT, Nconst, calo type)
579  histname = TString::Format("%s/DijetJetHistograms/hNConstVsPt", jets->GetArrayName().Data());
580  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents;type";
581  Int_t nbins4[4] = {20, nPtBins, 50, 3};
582  Double_t min4[4] = {0, 0, 0, -0.5};
583  Double_t max4[4] = {100, fMaxPt, fMaxPt, 2.5};
584  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins4, min4, max4);
585  }
586 }
587 
588 /*
589  * This function allocates the histograms for the momentum balance study.
590  */
592 {
593  AliJetContainer* jets = 0;
594  TIter nextJetColl(&fJetCollArray);
595  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
596 
597  // Allocate THnSparse
598  TString axisTitle[30]= {""};
599  Int_t nbins[30] = {0};
600  Double_t min[30] = {0.};
601  Double_t max[30] = {0.};
602  Double_t *binEdges[20] = {0};
603  Int_t dim = 0;
604 
605  axisTitle[dim] = "A_{J}";
606  nbins[dim] = 100;
607  min[dim] = 0;
608  max[dim] = 1;
609  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
610  dim++;
611 
612  axisTitle[dim] = "#Delta#phi";
613  nbins[dim] = 100;
614  min[dim] = -4;
615  max[dim] = 4;
616  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
617  dim++;
618 
619  axisTitle[dim] = "#it{p}_{T,particle} (GeV/#it{c})";
620  nbins[dim] = 9;
621  Double_t* pTParticleBins = new Double_t[nbins[dim]+1];
622  GenerateFixedBinArray(1, 0.15, 0.3, pTParticleBins);
623  GenerateFixedBinArray(1, 0.3, 0.5, pTParticleBins+1);
624  GenerateFixedBinArray(1, 0.5, 1, pTParticleBins+2);
625  GenerateFixedBinArray(2, 1, 5, pTParticleBins+3);
626  GenerateFixedBinArray(3, 5, 20, pTParticleBins+5);
627  GenerateFixedBinArray(1, 20, 150, pTParticleBins+8);
628  min[dim] = 0;
629  max[dim] = pTParticleBins[nbins[dim]];
630  binEdges[dim] = pTParticleBins;
631  dim++;
632 
633  axisTitle[dim] = "#it{p}_{T}#parallel (GeV/#it{c})";
634  nbins[dim] = 80;
635  Double_t* pTParallelBins = new Double_t[nbins[dim]+1];
636  GenerateFixedBinArray(20, 0, 2, pTParallelBins);
637  GenerateFixedBinArray(16, 2, 10, pTParallelBins+20);
638  GenerateFixedBinArray(10, 10, 20, pTParallelBins+36);
639  GenerateFixedBinArray(10, 20, 40, pTParallelBins+46);
640  GenerateFixedBinArray(24, 40, 150, pTParallelBins+56);
641  min[dim] = 0;
642  max[dim] = pTParallelBins[nbins[dim]];
643  binEdges[dim] = pTParallelBins;
644  dim++;
645 
646  TString thnname = TString::Format("%s/MomentumBalance", jets->GetArrayName().Data());
647  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
648  for (Int_t i = 0; i < dim; i++) {
649  hn->GetAxis(i)->SetTitle(axisTitle[i]);
650  hn->SetBinEdges(i, binEdges[i]);
651  }
652  }
653 }
654 
655 /*
656  * This function allocates the histograms for the constituent dijet study.
657  */
659 {
660  // Allocate geometrical matching THnSparse
661  TString axisTitle[30]= {""};
662  Int_t nbins[30] = {0};
663  Double_t min[30] = {0.};
664  Double_t max[30] = {0.};
665  Double_t *binEdges[20] = {0};
666  Int_t dim = 0;
667 
668  if (fForceBeamType != kpp) {
669  axisTitle[dim] = "Centrality (%)";
670  nbins[dim] = fNCentHistBins;
671  binEdges[dim] = fCentHistBins;
672  min[dim] = fCentHistBins[0];
673  max[dim] = fCentHistBins[fNCentHistBins];
674  dim++;
675  }
676 
677  axisTitle[dim] = "isSwitched";
678  nbins[dim] = 2;
679  min[dim] = -0.5;
680  max[dim] = 1.5;
681  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
682  dim++;
683 
684  axisTitle[dim] = "#DeltaR_{trig}";
685  nbins[dim] = 50;
686  min[dim] = 0;
687  max[dim] = 0.5;
688  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
689  dim++;
690 
691  axisTitle[dim] = "#DeltaR_{ass}";
692  nbins[dim] = 50;
693  min[dim] = 0;
694  max[dim] = 0.5;
695  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
696  dim++;
697 
698  axisTitle[dim] = "trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
699  nbins[dim] = 100;
700  min[dim] = -50;
701  max[dim] = 50;
702  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
703  dim++;
704 
705  axisTitle[dim] = "ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
706  nbins[dim] = 100;
707  min[dim] = -50;
708  max[dim] = 50;
709  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
710  dim++;
711 
712  axisTitle[dim] = "A_{J} low-threshold";
713  nbins[dim] = 100;
714  min[dim] = 0;
715  max[dim] = 1;
716  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
717  dim++;
718 
719  axisTitle[dim] = "A_{J} hard-core";
720  nbins[dim] = 100;
721  min[dim] = 0;
722  max[dim] = 1;
723  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
724  dim++;
725 
726  TString thnname = "GeometricalMatching";
727  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
728  for (Int_t i = 0; i < dim; i++) {
729  hn->GetAxis(i)->SetTitle(axisTitle[i]);
730  hn->SetBinEdges(i, binEdges[i]);
731  }
732 
733  // Allocate other histograms
734  TString histname;
735  TString title;
736  histname = "GeometricalMatchingEfficiency";
737  title = histname + ";isMatched;counts";
738  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 2, -0.5, 1.5);
739 }
740 
741 /*
742  * This function allocates the histograms for the calorimeter performance study.
743  */
745 {
746  TString histname;
747  TString htitle;
748  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
749 
750  Double_t* clusType = new Double_t[3+1];
751  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
752  const Int_t nRejBins = 32;
753  Double_t* rejReasonBins = new Double_t[nRejBins+1];
754  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
755  const Int_t nExBins = 200;
756  Double_t* exBins = new Double_t[nExBins+1];
757  GenerateFixedBinArray(nExBins, 0, 1, exBins);
758 
759  AliEmcalContainer* cont = 0;
760  TIter nextClusColl(&fClusterCollArray);
761  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
762 
763  // rejection reason plot, to make efficiency correction
764  histname = TString::Format("%s/hClusterRejectionReasonEMCal", cont->GetArrayName().Data());
765  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
766  TH2* hist = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
767  SetRejectionReasonLabels(hist->GetXaxis());
768 
769  histname = TString::Format("%s/hClusterRejectionReasonPHOS", cont->GetArrayName().Data());
770  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
771  TH2* histPhos = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
772  SetRejectionReasonLabels(histPhos->GetXaxis());
773 
774  // plot by SM
775  const Int_t nEmcalSM = 20;
776  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
777  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
778  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
779  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
780  }
781 
782  for (Int_t sm = 1; sm < 5; sm++) {
783  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
784  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
785  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
786  }
787 
788  // Plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
789  if (fPlotClusterTHnSparse) {
790  Int_t dim = 0;
791  TString title[20];
792  Int_t nbins[20] = {0};
793  Double_t min[30] = {0.};
794  Double_t max[30] = {0.};
795  Double_t *binEdges[20] = {0};
796 
798  title[dim] = "Centrality %";
799  nbins[dim] = fNCentHistBins;
800  binEdges[dim] = fCentHistBins;
801  min[dim] = fCentHistBins[0];
802  max[dim] = fCentHistBins[fNCentHistBins];
803  dim++;
804  }
805 
806  title[dim] = "#eta";
807  nbins[dim] = 28;
808  min[dim] = -0.7;
809  max[dim] = 0.7;
810  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
811  dim++;
812 
813  title[dim] = "#phi";
814  nbins[dim] = 100;
815  min[dim] = 1.;
816  max[dim] = 6.;
817  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
818  dim++;
819 
820  title[dim] = "#it{E}_{clus} (GeV)";
821  nbins[dim] = fNPtHistBins;
822  binEdges[dim] = fPtHistBins;
823  min[dim] = fPtHistBins[0];
824  max[dim] = fPtHistBins[fNPtHistBins];
825  dim++;
826 
827  title[dim] = "#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)";
828  nbins[dim] = fNPtHistBins;
829  binEdges[dim] = fPtHistBins;
830  min[dim] = fPtHistBins[0];
831  max[dim] = fPtHistBins[fNPtHistBins];
832  dim++;
833 
834  title[dim] = "Matched track";
835  nbins[dim] = 2;
836  min[dim] = -0.5;
837  max[dim] = 1.5;
838  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
839  dim++;
840 
841  title[dim] = "M02";
842  nbins[dim] = 50;
843  min[dim] = 0;
844  max[dim] = 5;
845  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
846  dim++;
847 
848  title[dim] = "Ncells";
849  nbins[dim] = 30;
850  min[dim] = -0.5;
851  max[dim] = 29.5;
852  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
853  dim++;
854 
855  title[dim] = "Dispersion cut";
856  nbins[dim] = 2;
857  min[dim] = -0.5;
858  max[dim] = 1.5;
859  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
860  dim++;
861 
862  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
863  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
864  for (Int_t i = 0; i < dim; i++) {
865  hn->GetAxis(i)->SetTitle(title[i]);
866  hn->SetBinEdges(i, binEdges[i]);
867  }
868  }
869 
870  if (fPlotExotics) {
871  histname = TString::Format("%s/hFcrossEMCal", cont->GetArrayName().Data());
872  htitle = histname + ";Centrality (%);Fcross;#it{E}_{clus} (GeV/)";
873  TH3* hist = fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNCentHistBins, fCentHistBins, nExBins, exBins, fNPtHistBins, fPtHistBins);
874  }
875  }
876 
877  AliJetContainer* jets = 0;
878  TIter nextJetColl(&fJetCollArray);
879  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
880 
881  // plot neutral jets
882  if (fPlotNeutralJets) {
883  TString axisTitle[30]= {""};
884  Int_t nbinsJet[30] = {0};
885  Double_t minJet[30] = {0.};
886  Double_t maxJet[30] = {0.};
887  Double_t *binEdgesJet[20] = {0};
888  Int_t dimJet = 0;
889 
890  if (fForceBeamType != kpp) {
891  axisTitle[dimJet] = "Centrality (%)";
892  nbinsJet[dimJet] = fNCentHistBins;
893  binEdgesJet[dimJet] = fCentHistBins;
894  minJet[dimJet] = fCentHistBins[0];
895  maxJet[dimJet] = fCentHistBins[fNCentHistBins];
896  dimJet++;
897  }
898 
899  axisTitle[dimJet] = "#eta_{jet}";
900  nbinsJet[dimJet] = 28;
901  minJet[dimJet] = -0.7;
902  maxJet[dimJet] = 0.7;
903  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
904  dimJet++;
905 
906  axisTitle[dimJet] = "#phi_{jet} (rad)";
907  nbinsJet[dimJet] = 100;
908  minJet[dimJet] = 1.;
909  maxJet[dimJet] = 6.;
910  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
911  dimJet++;
912 
913  axisTitle[dimJet] = "#it{E}_{T} (GeV)";
914  nbinsJet[dimJet] = fNPtHistBins;
915  binEdgesJet[dimJet] = fPtHistBins;
916  minJet[dimJet] = fPtHistBins[0];
917  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
918  dimJet++;
919 
920  axisTitle[dimJet] = "#rho (GeV/#it{c})";
921  nbinsJet[dimJet] = 100;
922  minJet[dimJet] = 0.;
923  maxJet[dimJet] = 1000.;
924  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
925  dimJet++;
926 
927  axisTitle[dimJet] = "N_{clusters}";
928  nbinsJet[dimJet] = 20;
929  minJet[dimJet] = -0.5;
930  maxJet[dimJet] = 19.5;
931  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
932  dimJet++;
933 
934  TString thnname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
935  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dimJet, nbinsJet, minJet, maxJet);
936  for (Int_t i = 0; i < dimJet; i++) {
937  hn->GetAxis(i)->SetTitle(axisTitle[i]);
938  hn->SetBinEdges(i, binEdgesJet[i]);
939  }
940  }
941 
942  // Plot cluster spectra within jets
943  // (centrality, cluster energy, jet pT, jet eta, jet phi)
944  if (fPlotClustersInJets) {
945  Int_t dim = 0;
946  TString title[20];
947  Int_t nbins[20] = {0};
948  Double_t min[30] = {0.};
949  Double_t max[30] = {0.};
950  Double_t *binEdges[20] = {0};
951 
953  title[dim] = "Centrality %";
954  nbins[dim] = fNCentHistBins;
955  binEdges[dim] = fCentHistBins;
956  min[dim] = fCentHistBins[0];
957  max[dim] = fCentHistBins[fNCentHistBins];
958  dim++;
959  }
960 
961  title[dim] = "#it{E}_{clus} (GeV)";
962  nbins[dim] = fNPtHistBins;
963  binEdges[dim] = fPtHistBins;
964  min[dim] = fPtHistBins[0];
965  max[dim] = fPtHistBins[fNPtHistBins];
966  dim++;
967 
968  title[dim] = "#it{p}_{T,jet}^{corr}";
969  nbins[dim] = nPtBins;
970  min[dim] = 0;
971  max[dim] = fMaxPt;
972  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
973  dim++;
974 
975  title[dim] = "#eta_{jet}";
976  nbins[dim] = fNEtaBins;
977  min[dim] = -0.5;
978  max[dim] = 0.5;
979  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
980  dim++;
981 
982  title[dim] = "#phi_{jet}";
983  nbins[dim] = fNPhiBins;
984  min[dim] = 1.;
985  max[dim] = 6.;
986  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
987  dim++;
988 
989  TString thnname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
990  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
991  for (Int_t i = 0; i < dim; i++) {
992  hn->GetAxis(i)->SetTitle(title[i]);
993  hn->SetBinEdges(i, binEdges[i]);
994  }
995 
996  // (jet type, jet pT, cluster shift)
997  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
998  htitle = histname + ";type;#it{p}_{T}^{corr} (GeV/#it{c});#Delta_{JES}";
999  fHistManager.CreateTH3(histname.Data(), htitle.Data(), 3, -0.5, 2.5, nPtBins, 0, fMaxPt, 100, 0, 20);
1000  }
1001 
1002  }
1003 
1004  // Plot cell histograms
1005 
1006  // centrality vs. cell energy vs. cell type (for all cells)
1007  histname = TString::Format("Cells/hCellEnergyAll");
1008  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
1009  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
1010 
1011  // centrality vs. cell energy vs. cell type (for cells in accepted clusters)
1012  histname = TString::Format("Cells/hCellEnergyAccepted");
1013  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
1014  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
1015 
1016  // centrality vs. cell energy vs. cell type (for leading cells in accepted clusters)
1017  histname = TString::Format("Cells/hCellEnergyLeading");
1018  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
1019  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
1020 
1021  // plot cell patches by SM
1022  const Int_t nEmcalSM = 20;
1023  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
1024  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
1025  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
1026  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
1027  }
1028 
1029  for (Int_t sm = 1; sm < 5; sm++) {
1030  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
1031  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
1032  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
1033  }
1034 
1035 }
1036 
1037 /*
1038  * This function allocates the histograms for single jets, when the "simulated" trigger has been fired.
1039  * A set of histograms is allocated per each jet container.
1040  */
1042 {
1043  TString histname;
1044  TString title;
1045  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
1046 
1047  //----------------------------------------------
1048  // Trigger patch histograms
1049 
1050  // patch eta vs. phi
1051  histname = "TriggerSimHistograms/hEtaVsPhi";
1052  title = histname + ";#eta_{patch} (rad);#phi_{patch} (rad)";
1053  fHistManager.CreateTH2(histname.Data(), title.Data(), 140, -0.7, 0.7, 500, 1., 6.);
1054 
1055  // N patches
1056  histname = "TriggerSimHistograms/hNPatches";
1057  title = histname + ";#it{N}_{patches};type";
1058  fHistManager.CreateTH2(histname.Data(), title.Data(), 200, 0, 200, 2, -0.5, 1.5);
1059 
1060  // patch E vs. centrality
1061  histname = "TriggerSimHistograms/hPatchE";
1062  title = histname + ";Centrality (%);#it{E}_{patch} (GeV)";
1063  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, nPtBins, 0, fMaxPt);
1064 
1065  // patch median vs. Centrality
1066  histname = "TriggerSimHistograms/hPatchMedianE";
1067  title = histname + ";Centrality (%);#it{E}_{patch,med} (GeV);type";
1068  fHistManager.CreateTH3(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 50, 2, -0.5, 1.5);
1069 
1070  // Median patch energy vs. centrality, for dijets
1071  histname = "TriggerSimHistograms/hMedPatchDijet";
1072  title = histname + ";Centrality (%);#it{p}_{T,trig}^{corr} (GeV/#it{c});#it{E}_{patch,med} (GeV);type";
1073  Int_t nbinsD[4] = {50, 40, 100, 2};
1074  Double_t minD[4] = {0, 0, 0, -0.5};
1075  Double_t maxD[4] = {100, 200, 50, 1.5};
1076  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbinsD, minD, maxD);
1077 
1078  //----------------------------------------------
1079  // Jet histograms for "triggered" events
1080  AliJetContainer* jets = 0;
1081  TIter nextJetColl(&fJetCollArray);
1082  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1083 
1084  // Jet rejection reason
1085  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1086  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
1087  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
1088  SetRejectionReasonLabels(hist->GetXaxis());
1089 
1090  // Rho vs. Centrality
1091  if (!jets->GetRhoName().IsNull()) {
1092  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
1093  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
1094  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
1095  }
1096 
1097  // (Centrality, pT, calo type)
1098  histname = TString::Format("%s/TriggerSimHistograms/hCentVsPt", jets->GetArrayName().Data());
1099  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});type";
1100  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, 3, -0.5, 2.5);
1101 
1102  // (Centrality, eta, phi, pT, NEF, calo type)
1103  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
1104  title = histname + ";Centrality (%);#eta_{jet} (rad);#phi_{jet} (rad);#it{p}_{T}^{corr} (GeV/#it{c});NEF;type";
1105  Int_t nbins1[6] = {20, fNEtaBins, fNPhiBins, nPtBins, 50, 3};
1106  Double_t min1[6] = {0, -0.5,1., 0, 0, -0.5};
1107  Double_t max1[6] = {100, 0.5,6., fMaxPt, 1., 2.5};
1108  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 6, nbins1, min1, max1);
1109 
1110  // pT-leading vs. pT
1111  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1112  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
1113  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
1114 
1115  // A vs. pT
1116  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
1117  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
1118  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/3, 0, 0.5);
1119 
1120  // (Centrality, pT, z-leading (charged), calo type)
1121  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
1122  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading};type";
1123  Int_t nbins2[4] = {20, nPtBins, 50, 3};
1124  Double_t min2[4] = {0, 0, 0, -0.5};
1125  Double_t max2[4] = {100, fMaxPt, 1., 2.5};
1126  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins2, min2, max2);
1127 
1128  // z (charged) vs. pT
1129  histname = TString::Format("%s/TriggerSimHistograms/hZVsPt", jets->GetArrayName().Data());
1130  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z};type";
1131  Int_t nbins3[4] = {20, nPtBins, 50, 3};
1132  Double_t min3[4] = {0, 0, 0, -0.5};
1133  Double_t max3[4] = {100, fMaxPt, 1., 2.5};
1134  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins3, min3, max3);
1135 
1136  // (Centrality, pT, Nconst, calo type)
1137  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPt", jets->GetArrayName().Data());
1138  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents;type";
1139  Int_t nbins4[4] = {20, nPtBins, 50, 3};
1140  Double_t min4[4] = {0, 0, 0, -0.5};
1141  Double_t max4[4] = {100, fMaxPt, fMaxPt, 2.5};
1142  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins4, min4, max4);
1143 
1144  }
1145 
1146 }
1147 
1151 void AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram(const char* path, const char* name1, const char* name2)
1152 {
1153 
1154  TString fname(path);
1155  if (fname.BeginsWith("alien://")) {
1156  TGrid::Connect("alien://");
1157  }
1158 
1159  TFile* file = TFile::Open(path);
1160 
1161  if (!file || file->IsZombie()) {
1162  ::Error("AliAnalysisTaskEmcalDijetImbalance", "Could not open background scaling histogram");
1163  return;
1164  }
1165 
1166  // Open background scale factor histogram
1167  TH2D* h1 = dynamic_cast<TH2D*>(file->Get(name1));
1168 
1169  if (h1) {
1170  ::Info("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s loaded from file %s.", name1, path);
1171  }
1172  else {
1173  ::Error("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s not found in file %s.", name1, path);
1174  return;
1175  }
1176 
1177  fBackgroundScalingWeights = static_cast<TH2D*>(h1->Clone());
1178 
1179  // Open jet pT scale factor histogram
1180  TH2D* h2 = dynamic_cast<TH2D*>(file->Get(name2));
1181 
1182  if (h2) {
1183  ::Info("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Jet pT scaling histogram %s loaded from file %s.", name2, path);
1184  }
1185  else {
1186  ::Error("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Jet pT scaling histogram %s not found in file %s.", name2, path);
1187  return;
1188  }
1189 
1190  fGapJetScalingWeights = static_cast<TH2D*>(h2->Clone());
1191 
1192  file->Close();
1193  delete file;
1194 
1195 }
1196 
1202 {
1203 
1204  // Configure base class to set fTriggerPatchInfo to array of trigger patches, each event
1205  // (Need to call this before base class ExecOnce)
1206  if (fDoTriggerSimulation) {
1207  this->SetCaloTriggerPatchInfoName("EmcalTriggers");
1208  }
1209 
1211 
1212  fNeedEmcalGeom = kTRUE;
1213 
1214  // Load the PHOS geometry
1215  if (fDoCaloStudy) {
1216  fPHOSGeo = AliPHOSGeometry::GetInstance();
1217  if (fPHOSGeo) {
1218  AliInfo("Found instance of PHOS geometry!");
1219  }
1220  else {
1221  AliInfo("Creating PHOS geometry!");
1222  Int_t runNum = InputEvent()->GetRunNumber();
1223  if(runNum<209122) //Run1
1224  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
1225  else
1226  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
1227 
1228  if (fPHOSGeo) {
1229  AliOADBContainer geomContainer("phosGeo");
1230  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
1231  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
1232  for(Int_t mod=0; mod<6; mod++) {
1233  if(!matrixes->At(mod)) continue;
1234  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
1235  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
1236  ((TGeoHMatrix*)matrixes->At(mod))->Print();
1237  }
1238  }
1239  }
1240  }
1241 
1242  // Check if trigger patches are loaded
1243  if (fDoTriggerSimulation) {
1244  if (fTriggerPatchInfo) {
1245  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1246  TClass cls(objname);
1247  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1248  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1249  GetName(), cls.GetName(), "EmcalTriggers"));
1250  fTriggerPatchInfo = 0;
1251  }
1252  }
1253  if (!fTriggerPatchInfo) {
1254  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1255  return;
1256  }
1257  }
1258 
1259  AliInfo(Form("Trigger jet threshold = %f, Associated jet threshold = %f", fMinTrigJetPt, fMinAssJetPt));
1260  AliInfo(Form("Leading hadron threshold (for dijet leading jet): %f GeV", fDijetLeadingHadronPt));
1261  AliInfo(Form("Momentum balance study: %d", fDoMomentumBalance));
1262  AliInfo(Form("Geometrical matching study: %d", fDoGeometricalMatching));
1263 
1264 }
1265 
1270 
1271  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
1272 
1274 
1275  // Get instance of the downscale factor helper class
1277  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
1278 
1279  // There are two possible min bias triggers for LHC15o
1280  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
1281  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
1282  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
1283 
1284  // Get the downscale factor for whichever MB trigger exists in the given run
1285  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
1286  Double_t downscalefactor;
1287  for (auto i : runtriggers) {
1288  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
1289  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
1290  break;
1291  }
1292  }
1293 
1294  // Store the inverse of the downscale factor, used later to weight the pT spectrum
1295  fMBUpscaleFactor = 1/downscalefactor;
1296 
1297  TString histname = "Trigger/hMBDownscaleFactor";
1298  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
1299 
1300  }
1301 
1302 }
1303 
1308 {
1309  if (fUseAliEventCuts) {
1310  if (!fEventCuts.AcceptEvent(InputEvent()))
1311  {
1312  PostData(1, fOutput);
1313  return kFALSE;
1314  }
1315  }
1316  else {
1318  }
1319  return kTRUE;
1320 }
1321 
1330 {
1331  TString histname;
1332  AliJetContainer* jetCont = 0;
1333  TIter next(&fJetCollArray);
1334  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1335  TString jetContName = jetCont->GetName();
1336  if (jetContName.Contains("HardCore")) continue;
1337 
1338  //-----------------------------------------------------------------------------
1339  // Find the leading di-jet candidate in each event, and if it satisfies the
1340  // trig jet pT threshold, then fill di-jet candidate histogram (regardless of ass jet).
1341  // The idea is to study the kinematic selections in post-processing.
1342 
1343  // Loop over leading hadron cut or not
1344  for (Int_t leadingHadronCutType=0; leadingHadronCutType<2; leadingHadronCutType++) {
1345 
1346  // Find the dijet candidate of the event and store its info in struct fDijet
1347  FindDijet(jetCont, leadingHadronCutType);
1348 
1349  // If we find a dijet candidate (i.e. acceptable trig jet; ass jet accepted or not), fill the di-jet candidate histogram
1351  FillDijetCandHistograms(jetCont);
1352  }
1353 
1354  }
1355 
1356  //---------------------------------------------------------------------------------------------------
1357  // Now, study the accepted dijet selection -- specified by the trig/ass jet pT conditions
1358 
1359  // Find the dijet candidate of the event and store its info in struct fDijet
1360  FindDijet(jetCont, 0);
1361 
1362  // If we find an accepted dijet, fill the dijet imbalance histogram
1365  }
1366  // If we find an accepted dijet, perform momentum-balance study (if requested)
1368  histname = TString::Format("%s/MomentumBalance", jetCont->GetArrayName().Data());
1369  DoMomentumBalance(histname);
1370  }
1371 
1372  //---------------------------------------------------------------------------
1373  // Do a simple trigger simulation (if requested)
1374  if (fDoTriggerSimulation) {
1376  }
1377 
1378  }
1379 
1380  //---------------------------------------------------------------------------
1381  // Do the constituent threshold and geometrical matching study (if requested)
1382  if (fDoGeometricalMatching) {
1384  }
1385 
1386  // Only fill the embedding qa plots if:
1387  // - We are using the embedding helper
1388  // - The class has been initialized
1389  // - Both jet collections are available
1390  if (fEmbeddingQA.IsInitialized()) {
1392  }
1393 
1394  return kTRUE;
1395 }
1396 
1405 {
1406  fDijet.clear();
1407  fDijet.leadingHadronCutType = leadingHadronCutBin;
1408 
1409  // Get trigger jet
1410  AliEmcalJet* trigJet = GetLeadingJet(jetCont);
1411  if(!trigJet) return;
1412 
1413  // Skip the event if the leading jet doesn't satisfy the pT threshold
1414  Double_t trigJetPt = GetJetPt(jetCont, trigJet);
1415  if ( trigJetPt < fMinTrigJetPt ) return;
1416 
1417  // Skip the event if the leading jet doesn't satisfy the leading hadron threshold
1418  if (jetCont->GetLeadingHadronPt(trigJet) < fDijetLeadingHadronPt*leadingHadronCutBin) return;
1419 
1420  // Fill the dijet struct
1421  fDijet.trigJet = trigJet;
1422  fDijet.trigJetPt = trigJetPt;
1423  fDijet.trigJetEta = trigJet->Eta();
1424  fDijet.trigJetPhi = trigJet->Phi();
1425 
1426  // Find the subleading jet in the opposite hemisphere
1427  AliEmcalJet *assJet = 0;
1428  for(auto assJetCand : jetCont->accepted()) {
1429  if (!assJetCand) continue;
1430  Double_t assJetCandPt = GetJetPt(jetCont, assJetCand);
1431  if ( TMath::Abs(trigJet->Phi() - assJetCand->Phi()) < fDeltaPhiMin ) continue;
1432  if (assJet) {
1433  Double_t assJetPt = GetJetPt(jetCont, assJet);
1434  if ( assJetCandPt < assJetPt ) continue;
1435  }
1436  assJet = assJetCand;
1437  }
1438  if (!assJet) return;
1439 
1440  // Fill the dijet struct
1441  fDijet.assJet = assJet;
1442  fDijet.assJetPt = GetJetPt(jetCont, assJet);
1443  fDijet.assJetPhi = assJet->Phi();
1444  fDijet.assJetEta = assJet->Eta();
1446 
1447  fDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
1448  fDijet.deltaEta = trigJet->Eta() - assJet->Eta();
1451  fDijet.kTy = TMath::Abs( fDijet.trigJetPt * TMath::Sin(fDijet.deltaPhi) );
1452 }
1453 
1459 {
1460 
1461  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1462 
1463  AliVTrack* track;
1464  for (auto trackIterator : trackCont->accepted_momentum() ) {
1465 
1466  track = trackIterator.second;
1467 
1468  // Compute the delta phi between the track and its nearest jet (of the two jets in the dijet),
1469  // as well as its pT-parallel projection onto the nearest jet's axis.
1470 
1471  Double_t trackPt = track->Pt();
1472  Double_t trackPhi = track->Phi();
1473  Double_t trigJetPhi = fDijet.trigJet->Phi();
1474  Double_t assJetPhi = fDijet.assJet->Phi();
1475 
1476  Double_t deltaPhiTrigJet = TMath::Abs(trackPhi - trigJetPhi);
1477  Double_t deltaPhiAssJet = TMath::Abs(trackPhi - assJetPhi);
1478  Bool_t isNearside = deltaPhiTrigJet < deltaPhiAssJet;
1479 
1480  Double_t deltaPhi;
1481  Double_t balancePt;
1482  if (isNearside) {
1483  deltaPhi = trackPhi - trigJetPhi;
1484  balancePt = trackPt * TMath::Cos(deltaPhi);
1485  }
1486  else {
1487  deltaPhi = trackPhi - assJetPhi;
1488  balancePt = -trackPt * TMath::Cos(deltaPhi);
1489  }
1490 
1491  FillMomentumBalanceHistograms(histname, deltaPhi, trackPt, balancePt);
1492 
1493  }
1494 }
1495 
1500 {
1501  // Get jet container with minimum constituent pT,E thresholds
1502  TString jetContAllName = Form("Jet_AKTFullR0%d0_tracks_pT0150_caloClusters_E0300_pt_scheme", (int) (fMatchingJetR*10) );
1503  AliJetContainer* jetContAll = GetJetContainer(jetContAllName.Data());
1504 
1505  // Get jet container with X GeV constituent pT,E thresholds
1506  Int_t trackThreshold = (int) (fTrackConstituentThreshold*1000); // in MeV
1507  Int_t clusThreshold = (int) (fClusterConstituentThreshold*1000); // in MeV
1508  TString jetContHardCoreName = Form("JetHardCore_AKTFullR0%d0_tracks_pT%d_caloClusters_E%d_pt_scheme", (int) (fMatchingJetR*10), trackThreshold, clusThreshold);
1509  AliJetContainer* jetContHardCore = GetJetContainer(jetContHardCoreName.Data());
1510 
1511  // Find the di-jet in the hard-core jet sample, then find the matching di-jet and fill histograms
1512  FindDijet(jetContHardCore, 0);
1513  if (fDijet.isAccepted) {
1514  FindMatchingDijet(jetContAll);
1516  }
1517 }
1518 
1523 {
1524  TString histname;
1525 
1526  // Check if trigger patches are loaded
1527  if (fTriggerPatchInfo) {
1528  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1529  TClass cls(objname);
1530  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1531  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1532  GetName(), cls.GetName(), "EmcalTriggers"));
1533  fTriggerPatchInfo = 0;
1534  }
1535  }
1536  if (!fTriggerPatchInfo) {
1537  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1538  return;
1539  }
1540 
1541  // Compute patches in EMCal, DCal (I want offline simple trigger patch, i.e. patch calculated using FEE energy)
1542  std::vector<Double_t> vecEMCal;
1543  std::vector<Double_t> vecDCal;
1544  for(auto p : *fTriggerPatchInfo){
1545  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1546  if (recpatch) {
1547 
1548  if(!recpatch->IsJetHighSimple()) continue;
1549 
1550  histname = "TriggerSimHistograms/hEtaVsPhi";
1551  fHistManager.FillTH2(histname.Data(), recpatch->GetEtaGeo(), recpatch->GetPhiGeo());
1552 
1553  histname = "TriggerSimHistograms/hPatchE";
1554  fHistManager.FillTH2(histname.Data(), fCent, recpatch->GetPatchE());
1555 
1556  if (recpatch->IsEMCal()) {
1557  vecEMCal.push_back(recpatch->GetPatchE());
1558  } else {
1559  vecDCal.push_back(recpatch->GetPatchE());
1560  }
1561 
1562  }
1563  }
1564 
1565  // Compute the median in each calorimeter
1566  const Int_t nBkgPatchesEMCal = vecEMCal.size(); // 6*8;
1567  const Int_t nBkgPatchesDCal = vecDCal.size(); // 4*5;
1568  fMedianEMCal = TMath::Median(nBkgPatchesEMCal, &vecEMCal[0]); // point to array used internally by vector
1569  fMedianDCal = TMath::Median(nBkgPatchesDCal, &vecDCal[0]);
1570 
1571  histname = "TriggerSimHistograms/hPatchMedianE";
1572  fHistManager.FillTH3(histname.Data(), fCent, fMedianEMCal, kEMCal);
1573  fHistManager.FillTH3(histname.Data(), fCent, fMedianDCal, kDCal);
1574 
1575  histname = "TriggerSimHistograms/hNPatches";
1576  fHistManager.FillTH2(histname.Data(), nBkgPatchesEMCal, kEMCal);
1577  fHistManager.FillTH2(histname.Data(), nBkgPatchesDCal, kDCal);
1578 
1579  // Median patch energy vs. pT for dijets
1580  if (fDijet.isAccepted) {
1581  histname = "TriggerSimHistograms/hMedPatchDijet";
1583  fHistManager.FillTHnSparse(histname, x);
1585  fHistManager.FillTHnSparse(histname, y);
1586  }
1587 
1588  // Then compute background subtracted patches, by subtracting from each patch the median patch E from the opposite hemisphere
1589  // If a patch is above threshold, the event is "triggered"
1590  Bool_t fkEMCEJE = kFALSE;
1591  Double_t threshold = 20;
1592  for(auto p : *fTriggerPatchInfo){
1593  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1594  if (recpatch) {
1595 
1596  if(!recpatch->IsJetHighSimple()) continue;
1597 
1598  if (recpatch->IsEMCal()) {
1599  if ((recpatch->GetPatchE() - fMedianDCal) > threshold) {
1600  fkEMCEJE = kTRUE;
1601  break;
1602  }
1603  } else {
1604  if ((recpatch->GetPatchE() - fMedianEMCal) > threshold) {
1605  fkEMCEJE = kTRUE;
1606  break;
1607  }
1608  }
1609  }
1610  }
1611 
1612  if (fkEMCEJE) {
1614  }
1615 
1616 }
1617 
1623 {
1625 
1626  // Loop over jets and find leading jet within R of fDijet.trigJet
1627  AliEmcalJet *matchingTrigJet = 0;
1628  for(auto matchingTrigJetCand : jetCont->accepted()) {
1629  if (!matchingTrigJetCand) continue;
1630  if ( GetDeltaR(matchingTrigJetCand, fDijet.trigJet) > fMatchingJetR ) continue;
1631  if (matchingTrigJet) {
1632  if ( GetJetPt(jetCont, matchingTrigJetCand) < GetJetPt(jetCont, matchingTrigJet) ) continue;
1633  }
1634  matchingTrigJet = matchingTrigJetCand;
1635  }
1636  if (!matchingTrigJet) return;
1637 
1638  // Loop over jets and find leading jet within R of fDijet.assJet
1639  AliEmcalJet *matchingAssJet = 0;
1640  for(auto matchingAssJetCand : jetCont->accepted()) {
1641  if (!matchingAssJetCand) continue;
1642  if ( GetDeltaR(matchingAssJetCand, fDijet.assJet) > fMatchingJetR ) continue;
1643  if (matchingAssJet) {
1644  if ( GetJetPt(jetCont, matchingAssJetCand) < GetJetPt(jetCont, matchingAssJet) ) continue;
1645  }
1646  matchingAssJet = matchingAssJetCand;
1647  }
1648 
1649  // Determine which matching jet is the leading jet (i.e. allow them to flip)
1650  if (matchingAssJet) {
1651  AliEmcalJet* trigJet = matchingTrigJet;
1652  AliEmcalJet* assJet = matchingAssJet;
1653  if ( GetJetPt(jetCont, matchingTrigJet) < GetJetPt(jetCont, matchingAssJet) ) {
1654  AliEmcalJet* trigJet = matchingAssJet;
1655  AliEmcalJet* assJet = matchingTrigJet;
1656  }
1657 
1658  // Fill the dijet struct
1659  fMatchingDijet.trigJet = trigJet;
1660  fMatchingDijet.trigJetPt = GetJetPt(jetCont, trigJet);
1661  fMatchingDijet.trigJetEta = trigJet->Eta();
1662  fMatchingDijet.trigJetPhi = trigJet->Phi();
1663 
1664  fMatchingDijet.assJet = assJet;
1665  fMatchingDijet.assJetPt = GetJetPt(jetCont, assJet);
1666  fMatchingDijet.assJetPhi = assJet->Phi();
1667  fMatchingDijet.assJetEta = assJet->Eta();
1669 
1670  fMatchingDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
1671  fMatchingDijet.deltaEta = trigJet->Eta() - assJet->Eta();
1674  fMatchingDijet.kTy = TMath::Abs( fMatchingDijet.trigJetPt * TMath::Sin(fMatchingDijet.deltaPhi) );
1675  }
1676 }
1677 
1682 {
1683  // Loop over tracks and clusters in order to:
1684  // (1) Compute scale factor for full jets
1685  // (2) Compute delta-pT for full jets, with the random cone method
1686  // For both the scale factor and delta-pT, we compute only one histogram each for EMCal.
1687  // And then we bin in eta-phi, in order to compute and perform a corretion to account for the DCal vs. PHOS vs. gap
1688 
1689  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
1690  Double_t etaTPC = 0.9;
1691  Double_t etaEMCal = 0.7;
1692  Double_t etaMinDCal = 0.22;
1693  Double_t etaMaxPHOS = 0.13;
1694  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
1695  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
1696  Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
1697  Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
1698  Double_t phiMinPHOS = 250 * TMath::DegToRad();
1699  Double_t phiMaxPHOS = 320 * TMath::DegToRad();
1700 
1701  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
1702  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
1703  Double_t accDCalRegion = 2 * etaEMCal * (phiMaxDCal - phiMinDCal);
1704 
1705  // Define fiducial acceptances, to be used to generate random cones
1706  TRandom3* r = new TRandom3(0);
1707  Double_t jetR = jetCont->GetJetRadius();
1708  Double_t accRC = TMath::Pi() * jetR * jetR;
1709  Double_t etaEMCalfid = etaEMCal - jetR;
1710  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
1711  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
1712  Double_t phiMinDCalRegionfid = phiMinDCal + jetR;
1713  Double_t phiMaxDCalRegionfid = phiMaxDCal - jetR;
1714 
1715  // Generate EMCal random cone eta-phi
1716  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
1717  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
1718 
1719  // For eta-phi correction, generate random eta, phi in each eta/phi bin, to be used as center of random cone
1720  Double_t etaDCalRC[fNEtaBins]; // array storing the RC eta values
1721  Double_t etaStep = 1./fNEtaBins;
1722  Double_t etaMin;
1723  Double_t etaMax;
1724  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1725  etaMin = -etaEMCalfid + etaBin*etaStep;
1726  etaMax = etaMin + etaStep;
1727  etaDCalRC[etaBin] = r->Uniform(etaMin, etaMax);
1728  }
1729 
1730  Double_t phiDCalRC[fNPhiBins]; // array storing the RC phi values
1731  Double_t phiStep = 5./fNPhiBins; // phi axis is [1,6] in order to have simple binning
1732  Double_t phiMin;
1733  Double_t phiMax;
1734  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1735  phiMin = 1 + phiBin*phiStep;
1736  phiMax = phiMin + phiStep;
1737  phiDCalRC[phiBin] = r->Uniform(phiMin, phiMax);
1738  }
1739 
1740  // Initialize the various sums to 0
1741  Double_t trackPtSumTPC = 0;
1742  Double_t trackPtSumEMCal = 0;
1743  Double_t trackPtSumEMCalRC = 0;
1744  Double_t clusESumEMCal = 0;
1745  Double_t clusESumEMCalRC = 0;
1746 
1747  // Define a 2D vector (initialized to 0) to store the sum of track pT, and another for cluster ET
1748  std::vector<std::vector<Double_t>> trackPtSumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1749  std::vector<std::vector<Double_t>> clusESumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1750 
1751  // Loop over tracks. Sum the track pT:
1752  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal random cone,
1753  // (4) in a random cone at each eta-phi
1754  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1755  AliTLorentzVector track;
1756  Double_t trackEta;
1757  Double_t trackPhi;
1758  Double_t trackPt;
1759  Double_t deltaR;
1760  for (auto trackIterator : trackCont->accepted_momentum() ) {
1761 
1762  track.Clear();
1763  track = trackIterator.first;
1764  trackEta = track.Eta();
1765  trackPhi = track.Phi_0_2pi();
1766  trackPt = track.Pt();
1767 
1768  // (1)
1769  if (TMath::Abs(trackEta) < etaTPC) {
1770  trackPtSumTPC += trackPt;
1771  }
1772 
1773  // (2)
1774  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
1775  trackPtSumEMCal += trackPt;
1776  }
1777 
1778  // (3)
1779  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
1780  if (deltaR < jetR) {
1781  trackPtSumEMCalRC += trackPt;
1782  }
1783 
1784  // (4)
1785  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1786  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1787  deltaR = GetDeltaR(&track, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1788  if (deltaR < jetR) {
1789  trackPtSumDCalRC[etaBin][phiBin] += trackPt;
1790  }
1791  }
1792  }
1793 
1794  }
1795 
1796  // Loop over clusters. Sum the cluster ET:
1797  // (1) in the EMCal, (2) in the EMCal random cone, (3) in a random cone at each eta-phi
1799  AliTLorentzVector clus;
1800  Double_t clusEta;
1801  Double_t clusPhi;
1802  Double_t clusE;
1803  for (auto clusIterator : clusCont->accepted_momentum() ) {
1804 
1805  clus.Clear();
1806  clus = clusIterator.first;
1807  clusEta = clus.Eta();
1808  clusPhi = clus.Phi_0_2pi();
1809  clusE = clus.E();
1810 
1811  // (1)
1812  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
1813  clusESumEMCal += clusE;
1814  }
1815 
1816  // (2)
1817  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
1818  if (deltaR < jetR) {
1819  clusESumEMCalRC += clusE;
1820  }
1821 
1822  // (3)
1823  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1824  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1825  deltaR = GetDeltaR(&clus, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1826  if (deltaR < jetR) {
1827  clusESumDCalRC[etaBin][phiBin] += clusE;
1828  }
1829  }
1830  }
1831 
1832  }
1833 
1834  // Compute the scale factor for EMCal, as a function of centrality
1835  Double_t numerator = (trackPtSumEMCal + clusESumEMCal) / accEMCal;
1836  Double_t denominator = trackPtSumTPC / accTPC;
1837  Double_t scaleFactor = numerator / denominator;
1838  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
1839  fHistManager.FillTH2(histname, fCent, scaleFactor);
1840 
1841  // Compute the scale factor in each eta-phi bin, as a function of centrality
1842  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1843  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1844  numerator = (trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin]) / accRC;
1845  scaleFactor = numerator / denominator;
1846  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jetCont->GetArrayName().Data());
1847  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, scaleFactor};
1848  fHistManager.FillTHnSparse(histname, x);
1849  }
1850  }
1851 
1852  // Compute delta pT for EMCal, as a function of centrality
1853  Double_t rho = jetCont->GetRhoVal();
1854  Double_t deltaPt = trackPtSumEMCalRC + clusESumEMCalRC - rho * TMath::Pi() * jetR * jetR;
1855  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
1856  fHistManager.FillTH2(histname, fCent, deltaPt);
1857 
1858  // Compute delta pT in each eta-phi bin, as a function of centrality
1859  Double_t sf;
1860  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1861  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1863  sf = fBackgroundScalingWeights->GetBinContent(fBackgroundScalingWeights->FindBin(etaDCalRC[etaBin], phiDCalRC[phiBin]));
1864  rho = sf * jetCont->GetRhoVal();
1865  }
1866  deltaPt = trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin] - rho * accRC;
1867  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jetCont->GetArrayName().Data());
1868  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, deltaPt};
1869  fHistManager.FillTHnSparse(histname, x);
1870  }
1871  }
1872 
1873  delete r;
1874 
1875 }
1876 
1881 {
1882 
1883  // Get eta-phi dependent jet pT scale factor
1884  Double_t jetPt = jet->Pt();
1885  if (fGapJetScalingWeights) {
1886  Double_t sf = fGapJetScalingWeights->GetBinContent(fGapJetScalingWeights->FindBin(jet->Eta(), jet->Phi_0_2pi()));
1887  jetPt = jetPt * (1 + sf * jet->NEF());
1888  }
1889 
1890  // Compute pTcorr
1891  Double_t rho = jetCont->GetRhoVal();
1892  Double_t pT = jetPt - rho * jet->Area();
1893 
1894  // If hard-core jet, don't subtract background
1895  TString jetContName = jetCont->GetName();
1896  if (jetContName.Contains("HardCore")) pT = jet->Pt();
1897 
1898  return pT;
1899 }
1900 
1905 {
1906  AliEmcalJet* leadingJet = 0;
1907 
1908  if (jetCont->GetRhoParameter()) {
1909  for(auto jetCand : jetCont->accepted()) {
1910  if (!jetCand) continue;
1911  Double_t jetCandPt = GetJetPt(jetCont, jetCand);
1912  if (leadingJet) {
1913  Double_t leadingJetPt = GetJetPt(jetCont, leadingJet);
1914  if ( jetCandPt < leadingJetPt ) continue;
1915  }
1916  leadingJet = jetCand;
1917  }
1918  }
1919  else {
1920  leadingJet = jetCont->GetLeadingJet();
1921  }
1922 
1923  return leadingJet;
1924 }
1925 
1930 {
1931  Double_t deltaPhi = TMath::Abs(jet1->Phi() - jet2->Phi());
1932  Double_t deltaEta = TMath::Abs(jet1->Eta() - jet2->Eta());
1933  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1934  return deltaR;
1935 }
1936 
1941 {
1942  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
1943  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
1944  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1945  return deltaR;
1946 }
1947 
1952 {
1953  UInt_t jetType = jet->GetJetAcceptanceType();
1954  Double_t type = -1;
1955  if (jetType & AliEmcalJet::kEMCAL) {
1956  type = kEMCal;
1957  }
1958  else if (jetType & AliEmcalJet::kDCALonly) {
1959  type = kDCal;
1960  }
1961  else if (jetType & AliEmcalJet::kPHOS) {
1962  type = kPHOS;
1963  }
1964 
1965  return type;
1966 }
1967 
1968 
1976 {
1979 
1980  return kTRUE;
1981 }
1982 
1988 {
1989  TString histname;
1990  AliJetContainer* jets = 0;
1991  TIter nextJetColl(&fJetCollArray);
1992  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1993  TString jetContName = jets->GetName();
1994  if (jetContName.Contains("HardCore")) continue;
1995  Double_t rhoVal = 0;
1996  if (jets->GetRhoParameter()) {
1997  rhoVal = jets->GetRhoVal();
1998  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
1999  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
2000  }
2001 
2002  for (auto jet : jets->all()) {
2003 
2004  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
2005  Float_t corrPt = GetJetPt(jets, jet);
2006 
2007  // A vs. pT (fill before area cut)
2008  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
2009  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
2010 
2011 
2012  // Rejection reason
2013  UInt_t rejectionReason = 0;
2014  if (!jets->AcceptJet(jet, rejectionReason)) {
2015  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
2016  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
2017  continue;
2018  }
2019 
2020  // compute jet acceptance type
2021  Double_t type = GetJetType(jet);
2022 
2023  // (Centrality, pT, calo type)
2024  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
2025  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type);
2026 
2027  // (Centrality, eta, phi, pT, NEF, calo type)
2028  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
2029  Double_t x[6] = {fCent, jet->Eta(), jet->Phi_0_2pi(), corrPt, jet->NEF(), type};
2030  fHistManager.FillTHnSparse(histname, x);
2031 
2032  // (Centrality, pT upscaled, calo type)
2033  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
2034  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type, fMBUpscaleFactor);
2035 
2036  // pT-leading vs. pT
2037  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
2038  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
2039 
2040  // (Centrality, pT, z-leading (charged), calo type)
2041  TLorentzVector leadPart;
2042  jets->GetLeadingHadronMomentum(leadPart, jet);
2043  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
2044  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
2045  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
2046  Double_t y[4] = {fCent, corrPt, z, type};
2047  fHistManager.FillTHnSparse(histname, y);
2048 
2049  // (Centrality, pT, z (charged), calo type)
2050  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
2051  AliVTrack* track;
2052  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
2053  track = static_cast<AliVTrack*>(jet->Track(i));
2054  z = track->Pt() / TMath::Abs(corrPt);
2055  Double_t y2[4] = {fCent, corrPt, z, type};
2056  fHistManager.FillTHnSparse(histname, y2);
2057  }
2058 
2059  // (Centrality, pT, Nconst, calo type)
2060  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
2061  Double_t a[4] = {fCent, corrPt, 1.*jet->GetNumberOfConstituents(), type};
2062  fHistManager.FillTHnSparse(histname, a);
2063 
2064  // (Median patch energy, calo type, jet pT, centrality)
2065  if (fDoTriggerSimulation) {
2066  histname = TString::Format("%s/JetHistograms/hMedPatchJet", jets->GetArrayName().Data());
2067  Double_t x[4] = {fMedianEMCal, kEMCal, corrPt, fCent};
2068  fHistManager.FillTHnSparse(histname, x);
2069  Double_t y[4] = {fMedianDCal, kDCal, corrPt, fCent};
2070  fHistManager.FillTHnSparse(histname, y);
2071  }
2072 
2073  } //jet loop
2074 
2075  //---------------------------------------------------------------------------
2076  // Do study of background (if requested)
2078  }
2079 }
2080 
2086 {
2087  TString histname;
2088  AliJetContainer* jets = 0;
2089  TIter nextJetColl(&fJetCollArray);
2090  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
2091  TString jetContName = jets->GetName();
2092  if (jetContName.Contains("HardCore")) continue;
2093  Double_t rhoVal = 0;
2094  if (jets->GetRhoParameter()) {
2095  rhoVal = jets->GetRhoVal();
2096  histname = TString::Format("%s/TriggerSimHistograms/hRhoVsCent", jets->GetArrayName().Data());
2097  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
2098  }
2099 
2100  for (auto jet : jets->all()) {
2101 
2102  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
2103  Float_t corrPt = GetJetPt(jets, jet);
2104 
2105  // A vs. pT (fill before area cut)
2106  histname = TString::Format("%s/TriggerSimHistograms/hAreaVsPt", jets->GetArrayName().Data());
2107  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
2108 
2109 
2110  // Rejection reason
2111  UInt_t rejectionReason = 0;
2112  if (!jets->AcceptJet(jet, rejectionReason)) {
2113  histname = TString::Format("%s/TriggerSimHistograms/hJetRejectionReason", jets->GetArrayName().Data());
2114  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
2115  continue;
2116  }
2117 
2118  // compute jet acceptance type
2119  Double_t type = GetJetType(jet);
2120 
2121  // Centrality vs. pT
2122  histname = TString::Format("%s/TriggerSimHistograms/hCentVsPt", jets->GetArrayName().Data());
2123  fHistManager.FillTH3(histname.Data(), fCent, corrPt, type);
2124 
2125  // (Centrality, eta, phi, pT, NEF, calo type)
2126  histname = TString::Format("%s/TriggerSimHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
2127  Double_t x[6] = {fCent, jet->Eta(), jet->Phi_0_2pi(), corrPt, jet->NEF(), type};
2128  fHistManager.FillTHnSparse(histname, x);
2129 
2130  // pT-leading vs. pT
2131  histname = TString::Format("%s/TriggerSimHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
2132  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
2133 
2134  // (Centrality, pT, z-leading (charged), calo type)
2135  TLorentzVector leadPart;
2136  jets->GetLeadingHadronMomentum(leadPart, jet);
2137  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
2138  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
2139  histname = TString::Format("%s/TriggerSimHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
2140  Double_t y[4] = {fCent, corrPt, z, type};
2141  fHistManager.FillTHnSparse(histname, y);
2142 
2143  // (Centrality, pT, z (charged), calo type)
2144  histname = TString::Format("%s/TriggerSimHistograms/hZVsPt", jets->GetArrayName().Data());
2145  AliVTrack* track;
2146  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
2147  track = static_cast<AliVTrack*>(jet->Track(i));
2148  z = track->Pt() / TMath::Abs(corrPt);
2149  Double_t y2[4] = {fCent, corrPt, z, type};
2150  fHistManager.FillTHnSparse(histname, y2);
2151  }
2152 
2153  // (Centrality, pT, Nconst, calo type)
2154  histname = TString::Format("%s/TriggerSimHistograms/hNConstVsPt", jets->GetArrayName().Data());
2155  Double_t a[4] = {fCent, corrPt, 1.*jet->GetNumberOfConstituents(), type};
2156  fHistManager.FillTHnSparse(histname, a);
2157 
2158  } //jet loop
2159 
2160  }
2161 }
2162 
2167 {
2168  Double_t contents[30]={0};
2169  TString histname = TString::Format("%s/DijetCandObservables", jets->GetArrayName().Data());
2170  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
2171  if (!histJetObservables) return;
2172  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
2173  TString title(histJetObservables->GetAxis(n)->GetTitle());
2174  if (title=="Centrality (%)")
2175  contents[n] = fCent;
2176  else if (title=="LeadingHadronRequired")
2177  contents[n] = fDijet.leadingHadronCutType;
2178  else if (title=="#it{p}_{T,trig jet} (GeV/#it{c})")
2179  contents[n] = fDijet.trigJetPt;
2180  else if (title=="#it{p}_{T,ass jet} (GeV/#it{c})")
2181  contents[n] = fDijet.assJetPt;
2182  else if (title=="#phi_{trig jet}")
2183  contents[n] = fDijet.trigJetPhi;
2184  else if (title=="#phi_{ass jet}")
2185  contents[n] = fDijet.assJetPhi;
2186  else if (title=="#eta_{trig jet}")
2187  contents[n] = fDijet.trigJetEta;
2188  else if (title=="#eta_{ass jet}")
2189  contents[n] = fDijet.assJetEta;
2190  else
2191  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
2192  }
2193  histJetObservables->Fill(contents);
2194 }
2195 
2200 {
2201  // Fill the dijet imbalance histogram (unless doing trigger simulation, in which case bypass this)
2202  if (!fDoTriggerSimulation) {
2203  Double_t contents[30]={0};
2204  TString histname = TString::Format("%s/DijetImbalanceObservables", jets->GetArrayName().Data());
2205  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
2206  if (!histJetObservables) return;
2207  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
2208  TString title(histJetObservables->GetAxis(n)->GetTitle());
2209  if (title=="Centrality (%)")
2210  contents[n] = fCent;
2211  else if (title=="#Delta#phi")
2212  contents[n] = fDijet.deltaPhi;
2213  else if (title=="#Delta#eta")
2214  contents[n] = fDijet.deltaEta;
2215  else if (title=="A_{J}")
2216  contents[n] = fDijet.AJ;
2217  else if (title=="x_{J}")
2218  contents[n] = fDijet.xJ;
2219  else if (title=="k_{Ty} (GeV)")
2220  contents[n] = fDijet.kTy;
2221  else if (title=="N_{tracks, trig jet}")
2222  contents[n] = fDijet.trigJet->GetNumberOfTracks();
2223  else if (title=="N_{tracks, ass jet}")
2224  contents[n] = fDijet.assJet->GetNumberOfTracks();
2225  else
2226  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
2227  }
2228  histJetObservables->Fill(contents);
2229  }
2230 
2231  // (Centrality, pT1, pT2) (upscaled)
2232  TString histname = TString::Format("%s/DijetJetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
2234 
2235  // Get jet acceptance type
2236  AliEmcalJet* trigJet = fDijet.trigJet;
2237  Double_t type = GetJetType(trigJet);
2238 
2239  // (Centrality, pT, NEF, calo type)
2240  histname = TString::Format("%s/DijetJetHistograms/hNEFVsPt", jets->GetArrayName().Data());
2241  Double_t x[4] = {fCent, fDijet.trigJetPt, trigJet->NEF(), type};
2242  fHistManager.FillTHnSparse(histname, x);
2243 
2244  // (Centrality, pT, z-leading (charged), calo type)
2245  histname = TString::Format("%s/DijetJetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
2246  TLorentzVector leadPart;
2247  jets->GetLeadingHadronMomentum(leadPart, trigJet);
2248  Double_t z = GetParallelFraction(leadPart.Vect(), trigJet);
2249  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
2250  Double_t y[4] = {fCent, fDijet.trigJetPt, z, type};
2251  fHistManager.FillTHnSparse(histname, y);
2252 
2253  // (Centrality, pT, z (charged), calo type)
2254  histname = TString::Format("%s/DijetJetHistograms/hZVsPt", jets->GetArrayName().Data());
2255  AliVTrack* track;
2256  for (Int_t i=0; i<trigJet->GetNumberOfTracks(); i++) {
2257  track = static_cast<AliVTrack*>(trigJet->Track(i));
2258  z = track->Pt() / TMath::Abs(fDijet.trigJetPt);
2259  Double_t y2[4] = {fCent, fDijet.trigJetPt, z, type};
2260  fHistManager.FillTHnSparse(histname, y2);
2261  }
2262 
2263  // (Centrality, pT, Nconst, calo type)
2264  histname = TString::Format("%s/DijetJetHistograms/hNConstVsPt", jets->GetArrayName().Data());
2265  Double_t a[4] = {fCent, fDijet.trigJetPt, 1.*trigJet->GetNumberOfConstituents(), type};
2266  fHistManager.FillTHnSparse(histname, a);
2267 
2268 }
2269 
2274 {
2275  Double_t contents[30]={0};
2276  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
2277  if (!histJetObservables) return;
2278  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
2279  TString title(histJetObservables->GetAxis(n)->GetTitle());
2280  if (title=="A_{J}")
2281  contents[n] = fDijet.AJ;
2282  else if (title=="#Delta#phi")
2283  contents[n] = deltaPhi;
2284  else if (title=="#it{p}_{T,particle} (GeV/#it{c})")
2285  contents[n] = trackPt;
2286  else if (title=="#it{p}_{T}#parallel (GeV/#it{c})")
2287  contents[n] = balancePt;
2288  else
2289  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
2290  }
2291  histJetObservables->Fill(contents);
2292 }
2293 
2298 {
2299  // Matching efficiency histogram
2300  TString histname = "GeometricalMatchingEfficiency";
2301 
2302  // If we have a matching di-jet, fill the geometrical matching histogram
2303  if (fMatchingDijet.assJet) {
2304  fHistManager.FillTH1(histname.Data(), 1);
2305 
2308  Bool_t isSwitched = trigDeltaR > fMatchingJetR;
2309 
2310  TString thnname = "GeometricalMatching";
2311  Double_t contents[30]={0};
2312  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(thnname.Data()));
2313  if (!histJetObservables) return;
2314  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
2315  TString title(histJetObservables->GetAxis(n)->GetTitle());
2316  if (title=="Centrality (%)")
2317  contents[n] = fCent;
2318  else if (title=="isSwitched")
2319  contents[n] = isSwitched;
2320  else if (title=="#DeltaR_{trig}")
2321  contents[n] = trigDeltaR;
2322  else if (title=="#DeltaR_{ass}")
2323  contents[n] = assDeltaR;
2324  else if (title=="trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
2325  contents[n] = fMatchingDijet.trigJetPt - fDijet.trigJetPt;
2326  else if (title=="ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
2327  contents[n] = fMatchingDijet.assJetPt - fDijet.assJetPt;
2328  else if (title=="A_{J} low-threshold")
2329  contents[n] = fMatchingDijet.AJ;
2330  else if (title=="A_{J} hard-core")
2331  contents[n] = fDijet.AJ;
2332  else
2333  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
2334  }
2335  histJetObservables->Fill(contents);
2336  }
2337  else {
2338  fHistManager.FillTH1(histname.Data(), 0.);
2339  }
2340 }
2341 
2342 /*
2343  * This function fills the histograms for the calorimeter performance study.
2344  */
2346 {
2347  // Define some vars
2348  TString histname;
2349  Double_t Enonlin;
2350  Double_t Ehadcorr;
2351  Int_t absId;
2352  Double_t ecell;
2353  Double_t leadEcell;
2354 
2355  // Get cells from event
2356  fCaloCells = InputEvent()->GetEMCALCells();
2357  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
2358 
2359  // Loop through clusters and plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
2360  AliClusterContainer* clusters = 0;
2361  TIter nextClusColl(&fClusterCollArray);
2362  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
2364  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
2365 
2366  // Determine cluster type (EMCal/DCal/Phos)
2367  ClusterType clusType = kNA;
2368  if (it->second->IsEMCAL()) {
2369  Double_t phi = it->first.Phi_0_2pi();
2370  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
2371  if (isDcal == 0) {
2372  clusType = kEMCal;
2373  } else if (isDcal == 1) {
2374  clusType = kDCal;
2375  }
2376  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
2377  clusType = kPHOS;
2378  }
2379 
2380  // rejection reason plots, to make efficiency correction
2381  if (it->second->IsEMCAL()) {
2382  histname = TString::Format("%s/hClusterRejectionReasonEMCal", clusters->GetArrayName().Data());
2383  UInt_t rejectionReason = 0;
2384  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
2385  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
2386  continue;
2387  }
2388  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
2389  histname = TString::Format("%s/hClusterRejectionReasonPHOS", clusters->GetArrayName().Data());
2390  UInt_t rejectionReason = 0;
2391  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
2392  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
2393  continue;
2394  }
2395  } else {
2396  continue;
2397  }
2398 
2399  // Fill cluster spectra by SM, and fill cell histograms
2400  Enonlin = 0;
2401  Ehadcorr = 0;
2402  if (it->second->IsEMCAL()) {
2403 
2404  Ehadcorr = it->second->GetHadCorrEnergy();
2405  Enonlin = it->second->GetNonLinCorrEnergy();
2407  Enonlin = it->second->E();
2408  }
2409 
2410  if (fPlotExotics) {
2411  histname = TString::Format("%s/hFcrossEMCal", clusters->GetArrayName().Data());
2412  Double_t Fcross = GetFcross(it->second, fCaloCells);
2413  fHistManager.FillTH3(histname, fCent, Fcross, it->second->E());
2414  }
2415 
2416  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
2417  if (sm >=0 && sm < 20) {
2418  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
2419  fHistManager.FillTH1(histname, it->second->E());
2420  }
2421  else {
2422  AliError(Form("Supermodule %d does not exist!", sm));
2423  }
2424 
2425  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
2426  histname = TString::Format("Cells/hCellEnergyAccepted");
2427  leadEcell = 0;
2428  for (Int_t iCell = 0; iCell < it->second->GetNCells(); iCell++){
2429  absId = it->second->GetCellAbsId(iCell);
2430  ecell = fCaloCells->GetCellAmplitude(absId);
2431  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
2432  if (ecell > leadEcell) {
2433  leadEcell = ecell;
2434  }
2435  }
2436  // Plot also the leading cell
2437  histname = TString::Format("Cells/hCellEnergyLeading");
2438  fHistManager.FillTH3(histname, leadEcell, fCent, kEMCal);
2439 
2440  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
2441 
2442  Ehadcorr = it->second->GetCoreEnergy();
2443  Enonlin = it->second->E();
2444 
2445  Int_t relid[4];
2446  if (fPHOSGeo) {
2447  fPHOSGeo->AbsToRelNumbering(it->second->GetCellAbsId(0), relid);
2448  Int_t sm = relid[0];
2449  if (sm >=1 && sm < 5) {
2450  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
2451  fHistManager.FillTH1(histname, it->second->E());
2452  }
2453  else {
2454  AliError(Form("Supermodule %d does not exist!", sm));
2455  }
2456  }
2457 
2458  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
2459  histname = TString::Format("Cells/hCellEnergyAccepted");
2460  leadEcell = 0;
2461  for (Int_t iCell = 0; iCell < it->second->GetNCells(); iCell++){
2462  absId = it->second->GetCellAbsId(iCell);
2463  ecell = phosCaloCells->GetCellAmplitude(absId);
2464  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
2465  if (ecell > leadEcell) {
2466  leadEcell = ecell;
2467  }
2468  }
2469  // Plot also the leading cell
2470  histname = TString::Format("Cells/hCellEnergyLeading");
2471  fHistManager.FillTH3(histname, leadEcell, fCent, kPHOS);
2472  }
2473 
2474  // Check if the cluster has a matched track
2475  Int_t hasMatchedTrack = -1;
2476  Int_t nMatchedTracks = it->second->GetNTracksMatched();
2477  if (nMatchedTracks == 0) {
2478  hasMatchedTrack = 0;
2479  } else if (nMatchedTracks > 0) {
2480  hasMatchedTrack = 1;
2481  }
2482 
2483  // Check if the cluster passes the dispersion cut for photon-like cluster (meaningful only for PHOS)
2484  Int_t passedDispersionCut = 0;
2485  if (it->second->Chi2() < 2.5*2.5) {
2486  passedDispersionCut = 1;
2487  }
2488 
2489  if (fPlotClusterTHnSparse) {
2490  Double_t contents[30]={0};
2491  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
2492  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
2493  if (!histClusterObservables) return;
2494  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
2495  TString title(histClusterObservables->GetAxis(i)->GetTitle());
2496  if (title=="Centrality %")
2497  contents[i] = fCent;
2498  else if (title=="#eta")
2499  contents[i] = it->first.Eta();
2500  else if (title=="#phi")
2501  contents[i] = it->first.Phi_0_2pi();
2502  else if (title=="#it{E}_{clus} (GeV)")
2503  contents[i] = Enonlin;
2504  else if (title=="#it{E}_{clus, hadcorr} or #it{E}_{core} (GeV)")
2505  contents[i] = Ehadcorr;
2506  else if (title=="Matched track")
2507  contents[i] = hasMatchedTrack;
2508  else if (title=="M02")
2509  contents[i] = it->second->GetM02();
2510  else if (title=="Ncells")
2511  contents[i] = it->second->GetNCells();
2512  else if (title=="Dispersion cut")
2513  contents[i] = passedDispersionCut;
2514  else
2515  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
2516  }
2517  histClusterObservables->Fill(contents);
2518  }
2519 
2520  }
2521 
2522  }
2523 
2524  AliJetContainer* jets = 0;
2525  TIter nextJetColl(&fJetCollArray);
2526  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
2527 
2528  for (auto jet : jets->accepted()) {
2529 
2530  // plot neutral jets THnSparse (centrality, eta, phi, E, Nclusters)
2531  if (fPlotNeutralJets) {
2532  Double_t contents[30]={0};
2533  histname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
2534  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
2535  if (!histJetObservables) return;
2536  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
2537  TString title(histJetObservables->GetAxis(i)->GetTitle());
2538  if (title=="Centrality (%)")
2539  contents[i] = fCent;
2540  else if (title=="#eta_{jet}")
2541  contents[i] = jet->Eta();
2542  else if (title=="#phi_{jet} (rad)")
2543  contents[i] = jet->Phi_0_2pi();
2544  else if (title=="#it{E}_{T} (GeV)")
2545  contents[i] = jet->Pt();
2546  else if (title=="#rho (GeV/#it{c})")
2547  contents[i] = jet->Pt() / jet->Area();
2548  else if (title=="N_{clusters}")
2549  contents[i] = jet->GetNumberOfClusters();
2550  else
2551  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
2552  }
2553  histJetObservables->Fill(contents);
2554  }
2555 
2556  // Fill cluster spectra of clusters within jets
2557  //(centrality, cluster energy, jet pT, jet eta, jet phi)
2558  if (fPlotClustersInJets) {
2559  histname = TString::Format("%s/hClustersInJets", jets->GetArrayName().Data());
2560  Int_t nClusters = jet->GetNumberOfClusters();
2561  AliVCluster* clus;
2562  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
2563  clus = jet->Cluster(iClus);
2564  Double_t x[5] = {fCent, clus->E(), GetJetPt(jets, jet), jet->Eta(), jet->Phi_0_2pi()};
2565  fHistManager.FillTHnSparse(histname, x);
2566  }
2567 
2568  // Loop through clusters, and plot estimated shift in JES due to cluster bump
2569  // Only do for 0-10% centrality, and for EMCal/DCal
2570  Double_t eclus;
2571  Double_t shift;
2572  Double_t shiftSum = 0;
2573  if (fCent < 10.) {
2574  if (GetJetType(jet) > -0.5 && GetJetType(jet) < 1.5) {
2575  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
2576  clus = jet->Cluster(iClus);
2577  eclus = clus->E();
2578  if (eclus > 0.5) {
2579  shift = 0.79 * TMath::Exp(-0.5 * ((eclus - 3.81) / 1.50)*((eclus - 3.81) / 1.50) );
2580  shiftSum += shift;
2581  }
2582  }
2583  histname = TString::Format("%s/hCaloJESshift", jets->GetArrayName().Data());
2584  fHistManager.FillTH3(histname, GetJetType(jet), GetJetPt(jets, jet), shiftSum);
2585  }
2586  }
2587  }
2588 
2589  }
2590 
2591  }
2592 
2593  // Loop through all cells and fill histos
2594  Int_t sm;
2595  Int_t relid[4];
2596  Double_t patchSumEMCal[20] = {0.};
2597  Double_t patchSumPHOS[4] = {0.};
2598  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
2599 
2600  absId = fCaloCells->GetCellNumber(i);
2601  ecell = fCaloCells->GetCellAmplitude(absId);
2602 
2603  // Fill cell histo
2604  histname = TString::Format("Cells/hCellEnergyAll");
2605  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
2606 
2607  // Fill cell patch histo, per SM
2608  sm = fGeom->GetSuperModuleNumber(absId);
2609  if (sm >=0 && sm < 20) {
2610  patchSumEMCal[sm] += ecell;
2611  }
2612 
2613  }
2614 
2615  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
2616 
2617  absId = phosCaloCells->GetCellNumber(i);
2618  ecell = phosCaloCells->GetCellAmplitude(absId);
2619 
2620  // Fill cell histo
2621  histname = TString::Format("Cells/hCellEnergyAll");
2622  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
2623 
2624  // Fill cell patch histo, per SM
2625  fPHOSGeo->AbsToRelNumbering(absId, relid);
2626  sm = relid[0];
2627  if (sm >=1 && sm < 5) {
2628  patchSumPHOS[sm-1] += ecell;
2629  }
2630 
2631  }
2632 
2633  for (Int_t sm = 0; sm < 20; sm++) {
2634  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
2635  fHistManager.FillTH2(histname, patchSumEMCal[sm], fCent);
2636  }
2637 
2638  for (Int_t sm = 1; sm < 5; sm++) {
2639  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
2640  fHistManager.FillTH2(histname, patchSumPHOS[sm-1], fCent);
2641  }
2642 
2643 }
2644 
2645 //________________________________________________________________________
2646 Double_t AliAnalysisTaskEmcalDijetImbalance::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
2647 {
2648  Int_t AbsIdseed = -1;
2649  Double_t Eseed = 0;
2650  for (Int_t i = 0; i < cluster->GetNCells(); i++) {
2651  if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > Eseed) {
2652  Eseed = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
2653  AbsIdseed = cluster->GetCellAbsId(i);
2654  }
2655  }
2656 
2657  if (Eseed < 1e-9) {
2658  return 100;
2659  }
2660 
2661  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
2662  fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
2663  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
2664 
2665  //Get close cells index and energy, not in corners
2666 
2667  Int_t absID1 = -1;
2668  Int_t absID2 = -1;
2669 
2670  if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) {
2671  absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
2672  }
2673  if (iphi > 0) {
2674  absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
2675  }
2676 
2677  // In case of cell in eta = 0 border, depending on SM shift the cross cell index
2678 
2679  Int_t absID3 = -1;
2680  Int_t absID4 = -1;
2681 
2682  if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
2683  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
2684  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
2685  }
2686  else if (ieta == 0 && imod%2) {
2687  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
2688  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
2689  }
2690  else {
2691  if (ieta < AliEMCALGeoParams::fgkEMCALCols-1) {
2692  absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
2693  }
2694  if (ieta > 0) {
2695  absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
2696  }
2697  }
2698 
2699  Double_t ecell1 = cells->GetCellAmplitude(absID1);
2700  Double_t ecell2 = cells->GetCellAmplitude(absID2);
2701  Double_t ecell3 = cells->GetCellAmplitude(absID3);
2702  Double_t ecell4 = cells->GetCellAmplitude(absID4);
2703 
2704  Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
2705 
2706  Double_t Fcross = 1 - Ecross/Eseed;
2707 
2708  return Fcross;
2709 }
2710 
AliEmcalJet * GetLeadingJet(AliJetContainer *jetCont)
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
Double_t Area() const
Definition: AliEmcalJet.h:123
TObjArray fClusterCollArray
cluster collection array
Double_t GetRhoVal() const
Int_t fNPhiBins
Number of phi bins in DCal region (for background/correction)
const TString & GetRhoName() const
Int_t fNPtHistBins
! number of variable pt bins
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
PHOS acceptance.
Definition: AliEmcalJet.h:68
Bool_t fUseManualEventCuts
Flag to use manual event cuts.
UInt_t fOffTrigger
offline trigger for event selection
Bool_t fUseAliEventCuts
Flag to use AliEventCuts (otherwise AliAnalysisTaskEmcal will be used)
Double_t Eta() const
Definition: AliEmcalJet.h:114
static AliEmcalDownscaleFactorsOCDB * Instance()
const AliClusterIterableMomentumContainer accepted_momentum() const
Int_t fNEtaBins
Number of eta bins in DCal region (for background/correction)
bidirectional stl iterator over the EMCAL iterable container
Double_t Phi() const
Definition: AliEmcalJet.h:110
Container with name, TClonesArray and cuts for particles.
Bool_t fPlotJetHistograms
Set whether to enable inclusive jet histograms.
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.
TList * fEventCutList
! Output list for event cut histograms
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.
TH2D * fGapJetScalingWeights
Histogram storing eta-phi weights scaling jets near the gap region.
AliVParticle * Track(Int_t idx) const
Double_t fDijetLeadingHadronPt
leading hadron pT threshold for leading jet in dijet
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Bool_t fkEMCEJE
! flag telling whether the event is "triggered" or not in "simulation"
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
Double_t fMedianEMCal
! median patch energy in EMCal, per event
Double_t fMedianDCal
! median patch energy in DCal, per event
Double_t fMinAssJetPt
subleading jet min pT in a dijet pair, for it to be accepted
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
void SetCaloTriggerPatchInfoName(const char *n)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Bool_t fDoMomentumBalance
Set whether to enable momentum balance study.
void FindDijet(AliJetContainer *jetCont, Int_t leadingHadronCutBin)
Bool_t fPlotDijetImbalanceHistograms
Set whether to enable dijet imbalance histograms.
bool AddQAPlotsToList(TList *list)
AliEventCuts fEventCuts
event selection utility
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Bool_t fPlotClusWithoutNonLinCorr
If true, use pre-nonlincorr energy in cluster thnsparse.
const Int_t nPtBins
AliEmcalJet * GetLeadingJet(const char *opt="")
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
Dijet_t fDijet
! dijet candidate (per event)
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.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Double_t GetDeltaR(AliEmcalJet *jet1, AliEmcalJet *jet2)
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
float Float_t
Definition: External.C:68
Double_t fClusterConstituentThreshold
constituent threshold for matching study
Di-jet imbalance analysis with full jets.
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 fDoTriggerSimulation
Set whether to perform a simple trigger simulation.
BeamType fForceBeamType
forced beam type
Definition: External.C:228
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
virtual Bool_t AcceptCluster(Int_t i, UInt_t &rejectionReason) const
Double_t fDeltaPhiMin
minimum delta phi between di-jets
Double_t fMBUpscaleFactor
! inverse of downscale factor, for MB trigger
Bool_t fComputeBackground
Set whether to enable study of background.
Double_t Phi_0_2pi() const
Definition: AliEmcalJet.h:122
std::vector< TString > GetTriggerClasses() const
Double_t fCent
!event centrality
void LoadBackgroundScalingHistogram(const char *path="alien:///alice/cern.ch/user/j/jmulliga/scaleFactorEMCalLHC15o.root", const char *name1="hEtaPhiSFCorrection", const char *name2="hEtaPhiJetPtCorrection")
Double_t GetJetPt(AliJetContainer *jetCont, AliEmcalJet *jet)
Double_t GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
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 FillMomentumBalanceHistograms(TString histname, Double_t deltaPhi, Double_t trackPt, Double_t balancePt)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
TObjArray fJetCollArray
jet collection array
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
Bool_t fComputeMBDownscaling
Set whether to compute and plot MB downscaling factors.
Double_t Pt() const
Definition: AliEmcalJet.h:102
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
const AliClusterIterableMomentumContainer all_momentum() const
virtual Bool_t IsEventSelected()
Performing event selection.
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
Definition: External.C:220
Bool_t fLoadBackgroundScalingWeights
Flag to load eta-phi weights for full-jet background scale factors.
Bool_t fPlotClusterTHnSparse
Set whether to plot cluster thnsparse in calo studies.
Handler for downscale factors for various triggers obtained from the OCDB.
TFile * file
TList with histograms for a given trigger.
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Base task in the EMCAL jet framework.
Bool_t fDoGeometricalMatching
Set whether to enable constituent study with geometrical matching.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Bool_t fPlotClustersInJets
Set whether to plot histogram of clusters within jets.
const AliTrackIterableMomentumContainer accepted_momentum() const
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
Dijet_t fMatchingDijet
! low-threshold matching dijet, for matching study
TH2D * fBackgroundScalingWeights
Histogram storing eta-phi weights for full-jet background scale factors.
void SetRejectionReasonLabels(TAxis *axis)
void UserCreateOutputObjects()
Main initialization function on the worker.
AliEmcalEmbeddingQA fEmbeddingQA
! QA hists for embedding (will only be added if embedding)
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Bool_t fPlotExotics
Set whether to plot exotic cluster study.
Double_t fTrackConstituentThreshold
constituent threshold for matching study
DCal acceptance – spans ONLY DCal (no PHOS or gap)
Definition: AliEmcalJet.h:66
Double_t NEF() const
Definition: AliEmcalJet.h:141
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 fMinTrigJetPt
leading jet min pT in a dijet pair
Container structure for EMCAL clusters.
Bool_t fPlotNeutralJets
Set whether to plot neutral jet histo.
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Bool_t fPlotDijetCandHistograms
Set whether to enable dijet pair histograms.
TH3 * CreateTH3(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, int nbinsz, double zmin, double zmax, Option_t *opt="")
Create a new TH2 within the container.
Bool_t fDoCaloStudy
Set whether to perform calorimeter detector study.
const AliJetIterableContainer all() const
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal