[u/mrichter/AliRoot.git] / PWG / EMCAL / AliEmcalTrackingQATask.cxx

// Track QA task (efficiency and pt resolution)
//
// Author: S.Aiola

#include <TH3F.h>
#include <THnSparse.h>
#include <TMath.h>
#include <TString.h>
#include <Riostream.h>

#include "AliPicoTrack.h"
#include "AliAODMCParticle.h"
#include "AliParticleContainer.h"
#include "AliLog.h"

#include "AliEmcalTrackingQATask.h"

ClassImp(AliEmcalTrackingQATask)

//________________________________________________________________________
AliEmcalTrackingQATask::AliEmcalTrackingQATask() : 
  AliAnalysisTaskEmcal("AliEmcalTrackingQA", kTRUE),
  fSelectHIJING(kTRUE),
  fGeneratorLevel(0),
  fDetectorLevel(0),
  fNPtHistBins(0),
  fPtHistBins(0),
  fNEtaHistBins(0),
  fEtaHistBins(0),
  fNPhiHistBins(0),
  fPhiHistBins(0),
  fNCentHistBins(0),
  fCentHistBins(0),
  fNPtResHistBins(0),
  fPtResHistBins(0),
  fNIntegerHistBins(0),
  fIntegerHistBins(0),
  fTracksAll(0),
  fTracksSelected(0),
  fParticlesAllPhysPrim(0),
  fParticlesSelected(0),
  fParticlesFindable(0),
  fParticlesMatched(0)
{
  // Default constructor.

  SetMakeGeneralHistograms(kTRUE);

  GenerateHistoBins();
}

//________________________________________________________________________
AliEmcalTrackingQATask::AliEmcalTrackingQATask(const char *name) : 
  AliAnalysisTaskEmcal(name, kTRUE),
  fSelectHIJING(kTRUE),
  fGeneratorLevel(0),
  fDetectorLevel(0),
  fNPtHistBins(0),
  fPtHistBins(0),
  fNEtaHistBins(0),
  fEtaHistBins(0),
  fNPhiHistBins(0),
  fPhiHistBins(0),
  fNCentHistBins(0),
  fCentHistBins(0),
  fNPtResHistBins(0),
  fPtResHistBins(0),
  fNIntegerHistBins(0),
  fIntegerHistBins(0),
  fTracksAll(0),
  fTracksSelected(0),
  fParticlesAllPhysPrim(0),
  fParticlesSelected(0),
  fParticlesFindable(0),
  fParticlesMatched(0)
{
  // Standard constructor.

  SetMakeGeneralHistograms(kTRUE);

  GenerateHistoBins();
}

//________________________________________________________________________
AliEmcalTrackingQATask::~AliEmcalTrackingQATask()
{
  // Destructor.
}

//________________________________________________________________________
void AliEmcalTrackingQATask::GenerateHistoBins()
{
  fNPtHistBins = 79;
  fPtHistBins = new Double_t[fNPtHistBins+1];
  GenerateFixedBinArray(10, 0, 1, fPtHistBins);
  GenerateFixedBinArray(10, 1, 3, fPtHistBins+10);
  GenerateFixedBinArray(14, 3, 10, fPtHistBins+20);
  GenerateFixedBinArray(10, 10, 20, fPtHistBins+34);
  GenerateFixedBinArray(15, 20, 50, fPtHistBins+44);
  GenerateFixedBinArray(20, 50, 150, fPtHistBins+59);

  fNEtaHistBins = 100;
  fEtaHistBins = new Double_t[fNEtaHistBins+1];
  GenerateFixedBinArray(fNEtaHistBins, -1, 1, fEtaHistBins);

  fNPhiHistBins = 101;
  fPhiHistBins = new Double_t[fNPhiHistBins+1];
  GenerateFixedBinArray(fNPhiHistBins, 0, TMath::Pi() * 2.02, fPhiHistBins);

  fNCentHistBins = 4;
  fCentHistBins = new Double_t[fNCentHistBins+1];
  fCentHistBins[0] = 0;
  fCentHistBins[1] = 10;
  fCentHistBins[2] = 30;
  fCentHistBins[3] = 50;
  fCentHistBins[4] = 90;

  fNPtResHistBins = 200;
  fPtResHistBins = new Double_t[fNPtResHistBins+1];
  GenerateFixedBinArray(fNPtResHistBins, -2, 2, fPtResHistBins);

  fNIntegerHistBins = 10;
  fIntegerHistBins = new Double_t[fNIntegerHistBins+1];
  GenerateFixedBinArray(fNIntegerHistBins, -0.5, 9.5, fIntegerHistBins);
}

