AliPhysics  5eaf189 (5eaf189)
AliAnalysisTaskPWGJEQA.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 <cstring>
17 
18 #include <TClonesArray.h>
19 #include <TH1F.h>
20 #include <TH2F.h>
21 #include <THnSparse.h>
22 #include <TMath.h>
23 #include <TString.h>
24 #include <Riostream.h>
25 #include <TChain.h>
26 #include <TFile.h>
27 
28 #include <AliVEventHandler.h>
29 #include <AliAnalysisManager.h>
30 #include "AliParticleContainer.h"
31 #include "AliClusterContainer.h"
32 #include "AliTrackContainer.h"
33 #include "AliCentrality.h"
34 #include "AliVCluster.h"
35 #include "AliVParticle.h"
36 #include "AliVTrack.h"
37 #include "AliLog.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliEMCALGeoParams.h"
40 #include "AliESDtrack.h"
41 #include "AliAODTrack.h"
42 #include "AliAODMCParticle.h"
43 #include "AliMCParticleContainer.h"
44 #include "AliEmcalJet.h"
45 #include "AliJetContainer.h"
46 #include "AliMCEvent.h"
47 #include "AliGenPythiaEventHeader.h"
48 #include "AliAODMCHeader.h"
49 #include "AliPHOSGeometry.h"
50 #include "AliOADBContainer.h"
51 
52 #include "AliAnalysisTaskPWGJEQA.h"
53 
54 using std::cout;
55 using std::endl;
56 
58 ClassImp(AliAnalysisTaskPWGJEQA);
60 
61 //________________________________________________________________________
64  AliAnalysisTaskEmcalJet("AliAnalysisTaskPWGJEQA", kTRUE),
65  fCellEnergyCut(0.05),
66  fMaxPt(250),
67  fNTotTracks(0),
68  fLeadingTrack(),
69  fDoTrackQA(kTRUE),
70  fDoCaloQA(kTRUE),
71  fDoJetQA(kTRUE),
72  fDoEventQA(kTRUE),
73  fGeneratorLevelName(),
74  fDetectorLevelName(),
75  fRejectOutlierEvents(kFALSE),
76  fIsPtHard(kFALSE),
77  fGeneratorLevel(0),
78  fDetectorLevel(0),
79  fNPtHistBins(0),
80  fPtHistBins(0),
81  fNEtaHistBins(0),
82  fEtaHistBins(0),
83  fNPhiHistBins(0),
84  fPhiHistBins(0),
85  fNCentHistBins(0),
86  fCentHistBins(0),
87  fNPtRelDiffHistBins(0),
88  fPtRelDiffHistBins(0),
89  fNPtResHistBins(0),
90  fPtResHistBins(0),
91  fNIntegerHistBins(0),
92  fIntegerHistBins(0),
93  fPHOSGeo(nullptr),
94  fHistManager("AliAnalysisTaskPWGJEQA")
95 {
96  // Default constructor.
97  memset(fNTotClusters, 0, sizeof(Int_t)*3);
100 }
101 
102 //________________________________________________________________________
107  AliAnalysisTaskEmcalJet(name, kTRUE),
108  fCellEnergyCut(0.05),
109  fMaxPt(250),
110  fNTotTracks(0),
111  fLeadingTrack(),
112  fDoTrackQA(kTRUE),
113  fDoCaloQA(kTRUE),
114  fDoJetQA(kTRUE),
115  fDoEventQA(kTRUE),
118  fRejectOutlierEvents(kFALSE),
119  fIsPtHard(kFALSE),
120  fGeneratorLevel(0),
121  fDetectorLevel(0),
122  fNPtHistBins(0),
123  fPtHistBins(0),
124  fNEtaHistBins(0),
125  fEtaHistBins(0),
126  fNPhiHistBins(0),
127  fPhiHistBins(0),
128  fNCentHistBins(0),
129  fCentHistBins(0),
132  fNPtResHistBins(0),
133  fPtResHistBins(0),
135  fIntegerHistBins(0),
136  fPHOSGeo(nullptr),
137  fHistManager(name)
138 {
139  // Standard
140  memset(fNTotClusters, 0, sizeof(Int_t)*3);
143 }
144 
145 //________________________________________________________________________
147 {
148  // Destructor
149 }
150 
151 //________________________________________________________________________
153 {
154  // Create histograms
156 
157  // Set track container pointers
160 
161  // Allocate histograms for tracks, cells, clusters, jets
167 
168  TIter next(fHistManager.GetListOfHistograms());
169  TObject* obj = 0;
170  while ((obj = next())) {
171  fOutput->Add(obj);
172  }
173 
174  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
175 }
176 
177 //________________________________________________________________________
179 
181  if (fGeneratorLevel) {
184  }
185 }
186 
187 //________________________________________________________________________
189 
190  TString histname;
191  TString title;
192 
193  if (!fCaloCellsName.IsNull()) {
194 
195  histname = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
196  title = histname + ";#it{E}_{cell} (GeV);counts";
197  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,0, 100);
198 
199  histname = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
200  title = histname + ";cell absId;<#it{E}_{cell}> (GeV)";
201  fHistManager.CreateTProfile(histname.Data(), title.Data(), 18000,0,18000);
202 
203  histname = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
204  title = histname + ";#it{t}_{cell} (s);counts";
205  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,-5e-6, 5e-6);
206 
207  histname = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
208  title = histname + ";cell absId;<#it{t}_{cell}> (s)";
209  fHistManager.CreateTProfile(histname.Data(), title.Data(), 18000,0,18000);
210 
211  histname = TString::Format("%s/fHistCellEvsTime", fCaloCellsName.Data());
212  title = histname + ";<#it{t}_{cell}> (s);<#it{E}_{cell}> (GeV)";
213  fHistManager.CreateTH2(histname.Data(), title.Data(), 1000, -1e-6, 1e-6, 200, 0, 100);
214 
215  }
216 }
217 
218 //________________________________________________________________________
220 
221  AliEmcalContainer* cont = 0;
222  TString histname;
223  TString title;
224  Int_t nPtBins = 2*fMaxPt;
225 
226  TIter nextClusColl(&fClusterCollArray);
227  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
228 
229  histname = TString::Format("%s/fHistClusterRejectionReason", cont->GetArrayName().Data());
230  title = histname + ";Rejection reason;#it{E}_{clus} (GeV/);counts";
231  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
232  SetRejectionReasonLabels(hist->GetXaxis());
233 
234  histname = TString::Format("%s/hPhosNCellsVsEnergy", cont->GetArrayName().Data());
235  title = histname + ";N cells;#it{E}_{clus} (GeV/);counts";
236  fHistManager.CreateTH2(histname.Data(), title.Data(), 40, 0, 40, nPtBins, 0, fMaxPt);
237 
238  histname = TString::Format("%s/hPhosM02VsEnergy", cont->GetArrayName().Data());
239  title = histname + ";M02;#it{E}_{clus} (GeV/);counts";
240  fHistManager.CreateTH2(histname.Data(), title.Data(), 100, 0, 20, nPtBins, 0, fMaxPt);
241 
242  histname = TString::Format("%s/hPhosCellIdVsEnergy", cont->GetArrayName().Data());
243  title = histname + ";Cell ID;#it{E}_{clus} (GeV/);counts";
244  fHistManager.CreateTH2(histname.Data(), title.Data(), 20000, 0, 20000, nPtBins, 0, fMaxPt);
245 
246  const Int_t nEmcalSM = 20;
247  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
248  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
249  title = histname + ";#it{E}_{cluster} (GeV);counts";
250  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
251  }
252 
253  for (Int_t sm = 1; sm < 5; sm++) {
254  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
255  title = histname + ";#it{E}_{cluster} (GeV);counts";
256  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
257  }
258 
259  // Cluster THnSparse
260  Int_t dim = 0;
261  TString title[20];
262  Int_t nbins[20] = {0};
263  Double_t min[30] = {0.};
264  Double_t max[30] = {0.};
265  Double_t *binEdges[20] = {0};
266 
268  title[dim] = "Centrality %";
269  nbins[dim] = fNCentHistBins;
270  binEdges[dim] = fCentHistBins;
271  min[dim] = fCentHistBins[0];
272  max[dim] = fCentHistBins[fNCentHistBins];
273  dim++;
274  }
275 
276  title[dim] = "#it{E}_{clus} (GeV)";
277  nbins[dim] = fNPtHistBins;
278  binEdges[dim] = fPtHistBins;
279  min[dim] = fPtHistBins[0];
280  max[dim] = fPtHistBins[fNPtHistBins];
281  dim++;
282 
283  title[dim] = "#eta";
284  nbins[dim] = fNEtaHistBins;
285  binEdges[dim] = fEtaHistBins;
286  min[dim] = fEtaHistBins[0];
287  max[dim] = fEtaHistBins[fNEtaHistBins];
288  dim++;
289 
290  title[dim] = "#phi";
291  nbins[dim] = fNPhiHistBins;
292  binEdges[dim] = fPhiHistBins;
293  min[dim] = fPhiHistBins[0];
294  max[dim] = fPhiHistBins[fNPhiHistBins];
295  dim++;
296 
297  title[dim] = "cluster type";
298  nbins[dim] = 3;
299  binEdges[dim] = fIntegerHistBins;
300  min[dim] = fIntegerHistBins[0];
301  max[dim] = fIntegerHistBins[3];
302  dim++;
303 
304  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
305  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
306  for (Int_t i = 0; i < dim; i++) {
307  hn->GetAxis(i)->SetTitle(title[i]);
308  hn->SetBinEdges(i, binEdges[i]);
309  }
310  }
311 }
312 
313 //________________________________________________________________________
315 
316  AliJetContainer* jets = 0;
317  TIter nextJetColl(&fJetCollArray);
318  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
319 
320  // Allocate THnSparse
321 
322  TString axisTitle[30]= {""};
323  Int_t nbins[30] = {0};
324  Double_t min[30] = {0.};
325  Double_t max[30] = {0.};
326  Double_t *binEdges[20] = {0};
327  Int_t dim = 0;
328 
329  if (fForceBeamType != kpp) {
330  axisTitle[dim] = "Centrality (%)";
331  nbins[dim] = fNCentHistBins;
332  binEdges[dim] = fCentHistBins;
333  min[dim] = fCentHistBins[0];
334  max[dim] = fCentHistBins[fNCentHistBins];
335  dim++;
336  }
337 
338  axisTitle[dim] = "#eta_{jet}";
339  nbins[dim] = fNEtaHistBins;
340  binEdges[dim] = fEtaHistBins;
341  min[dim] = fEtaHistBins[0];
342  max[dim] = fEtaHistBins[fNEtaHistBins];
343  dim++;
344 
345  axisTitle[dim] = "#phi_{jet} (rad)";
346  nbins[dim] = fNPhiHistBins;
347  binEdges[dim] = fPhiHistBins;
348  min[dim] = fPhiHistBins[0];
349  max[dim] = fPhiHistBins[fNPhiHistBins];
350  dim++;
351 
352  axisTitle[dim] = "#it{p}_{T} (GeV/#it{c})";
353  nbins[dim] = TMath::CeilNint(fMaxPt/3);
354  min[dim] = 0;
355  max[dim] = fMaxPt;
356  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
357  dim++;
358 
359  if (fForceBeamType != kpp) {
360  axisTitle[dim] = "#it{p}_{T}^{corr} (GeV/#it{c})";
361  nbins[dim] = TMath::CeilNint(fMaxPt/3);
362  min[dim] = -fMaxPt/2 + 25;
363  max[dim] = fMaxPt/2 + 25;
364  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
365  dim++;
366  }
367 
368  axisTitle[dim] = "#it{p}_{T,particle}^{leading} (GeV/#it{c})";
369  nbins[dim] = fMaxPt/5;
370  min[dim] = 0;
371  max[dim] = 150;
372  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
373  dim++;
374 
375  TString thnname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
376  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
377  for (Int_t i = 0; i < dim; i++) {
378  hn->GetAxis(i)->SetTitle(axisTitle[i]);
379  hn->SetBinEdges(i, binEdges[i]);
380  }
381 
382  // Allocate other jet histograms
383  TString histname;
384  TString title;
385  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
386  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
387  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
388  SetRejectionReasonLabels(hist->GetXaxis());
389 
390  if (!jets->GetRhoName().IsNull()) {
391  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
392  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
393  fHistManager.CreateTH2(histname.Data(), title.Data(), 101, 0, 101, 100, 0, 500);
394  }
395  if (fForceBeamType != kpp) {
396  histname = TString::Format("%s/hNEFVsPtEMC", jets->GetArrayName().Data());
397  title = histname + ";Centrality (%); #it{p}_{T}^{corr} (GeV/#it{c});NEF";
398  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, 250, 0, 250, 50, 0, 1);
399 
400  histname = TString::Format("%s/hNEFVsPtDCal", jets->GetArrayName().Data());
401  title = histname + ";Centrality (%); #it{p}_{T}^{corr} (GeV/#it{c});NEF";
402  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, 250, 0, 250, 50, 0, 1);
403 
404  histname = TString::Format("%s/hDeltaEHadCorr", jets->GetArrayName().Data());
405  title = histname + ";Centrality (%); #it{p}_{T}^{corr} (GeV/#it{c}); #sum #it{E_{nonlincorr}} - #it{E}_{hadcorr}";
406  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, 250, 0, 250, 250, 0, 50);
407  }
408  else{
409  histname = TString::Format("%s/hNEFVsPtEMC", jets->GetArrayName().Data());
410  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
411  fHistManager.CreateTH2(histname.Data(), title.Data(), 250, 0, 250, 50, 0, 1);
412 
413  histname = TString::Format("%s/hNEFVsPtDCal", jets->GetArrayName().Data());
414  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
415  fHistManager.CreateTH2(histname.Data(), title.Data(), 250, 0, 250, 50, 0, 1);
416 
417  histname = TString::Format("%s/hDeltaEHadCorr", jets->GetArrayName().Data());
418  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c}); #sum#it{E_{nonlincorr}} - #it{E}_{hadcorr}";
419  fHistManager.CreateTH2(histname.Data(), title.Data(), 250, 0, 250, 250, 0, 50);
420  }
421  }
422 }
423 
424 //________________________________________________________________________
426 
427  Int_t nPtBins = 2*fMaxPt;
428 
429  Int_t dim = 0;
430  TString axistitle[40];
431  Int_t nbins[40] = {0};
432  Double_t min[40] = {0};
433  Double_t max[40] = {0};
434  Double_t *binEdges[20] = {0};
435 
437  axistitle[dim] = "Centrality %";
438  nbins[dim] = fNCentHistBins;
439  binEdges[dim] = fCentHistBins;
440  min[dim] = fCentHistBins[0];
441  max[dim] = fCentHistBins[fNCentHistBins];
442  dim++;
443  }
444 
445  if (fParticleCollArray.GetEntriesFast()>0) {
446  axistitle[dim] = "No. of tracks";
448  nbins[dim] = 50;
449  min[dim] = -0.5;
450  max[dim] = 10000-0.5;
451  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
452  }
453  else {
454  nbins[dim] = 50;
455  min[dim] = 0;
456  max[dim] = 200;
457  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
458  }
459  dim++;
460 
461  axistitle[dim] = "#it{p}_{T,track}^{leading} (GeV/c)";
462  nbins[dim] = fMaxPt/5;
463  min[dim] = 0;
464  max[dim] = 150;
465  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
466  dim++;
467  }
468 
469  if (fClusterCollArray.GetEntriesFast()>0) {
470  axistitle[dim] = "No. of clusters";
471 
473  nbins[dim] = 50;
474  min[dim] = -0.5;
475  max[dim] = 4000-0.5;
476  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
477  }
478  else {
479  nbins[dim] = 50;
480  min[dim] = 0;
481  max[dim] = 200;
482  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
483  }
484  dim++;
485 
486  axistitle[dim] = "#it{E}_{cluster}^{leading} (GeV)";
487  nbins[dim] = fMaxPt/5;
488  min[dim] = 0;
489  max[dim] = 150;
490  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
491  dim++;
492  }
493 
494  THnSparse* hn = fHistManager.CreateTHnSparse("eventQA","eventQA",dim,nbins,min,max);
495  for (Int_t i = 0; i < dim; i++) {
496  hn->GetAxis(i)->SetTitle(axistitle[i]);
497  hn->SetBinEdges(i, binEdges[i]);
498  }
499 
500  if (fIsPtHard) {
501  TString histname = "hPtHard";
502  TString title = histname + ";#it{p}_{T,hard} (GeV/c);counts";
503  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
504 
505  TH1* hTrials = fHistManager.CreateTH1("hNtrials", "hNtrials", 1, 0, 1);
506  hTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
507 
508  TH1* hxsection = fHistManager.CreateTH1("hXsec", "hXsec", 1, 0, 1);
509  hxsection->GetXaxis()->SetBinLabel(1,"<#sigma>");
510  }
511 
512  fHistEventRejection->GetXaxis()->SetBinLabel(15,"PtHardBinJetOutlier");
513 
514 }
515 
516 //________________________________________________________________________
518 {
519  Int_t dim = 0;
520  TString title[20];
521  Int_t nbins[20] = {0};
522  Double_t min[30] = {0.};
523  Double_t max[30] = {0.};
524  Double_t *binEdges[20] = {0};
525 
527  title[dim] = "Centrality %";
528  nbins[dim] = fNCentHistBins;
529  binEdges[dim] = fCentHistBins;
530  min[dim] = fCentHistBins[0];
531  max[dim] = fCentHistBins[fNCentHistBins];
532  dim++;
533  }
534 
535  title[dim] = "#it{p}_{T} (GeV/#it{c})";
536  nbins[dim] = fNPtHistBins;
537  binEdges[dim] = fPtHistBins;
538  min[dim] = fPtHistBins[0];
539  max[dim] = fPtHistBins[fNPtHistBins];
540  dim++;
541 
542  title[dim] = "#eta";
543  nbins[dim] = fNEtaHistBins;
544  binEdges[dim] = fEtaHistBins;
545  min[dim] = fEtaHistBins[0];
546  max[dim] = fEtaHistBins[fNEtaHistBins];
547  dim++;
548 
549  title[dim] = "#phi";
550  nbins[dim] = fNPhiHistBins;
551  binEdges[dim] = fPhiHistBins;
552  min[dim] = fPhiHistBins[0];
553  max[dim] = fPhiHistBins[fNPhiHistBins];
554  dim++;
555 
556  title[dim] = "track type";
557  nbins[dim] = 3;
558  binEdges[dim] = fIntegerHistBins;
559  min[dim] = fIntegerHistBins[0];
560  max[dim] = fIntegerHistBins[3];
561  dim++;
562 
563  title[dim] = "#sigma(#it{p}_{T}) / #it{p}_{T}";
564  nbins[dim] = fNPtResHistBins;
565  binEdges[dim] = fPtResHistBins;
566  min[dim] = fPtResHistBins[0];
567  max[dim] = fPtResHistBins[fNPtResHistBins];
568  dim++;
569 
570  THnSparse* hn = fHistManager.CreateTHnSparse("tracks","tracks", dim, nbins, min, max);
571  for (Int_t i = 0; i < dim; i++) {
572  hn->GetAxis(i)->SetTitle(title[i]);
573  hn->SetBinEdges(i, binEdges[i]);
574  }
575 
576 }
577 
578 //________________________________________________________________________
580 {
581  Int_t dim = 0;
582  TString title[20];
583  Int_t nbins[20] = {0};
584  Double_t min[30] = {0.};
585  Double_t max[30] = {0.};
586  Double_t *binEdges[20] = {0};
587 
589  title[dim] = "Centrality %";
590  nbins[dim] = fNCentHistBins;
591  binEdges[dim] = fCentHistBins;
592  min[dim] = fCentHistBins[0];
593  max[dim] = fCentHistBins[fNCentHistBins];
594  dim++;
595  }
596 
597  title[dim] = "#it{p}_{T} (GeV/#it{c})";
598  nbins[dim] = fNPtHistBins;
599  binEdges[dim] = fPtHistBins;
600  min[dim] = fPtHistBins[0];
601  max[dim] = fPtHistBins[fNPtHistBins];
602  dim++;
603 
604  title[dim] = "#eta";
605  nbins[dim] = fNEtaHistBins;
606  binEdges[dim] = fEtaHistBins;
607  min[dim] = fEtaHistBins[0];
608  max[dim] = fEtaHistBins[fNEtaHistBins];
609  dim++;
610 
611  title[dim] = "#phi";
612  nbins[dim] = fNPhiHistBins;
613  binEdges[dim] = fPhiHistBins;
614  min[dim] = fPhiHistBins[0];
615  max[dim] = fPhiHistBins[fNPhiHistBins];
616  dim++;
617 
618  title[dim] = "Findable";
619  nbins[dim] = 2;
620  binEdges[dim] = fIntegerHistBins;
621  min[dim] = fIntegerHistBins[0];
622  max[dim] = fIntegerHistBins[2];
623  dim++;
624 
625  THnSparse* hn = fHistManager.CreateTHnSparse("tracks_PhysPrim","tracks_PhysPrim", dim, nbins, min, max);
626  for (Int_t i = 0; i < dim; i++) {
627  hn->GetAxis(i)->SetTitle(title[i]);
628  hn->SetBinEdges(i, binEdges[i]);
629  }
630 }
631 
632 //________________________________________________________________________
634 {
635  Int_t dim = 0;
636  TString title[20];
637  Int_t nbins[20] = {0};
638  Double_t min[30] = {0.};
639  Double_t max[30] = {0.};
640  Double_t *binEdges[20] = {0};
641 
643  title[dim] = "Centrality %";
644  nbins[dim] = fNCentHistBins;
645  binEdges[dim] = fCentHistBins;
646  min[dim] = fCentHistBins[0];
647  max[dim] = fCentHistBins[fNCentHistBins];
648  dim++;
649  }
650 
651  title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
652  nbins[dim] = fNPtHistBins;
653  binEdges[dim] = fPtHistBins;
654  min[dim] = fPtHistBins[0];
655  max[dim] = fPtHistBins[fNPtHistBins];
656  dim++;
657 
658  title[dim] = "#eta^{gen}";
659  nbins[dim] = fNEtaHistBins;
660  binEdges[dim] = fEtaHistBins;
661  min[dim] = fEtaHistBins[0];
662  max[dim] = fEtaHistBins[fNEtaHistBins];
663  dim++;
664 
665  title[dim] = "#phi^{gen}";
666  nbins[dim] = fNPhiHistBins;
667  binEdges[dim] = fPhiHistBins;
668  min[dim] = fPhiHistBins[0];
669  max[dim] = fPhiHistBins[fNPhiHistBins];
670  dim++;
671 
672  title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
673  nbins[dim] = fNPtHistBins;
674  binEdges[dim] = fPtHistBins;
675  min[dim] = fPtHistBins[0];
676  max[dim] = fPtHistBins[fNPtHistBins];
677  dim++;
678 
679  title[dim] = "#eta^{det}";
680  nbins[dim] = fNEtaHistBins;
681  binEdges[dim] = fEtaHistBins;
682  min[dim] = fEtaHistBins[0];
683  max[dim] = fEtaHistBins[fNEtaHistBins];
684  dim++;
685 
686  title[dim] = "#phi^{det}";
687  nbins[dim] = fNPhiHistBins;
688  binEdges[dim] = fPhiHistBins;
689  min[dim] = fPhiHistBins[0];
690  max[dim] = fPhiHistBins[fNPhiHistBins];
691  dim++;
692 
693  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
694  nbins[dim] = fNPtRelDiffHistBins;
695  binEdges[dim] = fPtRelDiffHistBins;
696  min[dim] = fPtRelDiffHistBins[0];
698  dim++;
699 
700  title[dim] = "track type";
701  nbins[dim] = 3;
702  binEdges[dim] = fIntegerHistBins;
703  min[dim] = fIntegerHistBins[0];
704  max[dim] = fIntegerHistBins[3];
705  dim++;
706 
707  THnSparse* hn = fHistManager.CreateTHnSparse("tracks_Matched","tracks_Matched", dim, nbins, min, max);
708  for (Int_t i = 0; i < dim; i++) {
709  hn->GetAxis(i)->SetTitle(title[i]);
710  hn->SetBinEdges(i, binEdges[i]);
711  }
712 }
713 
714 
715 //________________________________________________________________________
718 {
719  fNPtHistBins = 82;
722  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
723  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
724  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
725  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
726  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
727  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
728 
729  fNEtaHistBins = 70;
732 
733  fNPhiHistBins = 101;
735  GenerateFixedBinArray(fNPhiHistBins, 0, TMath::Pi() * 2.02, fPhiHistBins);
736 
737  fNCentHistBins = 4;
739  fCentHistBins[0] = 0;
740  fCentHistBins[1] = 10;
741  fCentHistBins[2] = 30;
742  fCentHistBins[3] = 50;
743  fCentHistBins[4] = 90;
744 
745  fNPtResHistBins = 100;
748  GenerateFixedBinArray(25, 0.05, 0.10, fPtResHistBins+50);
749  GenerateFixedBinArray(15, 0.10, 0.20, fPtResHistBins+75);
750  GenerateFixedBinArray(10, 0.20, 0.50, fPtResHistBins+90);
751 
752  fNPtRelDiffHistBins = 50;
755 
756  fNIntegerHistBins = 3;
759 }
760 
761 //________________________________________________________________________
763 {
764  if (fClusterCollArray.GetEntriesFast() == 0 && fCaloCellsName.IsNull()) {
765  fNeedEmcalGeom = kFALSE;
766  }
767  else {
768  fNeedEmcalGeom = kTRUE;
769  }
770 
772 
773  // Load the PHOS geometry
774  fPHOSGeo = AliPHOSGeometry::GetInstance();
775  if (fPHOSGeo) {
776  AliInfo("Found instance of PHOS geometry!");
777  }
778  else {
779  AliInfo("Creating PHOS geometry!");
780  Int_t runNum = InputEvent()->GetRunNumber();
781  if(runNum<209122) //Run1
782  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
783  else
784  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
785 
786  if (fPHOSGeo) {
787  AliOADBContainer geomContainer("phosGeo");
788  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
789  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
790  for(Int_t mod=0; mod<6; mod++) {
791  if(!matrixes->At(mod)) continue;
792  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
793  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
794  ((TGeoHMatrix*)matrixes->At(mod))->Print();
795  }
796  }
797  }
798 }
799 
800 //________________________________________________________________________
802 {
803  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
804 
805  // If Pt-hard production, get the Pt-hard of the event, and have possibility to reject the event for jet outliers
806  if (fIsPtHard) {
807  AliGenPythiaEventHeader* pygen = nullptr;
808  if (MCEvent()) {
809  pygen = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
810  if (!pygen) {
811  // Check if AOD
812  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
813  if (aodMCH) {
814  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
815  pygen = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
816  if (pygen) break;
817  }
818  }
819  }
820  }
821  if (pygen) {
822  fPtHard = pygen->GetPtHard();
823  //fNTrials = pygen->Trials();
824  //fXsection = pygen->GetXsection();
825 
826  // reject outlier events, where the jet Pt is much larger than Pt-hard
827  if (fRejectOutlierEvents) {
828  AliTLorentzVector jet;
829  Int_t nTriggerJets = pygen->NTriggerJets();
830  Float_t tmpjet[]={0,0,0,0};
831  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
832  pygen->TriggerJet(ijet, tmpjet);
833  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
834 
835  if (jet.Pt() > 4. * fPtHard) {
836  //AliInfo(Form("Reject jet event with: pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), 4.));
837  if (fGeneralHistograms) fHistEventRejection->Fill("PtHardBinJetOutlier",1);
838  return kFALSE;
839  }
840  }
841  }
842  }
843  }
844 
845  return kTRUE;
846 }
847 
855 //________________________________________________________________________
857 {
858  if (fIsPtHard) {
859  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
860  if (!tree) {
861  AliError(Form("%s - UserNotify: No current tree!",GetName()));
862  return kFALSE;
863  }
864 
865  Float_t xsection = 0;
866  Float_t trials = 0;
867  Int_t pthardbin = 0;
868 
869  TFile *curfile = tree->GetCurrentFile();
870  if (!curfile) {
871  AliError(Form("%s - UserNotify: No current file!",GetName()));
872  return kFALSE;
873  }
874 
875  TChain *chain = dynamic_cast<TChain*>(tree);
876  if (chain) tree = chain->GetTree();
877 
878  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
879 
880  fHistManager.FillTH1("hNtrials", "#sum{ntrials}", trials);
881  fHistManager.FillTH1("hXsec", "<#sigma>", xsection);
882  }
883  return kTRUE;
884 }
885 
886 //________________________________________________________________________
888 {
894 
895  return kTRUE;
896 }
897 
898 //________________________________________________________________________
900 
901  fNTotTracks = 0;
902  fLeadingTrack.SetPxPyPzE(0,0,0,0);
903 
904  AliVTrack* track;
905  for (auto trackIterator : fDetectorLevel->accepted_momentum() ) {
906 
907  fNTotTracks++;
908  if (fLeadingTrack.Pt() < trackIterator.first.Pt()) fLeadingTrack = trackIterator.first;
909 
910  track = trackIterator.second;
911  Byte_t type = fDetectorLevel->GetTrackType(track);
912  if (type <= 2) {
913  Double_t sigma = 0;
914 
915  if (fIsEsd) {
916 
917  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
918  if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
919 
920  } else { // AOD
921 
922  AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
923  if(!aodtrack) AliFatal("Not a standard AOD");
924 
925  AliExternalTrackParam exParam;
926 
927  //get covariance matrix
928  Double_t cov[21] = {0,};
929  aodtrack->GetCovMatrix(cov);
930  Double_t pxpypz[3] = {0,};
931  aodtrack->PxPyPz(pxpypz);
932  Double_t xyz[3] = {0,};
933  aodtrack->GetXYZ(xyz);
934  Short_t sign = aodtrack->Charge();
935  exParam.Set(xyz,pxpypz,cov,sign);
936 
937  sigma = TMath::Sqrt(exParam.GetSigma1Pt2());
938 
939  }
940 
941  Int_t label = TMath::Abs(track->GetLabel());
942 
943  FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, type);
944 
945  if (fGeneratorLevel && label > 0) {
946  AliAODMCParticle *part = fGeneratorLevel->GetAcceptMCParticleWithLabel(label);
947  if (part) {
948  if (part->GetGeneratorIndex() == 0) {
949  Int_t pdg = TMath::Abs(part->PdgCode());
950  // select charged pions, protons, kaons, electrons, muons
951  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
952  FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
953  }
954  }
955  }
956  }
957  }
958  else {
959  AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
960  }
961  }
962 
963  if (fGeneratorLevel) {
964  AliAODMCParticle* part;
965  for (auto partIterator : fGeneratorLevel->accepted_momentum() ) {
966  part = partIterator.second;
967 
968  Byte_t findable = 0;
969 
970  Int_t pdg = TMath::Abs(part->PdgCode());
971  // select charged pions, protons, kaons, electrons, muons
972  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
973 
974  FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), findable);
975  }
976  }
977 }
978 
979 //________________________________________________________________________
981 
982  if (!fCaloCells) return;
983 
984  TString histname_energy = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
985  TString histname_time = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
986  TString histname_energyProf = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
987  TString histname_timeProf = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
988  TString histname_EnergyvsTime = TString::Format("%s/fHistCellEvsTime", fCaloCellsName.Data());
989 
990  const Int_t ncells = fCaloCells->GetNumberOfCells();
991  for (Int_t pos = 0; pos < ncells; pos++) {
992  Float_t amp = fCaloCells->GetAmplitude(pos);
993  Float_t time = fCaloCells->GetTime(pos);
994  Int_t absId = fCaloCells->GetCellNumber(pos);
995 
996  if (amp < fCellEnergyCut) continue;
997 
998  fHistManager.FillTH1(histname_energy, amp);
999  fHistManager.FillTH1(histname_time, time);
1000 
1001  fHistManager.FillProfile(histname_energyProf, absId, amp);
1002  fHistManager.FillProfile(histname_timeProf, absId, time);
1003 
1004  fHistManager.FillTH2(histname_EnergyvsTime, time, amp);
1005  }
1006 }
1007 
1008 //________________________________________________________________________
1010 
1011  memset(fNTotClusters, 0, sizeof(Int_t)*3);
1012  for (Int_t i = 0; i < 3; i++) fLeadingCluster[i].SetPxPyPzE(0,0,0,0);
1013 
1014  TString histname;
1015  AliClusterContainer* clusters = 0;
1016  TIter nextClusColl(&fClusterCollArray);
1017  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
1018  // Cluster loop
1020  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1021 
1022  UInt_t rejectionReason = 0;
1023  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
1024  histname = TString::Format("%s/fHistClusterRejectionReason", clusters->GetArrayName().Data());
1025  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
1026  continue;
1027  }
1028 
1029  // Determine cluster type (EMCal/DCal/Phos) and record relevant eventQA info
1030  ClusterType clusType = kNA;
1031  if (it->second->IsEMCAL()) {
1032  Double_t phi = it->first.Phi_0_2pi();
1033  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
1034  if (isDcal == 0) {
1035  clusType = kEMCal;
1036  } else if (isDcal == 1) {
1037  clusType = kDCal;
1038  }
1039 
1040  if (fLeadingCluster[isDcal].E() < it->first.E()) fLeadingCluster[isDcal] = it->first;
1041  fNTotClusters[isDcal]++;
1042 
1043  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
1044  if (sm >=0 && sm < 20) {
1045  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1046  fHistManager.FillTH1(histname, it->second->E());
1047  }
1048  else {
1049  AliError(Form("Supermodule %d does not exist!", sm));
1050  }
1051 
1052  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
1053 
1054  Int_t nCells = it->second->GetNCells();
1055  Double_t M02 = it->second->GetM02();
1056  Double_t energy = it->second->E();
1057 
1058  histname = TString::Format("%s/hPhosNCellsVsEnergy", clusters->GetArrayName().Data());
1059  fHistManager.FillTH2(histname, nCells, energy);
1060 
1061  histname = TString::Format("%s/hPhosM02VsEnergy", clusters->GetArrayName().Data());
1062  fHistManager.FillTH2(histname, M02, energy);
1063 
1064  clusType = kPHOS;
1065  fNTotClusters[2]++;
1066  if (fLeadingCluster[2].E() < it->first.E()) fLeadingCluster[2] = it->first;
1067 
1068  // Fill Phos spectra by module
1069  Int_t relid[4];
1070  if (fPHOSGeo) {
1071  fPHOSGeo->AbsToRelNumbering(it->second->GetCellAbsId(0), relid);
1072  Int_t sm = relid[0];
1073  if (sm >=1 && sm < 5) {
1074  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1075  fHistManager.FillTH1(histname, energy);
1076  }
1077  else {
1078  AliError(Form("Supermodule %d does not exist!", sm));
1079  }
1080  }
1081 
1082  // Loop through cells in the cluster, and histogram their energy
1083  histname = TString::Format("%s/hPhosCellIdVsEnergy", clusters->GetArrayName().Data());
1084  for (Int_t i=0; i < nCells; i++) {
1085  fHistManager.FillTH2(histname, it->second->GetCellAbsId(i), energy);
1086  }
1087 
1088  }
1089 
1090  Double_t contents[30]={0};
1091  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
1092  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1093  if (!histClusterObservables) return;
1094  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1095  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1096  if (title=="Centrality %")
1097  contents[i] = fCent;
1098  else if (title=="#it{E}_{clus} (GeV)")
1099  contents[i] = it->first.E();
1100  else if (title=="#eta")
1101  contents[i] = it->first.Eta();
1102  else if (title=="#phi")
1103  contents[i] = it->first.Phi_0_2pi();
1104  else if (title=="cluster type")
1105  contents[i] = clusType;
1106  else
1107  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1108  }
1109  histClusterObservables->Fill(contents);
1110 
1111  }
1112  }
1113 }
1114 
1115 //________________________________________________________________________
1117 
1118  TString histname;
1119  AliJetContainer* jets = 0;
1120  TIter nextJetColl(&fJetCollArray);
1121  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1122  Double_t rhoVal = 0;
1123  if (jets->GetRhoParameter()) {
1124  rhoVal = jets->GetRhoVal();
1125  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
1126  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1127  }
1128 
1129  for (auto jet : jets->all()) {
1130 
1131  UInt_t rejectionReason = 0;
1132  if (!jets->AcceptJet(jet, rejectionReason)) {
1133  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
1134  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1135  continue;
1136  }
1137 
1138  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1139  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
1140 
1141  TLorentzVector leadPart;
1142  jets->GetLeadingHadronMomentum(leadPart, jet);
1143 
1144  Double_t contents[30]={0};
1145  histname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
1146  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1147  if (!histJetObservables) return;
1148  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
1149  TString title(histJetObservables->GetAxis(i)->GetTitle());
1150  if (title=="Centrality (%)")
1151  contents[i] = fCent;
1152  else if (title=="#eta_{jet}")
1153  contents[i] = jet->Eta();
1154  else if (title=="#phi_{jet} (rad)")
1155  contents[i] = jet->Phi_0_2pi();
1156  else if (title=="#it{p}_{T} (GeV/#it{c})")
1157  contents[i] = jet->Pt();
1158  else if (title=="#it{p}_{T}^{MC} (GeV/#it{c})")
1159  contents[i] = jet->MCPt();
1160  else if (title=="#it{p}_{T}^{corr} (GeV/#it{c})")
1161  contents[i] = corrPt;
1162  else if (title=="#it{p}_{T,particle}^{leading} (GeV/#it{c})")
1163  contents[i] = ptLeading;
1164  else
1165  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1166  }
1167  histJetObservables->Fill(contents);
1168 
1169  UInt_t jetType = jet->GetJetAcceptanceType();
1170  if(jetType & AliEmcalJet::kEMCALfid)
1171  {
1172  histname = TString::Format("%s/hNEFVsPtEMC", jets->GetArrayName().Data());
1173  if (fForceBeamType != kpp)fHistManager.FillTH3(histname.Data(), fCent, jet->Pt(), jet->NEF());
1174  else fHistManager.FillTH2(histname.Data(), jet->Pt(), jet->NEF());
1175  }
1176  else if(jetType & AliEmcalJet::kDCALonlyfid)
1177  {
1178  histname = TString::Format("%s/hNEFVsPtDCal", jets->GetArrayName().Data());
1179  if (fForceBeamType != kpp)fHistManager.FillTH3(histname.Data(), fCent, jet->Pt(), jet->NEF());
1180  else fHistManager.FillTH2(histname.Data(), jet->Pt(), jet->NEF());
1181  }
1182  Double_t deltaEhadcorr = 0;
1183  const AliVCluster* clus = nullptr;
1184  Int_t nClusters = jet->GetNumberOfClusters();
1185  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1186  clus = jet->Cluster(iClus);
1187  deltaEhadcorr += (clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy());
1188  }
1189 
1190  histname = TString::Format("%s/hDeltaEHadCorr", jets->GetArrayName().Data());
1191  if (fForceBeamType != kpp)fHistManager.FillTH3(histname, fCent, jet->Pt(), deltaEhadcorr);
1192  else fHistManager.FillTH2(histname, jet->Pt(), deltaEhadcorr);
1193  } //jet loop
1194  }
1195 }
1196 
1197 //________________________________________________________________________
1199 
1200  EventQA_t eventQA;
1201  for (Int_t i = 0; i < 3; i++) {
1202  eventQA.fMaxCluster[i] = fLeadingCluster[i];
1203  }
1204  eventQA.fCent = fCent;
1205  eventQA.fNTracks = fNTotTracks;
1206  eventQA.fNClusters[0] = fNTotClusters[0];
1207  eventQA.fNClusters[1] = fNTotClusters[1];
1208  eventQA.fNClusters[2] = fNTotClusters[2];
1209  eventQA.fMaxTrack = fLeadingTrack;
1210 
1211  Int_t globalNclusters = eventQA.fNClusters[0] + eventQA.fNClusters[1] + eventQA.fNClusters[2];
1212  AliTLorentzVector globalMaxCluster;
1213  for (Int_t i = 0; i < 3; i++) {
1214  if (globalMaxCluster.E() < eventQA.fMaxCluster[i].E()) globalMaxCluster = eventQA.fMaxCluster[i];
1215  }
1216 
1217  THnSparse* histEventQA = static_cast<THnSparse*>(fHistManager.FindObject("eventQA"));
1218  Double_t contents[40]={0};
1219  for (Int_t i = 0; i < histEventQA->GetNdimensions(); i++) {
1220  TString title(histEventQA->GetAxis(i)->GetTitle());
1221  if (title=="Centrality %")
1222  contents[i] = eventQA.fCent;
1223  else if (title=="No. of tracks")
1224  contents[i] = eventQA.fNTracks;
1225  else if (title=="No. of clusters")
1226  contents[i] = globalNclusters;
1227  else if (title=="#it{p}_{T,track}^{leading} (GeV/c)")
1228  contents[i] = eventQA.fMaxTrack.Pt();
1229  else if (title=="#it{E}_{cluster}^{leading} (GeV)")
1230  contents[i] = globalMaxCluster.E();
1231  else
1232  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1233  }
1234 
1235  histEventQA->Fill(contents);
1236 
1237  // Fill pythia pt hard histograms, if applicable
1238  if (fPtHard > 1e-6) {
1239  fHistManager.FillTH1("hPtHard", fPtHard);
1240  }
1241 }
1242 
1243 //________________________________________________________________________
1245  Double_t sigma1OverPt, Byte_t trackType)
1246 {
1247  Double_t contents[20]={0};
1248 
1249  THnSparse* thnTracks = static_cast<THnSparse*>(fHistManager.FindObject("tracks"));
1250  if (!thnTracks) return;
1251  for (Int_t i = 0; i < thnTracks->GetNdimensions(); i++) {
1252  TString title(thnTracks->GetAxis(i)->GetTitle());
1253  if (title=="Centrality %")
1254  contents[i] = cent;
1255  else if (title=="#it{p}_{T} (GeV/#it{c})")
1256  contents[i] = trackPt;
1257  else if (title=="#eta")
1258  contents[i] = trackEta;
1259  else if (title=="#phi")
1260  contents[i] = trackPhi;
1261  else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
1262  contents[i] = sigma1OverPt;
1263  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
1264  contents[i] = sigma1OverPt*trackPt;
1265  else if (title=="track type")
1266  contents[i] = trackType;
1267  else
1268  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks->GetName()));
1269  }
1270 
1271  thnTracks->Fill(contents);
1272 }
1273 
1274 //________________________________________________________________________
1275 void AliAnalysisTaskPWGJEQA::FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Byte_t findable)
1276 {
1277  Double_t contents[20]={0};
1278 
1279  THnSparse* thnTracks_PhysPrim = static_cast<THnSparse*>(fHistManager.FindObject("tracks_PhysPrim"));
1280  if (!thnTracks_PhysPrim) return;
1281  for (Int_t i = 0; i < thnTracks_PhysPrim->GetNdimensions(); i++) {
1282  TString title(thnTracks_PhysPrim->GetAxis(i)->GetTitle());
1283  if (title=="Centrality %")
1284  contents[i] = cent;
1285  else if (title=="#it{p}_{T} (GeV/#it{c})")
1286  contents[i] = partPt;
1287  else if (title=="#eta")
1288  contents[i] = partEta;
1289  else if (title=="#phi")
1290  contents[i] = partPhi;
1291  else if (title=="Findable")
1292  contents[i] = findable;
1293  else
1294  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_PhysPrim->GetName()));
1295  }
1296 
1297  thnTracks_PhysPrim->Fill(contents);
1298 }
1299 
1300 //________________________________________________________________________
1302  Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
1303 {
1304  Double_t contents[20]={0};
1305 
1306  THnSparse* thnTracks_Matched = static_cast<THnSparse*>(fHistManager.FindObject("tracks_Matched"));
1307  if (!thnTracks_Matched) return;
1308  for (Int_t i = 0; i < thnTracks_Matched->GetNdimensions(); i++) {
1309  TString title(thnTracks_Matched->GetAxis(i)->GetTitle());
1310  if (title=="Centrality %")
1311  contents[i] = cent;
1312  else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
1313  contents[i] = partPt;
1314  else if (title=="#eta^{gen}")
1315  contents[i] = partEta;
1316  else if (title=="#phi^{gen}")
1317  contents[i] = partPhi;
1318  else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
1319  contents[i] = trackPt;
1320  else if (title=="#eta^{det}")
1321  contents[i] = trackEta;
1322  else if (title=="#phi^{det}")
1323  contents[i] = trackPhi;
1324  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
1325  contents[i] = (partPt - trackPt) / partPt;
1326  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}")
1327  contents[i] = (partPt - trackPt) / trackPt;
1328  else if (title=="track type")
1329  contents[i] = (Double_t)trackType;
1330  else
1331  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_Matched->GetName()));
1332  }
1333 
1334  thnTracks_Matched->Fill(contents);
1335 }
1336 
1343  const char* ntracks,
1344  const char* nclusters,
1345  const char* ncells,
1346  const char *nGenLev,
1347  Bool_t doTrackQA,
1348  Bool_t doCaloQA,
1349  Bool_t doJetQA,
1350  Bool_t doEventQA,
1351  Double_t trackPtCut,
1352  Double_t clusECut,
1353  const char* suffix)
1354 {
1355  // Get the pointer to the existing analysis manager via the static access method
1356  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1357  if (!mgr) {
1358  ::Error("AddTaskPWGJEQA", "No analysis manager to connect to.");
1359  return NULL;
1360  }
1361 
1362  // Check the analysis type using the event handlers connected to the analysis manager
1363  AliVEventHandler* handler = mgr->GetInputEventHandler();
1364  if (!handler) {
1365  ::Error("AddTaskPWGJEQA", "This task requires an input event handler");
1366  return NULL;
1367  }
1368 
1369  enum EDataType_t {
1370  kUnknown,
1371  kESD,
1372  kAOD
1373  };
1374 
1375  EDataType_t dataType = kUnknown;
1376 
1377  if (handler->InheritsFrom("AliESDInputHandler")) {
1378  dataType = kESD;
1379  }
1380  else if (handler->InheritsFrom("AliAODInputHandler")) {
1381  dataType = kAOD;
1382  }
1383 
1384  // Init the task and do settings
1385  TString trackName(ntracks);
1386  TString clusName(nclusters);
1387  TString cellName(ncells);
1388 
1389  if (trackName == "usedefault") {
1390  if (dataType == kESD) {
1391  trackName = "Tracks";
1392  }
1393  else if (dataType == kAOD) {
1394  trackName = "tracks";
1395  }
1396  else {
1397  trackName = "";
1398  }
1399  }
1400 
1401  if (clusName == "usedefault") {
1402  if (dataType == kESD) {
1403  clusName = "CaloClusters";
1404  }
1405  else if (dataType == kAOD) {
1406  clusName = "caloClusters";
1407  }
1408  else {
1409  clusName = "";
1410  }
1411  }
1412 
1413  if (cellName == "usedefault") {
1414  if (dataType == kESD) {
1415  cellName = "EMCALCells";
1416  }
1417  else if (dataType == kAOD) {
1418  cellName = "emcalCells";
1419  }
1420  else {
1421  cellName = "";
1422  }
1423  }
1424 
1425  TString name("AliAnalysisTaskPWGJEQA");
1426  if (!trackName.IsNull()) {
1427  name += "_";
1428  name += trackName;
1429  }
1430  if (!clusName.IsNull()) {
1431  name += "_";
1432  name += clusName;
1433  }
1434 
1435  if (!cellName.IsNull()) {
1436  name += "_";
1437  name += cellName;
1438  }
1439 
1440  if (strcmp(suffix,"")) {
1441  name += "_";
1442  name += suffix;
1443  }
1444 
1445  if (nGenLev && strcmp(nGenLev,"")!=0) cout<<"MC Part: "<< nGenLev<<endl;
1446  else cout<<"No MC Part: "<< nGenLev<<endl;
1447 
1448  AliAnalysisTaskPWGJEQA* qaTask = new AliAnalysisTaskPWGJEQA(name);
1449  qaTask->SetVzRange(-10,10);
1450  qaTask->SetNeedEmcalGeom(kFALSE);
1451  qaTask->SetCaloCellsName(cellName);
1452  qaTask->SetDetectorLevelName(trackName);
1453  if (nGenLev && strcmp(nGenLev,"")!=0) qaTask->SetGeneratorLevelName(nGenLev);
1454  qaTask->SetDoTrackQA(doTrackQA);
1455  qaTask->SetDoCaloQA(doCaloQA);
1456  qaTask->SetDoJetQA(doJetQA);
1457  qaTask->SetDoEventQA(doEventQA);
1458 
1459  // Add the detector-level track container
1460  if (trackName == "tracks" || trackName == "Tracks") {
1461  AliTrackContainer* trackCont = qaTask->AddTrackContainer(trackName);
1462  trackCont->SetFilterHybridTracks(kTRUE);
1463  }
1464  else if (!trackName.IsNull()) {
1465  qaTask->AddParticleContainer(trackName);
1466  }
1467  AliParticleContainer *partCont = qaTask->GetParticleContainer(trackName);
1468  if (partCont) {
1469  partCont->SetParticlePtCut(trackPtCut);
1470  }
1471 
1472  // Add the generator-level container, if specified
1473  if (nGenLev && strcmp(nGenLev,"")!=0) {
1474  AliMCParticleContainer* mcpartCont = qaTask->AddMCParticleContainer(nGenLev);
1475  mcpartCont->SelectPhysicalPrimaries(kTRUE);
1476  mcpartCont->SetParticlePtCut(0);
1477  }
1478 
1479  // Add the cluster container
1480  AliClusterContainer *clusterCont = qaTask->AddClusterContainer(clusName);
1481  if (clusterCont) {
1482  clusterCont->SetClusECut(0.);
1483  clusterCont->SetClusPtCut(0.);
1484  clusterCont->SetClusHadCorrEnergyCut(clusECut);
1485  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1486  }
1487 
1488  // Final settings, pass to manager and set the containers
1489  mgr->AddTask(qaTask);
1490 
1491  // Create containers for input/output
1492  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
1493 
1494  TString contName = TString::Format("%s_histos", name.Data());
1495  TString commonoutput;
1496  commonoutput = mgr->GetCommonFileName();
1497  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(),
1498  TList::Class(),AliAnalysisManager::kOutputContainer,
1499  commonoutput);
1500  mgr->ConnectInput (qaTask, 0, cinput1 );
1501  mgr->ConnectOutput (qaTask, 1, coutput1 );
1502 
1503  return qaTask;
1504 
1505 }
Int_t pdg
Double_t * fEtaHistBins
! eta bins
Double_t * fPtResHistBins
! pt res bins
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
TObjArray fClusterCollArray
cluster collection array
void SetParticlePtCut(Double_t cut)
Bool_t fDoEventQA
Set whether to enable event QA.
Double_t GetRhoVal() const
const TString & GetRhoName() const
Int_t fNEtaHistBins
! number of eta bins
double Double_t
Definition: External.C:58
static AliAnalysisTaskPWGJEQA * AddTaskPWGJEQA(const char *ntracks="usedefault", const char *nclusters="usedefault", const char *ncells="usedefault", const char *nGenLev="mcparticles", Bool_t doTrackQA=kTRUE, Bool_t doCaloQA=kTRUE, Bool_t doJetQA=kTRUE, Bool_t doEventQA=kTRUE, Double_t trackPtCut=0.15, Double_t clusECut=0.30, const char *suffix="")
const char * title
Definition: MakeQAPdf.C:27
Int_t fNPtResHistBins
! number of pt res bins
AliTrackContainer * fDetectorLevel
! detector level container
Int_t fNPhiHistBins
! number of phi bins
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
const AliMCParticleIterableMomentumContainer accepted_momentum() const
bidirectional stl iterator over the EMCAL iterable container
Container with name, TClonesArray and cuts for particles.
energy
Definition: HFPtSpectrum.C:44
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
AliAnalysisTaskPWGJEQA()
Default constructor for ROOT I/O purposes.
Bool_t fGeneralHistograms
whether or not it should fill some general 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.
THistManager fHistManager
Histogram manager.
Bool_t fRejectOutlierEvents
flag to reject pythia pt-hard jet outlier events
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Int_t fNTotClusters[3]
!Total number of accepted clusters in current event (DCal/EMCal)
TString fDetectorLevelName
detector level container name
void SetVzRange(Double_t min, Double_t max)
Bool_t FillHistograms()
Function filling histograms.
Bool_t fDoTrackQA
Set whether to enable track QA.
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
Container for particles within the EMCAL framework.
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Double_t * sigma
const Int_t nPtBins
Double_t * fPtHistBins
! pt bins
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) 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
void CreateTProfile(const char *name, const char *title, int nbinsX, double xmin, double xmax, Option_t *opt="")
Create a new TProfile within the container.
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
float Float_t
Definition: External.C:68
AliEMCALGeometry * fGeom
!emcal geometry
AliPHOSGeometry * fPHOSGeo
! phos geometry
void FillDetectorLevelTHnSparse(Double_t cent, Double_t trackEta, Double_t trackPhi, Double_t trackPt, Double_t sigma1OverPt, Byte_t trackType)
virtual AliAODMCParticle * GetAcceptMCParticleWithLabel(Int_t lab)
void SetFilterHybridTracks(Bool_t f)
Int_t fNTotTracks
!Total number of accepted tracks in current event
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Float_t fCellEnergyCut
Energy cell cut.
BeamType fForceBeamType
forced beam type
virtual Bool_t AcceptCluster(Int_t i, UInt_t &rejectionReason) const
Double_t * fPtRelDiffHistBins
! pt relative difference bins
void FillProfile(const char *name, double x, double y, double weight=1.)
Declaration of class AliAnalysisTaskPWGJEQA.
Double_t fCent
!event centrality
Int_t fNPtHistBins
! number of pt bins
TString fCaloCellsName
name of calo cell collection
void SetGeneratorLevelName(const char *name)
Int_t fNIntegerHistBins
! number of integer bins
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
Loading PYTHIA information from external cross section file into the task.
TObjArray fJetCollArray
jet collection array
short Short_t
Definition: External.C:23
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
Int_t fNCentHistBins
! number of cent bins
void ExecOnce()
Perform steps needed to initialize the analysis.
Double_t * fIntegerHistBins
! integer bins
const AliClusterIterableMomentumContainer all_momentum() const
void SetDetectorLevelName(const char *name)
void FillMatchedParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
Float_t fPtHard
!event -hard
void FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Byte_t findable)
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
AliTLorentzVector fLeadingTrack
!Leading track in current event
AliEmcalList * fOutput
!output list
AliMCParticleContainer * GetMCParticleContainer(Int_t i=0) const
Char_t GetTrackType(const AliVTrack *track) const
Definition: External.C:220
Bool_t fIsEsd
!whether it&#39;s an ESD analysis
void GenerateHistoBins()
Generate histogram binning arrays.
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
AliTrackContainer * GetTrackContainer(Int_t i=0) const
void SelectPhysicalPrimaries(Bool_t s)
Bool_t fDoCaloQA
Set whether to enable cell/cluster QA.
TH1 * fHistEventRejection
!book keep reasons for rejecting event
void SetMakeGeneralHistograms(Bool_t g)
This is a task used to do basic PWGJE QA on tracks, clusters, and jets.
void SetNeedEmcalGeom(Bool_t n)
Base task in the EMCAL jet framework.
Bool_t fIsPtHard
flag to enable pt-hard histos and make available outlier cut
const AliTrackIterableMomentumContainer accepted_momentum() const
Bool_t fDoJetQA
Set whether to enable jet QA.
void SetClusPtCut(Double_t cut)
TString fGeneratorLevelName
generator level container name
void SetRejectionReasonLabels(TAxis *axis)
void UserCreateOutputObjects()
Main initialization function on the worker.
Float_t fMaxPt
Histogram pt limit.
AliMCParticleContainer * fGeneratorLevel
! generator level container
const Int_t nbins
bool Bool_t
Definition: External.C:53
Double_t * fCentHistBins
! cent bins
EDataType_t
Switch for the data type.
void SetCaloCellsName(const char *n)
void SetClusECut(Double_t cut)
void SetDefaultClusterEnergy(Int_t d)
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.
Container structure for EMCAL clusters.
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:74
Container for MC-true particles within the EMCAL framework.
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Int_t fNPtRelDiffHistBins
! number of pt relative difference bins
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
Container for jet within the EMCAL jet framework.
Definition: External.C:196
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.
Double_t * fPhiHistBins
! phi bins
AliTLorentzVector fLeadingCluster[3]
!Leading cluster in current event (EMCal/DCal)
void SetClusHadCorrEnergyCut(Double_t cut)
const AliJetIterableContainer all() const
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal