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