//________________________________________________________________________
void AliEmcalTrackingQATask::UserCreateOutputObjects()
{
  // Create my user objects.

  AliAnalysisTaskEmcal::UserCreateOutputObjects();

  if (!fCreateHisto) return;

  fTracksAll = new TH3**[fNcentBins];
  fTracksSelected = new TH3**[fNcentBins];
  fParticlesAllPhysPrim = new TH3*[fNcentBins];
  fParticlesSelected = new TH3*[fNcentBins];  

  TString histname;

  for (Int_t i = 0; i < fNcentBins; i++) {

    fTracksAll[i] = new TH3*[3];
    fTracksSelected[i] = new TH3*[3];
    for (Int_t j = 0; j < 3; j++) {
      histname = Form("fTracksAll_%d_%d",i,j);
      fTracksAll[i][j] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
      fTracksAll[i][j]->GetXaxis()->SetTitle("#eta");
      fTracksAll[i][j]->GetYaxis()->SetTitle("#phi");
      fTracksAll[i][j]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
      fOutput->Add(fTracksAll[i][j]);

      if (fSelectHIJING) {
	histname = Form("fTracksSelected_%d_%d",i,j);
	fTracksSelected[i][j] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
	fTracksSelected[i][j]->GetXaxis()->SetTitle("#eta");
	fTracksSelected[i][j]->GetYaxis()->SetTitle("#phi");
	fTracksSelected[i][j]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
	fOutput->Add(fTracksSelected[i][j]);
      }
    }

    histname = Form("fParticlesAllPhysPrim_%d",i);
    fParticlesAllPhysPrim[i] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
    fParticlesAllPhysPrim[i]->GetXaxis()->SetTitle("#eta");
    fParticlesAllPhysPrim[i]->GetYaxis()->SetTitle("#phi");
    fParticlesAllPhysPrim[i]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
    fOutput->Add(fParticlesAllPhysPrim[i]);

    if (fSelectHIJING) {
      histname = Form("fParticlesSelected_%d",i);
      fParticlesSelected[i] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
      fParticlesSelected[i]->GetXaxis()->SetTitle("#eta");
      fParticlesSelected[i]->GetYaxis()->SetTitle("#phi");
      fParticlesSelected[i]->GetZaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
      fOutput->Add(fParticlesSelected[i]);
    }
  }

  AllocateFindableParticlesTHnSparse();
  AllocateMatchedParticlesTHnSparse();

  PostData(1, fOutput);
}

//________________________________________________________________________
void AliEmcalTrackingQATask::AllocateFindableParticlesTHnSparse()
{
  Int_t dim = 0;
  TString title[20];
  Int_t nbins[20] = {0};
  Double_t *binEdges[20] = {0};
  
  if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
    title[dim] = "Centrality %";
    nbins[dim] = fNCentHistBins;
    binEdges[dim] = fCentHistBins;
    dim++;
  }

  title[dim] = "#it{p}_{T} (GeV/#it{c})";
  nbins[dim] = fNPtHistBins;
  binEdges[dim] = fPtHistBins;
  dim++;

  title[dim] = "#eta";
  nbins[dim] = fNEtaHistBins;
  binEdges[dim] = fEtaHistBins;
  dim++;

  title[dim] = "#phi";
  nbins[dim] = fNPhiHistBins;
  binEdges[dim] = fPhiHistBins;
  dim++;
 
  fParticlesFindable = new THnSparseF("fParticlesFindable","fParticlesFindable",dim,nbins);
  for (Int_t i = 0; i < dim; i++) {
    fParticlesFindable->GetAxis(i)->SetTitle(title[i]);
    fParticlesFindable->SetBinEdges(i, binEdges[i]);
  }

  fOutput->Add(fParticlesFindable);
}

