PWG/EMCAL/AliEmcalTrackingQATask.cxx

   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 "AliParticleContainer.h"
  15 #include "AliLog.h"
  16
  17 #include "AliEmcalTrackingQATask.h"
  18
  19 ClassImp(AliEmcalTrackingQATask)
  20
  21 //________________________________________________________________________
  22 AliEmcalTrackingQATask::AliEmcalTrackingQATask() :
  23   AliAnalysisTaskEmcal("AliEmcalTrackingQA", kTRUE),
  24   fSelectHIJING(kTRUE),
  25   fGeneratorLevel(0),
  26   fDetectorLevel(0),
  27   fNPtHistBins(0),
  28   fPtHistBins(0),
  29   fNEtaHistBins(0),
  30   fEtaHistBins(0),
  31   fNPhiHistBins(0),
  32   fPhiHistBins(0),
  33   fNCentHistBins(0),
  34   fCentHistBins(0),
  35   fNPtResHistBins(0),
  36   fPtResHistBins(0),
  37   f1OverPtResHistBins(0),
  38   fN1OverPtResHistBins(0),
  39   fNIntegerHistBins(0),
  40   fIntegerHistBins(0),
  41   fTracks(0),
  42   fParticlesPhysPrim(0),
  43   fParticlesMatched(0)
  44 {
  45   // Default constructor.
  46
  47   SetMakeGeneralHistograms(kTRUE);
  48
  49   GenerateHistoBins();
  50 }
  51
  52 //________________________________________________________________________
  53 AliEmcalTrackingQATask::AliEmcalTrackingQATask(const char *name) :
  54   AliAnalysisTaskEmcal(name, kTRUE),
  55   fSelectHIJING(kTRUE),
  56   fGeneratorLevel(0),
  57   fDetectorLevel(0),
  58   fNPtHistBins(0),
  59   fPtHistBins(0),
  60   fNEtaHistBins(0),
  61   fEtaHistBins(0),
  62   fNPhiHistBins(0),
  63   fPhiHistBins(0),
  64   fNCentHistBins(0),
  65   fCentHistBins(0),
  66   fNPtResHistBins(0),
  67   fPtResHistBins(0),
  68   f1OverPtResHistBins(0),
  69   fN1OverPtResHistBins(0),
  70   fNIntegerHistBins(0),
  71   fIntegerHistBins(0),
  72   fTracks(0),
  73   fParticlesPhysPrim(0),
  74   fParticlesMatched(0)
  75 {
  76   // Standard constructor.
  77
  78   SetMakeGeneralHistograms(kTRUE);
  79
  80   GenerateHistoBins();
  81 }
  82
  83 //________________________________________________________________________
  84 AliEmcalTrackingQATask::~AliEmcalTrackingQATask()
  85 {
  86   // Destructor.
  87 }
  88
  89 //________________________________________________________________________
  90 void AliEmcalTrackingQATask::GenerateHistoBins()
  91 {
  92   fNPtHistBins = 82;
  93   fPtHistBins = new Double_t[fNPtHistBins+1];
  94   GenerateFixedBinArray(6, 0, 0.3, fPtHistBins);
  95   GenerateFixedBinArray(7, 0.3, 1, fPtHistBins+6);
  96   GenerateFixedBinArray(10, 1, 3, fPtHistBins+13);
  97   GenerateFixedBinArray(14, 3, 10, fPtHistBins+23);
  98   GenerateFixedBinArray(10, 10, 20, fPtHistBins+37);
  99   GenerateFixedBinArray(15, 20, 50, fPtHistBins+47);
 100   GenerateFixedBinArray(20, 50, 150, fPtHistBins+62);
 101
 102   fNEtaHistBins = 100;
 103   fEtaHistBins = new Double_t[fNEtaHistBins+1];
 104   GenerateFixedBinArray(fNEtaHistBins, -1, 1, fEtaHistBins);
 105
 106   fNPhiHistBins = 101;
 107   fPhiHistBins = new Double_t[fNPhiHistBins+1];
 108   GenerateFixedBinArray(fNPhiHistBins, 0, TMath::Pi() * 2.02, fPhiHistBins);
 109
 110   fNCentHistBins = 4;
 111   fCentHistBins = new Double_t[fNCentHistBins+1];
 112   fCentHistBins[0] = 0;
 113   fCentHistBins[1] = 10;
 114   fCentHistBins[2] = 30;
 115   fCentHistBins[3] = 50;
 116   fCentHistBins[4] = 90;
 117
 118   fNPtResHistBins = 200;
 119   fPtResHistBins = new Double_t[fNPtResHistBins+1];
 120   GenerateFixedBinArray(fNPtResHistBins, -2, 2, fPtResHistBins);
 121
 122   fN1OverPtResHistBins = 325;
 123   f1OverPtResHistBins = new Double_t[fN1OverPtResHistBins+1];
 124   GenerateFixedBinArray(100, 0, 0.05, f1OverPtResHistBins);
 125   GenerateFixedBinArray(50, 0.05, 0.1, f1OverPtResHistBins+100);
 126   GenerateFixedBinArray(50, 0.1, 0.2, f1OverPtResHistBins+50);
 127   GenerateFixedBinArray(75, 0.2, 0.5, f1OverPtResHistBins+50);
 128   GenerateFixedBinArray(50, 0.5, 1.5, f1OverPtResHistBins+75);
 129
 130   fNIntegerHistBins = 10;
 131   fIntegerHistBins = new Double_t[fNIntegerHistBins+1];
 132   GenerateFixedBinArray(fNIntegerHistBins, -0.5, 9.5, fIntegerHistBins);
 133 }
 134
 135 //________________________________________________________________________
 136 void AliEmcalTrackingQATask::UserCreateOutputObjects()
 137 {
 138   // Create my user objects.
 139
 140   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 141
 142   if (fParticleCollArray.GetEntriesFast() < 1) {
 143     AliFatal("This task needs at least one particle container!");
 144   }
 145
 146   if (!fDetectorLevel) {
 147     fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 148     fDetectorLevel->SetClassName("AliPicoTrack");
 149   }
 150
 151   if (!fGeneratorLevel && fParticleCollArray.GetEntriesFast() > 1) {
 152     fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 153     fGeneratorLevel->SetClassName("AliAODMCParticle");
 154   }
 155
 156   AllocateDetectorLevelTHnSparse();
 157
 158   if (fGeneratorLevel) {
 159     AllocateGeneratorLevelTHnSparse();
 160     AllocateMatchedParticlesTHnSparse();
 161   }
 162 }
 163
 164 //________________________________________________________________________
 165 void AliEmcalTrackingQATask::AllocateDetectorLevelTHnSparse()
 166 {
 167   Int_t dim = 0;
 168   TString title[20];
 169   Int_t nbins[20] = {0};
 170   Double_t *binEdges[20] = {0};
 171
 172   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 173     title[dim] = "Centrality %";
 174     nbins[dim] = fNCentHistBins;
 175     binEdges[dim] = fCentHistBins;
 176     dim++;
 177   }
 178
 179   title[dim] = "#it{p}_{T} (GeV/#it{c})";
 180   nbins[dim] = fNPtHistBins;
 181   binEdges[dim] = fPtHistBins;
 182   dim++;
 183
 184   title[dim] = "#eta";
 185   nbins[dim] = fNEtaHistBins;
 186   binEdges[dim] = fEtaHistBins;
 187   dim++;
 188
 189   title[dim] = "#phi";
 190   nbins[dim] = fNPhiHistBins;
 191   binEdges[dim] = fPhiHistBins;
 192   dim++;
 193
 194   if (fSelectHIJING) {
 195     title[dim] = "MC Generator";
 196     nbins[dim] = 2;
 197     binEdges[dim] = fIntegerHistBins;
 198     dim++;
 199   }
 200
 201   title[dim] = "track type";
 202   nbins[dim] = 3;
 203   binEdges[dim] = fIntegerHistBins;
 204   dim++;
 205
 206   if (fIsEsd) {
 207     title[dim] = "#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}";
 208     nbins[dim] = fN1OverPtResHistBins;
 209     binEdges[dim] = f1OverPtResHistBins;
 210     dim++;
 211   }
 212
 213   fTracks = new THnSparseF("fTracks","fTracks",dim,nbins);
 214   for (Int_t i = 0; i < dim; i++) {
 215     fTracks->GetAxis(i)->SetTitle(title[i]);
 216     fTracks->SetBinEdges(i, binEdges[i]);
 217   }
 218
 219   fOutput->Add(fTracks);
 220 }
 221
 222 //________________________________________________________________________
 223 void AliEmcalTrackingQATask::AllocateGeneratorLevelTHnSparse()
 224 {
 225   Int_t dim = 0;
 226   TString title[20];
 227   Int_t nbins[20] = {0};
 228   Double_t *binEdges[20] = {0};
 229
 230   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 231     title[dim] = "Centrality %";
 232     nbins[dim] = fNCentHistBins;
 233     binEdges[dim] = fCentHistBins;
 234     dim++;
 235   }
 236
 237   title[dim] = "#it{p}_{T} (GeV/#it{c})";
 238   nbins[dim] = fNPtHistBins;
 239   binEdges[dim] = fPtHistBins;
 240   dim++;
 241
 242   title[dim] = "#eta";
 243   nbins[dim] = fNEtaHistBins;
 244   binEdges[dim] = fEtaHistBins;
 245   dim++;
 246
 247   title[dim] = "#phi";
 248   nbins[dim] = fNPhiHistBins;
 249   binEdges[dim] = fPhiHistBins;
 250   dim++;
 251
 252   if (fSelectHIJING) {
 253     title[dim] = "MC Generator";
 254     nbins[dim] = 2;
 255     binEdges[dim] = fIntegerHistBins;
 256     dim++;
 257   }
 258
 259   title[dim] = "Findable";
 260   nbins[dim] = 2;
 261   binEdges[dim] = fIntegerHistBins;
 262   dim++;
 263
 264   fParticlesPhysPrim = new THnSparseF("fParticlesPhysPrim","fParticlesPhysPrim",dim,nbins);
 265   for (Int_t i = 0; i < dim; i++) {
 266     fParticlesPhysPrim->GetAxis(i)->SetTitle(title[i]);
 267     fParticlesPhysPrim->SetBinEdges(i, binEdges[i]);
 268   }
 269
 270   fOutput->Add(fParticlesPhysPrim);
 271 }
 272
 273 //________________________________________________________________________
 274 void AliEmcalTrackingQATask::AllocateMatchedParticlesTHnSparse()
 275 {
 276   Int_t dim = 0;
 277   TString title[20];
 278   Int_t nbins[20] = {0};
 279   Double_t *binEdges[20] = {0};
 280
 281   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 282     title[dim] = "Centrality %";
 283     nbins[dim] = fNCentHistBins;
 284     binEdges[dim] = fCentHistBins;
 285     dim++;
 286   }
 287
 288   title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
 289   nbins[dim] = fNPtHistBins;
 290   binEdges[dim] = fPtHistBins;
 291   dim++;
 292
 293   title[dim] = "#eta^{gen}";
 294   nbins[dim] = fNEtaHistBins;
 295   binEdges[dim] = fEtaHistBins;
 296   dim++;
 297
 298   title[dim] = "#phi^{gen}";
 299   nbins[dim] = fNPhiHistBins;
 300   binEdges[dim] = fPhiHistBins;
 301   dim++;
 302
 303   title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
 304   nbins[dim] = fNPtHistBins;
 305   binEdges[dim] = fPtHistBins;
 306   dim++;
 307
 308   title[dim] = "#eta^{det}";
 309   nbins[dim] = fNEtaHistBins;
 310   binEdges[dim] = fEtaHistBins;
 311   dim++;
 312
 313   title[dim] = "#phi^{det}";
 314   nbins[dim] = fNPhiHistBins;
 315   binEdges[dim] = fPhiHistBins;
 316   dim++;
 317
 318   title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
 319   nbins[dim] = fNPtResHistBins;
 320   binEdges[dim] = fPtResHistBins;
 321   dim++;
 322
 323   title[dim] = "track type";
 324   nbins[dim] = 3;
 325   binEdges[dim] = fIntegerHistBins;
 326   dim++;
 327
 328   fParticlesMatched = new THnSparseF("fParticlesMatched","fParticlesMatched",dim,nbins);
 329   for (Int_t i = 0; i < dim; i++) {
 330     fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
 331     fParticlesMatched->SetBinEdges(i, binEdges[i]);
 332   }
 333
 334   fOutput->Add(fParticlesMatched);
 335 }
 336
 337 //________________________________________________________________________
 338 void AliEmcalTrackingQATask::SetGeneratorLevelName(const char* name)
 339 {
 340   if (!fDetectorLevel) {
 341     AliError("Please, first set the detector level array!");
 342     return;
 343   }
 344   if (!fGeneratorLevel) {  // first check if the generator level array is set
 345     fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 346     if (fGeneratorLevel) {  // now check if the first collection array has been added already
 347       fGeneratorLevel->SetArrayName(name);
 348     }
 349     else {
 350       fGeneratorLevel = AddParticleContainer(name);
 351     }
 352     fGeneratorLevel->SetClassName("AliAODMCParticle");
 353     fGeneratorLevel->SelectPhysicalPrimaries(kTRUE);
 354     fGeneratorLevel->SetParticlePtCut(0);
 355   }
 356   fGeneratorLevel->SetArrayName(name);
 357 }
 358
 359 //________________________________________________________________________
 360 void AliEmcalTrackingQATask::SetDetectorLevelName(const char* name)
 361 {
 362   if (!fDetectorLevel) {  // first check if the detector level array is set
 363     fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 364     if (fDetectorLevel) {  // now check if the second collection array has been added already
 365       fDetectorLevel->SetArrayName(name);
 366     }
 367     else {
 368       fDetectorLevel = AddParticleContainer(name);
 369     }
 370     fDetectorLevel->SetClassName("AliPicoTrack");
 371   }
 372   fDetectorLevel->SetArrayName(name);
 373 }
 374
 375 //________________________________________________________________________
 376 void AliEmcalTrackingQATask::ExecOnce()
 377 {
 378   // Init the analysis.
 379
 380   AliAnalysisTaskEmcal::ExecOnce();
 381 }
 382
 383 //________________________________________________________________________
 384 void AliEmcalTrackingQATask::FillDetectorLevelTHnSparse(Double_t cent, Double_t trackEta, Double_t trackPhi, Double_t trackPt,
 385                                                         Double_t sigma1OverPt, Int_t mcGen, Byte_t trackType)
 386 {
 387   Double_t contents[20]={0};
 388
 389   for (Int_t i = 0; i < fTracks->GetNdimensions(); i++) {
 390     TString title(fTracks->GetAxis(i)->GetTitle());
 391     if (title=="Centrality %")
 392       contents[i] = cent;
 393     else if (title=="#it{p}_{T} (GeV/#it{c})")
 394       contents[i] = trackPt;
 395     else if (title=="#eta")
 396       contents[i] = trackEta;
 397     else if (title=="#phi")
 398       contents[i] = trackPhi;
 399     else if (title=="#sigma(1/#it{p}_{T}) (GeV/#it{c})^{-1}")
 400       contents[i] = sigma1OverPt;
 401     else if (title=="MC Generator")
 402       contents[i] = mcGen;
 403     else if (title=="track type")
 404       contents[i] = trackType;
 405     else
 406       AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fTracks->GetName()));
 407   }
 408
 409   fTracks->Fill(contents);
 410 }
 411
 412 //________________________________________________________________________
 413 void AliEmcalTrackingQATask::FillGeneratorLevelTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt, Int_t mcGen, Byte_t findable)
 414 {
 415   Double_t contents[20]={0};
 416
 417   for (Int_t i = 0; i < fParticlesPhysPrim->GetNdimensions(); i++) {
 418     TString title(fParticlesPhysPrim->GetAxis(i)->GetTitle());
 419     if (title=="Centrality %")
 420       contents[i] = cent;
 421     else if (title=="#it{p}_{T} (GeV/#it{c})")
 422       contents[i] = partPt;
 423     else if (title=="#eta")
 424       contents[i] = partEta;
 425     else if (title=="#phi")
 426       contents[i] = partPhi;
 427     else if (title=="MC Generator")
 428       contents[i] = mcGen;
 429     else if (title=="Findable")
 430       contents[i] = findable;
 431     else
 432       AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesPhysPrim->GetName()));
 433   }
 434
 435   fParticlesPhysPrim->Fill(contents);
 436 }
 437
 438 //________________________________________________________________________
 439 void AliEmcalTrackingQATask::FillMatchedParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt,
 440                                                            Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
 441 {
 442   Double_t contents[20]={0};
 443
 444   for (Int_t i = 0; i < fParticlesMatched->GetNdimensions(); i++) {
 445     TString title(fParticlesMatched->GetAxis(i)->GetTitle());
 446     if (title=="Centrality %")
 447       contents[i] = cent;
 448     else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
 449       contents[i] = partPt;
 450     else if (title=="#eta^{gen}")
 451       contents[i] = partEta;
 452     else if (title=="#phi^{gen}")
 453       contents[i] = partPhi;
 454     else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
 455       contents[i] = trackPt;
 456     else if (title=="#eta^{det}")
 457       contents[i] = trackEta;
 458     else if (title=="#phi^{det}")
 459       contents[i] = trackPhi;
 460     else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
 461       contents[i] = (partPt - trackPt) / partPt;
 462     else if (title=="track type")
 463       contents[i] = (Double_t)trackType;
 464     else
 465       AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesMatched->GetName()));
 466   }
 467
 468   fParticlesMatched->Fill(contents);
 469 }
 470
 471 //________________________________________________________________________
 472 Bool_t AliEmcalTrackingQATask::FillHistograms()
 473 {
 474   // Fill the histograms.
 475
 476   AliPicoTrack *track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle(0));
 477   while (track != 0) {
 478     Byte_t type = track->GetTrackType();
 479     if (type <= 2) {
 480       Double_t sigma = 0;
 481       if (fIsEsd) {
 482         AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track->GetTrack());
 483         if (esdTrack) sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
 484       }
 485
 486       Int_t label = TMath::Abs(track->GetLabel());
 487       Int_t mcGen = 1;
 488       // reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
 489       if (fSelectHIJING && (label==0 || track->GetGeneratorIndex() == 0)) mcGen = 0;
 490
 491       FillDetectorLevelTHnSparse(fCent, track->Eta(), track->Phi(), track->Pt(), sigma, mcGen, type);
 492
 493       if (fGeneratorLevel && label > 0) {
 494         AliAODMCParticle *part =  static_cast<AliAODMCParticle*>(fGeneratorLevel->GetAcceptParticleWithLabel(label));
 495         if (part) {
 496           if (!fSelectHIJING || part->GetGeneratorIndex() == 0) {
 497             Int_t pdg = TMath::Abs(part->PdgCode());
 498             // select charged pions, protons, kaons , electrons, muons
 499             if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
 500               FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
 501             }
 502           }
 503         }
 504       }
 505     }
 506     else {
 507       AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
 508     }
 509
 510     track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle());
 511   }
 512
 513   if (fGeneratorLevel) {
 514     AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle(0));
 515     while (part != 0) {
 516       Int_t mcGen = 1;
 517       Byte_t findable = 0;
 518
 519       if (fSelectHIJING && part->GetGeneratorIndex() == 0) mcGen = 0;
 520
 521       Int_t pdg = TMath::Abs(part->PdgCode());
 522       // select charged pions, protons, kaons , electrons, muons
 523       if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) findable = 1;
 524
 525       FillGeneratorLevelTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), mcGen, findable);
 526       part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle());
 527     }
 528   }
 529
 530   return kTRUE;
 531 }