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