//________________________________________________________________________
void AliEmcalTrackingQATask::AllocateMatchedParticlesTHnSparse()
{
  Int_t dim = 0;
  TString title[20];
  Int_t nbins[20] = {0};
  Double_t *binEdges[20] = {0};
  
  if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
    title[dim] = "Centrality %";
    nbins[dim] = fNCentHistBins;
    binEdges[dim] = fCentHistBins;
    dim++;
  }

  title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
  nbins[dim] = fNPtHistBins;
  binEdges[dim] = fPtHistBins;
  dim++;

  title[dim] = "#eta^{gen}";
  nbins[dim] = fNEtaHistBins;
  binEdges[dim] = fEtaHistBins;
  dim++;

  title[dim] = "#phi^{gen}";
  nbins[dim] = fNPhiHistBins;
  binEdges[dim] = fPhiHistBins;
  dim++;

  title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
  nbins[dim] = fNPtHistBins;
  binEdges[dim] = fPtHistBins;
  dim++;

  title[dim] = "#eta^{det}";
  nbins[dim] = fNEtaHistBins;
  binEdges[dim] = fEtaHistBins;
  dim++;

  title[dim] = "#phi^{det}";
  nbins[dim] = fNPhiHistBins;
  binEdges[dim] = fPhiHistBins;
  dim++;

  title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
  nbins[dim] = fNPtResHistBins;
  binEdges[dim] = fPtResHistBins;
  dim++;

  title[dim] = "track type";
  nbins[dim] = 3;
  binEdges[dim] = fIntegerHistBins;
  dim++;

  fParticlesMatched = new THnSparseF("fParticlesMatched","fParticlesMatched",dim,nbins);
  for (Int_t i = 0; i < dim; i++) {
    fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
    fParticlesMatched->SetBinEdges(i, binEdges[i]);
  }

  fOutput->Add(fParticlesMatched);
}

//________________________________________________________________________
void AliEmcalTrackingQATask::SetGeneratorLevelName(const char* name)
{
  if (!fGeneratorLevel) {  // first check if the generator level array is set
    fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
    if (fGeneratorLevel) {  // now check if the first collection array has been added already
      fGeneratorLevel->SetArrayName(name);
    }
    else {
      fGeneratorLevel = AddParticleContainer(name);
    }
    fGeneratorLevel->SetClassName("AliAODMCParticle");
    fGeneratorLevel->SelectPhysicalPrimaries(kTRUE);
  }
  fGeneratorLevel->SetArrayName(name);
}

//________________________________________________________________________
void AliEmcalTrackingQATask::SetDetectorLevelName(const char* name)
{
  if (!fGeneratorLevel) {
    AliError("Please, first set the generatol level array!");
    return;
  }
  if (!fDetectorLevel) {  // first check if the detector level array is set
    fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
    if (fDetectorLevel) {  // now check if the second collection array has been added already
      fDetectorLevel->SetArrayName(name);
    }
    else {
      fDetectorLevel = AddParticleContainer(name);
    }
    fDetectorLevel->SetClassName("AliPicoTrack");
  }
  fDetectorLevel->SetArrayName(name);
}

