AliPhysics  c923f52 (c923f52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliEmcalTrackingQATask.cxx
Go to the documentation of this file.
1 // Track QA task (efficiency and pt resolution)
2 //
3 // Author: S.Aiola
4 
5 #include <TH3F.h>
6 #include <THnSparse.h>
7 #include <TMath.h>
8 #include <TString.h>
9 #include <Riostream.h>
10 
11 #include "AliPicoTrack.h"
12 #include "AliESDtrack.h"
13 #include "AliAODMCParticle.h"
14 #include "AliMCParticleContainer.h"
15 #include "AliLog.h"
16 
17 #include "AliEmcalTrackingQATask.h"
18 
20 
21 //________________________________________________________________________
23  AliAnalysisTaskEmcal("AliEmcalTrackingQA", kTRUE),
24  fDoSigma1OverPt(kFALSE),
25  fDoSigmaPtOverPtGen(kFALSE),
26  fGeneratorLevel(0),
27  fDetectorLevel(0),
28  fNPtHistBins(0),
29  fPtHistBins(0),
30  fNEtaHistBins(0),
31  fEtaHistBins(0),
32  fNPhiHistBins(0),
33  fPhiHistBins(0),
34  fNCentHistBins(0),
35  fCentHistBins(0),
36  fNPtRelDiffHistBins(0),
37  fPtRelDiffHistBins(0),
38  fNPtResHistBins(0),
39  fPtResHistBins(0),
40  f1OverPtResHistBins(0),
41  fN1OverPtResHistBins(0),
42  fNIntegerHistBins(0),
43  fIntegerHistBins(0),
44  fTracks(0),
45  fParticlesPhysPrim(0),
46  fParticlesMatched(0)
47 {
48  // Default constructor.
49 
50  SetMakeGeneralHistograms(kTRUE);
51 
52  GenerateHistoBins();
53 }
54 
55 //________________________________________________________________________
57  AliAnalysisTaskEmcal(name, kTRUE),
58  fDoSigma1OverPt(kFALSE),
59  fDoSigmaPtOverPtGen(kFALSE),
60  fGeneratorLevel(0),
61  fDetectorLevel(0),
62  fNPtHistBins(0),
63  fPtHistBins(0),
64  fNEtaHistBins(0),
65  fEtaHistBins(0),
66  fNPhiHistBins(0),
67  fPhiHistBins(0),
68  fNCentHistBins(0),
69  fCentHistBins(0),
70  fNPtRelDiffHistBins(0),
71  fPtRelDiffHistBins(0),
72  fNPtResHistBins(0),
73  fPtResHistBins(0),
74  f1OverPtResHistBins(0),
75  fN1OverPtResHistBins(0),
76  fNIntegerHistBins(0),
77  fIntegerHistBins(0),
78  fTracks(0),
79  fParticlesPhysPrim(0),
80  fParticlesMatched(0)
81 {
82  // Standard constructor.
83 
85 
87 }
88 
89 //________________________________________________________________________
91 {
92  // Destructor.
93 }
94 
95 //________________________________________________________________________
97 {
98  fNPtHistBins = 82;
101  GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
102  GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
103  GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
104  GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
105  GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
106  GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
107 
108  fNEtaHistBins = 100;
111 
112  fNPhiHistBins = 101;
114  GenerateFixedBinArray(fNPhiHistBins, 0, TMath::Pi() * 2.02, fPhiHistBins);
115 
116  fNCentHistBins = 4;
118  fCentHistBins[0] = 0;
119  fCentHistBins[1] = 10;
120  fCentHistBins[2] = 30;
121  fCentHistBins[3] = 50;
122  fCentHistBins[4] = 90;
123 
124  fNPtResHistBins = 175;
127  GenerateFixedBinArray(25, 0.05, 0.10, fPtResHistBins+50);
128  GenerateFixedBinArray(25, 0.10, 0.20, fPtResHistBins+75);
129  GenerateFixedBinArray(30, 0.20, 0.50, fPtResHistBins+100);
130  GenerateFixedBinArray(25, 0.50, 1.00, fPtResHistBins+130);
131  GenerateFixedBinArray(20, 1.00, 2.00, fPtResHistBins+155);
132 
133  fNPtRelDiffHistBins = 200;
136 
137  fN1OverPtResHistBins = 385;
140  GenerateFixedBinArray(60, 0.02, 0.05, f1OverPtResHistBins+100);
141  GenerateFixedBinArray(50, 0.05, 0.1, f1OverPtResHistBins+160);
142  GenerateFixedBinArray(50, 0.1, 0.2, f1OverPtResHistBins+210);
143  GenerateFixedBinArray(75, 0.2, 0.5, f1OverPtResHistBins+260);
144  GenerateFixedBinArray(50, 0.5, 1.5, f1OverPtResHistBins+335);
145 
146  fNIntegerHistBins = 10;
149 }
150 
151 //________________________________________________________________________
153 {
154  // Create my user objects.
155 
157 
158  if (fParticleCollArray.GetEntriesFast() < 1) {
159  AliFatal("This task needs at least one particle container!");
160  }
161 
162  if (!fDetectorLevel) {
163  fDetectorLevel = static_cast<AliTrackContainer*>(fParticleCollArray.At(0));
164  }
165 
166  if (!fGeneratorLevel && fParticleCollArray.GetEntriesFast() > 1) {
168  }
169 
171 
172  if (fGeneratorLevel) {
175  }
176 }
177 
178 //________________________________________________________________________
180 {
181  Int_t dim = 0;
182  TString title[20];
183  Int_t nbins[20] = {0};
184  Double_t *binEdges[20] = {0};
185 
187  title[dim] = "Centrality %";
188  nbins[dim] = fNCentHistBins;
189  binEdges[dim] = fCentHistBins;
190  dim++;
191  }
192 
193  title[dim] = "#it{p}_{T} (GeV/#it{c})";
194  nbins[dim] = fNPtHistBins;
195  binEdges[dim] = fPtHistBins;
196  dim++;
197 
198  title[dim] = "#eta";
199  nbins[dim] = fNEtaHistBins;
200  binEdges[dim] = fEtaHistBins;
201  dim++;
202 
203  title[dim] = "#phi";
204  nbins[dim] = fNPhiHistBins;
205  binEdges[dim] = fPhiHistBins;
206  dim++;
207 
208  title[dim] = "MC Generator";
209  nbins[dim] = 2;
210  binEdges[dim] = fIntegerHistBins;
211  dim++;
212 
213  title[dim] = "track type";
214  nbins[dim] = 3;
215  binEdges[dim] = fIntegerHistBins;
216  dim++;
217 
218  if (fIsEsd) {
219  if (fDoSigma1OverPt) {
220  title[dim] = "#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}";
221  nbins[dim] = fN1OverPtResHistBins;
222  binEdges[dim] = f1OverPtResHistBins;
223  dim++;
224  }
225  else {
226  title[dim] = "#sigma(#it{p}_{T}) / #it{p}_{T}";
227  nbins[dim] = fNPtResHistBins;
228  binEdges[dim] = fPtResHistBins;
229  dim++;
230  }
231  }
232 
233  fTracks = new THnSparseF("fTracks","fTracks",dim,nbins);
234  for (Int_t i = 0; i < dim; i++) {
235  fTracks->GetAxis(i)->SetTitle(title[i]);
236  fTracks->SetBinEdges(i, binEdges[i]);
237  }
238 
239  fOutput->Add(fTracks);
240 }
241 
242 //________________________________________________________________________
244 {
245  Int_t dim = 0;
246  TString title[20];
247  Int_t nbins[20] = {0};
248  Double_t *binEdges[20] = {0};
249 
251  title[dim] = "Centrality %";
252  nbins[dim] = fNCentHistBins;
253  binEdges[dim] = fCentHistBins;
254  dim++;
255  }
256 
257  title[dim] = "#it{p}_{T} (GeV/#it{c})";
258  nbins[dim] = fNPtHistBins;
259  binEdges[dim] = fPtHistBins;
260  dim++;
261 
262  title[dim] = "#eta";
263  nbins[dim] = fNEtaHistBins;
264  binEdges[dim] = fEtaHistBins;
265  dim++;
266 
267  title[dim] = "#phi";
268  nbins[dim] = fNPhiHistBins;
269  binEdges[dim] = fPhiHistBins;
270  dim++;
271 
272  title[dim] = "MC Generator";
273  nbins[dim] = 2;
274  binEdges[dim] = fIntegerHistBins;
275  dim++;
276 
277  title[dim] = "Findable";
278  nbins[dim] = 2;
279  binEdges[dim] = fIntegerHistBins;
280  dim++;
281 
282  fParticlesPhysPrim = new THnSparseF("fParticlesPhysPrim","fParticlesPhysPrim",dim,nbins);
283  for (Int_t i = 0; i < dim; i++) {
284  fParticlesPhysPrim->GetAxis(i)->SetTitle(title[i]);
285  fParticlesPhysPrim->SetBinEdges(i, binEdges[i]);
286  }
287 
289 }
290 
291 //________________________________________________________________________
293 {
294  Int_t dim = 0;
295  TString title[20];
296  Int_t nbins[20] = {0};
297  Double_t *binEdges[20] = {0};
298 
300  title[dim] = "Centrality %";
301  nbins[dim] = fNCentHistBins;
302  binEdges[dim] = fCentHistBins;
303  dim++;
304  }
305 
306  title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
307  nbins[dim] = fNPtHistBins;
308  binEdges[dim] = fPtHistBins;
309  dim++;
310 
311  title[dim] = "#eta^{gen}";
312  nbins[dim] = fNEtaHistBins;
313  binEdges[dim] = fEtaHistBins;
314  dim++;
315 
316  title[dim] = "#phi^{gen}";
317  nbins[dim] = fNPhiHistBins;
318  binEdges[dim] = fPhiHistBins;
319  dim++;
320 
321  title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
322  nbins[dim] = fNPtHistBins;
323  binEdges[dim] = fPtHistBins;
324  dim++;
325 
326  title[dim] = "#eta^{det}";
327  nbins[dim] = fNEtaHistBins;
328  binEdges[dim] = fEtaHistBins;
329  dim++;
330 
331  title[dim] = "#phi^{det}";
332  nbins[dim] = fNPhiHistBins;
333  binEdges[dim] = fPhiHistBins;
334  dim++;
335 
336  if (fDoSigmaPtOverPtGen) {
337  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
338  nbins[dim] = fNPtRelDiffHistBins;
339  binEdges[dim] = fPtRelDiffHistBins;
340  dim++;
341  }
342  else {
343  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}";
344  nbins[dim] = fNPtRelDiffHistBins;
345  binEdges[dim] = fPtRelDiffHistBins;
346  dim++;
347  }
348 
349  title[dim] = "track type";
350  nbins[dim] = 3;
351  binEdges[dim] = fIntegerHistBins;
352  dim++;
353 
354  fParticlesMatched = new THnSparseF("fParticlesMatched","fParticlesMatched",dim,nbins);
355  for (Int_t i = 0; i < dim; i++) {
356  fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
357  fParticlesMatched->SetBinEdges(i, binEdges[i]);
358  }
359 
361 }
362 
363 //________________________________________________________________________
365 {
366  if (!fDetectorLevel) {
367  AliError("Please, first set the detector level array!");
368  return;
369  }
370  if (!fGeneratorLevel) { // first check if the generator level array is set
372  if (fGeneratorLevel) { // now check if the first collection array has been added already
373  fGeneratorLevel->SetArrayName(name);
374  }
375  else {
377  }
380  }
381  fGeneratorLevel->SetArrayName(name);
382 }
383 
384 //________________________________________________________________________
386 {
387  if (!fDetectorLevel) { // first check if the detector level array is set
388  fDetectorLevel = static_cast<AliTrackContainer*>(fParticleCollArray.At(0));
389  if (fDetectorLevel) { // now check if the second collection array has been added already
390  fDetectorLevel->SetArrayName(name);
391  }
392  else {
394  }
395  }
396  fDetectorLevel->SetArrayName(name);
397 }
398 
399 //________________________________________________________________________
401 {
402  // Init the analysis.
403 
405 }
406 
407 //________________________________________________________________________
409  Double_t sigma1OverPt, Int_t mcGen, Byte_t trackType)
410 {
411  Double_t contents[20]={0};
412 
413  for (Int_t i = 0; i < fTracks->GetNdimensions(); i++) {
414  TString title(fTracks->GetAxis(i)->GetTitle());
415  if (title=="Centrality %")
416  contents[i] = cent;
417  else if (title=="#it{p}_{T} (GeV/#it{c})")
418  contents[i] = trackPt;
419  else if (title=="#eta")
420  contents[i] = trackEta;
421  else if (title=="#phi")
422  contents[i] = trackPhi;
423  else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
424  contents[i] = sigma1OverPt;
425  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
426  contents[i] = sigma1OverPt*trackPt;
427  else if (title=="MC Generator")
428  contents[i] = mcGen;
429  else if (title=="track type")
430  contents[i] = trackType;
431  else
432  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fTracks->GetName()));
433  }
434 
435  fTracks->Fill(contents);
436 }
437 
438 //________________________________________________________________________
439 void AliEmcalTrackingQATask::FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Int_t mcGen, Byte_t findable)
440 {
441  Double_t contents[20]={0};
442 
443  for (Int_t i = 0; i < fParticlesPhysPrim->GetNdimensions(); i++) {
444  TString title(fParticlesPhysPrim->GetAxis(i)->GetTitle());
445  if (title=="Centrality %")
446  contents[i] = cent;
447  else if (title=="#it{p}_{T} (GeV/#it{c})")
448  contents[i] = partPt;
449  else if (title=="#eta")
450  contents[i] = partEta;
451  else if (title=="#phi")
452  contents[i] = partPhi;
453  else if (title=="MC Generator")
454  contents[i] = mcGen;
455  else if (title=="Findable")
456  contents[i] = findable;
457  else
458  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesPhysPrim->GetName()));
459  }
460 
461  fParticlesPhysPrim->Fill(contents);
462 }
463 
464 //________________________________________________________________________
466  Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
467 {
468  Double_t contents[20]={0};
469 
470  for (Int_t i = 0; i < fParticlesMatched->GetNdimensions(); i++) {
471  TString title(fParticlesMatched->GetAxis(i)->GetTitle());
472  if (title=="Centrality %")
473  contents[i] = cent;
474  else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
475  contents[i] = partPt;
476  else if (title=="#eta^{gen}")
477  contents[i] = partEta;
478  else if (title=="#phi^{gen}")
479  contents[i] = partPhi;
480  else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
481  contents[i] = trackPt;
482  else if (title=="#eta^{det}")
483  contents[i] = trackEta;
484  else if (title=="#phi^{det}")
485  contents[i] = trackPhi;
486  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
487  contents[i] = (partPt - trackPt) / partPt;
488  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}")
489  contents[i] = (partPt - trackPt) / trackPt;
490  else if (title=="track type")
491  contents[i] = (Double_t)trackType;
492  else
493  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesMatched->GetName()));
494  }
495 
496  fParticlesMatched->Fill(contents);
497 }
498 
499 //________________________________________________________________________
501 {
502  // Fill the histograms.
503 
504  fDetectorLevel->ResetCurrentID();
505  AliVTrack *track = fDetectorLevel->GetNextAcceptTrack();
506  while (track != 0) {
507  Byte_t type = fDetectorLevel->GetTrackType(fDetectorLevel->GetCurrentID());
508  if (type <= 2) {
509  Double_t sigma = 0;
510  if (fIsEsd) {
511  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
512  if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
513  }
514 
515  Int_t label = TMath::Abs(track->GetLabel());
516  Int_t mcGen = 1;
517  // reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
518  if (label==0 || track->GetGeneratorIndex() == 0) mcGen = 0;
519 
520  FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, mcGen, type);
521 
522  if (fGeneratorLevel && label > 0) {
523  AliAODMCParticle *part = fGeneratorLevel->GetAcceptMCParticleWithLabel(label);
524  if (part) {
525  if (part->GetGeneratorIndex() == 0) {
526  Int_t pdg = TMath::Abs(part->PdgCode());
527  // select charged pions, protons, kaons , electrons, muons
528  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
529  FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
530  }
531  }
532  }
533  }
534  }
535  else {
536  AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
537  }
538 
540  }
541 
542  if (fGeneratorLevel) {
543  fGeneratorLevel->ResetCurrentID();
544  AliAODMCParticle *part = fGeneratorLevel->GetNextAcceptMCParticle();
545  while (part != 0) {
546  Int_t mcGen = 1;
547  Byte_t findable = 0;
548 
549  if (part->GetGeneratorIndex() == 0) mcGen = 0;
550 
551  Int_t pdg = TMath::Abs(part->PdgCode());
552  // select charged pions, protons, kaons , electrons, muons
553  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
554 
555  FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), mcGen, findable);
557  }
558  }
559 
560  return kTRUE;
561 }
Double_t * fPtRelDiffHistBins
number of pt relative difference bins
Int_t pdg
void SetGeneratorLevelName(const char *name)
void SetParticlePtCut(Double_t cut)
double Double_t
Definition: External.C:58
AliTrackContainer * fDetectorLevel
generator level container
const char * title
Definition: MakeQAPdf.C:26
Base task in the EMCAL framework.
Container with name, TClonesArray and cuts for particles.
void SetDetectorLevelName(const char *name)
Int_t fNPtResHistBins
pt relative difference bins
Int_t fNPtHistBins
detector level container
ClassImp(AliEmcalTrackingQATask) AliEmcalTrackingQATask
Double_t * fCentHistBins
number of cent bins
TObjArray fParticleCollArray
particle/track collection array
Double_t * fPtResHistBins
number of pt res bins
Double_t * fEtaHistBins
number of eta bins
Double_t * sigma
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 FillDetectorLevelTHnSparse(Double_t cent, Double_t trackEta, Double_t trackPhi, Double_t trackPt, Double_t sigma1OverPt, Int_t mcGen, Byte_t trackType)
int Int_t
Definition: External.C:63
virtual AliAODMCParticle * GetAcceptMCParticleWithLabel(Int_t lab)
BeamType fForceBeamType
forced beam type
virtual AliAODMCParticle * GetNextAcceptMCParticle()
Double_t fCent
!event centrality
Double_t * f1OverPtResHistBins
pt res bins
THnSparse * fTracks
integer bins
AliMCParticleContainer * AddMCParticleContainer(const char *n)
static Double_t * GenerateFixedBinArray(Int_t n, Double_t min, Double_t max)
void FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Int_t mcGen, Byte_t findable)
virtual AliVTrack * GetNextAcceptTrack()
AliEmcalList * fOutput
!output list
Int_t fNIntegerHistBins
number of 1/pt res bins
AliMCParticleContainer * fGeneratorLevel
Char_t GetTrackType(const AliVTrack *track) const
Double_t * fPtHistBins
number of pt bins
Bool_t fIsEsd
!whether it's an ESD analysis
AliTrackContainer * AddTrackContainer(const char *n)
Double_t * fIntegerHistBins
number of integer bins
void SelectPhysicalPrimaries(Bool_t s)
Int_t fNPtRelDiffHistBins
cent bins
void SetMakeGeneralHistograms(Bool_t g)
Double_t * fPhiHistBins
number of phi bins
const Int_t nbins
bool Bool_t
Definition: External.C:53
THnSparse * fParticlesPhysPrim
all tracks
Container for MC-true particles within the EMCAL framework.
Int_t fN1OverPtResHistBins
1/pt res bins
THnSparse * fParticlesMatched
all physical primary particles