AliPhysics  d9e9949 (d9e9949)
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  fUseAliEventCuts(kTRUE),
66  fEventCuts(0),
67  fEventCutList(0),
68  fUseManualEventCuts(kFALSE),
69  fCellEnergyCut(0.05),
70  fMaxPt(250),
71  fNTotTracks(0),
72  fLeadingTrack(),
73  fDoTrackQA(kTRUE),
74  fDoCaloQA(kTRUE),
75  fDoJetQA(kTRUE),
76  fDoEventQA(kTRUE),
77  fGeneratorLevelName(),
78  fDetectorLevelName(),
79  fRejectOutlierEvents(kFALSE),
80  fIsPtHard(kFALSE),
81  fGeneratorLevel(0),
82  fDetectorLevel(0),
83  fNPtHistBins(0),
84  fPtHistBins(0),
85  fNEtaHistBins(0),
86  fEtaHistBins(0),
87  fNPhiHistBins(0),
88  fPhiHistBins(0),
89  fNCentHistBins(0),
90  fCentHistBins(0),
91  fNPtRelDiffHistBins(0),
92  fPtRelDiffHistBins(0),
93  fNPtResHistBins(0),
94  fPtResHistBins(0),
95  fNIntegerHistBins(0),
96  fIntegerHistBins(0),
97  fPHOSGeo(nullptr),
98  fHistManager("AliAnalysisTaskPWGJEQA")
99 {
100  // Default constructor.
101  memset(fNTotClusters, 0, sizeof(Int_t)*3);
104 }
105 
106 //________________________________________________________________________
111  AliAnalysisTaskEmcalJet(name, kTRUE),
112  fUseAliEventCuts(kTRUE),
113  fEventCuts(0),
114  fEventCutList(0),
115  fUseManualEventCuts(kFALSE),
116  fCellEnergyCut(0.05),
117  fMaxPt(250),
118  fNTotTracks(0),
119  fLeadingTrack(),
120  fDoTrackQA(kTRUE),
121  fDoCaloQA(kTRUE),
122  fDoJetQA(kTRUE),
123  fDoEventQA(kTRUE),
126  fRejectOutlierEvents(kFALSE),
127  fIsPtHard(kFALSE),
128  fGeneratorLevel(0),
129  fDetectorLevel(0),
130  fNPtHistBins(0),
131  fPtHistBins(0),
132  fNEtaHistBins(0),
133  fEtaHistBins(0),
134  fNPhiHistBins(0),
135  fPhiHistBins(0),
136  fNCentHistBins(0),
137  fCentHistBins(0),
140  fNPtResHistBins(0),
141  fPtResHistBins(0),
143  fIntegerHistBins(0),
144  fPHOSGeo(nullptr),
145  fHistManager(name)
146 {
147  // Standard
148  memset(fNTotClusters, 0, sizeof(Int_t)*3);
151 }
152 
153 //________________________________________________________________________
155 {
156  // Destructor
157 }
158 
159 //________________________________________________________________________
161 {
162  // Create histograms
164 
165  // Intialize AliEventCuts
166  if (fUseAliEventCuts) {
167  fEventCutList = new TList();
168  fEventCutList ->SetOwner();
169  fEventCutList ->SetName("EventCutOutput");
170 
171  fEventCuts.OverrideAutomaticTriggerSelection(fOffTrigger);
172  if(fUseManualEventCuts==1)
173  {
174  fEventCuts.SetManualMode();
175  fEventCuts.fMC = MCEvent(); //before was= false
176  if(fForceBeamType != kpp)
177  {
178  fEventCuts.SetupLHC15o();
179  fEventCuts.fUseVariablesCorrelationCuts = true;
180  }
181  else if(fForceBeamType == kpp)
182  {
183  fEventCuts.SetupRun2pp();
184  //no other cuts known so far
185  }
186  else
187  {
188  printf("No implementation of manuel event cuts for pPb yet!");
189  }
190  }
191  fEventCuts.AddQAplotsToList(fEventCutList);
192  fOutput->Add(fEventCutList);
193  }
194 
195  // Set track container pointers
198 
199  // Allocate histograms for tracks, cells, clusters, jets
205 
206  TIter next(fHistManager.GetListOfHistograms());
207  TObject* obj = 0;
208  while ((obj = next())) {
209  fOutput->Add(obj);
210  }
211 
212  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
213 }
214 
215 //________________________________________________________________________
217 
219  if (fGeneratorLevel) {
222  }
223 }
224 
225 //________________________________________________________________________
227 
228  TString histname;
229  TString title;
230 
231  if (!fCaloCellsName.IsNull()) {
232 
233  histname = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
234  title = histname + ";#it{E}_{cell} (GeV);counts";
235  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,0, 100);
236 
237  histname = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
238  title = histname + ";cell absId;<#it{E}_{cell}> (GeV)";
239  fHistManager.CreateTProfile(histname.Data(), title.Data(), 18000,0,18000);
240 
241  histname = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
242  title = histname + ";#it{t}_{cell} (s);counts";
243  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,-3e-6, 3e-6);
244 
245  histname = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
246  title = histname + ";cell absId;<#it{t}_{cell}> (s)";
247  fHistManager.CreateTProfile(histname.Data(), title.Data(), 18000,0,18000);
248 
249  histname = TString::Format("%s/fHistCellEvsTime", fCaloCellsName.Data());
250  title = histname + ";<#it{t}_{cell}> (s);<#it{E}_{cell}> (GeV)";
251  fHistManager.CreateTH2(histname.Data(), title.Data(), 1000, -1e-6, 1e-6, 200, 0, 100);
252 
253  }
254 }
255 
256 //________________________________________________________________________
258 
259  AliEmcalContainer* cont = 0;
260  TString histname;
261  TString title;
262  Int_t nPtBins = 2*fMaxPt;
263 
264  TIter nextClusColl(&fClusterCollArray);
265  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
266 
267  histname = TString::Format("%s/fHistClusterRejectionReason", cont->GetArrayName().Data());
268  title = histname + ";Rejection reason;#it{E}_{clus} (GeV/);counts";
269  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
270  SetRejectionReasonLabels(hist->GetXaxis());
271 
272  histname = TString::Format("%s/hPhosNCellsVsEnergy", cont->GetArrayName().Data());
273  title = histname + ";N cells;#it{E}_{clus} (GeV/);counts";
274  fHistManager.CreateTH2(histname.Data(), title.Data(), 40, 0, 40, nPtBins, 0, fMaxPt);
275 
276  histname = TString::Format("%s/hPhosM02VsEnergy", cont->GetArrayName().Data());
277  title = histname + ";M02;#it{E}_{clus} (GeV/);counts";
278  fHistManager.CreateTH2(histname.Data(), title.Data(), 100, 0, 20, nPtBins, 0, fMaxPt);
279 
280  histname = TString::Format("%s/hPhosCellIdVsEnergy", cont->GetArrayName().Data());
281  title = histname + ";Cell ID;#it{E}_{clus} (GeV/);counts";
282  fHistManager.CreateTH2(histname.Data(), title.Data(), 20000, 0, 20000, nPtBins, 0, fMaxPt);
283 
284  const Int_t nEmcalSM = 20;
285  for (Int_t sm = 0; sm < nEmcalSM; sm++) {
286  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
287  title = histname + ";#it{E}_{cluster} (GeV);counts";
288  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
289  }
290 
291  for (Int_t sm = 1; sm < 5; sm++) {
292  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
293  title = histname + ";#it{E}_{cluster} (GeV);counts";
294  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
295  }
296 
297  // Cluster THnSparse
298  Int_t dim = 0;
299  TString title[20];
300  Int_t nbins[20] = {0};
301  Double_t min[30] = {0.};
302  Double_t max[30] = {0.};
303  Double_t *binEdges[20] = {0};
304 
306  title[dim] = "Centrality %";
307  nbins[dim] = fNCentHistBins;
308  binEdges[dim] = fCentHistBins;
309  min[dim] = fCentHistBins[0];
310  max[dim] = fCentHistBins[fNCentHistBins];
311  dim++;
312  }
313 
314  title[dim] = "#it{E}_{clus} (GeV)";
315  nbins[dim] = fNPtHistBins;
316  binEdges[dim] = fPtHistBins;
317  min[dim] = fPtHistBins[0];
318  max[dim] = fPtHistBins[fNPtHistBins];
319  dim++;
320 
321  title[dim] = "#eta";
322  nbins[dim] = fNEtaHistBins;
323  binEdges[dim] = fEtaHistBins;
324  min[dim] = fEtaHistBins[0];
325  max[dim] = fEtaHistBins[fNEtaHistBins];
326  dim++;
327 
328  title[dim] = "#phi";
329  nbins[dim] = fNPhiHistBins;
330  binEdges[dim] = fPhiHistBins;
331  min[dim] = fPhiHistBins[0];
332  max[dim] = fPhiHistBins[fNPhiHistBins];
333  dim++;
334 
335  title[dim] = "cluster type";
336  nbins[dim] = 3;
337  binEdges[dim] = fIntegerHistBins;
338  min[dim] = fIntegerHistBins[0];
339  max[dim] = fIntegerHistBins[3];
340  dim++;
341 
342  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
343  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
344  for (Int_t i = 0; i < dim; i++) {
345  hn->GetAxis(i)->SetTitle(title[i]);
346  hn->SetBinEdges(i, binEdges[i]);
347  }
348  }
349 }
350 
351 //________________________________________________________________________
353 
354  AliJetContainer* jets = 0;
355  TIter nextJetColl(&fJetCollArray);
356  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
357 
358  // Allocate THnSparse
359 
360  TString axisTitle[30]= {""};
361  Int_t nbins[30] = {0};
362  Double_t min[30] = {0.};
363  Double_t max[30] = {0.};
364  Double_t *binEdges[20] = {0};
365  Int_t dim = 0;
366 
367  if (fForceBeamType != kpp) {
368  axisTitle[dim] = "Centrality (%)";
369  nbins[dim] = fNCentHistBins;
370  binEdges[dim] = fCentHistBins;
371  min[dim] = fCentHistBins[0];
372  max[dim] = fCentHistBins[fNCentHistBins];
373  dim++;
374  }
375 
376  axisTitle[dim] = "#eta_{jet}";
377  nbins[dim] = fNEtaHistBins;
378  binEdges[dim] = fEtaHistBins;
379  min[dim] = fEtaHistBins[0];
380  max[dim] = fEtaHistBins[fNEtaHistBins];
381  dim++;
382 
383  axisTitle[dim] = "#phi_{jet} (rad)";
384  nbins[dim] = fNPhiHistBins;
385  binEdges[dim] = fPhiHistBins;
386  min[dim] = fPhiHistBins[0];
387  max[dim] = fPhiHistBins[fNPhiHistBins];
388  dim++;
389 
390  axisTitle[dim] = "#it{p}_{T} (GeV/#it{c})";
391  nbins[dim] = TMath::CeilNint(fMaxPt/3);
392  min[dim] = 0;
393  max[dim] = fMaxPt;
394  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
395  dim++;
396 
397  if (fForceBeamType != kpp) {
398  axisTitle[dim] = "#it{p}_{T}^{corr} (GeV/#it{c})";
399  nbins[dim] = TMath::CeilNint(fMaxPt/3);
400  min[dim] = -fMaxPt/2 + 25;
401  max[dim] = fMaxPt/2 + 25;
402  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
403  dim++;
404  }
405 
406  axisTitle[dim] = "#it{p}_{T,particle}^{leading} (GeV/#it{c})";
407  nbins[dim] = fMaxPt/5;
408  min[dim] = 0;
409  max[dim] = 150;
410  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
411  dim++;
412 
413  TString thnname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
414  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
415  for (Int_t i = 0; i < dim; i++) {
416  hn->GetAxis(i)->SetTitle(axisTitle[i]);
417  hn->SetBinEdges(i, binEdges[i]);
418  }
419 
420  // Allocate other jet histograms
421  TString histname;
422  TString title;
423  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
424  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
425  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
426  SetRejectionReasonLabels(hist->GetXaxis());
427 
428  if (!jets->GetRhoName().IsNull()) {
429  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
430  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
431  fHistManager.CreateTH2(histname.Data(), title.Data(), 101, 0, 101, 100, 0, 500);
432  }
433  if (fForceBeamType != kpp) {
434  histname = TString::Format("%s/hNEFVsPtEMC", jets->GetArrayName().Data());
435  title = histname + ";Centrality (%); #it{p}_{T}^{corr} (GeV/#it{c});NEF";
436  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, 250, 0, 250, 50, 0, 1);
437 
438  histname = TString::Format("%s/hNEFVsPtDCal", jets->GetArrayName().Data());
439  title = histname + ";Centrality (%); #it{p}_{T}^{corr} (GeV/#it{c});NEF";
440  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, 250, 0, 250, 50, 0, 1);
441 
442  histname = TString::Format("%s/hDeltaEHadCorr", jets->GetArrayName().Data());
443  title = histname + ";Centrality (%); #it{p}_{T}^{corr} (GeV/#it{c}); #sum #it{E_{nonlincorr}} - #it{E}_{hadcorr}";
444  fHistManager.CreateTH3(histname.Data(), title.Data(), 20, 0, 100, 250, 0, 250, 250, 0, 50);
445  }
446  else{
447  histname = TString::Format("%s/hNEFVsPtEMC", jets->GetArrayName().Data());
448  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
449  fHistManager.CreateTH2(histname.Data(), title.Data(), 250, 0, 250, 50, 0, 1);
450 
451  histname = TString::Format("%s/hNEFVsPtDCal", jets->GetArrayName().Data());
452  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c});NEF";
453  fHistManager.CreateTH2(histname.Data(), title.Data(), 250, 0, 250, 50, 0, 1);
454 
455  histname = TString::Format("%s/hDeltaEHadCorr", jets->GetArrayName().Data());
456  title = histname + ";#it{p}_{T}^{corr} (GeV/#it{c}); #sum#it{E_{nonlincorr}} - #it{E}_{hadcorr}";
457  fHistManager.CreateTH2(histname.Data(), title.Data(), 250, 0, 250, 250, 0, 50);
458  }
459  }
460 }
461 
462 //________________________________________________________________________
464 
465  Int_t nPtBins = 2*fMaxPt;
466 
467  Int_t dim = 0;
468  TString axistitle[40];
469  Int_t nbins[40] = {0};
470  Double_t min[40] = {0};
471  Double_t max[40] = {0};
472  Double_t *binEdges[20] = {0};
473 
475  axistitle[dim] = "Centrality %";
476  nbins[dim] = fNCentHistBins;
477  binEdges[dim] = fCentHistBins;
478  min[dim] = fCentHistBins[0];
479  max[dim] = fCentHistBins[fNCentHistBins];
480  dim++;
481  }
482 
483  if (fParticleCollArray.GetEntriesFast()>0) {
484  axistitle[dim] = "No. of tracks";
486  nbins[dim] = 50;
487  min[dim] = -0.5;
488  max[dim] = 10000-0.5;
489  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
490  }
491  else {
492  nbins[dim] = 50;
493  min[dim] = 0;
494  max[dim] = 200;
495  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
496  }
497  dim++;
498 
499  axistitle[dim] = "#it{p}_{T,track}^{leading} (GeV/c)";
500  nbins[dim] = fMaxPt/5;
501  min[dim] = 0;
502  max[dim] = 150;
503  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
504  dim++;
505  }
506 
507  if (fClusterCollArray.GetEntriesFast()>0) {
508  axistitle[dim] = "No. of clusters";
509 
511  nbins[dim] = 50;
512  min[dim] = -0.5;
513  max[dim] = 4000-0.5;
514  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
515  }
516  else {
517  nbins[dim] = 50;
518  min[dim] = 0;
519  max[dim] = 200;
520  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
521  }
522  dim++;
523 
524  axistitle[dim] = "#it{E}_{cluster}^{leading} (GeV)";
525  nbins[dim] = fMaxPt/5;
526  min[dim] = 0;
527  max[dim] = 150;
528  binEdges[dim] = GenerateFixedBinArray(nbins[dim], min[dim], max[dim]);
529  dim++;
530  }
531 
532  THnSparse* hn = fHistManager.CreateTHnSparse("eventQA","eventQA",dim,nbins,min,max);
533  for (Int_t i = 0; i < dim; i++) {
534  hn->GetAxis(i)->SetTitle(axistitle[i]);
535  hn->SetBinEdges(i, binEdges[i]);
536  }
537 
538  if (fIsPtHard) {
539  TString histname = "hPtHard";
540  TString title = histname + ";#it{p}_{T,hard} (GeV/c);counts";
541  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
542 
543  TH1* hTrials = fHistManager.CreateTH1("hNtrials", "hNtrials", 1, 0, 1);
544  hTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
545 
546  TH1* hxsection = fHistManager.CreateTH1("hXsec", "hXsec", 1, 0, 1);
547  hxsection->GetXaxis()->SetBinLabel(1,"<#sigma>");
548  }
549 
550  fHistEventRejection->GetXaxis()->SetBinLabel(15,"PtHardBinJetOutlier");
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] = "track type";
595  nbins[dim] = 3;
596  binEdges[dim] = fIntegerHistBins;
597  min[dim] = fIntegerHistBins[0];
598  max[dim] = fIntegerHistBins[3];
599  dim++;
600 
601  title[dim] = "#sigma(#it{p}_{T}) / #it{p}_{T}";
602  nbins[dim] = fNPtResHistBins;
603  binEdges[dim] = fPtResHistBins;
604  min[dim] = fPtResHistBins[0];
605  max[dim] = fPtResHistBins[fNPtResHistBins];
606  dim++;
607 
608  THnSparse* hn = fHistManager.CreateTHnSparse("tracks","tracks", dim, nbins, min, max);
609  for (Int_t i = 0; i < dim; i++) {
610  hn->GetAxis(i)->SetTitle(title[i]);
611  hn->SetBinEdges(i, binEdges[i]);
612  }
613 
614 }
615 
616 //________________________________________________________________________
618 {
619  Int_t dim = 0;
620  TString title[20];
621  Int_t nbins[20] = {0};
622  Double_t min[30] = {0.};
623  Double_t max[30] = {0.};
624  Double_t *binEdges[20] = {0};
625 
627  title[dim] = "Centrality %";
628  nbins[dim] = fNCentHistBins;
629  binEdges[dim] = fCentHistBins;
630  min[dim] = fCentHistBins[0];
631  max[dim] = fCentHistBins[fNCentHistBins];
632  dim++;
633  }
634 
635  title[dim] = "#it{p}_{T} (GeV/#it{c})";
636  nbins[dim] = fNPtHistBins;
637  binEdges[dim] = fPtHistBins;
638  min[dim] = fPtHistBins[0];
639  max[dim] = fPtHistBins[fNPtHistBins];
640  dim++;
641 
642  title[dim] = "#eta";
643  nbins[dim] = fNEtaHistBins;
644  binEdges[dim] = fEtaHistBins;
645  min[dim] = fEtaHistBins[0];
646  max[dim] = fEtaHistBins[fNEtaHistBins];
647  dim++;
648 
649  title[dim] = "#phi";
650  nbins[dim] = fNPhiHistBins;
651  binEdges[dim] = fPhiHistBins;
652  min[dim] = fPhiHistBins[0];
653  max[dim] = fPhiHistBins[fNPhiHistBins];
654  dim++;
655 
656  title[dim] = "Findable";
657  nbins[dim] = 2;
658  binEdges[dim] = fIntegerHistBins;
659  min[dim] = fIntegerHistBins[0];
660  max[dim] = fIntegerHistBins[2];
661  dim++;
662 
663  THnSparse* hn = fHistManager.CreateTHnSparse("tracks_PhysPrim","tracks_PhysPrim", dim, nbins, min, max);
664  for (Int_t i = 0; i < dim; i++) {
665  hn->GetAxis(i)->SetTitle(title[i]);
666  hn->SetBinEdges(i, binEdges[i]);
667  }
668 }
669 
670 //________________________________________________________________________
672 {
673  Int_t dim = 0;
674  TString title[20];
675  Int_t nbins[20] = {0};
676  Double_t min[30] = {0.};
677  Double_t max[30] = {0.};
678  Double_t *binEdges[20] = {0};
679 
681  title[dim] = "Centrality %";
682  nbins[dim] = fNCentHistBins;
683  binEdges[dim] = fCentHistBins;
684  min[dim] = fCentHistBins[0];
685  max[dim] = fCentHistBins[fNCentHistBins];
686  dim++;
687  }
688 
689  title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
690  nbins[dim] = fNPtHistBins;
691  binEdges[dim] = fPtHistBins;
692  min[dim] = fPtHistBins[0];
693  max[dim] = fPtHistBins[fNPtHistBins];
694  dim++;
695 
696  title[dim] = "#eta^{gen}";
697  nbins[dim] = fNEtaHistBins;
698  binEdges[dim] = fEtaHistBins;
699  min[dim] = fEtaHistBins[0];
700  max[dim] = fEtaHistBins[fNEtaHistBins];
701  dim++;
702 
703  title[dim] = "#phi^{gen}";
704  nbins[dim] = fNPhiHistBins;
705  binEdges[dim] = fPhiHistBins;
706  min[dim] = fPhiHistBins[0];
707  max[dim] = fPhiHistBins[fNPhiHistBins];
708  dim++;
709 
710  title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
711  nbins[dim] = fNPtHistBins;
712  binEdges[dim] = fPtHistBins;
713  min[dim] = fPtHistBins[0];
714  max[dim] = fPtHistBins[fNPtHistBins];
715  dim++;
716 
717  title[dim] = "#eta^{det}";
718  nbins[dim] = fNEtaHistBins;
719  binEdges[dim] = fEtaHistBins;
720  min[dim] = fEtaHistBins[0];
721  max[dim] = fEtaHistBins[fNEtaHistBins];
722  dim++;
723 
724  title[dim] = "#phi^{det}";
725  nbins[dim] = fNPhiHistBins;
726  binEdges[dim] = fPhiHistBins;
727  min[dim] = fPhiHistBins[0];
728  max[dim] = fPhiHistBins[fNPhiHistBins];
729  dim++;
730 
731  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
732  nbins[dim] = fNPtRelDiffHistBins;
733  binEdges[dim] = fPtRelDiffHistBins;
734  min[dim] = fPtRelDiffHistBins[0];
736  dim++;
737 
738  title[dim] = "track type";
739  nbins[dim] = 3;
740  binEdges[dim] = fIntegerHistBins;
741  min[dim] = fIntegerHistBins[0];
742  max[dim] = fIntegerHistBins[3];
743  dim++;
744 
745  THnSparse* hn = fHistManager.CreateTHnSparse("tracks_Matched","tracks_Matched", dim, nbins, min, max);
746  for (Int_t i = 0; i < dim; i++) {
747  hn->GetAxis(i)->SetTitle(title[i]);
748  hn->SetBinEdges(i, binEdges[i]);
749  }
750 }
751 
752 
753 //________________________________________________________________________
756 {
757  fNPtHistBins = 82;
760  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
761  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
762  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
763  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
764  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
765  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
766 
767  fNEtaHistBins = 70;
770 
771  fNPhiHistBins = 101;
773  GenerateFixedBinArray(fNPhiHistBins, 0, TMath::Pi() * 2.02, fPhiHistBins);
774 
775  fNCentHistBins = 4;
777  fCentHistBins[0] = 0;
778  fCentHistBins[1] = 10;
779  fCentHistBins[2] = 30;
780  fCentHistBins[3] = 50;
781  fCentHistBins[4] = 90;
782 
783  fNPtResHistBins = 100;
786  GenerateFixedBinArray(25, 0.05, 0.10, fPtResHistBins+50);
787  GenerateFixedBinArray(15, 0.10, 0.20, fPtResHistBins+75);
788  GenerateFixedBinArray(10, 0.20, 0.50, fPtResHistBins+90);
789 
790  fNPtRelDiffHistBins = 50;
793 
794  fNIntegerHistBins = 3;
797 }
798 
799 //________________________________________________________________________
801 {
802  if (fClusterCollArray.GetEntriesFast() == 0 && fCaloCellsName.IsNull()) {
803  fNeedEmcalGeom = kFALSE;
804  }
805  else {
806  fNeedEmcalGeom = kTRUE;
807  }
808 
810 
811  // Load the PHOS geometry
812  fPHOSGeo = AliPHOSGeometry::GetInstance();
813  if (fPHOSGeo) {
814  AliInfo("Found instance of PHOS geometry!");
815  }
816  else {
817  AliInfo("Creating PHOS geometry!");
818  Int_t runNum = InputEvent()->GetRunNumber();
819  if(runNum<209122) //Run1
820  fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP");
821  else
822  fPHOSGeo = AliPHOSGeometry::GetInstance("Run2");
823 
824  if (fPHOSGeo) {
825  AliOADBContainer geomContainer("phosGeo");
826  geomContainer.InitFromFile("$ALICE_PHYSICS/OADB/PHOS/PHOSMCGeometry.root","PHOSMCRotationMatrixes");
827  TObjArray* matrixes = (TObjArray*)geomContainer.GetObject(runNum,"PHOSRotationMatrixes");
828  for(Int_t mod=0; mod<6; mod++) {
829  if(!matrixes->At(mod)) continue;
830  fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
831  printf(".........Adding Matrix(%d), geo=%p\n",mod,fPHOSGeo);
832  ((TGeoHMatrix*)matrixes->At(mod))->Print();
833  }
834  }
835  }
836 }
837 
842 {
843  if (fUseAliEventCuts) {
844  if (!fEventCuts.AcceptEvent(InputEvent()))
845  {
846  PostData(1, fOutput);
847  return kFALSE;
848  }
849  }
850  else {
852  }
853  return kTRUE;
854 }
855 //________________________________________________________________________
857 {
858  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
859 
860  // If Pt-hard production, get the Pt-hard of the event, and have possibility to reject the event for jet outliers
861  if (fIsPtHard) {
862  AliGenPythiaEventHeader* pygen = nullptr;
863  if (MCEvent()) {
864  pygen = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
865  if (!pygen) {
866  // Check if AOD
867  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
868  if (aodMCH) {
869  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
870  pygen = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
871  if (pygen) break;
872  }
873  }
874  }
875  }
876  if (pygen) {
877  fPtHard = pygen->GetPtHard();
878  //fNTrials = pygen->Trials();
879  //fXsection = pygen->GetXsection();
880 
881  // reject outlier events, where the jet Pt is much larger than Pt-hard
882  if (fRejectOutlierEvents) {
883  AliTLorentzVector jet;
884  Int_t nTriggerJets = pygen->NTriggerJets();
885  Float_t tmpjet[]={0,0,0,0};
886  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
887  pygen->TriggerJet(ijet, tmpjet);
888  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
889 
890  if (jet.Pt() > 4. * fPtHard) {
891  //AliInfo(Form("Reject jet event with: pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), 4.));
892  if (fGeneralHistograms) fHistEventRejection->Fill("PtHardBinJetOutlier",1);
893  return kFALSE;
894  }
895  }
896  }
897  }
898  }
899 
900  return kTRUE;
901 }
902 
910 //________________________________________________________________________
912 {
913  if (fIsPtHard) {
914  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
915  if (!tree) {
916  AliError(Form("%s - UserNotify: No current tree!",GetName()));
917  return kFALSE;
918  }
919 
920  Float_t xsection = 0;
921  Float_t trials = 0;
922  Int_t pthardbin = 0;
923 
924  TFile *curfile = tree->GetCurrentFile();
925  if (!curfile) {
926  AliError(Form("%s - UserNotify: No current file!",GetName()));
927  return kFALSE;
928  }
929 
930  TChain *chain = dynamic_cast<TChain*>(tree);
931  if (chain) tree = chain->GetTree();
932 
933  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
934 
935  fHistManager.FillTH1("hNtrials", "#sum{ntrials}", trials);
936  fHistManager.FillTH1("hXsec", "<#sigma>", xsection);
937  }
938  return kTRUE;
939 }
940 
941 //________________________________________________________________________
943 {
949 
950  return kTRUE;
951 }
952 
953 //________________________________________________________________________
955 
956  fNTotTracks = 0;
957  fLeadingTrack.SetPxPyPzE(0,0,0,0);
958 
959  AliVTrack* track;
960  for (auto trackIterator : fDetectorLevel->accepted_momentum() ) {
961 
962  fNTotTracks++;
963  if (fLeadingTrack.Pt() < trackIterator.first.Pt()) fLeadingTrack = trackIterator.first;
964 
965  track = trackIterator.second;
966  Byte_t type = fDetectorLevel->GetTrackType(track);
967  if (type <= 2) {
968  Double_t sigma = 0;
969 
970  if (fIsEsd) {
971 
972  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
973  if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
974 
975  } else { // AOD
976 
977  AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
978  if(!aodtrack) AliFatal("Not a standard AOD");
979 
980  AliExternalTrackParam exParam;
981 
982  //get covariance matrix
983  Double_t cov[21] = {0,};
984  aodtrack->GetCovMatrix(cov);
985  Double_t pxpypz[3] = {0,};
986  aodtrack->PxPyPz(pxpypz);
987  Double_t xyz[3] = {0,};
988  aodtrack->GetXYZ(xyz);
989  Short_t sign = aodtrack->Charge();
990  exParam.Set(xyz,pxpypz,cov,sign);
991 
992  sigma = TMath::Sqrt(exParam.GetSigma1Pt2());
993 
994  }
995 
996  Int_t label = TMath::Abs(track->GetLabel());
997 
998  FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, type);
999 
1000  if (fGeneratorLevel && label > 0) {
1001  AliAODMCParticle *part = fGeneratorLevel->GetAcceptMCParticleWithLabel(label);
1002  if (part) {
1003  if (part->GetGeneratorIndex() == 0) {
1004  Int_t pdg = TMath::Abs(part->PdgCode());
1005  // select charged pions, protons, kaons, electrons, muons
1006  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
1007  FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
1008  }
1009  }
1010  }
1011  }
1012  }
1013  else {
1014  AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
1015  }
1016  }
1017 
1018  if (fGeneratorLevel) {
1019  AliAODMCParticle* part;
1020  for (auto partIterator : fGeneratorLevel->accepted_momentum() ) {
1021  part = partIterator.second;
1022 
1023  Byte_t findable = 0;
1024 
1025  Int_t pdg = TMath::Abs(part->PdgCode());
1026  // select charged pions, protons, kaons, electrons, muons
1027  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
1028 
1029  FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), findable);
1030  }
1031  }
1032 }
1033 
1034 //________________________________________________________________________
1036 
1037  if (!fCaloCells) return;
1038 
1039  TString histname_energy = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
1040  TString histname_time = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
1041  TString histname_energyProf = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
1042  TString histname_timeProf = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
1043  TString histname_EnergyvsTime = TString::Format("%s/fHistCellEvsTime", fCaloCellsName.Data());
1044 
1045  const Int_t ncells = fCaloCells->GetNumberOfCells();
1046  for (Int_t pos = 0; pos < ncells; pos++) {
1047  Float_t amp = fCaloCells->GetAmplitude(pos);
1048  Float_t time = fCaloCells->GetTime(pos);
1049  Int_t absId = fCaloCells->GetCellNumber(pos);
1050 
1051  if (amp < fCellEnergyCut) continue;
1052 
1053  fHistManager.FillTH1(histname_energy, amp);
1054  fHistManager.FillTH1(histname_time, time);
1055 
1056  fHistManager.FillProfile(histname_energyProf, absId, amp);
1057  fHistManager.FillProfile(histname_timeProf, absId, time);
1058 
1059  fHistManager.FillTH2(histname_EnergyvsTime, time, amp);
1060  }
1061 }
1062 
1063 //________________________________________________________________________
1065 
1066  memset(fNTotClusters, 0, sizeof(Int_t)*3);
1067  for (Int_t i = 0; i < 3; i++) fLeadingCluster[i].SetPxPyPzE(0,0,0,0);
1068 
1069  TString histname;
1070  AliClusterContainer* clusters = 0;
1071  TIter nextClusColl(&fClusterCollArray);
1072  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
1073  // Cluster loop
1075  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1076 
1077  UInt_t rejectionReason = 0;
1078  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
1079  histname = TString::Format("%s/fHistClusterRejectionReason", clusters->GetArrayName().Data());
1080  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
1081  continue;
1082  }
1083 
1084  // Determine cluster type (EMCal/DCal/Phos) and record relevant eventQA info
1085  ClusterType clusType = kNA;
1086  if (it->second->IsEMCAL()) {
1087  Double_t phi = it->first.Phi_0_2pi();
1088  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
1089  if (isDcal == 0) {
1090  clusType = kEMCal;
1091  } else if (isDcal == 1) {
1092  clusType = kDCal;
1093  }
1094 
1095  if (fLeadingCluster[isDcal].E() < it->first.E()) fLeadingCluster[isDcal] = it->first;
1096  fNTotClusters[isDcal]++;
1097 
1098  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
1099  if (sm >=0 && sm < 20) {
1100  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1101  fHistManager.FillTH1(histname, it->second->E());
1102  }
1103  else {
1104  AliError(Form("Supermodule %d does not exist!", sm));
1105  }
1106 
1107  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
1108 
1109  Int_t nCells = it->second->GetNCells();
1110  Double_t M02 = it->second->GetM02();
1111  Double_t energy = it->second->E();
1112 
1113  histname = TString::Format("%s/hPhosNCellsVsEnergy", clusters->GetArrayName().Data());
1114  fHistManager.FillTH2(histname, nCells, energy);
1115 
1116  histname = TString::Format("%s/hPhosM02VsEnergy", clusters->GetArrayName().Data());
1117  fHistManager.FillTH2(histname, M02, energy);
1118 
1119  clusType = kPHOS;
1120  fNTotClusters[2]++;
1121  if (fLeadingCluster[2].E() < it->first.E()) fLeadingCluster[2] = it->first;
1122 
1123  // Fill Phos spectra by module
1124  Int_t relid[4];
1125  if (fPHOSGeo) {
1126  fPHOSGeo->AbsToRelNumbering(it->second->GetCellAbsId(0), relid);
1127  Int_t sm = relid[0];
1128  if (sm >=1 && sm < 5) {
1129  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1130  fHistManager.FillTH1(histname, energy);
1131  }
1132  else {
1133  AliError(Form("Supermodule %d does not exist!", sm));
1134  }
1135  }
1136 
1137  // Loop through cells in the cluster, and histogram their energy
1138  histname = TString::Format("%s/hPhosCellIdVsEnergy", clusters->GetArrayName().Data());
1139  for (Int_t i=0; i < nCells; i++) {
1140  fHistManager.FillTH2(histname, it->second->GetCellAbsId(i), energy);
1141  }
1142 
1143  }
1144 
1145  Double_t contents[30]={0};
1146  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
1147  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1148  if (!histClusterObservables) return;
1149  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1150  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1151  if (title=="Centrality %")
1152  contents[i] = fCent;
1153  else if (title=="#it{E}_{clus} (GeV)")
1154  contents[i] = it->first.E();
1155  else if (title=="#eta")
1156  contents[i] = it->first.Eta();
1157  else if (title=="#phi")
1158  contents[i] = it->first.Phi_0_2pi();
1159  else if (title=="cluster type")
1160  contents[i] = clusType;
1161  else
1162  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1163  }
1164  histClusterObservables->Fill(contents);
1165 
1166  }
1167  }
1168 }
1169 
1170 //________________________________________________________________________
1172 
1173  TString histname;
1174  AliJetContainer* jets = 0;
1175  TIter nextJetColl(&fJetCollArray);
1176  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1177  Double_t rhoVal = 0;
1178  if (jets->GetRhoParameter()) {
1179  rhoVal = jets->GetRhoVal();
1180  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
1181  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1182  }
1183 
1184  for (auto jet : jets->all()) {
1185 
1186  UInt_t rejectionReason = 0;
1187  if (!jets->AcceptJet(jet, rejectionReason)) {
1188  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
1189  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1190  continue;
1191  }
1192 
1193  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1194  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
1195 
1196  TLorentzVector leadPart;
1197  jets->GetLeadingHadronMomentum(leadPart, jet);
1198 
1199  Double_t contents[30]={0};
1200  histname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
1201  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1202  if (!histJetObservables) return;
1203  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
1204  TString title(histJetObservables->GetAxis(i)->GetTitle());
1205  if (title=="Centrality (%)")
1206  contents[i] = fCent;
1207  else if (title=="#eta_{jet}")
1208  contents[i] = jet->Eta();
1209  else if (title=="#phi_{jet} (rad)")
1210  contents[i] = jet->Phi_0_2pi();
1211  else if (title=="#it{p}_{T} (GeV/#it{c})")
1212  contents[i] = jet->Pt();
1213  else if (title=="#it{p}_{T}^{MC} (GeV/#it{c})")
1214  contents[i] = jet->MCPt();
1215  else if (title=="#it{p}_{T}^{corr} (GeV/#it{c})")
1216  contents[i] = corrPt;
1217  else if (title=="#it{p}_{T,particle}^{leading} (GeV/#it{c})")
1218  contents[i] = ptLeading;
1219  else
1220  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1221  }
1222  histJetObservables->Fill(contents);
1223 
1224  UInt_t jetType = jet->GetJetAcceptanceType();
1225  if(jetType & AliEmcalJet::kEMCALfid)
1226  {
1227  histname = TString::Format("%s/hNEFVsPtEMC", jets->GetArrayName().Data());
1228  if (fForceBeamType != kpp){fHistManager.FillTH3(histname.Data(), fCent, jet->Pt(), jet->NEF());}
1229  else {fHistManager.FillTH2(histname.Data(), jet->Pt(), jet->NEF());}
1230  }
1231  else if(jetType & AliEmcalJet::kDCALonlyfid)
1232  {
1233  histname = TString::Format("%s/hNEFVsPtDCal", jets->GetArrayName().Data());
1234  if (fForceBeamType != kpp){
1235  fHistManager.FillTH3(histname.Data(), fCent, jet->Pt(), jet->NEF());
1236  }
1237  else{
1238  fHistManager.FillTH2(histname.Data(), jet->Pt(), jet->NEF());
1239  }
1240  }
1241  Double_t deltaEhadcorr = 0;
1242  const AliVCluster* clus = nullptr;
1243  Int_t nClusters = jet->GetNumberOfClusters();
1244  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1245  clus = jet->Cluster(iClus);
1246  deltaEhadcorr += (clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy());
1247  }
1248 
1249  histname = TString::Format("%s/hDeltaEHadCorr", jets->GetArrayName().Data());
1250  if (fForceBeamType != kpp){
1251  fHistManager.FillTH3(histname, fCent, jet->Pt(), deltaEhadcorr);
1252  }
1253  else{
1254  fHistManager.FillTH2(histname, jet->Pt(), deltaEhadcorr);
1255  }
1256  } //jet loop
1257  }
1258 }
1259 
1260 //________________________________________________________________________
1262 
1263  EventQA_t eventQA;
1264  for (Int_t i = 0; i < 3; i++) {
1265  eventQA.fMaxCluster[i] = fLeadingCluster[i];
1266  }
1267  eventQA.fCent = fCent;
1268  eventQA.fNTracks = fNTotTracks;
1269  eventQA.fNClusters[0] = fNTotClusters[0];
1270  eventQA.fNClusters[1] = fNTotClusters[1];
1271  eventQA.fNClusters[2] = fNTotClusters[2];
1272  eventQA.fMaxTrack = fLeadingTrack;
1273 
1274  Int_t globalNclusters = eventQA.fNClusters[0] + eventQA.fNClusters[1] + eventQA.fNClusters[2];
1275  AliTLorentzVector globalMaxCluster;
1276  for (Int_t i = 0; i < 3; i++) {
1277  if (globalMaxCluster.E() < eventQA.fMaxCluster[i].E()) globalMaxCluster = eventQA.fMaxCluster[i];
1278  }
1279 
1280  THnSparse* histEventQA = static_cast<THnSparse*>(fHistManager.FindObject("eventQA"));
1281  Double_t contents[40]={0};
1282  for (Int_t i = 0; i < histEventQA->GetNdimensions(); i++) {
1283  TString title(histEventQA->GetAxis(i)->GetTitle());
1284  if (title=="Centrality %")
1285  contents[i] = eventQA.fCent;
1286  else if (title=="No. of tracks")
1287  contents[i] = eventQA.fNTracks;
1288  else if (title=="No. of clusters")
1289  contents[i] = globalNclusters;
1290  else if (title=="#it{p}_{T,track}^{leading} (GeV/c)")
1291  contents[i] = eventQA.fMaxTrack.Pt();
1292  else if (title=="#it{E}_{cluster}^{leading} (GeV)")
1293  contents[i] = globalMaxCluster.E();
1294  else
1295  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1296  }
1297 
1298  histEventQA->Fill(contents);
1299 
1300  // Fill pythia pt hard histograms, if applicable
1301  if (fPtHard > 1e-6) {
1302  fHistManager.FillTH1("hPtHard", fPtHard);
1303  }
1304 }
1305 
1306 //________________________________________________________________________
1308  Double_t sigma1OverPt, Byte_t trackType)
1309 {
1310  Double_t contents[20]={0};
1311 
1312  THnSparse* thnTracks = static_cast<THnSparse*>(fHistManager.FindObject("tracks"));
1313  if (!thnTracks) return;
1314  for (Int_t i = 0; i < thnTracks->GetNdimensions(); i++) {
1315  TString title(thnTracks->GetAxis(i)->GetTitle());
1316  if (title=="Centrality %")
1317  contents[i] = cent;
1318  else if (title=="#it{p}_{T} (GeV/#it{c})")
1319  contents[i] = trackPt;
1320  else if (title=="#eta")
1321  contents[i] = trackEta;
1322  else if (title=="#phi")
1323  contents[i] = trackPhi;
1324  else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
1325  contents[i] = sigma1OverPt;
1326  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
1327  contents[i] = sigma1OverPt*trackPt;
1328  else if (title=="track type")
1329  contents[i] = trackType;
1330  else
1331  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks->GetName()));
1332  }
1333 
1334  thnTracks->Fill(contents);
1335 }
1336 
1337 //________________________________________________________________________
1338 void AliAnalysisTaskPWGJEQA::FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Byte_t findable)
1339 {
1340  Double_t contents[20]={0};
1341 
1342  THnSparse* thnTracks_PhysPrim = static_cast<THnSparse*>(fHistManager.FindObject("tracks_PhysPrim"));
1343  if (!thnTracks_PhysPrim) return;
1344  for (Int_t i = 0; i < thnTracks_PhysPrim->GetNdimensions(); i++) {
1345  TString title(thnTracks_PhysPrim->GetAxis(i)->GetTitle());
1346  if (title=="Centrality %")
1347  contents[i] = cent;
1348  else if (title=="#it{p}_{T} (GeV/#it{c})")
1349  contents[i] = partPt;
1350  else if (title=="#eta")
1351  contents[i] = partEta;
1352  else if (title=="#phi")
1353  contents[i] = partPhi;
1354  else if (title=="Findable")
1355  contents[i] = findable;
1356  else
1357  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_PhysPrim->GetName()));
1358  }
1359 
1360  thnTracks_PhysPrim->Fill(contents);
1361 }
1362 
1363 //________________________________________________________________________
1365  Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
1366 {
1367  Double_t contents[20]={0};
1368 
1369  THnSparse* thnTracks_Matched = static_cast<THnSparse*>(fHistManager.FindObject("tracks_Matched"));
1370  if (!thnTracks_Matched) return;
1371  for (Int_t i = 0; i < thnTracks_Matched->GetNdimensions(); i++) {
1372  TString title(thnTracks_Matched->GetAxis(i)->GetTitle());
1373  if (title=="Centrality %")
1374  contents[i] = cent;
1375  else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
1376  contents[i] = partPt;
1377  else if (title=="#eta^{gen}")
1378  contents[i] = partEta;
1379  else if (title=="#phi^{gen}")
1380  contents[i] = partPhi;
1381  else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
1382  contents[i] = trackPt;
1383  else if (title=="#eta^{det}")
1384  contents[i] = trackEta;
1385  else if (title=="#phi^{det}")
1386  contents[i] = trackPhi;
1387  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
1388  contents[i] = (partPt - trackPt) / partPt;
1389  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}")
1390  contents[i] = (partPt - trackPt) / trackPt;
1391  else if (title=="track type")
1392  contents[i] = (Double_t)trackType;
1393  else
1394  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_Matched->GetName()));
1395  }
1396 
1397  thnTracks_Matched->Fill(contents);
1398 }
1399 
1406  const char* ntracks,
1407  const char* nclusters,
1408  const char* ncells,
1409  const char *nGenLev,
1410  Bool_t doTrackQA,
1411  Bool_t doCaloQA,
1412  Bool_t doJetQA,
1413  Bool_t doEventQA,
1414  Double_t trackPtCut,
1415  Double_t clusECut,
1416  const char* suffix)
1417 {
1418  // Get the pointer to the existing analysis manager via the static access method
1419  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1420  if (!mgr) {
1421  ::Error("AddTaskPWGJEQA", "No analysis manager to connect to.");
1422  return NULL;
1423  }
1424 
1425  // Check the analysis type using the event handlers connected to the analysis manager
1426  AliVEventHandler* handler = mgr->GetInputEventHandler();
1427  if (!handler) {
1428  ::Error("AddTaskPWGJEQA", "This task requires an input event handler");
1429  return NULL;
1430  }
1431 
1432  enum EDataType_t {
1433  kUnknown,
1434  kESD,
1435  kAOD
1436  };
1437 
1438  EDataType_t dataType = kUnknown;
1439 
1440  if (handler->InheritsFrom("AliESDInputHandler")) {
1441  dataType = kESD;
1442  }
1443  else if (handler->InheritsFrom("AliAODInputHandler")) {
1444  dataType = kAOD;
1445  }
1446 
1447  // Init the task and do settings
1448  TString trackName(ntracks);
1449  TString clusName(nclusters);
1450  TString cellName(ncells);
1451 
1452  if (trackName == "usedefault") {
1453  if (dataType == kESD) {
1454  trackName = "Tracks";
1455  }
1456  else if (dataType == kAOD) {
1457  trackName = "tracks";
1458  }
1459  else {
1460  trackName = "";
1461  }
1462  }
1463 
1464  if (clusName == "usedefault") {
1465  if (dataType == kESD) {
1466  clusName = "CaloClusters";
1467  }
1468  else if (dataType == kAOD) {
1469  clusName = "caloClusters";
1470  }
1471  else {
1472  clusName = "";
1473  }
1474  }
1475 
1476  if (cellName == "usedefault") {
1477  if (dataType == kESD) {
1478  cellName = "EMCALCells";
1479  }
1480  else if (dataType == kAOD) {
1481  cellName = "emcalCells";
1482  }
1483  else {
1484  cellName = "";
1485  }
1486  }
1487 
1488  TString name("AliAnalysisTaskPWGJEQA");
1489  if (!trackName.IsNull()) {
1490  name += "_";
1491  name += trackName;
1492  }
1493  if (!clusName.IsNull()) {
1494  name += "_";
1495  name += clusName;
1496  }
1497 
1498  if (!cellName.IsNull()) {
1499  name += "_";
1500  name += cellName;
1501  }
1502 
1503  if (strcmp(suffix,"")) {
1504  name += "_";
1505  name += suffix;
1506  }
1507 
1508  if (nGenLev && strcmp(nGenLev,"")!=0) cout<<"MC Part: "<< nGenLev<<endl;
1509  else cout<<"No MC Part: "<< nGenLev<<endl;
1510 
1511  AliAnalysisTaskPWGJEQA* qaTask = new AliAnalysisTaskPWGJEQA(name);
1512  qaTask->SetVzRange(-10,10);
1513  qaTask->SetNeedEmcalGeom(kFALSE);
1514  qaTask->SetCaloCellsName(cellName);
1515  qaTask->SetDetectorLevelName(trackName);
1516  if (nGenLev && strcmp(nGenLev,"")!=0) qaTask->SetGeneratorLevelName(nGenLev);
1517  qaTask->SetDoTrackQA(doTrackQA);
1518  qaTask->SetDoCaloQA(doCaloQA);
1519  qaTask->SetDoJetQA(doJetQA);
1520  qaTask->SetDoEventQA(doEventQA);
1521 
1522  // Add the detector-level track container
1523  if (trackName == "tracks" || trackName == "Tracks") {
1524  AliTrackContainer* trackCont = qaTask->AddTrackContainer(trackName);
1525  trackCont->SetFilterHybridTracks(kTRUE);
1526  }
1527  else if (!trackName.IsNull()) {
1528  qaTask->AddParticleContainer(trackName);
1529  }
1530  AliParticleContainer *partCont = qaTask->GetParticleContainer(trackName);
1531  if (partCont) {
1532  partCont->SetParticlePtCut(trackPtCut);
1533  }
1534 
1535  // Add the generator-level container, if specified
1536  if (nGenLev && strcmp(nGenLev,"")!=0) {
1537  AliMCParticleContainer* mcpartCont = qaTask->AddMCParticleContainer(nGenLev);
1538  mcpartCont->SelectPhysicalPrimaries(kTRUE);
1539  mcpartCont->SetParticlePtCut(0);
1540  }
1541 
1542  // Add the cluster container
1543  AliClusterContainer *clusterCont = qaTask->AddClusterContainer(clusName);
1544  if (clusterCont) {
1545  clusterCont->SetClusECut(0.);
1546  clusterCont->SetClusPtCut(0.);
1547  clusterCont->SetClusHadCorrEnergyCut(clusECut);
1548  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1549  }
1550 
1551  // Final settings, pass to manager and set the containers
1552  mgr->AddTask(qaTask);
1553 
1554  // Create containers for input/output
1555  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
1556 
1557  TString contName = TString::Format("%s_histos", name.Data());
1558  TString commonoutput;
1559  commonoutput = mgr->GetCommonFileName();
1560  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(),
1561  TList::Class(),AliAnalysisManager::kOutputContainer,
1562  commonoutput);
1563  mgr->ConnectInput (qaTask, 0, cinput1 );
1564  mgr->ConnectOutput (qaTask, 1, coutput1 );
1565 
1566  return qaTask;
1567 
1568 }
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
UInt_t fOffTrigger
offline trigger for event selection
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.
TList * fEventCutList
! Output list for event cut histograms
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.
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
AliEventCuts fEventCuts
event selection utility
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.
Bool_t fUseAliEventCuts
Flag to use AliEventCuts (otherwise AliAnalysisTaskEmcal will be used)
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)
virtual Bool_t IsEventSelected()
Performing event selection.
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 fUseManualEventCuts
Flag to use manual event cuts.
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