//________________________________________________________________________
void AliEmcalTrackingQATask::ExecOnce()
{
  // Init the analysis.

  if (fParticleCollArray.GetEntriesFast() < 2) {
    AliFatal("This task needs at least two particle containers!");
  }

  if (!fGeneratorLevel) {
    fGeneratorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
    fGeneratorLevel->SetClassName("AliAODMCParticle");
  }

  if (!fDetectorLevel) {
    fDetectorLevel = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
    fDetectorLevel->SetClassName("AliPicoTrack");
  }

  AliAnalysisTaskEmcal::ExecOnce();
}

//________________________________________________________________________
void AliEmcalTrackingQATask::FillFindableParticlesTHnSparse(Double_t cent, Double_t partEta, Double_t partPhi, Double_t partPt)
{
  Double_t contents[20]={0};

  for (Int_t i = 0; i < fParticlesFindable->GetNdimensions(); i++) {
    TString title(fParticlesFindable->GetAxis(i)->GetTitle());
    if (title=="Centrality %")
      contents[i] = cent;
    else if (title=="#it{p}_{T} (GeV/#it{c})")
      contents[i] = partPt;
    else if (title=="#eta")
      contents[i] = partEta;
    else if (title=="#phi")
      contents[i] = partPhi;
    else 
      AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesFindable->GetName()));
  }

  fParticlesFindable->Fill(contents);
}

//________________________________________________________________________
void AliEmcalTrackingQATask::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)
{
  Double_t contents[20]={0};

  for (Int_t i = 0; i < fParticlesMatched->GetNdimensions(); i++) {
    TString title(fParticlesMatched->GetAxis(i)->GetTitle());
    if (title=="Centrality %")
      contents[i] = cent;
    else if (title=="#it{p}_{T}^{gen} (GeV/#it{c})")
      contents[i] = partPt;
    else if (title=="#eta^{gen}")
      contents[i] = partEta;
    else if (title=="#phi^{gen}")
      contents[i] = partPhi;
    else if (title=="#it{p}_{T}^{det} (GeV/#it{c})")
      contents[i] = trackPt;
    else if (title=="#eta^{det}")
      contents[i] = trackEta;
    else if (title=="#phi^{det}")
      contents[i] = trackPhi;
    else if (title=="(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}")
      contents[i] = (partPt - trackPt) / partPt;
    else if (title=="track type")
      contents[i] = (Double_t)trackType;
    else 
      AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fParticlesMatched->GetName()));
  }

  fParticlesMatched->Fill(contents);
}

//________________________________________________________________________
Bool_t AliEmcalTrackingQATask::FillHistograms()
{
  // Fill the histograms.

  AliPicoTrack *track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle(0));
  while (track != 0) {
    Byte_t type = track->GetTrackType();
    if (type<= 2 && type >= 0) {
      fTracksAll[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());

      Int_t label = TMath::Abs(track->GetLabel());

      if (fSelectHIJING && (label==0 || track->GetGeneratorIndex() == 0)) {  
	// reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
	fTracksSelected[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
      }
      
      if (label > 0) {
	AliAODMCParticle *part =  static_cast<AliAODMCParticle*>(fGeneratorLevel->GetAcceptParticleWithLabel(label));
	if (part) {
	  if (!fSelectHIJING || part->GetGeneratorIndex() == 0) {
	    Int_t pdg = TMath::Abs(part->PdgCode());
	    // select charged pions, protons, kaons , electrons, muons
	    if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
	      FillMatchedParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt(), track->Eta(), track->Phi(), track->Pt(), type);
	    }
	  }
	}
      }
    }
    else {
      AliError(Form("Track %d has type %d not recognized!", fDetectorLevel->GetCurrentID(), type));
    }

    track = static_cast<AliPicoTrack*>(fDetectorLevel->GetNextAcceptParticle());
  }

  AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle(0));
  while (part != 0) {
    fParticlesAllPhysPrim[fCentBin]->Fill(part->Eta(), part->Phi(), part->Pt());
    if (!fSelectHIJING || part->GetGeneratorIndex() == 0) {
      if (fSelectHIJING) fParticlesSelected[fCentBin]->Fill(part->Eta(), part->Phi(), part->Pt());

      Int_t pdg = TMath::Abs(part->PdgCode());
      // select charged pions, protons, kaons , electrons, muons
      if (pdg == 211 || pdg == 2212 || pdg == 321 || pdg == 11 || pdg == 13) {
	FillFindableParticlesTHnSparse(fCent, part->Eta(), part->Phi(), part->Pt());
      }
    }
    
    part = static_cast<AliAODMCParticle*>(fGeneratorLevel->GetNextAcceptParticle());
  }

  return kTRUE;
}
Commit	Line	Data
eec7bbb0	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),
d7ab1a93	23	fSelectHIJING(kTRUE),
eec7bbb0	24	fGeneratorLevel(0),
eec7bbb0	25	fDetectorLevel(0),
427d6cf5	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),
eec7bbb0	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);
427d6cf5	48
427d6cf5	49	GenerateHistoBins();
eec7bbb0	50	}
	51
	52	//________________________________________________________________________
	53	AliEmcalTrackingQATask::AliEmcalTrackingQATask(const char *name) :
	54	AliAnalysisTaskEmcal(name, kTRUE),
