AliPhysics  cda3415 (cda3415)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskPWGJEQA.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <cstring>
17 
18 #include <TClonesArray.h>
19 #include <TH1F.h>
20 #include <TH2F.h>
21 #include <TH3F.h>
22 #include <THnSparse.h>
23 #include <TMath.h>
24 #include <TString.h>
25 #include <Riostream.h>
26 
27 #include <AliVEventHandler.h>
28 #include <AliAnalysisManager.h>
29 #include "AliParticleContainer.h"
30 #include "AliClusterContainer.h"
31 #include "AliTrackContainer.h"
32 #include "AliCentrality.h"
33 #include "AliVCluster.h"
34 #include "AliVParticle.h"
35 #include "AliVTrack.h"
36 #include "AliLog.h"
37 #include "AliEMCALGeometry.h"
38 #include "AliEMCALGeoParams.h"
39 #include "AliESDtrack.h"
40 #include "AliAODTrack.h"
41 #include "AliAODMCParticle.h"
42 #include "AliMCParticleContainer.h"
43 #include "AliEmcalJet.h"
44 #include "AliJetContainer.h"
45 
46 #include "AliAnalysisTaskPWGJEQA.h"
47 
51 
52 //________________________________________________________________________
55  AliAnalysisTaskEmcalJet("AliAnalysisTaskPWGJEQA", kTRUE),
56  fCellEnergyCut(0.05),
57  fPtBinWidth(0.5),
58  fMaxPt(250),
59  fSeparateEMCalDCal(kTRUE),
60  fNTotTracks(0),
61  fLeadingTrack(),
62  fGeneratorLevel(0),
63  fGeneratorLevelName(),
64  fDetectorLevel(0),
65  fDetectorLevelName(),
66  fNPtHistBins(0),
67  fPtHistBins(0),
68  fNEtaHistBins(0),
69  fEtaHistBins(0),
70  fNPhiHistBins(0),
71  fPhiHistBins(0),
72  fNCentHistBins(0),
73  fCentHistBins(0),
74  fNPtRelDiffHistBins(0),
75  fPtRelDiffHistBins(0),
76  fNPtResHistBins(0),
77  fPtResHistBins(0),
78  f1OverPtResHistBins(0),
79  fN1OverPtResHistBins(0),
80  fNIntegerHistBins(0),
81  fIntegerHistBins(0),
82  fTracks(0),
83  fParticlesPhysPrim(0),
84  fParticlesMatched(0),
85  fHistManager("AliAnalysisTaskPWGJEQA"),
86  fDoTrackQA(kTRUE),
87  fDoEmcalQA(kTRUE),
88  fDoJetQA(kTRUE),
89  fDoEventQA(kTRUE)
90 {
91  // Default constructor.
92  memset(fNTotClusters, 0, sizeof(Int_t)*3);
95 }
96 
97 //________________________________________________________________________
102  AliAnalysisTaskEmcalJet(name, kTRUE),
103  fCellEnergyCut(0.05),
104  fPtBinWidth(0.5),
105  fMaxPt(250),
106  fSeparateEMCalDCal(kTRUE),
107  fNTotTracks(0),
108  fLeadingTrack(),
109  fGeneratorLevel(0),
110  fGeneratorLevelName(),
111  fDetectorLevel(0),
112  fDetectorLevelName(),
113  fNPtHistBins(0),
114  fPtHistBins(0),
115  fNEtaHistBins(0),
116  fEtaHistBins(0),
117  fNPhiHistBins(0),
118  fPhiHistBins(0),
119  fNCentHistBins(0),
120  fCentHistBins(0),
121  fNPtRelDiffHistBins(0),
122  fPtRelDiffHistBins(0),
123  fNPtResHistBins(0),
124  fPtResHistBins(0),
125  f1OverPtResHistBins(0),
126  fN1OverPtResHistBins(0),
127  fNIntegerHistBins(0),
128  fIntegerHistBins(0),
129  fTracks(0),
130  fParticlesPhysPrim(0),
131  fParticlesMatched(0),
132  fHistManager(name),
133  fDoTrackQA(kTRUE),
134  fDoEmcalQA(kTRUE),
135  fDoJetQA(kTRUE),
136  fDoEventQA(kTRUE)
137 {
138  // Standard
139  memset(fNTotClusters, 0, sizeof(Int_t)*3);
142 }
143 
144 //________________________________________________________________________
146 {
147  // Destructor
148 }
149 
150 //________________________________________________________________________
152 {
153  // Create histograms
155 
156  // Set track container pointers
159 
160  // Allocate histograms for tracks, cells, clusters, jets
166 
167  TIter next(fHistManager.GetListOfHistograms());
168  TObject* obj = 0;
169  while ((obj = next())) {
170  if (!strcmp(obj->GetName(),"fHistCellsAbsIdEnergy") || !strcmp(obj->GetName(),"fHistCellsAbsIdTime")) {
171  continue;
172  }
173  fOutput->Add(obj);
174  }
175 
176  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
177 }
178 
179 //________________________________________________________________________
181 
183  if (fGeneratorLevel) {
186  }
187 }
188 
189 //________________________________________________________________________
191 
192  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
193  TString histname;
194  TString title;
195 
196  if (!fCaloCellsName.IsNull()) {
197 
198  // TH2s used to create the cell histograms (we do not write these to the output)
199  histname = "fHistCellsAbsIdEnergy";
200  title = histname + ";cell abs. Id;#it{E}_{cell} (GeV);counts";
201  fHistManager.CreateTH2(histname.Data(), title.Data(), 20000,0,20000,(Int_t)(nPtBins / 2), 0, fMaxPt / 2);
202 
203  histname = "fHistCellsAbsIdTime";
204  title = histname + ";cell abs. Id;#it{t}_{cell} (s);counts";
205  fHistManager.CreateTH2(histname.Data(), title.Data(), 20000,0,20000,(Int_t)(nPtBins / 2), -5e-6, 5e-6);
206 
207  // TH1s that we write to the output
208 
209  histname = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
210  title = histname + ";#it{E}_{cell} (GeV);counts";
211  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,0, 100);
212 
213  /*
214  histname = TString::Format("%s/fProfCellAbsIdEnergy", fCaloCellsName.Data());
215  title = histname;
216  fHistManager.CreateTH1(histname.Data(), title.Data(), 20000,0,20000);
217  */
218 
219  histname = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
220  title = histname + ";#it{t}_{cell} (s);counts";
221  fHistManager.CreateTH1(histname.Data(), title.Data(), 500,-5e-6, 5e-6);
222 
223  /*
224  histname = TString::Format("%s/fProfCellAbsIdTime", fCaloCellsName.Data());
225  title = histname + ";#it{t}_{cell} (s);counts";
226  fHistManager.CreateTH1(histname.Data(), title.Data(), 20000,0,20000);
227  */
228  }
229 
230 }
231 
232 //________________________________________________________________________
234 
235  AliEmcalContainer* cont = 0;
236  TString histname;
237  TString title;
238  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
239 
240  TIter nextClusColl(&fClusterCollArray);
241  while ((cont = static_cast<AliEmcalContainer*>(nextClusColl()))) {
242 
243  histname = TString::Format("%s/fHistClusterRejectionReason", cont->GetArrayName().Data());
244  title = histname + ";Rejection reason;#it{E}_{clus} (GeV/);counts";
245  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
246  SetRejectionReasonLabels(hist->GetXaxis());
247 
248  const Int_t nSM = 20;
249  for (Int_t sm = 0; sm < nSM; sm++) {
250  histname = TString::Format("%s/BySM/fHistClusEnergy_SM%d", cont->GetArrayName().Data(), sm);
251  title = histname + ";#it{E}_{cluster} (GeV);counts";
252  fHistManager.CreateTH1(histname.Data(), title.Data(), nPtBins, 0, fMaxPt);
253  }
254 
255  // Cluster THnSparse
256  Int_t dim = 0;
257  TString title[20];
258  Int_t nbins[20] = {0};
259  Double_t min[30] = {0.};
260  Double_t max[30] = {0.};
261 
263  title[dim] = "Centrality %";
264  nbins[dim] = 5;
265  min[dim] = 0;
266  max[dim]= 100;
267  dim++;
268  }
269 
270  title[dim] = "#it{E}_{clus} (GeV)";
271  nbins[dim] = fNPtHistBins;
272  min[dim] = 0;
273  max[dim] = 250;
274  dim++;
275 
276  title[dim] = "#eta";
277  nbins[dim] = fNEtaHistBins;
278  min[dim] = -1;
279  max[dim] = 1;
280  dim++;
281 
282  title[dim] = "#phi";
283  nbins[dim] = fNPhiHistBins;
284  min[dim] = 0;
285  max[dim] = 2*TMath::Pi();
286  dim++;
287 
288  title[dim] = "cluster type";
289  nbins[dim] = 3;
290  min[dim] = -0.5;
291  max[dim] = 2.5;
292  dim++;
293 
294  TString thnname = TString::Format("%s/clusterObservables", cont->GetArrayName().Data());
295  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
296  for (Int_t i = 0; i < dim; i++) {
297  hn->GetAxis(i)->SetTitle(title[i]);
298  }
299  }
300 }
301 
302 //________________________________________________________________________
304 
305  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
306 
307  AliJetContainer* jets = 0;
308  TIter nextJetColl(&fJetCollArray);
309  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
310 
311  // Allocate THnSparse
312  Double_t jetRadius = jets->GetJetRadius();
313 
314  TString axisTitle[30]= {""};
315  Int_t nbins[30] = {0};
316  Double_t min[30] = {0.};
317  Double_t max[30] = {0.};
318  Int_t dim = 0;
319 
320  if (fForceBeamType != kpp) {
321  axisTitle[dim] = "Centrality (%)";
322  nbins[dim] = 5;
323  min[dim] = 0;
324  max[dim] = 100;
325  dim++;
326  }
327 
328  axisTitle[dim] = "#eta_{jet}";
329  nbins[dim] = nPtBins/10;
330  min[dim] = -1;
331  max[dim] = 1;
332  dim++;
333 
334  axisTitle[dim] = "#phi_{jet} (rad)";
335  nbins[dim] = nPtBins/10*3;
336  min[dim] = 0;
337  max[dim] = 2*TMath::Pi();
338  dim++;
339 
340  axisTitle[dim] = "#it{p}_{T} (GeV/#it{c})";
341  nbins[dim] = nPtBins/2;
342  min[dim] = 0;
343  max[dim] = fMaxPt;
344  dim++;
345 
346  if (fForceBeamType != kpp) {
347  axisTitle[dim] = "#it{p}_{T}^{corr} (GeV/#it{c})";
348  nbins[dim] = nPtBins/2;
349  min[dim] = -fMaxPt/2;
350  max[dim] = fMaxPt/2;
351  dim++;
352  }
353 
354  if (fClusterCollArray.GetEntriesFast() > 0 && fParticleCollArray.GetEntriesFast() > 0) {
355  axisTitle[dim] = "NEF";
356  nbins[dim] = nPtBins/20;
357  min[dim] = 0;
358  max[dim] = 1.0;
359  dim++;
360  }
361 
362  if (fForceBeamType != kpp) {
363  axisTitle[dim] = "No. of constituents";
364  nbins[dim] = 125;
365  min[dim] = 0;
366  max[dim] = 250;
367  dim++;
368  }
369  else {
370  axisTitle[dim] = "No. of constituents";
371  nbins[dim] = 50;
372  min[dim] = -0.5;
373  max[dim] = 49.5;
374  dim++;
375  }
376 
377  axisTitle[dim] = "#it{p}_{T,particle}^{leading} (GeV/#it{c})";
378  nbins[dim] = nPtBins/10*3;
379  min[dim] = 0;
380  max[dim] = 150;
381  dim++;
382 
383  TString thnname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
384  THnSparse* hn = fHistManager.CreateTHnSparse(thnname.Data(), thnname.Data(), dim, nbins, min, max);
385  for (Int_t i = 0; i < dim; i++) {
386  hn->GetAxis(i)->SetTitle(axisTitle[i]);
387  }
388 
389  // Allocate other jet histograms
390  TString histname;
391  TString title;
392  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
393  title = histname + ";Rejection reason;#it{p}_{T,jet} (GeV/#it{c});counts";
394  TH2* hist = fHistManager.CreateTH2(histname.Data(), title.Data(), 32, 0, 32, 50, 0, 250);
395  SetRejectionReasonLabels(hist->GetXaxis());
396 
397  if (!jets->GetRhoName().IsNull()) {
398  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
399  title = histname + ";Centrality (%);#rho (GeV/#it{c});counts";
400  fHistManager.CreateTH2(histname.Data(), title.Data(), 101, 0, 101, 100, 0, 500);
401  }
402  }
403 }
404 
405 //________________________________________________________________________
407 
408  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
409 
410  Int_t dim = 0;
411  TString axistitle[40];
412  Int_t nbins[40] = {0};
413  Double_t min[40] = {0};
414  Double_t max[40] = {0};
415 
417  axistitle[dim] = "Centrality %";
418  nbins[dim] = 5;
419  min[dim] = 0;
420  max[dim] = 100;
421  dim++;
422  }
423 
424  if (fParticleCollArray.GetEntriesFast()>0) {
425  axistitle[dim] = "No. of tracks";
427  nbins[dim] = 100;
428  min[dim] = -0.5;
429  max[dim] = 6000-0.5;
430  }
431  else {
432  nbins[dim] = 100;
433  min[dim] = 0;
434  max[dim] = 200;
435  }
436  dim++;
437 
438  axistitle[dim] = "#it{p}_{T,track}^{leading} (GeV/c)";
439  nbins[dim] = nPtBins/2;
440  min[dim] = 0;
441  max[dim] = fMaxPt;
442  dim++;
443  }
444 
445  if (fClusterCollArray.GetEntriesFast()>0) {
446  axistitle[dim] = "No. of clusters";
447 
449  nbins[dim] = 100;
450  min[dim] = -0.5;
451  max[dim] = 4000-0.5;
452  }
453  else {
454  nbins[dim] = 100;
455  min[dim] = 0;
456  max[dim] = 200;
457  }
458  dim++;
459 
460  if (fSeparateEMCalDCal) {
461  axistitle[dim] = "#it{E}_{EMCal cluster}^{leading} (GeV)";
462  nbins[dim] = nPtBins/2;
463  min[dim] = 0;
464  max[dim] = fMaxPt;
465  dim++;
466 
467  axistitle[dim] = "#it{E}_{DCal cluster}^{leading} (GeV)";
468  nbins[dim] = nPtBins/2;
469  min[dim] = 0;
470  max[dim] = fMaxPt;
471  dim++;
472 
473  axistitle[dim] = "#it{E}_{PHOS cluster}^{leading} (GeV)";
474  nbins[dim] = nPtBins/2;
475  min[dim] = 0;
476  max[dim] = fMaxPt;
477  dim++;
478  }
479  else {
480  axistitle[dim] = "#it{E}_{cluster}^{leading} (GeV)";
481  nbins[dim] = nPtBins/2;
482  min[dim] = 0;
483  max[dim] = fMaxPt;
484  dim++;
485  }
486  }
487 
488  THnSparse* hn = fHistManager.CreateTHnSparse("eventQA","eventQA",dim,nbins,min,max);
489  for (Int_t i = 0; i < dim; i++)
490  hn->GetAxis(i)->SetTitle(axistitle[i]);
491 
492 }
493 
494 //________________________________________________________________________
497 {
498  fNPtHistBins = 82;
501  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
502  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
503  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
504  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
505  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
506  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
507 
508  fNEtaHistBins = 100;
511 
512  fNPhiHistBins = 101;
514  GenerateFixedBinArray(fNPhiHistBins, 0, TMath::Pi() * 2.02, fPhiHistBins);
515 
516  fNCentHistBins = 4;
518  fCentHistBins[0] = 0;
519  fCentHistBins[1] = 10;
520  fCentHistBins[2] = 30;
521  fCentHistBins[3] = 50;
522  fCentHistBins[4] = 90;
523 
524  fNPtResHistBins = 175;
527  GenerateFixedBinArray(25, 0.05, 0.10, fPtResHistBins+50);
528  GenerateFixedBinArray(25, 0.10, 0.20, fPtResHistBins+75);
529  GenerateFixedBinArray(30, 0.20, 0.50, fPtResHistBins+100);
530  GenerateFixedBinArray(25, 0.50, 1.00, fPtResHistBins+130);
531  GenerateFixedBinArray(20, 1.00, 2.00, fPtResHistBins+155);
532 
533  fNPtRelDiffHistBins = 200;
536 
537  fN1OverPtResHistBins = 385;
540  GenerateFixedBinArray(60, 0.02, 0.05, f1OverPtResHistBins+100);
541  GenerateFixedBinArray(50, 0.05, 0.1, f1OverPtResHistBins+160);
542  GenerateFixedBinArray(50, 0.1, 0.2, f1OverPtResHistBins+210);
543  GenerateFixedBinArray(75, 0.2, 0.5, f1OverPtResHistBins+260);
544  GenerateFixedBinArray(50, 0.5, 1.5, f1OverPtResHistBins+335);
545 
546  fNIntegerHistBins = 10;
549 }
550 
551 //________________________________________________________________________
553 {
554  Int_t dim = 0;
555  TString title[20];
556  Int_t nbins[20] = {0};
557  Double_t *binEdges[20] = {0};
558 
560  title[dim] = "Centrality %";
561  nbins[dim] = fNCentHistBins;
562  binEdges[dim] = fCentHistBins;
563  dim++;
564  }
565 
566  title[dim] = "#it{p}_{T} (GeV/#it{c})";
567  nbins[dim] = fNPtHistBins;
568  binEdges[dim] = fPtHistBins;
569  dim++;
570 
571  title[dim] = "#eta";
572  nbins[dim] = fNEtaHistBins;
573  binEdges[dim] = fEtaHistBins;
574  dim++;
575 
576  title[dim] = "#phi";
577  nbins[dim] = fNPhiHistBins;
578  binEdges[dim] = fPhiHistBins;
579  dim++;
580 
581  title[dim] = "track type";
582  nbins[dim] = 3;
583  binEdges[dim] = fIntegerHistBins;
584  dim++;
585 
586  title[dim] = "#sigma(#it{p}_{T}) / #it{p}_{T}";
587  nbins[dim] = fNPtResHistBins;
588  binEdges[dim] = fPtResHistBins;
589  dim++;
590 
591  fTracks = new THnSparseF("tracks","tracks",dim,nbins);
592  for (Int_t i = 0; i < dim; i++) {
593  fTracks->GetAxis(i)->SetTitle(title[i]);
594  fTracks->SetBinEdges(i, binEdges[i]);
595  }
596 
597  fOutput->Add(fTracks);
598 
599 }
600 
601 //________________________________________________________________________
603 {
604  Int_t dim = 0;
605  TString title[20];
606  Int_t nbins[20] = {0};
607  Double_t *binEdges[20] = {0};
608 
610  title[dim] = "Centrality %";
611  nbins[dim] = fNCentHistBins;
612  binEdges[dim] = fCentHistBins;
613  dim++;
614  }
615 
616  title[dim] = "#it{p}_{T} (GeV/#it{c})";
617  nbins[dim] = fNPtHistBins;
618  binEdges[dim] = fPtHistBins;
619  dim++;
620 
621  title[dim] = "#eta";
622  nbins[dim] = fNEtaHistBins;
623  binEdges[dim] = fEtaHistBins;
624  dim++;
625 
626  title[dim] = "#phi";
627  nbins[dim] = fNPhiHistBins;
628  binEdges[dim] = fPhiHistBins;
629  dim++;
630 
631  fParticlesPhysPrim = new THnSparseF("tracks_PhysPrim","tracks_PhysPrim",dim,nbins);
632  for (Int_t i = 0; i < dim; i++) {
633  fParticlesPhysPrim->GetAxis(i)->SetTitle(title[i]);
634  fParticlesPhysPrim->SetBinEdges(i, binEdges[i]);
635  }
636 
638 }
639 
640 //________________________________________________________________________
642 {
643  Int_t dim = 0;
644  TString title[20];
645  Int_t nbins[20] = {0};
646  Double_t *binEdges[20] = {0};
647 
649  title[dim] = "Centrality %";
650  nbins[dim] = fNCentHistBins;
651  binEdges[dim] = fCentHistBins;
652  dim++;
653  }
654 
655  title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
656  nbins[dim] = fNPtHistBins;
657  binEdges[dim] = fPtHistBins;
658  dim++;
659 
660  title[dim] = "#eta^{gen}";
661  nbins[dim] = fNEtaHistBins;
662  binEdges[dim] = fEtaHistBins;
663  dim++;
664 
665  title[dim] = "#phi^{gen}";
666  nbins[dim] = fNPhiHistBins;
667  binEdges[dim] = fPhiHistBins;
668  dim++;
669 
670  title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
671  nbins[dim] = fNPtHistBins;
672  binEdges[dim] = fPtHistBins;
673  dim++;
674 
675  title[dim] = "#eta^{det}";
676  nbins[dim] = fNEtaHistBins;
677  binEdges[dim] = fEtaHistBins;
678  dim++;
679 
680  title[dim] = "#phi^{det}";
681  nbins[dim] = fNPhiHistBins;
682  binEdges[dim] = fPhiHistBins;
683  dim++;
684 
685  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
686  nbins[dim] = fNPtRelDiffHistBins;
687  binEdges[dim] = fPtRelDiffHistBins;
688  dim++;
689 
690  title[dim] = "track type";
691  nbins[dim] = 3;
692  binEdges[dim] = fIntegerHistBins;
693  dim++;
694 
695  fParticlesMatched = new THnSparseF("tracks_Matched","tracks_Matched",dim,nbins);
696  for (Int_t i = 0; i < dim; i++) {
697  fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
698  fParticlesMatched->SetBinEdges(i, binEdges[i]);
699  }
700 
702 }
703 
704 //________________________________________________________________________
706 {
707  if (fClusterCollArray.GetEntriesFast() == 0 && fCaloCellsName.IsNull()) {
708  fNeedEmcalGeom = kFALSE;
709  }
710  else {
711  fNeedEmcalGeom = kTRUE;
712  }
713 
715 }
716 
717 //________________________________________________________________________
719 {
720  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
721 
722  return kTRUE;
723 }
724 
725 //________________________________________________________________________
727 {
733 
734  return kTRUE;
735 }
736 
737 //________________________________________________________________________
739 
740  fNTotTracks = 0;
741  fLeadingTrack.SetPxPyPzE(0,0,0,0);
742 
743  fDetectorLevel->ResetCurrentID();
744  AliVTrack* track;
745  for (auto trackIterator : fDetectorLevel->accepted_momentum() ) {
746 
747  fNTotTracks++;
748  if (fLeadingTrack.Pt() < trackIterator.first.Pt()) fLeadingTrack = trackIterator.first;
749 
750  track = trackIterator.second;
751  Byte_t type = fDetectorLevel->GetTrackType(track);
752  if (type <= 2) {
753  Double_t sigma = 0;
754 
755  if (fIsEsd) {
756 
757  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
758  if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
759 
760  } else { // AOD
761 
762  AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
763  if(!aodtrack) AliFatal("Not a standard AOD");
764 
765  AliExternalTrackParam exParam;
766 
767  //get covariance matrix
768  Double_t cov[21] = {0,};
769  aodtrack->GetCovMatrix(cov);
770  Double_t pxpypz[3] = {0,};
771  aodtrack->PxPyPz(pxpypz);
772  Double_t xyz[3] = {0,};
773  aodtrack->GetXYZ(xyz);
774  Short_t sign = aodtrack->Charge();
775  exParam.Set(xyz,pxpypz,cov,sign);
776 
777  sigma = TMath::Sqrt(exParam.GetSigma1Pt2());
778 
779  }
780 
781  Int_t label = TMath::Abs(track->GetLabel());
782  Int_t mcGen = 1;
783  // reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
784  if (label==0 || track->GetGeneratorIndex() == 0) mcGen = 0;
785 
786  FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, type);
787 
788  if (fGeneratorLevel && label > 0) {
789  AliAODMCParticle *part = fGeneratorLevel->GetAcceptMCParticleWithLabel(label);
790  if (part) {
791  if (part->GetGeneratorIndex() == 0) {
792  Int_t pdg = TMath::Abs(part->PdgCode());
793  // select charged pions, protons, kaons, electrons, muons
794  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
795  FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
796  }
797  }
798  }
799  }
800  }
801  else {
802  AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
803  }
804  }
805 
806  if (fGeneratorLevel) {
807  fGeneratorLevel->ResetCurrentID();
808  AliAODMCParticle* part;
809  for (auto partIterator : fGeneratorLevel->accepted_momentum() ) {
810  part = partIterator.second;
811 
812  Int_t mcGen = 1;
813  Byte_t findable = 0;
814 
815  if (part->GetGeneratorIndex() == 0) mcGen = 0;
816 
817  Int_t pdg = TMath::Abs(part->PdgCode());
818  // select charged pions, protons, kaons, electrons, muons
819  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
820 
821  FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt());
822  }
823  }
824 }
825 
826 //________________________________________________________________________
828 
829  if (!fCaloCells) return;
830 
831  TString histname_en2D = "fHistCellsAbsIdEnergy";
832  TString histname_tm2D = "fHistCellsAbsIdTime";
833  TString histname_energy = TString::Format("%s/fHistCellEnergy", fCaloCellsName.Data());
834  TString histname_time = TString::Format("%s/fHistCellTime", fCaloCellsName.Data());
835 
836  const Int_t ncells = fCaloCells->GetNumberOfCells();
837  for (Int_t pos = 0; pos < ncells; pos++) {
838  Float_t amp = fCaloCells->GetAmplitude(pos);
839  Float_t time = fCaloCells->GetTime(pos);
840  Int_t absId = fCaloCells->GetCellNumber(pos);
841 
842  if (amp < fCellEnergyCut) continue;
843 
844  fHistManager.FillTH2(histname_en2D, absId, amp);
845  fHistManager.FillTH2(histname_tm2D, absId, time);
846 
847  fHistManager.FillTH1(histname_energy, amp);
848  fHistManager.FillTH1(histname_time, time);
849  }
850 }
851 
852 //________________________________________________________________________
854 
855  memset(fNTotClusters, 0, sizeof(Int_t)*3);
856  for (Int_t i = 0; i < 3; i++) fLeadingCluster[i].SetPxPyPzE(0,0,0,0);
857 
858  TString histname;
859  AliClusterContainer* clusters = 0;
860  TIter nextClusColl(&fClusterCollArray);
861  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
862  // Cluster loop
864  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
865 
866  UInt_t rejectionReason = 0;
867  if (!clusters->AcceptCluster(it.current_index(), rejectionReason)) {
868  histname = TString::Format("%s/fHistClusterRejectionReason", clusters->GetArrayName().Data());
869  fHistManager.FillTH2(histname, clusters->GetRejectionReasonBitPosition(rejectionReason), it->first.E());
870  continue;
871  }
872 
873  // Determine cluster type (EMCal/DCal/Phos) and record relevant eventQA info
874  ClusterType clusType = kNA;
875  if (it->second->IsEMCAL()) {
876  Double_t phi = it->first.Phi_0_2pi();
877  Int_t isDcal = Int_t(phi > fgkEMCalDCalPhiDivide);
878  if (isDcal == 0) {
879  clusType = kEMCal;
880  } else if (isDcal == 1) {
881  clusType = kDCal;
882  }
883 
884  if (fLeadingCluster[isDcal].E() < it->first.E()) fLeadingCluster[isDcal] = it->first;
885  fNTotClusters[isDcal]++;
886 
887  Int_t sm = fGeom->GetSuperModuleNumber(it->second->GetCellAbsId(0));
888  if (sm >=0 && sm < 20) {
889  histname = TString::Format("%s/BySM/fHistClusEnergy_SM%d", clusters->GetArrayName().Data(), sm);
890  fHistManager.FillTH1(histname, it->second->E());
891  }
892  else {
893  AliError(Form("Supermodule %d does not exist!", sm));
894  }
895 
896  } else if (it->second->IsPHOS()){
897  clusType = kPHOS;
898  fNTotClusters[2]++;
899  if (fLeadingCluster[2].E() < it->first.E()) fLeadingCluster[2] = it->first;
900  }
901 
902  Double_t contents[30]={0};
903  histname = TString::Format("%s/clusterObservables", clusters->GetArrayName().Data());
904  THnSparse* histClusterObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
905  if (!histClusterObservables) return;
906  for (Int_t i = 0; i < histClusterObservables->GetNdimensions(); i++) {
907  TString title(histClusterObservables->GetAxis(i)->GetTitle());
908  if (title=="Centrality %")
909  contents[i] = fCent;
910  else if (title=="#it{E}_{clus} (GeV)")
911  contents[i] = it->first.E();
912  else if (title=="#eta")
913  contents[i] = it->first.Eta();
914  else if (title=="#phi")
915  contents[i] = it->first.Phi_0_2pi();
916  else if (title=="cluster type")
917  contents[i] = clusType;
918  else
919  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
920  }
921  histClusterObservables->Fill(contents);
922 
923  }
924  }
925 }
926 
927 //________________________________________________________________________
929 
930  TString histname;
931  AliJetContainer* jets = 0;
932  TIter nextJetColl(&fJetCollArray);
933  while ((jets = static_cast<AliJetContainer*>(nextJetColl()))) {
934  Double_t rhoVal = 0;
935  if (jets->GetRhoParameter()) {
936  rhoVal = jets->GetRhoVal();
937  histname = TString::Format("%s/fHistRhoVsCent", jets->GetArrayName().Data());
938  fHistManager.FillTH2(histname.Data(), fCent, rhoVal);
939  }
940 
941  for (auto jet : jets->all()) {
942 
943  UInt_t rejectionReason = 0;
944  if (!jets->AcceptJet(jet, rejectionReason)) {
945  histname = TString::Format("%s/fHistJetRejectionReason", jets->GetArrayName().Data());
946  fHistManager.FillTH2(histname.Data(), jets->GetRejectionReasonBitPosition(rejectionReason), jet->Pt());
947  continue;
948  }
949 
950  Float_t ptLeading = jets->GetLeadingHadronPt(jet);
951  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
952 
953  TLorentzVector leadPart;
954  jets->GetLeadingHadronMomentum(leadPart, jet);
955 
956  Double_t contents[30]={0};
957  histname = TString::Format("%s/fHistJetObservables", jets->GetArrayName().Data());
958  THnSparse* histJetObservables = static_cast<THnSparse*>(fHistManager.FindObject(histname));
959  if (!histJetObservables) return;
960  for (Int_t i = 0; i < histJetObservables->GetNdimensions(); i++) {
961  TString title(histJetObservables->GetAxis(i)->GetTitle());
962  if (title=="Centrality (%)")
963  contents[i] = fCent;
964  else if (title=="#eta_{jet}")
965  contents[i] = jet->Eta();
966  else if (title=="#phi_{jet} (rad)")
967  contents[i] = jet->Phi_0_2pi();
968  else if (title=="#it{p}_{T} (GeV/#it{c})")
969  contents[i] = jet->Pt();
970  else if (title=="#it{p}_{T}^{MC} (GeV/#it{c})")
971  contents[i] = jet->MCPt();
972  else if (title=="#it{p}_{T}^{corr} (GeV/#it{c})")
973  contents[i] = corrPt;
974  else if (title=="NEF")
975  contents[i] = jet->NEF();
976  else if (title=="No. of constituents")
977  contents[i] = jet->GetNumberOfConstituents();
978  else if (title=="#it{p}_{T,particle}^{leading} (GeV/#it{c})")
979  contents[i] = ptLeading;
980  else
981  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
982  }
983  histJetObservables->Fill(contents);
984 
985  } //jet loop
986  }
987 }
988 
989 //________________________________________________________________________
991 
992  EventQA_t eventQA;
993  for (Int_t i = 0; i < 3; i++) {
994  eventQA.fMaxCluster[i] = fLeadingCluster[i];
995  }
996  eventQA.fCent = fCent;
997  eventQA.fNTracks = fNTotTracks;
998  eventQA.fNClusters[0] = fNTotClusters[0];
999  eventQA.fNClusters[1] = fNTotClusters[1];
1000  eventQA.fNClusters[2] = fNTotClusters[2];
1001  eventQA.fMaxTrack = fLeadingTrack;
1002 
1003  Int_t globalNclusters = eventQA.fNClusters[0] + eventQA.fNClusters[1] + eventQA.fNClusters[2];
1004  AliTLorentzVector globalMaxCluster;
1005  for (Int_t i = 0; i < 3; i++) {
1006  if (globalMaxCluster.E() < eventQA.fMaxCluster[i].E()) globalMaxCluster = eventQA.fMaxCluster[i];
1007  }
1008 
1009  THnSparse* histEventQA = static_cast<THnSparse*>(fHistManager.FindObject("eventQA"));
1010  Double_t contents[40]={0};
1011  for (Int_t i = 0; i < histEventQA->GetNdimensions(); i++) {
1012  TString title(histEventQA->GetAxis(i)->GetTitle());
1013  if (title=="Centrality %")
1014  contents[i] = eventQA.fCent;
1015  else if (title=="No. of tracks")
1016  contents[i] = eventQA.fNTracks;
1017  else if (title=="No. of clusters")
1018  contents[i] = globalNclusters;
1019  else if (title=="#it{p}_{T,track}^{leading} (GeV/c)")
1020  contents[i] = eventQA.fMaxTrack.Pt();
1021  else if (title=="#it{E}_{cluster}^{leading} (GeV)")
1022  contents[i] = globalMaxCluster.E();
1023  else if (title=="#it{E}_{EMCal cluster}^{leading} (GeV)")
1024  contents[i] = eventQA.fMaxCluster[0].E();
1025  else if (title=="#it{E}_{DCal cluster}^{leading} (GeV)")
1026  contents[i] = eventQA.fMaxCluster[1].E();
1027  else if (title=="#it{E}_{PHOS cluster}^{leading} (GeV)")
1028  contents[i] = eventQA.fMaxCluster[2].E();
1029  else
1030  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1031  }
1032 
1033  histEventQA->Fill(contents);
1034 }
1035 
1036 //________________________________________________________________________
1038  Double_t sigma1OverPt, Byte_t trackType)
1039 {
1040  Double_t contents[20]={0};
1041 
1042  for (Int_t i = 0; i < fTracks->GetNdimensions(); i++) {
1043  TString title(fTracks->GetAxis(i)->GetTitle());
1044  if (title=="Centrality %")
1045  contents[i] = cent;
1046  else if (title=="#it{p}_{T} (GeV/#it{c})")
1047  contents[i] = trackPt;
1048  else if (title=="#eta")
1049  contents[i] = trackEta;
1050  else if (title=="#phi")
1051  contents[i] = trackPhi;
1052  else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
1053  contents[i] = sigma1OverPt;
1054  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
1055  contents[i] = sigma1OverPt*trackPt;
1056  else if (title=="track type")
1057  contents[i] = trackType;
1058  else
1059  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fTracks->GetName()));
1060  }
1061 
1062  fTracks->Fill(contents);
1063 }
1064 
1065 //________________________________________________________________________
1067 {
1068  Double_t contents[20]={0};
1069 
1070  for (Int_t i = 0; i < fParticlesPhysPrim->GetNdimensions(); i++) {
1071  TString title(fParticlesPhysPrim->GetAxis(i)->GetTitle());
1072  if (title=="Centrality %")
1073  contents[i] = cent;
1074  else if (title=="#it{p}_{T} (GeV/#it{c})")
1075  contents[i] = partPt;
1076  else if (title=="#eta")
1077  contents[i] = partEta;
1078  else if (title=="#phi")
1079  contents[i] = partPhi;
1080  else
1081  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesPhysPrim->GetName()));
1082  }
1083 
1084  fParticlesPhysPrim->Fill(contents);
1085 }
1086 
1087 //________________________________________________________________________
1089  Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
1090 {
1091  Double_t contents[20]={0};
1092 
1093  for (Int_t i = 0; i < fParticlesMatched->GetNdimensions(); i++) {
1094  TString title(fParticlesMatched->GetAxis(i)->GetTitle());
1095  if (title=="Centrality %")
1096  contents[i] = cent;
1097  else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
1098  contents[i] = partPt;
1099  else if (title=="#eta^{gen}")
1100  contents[i] = partEta;
1101  else if (title=="#phi^{gen}")
1102  contents[i] = partPhi;
1103  else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
1104  contents[i] = trackPt;
1105  else if (title=="#eta^{det}")
1106  contents[i] = trackEta;
1107  else if (title=="#phi^{det}")
1108  contents[i] = trackPhi;
1109  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
1110  contents[i] = (partPt - trackPt) / partPt;
1111  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}")
1112  contents[i] = (partPt - trackPt) / trackPt;
1113  else if (title=="track type")
1114  contents[i] = (Double_t)trackType;
1115  else
1116  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesMatched->GetName()));
1117  }
1118 
1119  fParticlesMatched->Fill(contents);
1120 }
Int_t pdg
Bool_t fDoEmcalQA
Set whether to enable cell/cluster QA.
Double_t * fEtaHistBins
number of eta bins
Double_t * fPtResHistBins
number of pt res bins
TObjArray fClusterCollArray
cluster collection array
Bool_t fDoEventQA
Set whether to enable event QA.
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Int_t fN1OverPtResHistBins
1/pt res bins
const char * title
Definition: MakeQAPdf.C:26
Float_t fPtBinWidth
Histogram pt bin width.
Int_t fNPtResHistBins
pt relative difference bins
AliTrackContainer * fDetectorLevel
generator level container
const AliMCParticleIterableMomentumContainer accepted_momentum() const
bidirectional stl iterator over the EMCAL iterable container
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
AliAnalysisTaskPWGJEQA()
Default constructor for ROOT I/O purposes.
THnSparse * fTracks
integer bins
THistManager fHistManager
primary particles matched to detector level tracks
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
Bool_t fDoTrackQA
Set whether to enable track QA.
THnSparse * fParticlesPhysPrim
all tracks
TObjArray fParticleCollArray
particle/track collection array
Double_t * sigma
const Int_t nPtBins
Double_t * fPtHistBins
number of pt bins
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
TObject * FindObject(const char *name) const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
float Float_t
Definition: External.C:68
AliEMCALGeometry * fGeom
!emcal geometry
void FillDetectorLevelTHnSparse(Double_t cent, Double_t trackEta, Double_t trackPhi, Double_t trackPt, Double_t sigma1OverPt, Byte_t trackType)
Bool_t fSeparateEMCalDCal
Separate EMCal from DCal in QA plots.
virtual AliAODMCParticle * GetAcceptMCParticleWithLabel(Int_t lab)
Int_t fNTotTracks
!Total number of accepted tracks in current event
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Float_t fCellEnergyCut
Energy cell cut.
BeamType fForceBeamType
forced beam type
virtual Bool_t AcceptCluster(Int_t i, UInt_t &rejectionReason) const
Double_t * fPtRelDiffHistBins
number of pt relative difference bins
Declaration of class AliAnalysisTaskPWGJEQA.
Double_t fCent
!event centrality
Int_t fNPtHistBins
detector level container
TString fCaloCellsName
name of calo cell collection
Int_t fNIntegerHistBins
number of 1/pt res bins
Double_t * f1OverPtResHistBins
pt res bins
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
TObjArray fJetCollArray
jet collection array
short Short_t
Definition: External.C:23
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
Double_t * fIntegerHistBins
number of integer bins
const AliClusterIterableMomentumContainer all_momentum() const
void FillMatchedParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
void FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt)
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
AliTLorentzVector fLeadingTrack
!Leading track in current event
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
AliMCParticleContainer * GetMCParticleContainer(Int_t i=0) const
Char_t GetTrackType(const AliVTrack *track) const
Definition: External.C:220
Bool_t fIsEsd
!whether it's an ESD analysis
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
void GenerateHistoBins()
Generate histogram binning arrays for track histograms.
AliTrackContainer * GetTrackContainer(Int_t i=0) const
void SetMakeGeneralHistograms(Bool_t g)
This is a task used to do basic PWGJE QA on tracks, clusters, and jets. Based on code from Salvatore ...
Base task in the EMCAL jet framework.
const AliTrackIterableMomentumContainer accepted_momentum() const
Bool_t fDoJetQA
Set whether to enable jet QA.
TString fGeneratorLevelName
generator level container name
void SetRejectionReasonLabels(TAxis *axis)
Float_t fMaxPt
Histogram pt limit.
AliMCParticleContainer * fGeneratorLevel
const Int_t nbins
bool Bool_t
Definition: External.C:53
Double_t * fCentHistBins
number of cent bins
THnSparse * fParticlesMatched
all physical primary particles
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Container structure for EMCAL clusters.
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Container for jet within the EMCAL jet framework.
Double_t * fPhiHistBins
number of phi bins
AliTLorentzVector fLeadingCluster[3]
!Leading cluster in current event (EMCal/DCal)
const AliJetIterableContainer all() const
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal