AliPhysics  15d9304 (15d9304)
 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 
44 
48 
54  fHistManager(),
55  fEventCuts(0),
56  fEventCutList(0),
57  fUseManualEventCuts(kFALSE),
58  fUseAliEventCuts(kTRUE),
59  fDeltaPhiMin(0),
60  fMinTrigJetPt(0),
61  fMinAssJetPt(0),
62  fDijetLeadingHadronPt(0),
63  fMaxPt(200),
64  fNCentHistBins(0),
65  fCentHistBins(0),
66  fNPtHistBins(0),
67  fPtHistBins(0),
68  fPlotJetHistograms(kFALSE),
69  fPlotDijetCandHistograms(kFALSE),
70  fPlotDijetImbalanceHistograms(kFALSE),
71  fComputeBackground(kFALSE),
72  fDoMomentumBalance(kFALSE),
73  fDoGeometricalMatching(kFALSE),
74  fMatchingJetR(0.2),
75  fTrackConstituentThreshold(0),
76  fClusterConstituentThreshold(0),
77  fMBUpscaleFactor(1.),
78  fNEtaBins(40),
79  fNPhiBins(200),
80  fLoadBackgroundScalingWeights(kTRUE),
81  fBackgroundScalingWeights(0),
82  fGapJetScalingWeights(0),
83  fComputeMBDownscaling(kFALSE),
84  fDoCaloStudy(kFALSE),
85  fPHOSGeo(nullptr)
86 {
90 }
91 
98  AliAnalysisTaskEmcalJet(name, kTRUE),
99  fHistManager(name),
100  fEventCuts(0),
101  fEventCutList(0),
102  fUseManualEventCuts(kFALSE),
103  fUseAliEventCuts(kTRUE),
104  fDeltaPhiMin(0),
105  fMinTrigJetPt(0),
106  fMinAssJetPt(0),
107  fDijetLeadingHadronPt(0),
108  fMaxPt(200),
109  fNCentHistBins(0),
110  fCentHistBins(0),
111  fNPtHistBins(0),
112  fPtHistBins(0),
113  fPlotJetHistograms(kFALSE),
114  fPlotDijetCandHistograms(kFALSE),
115  fPlotDijetImbalanceHistograms(kFALSE),
116  fComputeBackground(kFALSE),
117  fDoMomentumBalance(kFALSE),
118  fDoGeometricalMatching(kFALSE),
119  fMatchingJetR(0.2),
120  fTrackConstituentThreshold(0),
121  fClusterConstituentThreshold(0),
122  fMBUpscaleFactor(1.),
123  fNEtaBins(40),
124  fNPhiBins(200),
125  fLoadBackgroundScalingWeights(kTRUE),
126  fBackgroundScalingWeights(0),
127  fGapJetScalingWeights(0),
128  fComputeMBDownscaling(kFALSE),
129  fDoCaloStudy(kFALSE),
130  fPHOSGeo(nullptr)
131 {
133  Dijet_t fDijet;
135 }
136 
141 {
142 }
143 
148 {
149  fNCentHistBins = 4;
151  fCentHistBins[0] = 0;
152  fCentHistBins[1] = 10;
153  fCentHistBins[2] = 30;
154  fCentHistBins[3] = 50;
155  fCentHistBins[4] = 90;
156 
157  fNPtHistBins = 82;
160  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
161  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
162  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
163  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
164  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
165  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
166 }
167 
173 {
175 
182 
183  TIter next(fHistManager.GetListOfHistograms());
184  TObject* obj = 0;
185  while ((obj = next())) {
186  fOutput->Add(obj);
187  }
188 
189  // Intialize AliEventCuts
190  if (fUseAliEventCuts) {
191  fEventCutList = new TList();
192  fEventCutList ->SetOwner();
193  fEventCutList ->SetName("EventCutOutput");
194 
195  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
196  if(fUseManualEventCuts==1)
197  {
198  fEventCuts.SetManualMode();
199  // Configure manual settings here
200  // ...
201  }
202  fEventCuts.AddQAplotsToList(fEventCutList);
203  fOutput->Add(fEventCutList);
204  }
205 
206  // Load eta-phi background scale factors from histogram on AliEn
209  }
210 
211  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
212 }
213 
214 /*
215  * This function allocates the histograms for single jets.
216  * A set of histograms is allocated per each jet container.
217  */
219 {
220  TString histname;
221  TString title;
222 
223  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
224 
225  AliJetContainer* jets = 0;
226  TIter nextJetColl(&fJetCollArray);
227  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
228 
229  // Jet rejection reason
230  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
231  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
232  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
233  SetRejectionReasonLabels(hist->GetXaxis());
234 
235  // Rho vs. Centrality
236  if (!jets->GetRhoName().IsNull()) {
237  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
238  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
239  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
240  }
241 
242  // Centrality vs. pT
243  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
244  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});Centrality (%);counts";
245  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 10, 0, 100);
246 
247  // NEF vs. pT vs. eta vs. phi
248  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
249  title = histname + ";#eta_{jet} (rad);#phi_{jet} (rad);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
250  Int_t nbins[4] = {fNEtaBins, fNPhiBins, nPtBins, 50};
251  Double_t min[4] = {-0.5,1., 0, 0};
252  Double_t max[4] = {0.5,6., fMaxPt, 1.};
253  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins, min, max);
254 
255  // pT upscaled
256  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
257  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c})";
258  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, "s");
259 
260  // pT-leading vs. pT
261  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
262  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
263  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
264 
265  // A vs. pT
266  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
267  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
268  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/3, 0, 0.5);
269 
270  // z-leading (charged) vs. pT
271  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
272  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
273  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/5, 0, 1.0);
274 
275  // z (charged) vs. pT
276  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
277  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
278  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/5, 0, 1.0);
279 
280  // Nconst vs. pT
281  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
282  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
283  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/5, 0, fMaxPt);
284 
285  // Allocate background subtraction histograms, if enabled
286  if (fComputeBackground) {
287 
288  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
289  title = histname + ";Centrality;Scale factor;counts";
290  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
291 
292  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
293  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
294  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
295 
296  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jets->GetArrayName().Data());
297  title = histname + ";#eta;#phi;Centrality;Scale factor;";
298  Int_t nbins[4] = {fNEtaBins, fNPhiBins, 50, 400};
299  Double_t min[4] = {-0.5,1., 0, 0};
300  Double_t max[4] = {0.5,6., 100, 20};
301  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins, min, max);
302 
303  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jets->GetArrayName().Data());
304  title = histname + ";#eta;#phi;Centrality;#delta#it{p}_{T} (GeV/#it{c})";
305  Int_t nbinsDpT[4] = {fNEtaBins, fNPhiBins, 10, 400};
306  Double_t minDpT[4] = {-0.5,1., 0, -50};
307  Double_t maxDpT[4] = {0.5,6., 100, 150};
308  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbinsDpT, minDpT, maxDpT);
309 
310  }
311 
312  }
313 
314  // MB downscale factor histogram
315  if (fComputeMBDownscaling) {
316  histname = "Trigger/hMBDownscaleFactor";
317  title = histname + ";Downscale factor;counts";
318  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 200);
319  }
320 
321 }
322 
323 /*
324  * This function allocates the histograms for dijet candidates, i.e. dijet pairs with the leading jet
325  * passing the trigger condition, but no condition on the subleading jet. In particular this histogram is
326  * designed to study the kinematic selection of dijets.
327  */
329 {
330  // Allocate dijet THnSparse
331  AliJetContainer* jets = 0;
332  TIter nextJetColl(&fJetCollArray);
333  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
334  TString axisTitle[30]= {""};
335  Int_t nbins[30] = {0};
336  Double_t min[30] = {0.};
337  Double_t max[30] = {0.};
338  Double_t *binEdges[20] = {0};
339  Int_t dim = 0;
340 
341  if (fForceBeamType != kpp) {
342  axisTitle[dim] = "Centrality (%)";
343  nbins[dim] = fNCentHistBins;
344  binEdges[dim] = fCentHistBins;
345  min[dim] = fCentHistBins[0];
346  max[dim] = fCentHistBins[fNCentHistBins];
347  dim++;
348  }
349 
350  axisTitle[dim] = "LeadingHadronRequired";
351  nbins[dim] = 2;
352  min[dim] = -0.5;
353  max[dim] = 1.5;
354  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
355  dim++;
356 
357  axisTitle[dim] = "#it{p}_{T,trig jet} (GeV/#it{c})";
358  nbins[dim] = TMath::CeilNint(fMaxPt/2);
359  min[dim] = 0;
360  max[dim] = fMaxPt;
361  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
362  dim++;
363 
364  axisTitle[dim] = "#it{p}_{T,ass jet} (GeV/#it{c})";
365  nbins[dim] = TMath::CeilNint(fMaxPt/2);
366  min[dim] = 0;
367  max[dim] = fMaxPt;
368  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
369  dim++;
370 
371  axisTitle[dim] = "#phi_{trig jet}";
372  nbins[dim] = TMath::CeilNint(fMaxPt/2);
373  min[dim] = 0;
374  max[dim] = TMath::TwoPi();
375  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
376  dim++;
377 
378  axisTitle[dim] = "#phi_{ass jet}";
379  nbins[dim] = TMath::CeilNint(fMaxPt/2);
380  min[dim] = 0;
381  max[dim] = TMath::TwoPi();
382  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
383  dim++;
384 
385  axisTitle[dim] = "#eta_{trig jet}";
386  nbins[dim] = TMath::CeilNint(fMaxPt/2);
387  min[dim] = -1;
388  max[dim] = 1;
389  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
390  dim++;
391 
392  axisTitle[dim] = "#eta_{ass jet}";
393  nbins[dim] = TMath::CeilNint(fMaxPt/2);
394  min[dim] = -1;
395  max[dim] = 1;
396  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
397  dim++;
398 
399  TString thnname = TString::Format("%s/DijetCandObservables", jets->GetArrayName().Data());
400  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
401  for (Int_t i = 0; i < dim; i++) {
402  hn->GetAxis(i)->SetTitle(axisTitle[i]);
403  hn->SetBinEdges(i, binEdges[i]);
404  }
405  }
406 }
407 
408 /*
409  * This function allocates the histograms for accepted dijets.
410  * The purpose is to study in detail the imbalance properties of the accepted dijets.
411  */
413 {
414  // Allocate dijet imbalance THnSparse
415  AliJetContainer* jets = 0;
416  TIter nextJetColl(&fJetCollArray);
417  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
418  TString axisTitle[30]= {""};
419  Int_t nbins[30] = {0};
420  Double_t min[30] = {0.};
421  Double_t max[30] = {0.};
422  Double_t *binEdges[20] = {0};
423  Int_t dim = 0;
424 
425  if (fForceBeamType != kpp) {
426  axisTitle[dim] = "Centrality (%)";
427  nbins[dim] = fNCentHistBins;
428  binEdges[dim] = fCentHistBins;
429  min[dim] = fCentHistBins[0];
430  max[dim] = fCentHistBins[fNCentHistBins];
431  dim++;
432  }
433 
434  axisTitle[dim] = "#Delta#phi";
435  nbins[dim] = 100;
436  min[dim] = 0;
437  max[dim] = 4;
438  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
439  dim++;
440 
441  axisTitle[dim] = "#Delta#eta";
442  nbins[dim] = 100;
443  min[dim] = -2;
444  max[dim] = 2;
445  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
446  dim++;
447 
448  axisTitle[dim] = "A_{J}";
449  nbins[dim] = 100;
450  min[dim] = 0;
451  max[dim] = 1;
452  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
453  dim++;
454 
455  axisTitle[dim] = "x_{J}";
456  nbins[dim] = 100;
457  min[dim] = 0;
458  max[dim] = 1;
459  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
460  dim++;
461 
462  axisTitle[dim] = "k_{Ty} (GeV)";
463  nbins[dim] = 100;
464  min[dim] = 0;
465  max[dim] = 100;
466  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
467  dim++;
468 
469  axisTitle[dim] = "N_{tracks, trig jet}";
470  nbins[dim] = fMaxPt/5;
471  min[dim] = 0;
472  max[dim] = 100;
473  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
474  dim++;
475 
476  axisTitle[dim] = "N_{tracks, ass jet}";
477  nbins[dim] = fMaxPt/5;
478  min[dim] = 0;
479  max[dim] = 100;
480  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
481  dim++;
482 
483  TString thnname = TString::Format("%s/DijetImbalanceObservables", jets->GetArrayName().Data());
484  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
485  for (Int_t i = 0; i < dim; i++) {
486  hn->GetAxis(i)->SetTitle(axisTitle[i]);
487  hn->SetBinEdges(i, binEdges[i]);
488  }
489  }
490 }
491 
492 /*
493  * This function allocates the histograms for the momentum balance study.
494  */
496 {
497  AliJetContainer* jets = 0;
498  TIter nextJetColl(&fJetCollArray);
499  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
500 
501  // Allocate THnSparse
502  TString axisTitle[30]= {""};
503  Int_t nbins[30] = {0};
504  Double_t min[30] = {0.};
505  Double_t max[30] = {0.};
506  Double_t *binEdges[20] = {0};
507  Int_t dim = 0;
508 
509  axisTitle[dim] = "A_{J}";
510  nbins[dim] = 100;
511  min[dim] = 0;
512  max[dim] = 1;
513  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
514  dim++;
515 
516  axisTitle[dim] = "#Delta#phi";
517  nbins[dim] = 100;
518  min[dim] = -4;
519  max[dim] = 4;
520  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
521  dim++;
522 
523  axisTitle[dim] = "#it{p}_{T,particle} (GeV/#it{c})";
524  nbins[dim] = 9;
525  Double_t* pTParticleBins = new Double_t[nbins[dim]+1];
526  GenerateFixedBinArray(1, 0.15, 0.3, pTParticleBins);
527  GenerateFixedBinArray(1, 0.3, 0.5, pTParticleBins+1);
528  GenerateFixedBinArray(1, 0.5, 1, pTParticleBins+2);
529  GenerateFixedBinArray(2, 1, 5, pTParticleBins+3);
530  GenerateFixedBinArray(3, 5, 20, pTParticleBins+5);
531  GenerateFixedBinArray(1, 20, 150, pTParticleBins+8);
532  min[dim] = 0;
533  max[dim] = pTParticleBins[nbins[dim]];
534  binEdges[dim] = pTParticleBins;
535  dim++;
536 
537  axisTitle[dim] = "#it{p}_{T}#parallel (GeV/#it{c})";
538  nbins[dim] = 80;
539  Double_t* pTParallelBins = new Double_t[nbins[dim]+1];
540  GenerateFixedBinArray(20, 0, 2, pTParallelBins);
541  GenerateFixedBinArray(16, 2, 10, pTParallelBins+20);
542  GenerateFixedBinArray(10, 10, 20, pTParallelBins+36);
543  GenerateFixedBinArray(10, 20, 40, pTParallelBins+46);
544  GenerateFixedBinArray(24, 40, 150, pTParallelBins+56);
545  min[dim] = 0;
546  max[dim] = pTParallelBins[nbins[dim]];
547  binEdges[dim] = pTParallelBins;
548  dim++;
549 
550  TString thnname = TString::Format("%s/MomentumBalance", jets->GetArrayName().Data());
551  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
552  for (Int_t i = 0; i < dim; i++) {
553  hn->GetAxis(i)->SetTitle(axisTitle[i]);
554  hn->SetBinEdges(i, binEdges[i]);
555  }
556  }
557 }
558 
559 /*
560  * This function allocates the histograms for the constituent dijet study.
561  */
563 {
564  // Allocate geometrical matching THnSparse
565  TString axisTitle[30]= {""};
566  Int_t nbins[30] = {0};
567  Double_t min[30] = {0.};
568  Double_t max[30] = {0.};
569  Double_t *binEdges[20] = {0};
570  Int_t dim = 0;
571 
572  if (fForceBeamType != kpp) {
573  axisTitle[dim] = "Centrality (%)";
574  nbins[dim] = fNCentHistBins;
575  binEdges[dim] = fCentHistBins;
576  min[dim] = fCentHistBins[0];
577  max[dim] = fCentHistBins[fNCentHistBins];
578  dim++;
579  }
580 
581  axisTitle[dim] = "isSwitched";
582  nbins[dim] = 2;
583  min[dim] = -0.5;
584  max[dim] = 1.5;
585  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
586  dim++;
587 
588  axisTitle[dim] = "#DeltaR_{trig}";
589  nbins[dim] = 50;
590  min[dim] = 0;
591  max[dim] = 0.5;
592  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
593  dim++;
594 
595  axisTitle[dim] = "#DeltaR_{ass}";
596  nbins[dim] = 50;
597  min[dim] = 0;
598  max[dim] = 0.5;
599  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
600  dim++;
601 
602  axisTitle[dim] = "trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
603  nbins[dim] = 100;
604  min[dim] = -50;
605  max[dim] = 50;
606  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
607  dim++;
608 
609  axisTitle[dim] = "ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
610  nbins[dim] = 100;
611  min[dim] = -50;
612  max[dim] = 50;
613  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
614  dim++;
615 
616  axisTitle[dim] = "A_{J} low-threshold";
617  nbins[dim] = 100;
618  min[dim] = 0;
619  max[dim] = 1;
620  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
621  dim++;
622 
623  axisTitle[dim] = "A_{J} hard-core";
624  nbins[dim] = 100;
625  min[dim] = 0;
626  max[dim] = 1;
627  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
628  dim++;
629 
630  TString thnname = "GeometricalMatching";
631  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
632  for (Int_t i = 0; i < dim; i++) {
633  hn->GetAxis(i)->SetTitle(axisTitle[i]);
634  hn->SetBinEdges(i, binEdges[i]);
635  }
636 
637  // Allocate other histograms
638  TString histname;
639  TString title;
640  histname = "GeometricalMatchingEfficiency";
641  title = histname + ";isMatched;counts";
642  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 2, -0.5, 1.5);
643 }
644 
645 /*
646  * This function allocates the histograms for the calorimeter performance study.
647  */
649 {
650  TString histname;
651  TString htitle;
652 
653  Double_t* clusType = new Double_t[3+1];
654  GenerateFixedBinArray(3, -0.5, 2.5, clusType);
655  const Int_t nRejBins = 32;
656  Double_t* rejReasonBins = new Double_t[nRejBins+1];
657  GenerateFixedBinArray(nRejBins, 0, nRejBins, rejReasonBins);
658 
659  AliEmcalContainer* cont = 0;
660  TIter nextClusColl(&fClusterCollArray);
661  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
662 
663  // rejection reason plot, to make efficiency correction
664  histname = TString::Format("%s/hClusterRejectionReasonEMCal", cont->GetArrayName().Data());
665  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
666  TH2* hist = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
667  SetRejectionReasonLabels(hist->GetXaxis());
668 
669  histname = TString::Format("%s/hClusterRejectionReasonPHOS", cont->GetArrayName().Data());
670  htitle = histname + ";Rejection reason;#it{E}_{clus} (GeV/)";
671  TH2* histPhos = fHistManager.CreateTH2(histname.Data(), htitle.Data(), nRejBins, rejReasonBins, fNPtHistBins, fPtHistBins);
672  SetRejectionReasonLabels(histPhos->GetXaxis());
673 
674  // plot by SM
675  const Int_t nEmcalSM = 20;
676  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
677  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
678  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
679  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
680  }
681 
682  for (Int_t sm = 1; sm < 5; sm++) {
683  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
684  htitle = histname + ";#it{E}_{cluster} (GeV);counts";
685  fHistManager.CreateTH1(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins);
686  }
687 
688  // Plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
689  Int_t dim = 0;
690  TString title[20];
691  Int_t nbins[20] = {0};
692  Double_t min[30] = {0.};
693  Double_t max[30] = {0.};
694  Double_t *binEdges[20] = {0};
695 
697  title[dim] = "Centrality %";
698  nbins[dim] = fNCentHistBins;
699  binEdges[dim] = fCentHistBins;
700  min[dim] = fCentHistBins[0];
701  max[dim] = fCentHistBins[fNCentHistBins];
702  dim++;
703  }
704 
705  title[dim] = "#eta";
706  nbins[dim] = 28;
707  min[dim] = -0.7;
708  max[dim] = 0.7;
709  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
710  dim++;
711 
712  title[dim] = "#phi";
713  nbins[dim] = 100;
714  min[dim] = 1.;
715  max[dim] = 6.;
716  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
717  dim++;
718 
719  title[dim] = "#it{E}_{clus} (GeV)";
720  nbins[dim] = fNPtHistBins;
721  binEdges[dim] = fPtHistBins;
722  min[dim] = fPtHistBins[0];
723  max[dim] = fPtHistBins[fNPtHistBins];
724  dim++;
725 
726  title[dim] = "#it{E}_{clus, hadcorr} (GeV)";
727  nbins[dim] = fNPtHistBins;
728  binEdges[dim] = fPtHistBins;
729  min[dim] = fPtHistBins[0];
730  max[dim] = fPtHistBins[fNPtHistBins];
731  dim++;
732 
733  title[dim] = "#it{E}_{core} (GeV)";
734  nbins[dim] = fNPtHistBins;
735  binEdges[dim] = fPtHistBins;
736  min[dim] = fPtHistBins[0];
737  max[dim] = fPtHistBins[fNPtHistBins];
738  dim++;
739 
740  title[dim] = "Matched track";
741  nbins[dim] = 2;
742  min[dim] = -0.5;
743  max[dim] = 1.5;
744  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
745  dim++;
746 
747  title[dim] = "M02";
748  nbins[dim] = 100;
749  min[dim] = 0;
750  max[dim] = 10;
751  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
752  dim++;
753 
754  title[dim] = "Ncells";
755  nbins[dim] = 40;
756  min[dim] = -0.5;
757  max[dim] = 49.5;
758  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
759  dim++;
760 
761  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
762  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
763  for (Int_t i = 0; i < dim; i++) {
764  hn->GetAxis(i)->SetTitle(title[i]);
765  hn->SetBinEdges(i, binEdges[i]);
766  }
767  }
768 
769  // plot neutral jets
770  AliJetContainer* jets = 0;
771  TIter nextJetColl(&fJetCollArray);
772  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
773 
774  TString axisTitle[30]= {""};
775  Int_t nbinsJet[30] = {0};
776  Double_t minJet[30] = {0.};
777  Double_t maxJet[30] = {0.};
778  Double_t *binEdgesJet[20] = {0};
779  Int_t dimJet = 0;
780 
781  if (fForceBeamType != kpp) {
782  axisTitle[dimJet] = "Centrality (%)";
783  nbinsJet[dimJet] = fNCentHistBins;
784  binEdgesJet[dimJet] = fCentHistBins;
785  minJet[dimJet] = fCentHistBins[0];
786  maxJet[dimJet] = fCentHistBins[fNCentHistBins];
787  dimJet++;
788  }
789 
790  axisTitle[dimJet] = "#eta_{jet}";
791  nbinsJet[dimJet] = 28;
792  minJet[dimJet] = -0.7;
793  maxJet[dimJet] = 0.7;
794  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
795  dimJet++;
796 
797  axisTitle[dimJet] = "#phi_{jet} (rad)";
798  nbinsJet[dimJet] = 100;
799  minJet[dimJet] = 1.;
800  maxJet[dimJet] = 6.;
801  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
802  dimJet++;
803 
804  axisTitle[dimJet] = "#it{E}_{T} (GeV)";
805  nbinsJet[dimJet] = fNPtHistBins;
806  binEdgesJet[dimJet] = fPtHistBins;
807  minJet[dimJet] = fPtHistBins[0];
808  maxJet[dimJet] = fPtHistBins[fNPtHistBins];
809  dimJet++;
810 
811  axisTitle[dimJet] = "#rho (GeV/#it{c})";
812  nbinsJet[dimJet] = 100;
813  minJet[dimJet] = 0.;
814  maxJet[dimJet] = 1000.;
815  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
816  dimJet++;
817 
818  axisTitle[dimJet] = "N_{clusters}";
819  nbinsJet[dimJet] = 20;
820  minJet[dimJet] = -0.5;
821  maxJet[dimJet] = 19.5;
822  binEdgesJet[dimJet] = GenerateFixedBinArray(nbinsJet[dimJet], minJet[dimJet], maxJet[dimJet]);
823  dimJet++;
824 
825  TString thnname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
826  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dimJet, nbinsJet, minJet, maxJet);
827  for (Int_t i = 0; i < dimJet; i++) {
828  hn->GetAxis(i)->SetTitle(axisTitle[i]);
829  hn->SetBinEdges(i, binEdgesJet[i]);
830  }
831 
832  }
833 
834  // Plot cell histograms
835 
836  // centrality vs. cell energy vs. cell type (for all cells)
837  histname = TString::Format("Cells/hCellEnergyAll");
838  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
839  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
840 
841  // centrality vs. cell energy vs. cell type (for cells in accepted clusters)
842  histname = TString::Format("Cells/hCellEnergyAccepted");
843  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
844  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
845 
846  // centrality vs. cell energy vs. cell type (for leading cells in accepted clusters)
847  histname = TString::Format("Cells/hCellEnergyLeading");
848  htitle = histname + ";#it{E}_{cell} (GeV);Centrality (%); Cluster type";
849  fHistManager.CreateTH3(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins, 3, clusType);
850 
851  // plot cell patches by SM
852  const Int_t nEmcalSM = 20;
853  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
854  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
855  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
856  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
857  }
858 
859  for (Int_t sm = 1; sm < 5; sm++) {
860  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
861  htitle = histname + ";#it{E}_{cell patch} (GeV);Centrality (%)";
862  fHistManager.CreateTH2(histname.Data(), htitle.Data(), fNPtHistBins, fPtHistBins, fNCentHistBins, fCentHistBins);
863  }
864 
865 }
866 
871 {
872 
873  TString fname(path);
874  if (fname.BeginsWith("alien://")) {
875  TGrid::Connect("alien://");
876  }
877 
878  TFile* file = TFile::Open(path);
879 
880  if (!file || file->IsZombie()) {
881  ::Error("AliAnalysisTaskEmcalDijetImbalance", "Could not open background scaling histogram");
882  return;
883  }
884 
885  TH1D* h = dynamic_cast<TH1D*>(file->Get(name));
886 
887  if (h) {
888  ::Info("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s loaded from file %s.", name, path);
889  }
890  else {
891  ::Error("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s not found in file %s.", name, path);
892  return;
893  }
894 
895  fBackgroundScalingWeights = static_cast<TH1D*>(h->Clone());
896 
897  file->Close();
898  delete file;
899 
900 }
901 
907 {
909 
910  fNeedEmcalGeom = kTRUE;
911 
912  // Load the PHOS geometry
913  if (fDoCaloStudy) {
914  fPHOSGeo = AliPHOSGeometry::GetInstance();
915  if (fPHOSGeo) {
916  AliInfo("Found instance of PHOS geometry!");
917  }
918  else {
919  AliInfo("Creating PHOS geometry!");
920  Int_t runNum = InputEvent()->GetRunNumber();
921  if(runNum<209122) //Run1
922  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
923  else
924  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
925 
926  if (fPHOSGeo) {
927  AliOADBContainer geomContainer("phosGeo");
928  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
929  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
930  for(Int_t mod=0; mod<6; mod++) {
931  if(!matrixes->At(mod)) continue;
932  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
933  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
934  ((TGeoHMatrix*)matrixes->At(mod))->Print();
935  }
936  }
937  }
938  }
939 
940  AliInfo(Form("Trigger jet threshold = %f, Associated jet threshold = %f", fMinTrigJetPt, fMinAssJetPt));
941  AliInfo(Form("Leading hadron threshold (for dijet leading jet): %f GeV", fDijetLeadingHadronPt));
942  AliInfo(Form("Momentum balance study: %d", fDoMomentumBalance));
943  AliInfo(Form("Geometrical matching study: %d", fDoGeometricalMatching));
944 
945 }
946 
951 
952  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
953 
955 
956  // Get instance of the downscale factor helper class
958  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
959 
960  // There are two possible min bias triggers for LHC15o
961  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
962  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
963  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
964 
965  // Get the downscale factor for whichever MB trigger exists in the given run
966  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
967  Double_t downscalefactor;
968  for (auto i : runtriggers) {
969  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
970  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
971  break;
972  }
973  }
974 
975  // Store the inverse of the downscale factor, used later to weight the pT spectrum
976  fMBUpscaleFactor = 1/downscalefactor;
977 
978  TString histname = "Trigger/hMBDownscaleFactor";
979  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
980 
981  }
982 
983 }
984 
989 {
990  if (fUseAliEventCuts) {
991  if (!fEventCuts.AcceptEvent(InputEvent()))
992  {
993  PostData(1, fOutput);
994  return kFALSE;
995  }
996  }
997  else {
999  }
1000  return kTRUE;
1001 }
1002 
1011 {
1012  TString histname;
1013  AliJetContainer* jetCont = 0;
1014  TIter next(&fJetCollArray);
1015  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1016  TString jetContName = jetCont->GetName();
1017  if (jetContName.Contains("HardCore")) continue;
1018 
1019  //-----------------------------------------------------------------------------
1020  // Find the leading di-jet candidate in each event, and if it satisfies the
1021  // trig jet pT threshold, then fill di-jet candidate histogram (regardless of ass jet).
1022  // The idea is to study the kinematic selections in post-processing.
1023 
1024  // Loop over leading hadron cut or not
1025  for (Int_t leadingHadronCutType=0; leadingHadronCutType<2; leadingHadronCutType++) {
1026 
1027  // Find the dijet candidate of the event and store its info in struct fDijet
1028  FindDijet(jetCont, leadingHadronCutType);
1029 
1030  // If we find a dijet candidate (i.e. acceptable trig jet; ass jet accepted or not), fill the di-jet candidate histogram
1032  TString histname = TString::Format("%s/DijetCandObservables", jetCont->GetArrayName().Data());
1033  FillDijetCandHistograms(histname);
1034  }
1035 
1036  }
1037 
1038  //---------------------------------------------------------------------------------------------------
1039  // Now, study the accepted dijet selection -- specified by the trig/ass jet pT conditions
1040 
1041  // Find the dijet candidate of the event and store its info in struct fDijet
1042  FindDijet(jetCont, 0);
1043 
1044  // If we find an accepted dijet, fill the dijet imbalance histogram
1046  TString histname = TString::Format("%s/DijetImbalanceObservables", jetCont->GetArrayName().Data());
1047  FillDijetImbalanceHistograms(histname);
1048  }
1049 
1050  // If we find an accepted dijet, perform momentum-balance study (if requested)
1052  histname = TString::Format("%s/MomentumBalance", jetCont->GetArrayName().Data());
1053  DoMomentumBalance(histname);
1054  }
1055 
1056  }
1057 
1058  //---------------------------------------------------------------------------
1059  // Do the constituent threshold and geometrical matching study (if requested)
1062 
1063  return kTRUE;
1064 }
1065 
1074 {
1075  fDijet.clear();
1076  fDijet.leadingHadronCutType = leadingHadronCutBin;
1077 
1078  // Get trigger jet
1079  AliEmcalJet* trigJet = GetLeadingJet(jetCont);
1080  if(!trigJet) return;
1081 
1082  // Skip the event if the leading jet doesn't satisfy the pT threshold
1083  Double_t trigJetPt = GetJetPt(jetCont, trigJet);
1084  if ( trigJetPt < fMinTrigJetPt ) return;
1085 
1086  // Skip the event if the leading jet doesn't satisfy the leading hadron threshold
1087  if (jetCont->GetLeadingHadronPt(trigJet) < fDijetLeadingHadronPt*leadingHadronCutBin) return;
1088 
1089  // Fill the dijet struct
1090  fDijet.trigJet = trigJet;
1091  fDijet.trigJetPt = trigJetPt;
1092  fDijet.trigJetEta = trigJet->Eta();
1093  fDijet.trigJetPhi = trigJet->Phi();
1094 
1095  // Find the subleading jet in the opposite hemisphere
1096  AliEmcalJet *assJet = 0;
1097  for(auto assJetCand : jetCont->accepted()) {
1098  if (!assJetCand) continue;
1099  Double_t assJetCandPt = GetJetPt(jetCont, assJetCand);
1100  if ( TMath::Abs(trigJet->Phi() - assJetCand->Phi()) < fDeltaPhiMin ) continue;
1101  if (assJet) {
1102  Double_t assJetPt = GetJetPt(jetCont, assJet);
1103  if ( assJetCandPt < assJetPt ) continue;
1104  }
1105  assJet = assJetCand;
1106  }
1107  if (!assJet) return;
1108 
1109  // Fill the dijet struct
1110  fDijet.assJet = assJet;
1111  fDijet.assJetPt = GetJetPt(jetCont, assJet);
1112  fDijet.assJetPhi = assJet->Phi();
1113  fDijet.assJetEta = assJet->Eta();
1115 
1116  fDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
1117  fDijet.deltaEta = trigJet->Eta() - assJet->Eta();
1120  fDijet.kTy = TMath::Abs( fDijet.trigJetPt * TMath::Sin(fDijet.deltaPhi) );
1121 }
1122 
1128 {
1129 
1130  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1131 
1132  AliVTrack* track;
1133  for (auto trackIterator : trackCont->accepted_momentum() ) {
1134 
1135  track = trackIterator.second;
1136 
1137  // Compute the delta phi between the track and its nearest jet (of the two jets in the dijet),
1138  // as well as its pT-parallel projection onto the nearest jet's axis.
1139 
1140  Double_t trackPt = track->Pt();
1141  Double_t trackPhi = track->Phi();
1142  Double_t trigJetPhi = fDijet.trigJet->Phi();
1143  Double_t assJetPhi = fDijet.assJet->Phi();
1144 
1145  Double_t deltaPhiTrigJet = TMath::Abs(trackPhi - trigJetPhi);
1146  Double_t deltaPhiAssJet = TMath::Abs(trackPhi - assJetPhi);
1147  Bool_t isNearside = deltaPhiTrigJet < deltaPhiAssJet;
1148 
1149  Double_t deltaPhi;
1150  Double_t balancePt;
1151  if (isNearside) {
1152  deltaPhi = trackPhi - trigJetPhi;
1153  balancePt = trackPt * TMath::Cos(deltaPhi);
1154  }
1155  else {
1156  deltaPhi = trackPhi - assJetPhi;
1157  balancePt = -trackPt * TMath::Cos(deltaPhi);
1158  }
1159 
1160  FillMomentumBalanceHistograms(histname, deltaPhi, trackPt, balancePt);
1161 
1162  }
1163 }
1164 
1169 {
1170  // Get jet container with minimum constituent pT,E thresholds
1171  TString jetContAllName = Form("Jet_AKTFullR0%d0_tracks_pT0150_caloClusters_E0300_pt_scheme", (int) (fMatchingJetR*10) );
1172  AliJetContainer* jetContAll = GetJetContainer(jetContAllName.Data());
1173 
1174  // Get jet container with X GeV constituent pT,E thresholds
1175  Int_t trackThreshold = (int) (fTrackConstituentThreshold*1000); // in MeV
1176  Int_t clusThreshold = (int) (fClusterConstituentThreshold*1000); // in MeV
1177  TString jetContHardCoreName = Form("JetHardCore_AKTFullR0%d0_tracks_pT%d_caloClusters_E%d_pt_scheme", (int) (fMatchingJetR*10), trackThreshold, clusThreshold);
1178  AliJetContainer* jetContHardCore = GetJetContainer(jetContHardCoreName.Data());
1179 
1180  // Find the di-jet in the hard-core jet sample, then find the matching di-jet and fill histograms
1181  FindDijet(jetContHardCore, 0);
1182  if (fDijet.isAccepted) {
1183  FindMatchingDijet(jetContAll);
1185  }
1186 }
1187 
1193 {
1195 
1196  // Loop over jets and find leading jet within R of fDijet.trigJet
1197  AliEmcalJet *matchingTrigJet = 0;
1198  for(auto matchingTrigJetCand : jetCont->accepted()) {
1199  if (!matchingTrigJetCand) continue;
1200  if ( GetDeltaR(matchingTrigJetCand, fDijet.trigJet) > fMatchingJetR ) continue;
1201  if (matchingTrigJet) {
1202  if ( GetJetPt(jetCont, matchingTrigJetCand) < GetJetPt(jetCont, matchingTrigJet) ) continue;
1203  }
1204  matchingTrigJet = matchingTrigJetCand;
1205  }
1206  if (!matchingTrigJet) return;
1207 
1208  // Loop over jets and find leading jet within R of fDijet.assJet
1209  AliEmcalJet *matchingAssJet = 0;
1210  for(auto matchingAssJetCand : jetCont->accepted()) {
1211  if (!matchingAssJetCand) continue;
1212  if ( GetDeltaR(matchingAssJetCand, fDijet.assJet) > fMatchingJetR ) continue;
1213  if (matchingAssJet) {
1214  if ( GetJetPt(jetCont, matchingAssJetCand) < GetJetPt(jetCont, matchingAssJet) ) continue;
1215  }
1216  matchingAssJet = matchingAssJetCand;
1217  }
1218 
1219  // Determine which matching jet is the leading jet (i.e. allow them to flip)
1220  if (matchingAssJet) {
1221  AliEmcalJet* trigJet = matchingTrigJet;
1222  AliEmcalJet* assJet = matchingAssJet;
1223  if ( GetJetPt(jetCont, matchingTrigJet) < GetJetPt(jetCont, matchingAssJet) ) {
1224  AliEmcalJet* trigJet = matchingAssJet;
1225  AliEmcalJet* assJet = matchingTrigJet;
1226  }
1227 
1228  // Fill the dijet struct
1229  fMatchingDijet.trigJet = trigJet;
1230  fMatchingDijet.trigJetPt = GetJetPt(jetCont, trigJet);
1231  fMatchingDijet.trigJetEta = trigJet->Eta();
1232  fMatchingDijet.trigJetPhi = trigJet->Phi();
1233 
1234  fMatchingDijet.assJet = assJet;
1235  fMatchingDijet.assJetPt = GetJetPt(jetCont, assJet);
1236  fMatchingDijet.assJetPhi = assJet->Phi();
1237  fMatchingDijet.assJetEta = assJet->Eta();
1239 
1240  fMatchingDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
1241  fMatchingDijet.deltaEta = trigJet->Eta() - assJet->Eta();
1244  fMatchingDijet.kTy = TMath::Abs( fMatchingDijet.trigJetPt * TMath::Sin(fMatchingDijet.deltaPhi) );
1245  }
1246 }
1247 
1252 {
1253  // Loop over tracks and clusters in order to:
1254  // (1) Compute scale factor for full jets
1255  // (2) Compute delta-pT for full jets, with the random cone method
1256  // For both the scale factor and delta-pT, we compute only one histogram each for EMCal.
1257  // But for DCal, we bin in eta-phi, in order to account for the DCal vs. PHOS vs. gap
1258 
1259  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
1260  Double_t etaTPC = 0.9;
1261  Double_t etaEMCal = 0.7;
1262  Double_t etaMinDCal = 0.22;
1263  Double_t etaMaxPHOS = 0.13;
1264  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
1265  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
1266  Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
1267  Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
1268  Double_t phiMinPHOS = 250 * TMath::DegToRad();
1269  Double_t phiMaxPHOS = 320 * TMath::DegToRad();
1270 
1271  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
1272  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
1273  Double_t accDCalRegion = 2 * etaEMCal * (phiMaxDCal - phiMinDCal);
1274 
1275  // Define fiducial acceptances, to be used to generate random cones
1276  TRandom3* r = new TRandom3(0);
1277  Double_t jetR = jetCont->GetJetRadius();
1278  Double_t accRC = TMath::Pi() * jetR * jetR;
1279  Double_t etaEMCalfid = etaEMCal - jetR;
1280  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
1281  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
1282  Double_t phiMinDCalRegionfid = phiMinDCal + jetR;
1283  Double_t phiMaxDCalRegionfid = phiMaxDCal - jetR;
1284 
1285  // Generate EMCal random cone eta-phi
1286  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
1287  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
1288 
1289  // For DCalRegion, generate random eta, phi in each eta/phi bin, to be used as center of random cone
1290  Double_t etaDCalRC[fNEtaBins]; // array storing the RC eta values
1291  Double_t etaStep = 1./fNEtaBins;
1292  Double_t etaMin;
1293  Double_t etaMax;
1294  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1295  etaMin = -etaEMCalfid + etaBin*etaStep;
1296  etaMax = etaMin + etaStep;
1297  etaDCalRC[etaBin] = r->Uniform(etaMin, etaMax);
1298  }
1299 
1300  Double_t phiDCalRC[fNPhiBins]; // array storing the RC phi values
1301  Double_t phiStep = 5./fNPhiBins; // phi axis is [1,6] in order to have simple binning
1302  Double_t phiMin;
1303  Double_t phiMax;
1304  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1305  phiMin = 1 + phiBin*phiStep;
1306  phiMax = phiMin + phiStep;
1307  phiDCalRC[phiBin] = r->Uniform(phiMin, phiMax);
1308  }
1309 
1310  // Initialize the various sums to 0
1311  Double_t trackPtSumTPC = 0;
1312  Double_t trackPtSumEMCal = 0;
1313  Double_t trackPtSumEMCalRC = 0;
1314  Double_t clusESumEMCal = 0;
1315  Double_t clusESumEMCalRC = 0;
1316 
1317  // Define a 2D vector (initialized to 0) to store the sum of track pT, and another for cluster ET
1318  std::vector<std::vector<Double_t>> trackPtSumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1319  std::vector<std::vector<Double_t>> clusESumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1320 
1321  // Loop over tracks. Sum the track pT:
1322  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal random cone,
1323  // (4) in the DCalRegion in a random cone at each eta-phi
1324  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1325  AliTLorentzVector track;
1326  Double_t trackEta;
1327  Double_t trackPhi;
1328  Double_t trackPt;
1329  Double_t deltaR;
1330  for (auto trackIterator : trackCont->accepted_momentum() ) {
1331 
1332  track.Clear();
1333  track = trackIterator.first;
1334  trackEta = track.Eta();
1335  trackPhi = track.Phi_0_2pi();
1336  trackPt = track.Pt();
1337 
1338  // (1)
1339  if (TMath::Abs(trackEta) < etaTPC) {
1340  trackPtSumTPC += trackPt;
1341  }
1342 
1343  // (2)
1344  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
1345  trackPtSumEMCal += trackPt;
1346  }
1347 
1348  // (3)
1349  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
1350  if (deltaR < jetR) {
1351  trackPtSumEMCalRC += trackPt;
1352  }
1353 
1354  // (4)
1355  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1356  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1357  deltaR = GetDeltaR(&track, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1358  if (deltaR < jetR) {
1359  trackPtSumDCalRC[etaBin][phiBin] += trackPt;
1360  }
1361  }
1362  }
1363 
1364  }
1365 
1366  // Loop over clusters. Sum the cluster ET:
1367  // (1) in the EMCal, (2) in the EMCal random cone, (3) in the DCalRegion in a random cone at each eta-phi
1369  AliTLorentzVector clus;
1370  Double_t clusEta;
1371  Double_t clusPhi;
1372  Double_t clusE;
1373  for (auto clusIterator : clusCont->accepted_momentum() ) {
1374 
1375  clus.Clear();
1376  clus = clusIterator.first;
1377  clusEta = clus.Eta();
1378  clusPhi = clus.Phi_0_2pi();
1379  clusE = clus.E();
1380 
1381  // (1)
1382  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
1383  clusESumEMCal += clusE;
1384  }
1385 
1386  // (2)
1387  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
1388  if (deltaR < jetR) {
1389  clusESumEMCalRC += clusE;
1390  }
1391 
1392  // (3)
1393  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1394  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1395  deltaR = GetDeltaR(&clus, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1396  if (deltaR < jetR) {
1397  clusESumDCalRC[etaBin][phiBin] += clusE;
1398  }
1399  }
1400  }
1401 
1402  }
1403 
1404  // Compute the scale factor for EMCal, as a function of centrality
1405  Double_t numerator = (trackPtSumEMCal + clusESumEMCal) / accEMCal;
1406  Double_t denominator = trackPtSumTPC / accTPC;
1407  Double_t scaleFactor = numerator / denominator;
1408  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
1409  fHistManager.FillTH2(histname, fCent, scaleFactor);
1410 
1411  // Compute the scale factor for DCalRegion in each eta-phi bin, as a function of centrality
1412  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1413  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1414  numerator = (trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin]) / accRC;
1415  scaleFactor = numerator / denominator;
1416  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jetCont->GetArrayName().Data());
1417  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, scaleFactor};
1418  fHistManager.FillTHnSparse(histname, x);
1419  }
1420  }
1421 
1422  // Compute delta pT for EMCal, as a function of centrality
1423  Double_t rho = jetCont->GetRhoVal();
1424  Double_t deltaPt = trackPtSumEMCalRC + clusESumEMCalRC - rho * TMath::Pi() * jetR * jetR;
1425  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
1426  fHistManager.FillTH2(histname, fCent, deltaPt);
1427 
1428  // Compute delta pT for DCalRegion in each eta-phi bin, as a function of centrality
1429  Double_t sf;
1430  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1431  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1433  sf = fBackgroundScalingWeights->GetBinContent(etaDCalRC[etaBin], phiDCalRC[phiBin]);
1434  rho = sf * rho;
1435  }
1436  deltaPt = trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin] - rho * accRC;
1437  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jetCont->GetArrayName().Data());
1438  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, deltaPt};
1439  fHistManager.FillTHnSparse(histname, x);
1440  }
1441  }
1442 
1443  delete r;
1444 
1445 }
1446 
1451 {
1452 
1453  Double_t rho = jetCont->GetRhoVal();
1454 
1455  // Get eta-phi dependent scale factor
1457  Double_t sf = fBackgroundScalingWeights->GetBinContent(jet->Eta(), jet->Phi_0_2pi());
1458  rho = sf * rho;
1459  }
1460 
1461  // Compute pT
1462  Double_t pT = jet->Pt() - rho * jet->Area();
1463 
1464  // If hard-core jet, don't subtract background
1465  TString jetContName = jetCont->GetName();
1466  if (jetContName.Contains("HardCore")) pT = jet->Pt();
1467 
1468  return pT;
1469 }
1470 
1475 {
1476  AliEmcalJet* leadingJet = 0;
1477 
1478  if (jetCont->GetRhoParameter()) {
1479  for(auto jetCand : jetCont->accepted()) {
1480  if (!jetCand) continue;
1481  Double_t jetCandPt = GetJetPt(jetCont, jetCand);
1482  if (leadingJet) {
1483  Double_t leadingJetPt = GetJetPt(jetCont, leadingJet);
1484  if ( jetCandPt < leadingJetPt ) continue;
1485  }
1486  leadingJet = jetCand;
1487  }
1488  }
1489  else {
1490  leadingJet = jetCont->GetLeadingJet();
1491  }
1492 
1493  return leadingJet;
1494 }
1495 
1500 {
1501  Double_t deltaPhi = TMath::Abs(jet1->Phi() - jet2->Phi());
1502  Double_t deltaEta = TMath::Abs(jet1->Eta() - jet2->Eta());
1503  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1504  return deltaR;
1505 }
1506 
1511 {
1512  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
1513  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
1514  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1515  return deltaR;
1516 }
1517 
1525 {
1528 
1529  return kTRUE;
1530 }
1531 
1537 {
1538  TString histname;
1539  AliJetContainer* jets = 0;
1540  TIter nextJetColl(&fJetCollArray);
1541  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1542  TString jetContName = jets->GetName();
1543  if (jetContName.Contains("HardCore")) continue;
1544  Double_t rhoVal = 0;
1545  if (jets->GetRhoParameter()) {
1546  rhoVal = jets->GetRhoVal();
1547  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
1548  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1549  }
1550 
1551  for (auto jet : jets->all()) {
1552 
1553  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1554  Float_t corrPt = GetJetPt(jets, jet);
1555 
1556  // A vs. pT (fill before area cut)
1557  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
1558  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
1559 
1560 
1561  // Rejection reason
1562  UInt_t rejectionReason = 0;
1563  if (!jets->AcceptJet(jet, rejectionReason)) {
1564  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1565  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1566  continue;
1567  }
1568 
1569  // Centrality vs. pT
1570  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
1571  fHistManager.FillTH2(histname.Data(), corrPt, fCent);
1572 
1573  // NEF vs. pT vs. eta vs. phi
1574  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
1575  Double_t x[4] = {jet->Eta(), jet->Phi_0_2pi(), corrPt, jet->NEF()};
1576  fHistManager.FillTHnSparse(histname, x);
1577 
1578  // pT un-downscaled
1579  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
1580  fHistManager.FillTH1(histname.Data(), corrPt, fMBUpscaleFactor);
1581 
1582  // pT-leading vs. pT
1583  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1584  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1585 
1586  // z-leading (charged) vs. pT
1587  TLorentzVector leadPart;
1588  jets->GetLeadingHadronMomentum(leadPart, jet);
1589  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1590  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
1591  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
1592  fHistManager.FillTH2(histname.Data(), corrPt, z);
1593 
1594  // z (charged) vs. pT
1595  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
1596  AliVTrack* track;
1597  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1598  track = static_cast<AliVTrack*>(jet->Track(i));
1599  z = track->Pt() / TMath::Abs(corrPt);
1600  fHistManager.FillTH2(histname.Data(), corrPt, z);
1601  }
1602 
1603  // Nconst vs. pT
1604  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
1605  fHistManager.FillTH2(histname.Data(), corrPt, jet->GetNumberOfConstituents());
1606 
1607  } //jet loop
1608 
1609  //---------------------------------------------------------------------------
1610  // Do study of background (if requested)
1612  }
1613 }
1614 
1619 {
1620  Double_t contents[30]={0};
1621  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1622  if (!histJetObservables) return;
1623  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1624  TString title(histJetObservables->GetAxis(n)->GetTitle());
1625  if (title=="Centrality (%)")
1626  contents[n] = fCent;
1627  else if (title=="LeadingHadronRequired")
1628  contents[n] = fDijet.leadingHadronCutType;
1629  else if (title=="#it{p}_{T,trig jet} (GeV/#it{c})")
1630  contents[n] = fDijet.trigJetPt;
1631  else if (title=="#it{p}_{T,ass jet} (GeV/#it{c})")
1632  contents[n] = fDijet.assJetPt;
1633  else if (title=="#phi_{trig jet}")
1634  contents[n] = fDijet.trigJetPhi;
1635  else if (title=="#phi_{ass jet}")
1636  contents[n] = fDijet.assJetPhi;
1637  else if (title=="#eta_{trig jet}")
1638  contents[n] = fDijet.trigJetEta;
1639  else if (title=="#eta_{ass jet}")
1640  contents[n] = fDijet.assJetEta;
1641  else
1642  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1643  }
1644  histJetObservables->Fill(contents);
1645 }
1646 
1651 {
1652  Double_t contents[30]={0};
1653  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1654  if (!histJetObservables) return;
1655  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1656  TString title(histJetObservables->GetAxis(n)->GetTitle());
1657  if (title=="Centrality (%)")
1658  contents[n] = fCent;
1659  else if (title=="#Delta#phi")
1660  contents[n] = fDijet.deltaPhi;
1661  else if (title=="#Delta#eta")
1662  contents[n] = fDijet.deltaEta;
1663  else if (title=="A_{J}")
1664  contents[n] = fDijet.AJ;
1665  else if (title=="x_{J}")
1666  contents[n] = fDijet.xJ;
1667  else if (title=="k_{Ty} (GeV)")
1668  contents[n] = fDijet.kTy;
1669  else if (title=="N_{tracks, trig jet}")
1670  contents[n] = fDijet.trigJet->GetNumberOfTracks();
1671  else if (title=="N_{tracks, ass jet}")
1672  contents[n] = fDijet.assJet->GetNumberOfTracks();
1673  else
1674  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1675  }
1676  histJetObservables->Fill(contents);
1677 }
1678 
1683 {
1684  Double_t contents[30]={0};
1685  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1686  if (!histJetObservables) return;
1687  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1688  TString title(histJetObservables->GetAxis(n)->GetTitle());
1689  if (title=="A_{J}")
1690  contents[n] = fDijet.AJ;
1691  else if (title=="#Delta#phi")
1692  contents[n] = deltaPhi;
1693  else if (title=="#it{p}_{T,particle} (GeV/#it{c})")
1694  contents[n] = trackPt;
1695  else if (title=="#it{p}_{T}#parallel (GeV/#it{c})")
1696  contents[n] = balancePt;
1697  else
1698  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1699  }
1700  histJetObservables->Fill(contents);
1701 }
1702 
1707 {
1708  // Matching efficiency histogram
1709  TString histname = "GeometricalMatchingEfficiency";
1710 
1711  // If we have a matching di-jet, fill the geometrical matching histogram
1712  if (fMatchingDijet.assJet) {
1713  fHistManager.FillTH1(histname.Data(), 1);
1714 
1717  Bool_t isSwitched = trigDeltaR > fMatchingJetR;
1718 
1719  TString thnname = "GeometricalMatching";
1720  Double_t contents[30]={0};
1721  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(thnname.Data()));
1722  if (!histJetObservables) return;
1723  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1724  TString title(histJetObservables->GetAxis(n)->GetTitle());
1725  if (title=="Centrality (%)")
1726  contents[n] = fCent;
1727  else if (title=="isSwitched")
1728  contents[n] = isSwitched;
1729  else if (title=="#DeltaR_{trig}")
1730  contents[n] = trigDeltaR;
1731  else if (title=="#DeltaR_{ass}")
1732  contents[n] = assDeltaR;
1733  else if (title=="trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1734  contents[n] = fMatchingDijet.trigJetPt - fDijet.trigJetPt;
1735  else if (title=="ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1736  contents[n] = fMatchingDijet.assJetPt - fDijet.assJetPt;
1737  else if (title=="A_{J} low-threshold")
1738  contents[n] = fMatchingDijet.AJ;
1739  else if (title=="A_{J} hard-core")
1740  contents[n] = fDijet.AJ;
1741  else
1742  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1743  }
1744  histJetObservables->Fill(contents);
1745  }
1746  else {
1747  fHistManager.FillTH1(histname.Data(), 0.);
1748  }
1749 }
1750 
1751 /*
1752  * This function fills the histograms for the calorimeter performance study.
1753  */
1755 {
1756  // Define some vars
1757  TString histname;
1758  Double_t Enonlin;
1759  Double_t Ehadcorr;
1760  Double_t Ecore;
1761  Int_t absId;
1762  Double_t ecell;
1763  Double_t leadEcell;
1764 
1765  // Get cells from event
1766  fCaloCells = InputEvent()->GetEMCALCells();
1767  AliVCaloCells* phosCaloCells = InputEvent()->GetPHOSCells();
1768 
1769  // Loop through clusters and plot cluster THnSparse (centrality, cluster type, E, E-hadcorr, has matched track, M02, Ncells)
1770  AliClusterContainer* clusters = 0;
1771  TIter nextClusColl(&fClusterCollArray);
1772  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
1774  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1775 
1776  // Determine cluster type (EMCal/DCal/Phos)
1777  ClusterType clusType = kNA;
1778  if (it->second->IsEMCAL()) {
1779  Double_t phi = it->first.Phi_0_2pi();
1780  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
1781  if (isDcal == 0) {
1782  clusType = kEMCal;
1783  } else if (isDcal == 1) {
1784  clusType = kDCal;
1785  }
1786  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
1787  clusType = kPHOS;
1788  }
1789 
1790  // rejection reason plots, to make efficiency correction
1791  if (it->second->IsEMCAL()) {
1792  histname = TString::Format("%s/hClusterRejectionReasonEMCal", clusters->GetArrayName().Data());
1793  UInt_t rejectionReason = 0;
1794  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
1795  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
1796  continue;
1797  }
1798  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
1799  histname = TString::Format("%s/hClusterRejectionReasonPHOS", clusters->GetArrayName().Data());
1800  UInt_t rejectionReason = 0;
1801  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
1802  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
1803  continue;
1804  }
1805  } else {
1806  continue;
1807  }
1808 
1809  // Fill cluster spectra by SM, and fill cell histograms
1810  Enonlin = 0;
1811  Ehadcorr = 0;
1812  Ecore = 0;
1813  if (it->second->IsEMCAL()) {
1814 
1815  Ehadcorr = it->second->GetHadCorrEnergy();
1816  Enonlin = it->second->GetNonLinCorrEnergy();
1817  Ecore = Ehadcorr;
1818 
1819  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
1820  if (sm >=0 && sm < 20) {
1821  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1822  fHistManager.FillTH1(histname, it->second->E());
1823  }
1824  else {
1825  AliError(Form("Supermodule %d does not exist!", sm));
1826  }
1827 
1828  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
1829  histname = TString::Format("Cells/hCellEnergyAccepted");
1830  leadEcell = 0;
1831  for (Int_t iCell = 0; iCell < it->second->GetNCells(); iCell++){
1832  absId = it->second->GetCellAbsId(iCell);
1833  ecell = fCaloCells->GetCellAmplitude(absId);
1834  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
1835  if (ecell > leadEcell) {
1836  leadEcell = ecell;
1837  }
1838  }
1839  // Plot also the leading cell
1840  histname = TString::Format("Cells/hCellEnergyLeading");
1841  fHistManager.FillTH3(histname, leadEcell, fCent, kEMCal);
1842 
1843  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
1844 
1845  Ehadcorr = it->second->E();
1846  Enonlin = it->second->E();
1847  Ecore = it->second->GetCoreEnergy();
1848 
1849  Int_t relid[4];
1850  if (fPHOSGeo) {
1851  fPHOSGeo->AbsToRelNumbering(it->second->GetCellAbsId(0), relid);
1852  Int_t sm = relid[0];
1853  if (sm >=1 && sm < 5) {
1854  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1855  fHistManager.FillTH1(histname, it->second->E());
1856  }
1857  else {
1858  AliError(Form("Supermodule %d does not exist!", sm));
1859  }
1860  }
1861 
1862  // Get cells from each accepted cluster, and plot centrality vs. cell energy vs. cell type
1863  histname = TString::Format("Cells/hCellEnergyAccepted");
1864  leadEcell = 0;
1865  for (Int_t iCell = 0; iCell < it->second->GetNCells(); iCell++){
1866  absId = it->second->GetCellAbsId(iCell);
1867  ecell = phosCaloCells->GetCellAmplitude(absId);
1868  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
1869  if (ecell > leadEcell) {
1870  leadEcell = ecell;
1871  }
1872  }
1873  // Plot also the leading cell
1874  histname = TString::Format("Cells/hCellEnergyLeading");
1875  fHistManager.FillTH3(histname, leadEcell, fCent, kPHOS);
1876  }
1877 
1878  // Check if the cluster has a matched track
1879  Int_t hasMatchedTrack = -1;
1880  Int_t nMatchedTracks = it->second->GetNTracksMatched();
1881  if (nMatchedTracks == 0) {
1882  hasMatchedTrack = 0;
1883  } else if (nMatchedTracks > 0) {
1884  hasMatchedTrack = 1;
1885  }
1886 
1887  Double_t contents[30]={0};
1888  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
1889  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1890  if (!histClusterObservables) return;
1891  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1892  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1893  if (title=="Centrality %")
1894  contents[i] = fCent;
1895  else if (title=="#eta")
1896  contents[i] = it->first.Eta();
1897  else if (title=="#phi")
1898  contents[i] = it->first.Phi_0_2pi();
1899  else if (title=="#it{E}_{clus} (GeV)")
1900  contents[i] = Enonlin;
1901  else if (title=="#it{E}_{clus, hadcorr} (GeV)")
1902  contents[i] = Ehadcorr;
1903  else if (title=="#it{E}_{core} (GeV)")
1904  contents[i] = Ecore;
1905  else if (title=="Matched track")
1906  contents[i] = hasMatchedTrack;
1907  else if (title=="M02")
1908  contents[i] = it->second->GetM02();
1909  else if (title=="Ncells")
1910  contents[i] = it->second->GetNCells();
1911  else
1912  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1913  }
1914  histClusterObservables->Fill(contents);
1915 
1916  }
1917 
1918  }
1919 
1920  // plot neutral jets THnSparse (centrality, eta, phi, E, Nclusters)
1921  AliJetContainer* jets = 0;
1922  TIter nextJetColl(&fJetCollArray);
1923  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1924 
1925  for (auto jet : jets->accepted()) {
1926 
1927  Double_t contents[30]={0};
1928  histname = TString::Format("%s/hClusterJetObservables", jets->GetArrayName().Data());
1929  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1930  if (!histJetObservables) return;
1931  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
1932  TString title(histJetObservables->GetAxis(i)->GetTitle());
1933  if (title=="Centrality (%)")
1934  contents[i] = fCent;
1935  else if (title=="#eta_{jet}")
1936  contents[i] = jet->Eta();
1937  else if (title=="#phi_{jet} (rad)")
1938  contents[i] = jet->Phi_0_2pi();
1939  else if (title=="#it{E}_{T} (GeV)")
1940  contents[i] = jet->Pt();
1941  else if (title=="#rho (GeV/#it{c})")
1942  contents[i] = jet->Pt() / jet->Area();
1943  else if (title=="N_{clusters}")
1944  contents[i] = jet->GetNumberOfClusters();
1945  else
1946  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1947  }
1948  histJetObservables->Fill(contents);
1949 
1950  }
1951 
1952  }
1953 
1954  // Loop through all cells and fill histos
1955  Int_t sm;
1956  Int_t relid[4];
1957  Double_t patchSumEMCal[20] = {0.};
1958  Double_t patchSumPHOS[4] = {0.};
1959  for (Int_t i=0; i<fCaloCells->GetNumberOfCells(); i++) {
1960 
1961  absId = fCaloCells->GetCellNumber(i);
1962  ecell = fCaloCells->GetCellAmplitude(absId);
1963 
1964  // Fill cell histo
1965  histname = TString::Format("Cells/hCellEnergyAll");
1966  fHistManager.FillTH3(histname, ecell, fCent, kEMCal); // Note: I don't distinguish EMCal from DCal cells
1967 
1968  // Fill cell patch histo, per SM
1969  sm = fGeom->GetSuperModuleNumber(absId);
1970  if (sm >=0 && sm < 20) {
1971  patchSumEMCal[sm] += ecell;
1972  }
1973 
1974  }
1975 
1976  for (Int_t i=0; i<phosCaloCells->GetNumberOfCells(); i++) {
1977 
1978  absId = phosCaloCells->GetCellNumber(i);
1979  ecell = phosCaloCells->GetCellAmplitude(absId);
1980 
1981  // Fill cell histo
1982  histname = TString::Format("Cells/hCellEnergyAll");
1983  fHistManager.FillTH3(histname, ecell, fCent, kPHOS);
1984 
1985  // Fill cell patch histo, per SM
1986  fPHOSGeo->AbsToRelNumbering(absId, relid);
1987  sm = relid[0];
1988  if (sm >=1 && sm < 5) {
1989  patchSumPHOS[sm-1] += ecell;
1990  }
1991 
1992  }
1993 
1994  for (Int_t sm = 0; sm < 20; sm++) {
1995  histname = TString::Format("Cells/BySM/hEmcalPatchEnergy_SM%d", sm);
1996  fHistManager.FillTH2(histname, patchSumEMCal[sm], fCent);
1997  }
1998 
1999  for (Int_t sm = 1; sm < 5; sm++) {
2000  histname = TString::Format("Cells/BySM/hPhosPatchEnergy_SM%d", sm);
2001  fHistManager.FillTH2(histname, patchSumPHOS[sm-1], fCent);
2002  }
2003 
2004 }
2005 
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:26
AliJetContainer * GetJetContainer(Int_t i=0) const
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.
Double_t fDijetLeadingHadronPt
leading hadron pT threshold for leading jet in dijet
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
Double_t fMinAssJetPt
subleading jet min pT in a dijet pair, for it to be accepted
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.
AliEventCuts fEventCuts
event selection utility
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
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="")
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
Double_t Phi_0_2pi() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
BeamType fForceBeamType
forced beam type
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Definition: External.C:212
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
Double_t GetJetPt(AliJetContainer *jetCont, AliEmcalJet *jet)
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
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.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Handler for downscale factors for various triggers obtained from the OCDB.
TFile * file
TList with histograms for a given trigger.
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
const AliTrackIterableMomentumContainer accepted_momentum() const
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
Dijet_t fMatchingDijet
! low-threshold matching dijet, for matching study
void SetRejectionReasonLabels(TAxis *axis)
void UserCreateOutputObjects()
Main initialization function on the worker.
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Double_t fTrackConstituentThreshold
constituent threshold for matching study
TH1D * fBackgroundScalingWeights
Histogram storing eta-phi weights for full-jet background scale factors.
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 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.
void LoadBackgroundScalingHistogram(const char *path="alien:///alice/cern.ch/user/j/jmulliga/BackgroundScalingWeights.root", const char *name="hBackgroundScalingWeights")
Bool_t fDoCaloStudy
Set whether to perform calorimeter detector study.
const AliJetIterableContainer all() const
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal