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   fNPtHistBins(0),
  27   fPtHistBins(0),
  28   fNEtaHistBins(0),
  29   fEtaHistBins(0),
  30   fNPhiHistBins(0),
  31   fPhiHistBins(0),
  32   fNCentHistBins(0),
  33   fCentHistBins(0),
  34   fNPtResHistBins(0),
  35   fPtResHistBins(0),
  36   fNIntegerHistBins(0),
  37   fIntegerHistBins(0),
  38   fTracksAll(0),
  39   fTracksSelected(0),
  40   fParticlesAllPhysPrim(0),
  41   fParticlesSelected(0),
  42   fParticlesFindable(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   fNIntegerHistBins(0),
  69   fIntegerHistBins(0),
  70   fTracksAll(0),
  71   fTracksSelected(0),
  72   fParticlesAllPhysPrim(0),
  73   fParticlesSelected(0),
  74   fParticlesFindable(0),
  75   fParticlesMatched(0)
  76 {
  77   // Standard constructor.
  78
  79   SetMakeGeneralHistograms(kTRUE);
  80
  81   GenerateHistoBins();
  82 }
  83
  84 //________________________________________________________________________
  85 AliEmcalTrackingQATask::~AliEmcalTrackingQATask()
  86 {
  87   // Destructor.
  88 }
  89
  90 //________________________________________________________________________
  91 void AliEmcalTrackingQATask::GenerateHistoBins()
  92 {
  93   fNPtHistBins = 79;
  94   fPtHistBins = new Double_t[fNPtHistBins+1];
  95   GenerateFixedBinArray(10, 0, 1, fPtHistBins);
  96   GenerateFixedBinArray(10, 1, 3, fPtHistBins+10);
  97   GenerateFixedBinArray(14, 3, 10, fPtHistBins+20);
  98   GenerateFixedBinArray(10, 10, 20, fPtHistBins+34);
  99   GenerateFixedBinArray(15, 20, 50, fPtHistBins+44);
 100   GenerateFixedBinArray(20, 50, 150, fPtHistBins+59);
 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   fNIntegerHistBins = 10;
 123   fIntegerHistBins = new Double_t[fNIntegerHistBins+1];
 124   GenerateFixedBinArray(fNIntegerHistBins, -0.5, 9.5, fIntegerHistBins);
 125 }
 126
 127 //________________________________________________________________________
 128 void AliEmcalTrackingQATask::UserCreateOutputObjects()
 129 {
 130   // Create my user objects.
 131
 132   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 133
 134   if (!fCreateHisto) return;
 135
 136   fTracksAll = new TH3**[fNcentBins];
 137   fTracksSelected = new TH3**[fNcentBins];
 138   fParticlesAllPhysPrim = new TH3*[fNcentBins];
 139   fParticlesSelected = new TH3*[fNcentBins];
 140
 141   TString histname;
 142
 143   for (Int_t i = 0; i < fNcentBins; i++) {
 144
 145     fTracksAll[i] = new TH3*[3];
 146     fTracksSelected[i] = new TH3*[3];
 147     for (Int_t j = 0; j < 3; j++) {
 148       histname = Form("fTracksAll_%d_%d",i,j);
 149       fTracksAll[i][j] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
 150       fTracksAll[i][j]->GetXaxis()->SetTitle("#eta");
 151       fTracksAll[i][j]->GetYaxis()->SetTitle("#phi");
 152       fTracksAll[i][j]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 153       fOutput->Add(fTracksAll[i][j]);
 154
 155       if (fSelectHIJING) {
 156         histname = Form("fTracksSelected_%d_%d",i,j);
 157         fTracksSelected[i][j] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
 158         fTracksSelected[i][j]->GetXaxis()->SetTitle("#eta");
 159         fTracksSelected[i][j]->GetYaxis()->SetTitle("#phi");
 160         fTracksSelected[i][j]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 161         fOutput->Add(fTracksSelected[i][j]);
 162       }
 163     }
 164
 165     histname = Form("fParticlesAllPhysPrim_%d",i);
 166     fParticlesAllPhysPrim[i] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
 167     fParticlesAllPhysPrim[i]->GetXaxis()->SetTitle("#eta");
 168     fParticlesAllPhysPrim[i]->GetYaxis()->SetTitle("#phi");
 169     fParticlesAllPhysPrim[i]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 170     fOutput->Add(fParticlesAllPhysPrim[i]);
 171
 172     if (fSelectHIJING) {
 173       histname = Form("fParticlesSelected_%d",i);
 174       fParticlesSelected[i] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
 175       fParticlesSelected[i]->GetXaxis()->SetTitle("#eta");
 176       fParticlesSelected[i]->GetYaxis()->SetTitle("#phi");
 177       fParticlesSelected[i]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 178       fOutput->Add(fParticlesSelected[i]);
 179     }
 180   }
 181
 182   AllocateFindableParticlesTHnSparse();
 183   AllocateMatchedParticlesTHnSparse();
 184
 185   PostData(1, fOutput);
 186 }
 187
 188 //________________________________________________________________________
 189 void AliEmcalTrackingQATask::AllocateFindableParticlesTHnSparse()
 190 {
 191   Int_t dim = 0;
 192   TString title[20];
 193   Int_t nbins[20] = {0};
 194   Double_t *binEdges[20] = {0};
 195
 196   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 197     title[dim] = "Centrality %";
 198     nbins[dim] = fNCentHistBins;
 199     binEdges[dim] = fCentHistBins;
 200     dim++;
 201   }
 202
 203   title[dim] = "#it{p}_{T} (GeV/#it{c})";
 204   nbins[dim] = fNPtHistBins;
 205   binEdges[dim] = fPtHistBins;
 206   dim++;
 207
 208   title[dim] = "#eta";
 209   nbins[dim] = fNEtaHistBins;
 210   binEdges[dim] = fEtaHistBins;
 211   dim++;
 212
 213   title[dim] = "#phi";
 214   nbins[dim] = fNPhiHistBins;
 215   binEdges[dim] = fPhiHistBins;
 216   dim++;
 217
 218   fParticlesFindable = new THnSparseF("fParticlesFindable","fParticlesFindable",dim,nbins);
 219   for (Int_t i = 0; i < dim; i++) {
 220     fParticlesFindable->GetAxis(i)->SetTitle(title[i]);
 221     fParticlesFindable->SetBinEdges(i, binEdges[i]);
 222   }
 223
 224   fOutput->Add(fParticlesFindable);
 225 }
 226
 227 //________________________________________________________________________
 228 void AliEmcalTrackingQATask::AllocateMatchedParticlesTHnSparse()
 229 {
 230   Int_t dim = 0;
 231   TString title[20];
 232   Int_t nbins[20] = {0};
 233   Double_t *binEdges[20] = {0};
 234
 235   if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
 236     title[dim] = "Centrality %";
 237     nbins[dim] = fNCentHistBins;
 238     binEdges[dim] = fCentHistBins;
 239     dim++;
 240   }
 241
 242   title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
 243   nbins[dim] = fNPtHistBins;
 244   binEdges[dim] = fPtHistBins;
 245   dim++;
 246
 247   title[dim] = "#eta^{gen}";
 248   nbins[dim] = fNEtaHistBins;
 249   binEdges[dim] = fEtaHistBins;
 250   dim++;
 251
 252   title[dim] = "#phi^{gen}";
 253   nbins[dim] = fNPhiHistBins;
 254   binEdges[dim] = fPhiHistBins;
 255   dim++;
 256
 257   title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
 258   nbins[dim] = fNPtHistBins;
 259   binEdges[dim] = fPtHistBins;
 260   dim++;
 261
 262   title[dim] = "#eta^{det}";
 263   nbins[dim] = fNEtaHistBins;
 264   binEdges[dim] = fEtaHistBins;
 265   dim++;
 266
 267   title[dim] = "#phi^{det}";
 268   nbins[dim] = fNPhiHistBins;
 269   binEdges[dim] = fPhiHistBins;
 270   dim++;
 271
 272   title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
 273   nbins[dim] = fNPtResHistBins;
 274   binEdges[dim] = fPtResHistBins;
 275   dim++;
 276
 277   title[dim] = "track type";
 278   nbins[dim] = 3;
 279   binEdges[dim] = fIntegerHistBins;
 280   dim++;
 281
 282   fParticlesMatched = new THnSparseF("fParticlesMatched","fParticlesMatched",dim,nbins);
 283   for (Int_t i = 0; i < dim; i++) {
 284     fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
 285     fParticlesMatched->SetBinEdges(i, binEdges[i]);
 286   }
 287
 288   fOutput->Add(fParticlesMatched);
 289 }
 290
 291 //________________________________________________________________________
 292 void AliEmcalTrackingQATask::SetGeneratorLevelName(const char* name)
 293 {
 294   if (!fGeneratorLevel) {  // first check if the generator level array is set
 295     fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 296     if (fGeneratorLevel) {  // now check if the first collection array has been added already
 297       fGeneratorLevel->SetArrayName(name);
 298     }
 299     else {
 300       fGeneratorLevel = AddParticleContainer(name);
 301     }
 302     fGeneratorLevel->SetClassName("AliAODMCParticle");
 303     fGeneratorLevel->SelectPhysicalPrimaries(kTRUE);
 304   }
 305   fGeneratorLevel->SetArrayName(name);
 306 }
 307
 308 //________________________________________________________________________
 309 void AliEmcalTrackingQATask::SetDetectorLevelName(const char* name)
 310 {
 311   if (!fGeneratorLevel) {
 312     AliError("Please, first set the generatol level array!");
 313     return;
 314   }
 315   if (!fDetectorLevel) {  // first check if the detector level array is set
 316     fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 317     if (fDetectorLevel) {  // now check if the second collection array has been added already
 318       fDetectorLevel->SetArrayName(name);
 319     }
 320     else {
 321       fDetectorLevel = AddParticleContainer(name);
 322     }
 323     fDetectorLevel->SetClassName("AliPicoTrack");
 324   }
 325   fDetectorLevel->SetArrayName(name);
 326 }
 327
 328 //________________________________________________________________________
 329 void AliEmcalTrackingQATask::ExecOnce()
 330 {
 331   // Init the analysis.
 332
 333   if (fParticleCollArray.GetEntriesFast() < 2) {
 334     AliFatal("This task needs at least two particle containers!");
 335   }
 336
 337   if (!fGeneratorLevel) {
 338     fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 339     fGeneratorLevel->SetClassName("AliAODMCParticle");
 340   }
 341
 342   if (!fDetectorLevel) {
 343     fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 344     fDetectorLevel->SetClassName("AliPicoTrack");
 345   }
 346
 347   AliAnalysisTaskEmcal::ExecOnce();
 348 }
 349
 350 //________________________________________________________________________
 351 void AliEmcalTrackingQATask::FillFindableParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt)
 352 {
 353   Double_t contents[20]={0};
 354
 355   for (Int_t i = 0; i < fParticlesFindable->GetNdimensions(); i++) {
 356     TString title(fParticlesFindable->GetAxis(i)->GetTitle());
 357     if (title=="Centrality %")
 358       contents[i] = cent;
 359     else if (title=="#it{p}_{T} (GeV/#it{c})")
 360       contents[i] = partPt;
 361     else if (title=="#eta")
 362       contents[i] = partEta;
 363     else if (title=="#phi")
 364       contents[i] = partPhi;
 365     else
 366       AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesFindable->GetName()));
 367   }
 368
 369   fParticlesFindable->Fill(contents);
 370 }
 371
 372 //________________________________________________________________________
 373 void AliEmcalTrackingQATask::FillMatchedParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt,
 374                                                            Double_t trackEta, Double_t trackPhi, Double_t trackPt, Byte_t trackType)
 375 {
 376   Double_t contents[20]={0};
 377
 378   for (Int_t i = 0; i < fParticlesMatched->GetNdimensions(); i++) {
 379     TString title(fParticlesMatched->GetAxis(i)->GetTitle());
 380     if (title=="Centrality %")
 381       contents[i] = cent;
 382     else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
 383       contents[i] = partPt;
 384     else if (title=="#eta^{gen}")
 385       contents[i] = partEta;
 386     else if (title=="#phi^{gen}")
 387       contents[i] = partPhi;
 388     else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
 389       contents[i] = trackPt;
 390     else if (title=="#eta^{det}")
 391       contents[i] = trackEta;
 392     else if (title=="#phi^{det}")
 393       contents[i] = trackPhi;
 394     else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
 395       contents[i] = (partPt - trackPt) / partPt;
 396     else if (title=="track type")
 397       contents[i] = (Double_t)trackType;
 398     else
 399       AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesMatched->GetName()));
 400   }
 401
 402   fParticlesMatched->Fill(contents);
 403 }
 404
 405 //________________________________________________________________________
 406 Bool_t AliEmcalTrackingQATask::FillHistograms()
 407 {
 408   // Fill the histograms.
 409
 410   AliPicoTrack *track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle(0));
 411   while (track != 0) {
 412     Byte_t type = track->GetTrackType();
 413     if (type <= 2) {
 414       fTracksAll[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
 415
 416       Int_t label = TMath::Abs(track->GetLabel());
 417
 418       if (fSelectHIJING && (label==0 || track->GetGeneratorIndex() == 0)) {
 419         // reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
 420         fTracksSelected[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
 421       }
 422
 423       if (label > 0) {
 424         AliAODMCParticle *part =  static_cast<AliAODMCParticle*>(fGeneratorLevel->GetAcceptParticleWithLabel(label));
 425         if (part) {
 426           if (!fSelectHIJING || part->GetGeneratorIndex() == 0) {
 427             Int_t pdg = TMath::Abs(part->PdgCode());
 428             // select charged pions, protons, kaons , electrons, muons
 429             if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
 430               FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
 431             }
 432           }
 433         }
 434       }
 435     }
 436     else {
 437       AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
 438     }
 439
 440     track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle());
 441   }
 442
 443   AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle(0));
 444   while (part != 0) {
 445     fParticlesAllPhysPrim[fCentBin]->Fill(part->Eta(), part->Phi(), part->Pt());
 446     if (!fSelectHIJING || part->GetGeneratorIndex() == 0) {
 447       if (fSelectHIJING) fParticlesSelected[fCentBin]->Fill(part->Eta(), part->Phi(), part->Pt());
 448
 449       Int_t pdg = TMath::Abs(part->PdgCode());
 450       // select charged pions, protons, kaons , electrons, muons
 451       if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
 452         FillFindableParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt());
 453       }
 454     }
 455
 456     part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle());
 457   }
 458
 459   return kTRUE;
 460 }