AliPhysics  master (3d17d9d)
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  return answer;
853  }
854  return kTRUE;
855 }
856 //________________________________________________________________________
858 {
859  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
860 
861  // If Pt-hard production, get the Pt-hard of the event, and have possibility to reject the event for jet outliers
862  if (fIsPtHard) {
863  AliGenPythiaEventHeader* pygen = nullptr;
864  if (MCEvent()) {
865  pygen = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
866  if (!pygen) {
867  // Check if AOD
868  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
869  if (aodMCH) {
870  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
871  pygen = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
872  if (pygen) break;
873  }
874  }
875  }
876  }
877  if (pygen) {
878  fPtHard = pygen->GetPtHard();
879  //fNTrials = pygen->Trials();
880  //fXsection = pygen->GetXsection();
881 
882  // reject outlier events, where the jet Pt is much larger than Pt-hard
883  if (fRejectOutlierEvents) {
884  AliTLorentzVector jet;
885  Int_t nTriggerJets = pygen->NTriggerJets();
886  Float_t tmpjet[]={0,0,0,0};
887  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
888  pygen->TriggerJet(ijet, tmpjet);
889  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
890 
891  if (jet.Pt() > 4. * fPtHard) {
892  //AliInfo(Form("Reject jet event with: pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), 4.));
893  if (fGeneralHistograms) fHistEventRejection->Fill("PtHardBinJetOutlier",1);
894  return kFALSE;
895  }
896  }
897  }
898  }
899  }
900 
901  return kTRUE;
902 }
903 
911 //________________________________________________________________________
913 {
914  if (fIsPtHard) {
915  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
916  if (!tree) {
917  AliError(Form("%s - UserNotify: No current tree!",GetName()));
918  return kFALSE;
919  }
920 
921  Float_t xsection = 0;
922  Float_t trials = 0;
923  Int_t pthardbin = 0;
924 
925  TFile *curfile = tree->GetCurrentFile();
926  if (!curfile) {
927  AliError(Form("%s - UserNotify: No current file!",GetName()));
928  return kFALSE;
929  }
930 
931  TChain *chain = dynamic_cast<TChain*>(tree);
932  if (chain) tree = chain->GetTree();
933 
934  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
935 
936  fHistManager.FillTH1("hNtrials", "#sum{ntrials}", trials);
937  fHistManager.FillTH1("hXsec", "<#sigma>", xsection);
938  }
939  return kTRUE;
940 }
941 
942 //________________________________________________________________________
944 {
950 
951  return kTRUE;
952 }
953 
954 //________________________________________________________________________
956 
957  fNTotTracks = 0;
958  fLeadingTrack.SetPxPyPzE(0,0,0,0);
959 
960  AliVTrack* track;
961  for (auto trackIterator : fDetectorLevel->accepted_momentum() ) {
962 
963  fNTotTracks++;
964  if (fLeadingTrack.Pt() < trackIterator.first.Pt()) fLeadingTrack = trackIterator.first;
965 
966  track = trackIterator.second;
967  Byte_t type = fDetectorLevel->GetTrackType(track);
968  if (type <= 2) {
969  Double_t sigma = 0;
970 
971  if (fIsEsd) {
972 
973  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
974  if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
975 
976  } else { // AOD
977 
978  AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
979  if(!aodtrack) AliFatal("Not a standard AOD");
980 
981  AliExternalTrackParam exParam;
982 
983  //get covariance matrix
984  Double_t cov[21] = {0,};
985  aodtrack->GetCovMatrix(cov);
986  Double_t pxpypz[3] = {0,};
987  aodtrack->PxPyPz(pxpypz);
988  Double_t xyz[3] = {0,};
989  aodtrack->GetXYZ(xyz);
990  Short_t sign = aodtrack->Charge();
991  exParam.Set(xyz,pxpypz,cov,sign);
992 
993  sigma = TMath::Sqrt(exParam.GetSigma1Pt2());
994 
995  }
996 
997  Int_t label = TMath::Abs(track->GetLabel());
998 
999  FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, type);
1000 
1001  if (fGeneratorLevel && label > 0) {
1002  AliAODMCParticle *part = fGeneratorLevel->GetAcceptMCParticleWithLabel(label);
1003  if (part) {
1004  if (part->GetGeneratorIndex() == 0) {
1005  Int_t pdg = TMath::Abs(part->PdgCode());
1006  // select charged pions, protons, kaons, electrons, muons
1007  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
1008  FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
1009  }
1010  }
1011  }
1012  }
1013  }
1014  else {
1015  AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
1016  }
1017  }
1018 
1019  if (fGeneratorLevel) {
1020  AliAODMCParticle* part;
1021  for (auto partIterator : fGeneratorLevel->accepted_momentum() ) {
1022  part = partIterator.second;
1023 
1024  Byte_t findable = 0;
1025 
1026  Int_t pdg = TMath::Abs(part->PdgCode());
1027  // select charged pions, protons, kaons, electrons, muons
1028  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
1029 
1030  FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), findable);
1031  }
1032  }
1033 }
1034 
1035 //________________________________________________________________________
1037 
1038  if (!fCaloCells) return;
1039 
1040  TString histname_energy = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
1041  TString histname_time = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
1042  TString histname_energyProf = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
1043  TString histname_timeProf = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
1044  TString histname_EnergyvsTime = TString::Format("%s/fHistCellEvsTime", fCaloCellsName.Data());
1045 
1046  const Int_t ncells = fCaloCells->GetNumberOfCells();
1047  for (Int_t pos = 0; pos < ncells; pos++) {
1048  Float_t amp = fCaloCells->GetAmplitude(pos);
1049  Float_t time = fCaloCells->GetTime(pos);
1050  Int_t absId = fCaloCells->GetCellNumber(pos);
1051 
1052  if (amp < fCellEnergyCut) continue;
1053 
1054  fHistManager.FillTH1(histname_energy, amp);
1055  fHistManager.FillTH1(histname_time, time);
1056 
1057  fHistManager.FillProfile(histname_energyProf, absId, amp);
1058  fHistManager.FillProfile(histname_timeProf, absId, time);
1059 
1060  fHistManager.FillTH2(histname_EnergyvsTime, time, amp);
1061  }
1062 }
1063 
1064 //________________________________________________________________________
1066 
1067  memset(fNTotClusters, 0, sizeof(Int_t)*3);
1068  for (Int_t i = 0; i < 3; i++) fLeadingCluster[i].SetPxPyPzE(0,0,0,0);
1069 
1070  TString histname;
1071  AliClusterContainer* clusters = 0;
1072  TIter nextClusColl(&fClusterCollArray);
1073  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
1074  // Cluster loop
1076  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1077 
1078  UInt_t rejectionReason = 0;
1079  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
1080  histname = TString::Format("%s/fHistClusterRejectionReason", clusters->GetArrayName().Data());
1081  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
1082  continue;
1083  }
1084 
1085  // Determine cluster type (EMCal/DCal/Phos) and record relevant eventQA info
1086  ClusterType clusType = kNA;
1087  if (it->second->IsEMCAL()) {
1088  Double_t phi = it->first.Phi_0_2pi();
1089  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
1090  if (isDcal == 0) {
1091  clusType = kEMCal;
1092  } else if (isDcal == 1) {
1093  clusType = kDCal;
1094  }
1095 
1096  if (fLeadingCluster[isDcal].E() < it->first.E()) fLeadingCluster[isDcal] = it->first;
1097  fNTotClusters[isDcal]++;
1098 
1099  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
1100  if (sm >=0 && sm < 20) {
1101  histname = TString::Format("%s/BySM/hEmcalClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1102  fHistManager.FillTH1(histname, it->second->E());
1103  }
1104  else {
1105  AliError(Form("Supermodule %d does not exist!", sm));
1106  }
1107 
1108  } else if (it->second->GetType() == AliVCluster::kPHOSNeutral){
1109 
1110  Int_t nCells = it->second->GetNCells();
1111  Double_t M02 = it->second->GetM02();
1112  Double_t energy = it->second->E();
1113 
1114  histname = TString::Format("%s/hPhosNCellsVsEnergy", clusters->GetArrayName().Data());
1115  fHistManager.FillTH2(histname, nCells, energy);
1116 
1117  histname = TString::Format("%s/hPhosM02VsEnergy", clusters->GetArrayName().Data());
1118  fHistManager.FillTH2(histname, M02, energy);
1119 
1120  clusType = kPHOS;
1121  fNTotClusters[2]++;
1122  if (fLeadingCluster[2].E() < it->first.E()) fLeadingCluster[2] = it->first;
1123 
1124  // Fill Phos spectra by module
1125  Int_t relid[4];
1126  if (fPHOSGeo) {
1127  fPHOSGeo->AbsToRelNumbering(it->second->GetCellAbsId(0), relid);
1128  Int_t sm = relid[0];
1129  if (sm >=1 && sm < 5) {
1130  histname = TString::Format("%s/BySM/hPhosClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
1131  fHistManager.FillTH1(histname, energy);
1132  }
1133  else {
1134  AliError(Form("Supermodule %d does not exist!", sm));
1135  }
1136  }
1137 
1138  // Loop through cells in the cluster, and histogram their energy
1139  histname = TString::Format("%s/hPhosCellIdVsEnergy", clusters->GetArrayName().Data());
1140  for (Int_t i=0; i < nCells; i++) {
1141  fHistManager.FillTH2(histname, it->second->GetCellAbsId(i), energy);
1142  }
1143 
1144  }
1145 
1146  Double_t contents[30]={0};
1147  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
1148  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1149  if (!histClusterObservables) return;
1150  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
1151  TString title(histClusterObservables->GetAxis(i)->GetTitle());
1152  if (title=="Centrality %")
1153  contents[i] = fCent;
1154  else if (title=="#it{E}_{clus} (GeV)")
1155  contents[i] = it->first.E();
1156  else if (title=="#eta")
1157  contents[i] = it->first.Eta();
1158  else if (title=="#phi")
1159  contents[i] = it->first.Phi_0_2pi();
1160  else if (title=="cluster type")
1161  contents[i] = clusType;
1162  else
1163  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1164  }
1165  histClusterObservables->Fill(contents);
1166 
1167  }
1168  }
1169 }
1170 
1171 //________________________________________________________________________
1173 
1174  TString histname;
1175  AliJetContainer* jets = 0;
1176  TIter nextJetColl(&fJetCollArray);
1177  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
1178  Double_t rhoVal = 0;
1179  if (jets->GetRhoParameter()) {
1180  rhoVal = jets->GetRhoVal();
1181  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
1182  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
1183  }
1184 
1185  for (auto jet : jets->all()) {
1186 
1187  UInt_t rejectionReason = 0;
1188  if (!jets->AcceptJet(jet, rejectionReason)) {
1189  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
1190  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
1191  continue;
1192  }
1193 
1194  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
1195  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
1196 
1197  TLorentzVector leadPart;
1198  jets->GetLeadingHadronMomentum(leadPart, jet);
1199 
1200  Double_t contents[30]={0};
1201  histname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
1202  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
1203  if (!histJetObservables) return;
1204  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
1205  TString title(histJetObservables->GetAxis(i)->GetTitle());
1206  if (title=="Centrality (%)")
1207  contents[i] = fCent;
1208  else if (title=="#eta_{jet}")
1209  contents[i] = jet->Eta();
1210  else if (title=="#phi_{jet} (rad)")
1211  contents[i] = jet->Phi_0_2pi();
1212  else if (title=="#it{p}_{T} (GeV/#it{c})")
1213  contents[i] = jet->Pt();
1214  else if (title=="#it{p}_{T}^{MC} (GeV/#it{c})")
1215  contents[i] = jet->MCPt();
1216  else if (title=="#it{p}_{T}^{corr} (GeV/#it{c})")
1217  contents[i] = corrPt;
1218  else if (title=="#it{p}_{T,particle}^{leading} (GeV/#it{c})")
1219  contents[i] = ptLeading;
1220  else
1221  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1222  }
1223  histJetObservables->Fill(contents);
1224 
1225  UInt_t jetType = jet->GetJetAcceptanceType();
1226  if(jetType & AliEmcalJet::kEMCALfid)
1227  {
1228  histname = TString::Format("%s/hNEFVsPtEMC", jets->GetArrayName().Data());
1229  if (fForceBeamType != kpp){fHistManager.FillTH3(histname.Data(), fCent, jet->Pt(), jet->NEF());}
1230  else {fHistManager.FillTH2(histname.Data(), jet->Pt(), jet->NEF());}
1231  }
1232  else if(jetType & AliEmcalJet::kDCALonlyfid)
1233  {
1234  histname = TString::Format("%s/hNEFVsPtDCal", jets->GetArrayName().Data());
1235  if (fForceBeamType != kpp){
1236  fHistManager.FillTH3(histname.Data(), fCent, jet->Pt(), jet->NEF());
1237  }
1238  else{
1239  fHistManager.FillTH2(histname.Data(), jet->Pt(), jet->NEF());
1240  }
1241  }
1242  Double_t deltaEhadcorr = 0;
1243  const AliVCluster* clus = nullptr;
1244  Int_t nClusters = jet->GetNumberOfClusters();
1245  for (Int_t iClus = 0; iClus < nClusters; iClus++) {
1246  clus = jet->Cluster(iClus);
1247  deltaEhadcorr += (clus->GetNonLinCorrEnergy() - clus->GetHadCorrEnergy());
1248  }
1249 
1250  histname = TString::Format("%s/hDeltaEHadCorr", jets->GetArrayName().Data());
1251  if (fForceBeamType != kpp){
1252  fHistManager.FillTH3(histname, fCent, jet->Pt(), deltaEhadcorr);
1253  }
1254  else{
1255  fHistManager.FillTH2(histname, jet->Pt(), deltaEhadcorr);
1256  }
1257  } //jet loop
1258  }
1259 }
1260 
1261 //________________________________________________________________________
1263 
1264  EventQA_t eventQA;
1265  for (Int_t i = 0; i < 3; i++) {
1266  eventQA.fMaxCluster[i] = fLeadingCluster[i];
1267  }
1268  eventQA.fCent = fCent;
1269  eventQA.fNTracks = fNTotTracks;
1270  eventQA.fNClusters[0] = fNTotClusters[0];
1271  eventQA.fNClusters[1] = fNTotClusters[1];
1272  eventQA.fNClusters[2] = fNTotClusters[2];
1273  eventQA.fMaxTrack = fLeadingTrack;
1274 
1275  Int_t globalNclusters = eventQA.fNClusters[0] + eventQA.fNClusters[1] + eventQA.fNClusters[2];
1276  AliTLorentzVector globalMaxCluster;
1277  for (Int_t i = 0; i < 3; i++) {
1278  if (globalMaxCluster.E() < eventQA.fMaxCluster[i].E()) globalMaxCluster = eventQA.fMaxCluster[i];
1279  }
1280 
1281  THnSparse* histEventQA = static_cast<THnSparse*>(fHistManager.FindObject("eventQA"));
1282  Double_t contents[40]={0};
1283  for (Int_t i = 0; i < histEventQA->GetNdimensions(); i++) {
1284  TString title(histEventQA->GetAxis(i)->GetTitle());
1285  if (title=="Centrality %")
1286  contents[i] = eventQA.fCent;
1287  else if (title=="No. of tracks")
1288  contents[i] = eventQA.fNTracks;
1289  else if (title=="No. of clusters")
1290  contents[i] = globalNclusters;
1291  else if (title=="#it{p}_{T,track}^{leading} (GeV/c)")
1292  contents[i] = eventQA.fMaxTrack.Pt();
1293  else if (title=="#it{E}_{cluster}^{leading} (GeV)")
1294  contents[i] = globalMaxCluster.E();
1295  else
1296  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1297  }
1298 
1299  histEventQA->Fill(contents);
1300 
1301  // Fill pythia pt hard histograms, if applicable
1302  if (fPtHard > 1e-6) {
1303  fHistManager.FillTH1("hPtHard", fPtHard);
1304  }
1305 }
1306 
1307 //________________________________________________________________________
1309  Double_t sigma1OverPt, Byte_t trackType)
1310 {
1311  Double_t contents[20]={0};
1312 
1313  THnSparse* thnTracks = static_cast<THnSparse*>(fHistManager.FindObject("tracks"));
1314  if (!thnTracks) return;
1315  for (Int_t i = 0; i < thnTracks->GetNdimensions(); i++) {
1316  TString title(thnTracks->GetAxis(i)->GetTitle());
1317  if (title=="Centrality %")
1318  contents[i] = cent;
1319  else if (title=="#it{p}_{T} (GeV/#it{c})")
1320  contents[i] = trackPt;
1321  else if (title=="#eta")
1322  contents[i] = trackEta;
1323  else if (title=="#phi")
1324  contents[i] = trackPhi;
1325  else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
1326  contents[i] = sigma1OverPt;
1327  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
1328  contents[i] = sigma1OverPt*trackPt;
1329  else if (title=="track type")
1330  contents[i] = trackType;
1331  else
1332  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks->GetName()));
1333  }
1334 
1335  thnTracks->Fill(contents);
1336 }
1337 
1338 //________________________________________________________________________
1339 void AliAnalysisTaskPWGJEQA::FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Byte_t findable)
1340 {
1341  Double_t contents[20]={0};
1342 
1343  THnSparse* thnTracks_PhysPrim = static_cast<THnSparse*>(fHistManager.FindObject("tracks_PhysPrim"));
1344  if (!thnTracks_PhysPrim) return;
1345  for (Int_t i = 0; i < thnTracks_PhysPrim->GetNdimensions(); i++) {
1346  TString title(thnTracks_PhysPrim->GetAxis(i)->GetTitle());
1347  if (title=="Centrality %")
1348  contents[i] = cent;
1349  else if (title=="#it{p}_{T} (GeV/#it{c})")
1350  contents[i] = partPt;
1351  else if (title=="#eta")
1352  contents[i] = partEta;
1353  else if (title=="#phi")
1354  contents[i] = partPhi;
1355  else if (title=="Findable")
1356  contents[i] = findable;
1357  else
1358  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_PhysPrim->GetName()));
1359  }
1360 
1361  thnTracks_PhysPrim->Fill(contents);
1362 }
1363 
1364 //________________________________________________________________________
1366  Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
1367 {
1368  Double_t contents[20]={0};
1369 
1370  THnSparse* thnTracks_Matched = static_cast<THnSparse*>(fHistManager.FindObject("tracks_Matched"));
1371  if (!thnTracks_Matched) return;
1372  for (Int_t i = 0; i < thnTracks_Matched->GetNdimensions(); i++) {
1373  TString title(thnTracks_Matched->GetAxis(i)->GetTitle());
1374  if (title=="Centrality %")
1375  contents[i] = cent;
1376  else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
1377  contents[i] = partPt;
1378  else if (title=="#eta^{gen}")
1379  contents[i] = partEta;
1380  else if (title=="#phi^{gen}")
1381  contents[i] = partPhi;
1382  else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
1383  contents[i] = trackPt;
1384  else if (title=="#eta^{det}")
1385  contents[i] = trackEta;
1386  else if (title=="#phi^{det}")
1387  contents[i] = trackPhi;
1388  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
1389  contents[i] = (partPt - trackPt) / partPt;
1390  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}")
1391  contents[i] = (partPt - trackPt) / trackPt;
1392  else if (title=="track type")
1393  contents[i] = (Double_t)trackType;
1394  else
1395  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), thnTracks_Matched->GetName()));
1396  }
1397 
1398  thnTracks_Matched->Fill(contents);
1399 }
1400 
1407  const char* ntracks,
1408  const char* nclusters,
1409  const char* ncells,
1410  const char *nGenLev,
1411  Bool_t doTrackQA,
1412  Bool_t doCaloQA,
1413  Bool_t doJetQA,
1414  Bool_t doEventQA,
1415  Double_t trackPtCut,
1416  Double_t clusECut,
1417  const char* suffix)
1418 {
1419  // Get the pointer to the existing analysis manager via the static access method
1420  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1421  if (!mgr) {
1422  ::Error("AddTaskPWGJEQA", "No analysis manager to connect to.");
1423  return NULL;
1424  }
1425 
1426  // Check the analysis type using the event handlers connected to the analysis manager
1427  AliVEventHandler* handler = mgr->GetInputEventHandler();
1428  if (!handler) {
1429  ::Error("AddTaskPWGJEQA", "This task requires an input event handler");
1430  return NULL;
1431  }
1432 
1433  enum EDataType_t {
1434  kUnknown,
1435  kESD,
1436  kAOD
1437  };
1438 
1439  EDataType_t dataType = kUnknown;
1440 
1441  if (handler->InheritsFrom("AliESDInputHandler")) {
1442  dataType = kESD;
1443  }
1444  else if (handler->InheritsFrom("AliAODInputHandler")) {
1445  dataType = kAOD;
1446  }
1447 
1448  // Init the task and do settings
1449  TString trackName(ntracks);
1450  TString clusName(nclusters);
1451  TString cellName(ncells);
1452 
1453  if (trackName == "usedefault") {
1454  if (dataType == kESD) {
1455  trackName = "Tracks";
1456  }
1457  else if (dataType == kAOD) {
1458  trackName = "tracks";
1459  }
1460  else {
1461  trackName = "";
1462  }
1463  }
1464 
1465  if (clusName == "usedefault") {
1466  if (dataType == kESD) {
1467  clusName = "CaloClusters";
1468  }
1469  else if (dataType == kAOD) {
1470  clusName = "caloClusters";
1471  }
1472  else {
1473  clusName = "";
1474  }
1475  }
1476 
1477  if (cellName == "usedefault") {
1478  if (dataType == kESD) {
1479  cellName = "EMCALCells";
1480  }
1481  else if (dataType == kAOD) {
1482  cellName = "emcalCells";
1483  }
1484  else {
1485  cellName = "";
1486  }
1487  }
1488 
1489  TString name("AliAnalysisTaskPWGJEQA");
1490  if (!trackName.IsNull()) {
1491  name += "_";
1492  name += trackName;
1493  }
1494  if (!clusName.IsNull()) {
1495  name += "_";
1496  name += clusName;
1497  }
1498 
1499  if (!cellName.IsNull()) {
1500  name += "_";
1501  name += cellName;
1502  }
1503 
1504  if (strcmp(suffix,"")) {
1505  name += "_";
1506  name += suffix;
1507  }
1508 
1509  if (nGenLev && strcmp(nGenLev,"")!=0) cout<<"MC Part: "<< nGenLev<<endl;
1510  else cout<<"No MC Part: "<< nGenLev<<endl;
1511 
1512  AliAnalysisTaskPWGJEQA* qaTask = new AliAnalysisTaskPWGJEQA(name);
1513  qaTask->SetVzRange(-10,10);
1514  qaTask->SetNeedEmcalGeom(kFALSE);
1515  qaTask->SetCaloCellsName(cellName);
1516  qaTask->SetDetectorLevelName(trackName);
1517  if (nGenLev && strcmp(nGenLev,"")!=0) qaTask->SetGeneratorLevelName(nGenLev);
1518  qaTask->SetDoTrackQA(doTrackQA);
1519  qaTask->SetDoCaloQA(doCaloQA);
1520  qaTask->SetDoJetQA(doJetQA);
1521  qaTask->SetDoEventQA(doEventQA);
1522 
1523  // Add the detector-level track container
1524  if (trackName == "tracks" || trackName == "Tracks") {
1525  AliTrackContainer* trackCont = qaTask->AddTrackContainer(trackName);
1526  trackCont->SetFilterHybridTracks(kTRUE);
1527  }
1528  else if (!trackName.IsNull()) {
1529  qaTask->AddParticleContainer(trackName);
1530  }
1531  AliParticleContainer *partCont = qaTask->GetParticleContainer(trackName);
1532  if (partCont) {
1533  partCont->SetParticlePtCut(trackPtCut);
1534  }
1535 
1536  // Add the generator-level container, if specified
1537  if (nGenLev && strcmp(nGenLev,"")!=0) {
1538  AliMCParticleContainer* mcpartCont = qaTask->AddMCParticleContainer(nGenLev);
1539  mcpartCont->SelectPhysicalPrimaries(kTRUE);
1540  mcpartCont->SetParticlePtCut(0);
1541  }
1542 
1543  // Add the cluster container
1544  AliClusterContainer *clusterCont = qaTask->AddClusterContainer(clusName);
1545  if (clusterCont) {
1546  clusterCont->SetClusECut(0.);
1547  clusterCont->SetClusPtCut(0.);
1548  clusterCont->SetClusHadCorrEnergyCut(clusECut);
1549  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1550  }
1551 
1552  // Final settings, pass to manager and set the containers
1553  mgr->AddTask(qaTask);
1554 
1555  // Create containers for input/output
1556  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
1557 
1558  TString contName = TString::Format("%s_histos", name.Data());
1559  TString commonoutput;
1560  commonoutput = mgr->GetCommonFileName();
1561  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(),
1562  TList::Class(),AliAnalysisManager::kOutputContainer,
1563  commonoutput);
1564  mgr->ConnectInput (qaTask, 0, cinput1 );
1565  mgr->ConnectOutput (qaTask, 1, coutput1 );
1566 
1567  return qaTask;
1568 
1569 }
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:45
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)
static UShort_t GetRejectionReasonBitPosition(UInt_t rejectionReason)
Returns the highest bit in the rejection map as reason why the object was rejected.
TString fDetectorLevelName
detector level container name
void SetVzRange(Double_t min, Double_t max)
Set pre-configured event cut object.
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
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
Int_t GetCurrentID() const
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
Base class for container structures within the EMCAL framework.
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 TString & GetArrayName() const
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)
Enable general histograms.
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)
TTree * tree
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