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 "AliAODMCParticle.h"
  13 #include "AliParticleContainer.h"
  14 #include "AliLog.h"
  15
  16 #include "AliEmcalTrackingQATask.h"
  17
  18 ClassImp(AliEmcalTrackingQATask)
  19
  20 //________________________________________________________________________
  21 AliEmcalTrackingQATask::AliEmcalTrackingQATask() :
  22   AliAnalysisTaskEmcal("AliEmcalTrackingQA", kTRUE),
  23   fSelectHIJING(kTRUE),
  24   fGeneratorLevel(0),
  25   fDetectorLevel(0),
  26   fTracksAll(0),
  27   fTracksSelected(0),
  28   fParticlesAllPhysPrim(0),
  29   fParticlesSelected(0),
  30   fParticlesFindable(0),
  31   fParticlesMatched(0)
  32 {
  33   // Default constructor.
  34
  35   SetMakeGeneralHistograms(kTRUE);
  36 }
  37
  38 //________________________________________________________________________
  39 AliEmcalTrackingQATask::AliEmcalTrackingQATask(const char *name) :
  40   AliAnalysisTaskEmcal(name, kTRUE),
  41   fSelectHIJING(kTRUE),
  42   fGeneratorLevel(0),
  43   fDetectorLevel(0),
  44   fTracksAll(0),
  45   fTracksSelected(0),
  46   fParticlesAllPhysPrim(0),
  47   fParticlesSelected(0),
  48   fParticlesFindable(0),
  49   fParticlesMatched(0)
  50 {
  51   // Standard constructor.
  52
  53   SetMakeGeneralHistograms(kTRUE);
  54 }
  55
  56 //________________________________________________________________________
  57 AliEmcalTrackingQATask::~AliEmcalTrackingQATask()
  58 {
  59   // Destructor.
  60 }
  61
  62 //________________________________________________________________________
  63 void AliEmcalTrackingQATask::UserCreateOutputObjects()
  64 {
  65   // Create my user objects.
  66
  67   AliAnalysisTaskEmcal::UserCreateOutputObjects();
  68
  69   if (!fCreateHisto) return;
  70
  71   fTracksAll = new TH3**[fNcentBins];
  72   fTracksSelected = new TH3**[fNcentBins];
  73   fParticlesAllPhysPrim = new TH3*[fNcentBins];
  74   fParticlesSelected = new TH3*[fNcentBins];
  75
  76   TString histname;
  77
  78   for (Int_t i = 0; i < fNcentBins; i++) {
  79
  80     fTracksAll[i] = new TH3*[3];
  81     fTracksSelected[i] = new TH3*[3];
  82     for (Int_t j = 0; j < 3; j++) {
  83       histname = Form("fTracksAll_%d_%d",i,j);
  84       fTracksAll[i][j] = new TH3F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
  85       fTracksAll[i][j]->GetXaxis()->SetTitle("#eta");
  86       fTracksAll[i][j]->GetYaxis()->SetTitle("#phi");
  87       fTracksAll[i][j]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
  88       fOutput->Add(fTracksAll[i][j]);
  89
  90       if (fSelectHIJING) {
  91         histname = Form("fTracksSelected_%d_%d",i,j);
  92         fTracksSelected[i][j] = new TH3F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
  93         fTracksSelected[i][j]->GetXaxis()->SetTitle("#eta");
  94         fTracksSelected[i][j]->GetYaxis()->SetTitle("#phi");
  95         fTracksSelected[i][j]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
  96         fOutput->Add(fTracksSelected[i][j]);
  97       }
  98     }
  99
 100     histname = Form("fParticlesAllPhysPrim_%d",i);
 101     fParticlesAllPhysPrim[i] = new TH3F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
 102     fParticlesAllPhysPrim[i]->GetXaxis()->SetTitle("#eta");
 103     fParticlesAllPhysPrim[i]->GetYaxis()->SetTitle("#phi");
 104     fParticlesAllPhysPrim[i]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 105     fOutput->Add(fParticlesAllPhysPrim[i]);
 106
 107     if (fSelectHIJING) {
 108       histname = Form("fParticlesSelected_%d",i);
 109       fParticlesSelected[i] = new TH3F(histname,histname, 100, -1, 1, 101, 0, TMath::Pi() * 2.02, fNbins, fMinBinPt, fMaxBinPt);
 110       fParticlesSelected[i]->GetXaxis()->SetTitle("#eta");
 111       fParticlesSelected[i]->GetYaxis()->SetTitle("#phi");
 112       fParticlesSelected[i]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 113       fOutput->Add(fParticlesSelected[i]);
 114     }
 115   }
 116
 117   AllocateFindableParticlesTHnSparse();
 118   AllocateMatchedParticlesTHnSparse();
 119
 120   PostData(1, fOutput);
 121 }
 122
 123 //________________________________________________________________________
 124 void AliEmcalTrackingQATask::AllocateFindableParticlesTHnSparse()
 125 {
 126   Int_t dim = 0;
 127   TString title[20];
 128   Int_t nbins[20] = {0};
 129   Double_t min[20] = {0};
 130   Double_t max[20] = {0};
 131
 132   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 133     title[dim] = "Centrality %";
 134     nbins[dim] = 101;
 135     min[dim] = 0;
 136     max[dim] = 101;
 137     dim++;
 138   }
 139
 140   title[dim] = "#it{p}_{T} (GeV/#it{c})";
 141   nbins[dim] = fNbins;
 142   min[dim] = fMinBinPt;
 143   max[dim] = fMaxBinPt;
 144   dim++;
 145
 146   title[dim] = "#eta";
 147   nbins[dim] = 100;
 148   min[dim] = -1;
 149   max[dim] = 1;
 150   dim++;
 151
 152   title[dim] = "#phi";
 153   nbins[dim] = 101;
 154   min[dim] = 0;
 155   max[dim] = TMath::Pi() * 2.02;
 156   dim++;
 157
 158   fParticlesFindable = new THnSparseF("fParticlesFindable","fParticlesFindable",dim,nbins,min,max);
 159   for (Int_t i = 0; i < dim; i++)
 160     fParticlesFindable->GetAxis(i)->SetTitle(title[i]);
 161   fOutput->Add(fParticlesFindable);
 162 }
 163
 164 //________________________________________________________________________
 165 void AliEmcalTrackingQATask::AllocateMatchedParticlesTHnSparse()
 166 {
 167   Int_t dim = 0;
 168   TString title[20];
 169   Int_t nbins[20] = {0};
 170   Double_t min[20] = {0};
 171   Double_t max[20] = {0};
 172
 173   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 174     title[dim] = "Centrality %";
 175     nbins[dim] = 101;
 176     min[dim] = 0;
 177     max[dim] = 101;
 178     dim++;
 179   }
 180
 181   title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
 182   nbins[dim] = fNbins;
 183   min[dim] = fMinBinPt;
 184   max[dim] = fMaxBinPt;
 185   dim++;
 186
 187   title[dim] = "#eta^{gen}";
 188   nbins[dim] = 100;
 189   min[dim] = -1;
 190   max[dim] = 1;
 191   dim++;
 192
 193   title[dim] = "#phi^{gen}";
 194   nbins[dim] = 101;
 195   min[dim] = 0;
 196   max[dim] = TMath::Pi() * 2.02;
 197   dim++;
 198
 199   title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
 200   nbins[dim] = fNbins;
 201   min[dim] = fMinBinPt;
 202   max[dim] = fMaxBinPt;
 203   dim++;
 204
 205   title[dim] = "#eta^{det}";
 206   nbins[dim] = 100;
 207   min[dim] = -1;
 208   max[dim] = 1;
 209   dim++;
 210
 211   title[dim] = "#phi^{det}";
 212   nbins[dim] = 101;
 213   min[dim] = 0;
 214   max[dim] = TMath::Pi() * 2.02;
 215   dim++;
 216
 217   title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
 218   nbins[dim] = fNbins;
 219   min[dim] = -1;
 220   max[dim] = 1;
 221   dim++;
 222
 223   title[dim] = "track type";
 224   nbins[dim] = 3;
 225   min[dim] = -0.5;
 226   max[dim] = 2.5;
 227   dim++;
 228
 229   fParticlesMatched = new THnSparseF("fParticlesMatched","fParticlesMatched",dim,nbins,min,max);
 230   for (Int_t i = 0; i < dim; i++)
 231     fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
 232   fOutput->Add(fParticlesMatched);
 233 }
 234
 235 //________________________________________________________________________
 236 void AliEmcalTrackingQATask::SetGeneratorLevelName(const char* name)
 237 {
 238   if (!fGeneratorLevel) {  // first check if the generator level array is set
 239     fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 240     if (fGeneratorLevel) {  // now check if the first collection array has been added already
 241       fGeneratorLevel->SetArrayName(name);
 242     }
 243     else {
 244       fGeneratorLevel = AddParticleContainer(name);
 245     }
 246     fGeneratorLevel->SetClassName("AliAODMCParticle");
 247     fGeneratorLevel->SelectPhysicalPrimaries(kTRUE);
 248   }
 249   fGeneratorLevel->SetArrayName(name);
 250 }
 251
 252 //________________________________________________________________________
 253 void AliEmcalTrackingQATask::SetDetectorLevelName(const char* name)
 254 {
 255   if (!fGeneratorLevel) {
 256     AliError("Please, first set the generatol level array!");
 257     return;
 258   }
 259   if (!fDetectorLevel) {  // first check if the detector level array is set
 260     fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 261     if (fDetectorLevel) {  // now check if the second collection array has been added already
 262       fDetectorLevel->SetArrayName(name);
 263     }
 264     else {
 265       fDetectorLevel = AddParticleContainer(name);
 266     }
 267     fDetectorLevel->SetClassName("AliPicoTrack");
 268   }
 269   fDetectorLevel->SetArrayName(name);
 270 }
 271
 272 //________________________________________________________________________
 273 void AliEmcalTrackingQATask::ExecOnce()
 274 {
 275   // Init the analysis.
 276
 277   if (fParticleCollArray.GetEntriesFast() < 2) {
 278     AliFatal("This task needs at least two particle containers!");
 279   }
 280
 281   if (!fGeneratorLevel) {
 282     fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 283     fGeneratorLevel->SetClassName("AliAODMCParticle");
 284   }
 285
 286   if (!fDetectorLevel) {
 287     fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 288     fDetectorLevel->SetClassName("AliPicoTrack");
 289   }
 290
 291   AliAnalysisTaskEmcal::ExecOnce();
 292 }
 293
 294 //________________________________________________________________________
 295 void AliEmcalTrackingQATask::FillFindableParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt)
 296 {
 297   Double_t contents[20]={0};
 298
 299   for (Int_t i = 0; i < fParticlesFindable->GetNdimensions(); i++) {
 300     TString title(fParticlesFindable->GetAxis(i)->GetTitle());
 301     if (title=="Centrality %")
 302       contents[i] = cent;
 303     else if (title=="#it{p}_{T} (GeV/#it{c})")
 304       contents[i] = partPt;
 305     else if (title=="#eta")
 306       contents[i] = partEta;
 307     else if (title=="#phi")
 308       contents[i] = partPhi;
 309     else
 310       AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesFindable->GetName()));
 311   }
 312
 313   fParticlesFindable->Fill(contents);
 314 }
 315
 316 //________________________________________________________________________
 317 void AliEmcalTrackingQATask::FillMatchedParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt,
 318                                                            Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
 319 {
 320   Double_t contents[20]={0};
 321
 322   for (Int_t i = 0; i < fParticlesMatched->GetNdimensions(); i++) {
 323     TString title(fParticlesMatched->GetAxis(i)->GetTitle());
 324     if (title=="Centrality %")
 325       contents[i] = cent;
 326     else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
 327       contents[i] = partPt;
 328     else if (title=="#eta^{gen}")
 329       contents[i] = partEta;
 330     else if (title=="#phi^{gen}")
 331       contents[i] = partPhi;
 332     else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
 333       contents[i] = trackPt;
 334     else if (title=="#eta^{det}")
 335       contents[i] = trackEta;
 336     else if (title=="#phi^{det}")
 337       contents[i] = trackPhi;
 338     else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
 339       contents[i] = (partPt - trackPt) / partPt;
 340     else if (title=="track type")
 341       contents[i] = (Double_t)trackType;
 342     else
 343       AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesMatched->GetName()));
 344   }
 345
 346   fParticlesMatched->Fill(contents);
 347 }
 348
 349 //________________________________________________________________________
 350 Bool_t AliEmcalTrackingQATask::FillHistograms()
 351 {
 352   // Fill the histograms.
 353
 354   AliPicoTrack *track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle(0));
 355   while (track != 0) {
 356     Byte_t type = track->GetTrackType();
 357     if (type<= 2 && type >= 0) {
 358       fTracksAll[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
 359
 360       Int_t label = TMath::Abs(track->GetLabel());
 361
 362       if (fSelectHIJING && (label==0 || track->GetGeneratorIndex() == 0)) {
 363         // reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
 364         fTracksSelected[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
 365       }
 366
 367       if (label > 0) {
 368         AliAODMCParticle *part =  static_cast<AliAODMCParticle*>(fGeneratorLevel->GetAcceptParticleWithLabel(label));
 369         if (part) {
 370           if (!fSelectHIJING || part->GetGeneratorIndex() == 0) {
 371             Int_t pdg = TMath::Abs(part->PdgCode());
 372             // select charged pions, protons, kaons , electrons, muons
 373             if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
 374               FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
 375             }
 376           }
 377         }
 378       }
 379     }
 380     else {
 381       AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
 382     }
 383
 384     track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle());
 385   }
 386
 387   AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle(0));
 388   while (part != 0) {
 389     fParticlesAllPhysPrim[fCentBin]->Fill(part->Eta(), part->Phi(), part->Pt());
 390     if (!fSelectHIJING || part->GetGeneratorIndex() == 0) {
 391       if (fSelectHIJING) fParticlesSelected[fCentBin]->Fill(part->Eta(), part->Phi(), part->Pt());
 392
 393       Int_t pdg = TMath::Abs(part->PdgCode());
 394       // select charged pions, protons, kaons , electrons, muons
 395       if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
 396         FillFindableParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt());
 397       }
 398     }
 399
 400     part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle());
 401   }
 402
 403   return kTRUE;
 404 }