d7ab1a93	55	fSelectHIJING(kTRUE),
eec7bbb0	56	fGeneratorLevel(0),
eec7bbb0	57	fDetectorLevel(0),
427d6cf5	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),
eec7bbb0	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);
427d6cf5	80
427d6cf5	81	GenerateHistoBins();
eec7bbb0	82	}
	83
	84	//________________________________________________________________________
	85	AliEmcalTrackingQATask::~AliEmcalTrackingQATask()
	86	{
	87	// Destructor.
	88	}
	89
427d6cf5	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
eec7bbb0	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);
427d6cf5	149	fTracksAll[i][j] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
eec7bbb0	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
d7ab1a93	155	if (fSelectHIJING) {
d7ab1a93	156	histname = Form("fTracksSelected_%d_%d",i,j);
427d6cf5	157	fTracksSelected[i][j] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
d7ab1a93	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	}
eec7bbb0	163	}
	164
	165	histname = Form("fParticlesAllPhysPrim_%d",i);
427d6cf5	166	fParticlesAllPhysPrim[i] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
eec7bbb0	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
d7ab1a93	172	if (fSelectHIJING) {
d7ab1a93	173	histname = Form("fParticlesSelected_%d",i);
427d6cf5	174	fParticlesSelected[i] = new TH3F(histname, histname, fNEtaHistBins, fEtaHistBins, fNPhiHistBins, fPhiHistBins, fNPtHistBins, fPtHistBins);
d7ab1a93	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	}
eec7bbb0	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};
427d6cf5	194	Double_t *binEdges[20] = {0};
eec7bbb0	195
	196	if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
	197	title[dim] = "Centrality %";
427d6cf5	198	nbins[dim] = fNCentHistBins;
427d6cf5	199	binEdges[dim] = fCentHistBins;
eec7bbb0	200	dim++;
	201	}
	202
	203	title[dim] = "#it{p}_{T} (GeV/#it{c})";
427d6cf5	204	nbins[dim] = fNPtHistBins;
427d6cf5	205	binEdges[dim] = fPtHistBins;
eec7bbb0	206	dim++;
	207
	208	title[dim] = "#eta";
427d6cf5	209	nbins[dim] = fNEtaHistBins;
427d6cf5	210	binEdges[dim] = fEtaHistBins;
eec7bbb0	211	dim++;
	212
	213	title[dim] = "#phi";
427d6cf5	214	nbins[dim] = fNPhiHistBins;
427d6cf5	215	binEdges[dim] = fPhiHistBins;
eec7bbb0	216	dim++;
eec7bbb0	217
427d6cf5	218	fParticlesFindable = new THnSparseF("fParticlesFindable","fParticlesFindable",dim,nbins);
427d6cf5	219	for (Int_t i = 0; i < dim; i++) {
eec7bbb0	220	fParticlesFindable->GetAxis(i)->SetTitle(title[i]);
427d6cf5	221	fParticlesFindable->SetBinEdges(i, binEdges[i]);
	222	}
	223
