AliPhysics  ec7afe5 (ec7afe5)
 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 <TClonesArray.h>
17 #include <TH1F.h>
18 #include <TH2F.h>
19 #include <TList.h>
20 #include <THnSparse.h>
21 #include <TRandom3.h>
22 
23 #include <AliVCluster.h>
24 #include <AliVParticle.h>
25 #include <AliLog.h>
26 
27 #include "AliTLorentzVector.h"
28 #include "AliEmcalJet.h"
29 #include "AliRhoParameter.h"
30 #include "AliJetContainer.h"
31 #include "AliParticleContainer.h"
32 #include "AliClusterContainer.h"
33 #include "AliEMCALGeometry.h"
35 
37 
41 
47  fHistManager(),
48  fEventCuts(0),
49  fEventCutList(0),
50  fUseManualEventCuts(kFALSE),
51  fDeltaPhiMin(0),
52  fMinTrigJetPt(0),
53  fMinAssJetPt(0),
54  fDijetLeadingHadronPt(0),
55  fMaxPt(250),
56  fNCentHistBins(0),
57  fCentHistBins(0),
58  fPlotJetHistograms(kFALSE),
59  fPlotDijetCandHistograms(kFALSE),
60  fPlotDijetImbalanceHistograms(kFALSE),
61  fComputeBackground(kFALSE),
62  fDoMomentumBalance(kFALSE),
63  fDoGeometricalMatching(kFALSE),
64  fMatchingJetR(0.2),
65  fTrackConstituentThreshold(0),
66  fClusterConstituentThreshold(0),
67  fMBUpscaleFactor(1.)
68 {
72 }
73 
80  AliAnalysisTaskEmcalJet(name, kTRUE),
81  fHistManager(name),
82  fEventCuts(0),
83  fEventCutList(0),
84  fUseManualEventCuts(kFALSE),
85  fDeltaPhiMin(0),
86  fMinTrigJetPt(0),
87  fMinAssJetPt(0),
88  fDijetLeadingHadronPt(0),
89  fMaxPt(250),
90  fNCentHistBins(0),
91  fCentHistBins(0),
92  fPlotJetHistograms(kFALSE),
93  fPlotDijetCandHistograms(kFALSE),
94  fPlotDijetImbalanceHistograms(kFALSE),
95  fComputeBackground(kFALSE),
96  fDoMomentumBalance(kFALSE),
97  fDoGeometricalMatching(kFALSE),
98  fMatchingJetR(0.2),
99  fTrackConstituentThreshold(0),
100  fClusterConstituentThreshold(0),
101  fMBUpscaleFactor(1.)
102 {
104  Dijet_t fDijet;
106 }
107 
112 {
113 }
114 
119 {
120  fNCentHistBins = 4;
122  fCentHistBins[0] = 0;
123  fCentHistBins[1] = 10;
124  fCentHistBins[2] = 30;
125  fCentHistBins[3] = 50;
126  fCentHistBins[4] = 90;
127 }
128 
134 {
136 
142 
143  TIter next(fHistManager.GetListOfHistograms());
144  TObject* obj = 0;
145  while ((obj = next())) {
146  fOutput->Add(obj);
147  }
148 
149  // Intialize AliEventCuts
150  fEventCutList = new TList();
151  fEventCutList ->SetOwner();
152  fEventCutList ->SetName("EventCutOutput");
153 
154  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
155  if(fUseManualEventCuts==1)
156  {
157  fEventCuts.SetManualMode();
158  // Configure manual settings here
159  // ...
160  }
161  fEventCuts.AddQAplotsToList(fEventCutList);
162  fOutput->Add(fEventCutList);
163 
164  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
165 }
166 
167 /*
168  * This function allocates the histograms for single jets.
169  * A set of histograms is allocated per each jet container.
170  */
172 {
173  TString histname;
174  TString title;
175 
176  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
177  Int_t minPtBin = -fMaxPt/2 + 25;
178  Int_t maxPtBin = fMaxPt/2 + 25;
179 
180  AliJetContainer* jets = 0;
181  TIter nextJetColl(&fJetCollArray);
182  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
183 
184  // Jet rejection reason
185  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
186  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
187  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
188  SetRejectionReasonLabels(hist->GetXaxis());
189 
190  // Rho vs. Centrality
191  if (!jets->GetRhoName().IsNull()) {
192  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
193  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
194  fHistManager.CreateTH2(histname.Data(), title.Data(), 101, 0, 101, 100, 0, 500);
195  }
196 
197  // Centrality vs. pT
198  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
199  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});Centrality (%);counts";
200  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, minPtBin, maxPtBin, 10, 0, 100);
201 
202  // pT vs. eta vs. phi
203  histname = TString::Format("%s/JetHistograms/hPtVsEtaVsPhi", jets->GetArrayName().Data());
204  title = histname + ";#eta_{jet} (rad);#phi_{jet} (rad);#it{p}_{T}^{corr} (GeV/#it{c})";
205  fHistManager.CreateTH3(histname.Data(), title.Data(), 50, -0.5, 0.5, 101, 0, TMath::Pi() * 2.02, 75, 0, maxPtBin);
206 
207  // pT upscaled
208  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
209  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c})";
210  fHistManager.CreateTH1(histname.Data(), title.Data(), 75, 0, maxPtBin, "s");
211 
212  // pT-leading vs. pT
213  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
214  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
215  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, minPtBin, maxPtBin, 150, 0, 150);
216 
217  // A vs. pT
218  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
219  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
220  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, minPtBin, maxPtBin, fMaxPt/3, 0, 1.5);
221 
222  // Cent vs. NEF vs. pT
223  histname = TString::Format("%s/JetHistograms/hCentVsNEFVsPt", jets->GetArrayName().Data());
224  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF;Centrality (%)";
225  fHistManager.CreateTH3(histname.Data(), title.Data(), nPtBins, minPtBin, maxPtBin, fMaxPt/5, 0, 1.0, 50, 0, 100);
226 
227  // Cent vs. z-leading (charged) vs. pT
228  histname = TString::Format("%s/JetHistograms/hCentVsZLeadingVsPt", jets->GetArrayName().Data());
229  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading};Centrality (%)";
230  fHistManager.CreateTH3(histname.Data(), title.Data(), nPtBins, minPtBin, maxPtBin, fMaxPt/5, 0, 1.0, 50, 0, 100);
231 
232  // Cent vs. z (charged) vs. pT
233  histname = TString::Format("%s/JetHistograms/hCentVsZVsPt", jets->GetArrayName().Data());
234  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z};Centrality (%)";
235  fHistManager.CreateTH3(histname.Data(), title.Data(), nPtBins, minPtBin, maxPtBin, fMaxPt/5, 0, 1.0, 50, 0, 100);
236 
237  // Cent vs. Nconst vs. pT
238  histname = TString::Format("%s/JetHistograms/hCentVsNConstVsPt", jets->GetArrayName().Data());
239  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents;Centrality (%)";
240  fHistManager.CreateTH3(histname.Data(), title.Data(), nPtBins, minPtBin, maxPtBin, fMaxPt/5, 0, fMaxPt, 50, 0, 100);
241 
242  // Allocate background subtraction histograms, if enabled
243  if (fComputeBackground) {
244 
245  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
246  title = histname + ";Centrality;Scale factor;counts";
247  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
248 
249  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
250  title = histname + ";#delta#it{p}_{T} (GeV/#it{c});counts";
251  fHistManager.CreateTH1(histname.Data(), title.Data(), 400, -100, 100);
252 
253  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCalRegion", jets->GetArrayName().Data());
254  title = histname + ";Centrality;Scale factor;Eta bin";
255  fHistManager.CreateTH3(histname.Data(), title.Data(), 50, 0, 100, 200, 0, 10, 28, -0.5, 27.5);
256 
257  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtDCalRegion", jets->GetArrayName().Data());
258  title = histname + ";#delta#it{p}_{T} (GeV/#it{c});Eta bin;counts";
259  fHistManager.CreateTH2(histname.Data(), title.Data(), 400, -100, 100, 20, -0.5, 19.5);
260 
261  }
262 
263  }
264 
265  // MB downscale factor histogram
266  histname = "Trigger/hMBDownscaleFactor";
267  title = histname + ";Downscale factor;counts";
268  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 200);
269 
270 }
271 
272 /*
273  * This function allocates the histograms for dijet candidates, i.e. dijet pairs with the leading jet
274  * passing the trigger condition, but no condition on the subleading jet. In particular this histogram is
275  * designed to study the kinematic selection of dijets.
276  */
278 {
279  // Allocate dijet THnSparse
280  AliJetContainer* jets = 0;
281  TIter nextJetColl(&fJetCollArray);
282  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
283  TString axisTitle[30]= {""};
284  Int_t nbins[30] = {0};
285  Double_t min[30] = {0.};
286  Double_t max[30] = {0.};
287  Double_t *binEdges[20] = {0};
288  Int_t dim = 0;
289 
290  if (fForceBeamType != kpp) {
291  axisTitle[dim] = "Centrality (%)";
292  nbins[dim] = fNCentHistBins;
293  binEdges[dim] = fCentHistBins;
294  min[dim] = fCentHistBins[0];
295  max[dim] = fCentHistBins[fNCentHistBins];
296  dim++;
297  }
298 
299  axisTitle[dim] = "LeadingHadronRequired";
300  nbins[dim] = 2;
301  min[dim] = -0.5;
302  max[dim] = 1.5;
303  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
304  dim++;
305 
306  axisTitle[dim] = "#it{p}_{T,trig jet} (GeV/#it{c})";
307  nbins[dim] = fMaxPt/3;
308  min[dim] = -fMaxPt/2 + 25;
309  max[dim] = fMaxPt/2 + 25;
310  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
311  dim++;
312 
313  axisTitle[dim] = "#it{p}_{T,ass jet} (GeV/#it{c})";
314  nbins[dim] = fMaxPt/3;
315  min[dim] = -fMaxPt/2 + 25;
316  max[dim] = fMaxPt/2 + 25;
317  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
318  dim++;
319 
320  axisTitle[dim] = "#phi_{trig jet}";
321  nbins[dim] = fMaxPt/3;
322  min[dim] = 0;
323  max[dim] = TMath::TwoPi();
324  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
325  dim++;
326 
327  axisTitle[dim] = "#phi_{ass jet}";
328  nbins[dim] = fMaxPt/3;
329  min[dim] = 0;
330  max[dim] = TMath::TwoPi();
331  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
332  dim++;
333 
334  axisTitle[dim] = "#eta_{trig jet}";
335  nbins[dim] = fMaxPt/3;
336  min[dim] = -1;
337  max[dim] = 1;
338  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
339  dim++;
340 
341  axisTitle[dim] = "#eta_{ass jet}";
342  nbins[dim] = fMaxPt/3;
343  min[dim] = -1;
344  max[dim] = 1;
345  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
346  dim++;
347 
348  TString thnname = TString::Format("%s/DijetCandObservables", jets->GetArrayName().Data());
349  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
350  for (Int_t i = 0; i < dim; i++) {
351  hn->GetAxis(i)->SetTitle(axisTitle[i]);
352  hn->SetBinEdges(i, binEdges[i]);
353  }
354  }
355 }
356 
357 /*
358  * This function allocates the histograms for accepted dijets.
359  * The purpose is to study in detail the imbalance properties of the accepted dijets.
360  */
362 {
363  // Allocate dijet imbalance THnSparse
364  AliJetContainer* jets = 0;
365  TIter nextJetColl(&fJetCollArray);
366  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
367  TString axisTitle[30]= {""};
368  Int_t nbins[30] = {0};
369  Double_t min[30] = {0.};
370  Double_t max[30] = {0.};
371  Double_t *binEdges[20] = {0};
372  Int_t dim = 0;
373 
374  if (fForceBeamType != kpp) {
375  axisTitle[dim] = "Centrality (%)";
376  nbins[dim] = fNCentHistBins;
377  binEdges[dim] = fCentHistBins;
378  min[dim] = fCentHistBins[0];
379  max[dim] = fCentHistBins[fNCentHistBins];
380  dim++;
381  }
382 
383  axisTitle[dim] = "#Delta#phi";
384  nbins[dim] = 100;
385  min[dim] = 0;
386  max[dim] = 4;
387  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
388  dim++;
389 
390  axisTitle[dim] = "#Delta#eta";
391  nbins[dim] = 100;
392  min[dim] = -2;
393  max[dim] = 2;
394  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
395  dim++;
396 
397  axisTitle[dim] = "A_{J}";
398  nbins[dim] = 100;
399  min[dim] = 0;
400  max[dim] = 1;
401  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
402  dim++;
403 
404  axisTitle[dim] = "x_{J}";
405  nbins[dim] = 100;
406  min[dim] = 0;
407  max[dim] = 1;
408  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
409  dim++;
410 
411  axisTitle[dim] = "k_{Ty} (GeV)";
412  nbins[dim] = 100;
413  min[dim] = 0;
414  max[dim] = 100;
415  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
416  dim++;
417 
418  axisTitle[dim] = "N_{tracks, trig jet}";
419  nbins[dim] = fMaxPt/5;
420  min[dim] = 0;
421  max[dim] = 100;
422  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
423  dim++;
424 
425  axisTitle[dim] = "N_{tracks, ass jet}";
426  nbins[dim] = fMaxPt/5;
427  min[dim] = 0;
428  max[dim] = 100;
429  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
430  dim++;
431 
432  TString thnname = TString::Format("%s/DijetImbalanceObservables", jets->GetArrayName().Data());
433  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
434  for (Int_t i = 0; i < dim; i++) {
435  hn->GetAxis(i)->SetTitle(axisTitle[i]);
436  hn->SetBinEdges(i, binEdges[i]);
437  }
438  }
439 }
440 
441 /*
442  * This function allocates the histograms for the momentum balance study.
443  */
445 {
446  AliJetContainer* jets = 0;
447  TIter nextJetColl(&fJetCollArray);
448  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
449 
450  // Allocate THnSparse
451  TString axisTitle[30]= {""};
452  Int_t nbins[30] = {0};
453  Double_t min[30] = {0.};
454  Double_t max[30] = {0.};
455  Double_t *binEdges[20] = {0};
456  Int_t dim = 0;
457 
458  axisTitle[dim] = "A_{J}";
459  nbins[dim] = 100;
460  min[dim] = 0;
461  max[dim] = 1;
462  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
463  dim++;
464 
465  axisTitle[dim] = "#Delta#phi";
466  nbins[dim] = 100;
467  min[dim] = -4;
468  max[dim] = 4;
469  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
470  dim++;
471 
472  axisTitle[dim] = "#it{p}_{T,particle} (GeV/#it{c})";
473  nbins[dim] = 9;
474  Double_t* pTParticleBins = new Double_t[nbins[dim]+1];
475  GenerateFixedBinArray(1, 0.15, 0.3, pTParticleBins);
476  GenerateFixedBinArray(1, 0.3, 0.5, pTParticleBins+1);
477  GenerateFixedBinArray(1, 0.5, 1, pTParticleBins+2);
478  GenerateFixedBinArray(2, 1, 5, pTParticleBins+3);
479  GenerateFixedBinArray(3, 5, 20, pTParticleBins+5);
480  GenerateFixedBinArray(1, 20, 150, pTParticleBins+8);
481  min[dim] = 0;
482  max[dim] = pTParticleBins[nbins[dim]];
483  binEdges[dim] = pTParticleBins;
484  dim++;
485 
486  axisTitle[dim] = "#it{p}_{T}#parallel (GeV/#it{c})";
487  nbins[dim] = 80;
488  Double_t* pTParallelBins = new Double_t[nbins[dim]+1];
489  GenerateFixedBinArray(20, 0, 2, pTParallelBins);
490  GenerateFixedBinArray(16, 2, 10, pTParallelBins+20);
491  GenerateFixedBinArray(10, 10, 20, pTParallelBins+36);
492  GenerateFixedBinArray(10, 20, 40, pTParallelBins+46);
493  GenerateFixedBinArray(24, 40, 150, pTParallelBins+56);
494  min[dim] = 0;
495  max[dim] = pTParallelBins[nbins[dim]];
496  binEdges[dim] = pTParallelBins;
497  dim++;
498 
499  TString thnname = TString::Format("%s/MomentumBalance", jets->GetArrayName().Data());
500  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
501  for (Int_t i = 0; i < dim; i++) {
502  hn->GetAxis(i)->SetTitle(axisTitle[i]);
503  hn->SetBinEdges(i, binEdges[i]);
504  }
505  }
506 }
507 
508 /*
509  * This function allocates the histograms for the constituent dijet study.
510  */
512 {
513  // Allocate geometrical matching THnSparse
514  TString axisTitle[30]= {""};
515  Int_t nbins[30] = {0};
516  Double_t min[30] = {0.};
517  Double_t max[30] = {0.};
518  Double_t *binEdges[20] = {0};
519  Int_t dim = 0;
520 
521  if (fForceBeamType != kpp) {
522  axisTitle[dim] = "Centrality (%)";
523  nbins[dim] = fNCentHistBins;
524  binEdges[dim] = fCentHistBins;
525  min[dim] = fCentHistBins[0];
526  max[dim] = fCentHistBins[fNCentHistBins];
527  dim++;
528  }
529 
530  axisTitle[dim] = "isSwitched";
531  nbins[dim] = 2;
532  min[dim] = -0.5;
533  max[dim] = 1.5;
534  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
535  dim++;
536 
537  axisTitle[dim] = "#DeltaR_{trig}";
538  nbins[dim] = 50;
539  min[dim] = 0;
540  max[dim] = 0.5;
541  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
542  dim++;
543 
544  axisTitle[dim] = "#DeltaR_{ass}";
545  nbins[dim] = 50;
546  min[dim] = 0;
547  max[dim] = 0.5;
548  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
549  dim++;
550 
551  axisTitle[dim] = "trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
552  nbins[dim] = 100;
553  min[dim] = -50;
554  max[dim] = 50;
555  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
556  dim++;
557 
558  axisTitle[dim] = "ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
559  nbins[dim] = 100;
560  min[dim] = -50;
561  max[dim] = 50;
562  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
563  dim++;
564 
565  axisTitle[dim] = "A_{J} low-threshold";
566  nbins[dim] = 100;
567  min[dim] = 0;
568  max[dim] = 1;
569  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
570  dim++;
571 
572  axisTitle[dim] = "A_{J} hard-core";
573  nbins[dim] = 100;
574  min[dim] = 0;
575  max[dim] = 1;
576  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
577  dim++;
578 
579  TString thnname = "GeometricalMatching";
580  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
581  for (Int_t i = 0; i < dim; i++) {
582  hn->GetAxis(i)->SetTitle(axisTitle[i]);
583  hn->SetBinEdges(i, binEdges[i]);
584  }
585 
586  // Allocate other histograms
587  TString histname;
588  TString title;
589  histname = "GeometricalMatchingEfficiency";
590  title = histname + ";isMatched;counts";
591  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 2, -0.5, 1.5);
592 }
593 
599 {
601 
602  fNeedEmcalGeom = kTRUE;
603 
604  AliInfo(Form("Trigger jet threshold = %f, Associated jet threshold = %f", fMinTrigJetPt, fMinAssJetPt));
605  AliInfo(Form("Leading hadron threshold (for dijet leading jet): %f GeV", fDijetLeadingHadronPt));
606  AliInfo(Form("Momentum balance study: %d", fDoMomentumBalance));
607  AliInfo(Form("Geometrical matching study: %d", fDoGeometricalMatching));
608 
609 }
610 
615 
616  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
617 
618  // Get instance of the downscale factor helper class
620  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
621 
622  // There are two possible min bias triggers for LHC15o
623  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
624  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
625  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
626 
627  // Get the downscale factor for whichever MB trigger exists in the given run
628  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
629  Double_t downscalefactor;
630  for (auto i : runtriggers) {
631  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
632  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
633  break;
634  }
635  }
636 
637  // Store the inverse of the downscale factor, used later to weight the pT spectrum
638  fMBUpscaleFactor = 1/downscalefactor;
639 
640  TString histname = "Trigger/hMBDownscaleFactor";
641  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
642 
643 }
644 
649 {
650  if (!fEventCuts.AcceptEvent(InputEvent()))
651  {
652  PostData(1, fOutput);
653  return kFALSE;
654  }
655  return kTRUE;
656 }
657 
666 {
667  TString histname;
668  AliJetContainer* jetCont = 0;
669  TIter next(&fJetCollArray);
670  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
671  TString jetContName = jetCont->GetName();
672  if (jetContName.Contains("HardCore")) continue;
673 
674  //-----------------------------------------------------------------------------
675  // Find the leading di-jet candidate in each event, and if it satisfies the
676  // trig jet pT threshold, then fill di-jet candidate histogram (regardless of ass jet).
677  // The idea is to study the kinematic selections in post-processing.
678 
679  // Loop over leading hadron cut or not
680  for (Int_t leadingHadronCutType=0; leadingHadronCutType<2; leadingHadronCutType++) {
681 
682  // Find the dijet candidate of the event and store its info in struct fDijet
683  FindDijet(jetCont, leadingHadronCutType);
684 
685  // If we find a dijet candidate (i.e. acceptable trig jet; ass jet accepted or not), fill the di-jet candidate histogram
687  TString histname = TString::Format("%s/DijetCandObservables", jetCont->GetArrayName().Data());
688  FillDijetCandHistograms(histname);
689  }
690 
691  }
692 
693  //---------------------------------------------------------------------------------------------------
694  // Now, study the accepted dijet selection -- specified by the trig/ass jet pT conditions
695 
696  // Find the dijet candidate of the event and store its info in struct fDijet
697  FindDijet(jetCont, 0);
698 
699  // If we find an accepted dijet, fill the dijet imbalance histogram
701  TString histname = TString::Format("%s/DijetImbalanceObservables", jetCont->GetArrayName().Data());
703  }
704 
705  // If we find an accepted dijet, perform momentum-balance study (if requested)
707  histname = TString::Format("%s/MomentumBalance", jetCont->GetArrayName().Data());
708  DoMomentumBalance(histname);
709  }
710 
711  }
712 
713  //---------------------------------------------------------------------------
714  // Do the constituent threshold and geometrical matching study (if requested)
717 
718  return kTRUE;
719 }
720 
729 {
730  fDijet.clear();
731  fDijet.leadingHadronCutType = leadingHadronCutBin;
732 
733  // Get trigger jet
734  AliEmcalJet* trigJet = 0;
735  if (jetCont->GetRhoParameter())
736  trigJet = jetCont->GetLeadingJet("rho");
737  else
738  trigJet = jetCont->GetLeadingJet();
739  if(!trigJet) return;
740 
741  // Skip the event if the leading jet doesn't satisfy the pT threshold
742  Double_t trigJetPt = GetJetPt(jetCont, trigJet);
743  if ( trigJetPt < fMinTrigJetPt ) return;
744 
745  // Skip the event if the leading jet doesn't satisfy the leading hadron threshold
746  if (jetCont->GetLeadingHadronPt(trigJet) < fDijetLeadingHadronPt*leadingHadronCutBin) return;
747 
748  // Fill the dijet struct
749  fDijet.trigJet = trigJet;
750  fDijet.trigJetPt = trigJetPt;
751  fDijet.trigJetEta = trigJet->Eta();
752  fDijet.trigJetPhi = trigJet->Phi();
753 
754  // Find the subleading jet in the opposite hemisphere
755  AliEmcalJet *assJet = 0;
756  for(auto assJetCand : jetCont->accepted()) {
757  if (!assJetCand) continue;
758  Double_t assJetCandPt = GetJetPt(jetCont, assJetCand);
759  if ( TMath::Abs(trigJet->Phi() - assJetCand->Phi()) < fDeltaPhiMin ) continue;
760  if (assJet) {
761  Double_t assJetPt = GetJetPt(jetCont, assJet);
762  if ( assJetCandPt < assJetPt ) continue;
763  }
764  assJet = assJetCand;
765  }
766  if (!assJet) return;
767 
768  // Fill the dijet struct
769  fDijet.assJet = assJet;
770  fDijet.assJetPt = GetJetPt(jetCont, assJet);
771  fDijet.assJetPhi = assJet->Phi();
772  fDijet.assJetEta = assJet->Eta();
774 
775  fDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
776  fDijet.deltaEta = trigJet->Eta() - assJet->Eta();
779  fDijet.kTy = TMath::Abs( fDijet.trigJetPt * TMath::Sin(fDijet.deltaPhi) );
780 }
781 
787 {
788 
789  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
790 
791  AliVTrack* track;
792  for (auto trackIterator : trackCont->accepted_momentum() ) {
793 
794  track = trackIterator.second;
795 
796  // Compute the delta phi between the track and its nearest jet (of the two jets in the dijet),
797  // as well as its pT-parallel projection onto the nearest jet's axis.
798 
799  Double_t trackPt = track->Pt();
800  Double_t trackPhi = track->Phi();
801  Double_t trigJetPhi = fDijet.trigJet->Phi();
802  Double_t assJetPhi = fDijet.assJet->Phi();
803 
804  Double_t deltaPhiTrigJet = TMath::Abs(trackPhi - trigJetPhi);
805  Double_t deltaPhiAssJet = TMath::Abs(trackPhi - assJetPhi);
806  Bool_t isNearside = deltaPhiTrigJet < deltaPhiAssJet;
807 
808  Double_t deltaPhi;
809  Double_t balancePt;
810  if (isNearside) {
811  deltaPhi = trackPhi - trigJetPhi;
812  balancePt = trackPt * TMath::Cos(deltaPhi);
813  }
814  else {
815  deltaPhi = trackPhi - assJetPhi;
816  balancePt = -trackPt * TMath::Cos(deltaPhi);
817  }
818 
819  FillMomentumBalanceHistograms(histname, deltaPhi, trackPt, balancePt);
820 
821  }
822 }
823 
828 {
829  // Get jet container with minimum constituent pT,E thresholds
830  TString jetContAllName = Form("Jet_AKTFullR0%d0_tracks_pT0150_caloClusters_E0300_pt_scheme", (int) (fMatchingJetR*10) );
831  AliJetContainer* jetContAll = GetJetContainer(jetContAllName.Data());
832 
833  // Get jet container with X GeV constituent pT,E thresholds
834  Int_t trackThreshold = (int) (fTrackConstituentThreshold*1000); // in MeV
835  Int_t clusThreshold = (int) (fClusterConstituentThreshold*1000); // in MeV
836  TString jetContHardCoreName = Form("JetHardCore_AKTFullR0%d0_tracks_pT%d_caloClusters_E%d_pt_scheme", (int) (fMatchingJetR*10), trackThreshold, clusThreshold);
837  AliJetContainer* jetContHardCore = GetJetContainer(jetContHardCoreName.Data());
838 
839  // Find the di-jet in the hard-core jet sample, then find the matching di-jet and fill histograms
840  FindDijet(jetContHardCore, 0);
841  if (fDijet.isAccepted) {
842  FindMatchingDijet(jetContAll);
844  }
845 }
846 
852 {
854 
855  // Loop over jets and find leading jet within R of fDijet.trigJet
856  AliEmcalJet *matchingTrigJet = 0;
857  for(auto matchingTrigJetCand : jetCont->accepted()) {
858  if (!matchingTrigJetCand) continue;
859  if ( GetDeltaR(matchingTrigJetCand, fDijet.trigJet) > fMatchingJetR ) continue;
860  if (matchingTrigJet) {
861  if ( GetJetPt(jetCont, matchingTrigJetCand) < GetJetPt(jetCont, matchingTrigJet) ) continue;
862  }
863  matchingTrigJet = matchingTrigJetCand;
864  }
865  if (!matchingTrigJet) return;
866 
867  // Loop over jets and find leading jet within R of fDijet.assJet
868  AliEmcalJet *matchingAssJet = 0;
869  for(auto matchingAssJetCand : jetCont->accepted()) {
870  if (!matchingAssJetCand) continue;
871  if ( GetDeltaR(matchingAssJetCand, fDijet.assJet) > fMatchingJetR ) continue;
872  if (matchingAssJet) {
873  if ( GetJetPt(jetCont, matchingAssJetCand) < GetJetPt(jetCont, matchingAssJet) ) continue;
874  }
875  matchingAssJet = matchingAssJetCand;
876  }
877 
878  // Determine which matching jet is the leading jet (i.e. allow them to flip)
879  if (matchingAssJet) {
880  AliEmcalJet* trigJet = matchingTrigJet;
881  AliEmcalJet* assJet = matchingAssJet;
882  if ( GetJetPt(jetCont, matchingTrigJet) < GetJetPt(jetCont, matchingAssJet) ) {
883  AliEmcalJet* trigJet = matchingAssJet;
884  AliEmcalJet* assJet = matchingTrigJet;
885  }
886 
887  // Fill the dijet struct
888  fMatchingDijet.trigJet = trigJet;
889  fMatchingDijet.trigJetPt = GetJetPt(jetCont, trigJet);
890  fMatchingDijet.trigJetEta = trigJet->Eta();
891  fMatchingDijet.trigJetPhi = trigJet->Phi();
892 
893  fMatchingDijet.assJet = assJet;
894  fMatchingDijet.assJetPt = GetJetPt(jetCont, assJet);
895  fMatchingDijet.assJetPhi = assJet->Phi();
896  fMatchingDijet.assJetEta = assJet->Eta();
898 
899  fMatchingDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
900  fMatchingDijet.deltaEta = trigJet->Eta() - assJet->Eta();
903  fMatchingDijet.kTy = TMath::Abs( fMatchingDijet.trigJetPt * TMath::Sin(fMatchingDijet.deltaPhi) );
904  }
905 }
906 
911 {
912  // Loop over tracks and clusters in order to:
913  // (1) Compute scale factor for full jets
914  // (2) Compute delta-pT for full jets, with the random cone method
915  // For both the scale factor and delta-pT, we compute only one histogram each for EMCal.
916  // But for DCal, we bin in eta, in order to study DCal vs. PHOS vs. gap
917 
918  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
919  Double_t etaTPC = 0.9;
920  Double_t etaEMCal = 0.7;
921  Double_t etaMinDCal = 0.22;
922  Double_t etaMaxPHOS = 0.13;
923  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
924  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
925  Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
926  Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
927  Double_t phiMinPHOS = 250 * TMath::DegToRad();
928  Double_t phiMaxPHOS = 320 * TMath::DegToRad();
929 
930  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
931  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
932  Double_t accDCalRegion = 2 * etaEMCal * (phiMaxDCal - phiMinDCal);
933 
934  // Define fiducial acceptances, to be used to generate random cones
935  TRandom3* r = new TRandom3(0);
936  Double_t jetR = jetCont->GetJetRadius();
937  Double_t etaEMCalfid = etaEMCal - jetR;
938  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
939  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
940  Double_t phiMinDCalRegionfid = phiMinDCal + jetR;
941  Double_t phiMaxDCalRegionfid = phiMaxDCal - jetR;
942 
943  // Generate EMCal random cone eta-phi
944  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
945  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
946 
947  // Generate DCalRegion random cone eta-phi inside each eta slice (same phi used for all)
948  Double_t etaStep = 0.05;
949  const Int_t nEtaBinsSF = 28; // 2 * 0.7 / 0.05
950  const Int_t nEtaBinsRC = 20; // 2 * 0.5 / 0.05
951 
952  Double_t phiDCalRC = r->Uniform(phiMinDCalRegionfid, phiMaxDCalRegionfid);
953  Double_t etaDCalRC[nEtaBinsRC];
954  Double_t etaMin;
955  Double_t etaMax;
956  for (Int_t bin=0; bin < nEtaBinsRC; bin++) {
957  etaMin = -etaEMCalfid + bin*etaStep;
958  etaMax = etaMin + etaStep;
959  etaDCalRC[bin] = r->Uniform(etaMin, etaMax);
960  }
961 
962  // Initialize the various sums to 0
963  Double_t trackPtSumTPC = 0;
964  Double_t trackPtSumEMCal = 0;
965  Double_t trackPtSumEMCalRC = 0;
966  Double_t clusESumEMCal = 0;
967  Double_t clusESumEMCalRC = 0;
968  Double_t trackPtSumDCal[nEtaBinsSF] = {0.};
969  Double_t trackPtSumDCalRC[nEtaBinsRC] = {0.};
970  Double_t clusESumDCal[nEtaBinsSF] = {0.};
971  Double_t clusESumDCalRC[nEtaBinsRC] = {0.};
972 
973  // Loop over tracks. Sum the track pT:
974  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal random cone,
975  // (4) in the DCalRegion at each eta, (5) in the DCalRegion random cone at each eta
976  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
977  AliTLorentzVector track;
978  Double_t trackEta;
979  Double_t trackPhi;
980  Double_t trackPt;
981  Double_t deltaR;
982  for (auto trackIterator : trackCont->accepted_momentum() ) {
983 
984  track.Clear();
985  track = trackIterator.first;
986  trackEta = track.Eta();
987  trackPhi = track.Phi_0_2pi();
988  trackPt = track.Pt();
989 
990  // (1)
991  if (TMath::Abs(trackEta) < etaTPC) {
992  trackPtSumTPC += trackPt;
993  }
994 
995  // (2)
996  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
997  trackPtSumEMCal += trackPt;
998  }
999 
1000  // (3)
1001  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
1002  if (deltaR < jetR) {
1003  trackPtSumEMCalRC += trackPt;
1004  }
1005 
1006  // (4)
1007  for (Int_t bin=0; bin < nEtaBinsSF; bin++) {
1008  if (trackPhi > phiMinDCal && trackPhi < phiMaxDCal) {
1009  etaMin = -etaEMCal+ bin*etaStep;
1010  etaMax = etaMin + etaStep;
1011  if (trackEta > etaMin && trackEta < etaMax) {
1012  trackPtSumDCal[bin] += trackPt;
1013  }
1014  }
1015  }
1016 
1017  // (5)
1018  for (Int_t bin=0; bin < nEtaBinsRC; bin++) {
1019  deltaR = GetDeltaR(&track, etaDCalRC[bin], phiDCalRC);
1020  if (deltaR < jetR) {
1021  trackPtSumDCalRC[bin] += trackPt;
1022  }
1023  }
1024 
1025  }
1026 
1027  // Loop over clusters. Sum the cluster ET:
1028  // (1) in the EMCal, (2) in the EMCal random cone, (3) in the DCalRegion at each eta,
1029  // (4) in the DCalRegion random cone at each eta
1031  AliTLorentzVector clus;
1032  Double_t clusEta;
1033  Double_t clusPhi;
1034  Double_t clusE;
1035  for (auto clusIterator : clusCont->accepted_momentum() ) {
1036 
1037  clus.Clear();
1038  clus = clusIterator.first;
1039  clusEta = clus.Eta();
1040  clusPhi = clus.Phi_0_2pi();
1041  clusE = clus.E();
1042 
1043  // (1)
1044  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
1045  clusESumEMCal += clusE;
1046  }
1047 
1048  // (2)
1049  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
1050  if (deltaR < jetR) {
1051  clusESumEMCalRC += clusE;
1052  }
1053 
1054  // (3)
1055  for (Int_t bin=0; bin < nEtaBinsSF; bin++) {
1056  if (clusPhi > phiMinDCal && clusPhi < phiMaxDCal) {
1057  etaMin = -etaEMCal+ bin*etaStep;
1058  etaMax = etaMin + etaStep;
1059  if (clusEta > etaMin && clusEta < etaMax) {
1060  clusESumDCal[bin] += clusE;
1061  }
1062  }
1063  }
1064 
1065  // (4)
1066  for (Int_t bin=0; bin < nEtaBinsRC; bin++) {
1067  deltaR = GetDeltaR(&clus, etaDCalRC[bin], phiDCalRC);
1068  if (deltaR < jetR) {
1069  clusESumDCalRC[bin] += clusE;
1070  }
1071  }
1072 
1073  }
1074 
1075  // Compute the scale factor for EMCal
1076  Double_t numerator = (trackPtSumEMCal + clusESumEMCal) / accEMCal;
1077  Double_t denominator = trackPtSumTPC / accTPC;
1078  Double_t scaleFactor = numerator / denominator;
1079  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
1080  fHistManager.FillTH2(histname, fCent, scaleFactor);
1081 
1082  // Compute the scale factor for DCalRegion
1083  Double_t accDCalRegionBin = accDCalRegion / nEtaBinsSF;
1084  for (Int_t bin=0; bin < nEtaBinsSF; bin++) {
1085  numerator = (trackPtSumDCal[bin] + clusESumDCal[bin]) / accDCalRegionBin;
1086  scaleFactor = numerator / denominator;
1087  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorDCalRegion", jetCont->GetArrayName().Data());
1088  fHistManager.FillTH3(histname, fCent, scaleFactor, bin);
1089  }
1090 
1091  // Compute delta pT for EMCal
1092  Double_t rho = jetCont->GetRhoVal();
1093  Double_t deltaPt = trackPtSumEMCalRC + clusESumEMCalRC - rho * TMath::Pi() * jetR * jetR;
1094  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
1095  fHistManager.FillTH1(histname, deltaPt);
1096 
1097  // Compute delta pT for DCalRegion
1098  for (Int_t bin=0; bin < nEtaBinsRC; bin++) {
1099  deltaPt = trackPtSumDCalRC[bin] + clusESumDCalRC[bin] - rho * TMath::Pi() * jetR * jetR;
1100  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtDCalRegion", jetCont->GetArrayName().Data());
1101  fHistManager.FillTH2(histname, deltaPt, bin);
1102  }
1103 
1104 }
1105 
1110 {
1111  Double_t pT = jet->Pt() - jetCont->GetRhoVal() * jet->Area();
1112  TString jetContName = jetCont->GetName();
1113  if (jetContName.Contains("HardCore")) pT = jet->Pt();
1114  return pT;
1115 }
1116 
1121 {
1122  Double_t deltaPhi = TMath::Abs(jet1->Phi() - jet2->Phi());
1123  Double_t deltaEta = TMath::Abs(jet1->Eta() - jet2->Eta());
1124  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1125  return deltaR;
1126 }
1127 
1132 {
1133  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
1134  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
1135  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1136  return deltaR;
1137 }
1138 
1146 {
1148 
1149  return kTRUE;
1150 }
1151 
1157 {
1158  TString histname;
1159  AliJetContainer* jets = 0;
1160  TIter nextJetColl(&fJetCollArray);
1161  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1162  TString jetContName = jets->GetName();
1163  if (jetContName.Contains("HardCore")) continue;
1164  Double_t rhoVal = 0;
1165  if (jets->GetRhoParameter()) {
1166  rhoVal = jets->GetRhoVal();
1167  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
1168  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1169  }
1170 
1171  for (auto jet : jets->all()) {
1172 
1173  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1174  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
1175 
1176  // A vs. pT (fill before area cut)
1177  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
1178  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
1179 
1180 
1181  // Rejection reason
1182  UInt_t rejectionReason = 0;
1183  if (!jets->AcceptJet(jet, rejectionReason)) {
1184  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1185  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1186  continue;
1187  }
1188 
1189  // Centrality vs. pT
1190  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
1191  fHistManager.FillTH2(histname.Data(), corrPt, fCent);
1192 
1193  // pT vs. eta vs. phi
1194  histname = TString::Format("%s/JetHistograms/hPtVsEtaVsPhi", jets->GetArrayName().Data());
1195  fHistManager.FillTH3(histname.Data(), jet->Eta(), jet->Phi_0_2pi(), corrPt);
1196 
1197  // pT un-downscaled
1198  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
1199  fHistManager.FillTH1(histname.Data(), corrPt, fMBUpscaleFactor);
1200 
1201  // pT-leading vs. pT
1202  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1203  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1204 
1205  // Cent vs. NEF vs. pT
1206  histname = TString::Format("%s/JetHistograms/hCentVsNEFVsPt", jets->GetArrayName().Data());
1207  fHistManager.FillTH3(histname.Data(), corrPt, jet->NEF(), fCent);
1208 
1209  // Cent vs. z-leading (charged) vs. pT
1210  TLorentzVector leadPart;
1211  jets->GetLeadingHadronMomentum(leadPart, jet);
1212  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1213  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
1214  histname = TString::Format("%s/JetHistograms/hCentVsZLeadingVsPt", jets->GetArrayName().Data());
1215  fHistManager.FillTH3(histname.Data(), corrPt, z, fCent);
1216 
1217  // Cent vs. z (charged) vs. pT
1218  histname = TString::Format("%s/JetHistograms/hCentVsZVsPt", jets->GetArrayName().Data());
1219  AliVTrack* track;
1220  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1221  track = static_cast<AliVTrack*>(jet->Track(i));
1222  z = track->Pt() / corrPt;
1223  fHistManager.FillTH3(histname.Data(), corrPt, z, fCent);
1224  }
1225 
1226  // Cent vs. Nconst vs. pT
1227  histname = TString::Format("%s/JetHistograms/hCentVsNConstVsPt", jets->GetArrayName().Data());
1228  fHistManager.FillTH3(histname.Data(), corrPt, jet->GetNumberOfConstituents(), fCent);
1229 
1230  } //jet loop
1231 
1232  //---------------------------------------------------------------------------
1233  // Do study of background (if requested)
1235  }
1236 }
1237 
1242 {
1243  Double_t contents[30]={0};
1244  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1245  if (!histJetObservables) return;
1246  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1247  TString title(histJetObservables->GetAxis(n)->GetTitle());
1248  if (title=="Centrality (%)")
1249  contents[n] = fCent;
1250  else if (title=="LeadingHadronRequired")
1251  contents[n] = fDijet.leadingHadronCutType;
1252  else if (title=="#it{p}_{T,trig jet} (GeV/#it{c})")
1253  contents[n] = fDijet.trigJetPt;
1254  else if (title=="#it{p}_{T,ass jet} (GeV/#it{c})")
1255  contents[n] = fDijet.assJetPt;
1256  else if (title=="#phi_{trig jet}")
1257  contents[n] = fDijet.trigJetPhi;
1258  else if (title=="#phi_{ass jet}")
1259  contents[n] = fDijet.assJetPhi;
1260  else if (title=="#eta_{trig jet}")
1261  contents[n] = fDijet.trigJetEta;
1262  else if (title=="#eta_{ass jet}")
1263  contents[n] = fDijet.assJetEta;
1264  else
1265  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1266  }
1267  histJetObservables->Fill(contents);
1268 }
1269 
1274 {
1275  Double_t contents[30]={0};
1276  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1277  if (!histJetObservables) return;
1278  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1279  TString title(histJetObservables->GetAxis(n)->GetTitle());
1280  if (title=="Centrality (%)")
1281  contents[n] = fCent;
1282  else if (title=="#Delta#phi")
1283  contents[n] = fDijet.deltaPhi;
1284  else if (title=="#Delta#eta")
1285  contents[n] = fDijet.deltaEta;
1286  else if (title=="A_{J}")
1287  contents[n] = fDijet.AJ;
1288  else if (title=="x_{J}")
1289  contents[n] = fDijet.xJ;
1290  else if (title=="k_{Ty} (GeV)")
1291  contents[n] = fDijet.kTy;
1292  else if (title=="N_{tracks, trig jet}")
1293  contents[n] = fDijet.trigJet->GetNumberOfTracks();
1294  else if (title=="N_{tracks, ass jet}")
1295  contents[n] = fDijet.assJet->GetNumberOfTracks();
1296  else
1297  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1298  }
1299  histJetObservables->Fill(contents);
1300 }
1301 
1306 {
1307  Double_t contents[30]={0};
1308  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1309  if (!histJetObservables) return;
1310  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1311  TString title(histJetObservables->GetAxis(n)->GetTitle());
1312  if (title=="A_{J}")
1313  contents[n] = fDijet.AJ;
1314  else if (title=="#Delta#phi")
1315  contents[n] = deltaPhi;
1316  else if (title=="#it{p}_{T,particle} (GeV/#it{c})")
1317  contents[n] = trackPt;
1318  else if (title=="#it{p}_{T}#parallel (GeV/#it{c})")
1319  contents[n] = balancePt;
1320  else
1321  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1322  }
1323  histJetObservables->Fill(contents);
1324 }
1325 
1330 {
1331  // Matching efficiency histogram
1332  TString histname = "GeometricalMatchingEfficiency";
1333 
1334  // If we have a matching di-jet, fill the geometrical matching histogram
1335  if (fMatchingDijet.assJet) {
1336  fHistManager.FillTH1(histname.Data(), 1);
1337 
1340  Bool_t isSwitched = trigDeltaR > fMatchingJetR;
1341 
1342  TString thnname = "GeometricalMatching";
1343  Double_t contents[30]={0};
1344  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(thnname.Data()));
1345  if (!histJetObservables) return;
1346  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1347  TString title(histJetObservables->GetAxis(n)->GetTitle());
1348  if (title=="Centrality (%)")
1349  contents[n] = fCent;
1350  else if (title=="isSwitched")
1351  contents[n] = isSwitched;
1352  else if (title=="#DeltaR_{trig}")
1353  contents[n] = trigDeltaR;
1354  else if (title=="#DeltaR_{ass}")
1355  contents[n] = assDeltaR;
1356  else if (title=="trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1357  contents[n] = fMatchingDijet.trigJetPt - fDijet.trigJetPt;
1358  else if (title=="ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1359  contents[n] = fMatchingDijet.assJetPt - fDijet.assJetPt;
1360  else if (title=="A_{J} low-threshold")
1361  contents[n] = fMatchingDijet.AJ;
1362  else if (title=="A_{J} hard-core")
1363  contents[n] = fDijet.AJ;
1364  else
1365  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1366  }
1367  histJetObservables->Fill(contents);
1368  }
1369  else {
1370  fHistManager.FillTH1(histname.Data(), 0.);
1371  }
1372 }
1373 
Double_t Area() const
Definition: AliEmcalJet.h:123
Double_t GetRhoVal() const
const TString & GetRhoName() const
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
Double_t Eta() const
Definition: AliEmcalJet.h:114
static AliEmcalDownscaleFactorsOCDB * Instance()
const AliClusterIterableMomentumContainer accepted_momentum() const
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)
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:658
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.
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.
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
AliRhoParameter * GetRhoParameter()
Double_t Pt() const
Definition: AliEmcalJet.h:102
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
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Handler for downscale factors for various triggers obtained from the OCDB.
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
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.
const AliJetIterableContainer all() const