AliPhysics  cda3415 (cda3415)
 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 "AliESDtrack.h"
12 #include "AliAODTrack.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 (fDoSigma1OverPt) {
219  title[dim] = "#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}";
220  nbins[dim] = fN1OverPtResHistBins;
221  binEdges[dim] = f1OverPtResHistBins;
222  dim++;
223  }
224  else {
225  title[dim] = "#sigma(#it{p}_{T}) / #it{p}_{T}";
226  nbins[dim] = fNPtResHistBins;
227  binEdges[dim] = fPtResHistBins;
228  dim++;
229  }
230 
231  fTracks = new THnSparseF("fTracks","fTracks",dim,nbins);
232  for (Int_t i = 0; i < dim; i++) {
233  fTracks->GetAxis(i)->SetTitle(title[i]);
234  fTracks->SetBinEdges(i, binEdges[i]);
235  }
236 
237  fOutput->Add(fTracks);
238 }
239 
240 //________________________________________________________________________
242 {
243  Int_t dim = 0;
244  TString title[20];
245  Int_t nbins[20] = {0};
246  Double_t *binEdges[20] = {0};
247 
249  title[dim] = "Centrality %";
250  nbins[dim] = fNCentHistBins;
251  binEdges[dim] = fCentHistBins;
252  dim++;
253  }
254 
255  title[dim] = "#it{p}_{T} (GeV/#it{c})";
256  nbins[dim] = fNPtHistBins;
257  binEdges[dim] = fPtHistBins;
258  dim++;
259 
260  title[dim] = "#eta";
261  nbins[dim] = fNEtaHistBins;
262  binEdges[dim] = fEtaHistBins;
263  dim++;
264 
265  title[dim] = "#phi";
266  nbins[dim] = fNPhiHistBins;
267  binEdges[dim] = fPhiHistBins;
268  dim++;
269 
270  title[dim] = "MC Generator";
271  nbins[dim] = 2;
272  binEdges[dim] = fIntegerHistBins;
273  dim++;
274 
275  title[dim] = "Findable";
276  nbins[dim] = 2;
277  binEdges[dim] = fIntegerHistBins;
278  dim++;
279 
280  fParticlesPhysPrim = new THnSparseF("fParticlesPhysPrim","fParticlesPhysPrim",dim,nbins);
281  for (Int_t i = 0; i < dim; i++) {
282  fParticlesPhysPrim->GetAxis(i)->SetTitle(title[i]);
283  fParticlesPhysPrim->SetBinEdges(i, binEdges[i]);
284  }
285 
287 }
288 
289 //________________________________________________________________________
291 {
292  Int_t dim = 0;
293  TString title[20];
294  Int_t nbins[20] = {0};
295  Double_t *binEdges[20] = {0};
296 
298  title[dim] = "Centrality %";
299  nbins[dim] = fNCentHistBins;
300  binEdges[dim] = fCentHistBins;
301  dim++;
302  }
303 
304  title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
305  nbins[dim] = fNPtHistBins;
306  binEdges[dim] = fPtHistBins;
307  dim++;
308 
309  title[dim] = "#eta^{gen}";
310  nbins[dim] = fNEtaHistBins;
311  binEdges[dim] = fEtaHistBins;
312  dim++;
313 
314  title[dim] = "#phi^{gen}";
315  nbins[dim] = fNPhiHistBins;
316  binEdges[dim] = fPhiHistBins;
317  dim++;
318 
319  title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
320  nbins[dim] = fNPtHistBins;
321  binEdges[dim] = fPtHistBins;
322  dim++;
323 
324  title[dim] = "#eta^{det}";
325  nbins[dim] = fNEtaHistBins;
326  binEdges[dim] = fEtaHistBins;
327  dim++;
328 
329  title[dim] = "#phi^{det}";
330  nbins[dim] = fNPhiHistBins;
331  binEdges[dim] = fPhiHistBins;
332  dim++;
333 
334  if (fDoSigmaPtOverPtGen) {
335  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
336  nbins[dim] = fNPtRelDiffHistBins;
337  binEdges[dim] = fPtRelDiffHistBins;
338  dim++;
339  }
340  else {
341  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}";
342  nbins[dim] = fNPtRelDiffHistBins;
343  binEdges[dim] = fPtRelDiffHistBins;
344  dim++;
345  }
346 
347  title[dim] = "track type";
348  nbins[dim] = 3;
349  binEdges[dim] = fIntegerHistBins;
350  dim++;
351 
352  fParticlesMatched = new THnSparseF("fParticlesMatched","fParticlesMatched",dim,nbins);
353  for (Int_t i = 0; i < dim; i++) {
354  fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
355  fParticlesMatched->SetBinEdges(i, binEdges[i]);
356  }
357 
359 }
360 
361 //________________________________________________________________________
363 {
364  if (!fDetectorLevel) {
365  AliError("Please, first set the detector level array!");
366  return;
367  }
368  if (!fGeneratorLevel) { // first check if the generator level array is set
370  if (fGeneratorLevel) { // now check if the first collection array has been added already
371  fGeneratorLevel->SetArrayName(name);
372  }
373  else {
375  }
378  }
379  fGeneratorLevel->SetArrayName(name);
380 }
381 
382 //________________________________________________________________________
384 {
385  if (!fDetectorLevel) { // first check if the detector level array is set
386  fDetectorLevel = static_cast<AliTrackContainer*>(fParticleCollArray.At(0));
387  if (fDetectorLevel) { // now check if the second collection array has been added already
388  fDetectorLevel->SetArrayName(name);
389  }
390  else {
392  }
393  }
394  fDetectorLevel->SetArrayName(name);
395 }
396 
397 //________________________________________________________________________
399 {
400  // Init the analysis.
401 
403 }
404 
405 //________________________________________________________________________
407  Double_t sigma1OverPt, Int_t mcGen, Byte_t trackType)
408 {
409  Double_t contents[20]={0};
410 
411  for (Int_t i = 0; i < fTracks->GetNdimensions(); i++) {
412  TString title(fTracks->GetAxis(i)->GetTitle());
413  if (title=="Centrality %")
414  contents[i] = cent;
415  else if (title=="#it{p}_{T} (GeV/#it{c})")
416  contents[i] = trackPt;
417  else if (title=="#eta")
418  contents[i] = trackEta;
419  else if (title=="#phi")
420  contents[i] = trackPhi;
421  else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
422  contents[i] = sigma1OverPt;
423  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
424  contents[i] = sigma1OverPt*trackPt;
425  else if (title=="MC Generator")
426  contents[i] = mcGen;
427  else if (title=="track type")
428  contents[i] = trackType;
429  else
430  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fTracks->GetName()));
431  }
432 
433  fTracks->Fill(contents);
434 }
435 
436 //________________________________________________________________________
437 void AliEmcalTrackingQATask::FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Int_t mcGen, Byte_t findable)
438 {
439  Double_t contents[20]={0};
440 
441  for (Int_t i = 0; i < fParticlesPhysPrim->GetNdimensions(); i++) {
442  TString title(fParticlesPhysPrim->GetAxis(i)->GetTitle());
443  if (title=="Centrality %")
444  contents[i] = cent;
445  else if (title=="#it{p}_{T} (GeV/#it{c})")
446  contents[i] = partPt;
447  else if (title=="#eta")
448  contents[i] = partEta;
449  else if (title=="#phi")
450  contents[i] = partPhi;
451  else if (title=="MC Generator")
452  contents[i] = mcGen;
453  else if (title=="Findable")
454  contents[i] = findable;
455  else
456  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesPhysPrim->GetName()));
457  }
458 
459  fParticlesPhysPrim->Fill(contents);
460 }
461 
462 //________________________________________________________________________
464  Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
465 {
466  Double_t contents[20]={0};
467 
468  for (Int_t i = 0; i < fParticlesMatched->GetNdimensions(); i++) {
469  TString title(fParticlesMatched->GetAxis(i)->GetTitle());
470  if (title=="Centrality %")
471  contents[i] = cent;
472  else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
473  contents[i] = partPt;
474  else if (title=="#eta^{gen}")
475  contents[i] = partEta;
476  else if (title=="#phi^{gen}")
477  contents[i] = partPhi;
478  else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
479  contents[i] = trackPt;
480  else if (title=="#eta^{det}")
481  contents[i] = trackEta;
482  else if (title=="#phi^{det}")
483  contents[i] = trackPhi;
484  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
485  contents[i] = (partPt - trackPt) / partPt;
486  else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{det}")
487  contents[i] = (partPt - trackPt) / trackPt;
488  else if (title=="track type")
489  contents[i] = (Double_t)trackType;
490  else
491  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesMatched->GetName()));
492  }
493 
494  fParticlesMatched->Fill(contents);
495 }
496 
497 //________________________________________________________________________
499 {
500  // Fill the histograms.
501 
502  fDetectorLevel->ResetCurrentID();
503  AliVTrack* track;
504  for (auto trackIterator : fDetectorLevel->accepted_momentum() ) {
505  track = trackIterator.second;
506  Byte_t type = fDetectorLevel->GetTrackType(track);
507  if (type <= 2) {
508  Double_t sigma = 0;
509 
510  if (fIsEsd) {
511 
512  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
513  if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
514 
515  } else { // AOD
516 
517  AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
518  if(!aodtrack) AliFatal("Not a standard AOD");
519 
520  AliExternalTrackParam exParam;
521 
522  //get covariance matrix
523  Double_t cov[21] = {0,};
524  aodtrack->GetCovMatrix(cov);
525  Double_t pxpypz[3] = {0,};
526  aodtrack->PxPyPz(pxpypz);
527  Double_t xyz[3] = {0,};
528  aodtrack->GetXYZ(xyz);
529  Short_t sign = aodtrack->Charge();
530  exParam.Set(xyz,pxpypz,cov,sign);
531 
532  sigma = TMath::Sqrt(exParam.GetSigma1Pt2());
533 
534  }
535 
536  Int_t label = TMath::Abs(track->GetLabel());
537  Int_t mcGen = 1;
538  // reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
539  if (label==0 || track->GetGeneratorIndex() == 0) mcGen = 0;
540 
541  FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, mcGen, type);
542 
543  if (fGeneratorLevel && label > 0) {
544  AliAODMCParticle *part = fGeneratorLevel->GetAcceptMCParticleWithLabel(label);
545  if (part) {
546  if (part->GetGeneratorIndex() == 0) {
547  Int_t pdg = TMath::Abs(part->PdgCode());
548  // select charged pions, protons, kaons , electrons, muons
549  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
550  FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
551  }
552  }
553  }
554  }
555  }
556  else {
557  AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
558  }
559  }
560 
561  if (fGeneratorLevel) {
562  fGeneratorLevel->ResetCurrentID();
563  AliAODMCParticle* part;
564  for (auto partIterator : fGeneratorLevel->accepted_momentum() ) {
565  part = partIterator.second;
566 
567  Int_t mcGen = 1;
568  Byte_t findable = 0;
569 
570  if (part->GetGeneratorIndex() == 0) mcGen = 0;
571 
572  Int_t pdg = TMath::Abs(part->PdgCode());
573  // select charged pions, protons, kaons , electrons, muons
574  if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
575 
576  FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), mcGen, findable);
577  }
578  }
579 
580  return kTRUE;
581 }
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
const AliMCParticleIterableMomentumContainer accepted_momentum() const
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
Double_t fCent
!event centrality
Double_t * f1OverPtResHistBins
pt res bins
THnSparse * fTracks
integer bins
AliMCParticleContainer * AddMCParticleContainer(const char *n)
short Short_t
Definition: External.C:23
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)
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)
const AliTrackIterableMomentumContainer accepted_momentum() const
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