eec7bbb0	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};
427d6cf5	233	Double_t *binEdges[20] = {0};
eec7bbb0	234
	235	if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
	236	title[dim] = "Centrality %";
427d6cf5	237	nbins[dim] = fNCentHistBins;
427d6cf5	238	binEdges[dim] = fCentHistBins;
eec7bbb0	239	dim++;
	240	}
	241
	242	title[dim] = "#it{p}_{T}^{gen} (GeV/#it{c})";
427d6cf5	243	nbins[dim] = fNPtHistBins;
427d6cf5	244	binEdges[dim] = fPtHistBins;
eec7bbb0	245	dim++;
	246
	247	title[dim] = "#eta^{gen}";
427d6cf5	248	nbins[dim] = fNEtaHistBins;
427d6cf5	249	binEdges[dim] = fEtaHistBins;
eec7bbb0	250	dim++;
	251
	252	title[dim] = "#phi^{gen}";
427d6cf5	253	nbins[dim] = fNPhiHistBins;
427d6cf5	254	binEdges[dim] = fPhiHistBins;
eec7bbb0	255	dim++;
	256
	257	title[dim] = "#it{p}_{T}^{det} (GeV/#it{c})";
427d6cf5	258	nbins[dim] = fNPtHistBins;
427d6cf5	259	binEdges[dim] = fPtHistBins;
eec7bbb0	260	dim++;
	261
	262	title[dim] = "#eta^{det}";
427d6cf5	263	nbins[dim] = fNEtaHistBins;
427d6cf5	264	binEdges[dim] = fEtaHistBins;
eec7bbb0	265	dim++;
	266
	267	title[dim] = "#phi^{det}";
427d6cf5	268	nbins[dim] = fNPhiHistBins;
427d6cf5	269	binEdges[dim] = fPhiHistBins;
eec7bbb0	270	dim++;
	271
	272	title[dim] = "(#it{p}_{T}^{gen} - #it{p}_{T}^{det}) / #it{p}_{T}^{gen}";
427d6cf5	273	nbins[dim] = fNPtResHistBins;
427d6cf5	274	binEdges[dim] = fPtResHistBins;
eec7bbb0	275	dim++;
	276
	277	title[dim] = "track type";
	278	nbins[dim] = 3;
427d6cf5	279	binEdges[dim] = fIntegerHistBins;
eec7bbb0	280	dim++;
eec7bbb0	281
427d6cf5	282	fParticlesMatched = new THnSparseF("fParticlesMatched","fParticlesMatched",dim,nbins);
427d6cf5	283	for (Int_t i = 0; i < dim; i++) {
eec7bbb0	284	fParticlesMatched->GetAxis(i)->SetTitle(title[i]);
427d6cf5	285	fParticlesMatched->SetBinEdges(i, binEdges[i]);
	286	}
	287
eec7bbb0	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");
eec7bbb0	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")
b63add30	397	contents[i] = (Double_t)trackType;
eec7bbb0	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 && type >= 0) {
	414	fTracksAll[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
	415
	416	Int_t label = TMath::Abs(track->GetLabel());
	417
d7ab1a93	418	if (fSelectHIJING && (label==0 \|\| track->GetGeneratorIndex() == 0)) {
d7ab1a93	419	// reject particles generated from other generators in the cocktail but keep fake tracks (label == 0)
eec7bbb0	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) {
d7ab1a93	426	if (!fSelectHIJING \|\| part->GetGeneratorIndex() == 0) {
eec7bbb0	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());
d7ab1a93	446	if (!fSelectHIJING \|\| part->GetGeneratorIndex() == 0) {
d7ab1a93	447	if (fSelectHIJING) fParticlesSelected[fCentBin]->Fill(part->Eta(), part->Phi(), part->Pt());
eec7bbb0	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	}