AliPhysics  6f1d526 (6f1d526)
 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 <TList.h>
22 #include <THnSparse.h>
23 #include <TRandom3.h>
24 #include <TGrid.h>
25 #include <TFile.h>
26 
27 #include <AliVCluster.h>
28 #include <AliVParticle.h>
29 #include <AliLog.h>
30 
31 #include "AliTLorentzVector.h"
32 #include "AliEmcalJet.h"
33 #include "AliRhoParameter.h"
34 #include "AliJetContainer.h"
35 #include "AliParticleContainer.h"
36 #include "AliClusterContainer.h"
37 #include "AliEMCALGeometry.h"
39 
41 
45 
51  fHistManager(),
52  fEventCuts(0),
53  fEventCutList(0),
54  fUseManualEventCuts(kFALSE),
55  fUseAliEventCuts(kTRUE),
56  fDeltaPhiMin(0),
57  fMinTrigJetPt(0),
58  fMinAssJetPt(0),
59  fDijetLeadingHadronPt(0),
60  fMaxPt(200),
61  fNCentHistBins(0),
62  fCentHistBins(0),
63  fPlotJetHistograms(kFALSE),
64  fPlotDijetCandHistograms(kFALSE),
65  fPlotDijetImbalanceHistograms(kFALSE),
66  fComputeBackground(kFALSE),
67  fDoMomentumBalance(kFALSE),
68  fDoGeometricalMatching(kFALSE),
69  fMatchingJetR(0.2),
70  fTrackConstituentThreshold(0),
71  fClusterConstituentThreshold(0),
72  fMBUpscaleFactor(1.),
73  fNEtaBins(40),
74  fNPhiBins(200),
75  fLoadBackgroundScalingWeights(kTRUE),
76  fBackgroundScalingWeights(0),
77  fGapJetScalingWeights(0),
78  fComputeMBDownscaling(kFALSE)
79 {
83 }
84 
91  AliAnalysisTaskEmcalJet(name, kTRUE),
92  fHistManager(name),
93  fEventCuts(0),
94  fEventCutList(0),
95  fUseManualEventCuts(kFALSE),
96  fUseAliEventCuts(kTRUE),
97  fDeltaPhiMin(0),
98  fMinTrigJetPt(0),
99  fMinAssJetPt(0),
100  fDijetLeadingHadronPt(0),
101  fMaxPt(200),
102  fNCentHistBins(0),
103  fCentHistBins(0),
104  fPlotJetHistograms(kFALSE),
105  fPlotDijetCandHistograms(kFALSE),
106  fPlotDijetImbalanceHistograms(kFALSE),
107  fComputeBackground(kFALSE),
108  fDoMomentumBalance(kFALSE),
109  fDoGeometricalMatching(kFALSE),
110  fMatchingJetR(0.2),
111  fTrackConstituentThreshold(0),
112  fClusterConstituentThreshold(0),
113  fMBUpscaleFactor(1.),
114  fNEtaBins(40),
115  fNPhiBins(200),
116  fLoadBackgroundScalingWeights(kTRUE),
117  fBackgroundScalingWeights(0),
118  fGapJetScalingWeights(0),
119  fComputeMBDownscaling(kFALSE)
120 {
122  Dijet_t fDijet;
124 }
125 
130 {
131 }
132 
137 {
138  fNCentHistBins = 4;
140  fCentHistBins[0] = 0;
141  fCentHistBins[1] = 10;
142  fCentHistBins[2] = 30;
143  fCentHistBins[3] = 50;
144  fCentHistBins[4] = 90;
145 }
146 
152 {
154 
160 
161  TIter next(fHistManager.GetListOfHistograms());
162  TObject* obj = 0;
163  while ((obj = next())) {
164  fOutput->Add(obj);
165  }
166 
167  // Intialize AliEventCuts
168  if (fUseAliEventCuts) {
169  fEventCutList = new TList();
170  fEventCutList ->SetOwner();
171  fEventCutList ->SetName("EventCutOutput");
172 
173  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
174  if(fUseManualEventCuts==1)
175  {
176  fEventCuts.SetManualMode();
177  // Configure manual settings here
178  // ...
179  }
180  fEventCuts.AddQAplotsToList(fEventCutList);
181  fOutput->Add(fEventCutList);
182  }
183 
184  // Load eta-phi background scale factors from histogram on AliEn
187  }
188 
189  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
190 }
191 
192 /*
193  * This function allocates the histograms for single jets.
194  * A set of histograms is allocated per each jet container.
195  */
197 {
198  TString histname;
199  TString title;
200 
201  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
202 
203  AliJetContainer* jets = 0;
204  TIter nextJetColl(&fJetCollArray);
205  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
206 
207  // Jet rejection reason
208  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
209  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
210  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, fMaxPt);
211  SetRejectionReasonLabels(hist->GetXaxis());
212 
213  // Rho vs. Centrality
214  if (!jets->GetRhoName().IsNull()) {
215  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
216  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
217  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 500);
218  }
219 
220  // Centrality vs. pT
221  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
222  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});Centrality (%);counts";
223  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, 10, 0, 100);
224 
225  // NEF vs. pT vs. eta vs. phi
226  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
227  title = histname + ";#eta_{jet} (rad);#phi_{jet} (rad);#it{p}_{T}^{corr} (GeV/#it{c});NEF";
228  Int_t nbins[4] = {fNEtaBins, fNPhiBins, nPtBins, 50};
229  Double_t min[4] = {-0.5,1., 0, 0};
230  Double_t max[4] = {0.5,6., fMaxPt, 1.};
231  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins, min, max);
232 
233  // pT upscaled
234  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
235  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c})";
236  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, "s");
237 
238  // pT-leading vs. pT
239  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
240  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{p}_{T,particle}^{leading} (GeV/#it{c})";
241  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt);
242 
243  // A vs. pT
244  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
245  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{A}_{jet}";
246  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/3, 0, 0.5);
247 
248  // z-leading (charged) vs. pT
249  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
250  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading}";
251  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/5, 0, 1.0);
252 
253  // z (charged) vs. pT
254  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
255  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});#it{z}";
256  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/5, 0, 1.0);
257 
258  // Nconst vs. pT
259  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
260  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents";
261  fHistManager.CreateTH2(histname.Data(), title.Data(), nPtBins, 0, fMaxPt, fMaxPt/5, 0, fMaxPt);
262 
263  // Allocate background subtraction histograms, if enabled
264  if (fComputeBackground) {
265 
266  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
267  title = histname + ";Centrality;Scale factor;counts";
268  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
269 
270  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
271  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
272  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
273 
274  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jets->GetArrayName().Data());
275  title = histname + ";#eta;#phi;Centrality;Scale factor;";
276  Int_t nbins[4] = {fNEtaBins, fNPhiBins, 50, 400};
277  Double_t min[4] = {-0.5,1., 0, 0};
278  Double_t max[4] = {0.5,6., 100, 20};
279  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins, min, max);
280 
281  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jets->GetArrayName().Data());
282  title = histname + ";#eta;#phi;Centrality;#delta#it{p}_{T} (GeV/#it{c})";
283  Int_t nbinsDpT[4] = {fNEtaBins, fNPhiBins, 10, 400};
284  Double_t minDpT[4] = {-0.5,1., 0, -50};
285  Double_t maxDpT[4] = {0.5,6., 100, 150};
286  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbinsDpT, minDpT, maxDpT);
287 
288  }
289 
290  }
291 
292  // MB downscale factor histogram
293  if (fComputeMBDownscaling) {
294  histname = "Trigger/hMBDownscaleFactor";
295  title = histname + ";Downscale factor;counts";
296  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 200, 0, 200);
297  }
298 
299 }
300 
301 /*
302  * This function allocates the histograms for dijet candidates, i.e. dijet pairs with the leading jet
303  * passing the trigger condition, but no condition on the subleading jet. In particular this histogram is
304  * designed to study the kinematic selection of dijets.
305  */
307 {
308  // Allocate dijet THnSparse
309  AliJetContainer* jets = 0;
310  TIter nextJetColl(&fJetCollArray);
311  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
312  TString axisTitle[30]= {""};
313  Int_t nbins[30] = {0};
314  Double_t min[30] = {0.};
315  Double_t max[30] = {0.};
316  Double_t *binEdges[20] = {0};
317  Int_t dim = 0;
318 
319  if (fForceBeamType != kpp) {
320  axisTitle[dim] = "Centrality (%)";
321  nbins[dim] = fNCentHistBins;
322  binEdges[dim] = fCentHistBins;
323  min[dim] = fCentHistBins[0];
324  max[dim] = fCentHistBins[fNCentHistBins];
325  dim++;
326  }
327 
328  axisTitle[dim] = "LeadingHadronRequired";
329  nbins[dim] = 2;
330  min[dim] = -0.5;
331  max[dim] = 1.5;
332  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
333  dim++;
334 
335  axisTitle[dim] = "#it{p}_{T,trig jet} (GeV/#it{c})";
336  nbins[dim] = TMath::CeilNint(fMaxPt/2);
337  min[dim] = 0;
338  max[dim] = fMaxPt;
339  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
340  dim++;
341 
342  axisTitle[dim] = "#it{p}_{T,ass jet} (GeV/#it{c})";
343  nbins[dim] = TMath::CeilNint(fMaxPt/2);
344  min[dim] = 0;
345  max[dim] = fMaxPt;
346  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
347  dim++;
348 
349  axisTitle[dim] = "#phi_{trig jet}";
350  nbins[dim] = TMath::CeilNint(fMaxPt/2);
351  min[dim] = 0;
352  max[dim] = TMath::TwoPi();
353  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
354  dim++;
355 
356  axisTitle[dim] = "#phi_{ass jet}";
357  nbins[dim] = TMath::CeilNint(fMaxPt/2);
358  min[dim] = 0;
359  max[dim] = TMath::TwoPi();
360  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
361  dim++;
362 
363  axisTitle[dim] = "#eta_{trig jet}";
364  nbins[dim] = TMath::CeilNint(fMaxPt/2);
365  min[dim] = -1;
366  max[dim] = 1;
367  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
368  dim++;
369 
370  axisTitle[dim] = "#eta_{ass jet}";
371  nbins[dim] = TMath::CeilNint(fMaxPt/2);
372  min[dim] = -1;
373  max[dim] = 1;
374  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
375  dim++;
376 
377  TString thnname = TString::Format("%s/DijetCandObservables", jets->GetArrayName().Data());
378  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
379  for (Int_t i = 0; i < dim; i++) {
380  hn->GetAxis(i)->SetTitle(axisTitle[i]);
381  hn->SetBinEdges(i, binEdges[i]);
382  }
383  }
384 }
385 
386 /*
387  * This function allocates the histograms for accepted dijets.
388  * The purpose is to study in detail the imbalance properties of the accepted dijets.
389  */
391 {
392  // Allocate dijet imbalance THnSparse
393  AliJetContainer* jets = 0;
394  TIter nextJetColl(&fJetCollArray);
395  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
396  TString axisTitle[30]= {""};
397  Int_t nbins[30] = {0};
398  Double_t min[30] = {0.};
399  Double_t max[30] = {0.};
400  Double_t *binEdges[20] = {0};
401  Int_t dim = 0;
402 
403  if (fForceBeamType != kpp) {
404  axisTitle[dim] = "Centrality (%)";
405  nbins[dim] = fNCentHistBins;
406  binEdges[dim] = fCentHistBins;
407  min[dim] = fCentHistBins[0];
408  max[dim] = fCentHistBins[fNCentHistBins];
409  dim++;
410  }
411 
412  axisTitle[dim] = "#Delta#phi";
413  nbins[dim] = 100;
414  min[dim] = 0;
415  max[dim] = 4;
416  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
417  dim++;
418 
419  axisTitle[dim] = "#Delta#eta";
420  nbins[dim] = 100;
421  min[dim] = -2;
422  max[dim] = 2;
423  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
424  dim++;
425 
426  axisTitle[dim] = "A_{J}";
427  nbins[dim] = 100;
428  min[dim] = 0;
429  max[dim] = 1;
430  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
431  dim++;
432 
433  axisTitle[dim] = "x_{J}";
434  nbins[dim] = 100;
435  min[dim] = 0;
436  max[dim] = 1;
437  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
438  dim++;
439 
440  axisTitle[dim] = "k_{Ty} (GeV)";
441  nbins[dim] = 100;
442  min[dim] = 0;
443  max[dim] = 100;
444  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
445  dim++;
446 
447  axisTitle[dim] = "N_{tracks, trig jet}";
448  nbins[dim] = fMaxPt/5;
449  min[dim] = 0;
450  max[dim] = 100;
451  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
452  dim++;
453 
454  axisTitle[dim] = "N_{tracks, ass jet}";
455  nbins[dim] = fMaxPt/5;
456  min[dim] = 0;
457  max[dim] = 100;
458  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
459  dim++;
460 
461  TString thnname = TString::Format("%s/DijetImbalanceObservables", jets->GetArrayName().Data());
462  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
463  for (Int_t i = 0; i < dim; i++) {
464  hn->GetAxis(i)->SetTitle(axisTitle[i]);
465  hn->SetBinEdges(i, binEdges[i]);
466  }
467  }
468 }
469 
470 /*
471  * This function allocates the histograms for the momentum balance study.
472  */
474 {
475  AliJetContainer* jets = 0;
476  TIter nextJetColl(&fJetCollArray);
477  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
478 
479  // Allocate THnSparse
480  TString axisTitle[30]= {""};
481  Int_t nbins[30] = {0};
482  Double_t min[30] = {0.};
483  Double_t max[30] = {0.};
484  Double_t *binEdges[20] = {0};
485  Int_t dim = 0;
486 
487  axisTitle[dim] = "A_{J}";
488  nbins[dim] = 100;
489  min[dim] = 0;
490  max[dim] = 1;
491  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
492  dim++;
493 
494  axisTitle[dim] = "#Delta#phi";
495  nbins[dim] = 100;
496  min[dim] = -4;
497  max[dim] = 4;
498  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
499  dim++;
500 
501  axisTitle[dim] = "#it{p}_{T,particle} (GeV/#it{c})";
502  nbins[dim] = 9;
503  Double_t* pTParticleBins = new Double_t[nbins[dim]+1];
504  GenerateFixedBinArray(1, 0.15, 0.3, pTParticleBins);
505  GenerateFixedBinArray(1, 0.3, 0.5, pTParticleBins+1);
506  GenerateFixedBinArray(1, 0.5, 1, pTParticleBins+2);
507  GenerateFixedBinArray(2, 1, 5, pTParticleBins+3);
508  GenerateFixedBinArray(3, 5, 20, pTParticleBins+5);
509  GenerateFixedBinArray(1, 20, 150, pTParticleBins+8);
510  min[dim] = 0;
511  max[dim] = pTParticleBins[nbins[dim]];
512  binEdges[dim] = pTParticleBins;
513  dim++;
514 
515  axisTitle[dim] = "#it{p}_{T}#parallel (GeV/#it{c})";
516  nbins[dim] = 80;
517  Double_t* pTParallelBins = new Double_t[nbins[dim]+1];
518  GenerateFixedBinArray(20, 0, 2, pTParallelBins);
519  GenerateFixedBinArray(16, 2, 10, pTParallelBins+20);
520  GenerateFixedBinArray(10, 10, 20, pTParallelBins+36);
521  GenerateFixedBinArray(10, 20, 40, pTParallelBins+46);
522  GenerateFixedBinArray(24, 40, 150, pTParallelBins+56);
523  min[dim] = 0;
524  max[dim] = pTParallelBins[nbins[dim]];
525  binEdges[dim] = pTParallelBins;
526  dim++;
527 
528  TString thnname = TString::Format("%s/MomentumBalance", jets->GetArrayName().Data());
529  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
530  for (Int_t i = 0; i < dim; i++) {
531  hn->GetAxis(i)->SetTitle(axisTitle[i]);
532  hn->SetBinEdges(i, binEdges[i]);
533  }
534  }
535 }
536 
537 /*
538  * This function allocates the histograms for the constituent dijet study.
539  */
541 {
542  // Allocate geometrical matching THnSparse
543  TString axisTitle[30]= {""};
544  Int_t nbins[30] = {0};
545  Double_t min[30] = {0.};
546  Double_t max[30] = {0.};
547  Double_t *binEdges[20] = {0};
548  Int_t dim = 0;
549 
550  if (fForceBeamType != kpp) {
551  axisTitle[dim] = "Centrality (%)";
552  nbins[dim] = fNCentHistBins;
553  binEdges[dim] = fCentHistBins;
554  min[dim] = fCentHistBins[0];
555  max[dim] = fCentHistBins[fNCentHistBins];
556  dim++;
557  }
558 
559  axisTitle[dim] = "isSwitched";
560  nbins[dim] = 2;
561  min[dim] = -0.5;
562  max[dim] = 1.5;
563  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
564  dim++;
565 
566  axisTitle[dim] = "#DeltaR_{trig}";
567  nbins[dim] = 50;
568  min[dim] = 0;
569  max[dim] = 0.5;
570  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
571  dim++;
572 
573  axisTitle[dim] = "#DeltaR_{ass}";
574  nbins[dim] = 50;
575  min[dim] = 0;
576  max[dim] = 0.5;
577  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
578  dim++;
579 
580  axisTitle[dim] = "trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
581  nbins[dim] = 100;
582  min[dim] = -50;
583  max[dim] = 50;
584  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
585  dim++;
586 
587  axisTitle[dim] = "ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
588  nbins[dim] = 100;
589  min[dim] = -50;
590  max[dim] = 50;
591  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
592  dim++;
593 
594  axisTitle[dim] = "A_{J} low-threshold";
595  nbins[dim] = 100;
596  min[dim] = 0;
597  max[dim] = 1;
598  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
599  dim++;
600 
601  axisTitle[dim] = "A_{J} hard-core";
602  nbins[dim] = 100;
603  min[dim] = 0;
604  max[dim] = 1;
605  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
606  dim++;
607 
608  TString thnname = "GeometricalMatching";
609  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
610  for (Int_t i = 0; i < dim; i++) {
611  hn->GetAxis(i)->SetTitle(axisTitle[i]);
612  hn->SetBinEdges(i, binEdges[i]);
613  }
614 
615  // Allocate other histograms
616  TString histname;
617  TString title;
618  histname = "GeometricalMatchingEfficiency";
619  title = histname + ";isMatched;counts";
620  TH1* hist = fHistManager.CreateTH1(histname.Data(), title.Data(), 2, -0.5, 1.5);
621 }
622 
627 {
628 
629  TString fname(path);
630  if (fname.BeginsWith("alien://")) {
631  TGrid::Connect("alien://");
632  }
633 
634  TFile* file = TFile::Open(path);
635 
636  if (!file || file->IsZombie()) {
637  ::Error("AliAnalysisTaskEmcalDijetImbalance", "Could not open background scaling histogram");
638  return;
639  }
640 
641  TH1D* h = dynamic_cast<TH1D*>(file->Get(name));
642 
643  if (h) {
644  ::Info("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s loaded from file %s.", name, path);
645  }
646  else {
647  ::Error("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s not found in file %s.", name, path);
648  return;
649  }
650 
651  fBackgroundScalingWeights = static_cast<TH1D*>(h->Clone());
652 
653  file->Close();
654  delete file;
655 
656 }
657 
663 {
665 
666  fNeedEmcalGeom = kTRUE;
667 
668  AliInfo(Form("Trigger jet threshold = %f, Associated jet threshold = %f", fMinTrigJetPt, fMinAssJetPt));
669  AliInfo(Form("Leading hadron threshold (for dijet leading jet): %f GeV", fDijetLeadingHadronPt));
670  AliInfo(Form("Momentum balance study: %d", fDoMomentumBalance));
671  AliInfo(Form("Geometrical matching study: %d", fDoGeometricalMatching));
672 
673 }
674 
679 
680  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
681 
683 
684  // Get instance of the downscale factor helper class
686  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
687 
688  // There are two possible min bias triggers for LHC15o
689  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
690  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
691  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
692 
693  // Get the downscale factor for whichever MB trigger exists in the given run
694  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
695  Double_t downscalefactor;
696  for (auto i : runtriggers) {
697  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
698  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
699  break;
700  }
701  }
702 
703  // Store the inverse of the downscale factor, used later to weight the pT spectrum
704  fMBUpscaleFactor = 1/downscalefactor;
705 
706  TString histname = "Trigger/hMBDownscaleFactor";
707  fHistManager.FillTH1(histname.Data(), fMBUpscaleFactor);
708 
709  }
710 
711 }
712 
717 {
718  if (fUseAliEventCuts) {
719  if (!fEventCuts.AcceptEvent(InputEvent()))
720  {
721  PostData(1, fOutput);
722  return kFALSE;
723  }
724  }
725  else {
727  }
728  return kTRUE;
729 }
730 
739 {
740  TString histname;
741  AliJetContainer* jetCont = 0;
742  TIter next(&fJetCollArray);
743  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
744  TString jetContName = jetCont->GetName();
745  if (jetContName.Contains("HardCore")) continue;
746 
747  //-----------------------------------------------------------------------------
748  // Find the leading di-jet candidate in each event, and if it satisfies the
749  // trig jet pT threshold, then fill di-jet candidate histogram (regardless of ass jet).
750  // The idea is to study the kinematic selections in post-processing.
751 
752  // Loop over leading hadron cut or not
753  for (Int_t leadingHadronCutType=0; leadingHadronCutType<2; leadingHadronCutType++) {
754 
755  // Find the dijet candidate of the event and store its info in struct fDijet
756  FindDijet(jetCont, leadingHadronCutType);
757 
758  // If we find a dijet candidate (i.e. acceptable trig jet; ass jet accepted or not), fill the di-jet candidate histogram
760  TString histname = TString::Format("%s/DijetCandObservables", jetCont->GetArrayName().Data());
761  FillDijetCandHistograms(histname);
762  }
763 
764  }
765 
766  //---------------------------------------------------------------------------------------------------
767  // Now, study the accepted dijet selection -- specified by the trig/ass jet pT conditions
768 
769  // Find the dijet candidate of the event and store its info in struct fDijet
770  FindDijet(jetCont, 0);
771 
772  // If we find an accepted dijet, fill the dijet imbalance histogram
774  TString histname = TString::Format("%s/DijetImbalanceObservables", jetCont->GetArrayName().Data());
776  }
777 
778  // If we find an accepted dijet, perform momentum-balance study (if requested)
780  histname = TString::Format("%s/MomentumBalance", jetCont->GetArrayName().Data());
781  DoMomentumBalance(histname);
782  }
783 
784  }
785 
786  //---------------------------------------------------------------------------
787  // Do the constituent threshold and geometrical matching study (if requested)
790 
791  return kTRUE;
792 }
793 
802 {
803  fDijet.clear();
804  fDijet.leadingHadronCutType = leadingHadronCutBin;
805 
806  // Get trigger jet
807  AliEmcalJet* trigJet = GetLeadingJet(jetCont);
808  if(!trigJet) return;
809 
810  // Skip the event if the leading jet doesn't satisfy the pT threshold
811  Double_t trigJetPt = GetJetPt(jetCont, trigJet);
812  if ( trigJetPt < fMinTrigJetPt ) return;
813 
814  // Skip the event if the leading jet doesn't satisfy the leading hadron threshold
815  if (jetCont->GetLeadingHadronPt(trigJet) < fDijetLeadingHadronPt*leadingHadronCutBin) return;
816 
817  // Fill the dijet struct
818  fDijet.trigJet = trigJet;
819  fDijet.trigJetPt = trigJetPt;
820  fDijet.trigJetEta = trigJet->Eta();
821  fDijet.trigJetPhi = trigJet->Phi();
822 
823  // Find the subleading jet in the opposite hemisphere
824  AliEmcalJet *assJet = 0;
825  for(auto assJetCand : jetCont->accepted()) {
826  if (!assJetCand) continue;
827  Double_t assJetCandPt = GetJetPt(jetCont, assJetCand);
828  if ( TMath::Abs(trigJet->Phi() - assJetCand->Phi()) < fDeltaPhiMin ) continue;
829  if (assJet) {
830  Double_t assJetPt = GetJetPt(jetCont, assJet);
831  if ( assJetCandPt < assJetPt ) continue;
832  }
833  assJet = assJetCand;
834  }
835  if (!assJet) return;
836 
837  // Fill the dijet struct
838  fDijet.assJet = assJet;
839  fDijet.assJetPt = GetJetPt(jetCont, assJet);
840  fDijet.assJetPhi = assJet->Phi();
841  fDijet.assJetEta = assJet->Eta();
843 
844  fDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
845  fDijet.deltaEta = trigJet->Eta() - assJet->Eta();
848  fDijet.kTy = TMath::Abs( fDijet.trigJetPt * TMath::Sin(fDijet.deltaPhi) );
849 }
850 
856 {
857 
858  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
859 
860  AliVTrack* track;
861  for (auto trackIterator : trackCont->accepted_momentum() ) {
862 
863  track = trackIterator.second;
864 
865  // Compute the delta phi between the track and its nearest jet (of the two jets in the dijet),
866  // as well as its pT-parallel projection onto the nearest jet's axis.
867 
868  Double_t trackPt = track->Pt();
869  Double_t trackPhi = track->Phi();
870  Double_t trigJetPhi = fDijet.trigJet->Phi();
871  Double_t assJetPhi = fDijet.assJet->Phi();
872 
873  Double_t deltaPhiTrigJet = TMath::Abs(trackPhi - trigJetPhi);
874  Double_t deltaPhiAssJet = TMath::Abs(trackPhi - assJetPhi);
875  Bool_t isNearside = deltaPhiTrigJet < deltaPhiAssJet;
876 
877  Double_t deltaPhi;
878  Double_t balancePt;
879  if (isNearside) {
880  deltaPhi = trackPhi - trigJetPhi;
881  balancePt = trackPt * TMath::Cos(deltaPhi);
882  }
883  else {
884  deltaPhi = trackPhi - assJetPhi;
885  balancePt = -trackPt * TMath::Cos(deltaPhi);
886  }
887 
888  FillMomentumBalanceHistograms(histname, deltaPhi, trackPt, balancePt);
889 
890  }
891 }
892 
897 {
898  // Get jet container with minimum constituent pT,E thresholds
899  TString jetContAllName = Form("Jet_AKTFullR0%d0_tracks_pT0150_caloClusters_E0300_pt_scheme", (int) (fMatchingJetR*10) );
900  AliJetContainer* jetContAll = GetJetContainer(jetContAllName.Data());
901 
902  // Get jet container with X GeV constituent pT,E thresholds
903  Int_t trackThreshold = (int) (fTrackConstituentThreshold*1000); // in MeV
904  Int_t clusThreshold = (int) (fClusterConstituentThreshold*1000); // in MeV
905  TString jetContHardCoreName = Form("JetHardCore_AKTFullR0%d0_tracks_pT%d_caloClusters_E%d_pt_scheme", (int) (fMatchingJetR*10), trackThreshold, clusThreshold);
906  AliJetContainer* jetContHardCore = GetJetContainer(jetContHardCoreName.Data());
907 
908  // Find the di-jet in the hard-core jet sample, then find the matching di-jet and fill histograms
909  FindDijet(jetContHardCore, 0);
910  if (fDijet.isAccepted) {
911  FindMatchingDijet(jetContAll);
913  }
914 }
915 
921 {
923 
924  // Loop over jets and find leading jet within R of fDijet.trigJet
925  AliEmcalJet *matchingTrigJet = 0;
926  for(auto matchingTrigJetCand : jetCont->accepted()) {
927  if (!matchingTrigJetCand) continue;
928  if ( GetDeltaR(matchingTrigJetCand, fDijet.trigJet) > fMatchingJetR ) continue;
929  if (matchingTrigJet) {
930  if ( GetJetPt(jetCont, matchingTrigJetCand) < GetJetPt(jetCont, matchingTrigJet) ) continue;
931  }
932  matchingTrigJet = matchingTrigJetCand;
933  }
934  if (!matchingTrigJet) return;
935 
936  // Loop over jets and find leading jet within R of fDijet.assJet
937  AliEmcalJet *matchingAssJet = 0;
938  for(auto matchingAssJetCand : jetCont->accepted()) {
939  if (!matchingAssJetCand) continue;
940  if ( GetDeltaR(matchingAssJetCand, fDijet.assJet) > fMatchingJetR ) continue;
941  if (matchingAssJet) {
942  if ( GetJetPt(jetCont, matchingAssJetCand) < GetJetPt(jetCont, matchingAssJet) ) continue;
943  }
944  matchingAssJet = matchingAssJetCand;
945  }
946 
947  // Determine which matching jet is the leading jet (i.e. allow them to flip)
948  if (matchingAssJet) {
949  AliEmcalJet* trigJet = matchingTrigJet;
950  AliEmcalJet* assJet = matchingAssJet;
951  if ( GetJetPt(jetCont, matchingTrigJet) < GetJetPt(jetCont, matchingAssJet) ) {
952  AliEmcalJet* trigJet = matchingAssJet;
953  AliEmcalJet* assJet = matchingTrigJet;
954  }
955 
956  // Fill the dijet struct
957  fMatchingDijet.trigJet = trigJet;
958  fMatchingDijet.trigJetPt = GetJetPt(jetCont, trigJet);
959  fMatchingDijet.trigJetEta = trigJet->Eta();
960  fMatchingDijet.trigJetPhi = trigJet->Phi();
961 
962  fMatchingDijet.assJet = assJet;
963  fMatchingDijet.assJetPt = GetJetPt(jetCont, assJet);
964  fMatchingDijet.assJetPhi = assJet->Phi();
965  fMatchingDijet.assJetEta = assJet->Eta();
967 
968  fMatchingDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
969  fMatchingDijet.deltaEta = trigJet->Eta() - assJet->Eta();
972  fMatchingDijet.kTy = TMath::Abs( fMatchingDijet.trigJetPt * TMath::Sin(fMatchingDijet.deltaPhi) );
973  }
974 }
975 
980 {
981  // Loop over tracks and clusters in order to:
982  // (1) Compute scale factor for full jets
983  // (2) Compute delta-pT for full jets, with the random cone method
984  // For both the scale factor and delta-pT, we compute only one histogram each for EMCal.
985  // But for DCal, we bin in eta-phi, in order to account for the DCal vs. PHOS vs. gap
986 
987  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
988  Double_t etaTPC = 0.9;
989  Double_t etaEMCal = 0.7;
990  Double_t etaMinDCal = 0.22;
991  Double_t etaMaxPHOS = 0.13;
992  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
993  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
994  Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
995  Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
996  Double_t phiMinPHOS = 250 * TMath::DegToRad();
997  Double_t phiMaxPHOS = 320 * TMath::DegToRad();
998 
999  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
1000  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
1001  Double_t accDCalRegion = 2 * etaEMCal * (phiMaxDCal - phiMinDCal);
1002 
1003  // Define fiducial acceptances, to be used to generate random cones
1004  TRandom3* r = new TRandom3(0);
1005  Double_t jetR = jetCont->GetJetRadius();
1006  Double_t accRC = TMath::Pi() * jetR * jetR;
1007  Double_t etaEMCalfid = etaEMCal - jetR;
1008  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
1009  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
1010  Double_t phiMinDCalRegionfid = phiMinDCal + jetR;
1011  Double_t phiMaxDCalRegionfid = phiMaxDCal - jetR;
1012 
1013  // Generate EMCal random cone eta-phi
1014  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
1015  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
1016 
1017  // For DCalRegion, generate random eta, phi in each eta/phi bin, to be used as center of random cone
1018  Double_t etaDCalRC[fNEtaBins]; // array storing the RC eta values
1019  Double_t etaStep = 1./fNEtaBins;
1020  Double_t etaMin;
1021  Double_t etaMax;
1022  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1023  etaMin = -etaEMCalfid + etaBin*etaStep;
1024  etaMax = etaMin + etaStep;
1025  etaDCalRC[etaBin] = r->Uniform(etaMin, etaMax);
1026  }
1027 
1028  Double_t phiDCalRC[fNPhiBins]; // array storing the RC phi values
1029  Double_t phiStep = 5./fNPhiBins; // phi axis is [1,6] in order to have simple binning
1030  Double_t phiMin;
1031  Double_t phiMax;
1032  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1033  phiMin = 1 + phiBin*phiStep;
1034  phiMax = phiMin + phiStep;
1035  phiDCalRC[phiBin] = r->Uniform(phiMin, phiMax);
1036  }
1037 
1038  // Initialize the various sums to 0
1039  Double_t trackPtSumTPC = 0;
1040  Double_t trackPtSumEMCal = 0;
1041  Double_t trackPtSumEMCalRC = 0;
1042  Double_t clusESumEMCal = 0;
1043  Double_t clusESumEMCalRC = 0;
1044 
1045  // Define a 2D vector (initialized to 0) to store the sum of track pT, and another for cluster ET
1046  std::vector<std::vector<Double_t>> trackPtSumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1047  std::vector<std::vector<Double_t>> clusESumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1048 
1049  // Loop over tracks. Sum the track pT:
1050  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal random cone,
1051  // (4) in the DCalRegion in a random cone at each eta-phi
1052  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1053  AliTLorentzVector track;
1054  Double_t trackEta;
1055  Double_t trackPhi;
1056  Double_t trackPt;
1057  Double_t deltaR;
1058  for (auto trackIterator : trackCont->accepted_momentum() ) {
1059 
1060  track.Clear();
1061  track = trackIterator.first;
1062  trackEta = track.Eta();
1063  trackPhi = track.Phi_0_2pi();
1064  trackPt = track.Pt();
1065 
1066  // (1)
1067  if (TMath::Abs(trackEta) < etaTPC) {
1068  trackPtSumTPC += trackPt;
1069  }
1070 
1071  // (2)
1072  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
1073  trackPtSumEMCal += trackPt;
1074  }
1075 
1076  // (3)
1077  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
1078  if (deltaR < jetR) {
1079  trackPtSumEMCalRC += trackPt;
1080  }
1081 
1082  // (4)
1083  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1084  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1085  deltaR = GetDeltaR(&track, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1086  if (deltaR < jetR) {
1087  trackPtSumDCalRC[etaBin][phiBin] += trackPt;
1088  }
1089  }
1090  }
1091 
1092  }
1093 
1094  // Loop over clusters. Sum the cluster ET:
1095  // (1) in the EMCal, (2) in the EMCal random cone, (3) in the DCalRegion in a random cone at each eta-phi
1097  AliTLorentzVector clus;
1098  Double_t clusEta;
1099  Double_t clusPhi;
1100  Double_t clusE;
1101  for (auto clusIterator : clusCont->accepted_momentum() ) {
1102 
1103  clus.Clear();
1104  clus = clusIterator.first;
1105  clusEta = clus.Eta();
1106  clusPhi = clus.Phi_0_2pi();
1107  clusE = clus.E();
1108 
1109  // (1)
1110  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
1111  clusESumEMCal += clusE;
1112  }
1113 
1114  // (2)
1115  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
1116  if (deltaR < jetR) {
1117  clusESumEMCalRC += clusE;
1118  }
1119 
1120  // (3)
1121  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1122  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1123  deltaR = GetDeltaR(&clus, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1124  if (deltaR < jetR) {
1125  clusESumDCalRC[etaBin][phiBin] += clusE;
1126  }
1127  }
1128  }
1129 
1130  }
1131 
1132  // Compute the scale factor for EMCal, as a function of centrality
1133  Double_t numerator = (trackPtSumEMCal + clusESumEMCal) / accEMCal;
1134  Double_t denominator = trackPtSumTPC / accTPC;
1135  Double_t scaleFactor = numerator / denominator;
1136  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
1137  fHistManager.FillTH2(histname, fCent, scaleFactor);
1138 
1139  // Compute the scale factor for DCalRegion in each eta-phi bin, as a function of centrality
1140  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1141  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1142  numerator = (trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin]) / accRC;
1143  scaleFactor = numerator / denominator;
1144  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jetCont->GetArrayName().Data());
1145  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, scaleFactor};
1146  fHistManager.FillTHnSparse(histname, x);
1147  }
1148  }
1149 
1150  // Compute delta pT for EMCal, as a function of centrality
1151  Double_t rho = jetCont->GetRhoVal();
1152  Double_t deltaPt = trackPtSumEMCalRC + clusESumEMCalRC - rho * TMath::Pi() * jetR * jetR;
1153  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
1154  fHistManager.FillTH2(histname, fCent, deltaPt);
1155 
1156  // Compute delta pT for DCalRegion in each eta-phi bin, as a function of centrality
1157  Double_t sf;
1158  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1159  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1161  sf = fBackgroundScalingWeights->GetBinContent(etaDCalRC[etaBin], phiDCalRC[phiBin]);
1162  rho = sf * rho;
1163  }
1164  deltaPt = trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin] - rho * accRC;
1165  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jetCont->GetArrayName().Data());
1166  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, deltaPt};
1167  fHistManager.FillTHnSparse(histname, x);
1168  }
1169  }
1170 
1171  delete r;
1172 
1173 }
1174 
1179 {
1180 
1181  Double_t rho = jetCont->GetRhoVal();
1182 
1183  // Get eta-phi dependent scale factor
1185  Double_t sf = fBackgroundScalingWeights->GetBinContent(jet->Eta(), jet->Phi_0_2pi());
1186  rho = sf * rho;
1187  }
1188 
1189  // Compute pT
1190  Double_t pT = jet->Pt() - rho * jet->Area();
1191 
1192  // If hard-core jet, don't subtract background
1193  TString jetContName = jetCont->GetName();
1194  if (jetContName.Contains("HardCore")) pT = jet->Pt();
1195 
1196  return pT;
1197 }
1198 
1203 {
1204  AliEmcalJet* leadingJet = 0;
1205 
1206  if (jetCont->GetRhoParameter()) {
1207  for(auto jetCand : jetCont->accepted()) {
1208  if (!jetCand) continue;
1209  Double_t jetCandPt = GetJetPt(jetCont, jetCand);
1210  if (leadingJet) {
1211  Double_t leadingJetPt = GetJetPt(jetCont, leadingJet);
1212  if ( jetCandPt < leadingJetPt ) continue;
1213  }
1214  leadingJet = jetCand;
1215  }
1216  }
1217  else {
1218  leadingJet = jetCont->GetLeadingJet();
1219  }
1220 
1221  return leadingJet;
1222 }
1223 
1228 {
1229  Double_t deltaPhi = TMath::Abs(jet1->Phi() - jet2->Phi());
1230  Double_t deltaEta = TMath::Abs(jet1->Eta() - jet2->Eta());
1231  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1232  return deltaR;
1233 }
1234 
1239 {
1240  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
1241  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
1242  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1243  return deltaR;
1244 }
1245 
1253 {
1255 
1256  return kTRUE;
1257 }
1258 
1264 {
1265  TString histname;
1266  AliJetContainer* jets = 0;
1267  TIter nextJetColl(&fJetCollArray);
1268  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1269  TString jetContName = jets->GetName();
1270  if (jetContName.Contains("HardCore")) continue;
1271  Double_t rhoVal = 0;
1272  if (jets->GetRhoParameter()) {
1273  rhoVal = jets->GetRhoVal();
1274  histname = TString::Format("%s/JetHistograms/hRhoVsCent", jets->GetArrayName().Data());
1275  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1276  }
1277 
1278  for (auto jet : jets->all()) {
1279 
1280  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1281  Float_t corrPt = GetJetPt(jets, jet);
1282 
1283  // A vs. pT (fill before area cut)
1284  histname = TString::Format("%s/JetHistograms/hAreaVsPt", jets->GetArrayName().Data());
1285  fHistManager.FillTH2(histname.Data(), corrPt, jet->Area());
1286 
1287 
1288  // Rejection reason
1289  UInt_t rejectionReason = 0;
1290  if (!jets->AcceptJet(jet, rejectionReason)) {
1291  histname = TString::Format("%s/JetHistograms/hJetRejectionReason", jets->GetArrayName().Data());
1292  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1293  continue;
1294  }
1295 
1296  // Centrality vs. pT
1297  histname = TString::Format("%s/JetHistograms/hCentVsPt", jets->GetArrayName().Data());
1298  fHistManager.FillTH2(histname.Data(), corrPt, fCent);
1299 
1300  // NEF vs. pT vs. eta vs. phi
1301  histname = TString::Format("%s/JetHistograms/hNEFVsPtVsEtaVsPhi", jets->GetArrayName().Data());
1302  Double_t x[4] = {jet->Eta(), jet->Phi_0_2pi(), corrPt, jet->NEF()};
1303  fHistManager.FillTHnSparse(histname, x);
1304 
1305  // pT un-downscaled
1306  histname = TString::Format("%s/JetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
1307  fHistManager.FillTH1(histname.Data(), corrPt, fMBUpscaleFactor);
1308 
1309  // pT-leading vs. pT
1310  histname = TString::Format("%s/JetHistograms/hPtLeadingVsPt", jets->GetArrayName().Data());
1311  fHistManager.FillTH2(histname.Data(), corrPt, ptLeading);
1312 
1313  // z-leading (charged) vs. pT
1314  TLorentzVector leadPart;
1315  jets->GetLeadingHadronMomentum(leadPart, jet);
1316  Double_t z = GetParallelFraction(leadPart.Vect(), jet);
1317  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
1318  histname = TString::Format("%s/JetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
1319  fHistManager.FillTH2(histname.Data(), corrPt, z);
1320 
1321  // z (charged) vs. pT
1322  histname = TString::Format("%s/JetHistograms/hZVsPt", jets->GetArrayName().Data());
1323  AliVTrack* track;
1324  for (Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
1325  track = static_cast<AliVTrack*>(jet->Track(i));
1326  z = track->Pt() / TMath::Abs(corrPt);
1327  fHistManager.FillTH2(histname.Data(), corrPt, z);
1328  }
1329 
1330  // Nconst vs. pT
1331  histname = TString::Format("%s/JetHistograms/hNConstVsPt", jets->GetArrayName().Data());
1332  fHistManager.FillTH2(histname.Data(), corrPt, jet->GetNumberOfConstituents());
1333 
1334  } //jet loop
1335 
1336  //---------------------------------------------------------------------------
1337  // Do study of background (if requested)
1339  }
1340 }
1341 
1346 {
1347  Double_t contents[30]={0};
1348  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1349  if (!histJetObservables) return;
1350  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1351  TString title(histJetObservables->GetAxis(n)->GetTitle());
1352  if (title=="Centrality (%)")
1353  contents[n] = fCent;
1354  else if (title=="LeadingHadronRequired")
1355  contents[n] = fDijet.leadingHadronCutType;
1356  else if (title=="#it{p}_{T,trig jet} (GeV/#it{c})")
1357  contents[n] = fDijet.trigJetPt;
1358  else if (title=="#it{p}_{T,ass jet} (GeV/#it{c})")
1359  contents[n] = fDijet.assJetPt;
1360  else if (title=="#phi_{trig jet}")
1361  contents[n] = fDijet.trigJetPhi;
1362  else if (title=="#phi_{ass jet}")
1363  contents[n] = fDijet.assJetPhi;
1364  else if (title=="#eta_{trig jet}")
1365  contents[n] = fDijet.trigJetEta;
1366  else if (title=="#eta_{ass jet}")
1367  contents[n] = fDijet.assJetEta;
1368  else
1369  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1370  }
1371  histJetObservables->Fill(contents);
1372 }
1373 
1378 {
1379  Double_t contents[30]={0};
1380  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1381  if (!histJetObservables) return;
1382  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1383  TString title(histJetObservables->GetAxis(n)->GetTitle());
1384  if (title=="Centrality (%)")
1385  contents[n] = fCent;
1386  else if (title=="#Delta#phi")
1387  contents[n] = fDijet.deltaPhi;
1388  else if (title=="#Delta#eta")
1389  contents[n] = fDijet.deltaEta;
1390  else if (title=="A_{J}")
1391  contents[n] = fDijet.AJ;
1392  else if (title=="x_{J}")
1393  contents[n] = fDijet.xJ;
1394  else if (title=="k_{Ty} (GeV)")
1395  contents[n] = fDijet.kTy;
1396  else if (title=="N_{tracks, trig jet}")
1397  contents[n] = fDijet.trigJet->GetNumberOfTracks();
1398  else if (title=="N_{tracks, ass jet}")
1399  contents[n] = fDijet.assJet->GetNumberOfTracks();
1400  else
1401  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1402  }
1403  histJetObservables->Fill(contents);
1404 }
1405 
1410 {
1411  Double_t contents[30]={0};
1412  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1413  if (!histJetObservables) return;
1414  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1415  TString title(histJetObservables->GetAxis(n)->GetTitle());
1416  if (title=="A_{J}")
1417  contents[n] = fDijet.AJ;
1418  else if (title=="#Delta#phi")
1419  contents[n] = deltaPhi;
1420  else if (title=="#it{p}_{T,particle} (GeV/#it{c})")
1421  contents[n] = trackPt;
1422  else if (title=="#it{p}_{T}#parallel (GeV/#it{c})")
1423  contents[n] = balancePt;
1424  else
1425  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1426  }
1427  histJetObservables->Fill(contents);
1428 }
1429 
1434 {
1435  // Matching efficiency histogram
1436  TString histname = "GeometricalMatchingEfficiency";
1437 
1438  // If we have a matching di-jet, fill the geometrical matching histogram
1439  if (fMatchingDijet.assJet) {
1440  fHistManager.FillTH1(histname.Data(), 1);
1441 
1444  Bool_t isSwitched = trigDeltaR > fMatchingJetR;
1445 
1446  TString thnname = "GeometricalMatching";
1447  Double_t contents[30]={0};
1448  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(thnname.Data()));
1449  if (!histJetObservables) return;
1450  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1451  TString title(histJetObservables->GetAxis(n)->GetTitle());
1452  if (title=="Centrality (%)")
1453  contents[n] = fCent;
1454  else if (title=="isSwitched")
1455  contents[n] = isSwitched;
1456  else if (title=="#DeltaR_{trig}")
1457  contents[n] = trigDeltaR;
1458  else if (title=="#DeltaR_{ass}")
1459  contents[n] = assDeltaR;
1460  else if (title=="trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1461  contents[n] = fMatchingDijet.trigJetPt - fDijet.trigJetPt;
1462  else if (title=="ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1463  contents[n] = fMatchingDijet.assJetPt - fDijet.assJetPt;
1464  else if (title=="A_{J} low-threshold")
1465  contents[n] = fMatchingDijet.AJ;
1466  else if (title=="A_{J} hard-core")
1467  contents[n] = fDijet.AJ;
1468  else
1469  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1470  }
1471  histJetObservables->Fill(contents);
1472  }
1473  else {
1474  fHistManager.FillTH1(histname.Data(), 0.);
1475  }
1476 }
1477 
AliEmcalJet * GetLeadingJet(AliJetContainer *jetCont)
Double_t Area() const
Definition: AliEmcalJet.h:123
Double_t GetRhoVal() const
Int_t fNPhiBins
Number of phi bins in DCal region (for background/correction)
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
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)
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
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
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
AliRhoParameter * GetRhoParameter()
Bool_t fComputeMBDownscaling
Set whether to compute and plot MB downscaling factors.
Double_t Pt() const
Definition: AliEmcalJet.h:102
virtual Bool_t IsEventSelected()
Performing event selection.
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
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
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.
void LoadBackgroundScalingHistogram(const char *path="alien:///alice/cern.ch/user/j/jmulliga/BackgroundScalingWeights.root", const char *name="hBackgroundScalingWeights")
const AliJetIterableContainer all() const