AliPhysics  fe039ad (fe039ad)
 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 "AliEMCALTriggerPatchInfo.h"
42 
44 
48 
54  fDeltaPhiMin(0),
55  fMinTrigJetPt(0),
56  fMinAssJetPt(0),
57  fDijetLeadingHadronPt(0),
58  fMatchingJetR(0.2),
59  fTrackConstituentThreshold(0),
60  fClusterConstituentThreshold(0),
61  fNEtaBins(40),
62  fNPhiBins(200),
63  fBackgroundScalingWeights(0),
64  fGapJetScalingWeights(0),
65  fPlotDijetCandHistograms(kFALSE),
66  fPlotDijetImbalanceHistograms(kFALSE),
67  fComputeBackground(kFALSE),
68  fDoMomentumBalance(kFALSE),
69  fDoGeometricalMatching(kFALSE),
70  fLoadBackgroundScalingWeights(kTRUE),
71  fComputeMBDownscaling(kFALSE),
72  fDoTriggerSimulation(kFALSE),
73  fMaxPt(200),
74  fNCentHistBins(0),
75  fCentHistBins(0),
76  fNPtHistBins(0),
77  fPtHistBins(0),
78  fUseAliEventCuts(kTRUE),
79  fEventCuts(0),
80  fEventCutList(0),
81  fUseManualEventCuts(kFALSE),
82  fMBUpscaleFactor(1.),
83  fMedianEMCal(0),
84  fMedianDCal(0),
85  fkEMCEJE(kFALSE),
86  fEmbeddingQA(),
87  fHistManager()
88 {
92 }
93 
100  AliAnalysisTaskEmcalJet(name, kTRUE),
101  fDeltaPhiMin(0),
102  fMinTrigJetPt(0),
103  fMinAssJetPt(0),
104  fDijetLeadingHadronPt(0),
105  fMatchingJetR(0.2),
106  fTrackConstituentThreshold(0),
107  fClusterConstituentThreshold(0),
108  fNEtaBins(40),
109  fNPhiBins(200),
110  fBackgroundScalingWeights(0),
111  fGapJetScalingWeights(0),
112  fPlotDijetCandHistograms(kFALSE),
113  fPlotDijetImbalanceHistograms(kFALSE),
114  fComputeBackground(kFALSE),
115  fDoMomentumBalance(kFALSE),
116  fDoGeometricalMatching(kFALSE),
117  fLoadBackgroundScalingWeights(kTRUE),
118  fComputeMBDownscaling(kFALSE),
119  fDoTriggerSimulation(kFALSE),
120  fMaxPt(200),
121  fNCentHistBins(0),
122  fCentHistBins(0),
123  fNPtHistBins(0),
124  fPtHistBins(0),
125  fUseAliEventCuts(kTRUE),
126  fEventCuts(0),
127  fEventCutList(0),
128  fUseManualEventCuts(kFALSE),
129  fMBUpscaleFactor(1.),
130  fMedianEMCal(0),
131  fMedianDCal(0),
132  fkEMCEJE(kFALSE),
133  fEmbeddingQA(),
134  fHistManager(name)
135 {
137  Dijet_t fDijet;
139 }
140 
145 {
146 }
147 
152 {
153  fNCentHistBins = 4;
155  fCentHistBins[0] = 0;
156  fCentHistBins[1] = 10;
157  fCentHistBins[2] = 30;
158  fCentHistBins[3] = 50;
159  fCentHistBins[4] = 90;
160 
161  fNPtHistBins = 82;
164  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
165  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
166  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
167  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
168  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
169  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
170 }
171 
177 {
179 
180  if (fComputeBackground) {
182  }
185  }
188  }
189  if (fDoMomentumBalance) {
191  }
194  }
195  if (fDoTriggerSimulation) {
197  }
198 
199  TIter next(fHistManager.GetListOfHistograms());
200  TObject* obj = 0;
201  while ((obj = next())) {
202  fOutput->Add(obj);
203  }
204 
205  // Intialize AliEventCuts
206  if (fUseAliEventCuts) {
207  fEventCutList = new TList();
208  fEventCutList ->SetOwner();
209  fEventCutList ->SetName("EventCutOutput");
210 
211  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
212  if(fUseManualEventCuts==1)
213  {
214  fEventCuts.SetManualMode();
215  // Configure manual settings here
216  // ...
217  }
218  fEventCuts.AddQAplotsToList(fEventCutList);
219  fOutput->Add(fEventCutList);
220  }
221 
222  // Load eta-phi background scale factors from histogram on AliEn
225  }
226 
227  // Initialize embedding QA
229  if (embeddingHelper) {
230  bool res = fEmbeddingQA.Initialize();
231  if (res) {
233  }
234  }
235 
236  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
237 }
238 
239 /*
240  * This function allocates background subtraction histograms, if enabled.
241  * A set of histograms is allocated per each jet container.
242  */
244 {
245  TString histname;
246  TString title;
247 
248  AliJetContainer* jets = 0;
249  TIter nextJetColl(&fJetCollArray);
250  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
251 
252  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
253  title = histname + ";Centrality;Scale factor;counts";
254  fHistManager.CreateTH2(histname.Data(), title.Data(), 50, 0, 100, 100, 0, 5);
255 
256  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jets->GetArrayName().Data());
257  title = histname + ";Centrality (%);#delta#it{p}_{T} (GeV/#it{c});counts";
258  fHistManager.CreateTH2(histname.Data(), title.Data(), 10, 0, 100, 400, -50, 150);
259 
260  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jets->GetArrayName().Data());
261  title = histname + ";#eta;#phi;Centrality;Scale factor;";
262  Int_t nbins[4] = {fNEtaBins, fNPhiBins, 50, 400};
263  Double_t min[4] = {-0.5,1., 0, 0};
264  Double_t max[4] = {0.5,6., 100, 20};
265  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins, min, max);
266 
267  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jets->GetArrayName().Data());
268  title = histname + ";#eta;#phi;Centrality;#delta#it{p}_{T} (GeV/#it{c})";
269  Int_t nbinsDpT[4] = {fNEtaBins, fNPhiBins, 10, 400};
270  Double_t minDpT[4] = {-0.5,1., 0, -50};
271  Double_t maxDpT[4] = {0.5,6., 100, 150};
272  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbinsDpT, minDpT, maxDpT);
273 
274  }
275 }
276 
277 /*
278  * This function allocates the histograms for dijet candidates, i.e. dijet pairs with the leading jet
279  * passing the trigger condition, but no condition on the subleading jet. In particular this histogram is
280  * designed to study the kinematic selection of dijets.
281  */
283 {
284  // Allocate dijet THnSparse
285  AliJetContainer* jets = 0;
286  TIter nextJetColl(&fJetCollArray);
287  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
288  TString axisTitle[30]= {""};
289  Int_t nbins[30] = {0};
290  Double_t min[30] = {0.};
291  Double_t max[30] = {0.};
292  Double_t *binEdges[20] = {0};
293  Int_t dim = 0;
294 
295  if (fForceBeamType != kpp) {
296  axisTitle[dim] = "Centrality (%)";
297  nbins[dim] = fNCentHistBins;
298  binEdges[dim] = fCentHistBins;
299  min[dim] = fCentHistBins[0];
300  max[dim] = fCentHistBins[fNCentHistBins];
301  dim++;
302  }
303 
304  axisTitle[dim] = "LeadingHadronRequired";
305  nbins[dim] = 2;
306  min[dim] = -0.5;
307  max[dim] = 1.5;
308  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
309  dim++;
310 
311  axisTitle[dim] = "#it{p}_{T,trig jet} (GeV/#it{c})";
312  nbins[dim] = TMath::CeilNint(fMaxPt/2);
313  min[dim] = 0;
314  max[dim] = fMaxPt;
315  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
316  dim++;
317 
318  axisTitle[dim] = "#it{p}_{T,ass jet} (GeV/#it{c})";
319  nbins[dim] = TMath::CeilNint(fMaxPt/2);
320  min[dim] = 0;
321  max[dim] = fMaxPt;
322  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
323  dim++;
324 
325  axisTitle[dim] = "#phi_{trig jet}";
326  nbins[dim] = TMath::CeilNint(fMaxPt/2);
327  min[dim] = 0;
328  max[dim] = TMath::TwoPi();
329  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
330  dim++;
331 
332  axisTitle[dim] = "#phi_{ass jet}";
333  nbins[dim] = TMath::CeilNint(fMaxPt/2);
334  min[dim] = 0;
335  max[dim] = TMath::TwoPi();
336  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
337  dim++;
338 
339  axisTitle[dim] = "#eta_{trig jet}";
340  nbins[dim] = TMath::CeilNint(fMaxPt/2);
341  min[dim] = -1;
342  max[dim] = 1;
343  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
344  dim++;
345 
346  axisTitle[dim] = "#eta_{ass jet}";
347  nbins[dim] = TMath::CeilNint(fMaxPt/2);
348  min[dim] = -1;
349  max[dim] = 1;
350  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
351  dim++;
352 
353  TString thnname = TString::Format("%s/DijetCandObservables", jets->GetArrayName().Data());
354  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
355  for (Int_t i = 0; i < dim; i++) {
356  hn->GetAxis(i)->SetTitle(axisTitle[i]);
357  hn->SetBinEdges(i, binEdges[i]);
358  }
359  }
360 }
361 
362 /*
363  * This function allocates the histograms for accepted dijets.
364  * The purpose is to study in detail the imbalance properties of the accepted dijets.
365  */
367 {
368  AliJetContainer* jets = 0;
369  TIter nextJetColl(&fJetCollArray);
370  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
371 
372  // Allocate dijet imbalance THnSparse, unless doing trigger simulation
373  if (!fDoTriggerSimulation) {
374  TString axisTitle[30]= {""};
375  Int_t nbins[30] = {0};
376  Double_t min[30] = {0.};
377  Double_t max[30] = {0.};
378  Double_t *binEdges[20] = {0};
379  Int_t dim = 0;
380 
381  if (fForceBeamType != kpp) {
382  axisTitle[dim] = "Centrality (%)";
383  nbins[dim] = fNCentHistBins;
384  binEdges[dim] = fCentHistBins;
385  min[dim] = fCentHistBins[0];
386  max[dim] = fCentHistBins[fNCentHistBins];
387  dim++;
388  }
389 
390  axisTitle[dim] = "#Delta#phi";
391  nbins[dim] = 100;
392  min[dim] = 0;
393  max[dim] = 4;
394  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
395  dim++;
396 
397  axisTitle[dim] = "#Delta#eta";
398  nbins[dim] = 100;
399  min[dim] = -2;
400  max[dim] = 2;
401  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
402  dim++;
403 
404  axisTitle[dim] = "A_{J}";
405  nbins[dim] = 100;
406  min[dim] = 0;
407  max[dim] = 1;
408  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
409  dim++;
410 
411  axisTitle[dim] = "x_{J}";
412  nbins[dim] = 100;
413  min[dim] = 0;
414  max[dim] = 1;
415  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
416  dim++;
417 
418  axisTitle[dim] = "k_{Ty} (GeV)";
419  nbins[dim] = 100;
420  min[dim] = 0;
421  max[dim] = 100;
422  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
423  dim++;
424 
425  axisTitle[dim] = "N_{tracks, trig jet}";
426  nbins[dim] = fMaxPt/5;
427  min[dim] = 0;
428  max[dim] = 100;
429  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
430  dim++;
431 
432  axisTitle[dim] = "N_{tracks, ass jet}";
433  nbins[dim] = fMaxPt/5;
434  min[dim] = 0;
435  max[dim] = 100;
436  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
437  dim++;
438 
439  TString thnname = TString::Format("%s/DijetImbalanceObservables", jets->GetArrayName().Data());
440  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
441  for (Int_t i = 0; i < dim; i++) {
442  hn->GetAxis(i)->SetTitle(axisTitle[i]);
443  hn->SetBinEdges(i, binEdges[i]);
444  }
445  }
446 
447  // Now, allocate 2d pt spectrum, upscaled according to MB downscaling factors (for trigger efficiency)
448  Int_t nPtBins = TMath::CeilNint(fMaxPt/2);
449 
450  // (Centrality, pT1, pT2) (upscaled)
451  TString histname = TString::Format("%s/DijetJetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
452  TString title = histname + ";Centrality (%);#it{p}_{T,1}^{corr} (GeV/#it{c});#it{p}_{T,2}^{corr} (GeV/#it{c})";
453  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, nPtBins, 0, fMaxPt, nPtBins, 0, fMaxPt, "s");
454 
455  // Now, plot jet histograms for the leading jet within the dijet (for comparison to single jets, and triggered to MB)
456 
457  // (Centrality, pT, NEF, calo type)
458  histname = TString::Format("%s/DijetJetHistograms/hNEFVsPt", jets->GetArrayName().Data());
459  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});NEF;type";
460  Int_t nbins1[4] = {20, nPtBins, 50, 3};
461  Double_t min1[4] = {0, 0, 0, -0.5};
462  Double_t max1[4] = {100, fMaxPt, 1., 2.5};
463  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins1, min1, max1);
464 
465  // (Centrality, pT, z-leading (charged), calo type)
466  histname = TString::Format("%s/DijetJetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
467  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z}_{leading};type";
468  Int_t nbins2[4] = {20, nPtBins, 50, 3};
469  Double_t min2[4] = {0, 0, 0, -0.5};
470  Double_t max2[4] = {100, fMaxPt, 1., 2.5};
471  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins2, min2, max2);
472 
473  // (Centrality, pT, z (charged), calo type)
474  histname = TString::Format("%s/DijetJetHistograms/hZVsPt", jets->GetArrayName().Data());
475  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});#it{z};type";
476  Int_t nbins3[4] = {20, nPtBins, 50, 3};
477  Double_t min3[4] = {0, 0, 0, -0.5};
478  Double_t max3[4] = {100, fMaxPt, 1., 2.5};
479  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins3, min3, max3);
480 
481  // (Centrality, pT, Nconst, calo type)
482  histname = TString::Format("%s/DijetJetHistograms/hNConstVsPt", jets->GetArrayName().Data());
483  title = histname + ";Centrality (%);#it{p}_{T}^{corr} (GeV/#it{c});No. of constituents;type";
484  Int_t nbins4[4] = {20, nPtBins, 50, 3};
485  Double_t min4[4] = {0, 0, 0, -0.5};
486  Double_t max4[4] = {100, fMaxPt, fMaxPt, 2.5};
487  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbins4, min4, max4);
488  }
489 }
490 
491 /*
492  * This function allocates the histograms for the momentum balance study.
493  */
495 {
496  AliJetContainer* jets = 0;
497  TIter nextJetColl(&fJetCollArray);
498  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
499 
500  // Allocate THnSparse
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  axisTitle[dim] = "A_{J}";
509  nbins[dim] = 100;
510  min[dim] = 0;
511  max[dim] = 1;
512  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
513  dim++;
514 
515  axisTitle[dim] = "#Delta#phi";
516  nbins[dim] = 100;
517  min[dim] = -4;
518  max[dim] = 4;
519  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
520  dim++;
521 
522  axisTitle[dim] = "#it{p}_{T,particle} (GeV/#it{c})";
523  nbins[dim] = 9;
524  Double_t* pTParticleBins = new Double_t[nbins[dim]+1];
525  GenerateFixedBinArray(1, 0.15, 0.3, pTParticleBins);
526  GenerateFixedBinArray(1, 0.3, 0.5, pTParticleBins+1);
527  GenerateFixedBinArray(1, 0.5, 1, pTParticleBins+2);
528  GenerateFixedBinArray(2, 1, 5, pTParticleBins+3);
529  GenerateFixedBinArray(3, 5, 20, pTParticleBins+5);
530  GenerateFixedBinArray(1, 20, 150, pTParticleBins+8);
531  min[dim] = 0;
532  max[dim] = pTParticleBins[nbins[dim]];
533  binEdges[dim] = pTParticleBins;
534  dim++;
535 
536  axisTitle[dim] = "#it{p}_{T}#parallel (GeV/#it{c})";
537  nbins[dim] = 80;
538  Double_t* pTParallelBins = new Double_t[nbins[dim]+1];
539  GenerateFixedBinArray(20, 0, 2, pTParallelBins);
540  GenerateFixedBinArray(16, 2, 10, pTParallelBins+20);
541  GenerateFixedBinArray(10, 10, 20, pTParallelBins+36);
542  GenerateFixedBinArray(10, 20, 40, pTParallelBins+46);
543  GenerateFixedBinArray(24, 40, 150, pTParallelBins+56);
544  min[dim] = 0;
545  max[dim] = pTParallelBins[nbins[dim]];
546  binEdges[dim] = pTParallelBins;
547  dim++;
548 
549  TString thnname = TString::Format("%s/MomentumBalance", jets->GetArrayName().Data());
550  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
551  for (Int_t i = 0; i < dim; i++) {
552  hn->GetAxis(i)->SetTitle(axisTitle[i]);
553  hn->SetBinEdges(i, binEdges[i]);
554  }
555  }
556 }
557 
558 /*
559  * This function allocates the histograms for the constituent dijet study.
560  */
562 {
563  // Allocate geometrical matching THnSparse
564  TString axisTitle[30]= {""};
565  Int_t nbins[30] = {0};
566  Double_t min[30] = {0.};
567  Double_t max[30] = {0.};
568  Double_t *binEdges[20] = {0};
569  Int_t dim = 0;
570 
571  if (fForceBeamType != kpp) {
572  axisTitle[dim] = "Centrality (%)";
573  nbins[dim] = fNCentHistBins;
574  binEdges[dim] = fCentHistBins;
575  min[dim] = fCentHistBins[0];
576  max[dim] = fCentHistBins[fNCentHistBins];
577  dim++;
578  }
579 
580  axisTitle[dim] = "isSwitched";
581  nbins[dim] = 2;
582  min[dim] = -0.5;
583  max[dim] = 1.5;
584  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
585  dim++;
586 
587  axisTitle[dim] = "#DeltaR_{trig}";
588  nbins[dim] = 50;
589  min[dim] = 0;
590  max[dim] = 0.5;
591  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
592  dim++;
593 
594  axisTitle[dim] = "#DeltaR_{ass}";
595  nbins[dim] = 50;
596  min[dim] = 0;
597  max[dim] = 0.5;
598  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
599  dim++;
600 
601  axisTitle[dim] = "trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
602  nbins[dim] = 100;
603  min[dim] = -50;
604  max[dim] = 50;
605  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
606  dim++;
607 
608  axisTitle[dim] = "ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}";
609  nbins[dim] = 100;
610  min[dim] = -50;
611  max[dim] = 50;
612  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
613  dim++;
614 
615  axisTitle[dim] = "A_{J} low-threshold";
616  nbins[dim] = 100;
617  min[dim] = 0;
618  max[dim] = 1;
619  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
620  dim++;
621 
622  axisTitle[dim] = "A_{J} hard-core";
623  nbins[dim] = 100;
624  min[dim] = 0;
625  max[dim] = 1;
626  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
627  dim++;
628 
629  TString thnname = "GeometricalMatching";
630  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
631  for (Int_t i = 0; i < dim; i++) {
632  hn->GetAxis(i)->SetTitle(axisTitle[i]);
633  hn->SetBinEdges(i, binEdges[i]);
634  }
635 
636  // Allocate other histograms
637  TString histname;
638  TString title;
639  histname = "GeometricalMatchingEfficiency";
640  title = histname + ";isMatched;counts";
641  fHistManager.CreateTH1(histname.Data(), title.Data(), 2, -0.5, 1.5);
642 }
643 
644 /*
645  * This function allocates the histograms for single jets, when the "simulated" trigger has been fired.
646  * A set of histograms is allocated per each jet container.
647  */
649 {
650  TString histname;
651  TString title;
652 
653  //----------------------------------------------
654  // Trigger patch histograms
655 
656  // Median patch energy vs. centrality, for dijets
657  histname = "TriggerSimHistograms/hMedPatchDijet";
658  title = histname + ";Centrality (%);#it{p}_{T,trig}^{corr} (GeV/#it{c});#it{E}_{patch,med} (GeV);type";
659  Int_t nbinsD[4] = {50, 40, 100, 2};
660  Double_t minD[4] = {0, 0, 0, -0.5};
661  Double_t maxD[4] = {100, 200, 50, 1.5};
662  fHistManager.CreateTHnSparse(histname.Data(), title.Data(), 4, nbinsD, minD, maxD);
663 
664 }
665 
669 void AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram(const char* path, const char* name1, const char* name2)
670 {
671 
672  TString fname(path);
673  if (fname.BeginsWith("alien://")) {
674  TGrid::Connect("alien://");
675  }
676 
677  TFile* file = TFile::Open(path);
678 
679  if (!file || file->IsZombie()) {
680  ::Error("AliAnalysisTaskEmcalDijetImbalance", "Could not open background scaling histogram");
681  return;
682  }
683 
684  // Open background scale factor histogram
685  TH2D* h1 = dynamic_cast<TH2D*>(file->Get(name1));
686 
687  if (h1) {
688  ::Info("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s loaded from file %s.", name1, path);
689  }
690  else {
691  ::Error("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Background histogram %s not found in file %s.", name1, path);
692  return;
693  }
694 
695  fBackgroundScalingWeights = static_cast<TH2D*>(h1->Clone());
696 
697  // Open jet pT scale factor histogram
698  TH2D* h2 = dynamic_cast<TH2D*>(file->Get(name2));
699 
700  if (h2) {
701  ::Info("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Jet pT scaling histogram %s loaded from file %s.", name2, path);
702  }
703  else {
704  ::Error("AliAnalysisTaskEmcalDijetImbalance::LoadBackgroundScalingHistogram", "Jet pT scaling histogram %s not found in file %s.", name2, path);
705  return;
706  }
707 
708  fGapJetScalingWeights = static_cast<TH2D*>(h2->Clone());
709 
710  file->Close();
711  delete file;
712 
713 }
714 
720 {
721 
722  // Configure base class to set fTriggerPatchInfo to array of trigger patches, each event
723  // (Need to call this before base class ExecOnce)
724  if (fDoTriggerSimulation) {
725  this->SetCaloTriggerPatchInfoName("EmcalTriggers");
726  }
727 
729 
730  fNeedEmcalGeom = kTRUE;
731 
732  // Check if trigger patches are loaded
733  if (fDoTriggerSimulation) {
734  if (fTriggerPatchInfo) {
735  TString objname(fTriggerPatchInfo->GetClass()->GetName());
736  TClass cls(objname);
737  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
738  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
739  GetName(), cls.GetName(), "EmcalTriggers"));
740  fTriggerPatchInfo = 0;
741  }
742  }
743  if (!fTriggerPatchInfo) {
744  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
745  return;
746  }
747  }
748 
749  AliInfo(Form("Trigger jet threshold = %f, Associated jet threshold = %f", fMinTrigJetPt, fMinAssJetPt));
750  AliInfo(Form("Leading hadron threshold (for dijet leading jet): %f GeV", fDijetLeadingHadronPt));
751  AliInfo(Form("Momentum balance study: %d", fDoMomentumBalance));
752  AliInfo(Form("Geometrical matching study: %d", fDoGeometricalMatching));
753 
754 }
755 
760 
761  // Get the downscaling factors for MB triggers (to be used to calculate trigger efficiency)
762 
763  if (fComputeMBDownscaling) {
764 
765  // Get instance of the downscale factor helper class
767  downscaleOCDB->SetRun(InputEvent()->GetRunNumber());
768 
769  // There are two possible min bias triggers for LHC15o
770  TString triggerNameMB1 = "CINT7-B-NOPF-CENT";
771  TString triggerNameMB2 = "CV0L7-B-NOPF-CENT";
772  TString triggerNameJE = "CINT7EJ1-B-NOPF-CENTNOPMD";
773 
774  // Get the downscale factor for whichever MB trigger exists in the given run
775  std::vector<TString> runtriggers = downscaleOCDB->GetTriggerClasses();
776  Double_t downscalefactor;
777  for (auto i : runtriggers) {
778  if (i.EqualTo(triggerNameMB1) || i.EqualTo(triggerNameMB2)) {
779  downscalefactor = downscaleOCDB->GetDownscaleFactorForTriggerClass(i.Data());
780  break;
781  }
782  }
783 
784  // Store the inverse of the downscale factor, used later to weight the pT spectrum
785  fMBUpscaleFactor = 1/downscalefactor;
786 
787  }
788 
789 }
790 
795 {
796  if (fUseAliEventCuts) {
797  if (!fEventCuts.AcceptEvent(InputEvent()))
798  {
799  PostData(1, fOutput);
800  return kFALSE;
801  }
802  }
803  else {
805  }
806  return kTRUE;
807 }
808 
817 {
818  TString histname;
819  AliJetContainer* jetCont = 0;
820  TIter next(&fJetCollArray);
821  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
822  TString jetContName = jetCont->GetName();
823  if (jetContName.Contains("HardCore")) continue;
824 
825  //-----------------------------------------------------------------------------
826  // Find the leading di-jet candidate in each event, and if it satisfies the
827  // trig jet pT threshold, then fill di-jet candidate histogram (regardless of ass jet).
828  // The idea is to study the kinematic selections in post-processing.
829 
830  // Loop over leading hadron cut or not
831  for (Int_t leadingHadronCutType=0; leadingHadronCutType<2; leadingHadronCutType++) {
832 
833  // Find the dijet candidate of the event and store its info in struct fDijet
834  FindDijet(jetCont, leadingHadronCutType);
835 
836  // If we find a dijet candidate (i.e. acceptable trig jet; ass jet accepted or not), fill the di-jet candidate histogram
838  FillDijetCandHistograms(jetCont);
839  }
840 
841  }
842 
843  //---------------------------------------------------------------------------------------------------
844  // Now, study the accepted dijet selection -- specified by the trig/ass jet pT conditions
845 
846  // Find the dijet candidate of the event and store its info in struct fDijet
847  FindDijet(jetCont, 0);
848 
849  // If we find an accepted dijet, fill the dijet imbalance histogram
852  }
853  // If we find an accepted dijet, perform momentum-balance study (if requested)
855  histname = TString::Format("%s/MomentumBalance", jetCont->GetArrayName().Data());
856  DoMomentumBalance(histname);
857  }
858 
859  //---------------------------------------------------------------------------
860  // Do a simple trigger simulation (if requested)
861  if (fDoTriggerSimulation) {
863  }
864 
865  }
866 
867  //---------------------------------------------------------------------------
868  // Do the constituent threshold and geometrical matching study (if requested)
871  }
872 
873  // Only fill the embedding qa plots if:
874  // - We are using the embedding helper
875  // - The class has been initialized
876  // - Both jet collections are available
877  if (fEmbeddingQA.IsInitialized()) {
879  }
880 
881  return kTRUE;
882 }
883 
892 {
893  fDijet.clear();
894  fDijet.leadingHadronCutType = leadingHadronCutBin;
895 
896  // Get trigger jet
897  AliEmcalJet* trigJet = GetLeadingJet(jetCont);
898  if(!trigJet) return;
899 
900  // Skip the event if the leading jet doesn't satisfy the pT threshold
901  Double_t trigJetPt = GetJetPt(jetCont, trigJet);
902  if ( trigJetPt < fMinTrigJetPt ) return;
903 
904  // Skip the event if the leading jet doesn't satisfy the leading hadron threshold
905  if (jetCont->GetLeadingHadronPt(trigJet) < fDijetLeadingHadronPt*leadingHadronCutBin) return;
906 
907  // Fill the dijet struct
908  fDijet.trigJet = trigJet;
909  fDijet.trigJetPt = trigJetPt;
910  fDijet.trigJetEta = trigJet->Eta();
911  fDijet.trigJetPhi = trigJet->Phi();
912 
913  // Find the subleading jet in the opposite hemisphere
914  AliEmcalJet *assJet = 0;
915  for(auto assJetCand : jetCont->accepted()) {
916  if (!assJetCand) continue;
917  Double_t assJetCandPt = GetJetPt(jetCont, assJetCand);
918  if ( TMath::Abs(trigJet->Phi() - assJetCand->Phi()) < fDeltaPhiMin ) continue;
919  if (assJet) {
920  Double_t assJetPt = GetJetPt(jetCont, assJet);
921  if ( assJetCandPt < assJetPt ) continue;
922  }
923  assJet = assJetCand;
924  }
925  if (!assJet) return;
926 
927  // Fill the dijet struct
928  fDijet.assJet = assJet;
929  fDijet.assJetPt = GetJetPt(jetCont, assJet);
930  fDijet.assJetPhi = assJet->Phi();
931  fDijet.assJetEta = assJet->Eta();
933 
934  fDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
935  fDijet.deltaEta = trigJet->Eta() - assJet->Eta();
938  fDijet.kTy = TMath::Abs( fDijet.trigJetPt * TMath::Sin(fDijet.deltaPhi) );
939 }
940 
946 {
947 
948  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
949 
950  AliVTrack* track;
951  for (auto trackIterator : trackCont->accepted_momentum() ) {
952 
953  track = trackIterator.second;
954 
955  // Compute the delta phi between the track and its nearest jet (of the two jets in the dijet),
956  // as well as its pT-parallel projection onto the nearest jet's axis.
957 
958  Double_t trackPt = track->Pt();
959  Double_t trackPhi = track->Phi();
960  Double_t trigJetPhi = fDijet.trigJet->Phi();
961  Double_t assJetPhi = fDijet.assJet->Phi();
962 
963  Double_t deltaPhiTrigJet = TMath::Abs(trackPhi - trigJetPhi);
964  Double_t deltaPhiAssJet = TMath::Abs(trackPhi - assJetPhi);
965  Bool_t isNearside = deltaPhiTrigJet < deltaPhiAssJet;
966 
967  Double_t deltaPhi;
968  Double_t balancePt;
969  if (isNearside) {
970  deltaPhi = trackPhi - trigJetPhi;
971  balancePt = trackPt * TMath::Cos(deltaPhi);
972  }
973  else {
974  deltaPhi = trackPhi - assJetPhi;
975  balancePt = -trackPt * TMath::Cos(deltaPhi);
976  }
977 
978  FillMomentumBalanceHistograms(histname, deltaPhi, trackPt, balancePt);
979 
980  }
981 }
982 
987 {
988  // Get jet container with minimum constituent pT,E thresholds
989  TString jetContAllName = Form("Jet_AKTFullR0%d0_tracks_pT0150_caloClusters_E0300_pt_scheme", (int) (fMatchingJetR*10) );
990  AliJetContainer* jetContAll = GetJetContainer(jetContAllName.Data());
991 
992  // Get jet container with X GeV constituent pT,E thresholds
993  Int_t trackThreshold = (int) (fTrackConstituentThreshold*1000); // in MeV
994  Int_t clusThreshold = (int) (fClusterConstituentThreshold*1000); // in MeV
995  TString jetContHardCoreName = Form("JetHardCore_AKTFullR0%d0_tracks_pT%d_caloClusters_E%d_pt_scheme", (int) (fMatchingJetR*10), trackThreshold, clusThreshold);
996  AliJetContainer* jetContHardCore = GetJetContainer(jetContHardCoreName.Data());
997 
998  // Find the di-jet in the hard-core jet sample, then find the matching di-jet and fill histograms
999  FindDijet(jetContHardCore, 0);
1000  if (fDijet.isAccepted) {
1001  FindMatchingDijet(jetContAll);
1003  }
1004 }
1005 
1011 {
1012  TString histname;
1013 
1014  // Check if trigger patches are loaded
1015  if (fTriggerPatchInfo) {
1016  TString objname(fTriggerPatchInfo->GetClass()->GetName());
1017  TClass cls(objname);
1018  if (!cls.InheritsFrom("AliEMCALTriggerPatchInfo")) {
1019  AliError(Form("%s: Objects of type %s in %s are not inherited from AliEMCALTriggerPatchInfo!",
1020  GetName(), cls.GetName(), "EmcalTriggers"));
1021  fTriggerPatchInfo = 0;
1022  }
1023  }
1024  if (!fTriggerPatchInfo) {
1025  AliError(Form("%s: Unable to get trigger patch container with name %s. Aborting", GetName(), "EmcalTriggers"));
1026  return;
1027  }
1028 
1029  // Compute patches in EMCal, DCal (I want offline simple trigger patch, i.e. patch calculated using FEE energy)
1030  std::vector<Double_t> vecEMCal;
1031  std::vector<Double_t> vecDCal;
1032  for(auto p : *fTriggerPatchInfo){
1033  AliEMCALTriggerPatchInfo *recpatch = static_cast<AliEMCALTriggerPatchInfo *>(p);
1034  if (recpatch) {
1035 
1036  if(!recpatch->IsJetHighSimple()) continue;
1037 
1038  if (recpatch->IsEMCal()) {
1039  vecEMCal.push_back(recpatch->GetPatchE());
1040  } else {
1041  vecDCal.push_back(recpatch->GetPatchE());
1042  }
1043 
1044  }
1045  }
1046 
1047  // Compute the median in each calorimeter
1048  const Int_t nBkgPatchesEMCal = vecEMCal.size(); // 6*8;
1049  const Int_t nBkgPatchesDCal = vecDCal.size(); // 4*5;
1050  fMedianEMCal = TMath::Median(nBkgPatchesEMCal, &vecEMCal[0]); // point to array used internally by vector
1051  fMedianDCal = TMath::Median(nBkgPatchesDCal, &vecDCal[0]);
1052 
1053  // Median patch energy vs. pT for dijets
1054  if (fDijet.isAccepted) {
1055  histname = "TriggerSimHistograms/hMedPatchDijet";
1057  fHistManager.FillTHnSparse(histname, x);
1059  fHistManager.FillTHnSparse(histname, y);
1060  }
1061 
1062 }
1063 
1069 {
1071 
1072  // Loop over jets and find leading jet within R of fDijet.trigJet
1073  AliEmcalJet *matchingTrigJet = 0;
1074  for(auto matchingTrigJetCand : jetCont->accepted()) {
1075  if (!matchingTrigJetCand) continue;
1076  if ( GetDeltaR(matchingTrigJetCand, fDijet.trigJet) > fMatchingJetR ) continue;
1077  if (matchingTrigJet) {
1078  if ( GetJetPt(jetCont, matchingTrigJetCand) < GetJetPt(jetCont, matchingTrigJet) ) continue;
1079  }
1080  matchingTrigJet = matchingTrigJetCand;
1081  }
1082  if (!matchingTrigJet) return;
1083 
1084  // Loop over jets and find leading jet within R of fDijet.assJet
1085  AliEmcalJet *matchingAssJet = 0;
1086  for(auto matchingAssJetCand : jetCont->accepted()) {
1087  if (!matchingAssJetCand) continue;
1088  if ( GetDeltaR(matchingAssJetCand, fDijet.assJet) > fMatchingJetR ) continue;
1089  if (matchingAssJet) {
1090  if ( GetJetPt(jetCont, matchingAssJetCand) < GetJetPt(jetCont, matchingAssJet) ) continue;
1091  }
1092  matchingAssJet = matchingAssJetCand;
1093  }
1094 
1095  // Determine which matching jet is the leading jet (i.e. allow them to flip)
1096  if (matchingAssJet) {
1097  AliEmcalJet* trigJet = matchingTrigJet;
1098  AliEmcalJet* assJet = matchingAssJet;
1099  if ( GetJetPt(jetCont, matchingTrigJet) < GetJetPt(jetCont, matchingAssJet) ) {
1100  trigJet = matchingAssJet;
1101  assJet = matchingTrigJet;
1102  }
1103 
1104  // Fill the dijet struct
1105  fMatchingDijet.trigJet = trigJet;
1106  fMatchingDijet.trigJetPt = GetJetPt(jetCont, trigJet);
1107  fMatchingDijet.trigJetEta = trigJet->Eta();
1108  fMatchingDijet.trigJetPhi = trigJet->Phi();
1109 
1110  fMatchingDijet.assJet = assJet;
1111  fMatchingDijet.assJetPt = GetJetPt(jetCont, assJet);
1112  fMatchingDijet.assJetPhi = assJet->Phi();
1113  fMatchingDijet.assJetEta = assJet->Eta();
1115 
1116  fMatchingDijet.deltaPhi = TMath::Abs(trigJet->Phi() - assJet->Phi());
1117  fMatchingDijet.deltaEta = trigJet->Eta() - assJet->Eta();
1120  fMatchingDijet.kTy = TMath::Abs( fMatchingDijet.trigJetPt * TMath::Sin(fMatchingDijet.deltaPhi) );
1121  }
1122 }
1123 
1128 {
1129  // Loop over tracks and clusters in order to:
1130  // (1) Compute scale factor for full jets
1131  // (2) Compute delta-pT for full jets, with the random cone method
1132  // For both the scale factor and delta-pT, we compute only one histogram each for EMCal.
1133  // And then we bin in eta-phi, in order to compute and perform a corretion to account for the DCal vs. PHOS vs. gap
1134 
1135  // Define the acceptance boundaries for the TPC and EMCal/DCal/PHOS
1136  Double_t etaTPC = 0.9;
1137  Double_t etaEMCal = 0.7;
1138  //Double_t etaMinDCal = 0.22;
1139  //Double_t etaMaxPHOS = 0.13;
1140  Double_t phiMinEMCal = fGeom->GetArm1PhiMin() * TMath::DegToRad(); // 80
1141  Double_t phiMaxEMCal = fGeom->GetEMCALPhiMax() * TMath::DegToRad(); // ~188
1142  //Double_t phiMinDCal = fGeom->GetDCALPhiMin() * TMath::DegToRad(); // 260
1143  //Double_t phiMaxDCal = fGeom->GetDCALPhiMax() * TMath::DegToRad(); // ~327 (1/3 SMs start at 320)
1144  //Double_t phiMinPHOS = 250 * TMath::DegToRad();
1145  //Double_t phiMaxPHOS = 320 * TMath::DegToRad();
1146 
1147  Double_t accTPC = 2 * etaTPC * 2 * TMath::Pi();
1148  Double_t accEMCal = 2 * etaEMCal * (phiMaxEMCal - phiMinEMCal);
1149  //Double_t accDCalRegion = 2 * etaEMCal * (phiMaxDCal - phiMinDCal);
1150 
1151  // Loop over jet containers
1152  AliJetContainer* jetCont = 0;
1153  TIter nextJetColl(&fJetCollArray);
1154  while ((jetCont = static_cast<AliJetContainer*>(nextJetColl()))) {
1155 
1156  // Define fiducial acceptances, to be used to generate random cones
1157  TRandom3* r = new TRandom3(0);
1158  Double_t jetR = jetCont->GetJetRadius();
1159  Double_t accRC = TMath::Pi() * jetR * jetR;
1160  Double_t etaEMCalfid = etaEMCal - jetR;
1161  Double_t phiMinEMCalfid = phiMinEMCal + jetR;
1162  Double_t phiMaxEMCalfid = phiMaxEMCal - jetR;
1163 
1164  // Generate EMCal random cone eta-phi
1165  Double_t etaEMCalRC = r->Uniform(-etaEMCalfid, etaEMCalfid);
1166  Double_t phiEMCalRC = r->Uniform(phiMinEMCalfid, phiMaxEMCalfid);
1167 
1168  // For eta-phi correction, generate random eta, phi in each eta/phi bin, to be used as center of random cone
1169  Double_t etaDCalRC[fNEtaBins]; // array storing the RC eta values
1170  Double_t etaStep = 1./fNEtaBins;
1171  Double_t etaMin;
1172  Double_t etaMax;
1173  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1174  etaMin = -etaEMCalfid + etaBin*etaStep;
1175  etaMax = etaMin + etaStep;
1176  etaDCalRC[etaBin] = r->Uniform(etaMin, etaMax);
1177  }
1178 
1179  Double_t phiDCalRC[fNPhiBins]; // array storing the RC phi values
1180  Double_t phiStep = 5./fNPhiBins; // phi axis is [1,6] in order to have simple binning
1181  Double_t phiMin;
1182  Double_t phiMax;
1183  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1184  phiMin = 1 + phiBin*phiStep;
1185  phiMax = phiMin + phiStep;
1186  phiDCalRC[phiBin] = r->Uniform(phiMin, phiMax);
1187  }
1188 
1189  // Initialize the various sums to 0
1190  Double_t trackPtSumTPC = 0;
1191  Double_t trackPtSumEMCal = 0;
1192  Double_t trackPtSumEMCalRC = 0;
1193  Double_t clusESumEMCal = 0;
1194  Double_t clusESumEMCalRC = 0;
1195 
1196  // Define a 2D vector (initialized to 0) to store the sum of track pT, and another for cluster ET
1197  std::vector<std::vector<Double_t>> trackPtSumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1198  std::vector<std::vector<Double_t>> clusESumDCalRC(fNEtaBins, std::vector<Double_t>(fNPhiBins));
1199 
1200  // Loop over tracks. Sum the track pT:
1201  // (1) in the entire TPC, (2) in the EMCal, (3) in the EMCal random cone,
1202  // (4) in a random cone at each eta-phi
1203  AliTrackContainer* trackCont = dynamic_cast<AliTrackContainer*>(GetParticleContainer("tracks"));
1204  AliTLorentzVector track;
1205  Double_t trackEta;
1206  Double_t trackPhi;
1207  Double_t trackPt;
1208  Double_t deltaR;
1209  for (auto trackIterator : trackCont->accepted_momentum() ) {
1210 
1211  track.Clear();
1212  track = trackIterator.first;
1213  trackEta = track.Eta();
1214  trackPhi = track.Phi_0_2pi();
1215  trackPt = track.Pt();
1216 
1217  // (1)
1218  if (TMath::Abs(trackEta) < etaTPC) {
1219  trackPtSumTPC += trackPt;
1220  }
1221 
1222  // (2)
1223  if (TMath::Abs(trackEta) < etaEMCal && trackPhi > phiMinEMCal && trackPhi < phiMaxEMCal) {
1224  trackPtSumEMCal += trackPt;
1225  }
1226 
1227  // (3)
1228  deltaR = GetDeltaR(&track, etaEMCalRC, phiEMCalRC);
1229  if (deltaR < jetR) {
1230  trackPtSumEMCalRC += trackPt;
1231  }
1232 
1233  // (4)
1234  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1235  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1236  deltaR = GetDeltaR(&track, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1237  if (deltaR < jetR) {
1238  trackPtSumDCalRC[etaBin][phiBin] += trackPt;
1239  }
1240  }
1241  }
1242 
1243  }
1244 
1245  // Loop over clusters. Sum the cluster ET:
1246  // (1) in the EMCal, (2) in the EMCal random cone, (3) in a random cone at each eta-phi
1248  AliTLorentzVector clus;
1249  Double_t clusEta;
1250  Double_t clusPhi;
1251  Double_t clusE;
1252  for (auto clusIterator : clusCont->accepted_momentum() ) {
1253 
1254  clus.Clear();
1255  clus = clusIterator.first;
1256  clusEta = clus.Eta();
1257  clusPhi = clus.Phi_0_2pi();
1258  clusE = clus.E();
1259 
1260  // (1)
1261  if (TMath::Abs(clusEta) < etaEMCal && clusPhi > phiMinEMCal && clusPhi < phiMaxEMCal) {
1262  clusESumEMCal += clusE;
1263  }
1264 
1265  // (2)
1266  deltaR = GetDeltaR(&clus, etaEMCalRC, phiEMCalRC);
1267  if (deltaR < jetR) {
1268  clusESumEMCalRC += clusE;
1269  }
1270 
1271  // (3)
1272  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1273  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1274  deltaR = GetDeltaR(&clus, etaDCalRC[etaBin], phiDCalRC[phiBin]);
1275  if (deltaR < jetR) {
1276  clusESumDCalRC[etaBin][phiBin] += clusE;
1277  }
1278  }
1279  }
1280 
1281  }
1282 
1283  // Compute the scale factor for EMCal, as a function of centrality
1284  Double_t numerator = (trackPtSumEMCal + clusESumEMCal) / accEMCal;
1285  Double_t denominator = trackPtSumTPC / accTPC;
1286  Double_t scaleFactor = numerator / denominator;
1287  TString histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jetCont->GetArrayName().Data());
1288  fHistManager.FillTH2(histname, fCent, scaleFactor);
1289 
1290  // Compute the scale factor in each eta-phi bin, as a function of centrality
1291  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1292  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1293  numerator = (trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin]) / accRC;
1294  scaleFactor = numerator / denominator;
1295  histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEtaPhi", jetCont->GetArrayName().Data());
1296  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, scaleFactor};
1297  fHistManager.FillTHnSparse(histname, x);
1298  }
1299  }
1300 
1301  // Compute delta pT for EMCal, as a function of centrality
1302  Double_t rho = jetCont->GetRhoVal();
1303  Double_t deltaPt = trackPtSumEMCalRC + clusESumEMCalRC - rho * TMath::Pi() * jetR * jetR;
1304  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEMCal", jetCont->GetArrayName().Data());
1305  fHistManager.FillTH2(histname, fCent, deltaPt);
1306 
1307  // Compute delta pT in each eta-phi bin, as a function of centrality
1308  Double_t sf;
1309  for (Int_t etaBin=0; etaBin < fNEtaBins; etaBin++) {
1310  for (Int_t phiBin=0; phiBin < fNPhiBins; phiBin++) {
1312  sf = fBackgroundScalingWeights->GetBinContent(fBackgroundScalingWeights->FindBin(etaDCalRC[etaBin], phiDCalRC[phiBin]));
1313  rho = sf * jetCont->GetRhoVal();
1314  }
1315  deltaPt = trackPtSumDCalRC[etaBin][phiBin] + clusESumDCalRC[etaBin][phiBin] - rho * accRC;
1316  histname = TString::Format("%s/BackgroundHistograms/hDeltaPtEtaPhi", jetCont->GetArrayName().Data());
1317  Double_t x[4] = {etaDCalRC[etaBin], phiDCalRC[phiBin], fCent, deltaPt};
1318  fHistManager.FillTHnSparse(histname, x);
1319  }
1320  }
1321 
1322  delete r;
1323 
1324  }
1325 
1326 }
1327 
1332 {
1333 
1334  // Get eta-phi dependent jet pT scale factor
1335  Double_t jetPt = jet->Pt();
1336  if (fGapJetScalingWeights) {
1337  Double_t sf = fGapJetScalingWeights->GetBinContent(fGapJetScalingWeights->FindBin(jet->Eta(), jet->Phi_0_2pi()));
1338  jetPt = jetPt * (1 + sf * jet->NEF());
1339  }
1340 
1341  // Compute pTcorr
1342  Double_t rho = jetCont->GetRhoVal();
1343  Double_t pT = jetPt - rho * jet->Area();
1344 
1345  // If hard-core jet, don't subtract background
1346  TString jetContName = jetCont->GetName();
1347  if (jetContName.Contains("HardCore")) pT = jet->Pt();
1348 
1349  return pT;
1350 }
1351 
1356 {
1357  AliEmcalJet* leadingJet = 0;
1358 
1359  if (jetCont->GetRhoParameter()) {
1360  for(auto jetCand : jetCont->accepted()) {
1361  if (!jetCand) continue;
1362  Double_t jetCandPt = GetJetPt(jetCont, jetCand);
1363  if (leadingJet) {
1364  Double_t leadingJetPt = GetJetPt(jetCont, leadingJet);
1365  if ( jetCandPt < leadingJetPt ) continue;
1366  }
1367  leadingJet = jetCand;
1368  }
1369  }
1370  else {
1371  leadingJet = jetCont->GetLeadingJet();
1372  }
1373 
1374  return leadingJet;
1375 }
1376 
1381 {
1382  Double_t deltaPhi = TMath::Abs(jet1->Phi() - jet2->Phi());
1383  Double_t deltaEta = TMath::Abs(jet1->Eta() - jet2->Eta());
1384  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1385  return deltaR;
1386 }
1387 
1392 {
1393  Double_t deltaPhi = TMath::Abs(part->Phi_0_2pi() - phiRef);
1394  Double_t deltaEta = TMath::Abs(part->Eta() - etaRef);
1395  Double_t deltaR = TMath::Sqrt( deltaPhi*deltaPhi + deltaEta*deltaEta );
1396  return deltaR;
1397 }
1398 
1403 {
1404  UInt_t jetType = jet->GetJetAcceptanceType();
1405  Double_t type = -1;
1406  if (jetType & AliEmcalJet::kEMCAL) {
1407  type = kEMCal;
1408  }
1409  else if (jetType & AliEmcalJet::kDCALonly) {
1410  type = kDCal;
1411  }
1412  else if (jetType & AliEmcalJet::kPHOS) {
1413  type = kPHOS;
1414  }
1415 
1416  return type;
1417 }
1418 
1419 
1427 {
1428 
1429  if (fComputeBackground) {
1431  }
1432 
1433  return kTRUE;
1434 }
1435 
1440 {
1441  Double_t contents[30]={0};
1442  TString histname = TString::Format("%s/DijetCandObservables", jets->GetArrayName().Data());
1443  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1444  if (!histJetObservables) return;
1445  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1446  TString title(histJetObservables->GetAxis(n)->GetTitle());
1447  if (title=="Centrality (%)")
1448  contents[n] = fCent;
1449  else if (title=="LeadingHadronRequired")
1450  contents[n] = fDijet.leadingHadronCutType;
1451  else if (title=="#it{p}_{T,trig jet} (GeV/#it{c})")
1452  contents[n] = fDijet.trigJetPt;
1453  else if (title=="#it{p}_{T,ass jet} (GeV/#it{c})")
1454  contents[n] = fDijet.assJetPt;
1455  else if (title=="#phi_{trig jet}")
1456  contents[n] = fDijet.trigJetPhi;
1457  else if (title=="#phi_{ass jet}")
1458  contents[n] = fDijet.assJetPhi;
1459  else if (title=="#eta_{trig jet}")
1460  contents[n] = fDijet.trigJetEta;
1461  else if (title=="#eta_{ass jet}")
1462  contents[n] = fDijet.assJetEta;
1463  else
1464  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1465  }
1466  histJetObservables->Fill(contents);
1467 }
1468 
1473 {
1474  // Fill the dijet imbalance histogram (unless doing trigger simulation, in which case bypass this)
1475  if (!fDoTriggerSimulation) {
1476  Double_t contents[30]={0};
1477  TString histname = TString::Format("%s/DijetImbalanceObservables", jets->GetArrayName().Data());
1478  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1479  if (!histJetObservables) return;
1480  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1481  TString title(histJetObservables->GetAxis(n)->GetTitle());
1482  if (title=="Centrality (%)")
1483  contents[n] = fCent;
1484  else if (title=="#Delta#phi")
1485  contents[n] = fDijet.deltaPhi;
1486  else if (title=="#Delta#eta")
1487  contents[n] = fDijet.deltaEta;
1488  else if (title=="A_{J}")
1489  contents[n] = fDijet.AJ;
1490  else if (title=="x_{J}")
1491  contents[n] = fDijet.xJ;
1492  else if (title=="k_{Ty} (GeV)")
1493  contents[n] = fDijet.kTy;
1494  else if (title=="N_{tracks, trig jet}")
1495  contents[n] = fDijet.trigJet->GetNumberOfTracks();
1496  else if (title=="N_{tracks, ass jet}")
1497  contents[n] = fDijet.assJet->GetNumberOfTracks();
1498  else
1499  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1500  }
1501  histJetObservables->Fill(contents);
1502  }
1503 
1504  // (Centrality, pT1, pT2) (upscaled)
1505  TString histname = TString::Format("%s/DijetJetHistograms/hPtUpscaledMB", jets->GetArrayName().Data());
1507 
1508  // Get jet acceptance type
1509  AliEmcalJet* trigJet = fDijet.trigJet;
1510  Double_t type = GetJetType(trigJet);
1511 
1512  // (Centrality, pT, NEF, calo type)
1513  histname = TString::Format("%s/DijetJetHistograms/hNEFVsPt", jets->GetArrayName().Data());
1514  Double_t x[4] = {fCent, fDijet.trigJetPt, trigJet->NEF(), type};
1515  fHistManager.FillTHnSparse(histname, x);
1516 
1517  // (Centrality, pT, z-leading (charged), calo type)
1518  histname = TString::Format("%s/DijetJetHistograms/hZLeadingVsPt", jets->GetArrayName().Data());
1519  TLorentzVector leadPart;
1520  jets->GetLeadingHadronMomentum(leadPart, trigJet);
1521  Double_t z = GetParallelFraction(leadPart.Vect(), trigJet);
1522  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
1523  Double_t y[4] = {fCent, fDijet.trigJetPt, z, type};
1524  fHistManager.FillTHnSparse(histname, y);
1525 
1526  // (Centrality, pT, z (charged), calo type)
1527  histname = TString::Format("%s/DijetJetHistograms/hZVsPt", jets->GetArrayName().Data());
1528  AliVTrack* track;
1529  for (Int_t i=0; i<trigJet->GetNumberOfTracks(); i++) {
1530  track = static_cast<AliVTrack*>(trigJet->Track(i));
1531  z = track->Pt() / TMath::Abs(fDijet.trigJetPt);
1532  Double_t y2[4] = {fCent, fDijet.trigJetPt, z, type};
1533  fHistManager.FillTHnSparse(histname, y2);
1534  }
1535 
1536  // (Centrality, pT, Nconst, calo type)
1537  histname = TString::Format("%s/DijetJetHistograms/hNConstVsPt", jets->GetArrayName().Data());
1538  Double_t a[4] = {fCent, fDijet.trigJetPt, 1.*trigJet->GetNumberOfConstituents(), type};
1539  fHistManager.FillTHnSparse(histname, a);
1540 
1541 }
1542 
1547 {
1548  Double_t contents[30]={0};
1549  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1550  if (!histJetObservables) return;
1551  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1552  TString title(histJetObservables->GetAxis(n)->GetTitle());
1553  if (title=="A_{J}")
1554  contents[n] = fDijet.AJ;
1555  else if (title=="#Delta#phi")
1556  contents[n] = deltaPhi;
1557  else if (title=="#it{p}_{T,particle} (GeV/#it{c})")
1558  contents[n] = trackPt;
1559  else if (title=="#it{p}_{T}#parallel (GeV/#it{c})")
1560  contents[n] = balancePt;
1561  else
1562  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1563  }
1564  histJetObservables->Fill(contents);
1565 }
1566 
1571 {
1572  // Matching efficiency histogram
1573  TString histname = "GeometricalMatchingEfficiency";
1574 
1575  // If we have a matching di-jet, fill the geometrical matching histogram
1576  if (fMatchingDijet.assJet) {
1577  fHistManager.FillTH1(histname.Data(), 1);
1578 
1581  Bool_t isSwitched = trigDeltaR > fMatchingJetR;
1582 
1583  TString thnname = "GeometricalMatching";
1584  Double_t contents[30]={0};
1585  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(thnname.Data()));
1586  if (!histJetObservables) return;
1587  for (Int_t n = 0; n < histJetObservables->GetNdimensions(); n++) {
1588  TString title(histJetObservables->GetAxis(n)->GetTitle());
1589  if (title=="Centrality (%)")
1590  contents[n] = fCent;
1591  else if (title=="isSwitched")
1592  contents[n] = isSwitched;
1593  else if (title=="#DeltaR_{trig}")
1594  contents[n] = trigDeltaR;
1595  else if (title=="#DeltaR_{ass}")
1596  contents[n] = assDeltaR;
1597  else if (title=="trig #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1598  contents[n] = fMatchingDijet.trigJetPt - fDijet.trigJetPt;
1599  else if (title=="ass #it{p}_{T,low-thresh} - #it{p}_{T,hard-core}")
1600  contents[n] = fMatchingDijet.assJetPt - fDijet.assJetPt;
1601  else if (title=="A_{J} low-threshold")
1602  contents[n] = fMatchingDijet.AJ;
1603  else if (title=="A_{J} hard-core")
1604  contents[n] = fDijet.AJ;
1605  else
1606  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1607  }
1608  histJetObservables->Fill(contents);
1609  }
1610  else {
1611  fHistManager.FillTH1(histname.Data(), 0.);
1612  }
1613 }
AliEmcalJet * GetLeadingJet(AliJetContainer *jetCont)
Double_t Area() const
Definition: AliEmcalJet.h:130
Double_t GetRhoVal() const
Int_t fNPhiBins
Number of phi bins in DCal region (for background/correction)
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
PHOS acceptance.
Definition: AliEmcalJet.h:75
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:121
static AliEmcalDownscaleFactorsOCDB * Instance()
const AliClusterIterableMomentumContainer accepted_momentum() const
Int_t fNEtaBins
Number of eta bins in DCal region (for background/correction)
Double_t Phi() const
Definition: AliEmcalJet.h:117
Container with name, TClonesArray and cuts for particles.
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
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:140
void SetCaloTriggerPatchInfoName(const char *n)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
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.
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
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:367
Double_t Phi_0_2pi() const
EMCal acceptance.
Definition: AliEmcalJet.h:69
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.
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:129
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)
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
Calculates the fraction of momentum z of part 1 w.r.t. part 2 in the direction of part 2...
void FillMomentumBalanceHistograms(TString histname, Double_t deltaPhi, Double_t trackPt, Double_t balancePt)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
TObjArray fJetCollArray
jet collection array
AliRhoParameter * GetRhoParameter()
Bool_t fComputeMBDownscaling
Set whether to compute and plot MB downscaling factors.
Double_t Pt() const
Definition: AliEmcalJet.h:109
virtual Bool_t IsEventSelected()
Performing event selection.
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
AliEmcalList * fOutput
!output list
Bool_t fLoadBackgroundScalingWeights
Flag to load eta-phi weights for full-jet background scale factors.
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:51
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 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
Double_t fTrackConstituentThreshold
constituent threshold for matching study
DCal acceptance – spans ONLY DCal (no PHOS or gap)
Definition: AliEmcalJet.h:73
Double_t NEF() const
Definition: AliEmcalJet.h:148
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
Double_t fMinTrigJetPt
leading jet min pT in a dijet pair
Container structure for EMCAL clusters.
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Container for jet within the EMCAL jet framework.
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.
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()