AliPhysics  a4b41ad (a4b41ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
57 
58 //________________________________________________________________________
61  AliAnalysisTaskEmcalJet("AliAnalysisTaskPWGJEQA", kTRUE),
62  fCellEnergyCut(0.05),
63  fMaxPt(250),
64  fNTotTracks(0),
65  fLeadingTrack(),
66  fGeneratorLevel(0),
67  fGeneratorLevelName(),
68  fDetectorLevel(0),
69  fDetectorLevelName(),
70  fNPtHistBins(0),
71  fPtHistBins(0),
72  fNEtaHistBins(0),
73  fEtaHistBins(0),
74  fNPhiHistBins(0),
75  fPhiHistBins(0),
76  fNCentHistBins(0),
77  fCentHistBins(0),
78  fNPtRelDiffHistBins(0),
79  fPtRelDiffHistBins(0),
80  fNPtResHistBins(0),
81  fPtResHistBins(0),
82  fNIntegerHistBins(0),
83  fIntegerHistBins(0),
84  fHistManager("AliAnalysisTaskPWGJEQA"),
85  fDoTrackQA(kTRUE),
86  fDoCaloQA(kTRUE),
87  fDoJetQA(kTRUE),
88  fDoEventQA(kTRUE),
89  fRejectOutlierEvents(kFALSE),
90  fIsPtHard(kFALSE),
91  fPhosMinNcells(0),
92  fPhosMinM02(0),
93  fPHOSGeo(nullptr)
94 {
95  // Default constructor.
96  memset(fNTotClusters, 0, sizeof(Int_t)*3);
99 }
100 
101 //________________________________________________________________________
106  AliAnalysisTaskEmcalJet(name, kTRUE),
107  fCellEnergyCut(0.05),
108  fMaxPt(250),
109  fNTotTracks(0),
110  fLeadingTrack(),
111  fGeneratorLevel(0),
112  fGeneratorLevelName(),
113  fDetectorLevel(0),
114  fDetectorLevelName(),
115  fNPtHistBins(0),
116  fPtHistBins(0),
117  fNEtaHistBins(0),
118  fEtaHistBins(0),
119  fNPhiHistBins(0),
120  fPhiHistBins(0),
121  fNCentHistBins(0),
122  fCentHistBins(0),
123  fNPtRelDiffHistBins(0),
124  fPtRelDiffHistBins(0),
125  fNPtResHistBins(0),
126  fPtResHistBins(0),
127  fNIntegerHistBins(0),
128  fIntegerHistBins(0),
129  fHistManager(name),
130  fDoTrackQA(kTRUE),
131  fDoCaloQA(kTRUE),
132  fDoJetQA(kTRUE),
133  fDoEventQA(kTRUE),
134  fRejectOutlierEvents(kFALSE),
135  fIsPtHard(kFALSE),
136  fPhosMinNcells(0),
137  fPhosMinM02(0),
138  fPHOSGeo(nullptr)
139 {
140  // Standard
141  memset(fNTotClusters, 0, sizeof(Int_t)*3);
144 }
145 
146 //________________________________________________________________________
148 {
149  // Destructor
150 }
151 
152 //________________________________________________________________________
154 {
155  // Create histograms
157 
158  // Set track container pointers
161 
162  // Allocate histograms for tracks, cells, clusters, jets
168 
169  TIter next(fHistManager.GetListOfHistograms());
170  TObject* obj = 0;
171  while ((obj = next())) {
172  fOutput->Add(obj);
173  }
174 
175  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
176 }
177 
178 //________________________________________________________________________
180 
182  if (fGeneratorLevel) {
185  }
186 }
187 
188 //________________________________________________________________________
190 
191  TString histname;
192  TString title;
193 
194  if (!fCaloCellsName.IsNull()) {
195 
196  histname = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
197  title = histname + ";#it{E}_{cell} (GeV);counts";
198  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,0, 100);
199 
200  histname = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
201  title = histname + ";cell absId;<#it{E}_{cell}> (GeV)";
202  fHistManager.CreateTProfile(histname.Data(), title.Data(), 18000,0,18000);
203 
204  histname = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
205  title = histname + ";#it{t}_{cell} (s);counts";
206  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,-5e-6, 5e-6);
207 
208  histname = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
209  title = histname + ";cell absId;<#it{t}_{cell}> (s)";
210  fHistManager.CreateTProfile(histname.Data(), title.Data(), 18000,0,18000);
211 
212  histname = TString::Format("%s/fHistCellEvsTime", fCaloCellsName.Data());
213  title = histname + ";<#it{t}_{cell}> (s);<#it{E}_{cell}> (GeV)";
214  fHistManager.CreateTH2(histname.Data(), title.Data(), 1000, -1e-6, 1e-6, 200, 0, 100);
215 
216  }
217 }
218 
219 //________________________________________________________________________
221 
222  AliEmcalContainer* cont = 0;
223  TString histname;
224  TString title;
225  Int_t nPtBins = 2*fMaxPt;
226 
227  TIter nextClusColl(&fClusterCollArray);
228  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
229 
230  histname = TString::Format("%s/fHistClusterRejectionReason", cont->GetArrayName().Data());
231  title = histname + ";Rejection reason;#it{E}_{clus} (GeV/);counts";
232  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
233  SetRejectionReasonLabels(hist->GetXaxis());
234 
235  histname = TString::Format("%s/hPhosNCellsVsEnergy", cont->GetArrayName().Data());
236  title = histname + ";N cells;#it{E}_{clus} (GeV/);counts";
237  fHistManager.CreateTH2(histname.Data(), title.Data(), 40, 0, 40, nPtBins, 0, fMaxPt);
238 
239  histname = TString::Format("%s/hPhosM02VsEnergy", cont->GetArrayName().Data());
240  title = histname + ";M02 (cm);#it{E}_{clus} (GeV/);counts";
241  fHistManager.CreateTH2(histname.Data(), title.Data(), 100, 0, 20, nPtBins, 0, fMaxPt);
242 
243  histname = TString::Format("%s/hPhosCellIdVsEnergy", cont->GetArrayName().Data());
244  title = histname + ";Cell ID;#it{E}_{clus} (GeV/);counts";
245  fHistManager.CreateTH2(histname.Data(), title.Data(), 20000, 0, 20000, nPtBins, 0, fMaxPt);
246 
247  const Int_t nEmcalSM = 20;
248  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
249  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
250  title = histname + ";#it{E}_{cluster} (GeV);counts";
251  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
252  }
253 
254  for (Int_t sm = 1; sm < 5; sm++) {
255  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
256  title = histname + ";#it{E}_{cluster} (GeV);counts";
257  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
258  }
259 
260  // Cluster THnSparse
261  Int_t dim = 0;
262  TString title[20];
263  Int_t nbins[20] = {0};
264  Double_t min[30] = {0.};
265  Double_t max[30] = {0.};
266  Double_t *binEdges[20] = {0};
267 
269  title[dim] = "Centrality %";
270  nbins[dim] = fNCentHistBins;
271  binEdges[dim] = fCentHistBins;
272  min[dim] = fCentHistBins[0];
273  max[dim] = fCentHistBins[fNCentHistBins];
274  dim++;
275  }
276 
277  title[dim] = "#it{E}_{clus} (GeV)";
278  nbins[dim] = fNPtHistBins;
279  binEdges[dim] = fPtHistBins;
280  min[dim] = fPtHistBins[0];
281  max[dim] = fPtHistBins[fNPtHistBins];
282  dim++;
283 
284  title[dim] = "#eta";
285  nbins[dim] = fNEtaHistBins;
286  binEdges[dim] = fEtaHistBins;
287  min[dim] = fEtaHistBins[0];
288  max[dim] = fEtaHistBins[fNEtaHistBins];
289  dim++;
290 
291  title[dim] = "#phi";
292  nbins[dim] = fNPhiHistBins;
293  binEdges[dim] = fPhiHistBins;
294  min[dim] = fPhiHistBins[0];
295  max[dim] = fPhiHistBins[fNPhiHistBins];
296  dim++;
297 
298  title[dim] = "cluster type";
299  nbins[dim] = 3;
300  binEdges[dim] = fIntegerHistBins;
301  min[dim] = fIntegerHistBins[0];
302  max[dim] = fIntegerHistBins[3];
303  dim++;
304 
305  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
306  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
307  for (Int_t i = 0; i < dim; i++) {
308  hn->GetAxis(i)->SetTitle(title[i]);
309  hn->SetBinEdges(i, binEdges[i]);
310  }
311  }
312 }
313 
314 //________________________________________________________________________
316 
317  AliJetContainer* jets = 0;
318  TIter nextJetColl(&fJetCollArray);
319  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
320 
321  // Allocate THnSparse
322  Double_t jetRadius = jets->GetJetRadius();
323 
324  TString axisTitle[30]= {""};
325  Int_t nbins[30] = {0};
326  Double_t min[30] = {0.};
327  Double_t max[30] = {0.};
328  Double_t *binEdges[20] = {0};
329  Int_t dim = 0;
330 
331  if (fForceBeamType != kpp) {
332  axisTitle[dim] = "Centrality (%)";
333  nbins[dim] = fNCentHistBins;
334  binEdges[dim] = fCentHistBins;
335  min[dim] = fCentHistBins[0];
336  max[dim] = fCentHistBins[fNCentHistBins];
337  dim++;
338  }
339 
340  axisTitle[dim] = "#eta_{jet}";
341  nbins[dim] = fNEtaHistBins;
342  binEdges[dim] = fEtaHistBins;
343  min[dim] = fEtaHistBins[0];
344  max[dim] = fEtaHistBins[fNEtaHistBins];
345  dim++;
346 
347  axisTitle[dim] = "#phi_{jet} (rad)";
348  nbins[dim] = fNPhiHistBins;
349  binEdges[dim] = fPhiHistBins;
350  min[dim] = fPhiHistBins[0];
351  max[dim] = fPhiHistBins[fNPhiHistBins];
352  dim++;
353 
354  axisTitle[dim] = "#it{p}_{T} (GeV/#it{c})";
355  nbins[dim] = TMath::CeilNint(fMaxPt/3);
356  min[dim] = 0;
357  max[dim] = fMaxPt;
358  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
359  dim++;
360 
361  if (fForceBeamType != kpp) {
362  axisTitle[dim] = "#it{p}_{T}^{corr} (GeV/#it{c})";
363  nbins[dim] = TMath::CeilNint(fMaxPt/3);
364  min[dim] = -fMaxPt/2 + 25;
365  max[dim] = fMaxPt/2 + 25;
366  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
367  dim++;
368  }
369 
370  axisTitle[dim] = "#it{p}_{T,particle}^{leading} (GeV/#it{c})";
371  nbins[dim] = fMaxPt/5;
372  min[dim] = 0;
373  max[dim] = 150;
374  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
375  dim++;
376 
377  TString thnname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
378  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
379  for (Int_t i = 0; i < dim; i++) {
380  hn->GetAxis(i)->SetTitle(axisTitle[i]);
381  hn->SetBinEdges(i, binEdges[i]);
382  }
383 
384  // Allocate other jet histograms
385  TString histname;
386  TString title;
387  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
388  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
389  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
390  SetRejectionReasonLabels(hist->GetXaxis());
391 
392  if (!jets->GetRhoName().IsNull()) {
393  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
394  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
395  fHistManager.CreateTH2(histname.Data(), title.Data(), 101, 0, 101, 100, 0, 500);
396  }
397  }
398 }
399 
400 //________________________________________________________________________
402 
403  Int_t nPtBins = 2*fMaxPt;
404 
405  Int_t dim = 0;
406  TString axistitle[40];
407  Int_t nbins[40] = {0};
408  Double_t min[40] = {0};
409  Double_t max[40] = {0};
410  Double_t *binEdges[20] = {0};
411 
413  axistitle[dim] = "Centrality %";
414  nbins[dim] = fNCentHistBins;
415  binEdges[dim] = fCentHistBins;
416  min[dim] = fCentHistBins[0];
417  max[dim] = fCentHistBins[fNCentHistBins];
418  dim++;
419  }
420 
421  if (fParticleCollArray.GetEntriesFast()>0) {
422  axistitle[dim] = "No. of tracks";
424  nbins[dim] = 50;
425  min[dim] = -0.5;
426  max[dim] = 10000-0.5;
427  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
428  }
429  else {
430  nbins[dim] = 50;
431  min[dim] = 0;
432  max[dim] = 200;
433  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
434  }
435  dim++;
436 
437  axistitle[dim] = "#it{p}_{T,track}^{leading} (GeV/c)";
438  nbins[dim] = fMaxPt/5;
439  min[dim] = 0;
440  max[dim] = 150;
441  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
442  dim++;
443  }
444 
445  if (fClusterCollArray.GetEntriesFast()>0) {
446  axistitle[dim] = "No. of clusters";
447 
449  nbins[dim] = 50;
450  min[dim] = -0.5;
451  max[dim] = 4000-0.5;
452  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
453  }
454  else {
455  nbins[dim] = 50;
456  min[dim] = 0;
457  max[dim] = 200;
458  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
459  }
460  dim++;
461 
462  axistitle[dim] = "#it{E}_{cluster}^{leading} (GeV)";
463  nbins[dim] = fMaxPt/5;
464  min[dim] = 0;
465  max[dim] = 150;
466  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
467  dim++;
468  }
469 
470  THnSparse* hn = fHistManager.CreateTHnSparse("eventQA","eventQA",dim,nbins,min,max);
471  for (Int_t i = 0; i < dim; i++) {
472  hn->GetAxis(i)->SetTitle(axistitle[i]);
473  hn->SetBinEdges(i, binEdges[i]);
474  }
475 
476  if (fIsPtHard) {
477  TString histname = "hPtHard";
478  TString title = histname + ";#it{p}_{T,hard} (GeV/c);counts";
479  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
480 
481  TH1* hTrials = fHistManager.CreateTH1("hNtrials", "hNtrials", 1, 0, 1);
482  hTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
483 
484  TH1* hxsection = fHistManager.CreateTH1("hXsec", "hXsec", 1, 0, 1);
485  hxsection->GetXaxis()->SetBinLabel(1,"<#sigma>");
486  }
487 
488  fHistEventRejection->GetXaxis()->SetBinLabel(15,"PtHardBinJetOutlier");
489 
490 }
491 
492 //________________________________________________________________________
494 {
495  Int_t dim = 0;
496  TString title[20];
497  Int_t nbins[20] = {0};
498  Double_t min[30] = {0.};
499  Double_t max[30] = {0.};
500  Double_t *binEdges[20] = {0};
501 
503  title[dim] = "Centrality %";
504  nbins[dim] = fNCentHistBins;
505  binEdges[dim] = fCentHistBins;
506  min[dim] = fCentHistBins[0];
507  max[dim] = fCentHistBins[fNCentHistBins];
508  dim++;
509  }
510 
511  title[dim] = "#it{p}_{T} (GeV/#it{c})";
512  nbins[dim] = fNPtHistBins;
513  binEdges[dim] = fPtHistBins;
514  min[dim] = fPtHistBins[0];
515  max[dim] = fPtHistBins[fNPtHistBins];
516  dim++;
517 
518  title[dim] = "#eta";
519  nbins[dim] = fNEtaHistBins;
520  binEdges[dim] = fEtaHistBins;
521  min[dim] = fEtaHistBins[0];
522  max[dim] = fEtaHistBins[fNEtaHistBins];
523  dim++;
524 
525  title[dim] = "#phi";
526  nbins[dim] = fNPhiHistBins;
527  binEdges[dim] = fPhiHistBins;
528  min[dim] = fPhiHistBins[0];
529  max[dim] = fPhiHistBins[fNPhiHistBins];
530  dim++;
531 
532  title[dim] = "track type";
533  nbins[dim] = 3;
534  binEdges[dim] = fIntegerHistBins;
535  min[dim] = fIntegerHistBins[0];
536  max[dim] = fIntegerHistBins[3];
537  dim++;
538 
539  title[dim] = "#sigma(#it{p}_{T}) / #it{p}_{T}";
540  nbins[dim] = fNPtResHistBins;
541  binEdges[dim] = fPtResHistBins;
542  min[dim] = fPtResHistBins[0];
543  max[dim] = fPtResHistBins[fNPtResHistBins];
544  dim++;
545 
546  THnSparse* hn = fHistManager.CreateTHnSparse("tracks","tracks", dim, nbins, min, max);
547  for (Int_t i = 0; i < dim; i++) {
548  hn->GetAxis(i)->SetTitle(title[i]);
549  hn->SetBinEdges(i, binEdges[i]);
550  }
551 
552 }
553 
554 //________________________________________________________________________
556 {
557  Int_t dim = 0;
558  TString title[20];
559  Int_t nbins[20] = {0};
560  Double_t min[30] = {0.};
561  Double_t max[30] = {0.};
562  Double_t *binEdges[20] = {0};
563 
565  title[dim] = "Centrality %";
566  nbins[dim] = fNCentHistBins;
567  binEdges[dim] = fCentHistBins;
568  min[dim] = fCentHistBins[0];
569  max[dim] = fCentHistBins[fNCentHistBins];
570  dim++;
571  }
572 
573  title[dim] = "#it{p}_{T} (GeV/#it{c})";
574  nbins[dim] = fNPtHistBins;
575  binEdges[dim] = fPtHistBins;
576  min[dim] = fPtHistBins[0];
577  max[dim] = fPtHistBins[fNPtHistBins];
578  dim++;
579 
580  title[dim] = "#eta";
581  nbins[dim] = fNEtaHistBins;
582  binEdges[dim] = fEtaHistBins;
583  min[dim] = fEtaHistBins[0];
584  max[dim] = fEtaHistBins[fNEtaHistBins];
585  dim++;
586 
587  title[dim] = "#phi";
588  nbins[dim] = fNPhiHistBins;
589  binEdges[dim] = fPhiHistBins;
590  min[dim] = fPhiHistBins[0];
591  max[dim] = fPhiHistBins[fNPhiHistBins];
592  dim++;
593 
594  title[dim] = "Findable";
595  nbins[dim] = 2;
596  binEdges[dim] = fIntegerHistBins;
597  min[dim] = fIntegerHistBins[0];
598  max[dim] = fIntegerHistBins[2];
599  dim++;
600 
601  THnSparse* hn = fHistManager.CreateTHnSparse("tracks_PhysPrim","tracks_PhysPrim", dim, nbins, min, max);
602  for (Int_t i = 0; i < dim; i++) {
603  hn->GetAxis(i)->SetTitle(title[i]);
604  hn->SetBinEdges(i, binEdges[i]);
605  }
606 }
607 
608 //________________________________________________________________________
610 {
611  Int_t dim = 0;
612  TString title[20];
613  Int_t nbins[20] = {0};
614  Double_t min[30] = {0.};
615  Double_t max[30] = {0.};
616  Double_t *binEdges[20] = {0};
617 
619  title[dim] = "Centrality %";
620  nbins[dim] = fNCentHistBins;
621  binEdges[dim] = fCentHistBins;
622  min[dim] = fCentHistBins[0];
623  max[dim] = fCentHistBins[fNCentHistBins];
624  dim++;
625  }
626 
627  title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
628  nbins[dim] = fNPtHistBins;
629  binEdges[dim] = fPtHistBins;
630  min[dim] = fPtHistBins[0];
631  max[dim] = fPtHistBins[fNPtHistBins];
632  dim++;
633 
634  title[dim] = "#eta^{gen}";
635  nbins[dim] = fNEtaHistBins;
636  binEdges[dim] = fEtaHistBins;
637  min[dim] = fEtaHistBins[0];
638  max[dim] = fEtaHistBins[fNEtaHistBins];
639  dim++;
640 
641  title[dim] = "#phi^{gen}";
642  nbins[dim] = fNPhiHistBins;
643  binEdges[dim] = fPhiHistBins;
644  min[dim] = fPhiHistBins[0];
645  max[dim] = fPhiHistBins[fNPhiHistBins];
646  dim++;
647 
648  title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
649  nbins[dim] = fNPtHistBins;
650  binEdges[dim] = fPtHistBins;
651  min[dim] = fPtHistBins[0];
652  max[dim] = fPtHistBins[fNPtHistBins];
653  dim++;
654 
655  title[dim] = "#eta^{det}";
656  nbins[dim] = fNEtaHistBins;
657  binEdges[dim] = fEtaHistBins;
658  min[dim] = fEtaHistBins[0];
659  max[dim] = fEtaHistBins[fNEtaHistBins];
660  dim++;
661 
662  title[dim] = "#phi^{det}";
663  nbins[dim] = fNPhiHistBins;
664  binEdges[dim] = fPhiHistBins;
665  min[dim] = fPhiHistBins[0];
666  max[dim] = fPhiHistBins[fNPhiHistBins];
667  dim++;
668 
669  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
670  nbins[dim] = fNPtRelDiffHistBins;
671  binEdges[dim] = fPtRelDiffHistBins;
672  min[dim] = fPtRelDiffHistBins[0];
674  dim++;
675 
676  title[dim] = "track type";
677  nbins[dim] = 3;
678  binEdges[dim] = fIntegerHistBins;
679  min[dim] = fIntegerHistBins[0];
680  max[dim] = fIntegerHistBins[3];
681  dim++;
682 
683  THnSparse* hn = fHistManager.CreateTHnSparse("tracks_Matched","tracks_Matched", dim, nbins, min, max);
684  for (Int_t i = 0; i < dim; i++) {
685  hn->GetAxis(i)->SetTitle(title[i]);
686  hn->SetBinEdges(i, binEdges[i]);
687  }
688 }
689 
690 
691 //________________________________________________________________________
694 {
695  fNPtHistBins = 82;
698  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
699  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
700  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
701  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
702  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
703  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
704 
705  fNEtaHistBins = 70;
708 
709  fNPhiHistBins = 101;
711  GenerateFixedBinArray(fNPhiHistBins, 0, TMath::Pi() * 2.02, fPhiHistBins);
712 
713  fNCentHistBins = 4;
715  fCentHistBins[0] = 0;
716  fCentHistBins[1] = 10;
717  fCentHistBins[2] = 30;
718  fCentHistBins[3] = 50;
719  fCentHistBins[4] = 90;
720 
721  fNPtResHistBins = 100;
724  GenerateFixedBinArray(25, 0.05, 0.10, fPtResHistBins+50);
725  GenerateFixedBinArray(15, 0.10, 0.20, fPtResHistBins+75);
726  GenerateFixedBinArray(10, 0.20, 0.50, fPtResHistBins+90);
727 
728  fNPtRelDiffHistBins = 50;
731 
732  fNIntegerHistBins = 3;
735 }
736 
737 //________________________________________________________________________
739 {
740  if (fClusterCollArray.GetEntriesFast() == 0 && fCaloCellsName.IsNull()) {
741  fNeedEmcalGeom = kFALSE;
742  }
743  else {
744  fNeedEmcalGeom = kTRUE;
745  }
746 
748 
749  // Load the PHOS geometry
750  fPHOSGeo = AliPHOSGeometry::GetInstance();
751  if (fPHOSGeo) {
752  AliInfo("Found instance of PHOS geometry!");
753  }
754  else {
755  AliInfo("Creating PHOS geometry!");
756  Int_t runNum = InputEvent()->GetRunNumber();
757  if(runNum<209122) //Run1
758  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
759  else
760  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
761 
762  if (fPHOSGeo) {
763  AliOADBContainer geomContainer("phosGeo");
764  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
765  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
766  for(Int_t mod=0; mod<6; mod++) {
767  if(!matrixes->At(mod)) continue;
768  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
769  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
770  ((TGeoHMatrix*)matrixes->At(mod))->Print();
771  }
772  }
773  }
774 
775 }
776 
777 //________________________________________________________________________
779 {
780  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
781 
782  // If Pt-hard production, get the Pt-hard of the event, and have possibility to reject the event for jet outliers
783  if (fIsPtHard) {
784  AliGenPythiaEventHeader* pygen;
785  if (MCEvent()) {
786  pygen = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
787  if (!pygen) {
788  // Check if AOD
789  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
790  if (aodMCH) {
791  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
792  pygen = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
793  if (pygen) break;
794  }
795  }
796  }
797  }
798  if (pygen) {
799  fPtHard = pygen->GetPtHard();
800  //fNTrials = pygen->Trials();
801  //fXsection = pygen->GetXsection();
802 
803  // reject outlier events, where the jet Pt is much larger than Pt-hard
804  if (fRejectOutlierEvents) {
805  AliTLorentzVector jet;
806  Int_t nTriggerJets = pygen->NTriggerJets();
807  Float_t tmpjet[]={0,0,0,0};
808  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
809  pygen->TriggerJet(ijet, tmpjet);
810  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
811 
812  if (jet.Pt() > 4. * fPtHard) {
813  //AliInfo(Form("Reject jet event with: pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), 4.));
814  if (fGeneralHistograms) fHistEventRejection->Fill("PtHardBinJetOutlier",1);
815  return kFALSE;
816  }
817  }
818  }
819  }
820  }
821 
822  return kTRUE;
823 }
824 
832 //________________________________________________________________________
834 {
835  if (fIsPtHard) {
836  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
837  if (!tree) {
838  AliError(Form("%s - UserNotify: No current tree!",GetName()));
839  return kFALSE;
840  }
841 
842  Float_t xsection = 0;
843  Float_t trials = 0;
844  Int_t pthardbin = 0;
845 
846  TFile *curfile = tree->GetCurrentFile();
847  if (!curfile) {
848  AliError(Form("%s - UserNotify: No current file!",GetName()));
849  return kFALSE;
850  }
851 
852  TChain *chain = dynamic_cast<TChain*>(tree);
853  if (chain) tree = chain->GetTree();
854 
855  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
856 
857  fHistManager.FillTH1("hNtrials", "#sum{ntrials}", trials);
858  fHistManager.FillTH1("hXsec", "<#sigma>", xsection);
859  }
860  return kTRUE;
861 }
862 
863 //________________________________________________________________________
865 {
871 
872  return kTRUE;
873 }
874 
875 //________________________________________________________________________
877 
878  fNTotTracks = 0;
879  fLeadingTrack.SetPxPyPzE(0,0,0,0);
880 
881  AliVTrack* track;
882  for (auto trackIterator : fDetectorLevel->accepted_momentum() ) {
883 
884  fNTotTracks++;
885  if (fLeadingTrack.Pt() < trackIterator.first.Pt()) fLeadingTrack = trackIterator.first;
886 
887  track = trackIterator.second;
888  Byte_t type = fDetectorLevel->GetTrackType(track);
889  if (type <= 2) {
890  Double_t sigma = 0;
891 
892  if (fIsEsd) {
893 
894  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
895  if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
896 
897  } else { // AOD
898 
899  AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
900  if(!aodtrack) AliFatal("Not a standard AOD");
901 
902  AliExternalTrackParam exParam;
903 
904  //get covariance matrix
905  Double_t cov[21] = {0,};
906  aodtrack->GetCovMatrix(cov);
907  Double_t pxpypz[3] = {0,};
908  aodtrack->PxPyPz(pxpypz);
909  Double_t xyz[3] = {0,};
910  aodtrack->GetXYZ(xyz);
911  Short_t sign = aodtrack->Charge();
912  exParam.Set(xyz,pxpypz,cov,sign);
913 
914  sigma = TMath::Sqrt(exParam.GetSigma1Pt2());
915 
916  }
917 
918  Int_t label = TMath::Abs(track->GetLabel());
919 
920  FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, type);
921 
922  if (fGeneratorLevel && label > 0) {
923  AliAODMCParticle *part = fGeneratorLevel->GetAcceptMCParticleWithLabel(label);
924  if (part) {
925  if (part->GetGeneratorIndex() == 0) {
926  Int_t pdg = TMath::Abs(part->PdgCode());
927  // select charged pions, protons, kaons, electrons, muons
928  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
929  FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
930  }
931  }
932  }
933  }
934  }
935  else {
936  AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
937  }
938  }
939 
940  if (fGeneratorLevel) {
941  AliAODMCParticle* part;
942  for (auto partIterator : fGeneratorLevel->accepted_momentum() ) {
943  part = partIterator.second;
944 
945  Byte_t findable = 0;
946 
947  Int_t pdg = TMath::Abs(part->PdgCode());
948  // select charged pions, protons, kaons, electrons, muons
949  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
950 
951  FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), findable);
952  }
953  }
954 }
955 
956 //________________________________________________________________________
958 
959  if (!fCaloCells) return;
960 
961  TString histname_energy = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
962  TString histname_time = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
963  TString histname_energyProf = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
964  TString histname_timeProf = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
965  TString histname_EnergyvsTime = TString::Format("%s/fHistCellEvsTime", fCaloCellsName.Data());
966 
967  const Int_t ncells = fCaloCells->GetNumberOfCells();
968  for (Int_t pos = 0; pos < ncells; pos++) {
969  Float_t amp = fCaloCells->GetAmplitude(pos);
970  Float_t time = fCaloCells->GetTime(pos);
971  Int_t absId = fCaloCells->GetCellNumber(pos);
972 
973  if (amp < fCellEnergyCut) continue;
974 
975  fHistManager.FillTH1(histname_energy, amp);
976  fHistManager.FillTH1(histname_time, time);
977 
978  fHistManager.FillProfile(histname_energyProf, absId, amp);
979  fHistManager.FillProfile(histname_timeProf, absId, time);
980 
981  fHistManager.FillTH2(histname_EnergyvsTime, time, amp);
982  }
983 }
984 
985 //________________________________________________________________________
987 
988  memset(fNTotClusters, 0, sizeof(Int_t)*3);
989  for (Int_t i = 0; i < 3; i++) fLeadingCluster[i].SetPxPyPzE(0,0,0,0);
990 
991  TString histname;
992  AliClusterContainer* clusters = 0;
993  TIter nextClusColl(&fClusterCollArray);
994  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
995  // Cluster loop
997  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
998 
999  UInt_t rejectionReason = 0;
1000  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
1001  histname = TString::Format("%s/fHistClusterRejectionReason", clusters->GetArrayName().Data());
1002  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
1003  continue;
1004  }
1005 
1006  // Determine cluster type (EMCal/DCal/Phos) and record relevant eventQA info
1007  ClusterType clusType = kNA;
1008  if (it->second->IsEMCAL()) {
1009  Double_t phi = it->first.Phi_0_2pi();
1010  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
1011  if (isDcal == 0) {
1012  clusType = kEMCal;
1013  } else if (isDcal == 1) {
1014  clusType = kDCal;
1015  }
1016 
1017  if (fLeadingCluster[isDcal].E() < it->first.E()) fLeadingCluster[isDcal] = it->first;
1018  fNTotClusters[isDcal]++;
1019 
1020  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
1021  if (sm >=0 && sm < 20) {
1022  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1023  fHistManager.FillTH1(histname, it->second->E());
1024  }
1025  else {
1026  AliError(Form("Supermodule %d does not exist!", sm));
1027  }
1028 
1029  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
1030 
1031  Int_t nCells = it->second->GetNCells();
1032  Double_t M02 = it->second->GetM02();
1033  Double_t energy = it->second->E();
1034 
1035  histname = TString::Format("%s/hPhosNCellsVsEnergy", clusters->GetArrayName().Data());
1036  fHistManager.FillTH2(histname, nCells, energy);
1037 
1038  histname = TString::Format("%s/hPhosM02VsEnergy", clusters->GetArrayName().Data());
1039  fHistManager.FillTH2(histname, M02, energy);
1040 
1041  // Apply recommended exotic cuts
1042  if (nCells < fPhosMinNcells) continue;
1043  if (M02 < fPhosMinM02) continue;
1044 
1045  clusType = kPHOS;
1046  fNTotClusters[2]++;
1047  if (fLeadingCluster[2].E() < it->first.E()) fLeadingCluster[2] = it->first;
1048 
1049  // Fill Phos spectra by module
1050  Int_t relid[4];
1051  if (fPHOSGeo) {
1052  fPHOSGeo->AbsToRelNumbering(it->second->GetCellAbsId(0), relid);
1053  Int_t sm = relid[0];
1054  if (sm >=1 && sm < 5) {
1055  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1056  fHistManager.FillTH1(histname, energy);
1057  }
1058  else {
1059  AliError(Form("Supermodule %d does not exist!", sm));
1060  }
1061  }
1062 
1063  // Loop through cells in the cluster, and histogram their energy
1064  histname = TString::Format("%s/hPhosCellIdVsEnergy", clusters->GetArrayName().Data());
1065  for (Int_t i=0; i < nCells; i++) {
1066  fHistManager.FillTH2(histname, it->second->GetCellAbsId(i), energy);
1067  }
1068 
1069  }
1070 
1071  Double_t contents[30]={0};
1072  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
1073  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1074  if (!histClusterObservables) return;
1075  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1076  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1077  if (title=="Centrality %")
1078  contents[i] = fCent;
1079  else if (title=="#it{E}_{clus} (GeV)")
1080  contents[i] = it->first.E();
1081  else if (title=="#eta")
1082  contents[i] = it->first.Eta();
1083  else if (title=="#phi")
1084  contents[i] = it->first.Phi_0_2pi();
1085  else if (title=="cluster type")
1086  contents[i] = clusType;
1087  else
1088  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1089  }
1090  histClusterObservables->Fill(contents);
1091 
1092  }
1093  }
1094 }
1095 
1096 //________________________________________________________________________
1098 
1099  TString histname;
1100  AliJetContainer* jets = 0;
1101  TIter nextJetColl(&fJetCollArray);
1102  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1103  Double_t rhoVal = 0;
1104  if (jets->GetRhoParameter()) {
1105  rhoVal = jets->GetRhoVal();
1106  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
1107  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1108  }
1109 
1110  for (auto jet : jets->all()) {
1111 
1112  UInt_t rejectionReason = 0;
1113  if (!jets->AcceptJet(jet, rejectionReason)) {
1114  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
1115  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1116  continue;
1117  }
1118 
1119  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1120  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
1121 
1122  TLorentzVector leadPart;
1123  jets->GetLeadingHadronMomentum(leadPart, jet);
1124 
1125  Double_t contents[30]={0};
1126  histname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
1127  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1128  if (!histJetObservables) return;
1129  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
1130  TString title(histJetObservables->GetAxis(i)->GetTitle());
1131  if (title=="Centrality (%)")
1132  contents[i] = fCent;
1133  else if (title=="#eta_{jet}")
1134  contents[i] = jet->Eta();
1135  else if (title=="#phi_{jet} (rad)")
1136  contents[i] = jet->Phi_0_2pi();
1137  else if (title=="#it{p}_{T} (GeV/#it{c})")
1138  contents[i] = jet->Pt();
1139  else if (title=="#it{p}_{T}^{MC} (GeV/#it{c})")
1140  contents[i] = jet->MCPt();
1141  else if (title=="#it{p}_{T}^{corr} (GeV/#it{c})")
1142  contents[i] = corrPt;
1143  else if (title=="#it{p}_{T,particle}^{leading} (GeV/#it{c})")
1144  contents[i] = ptLeading;
1145  else
1146  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1147  }
1148  histJetObservables->Fill(contents);
1149 
1150  } //jet loop
1151  }
1152 }
1153 
1154 //________________________________________________________________________
1156 
1157  EventQA_t eventQA;
1158  for (Int_t i = 0; i < 3; i++) {
1159  eventQA.fMaxCluster[i] = fLeadingCluster[i];
1160  }
1161  eventQA.fCent = fCent;
1162  eventQA.fNTracks = fNTotTracks;
1163  eventQA.fNClusters[0] = fNTotClusters[0];
1164  eventQA.fNClusters[1] = fNTotClusters[1];
1165  eventQA.fNClusters[2] = fNTotClusters[2];
1166  eventQA.fMaxTrack = fLeadingTrack;
1167 
1168  Int_t globalNclusters = eventQA.fNClusters[0] + eventQA.fNClusters[1] + eventQA.fNClusters[2];
1169  AliTLorentzVector globalMaxCluster;
1170  for (Int_t i = 0; i < 3; i++) {
1171  if (globalMaxCluster.E() < eventQA.fMaxCluster[i].E()) globalMaxCluster = eventQA.fMaxCluster[i];
1172  }
1173 
1174  THnSparse* histEventQA = static_cast<THnSparse*>(fHistManager.FindObject("eventQA"));
1175  Double_t contents[40]={0};
1176  for (Int_t i = 0; i < histEventQA->GetNdimensions(); i++) {
1177  TString title(histEventQA->GetAxis(i)->GetTitle());
1178  if (title=="Centrality %")
1179  contents[i] = eventQA.fCent;
1180  else if (title=="No. of tracks")
1181  contents[i] = eventQA.fNTracks;
1182  else if (title=="No. of clusters")
1183  contents[i] = globalNclusters;
1184  else if (title=="#it{p}_{T,track}^{leading} (GeV/c)")
1185  contents[i] = eventQA.fMaxTrack.Pt();
1186  else if (title=="#it{E}_{cluster}^{leading} (GeV)")
1187  contents[i] = globalMaxCluster.E();
1188  else
1189  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1190  }
1191 
1192  histEventQA->Fill(contents);
1193 
1194  // Fill pythia pt hard histograms, if applicable
1195  if (fPtHard > 1e-6) {
1196  fHistManager.FillTH1("hPtHard", fPtHard);
1197  }
1198 }
1199 
1200 //________________________________________________________________________
1202  Double_t sigma1OverPt, Byte_t trackType)
1203 {
1204  Double_t contents[20]={0};
1205 
1206  THnSparse* thnTracks = static_cast<THnSparse*>(fHistManager.FindObject("tracks"));
1207  if (!thnTracks) return;
1208  for (Int_t i = 0; i < thnTracks->GetNdimensions(); i++) {
1209  TString title(thnTracks->GetAxis(i)->GetTitle());
1210  if (title=="Centrality %")
1211  contents[i] = cent;
1212  else if (title=="#it{p}_{T} (GeV/#it{c})")
1213  contents[i] = trackPt;
1214  else if (title=="#eta")
1215  contents[i] = trackEta;
1216  else if (title=="#phi")
1217  contents[i] = trackPhi;
1218  else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
1219  contents[i] = sigma1OverPt;
1220  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
1221  contents[i] = sigma1OverPt*trackPt;
1222  else if (title=="track type")
1223  contents[i] = trackType;
1224  else
1225  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks->GetName()));
1226  }
1227 
1228  thnTracks->Fill(contents);
1229 }
1230 
1231 //________________________________________________________________________
1232 void AliAnalysisTaskPWGJEQA::FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Byte_t findable)
1233 {
1234  Double_t contents[20]={0};
1235 
1236  THnSparse* thnTracks_PhysPrim = static_cast<THnSparse*>(fHistManager.FindObject("tracks_PhysPrim"));
1237  if (!thnTracks_PhysPrim) return;
1238  for (Int_t i = 0; i < thnTracks_PhysPrim->GetNdimensions(); i++) {
1239  TString title(thnTracks_PhysPrim->GetAxis(i)->GetTitle());
1240  if (title=="Centrality %")
1241  contents[i] = cent;
1242  else if (title=="#it{p}_{T} (GeV/#it{c})")
1243  contents[i] = partPt;
1244  else if (title=="#eta")
1245  contents[i] = partEta;
1246  else if (title=="#phi")
1247  contents[i] = partPhi;
1248  else if (title=="Findable")
1249  contents[i] = findable;
1250  else
1251  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_PhysPrim->GetName()));
1252  }
1253 
1254  thnTracks_PhysPrim->Fill(contents);
1255 }
1256 
1257 //________________________________________________________________________
1259  Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
1260 {
1261  Double_t contents[20]={0};
1262 
1263  THnSparse* thnTracks_Matched = static_cast<THnSparse*>(fHistManager.FindObject("tracks_Matched"));
1264  if (!thnTracks_Matched) return;
1265  for (Int_t i = 0; i < thnTracks_Matched->GetNdimensions(); i++) {
1266  TString title(thnTracks_Matched->GetAxis(i)->GetTitle());
1267  if (title=="Centrality %")
1268  contents[i] = cent;
1269  else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
1270  contents[i] = partPt;
1271  else if (title=="#eta^{gen}")
1272  contents[i] = partEta;
1273  else if (title=="#phi^{gen}")
1274  contents[i] = partPhi;
1275  else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
1276  contents[i] = trackPt;
1277  else if (title=="#eta^{det}")
1278  contents[i] = trackEta;
1279  else if (title=="#phi^{det}")
1280  contents[i] = trackPhi;
1281  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
1282  contents[i] = (partPt - trackPt) / partPt;
1283  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}")
1284  contents[i] = (partPt - trackPt) / trackPt;
1285  else if (title=="track type")
1286  contents[i] = (Double_t)trackType;
1287  else
1288  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_Matched->GetName()));
1289  }
1290 
1291  thnTracks_Matched->Fill(contents);
1292 }
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
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
const char * title
Definition: MakeQAPdf.C:26
Int_t fNPtResHistBins
! number of pt res bins
AliTrackContainer * fDetectorLevel
! detector level container
Int_t fNPhiHistBins
! number of phi bins
const AliMCParticleIterableMomentumContainer accepted_momentum() const
bidirectional stl iterator over the EMCAL iterable container
energy
Definition: HFPtSpectrum.C:44
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Double_t fPtHard
!event pt hard
AliAnalysisTaskPWGJEQA()
Default constructor for ROOT I/O purposes.
Bool_t fGeneralHistograms
whether or not it should fill some general histograms
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)
Int_t fPhosMinNcells
min number of phos cells per cluster
TString fDetectorLevelName
detector level container name
Bool_t fDoTrackQA
Set whether to enable track QA.
TObjArray fParticleCollArray
particle/track collection array
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="")
TObject * FindObject(const char *name) const
int Int_t
Definition: External.C:63
void CreateTProfile(const char *name, const char *title, int nbinsX, double xmin, double xmax, Option_t *opt="")
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
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)
Int_t fNTotTracks
!Total number of accepted tracks in current event
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
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
Int_t fNIntegerHistBins
! number of integer bins
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
TObjArray fJetCollArray
jet collection array
short Short_t
Definition: External.C:23
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
Int_t fNCentHistBins
! number of cent bins
Double_t * fIntegerHistBins
! integer bins
const AliClusterIterableMomentumContainer all_momentum() const
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)
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
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
AliMCParticleContainer * GetMCParticleContainer(Int_t i=0) const
Char_t GetTrackType(const AliVTrack *track) const
Definition: External.C:220
Double_t fPhosMinM02
min value of M02 for phos clusters
Bool_t fIsEsd
!whether it's an ESD analysis
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
void GenerateHistoBins()
Generate histogram binning arrays.
AliTrackContainer * GetTrackContainer(Int_t i=0) const
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.
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.
TString fGeneratorLevelName
generator level container name
void SetRejectionReasonLabels(TAxis *axis)
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
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Container structure for EMCAL clusters.
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Int_t fNPtRelDiffHistBins
! number of pt relative difference bins
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Double_t * fPhiHistBins
! phi bins
AliTLorentzVector fLeadingCluster[3]
!Leading cluster in current event (EMCal/DCal)
const AliJetIterableContainer all() const
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal