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