[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ParticleEfficiency / AliAnalysisTaskParticleEfficiency.cxx

#include "AliAnalysisTaskParticleEfficiency.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliInputEventHandler.h"
#include "AliAnalysisManager.h"
#include "AliESDVertex.h"
#include "AliCentrality.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "TParticle.h"
#include "TObjArray.h"
#include "TList.h"
#include "TH2F.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"

ClassImp(AliAnalysisTaskParticleEfficiency)

//_______________________________________________________

AliAnalysisTaskParticleEfficiency::AliAnalysisTaskParticleEfficiency(const Char_t *partName) :
  AliAnalysisTaskSE(partName),
  fParticlePdgCode(0),
  fTrackCuts(NULL)
{

  /* 
   * default constructor
   */

  /* set particle PDG code */
  TParticlePDG *ppdg = TDatabasePDG::Instance()->GetParticle(partName);
  if (ppdg) fParticlePdgCode = ppdg->PdgCode();

  /* init track cuts */
  fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
  fTrackCuts->SetEtaRange(-0.8, 0.8);
  fTrackCuts->SetPtRange(0.15, 1.e10);

  /* create output */
  fHistoList = new TList();
  fHistoList->SetOwner(kTRUE);
  DefineOutput(1, TList::Class());
}

//_______________________________________________________

AliAnalysisTaskParticleEfficiency::~AliAnalysisTaskParticleEfficiency()
{

  /* 
   * default destructor
   */
  
}

//_______________________________________________________

void
AliAnalysisTaskParticleEfficiency::UserCreateOutputObjects()
{

  /* 
   * user create output objects
   */

  fHistoEvents = new TH2F("hHistoEvents", ";centrality percentile;track multiplicity", 20, 0., 100., 200, 0., 2000.);
  fHistoList->Add(fHistoEvents);

  fHistoGenerated = new TH2F("hHistoGenerated", ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
  fHistoList->Add(fHistoGenerated);
  
  fHistoReconstructed = new TH2F("hHistoReconstructed", ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
  fHistoList->Add(fHistoReconstructed);
  
  for (Int_t idau = 0; idau < 2; idau++) {
    
    fHistoGeneratedDaughter[idau] = new TH2F(Form("hHistoGeneratedDaughter_%d", idau), ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
    fHistoList->Add(fHistoGeneratedDaughter[idau]);
    
    fHistoReconstructedDaughter[idau] = new TH2F(Form("hHistoReconstructedDaughter_%d", idau), ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
    fHistoList->Add(fHistoReconstructedDaughter[idau]);
    
  }  

  PostData(1, fHistoList);
}

//_______________________________________________________

void
AliAnalysisTaskParticleEfficiency::UserExec(Option_t *option)
{

  /* 
   * user exec
   */

  /*** INIT ***/

  /* get ESD event */
  AliESDEvent *esdEvent = dynamic_cast<AliESDEvent *>(InputEvent());
  if (!esdEvent) return;
  /* get MC event */
  AliMCEvent *mcEvent = dynamic_cast<AliMCEvent *>(MCEvent());
  if (!mcEvent) return;
  /* get stack */
  AliStack *mcStack = mcEvent->Stack();
  if (!mcStack) return;

  /*** EVENT SELECTION ***/

  /* collision candidate */
  if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
  /* vertex selection */
  const AliESDVertex *vertex = esdEvent->GetPrimaryVertexTracks();
  if (vertex->GetNContributors() < 1) {
    vertex = esdEvent->GetPrimaryVertexSPD();
    if (vertex->GetNContributors() < 1) return;
  }
  /* z-vertex cut */
  if (TMath::Abs(vertex->GetZ()) > 10.) return;
  /* centrality selection */
  AliCentrality *centrality = esdEvent->GetCentrality();
  if (centrality->GetQuality() != 0) return;
  Double_t cent = centrality->GetCentralityPercentileUnchecked("V0M");

  /* fill event histos */
  fHistoEvents->Fill(cent, fTrackCuts->CountAcceptedTracks(esdEvent));

  /*** RECONSTRUCTED TRACKS ***/

  TObjArray recoParticleArray;

  /* loop over ESD tracks */
  for (Int_t itrk = 0; itrk < esdEvent->GetNumberOfTracks(); itrk++) {

    /* get track */
    AliESDtrack *track = esdEvent->GetTrack(itrk);
    if (!track) continue;
    /* check accept track */
    if (!fTrackCuts->AcceptTrack(track)) continue;

    /* get coresponding MC particle */
    TParticle *particle = mcStack->Particle(TMath::Abs(track->GetLabel()));
    if (!particle) continue;

    /* add to reconstructed particle array */
    recoParticleArray.Add(particle);
  }

  /*** MONTECARLO PARTICLES ***/

  /* loop over MC stack */
  for (Int_t ipart = 0; ipart < mcStack->GetNtrack(); ipart++) {
    /* get particle */
    TParticle *particle = mcStack->Particle(ipart);
    if (!particle) continue;
    /* check PDG */
    if (particle->GetPdgCode() != fParticlePdgCode) continue;
    /* check rapidity */
    if (TMath::Abs(particle->Y()) > 0.5) continue;

    /* switch PDG */
    switch (fParticlePdgCode) {
      
    case 211: case -211: /* pions */
    case 321: case -321: /* kaons */
    case 2212: case -2212: /* protons */

      /* check physical primary */
      if (!mcStack->IsPhysicalPrimary(ipart)) continue;

      /* fill particle histos */
      fHistoGenerated->Fill(cent, particle->Pt());
      if (recoParticleArray.Contains(particle))
	fHistoReconstructed->Fill(cent, particle->Pt());
      break;

    case 333: /* phi */
    case 313: /* K*0 */
    
      /* check number of daughters */
      if (particle->GetNDaughters() != 2) continue;
      /* check daughter charge */
      TParticle *daughter0 =  mcStack->Particle(particle->GetDaughter(0));
      TParticle *daughter1 =  mcStack->Particle(particle->GetDaughter(1));
      if (TDatabasePDG::Instance()->GetParticle(daughter0->GetPdgCode())->Charge() == 0.) continue;
      if (TDatabasePDG::Instance()->GetParticle(daughter1->GetPdgCode())->Charge() == 0.) continue;

      /* fill daughter histos */
      fHistoGeneratedDaughter[0]->Fill(cent, daughter0->Pt());
      if (recoParticleArray.Contains(daughter0))
	fHistoReconstructedDaughter[0]->Fill(cent, daughter0->Pt());
      fHistoGeneratedDaughter[1]->Fill(cent, daughter1->Pt());
      if (recoParticleArray.Contains(daughter1))
	fHistoReconstructedDaughter[1]->Fill(cent, daughter1->Pt());

      /* fill particle histos */
      fHistoGenerated->Fill(cent, particle->Pt());
      if (recoParticleArray.Contains(daughter0) &&
	  recoParticleArray.Contains(daughter1)) 
	fHistoReconstructed->Fill(cent, particle->Pt());
      break;
    }
    
  }

  PostData(1, fHistoList);
}
Commit	Line	Data
52961962	1	#include "AliAnalysisTaskParticleEfficiency.h"
	2	#include "AliESDEvent.h"
	3	#include "AliMCEvent.h"
	4	#include "AliStack.h"
	5	#include "AliInputEventHandler.h"
	6	#include "AliAnalysisManager.h"
	7	#include "AliESDVertex.h"
	8	#include "AliCentrality.h"
	9	#include "AliESDtrack.h"
	10	#include "AliESDtrackCuts.h"
	11	#include "TParticle.h"
	12	#include "TObjArray.h"
	13	#include "TList.h"
	14	#include "TH2F.h"
	15	#include "TDatabasePDG.h"
	16	#include "TParticlePDG.h"
	17
	18	ClassImp(AliAnalysisTaskParticleEfficiency)
	19
	20	//_______________________________________________________
	21
	22	AliAnalysisTaskParticleEfficiency::AliAnalysisTaskParticleEfficiency(const Char_t *partName) :
	23	AliAnalysisTaskSE(partName),
	24	fParticlePdgCode(0),
	25	fTrackCuts(NULL)
	26	{
	27
	28	/*
	29	* default constructor
	30	*/
	31
	32	/* set particle PDG code */
	33	TParticlePDG *ppdg = TDatabasePDG::Instance()->GetParticle(partName);
	34	if (ppdg) fParticlePdgCode = ppdg->PdgCode();
	35
	36	/* init track cuts */
	37	fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
	38	fTrackCuts->SetEtaRange(-0.8, 0.8);
	39	fTrackCuts->SetPtRange(0.15, 1.e10);
	40
	41	/* create output */
	42	fHistoList = new TList();
	43	fHistoList->SetOwner(kTRUE);
	44	DefineOutput(1, TList::Class());
	45	}
	46
	47	//_______________________________________________________
	48
	49	AliAnalysisTaskParticleEfficiency::~AliAnalysisTaskParticleEfficiency()
	50	{
	51
	52	/*
	53	* default destructor
	54	*/
	55
	56	}
	57
	58	//_______________________________________________________
	59
	60	void
	61	AliAnalysisTaskParticleEfficiency::UserCreateOutputObjects()
	62	{
	63
	64	/*
65	* user create output objects
66	*/
67
68	fHistoEvents = new TH2F("hHistoEvents", ";centrality percentile;track multiplicity", 20, 0., 100., 200, 0., 2000.);
69	fHistoList->Add(fHistoEvents);
70
71	fHistoGenerated = new TH2F("hHistoGenerated", ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
72	fHistoList->Add(fHistoGenerated);
73
74	fHistoReconstructed = new TH2F("hHistoReconstructed", ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
75	fHistoList->Add(fHistoReconstructed);
76
77	for (Int_t idau = 0; idau < 2; idau++) {
78
79	fHistoGeneratedDaughter[idau] = new TH2F(Form("hHistoGeneratedDaughter_%d", idau), ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
80	fHistoList->Add(fHistoGeneratedDaughter[idau]);
81
82	fHistoReconstructedDaughter[idau] = new TH2F(Form("hHistoReconstructedDaughter_%d", idau), ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
83	fHistoList->Add(fHistoReconstructedDaughter[idau]);
84
85	}
86
87	PostData(1, fHistoList);
88	}
89
90	//_______________________________________________________
91
92	void
93	AliAnalysisTaskParticleEfficiency::UserExec(Option_t *option)
94	{
95
96	/*
97	* user exec
98	*/
99
100	/* INIT */
101
102	/* get ESD event */
103	AliESDEvent esdEvent = dynamic_cast<AliESDEvent >(InputEvent());
104	if (!esdEvent) return;
105	/* get MC event */
106	AliMCEvent mcEvent = dynamic_cast<AliMCEvent >(MCEvent());
107	if (!mcEvent) return;
108	/* get stack */
109	AliStack *mcStack = mcEvent->Stack();
110	if (!mcStack) return;
111
112	/* EVENT SELECTION */
113
114	/* collision candidate */
115	if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
116	/* vertex selection */
117	const AliESDVertex *vertex = esdEvent->GetPrimaryVertexTracks();
118	if (vertex->GetNContributors() < 1) {
119	vertex = esdEvent->GetPrimaryVertexSPD();
120	if (vertex->GetNContributors() < 1) return;
121	}
122	/* z-vertex cut */
123	if (TMath::Abs(vertex->GetZ()) > 10.) return;
124	/* centrality selection */
125	AliCentrality *centrality = esdEvent->GetCentrality();
126	if (centrality->GetQuality() != 0) return;
127	Double_t cent = centrality->GetCentralityPercentileUnchecked("V0M");
128
129	/* fill event histos */
130	fHistoEvents->Fill(cent, fTrackCuts->CountAcceptedTracks(esdEvent));
131
132	/* RECONSTRUCTED TRACKS */
133
134	TObjArray recoParticleArray;
135
136	/* loop over ESD tracks */
137	for (Int_t itrk = 0; itrk < esdEvent->GetNumberOfTracks(); itrk++) {
138
139	/* get track */
140	AliESDtrack *track = esdEvent->GetTrack(itrk);
141	if (!track) continue;
142	/* check accept track */
143	if (!fTrackCuts->AcceptTrack(track)) continue;
144
145	/* get coresponding MC particle */
146	TParticle *particle = mcStack->Particle(TMath::Abs(track->GetLabel()));
147	if (!particle) continue;
148
149	/* add to reconstructed particle array */
150	recoParticleArray.Add(particle);
151	}
152
153	/* MONTECARLO PARTICLES */
154
155	/* loop over MC stack */
156	for (Int_t ipart = 0; ipart < mcStack->GetNtrack(); ipart++) {
157	/* get particle */
158	TParticle *particle = mcStack->Particle(ipart);
159	if (!particle) continue;
160	/* check PDG */
161	if (particle->GetPdgCode() != fParticlePdgCode) continue;
162	/* check rapidity */
163	if (TMath::Abs(particle->Y()) > 0.5) continue;
164
165	/* switch PDG */
166	switch (fParticlePdgCode) {
167
168	case 211: case -211: /* pions */
169	case 321: case -321: /* kaons */
170	case 2212: case -2212: /* protons */
171
172	/* check physical primary */
173	if (!mcStack->IsPhysicalPrimary(ipart)) continue;
174
175	/* fill particle histos */
176	fHistoGenerated->Fill(cent, particle->Pt());
177	if (recoParticleArray.Contains(particle))
178	fHistoReconstructed->Fill(cent, particle->Pt());
179	break;
180
181	case 333: /* phi */
182	case 313: /* K0 /
183
184	/* check number of daughters */
185	if (particle->GetNDaughters() != 2) continue;
186	/* check daughter charge */
187	TParticle *daughter0 = mcStack->Particle(particle->GetDaughter(0));
188	TParticle *daughter1 = mcStack->Particle(particle->GetDaughter(1));
189	if (TDatabasePDG::Instance()->GetParticle(daughter0->GetPdgCode())->Charge() == 0.) continue;
190	if (TDatabasePDG::Instance()->GetParticle(daughter1->GetPdgCode())->Charge() == 0.) continue;
191
192	/* fill daughter histos */
193	fHistoGeneratedDaughter[0]->Fill(cent, daughter0->Pt());
194	if (recoParticleArray.Contains(daughter0))
195	fHistoReconstructedDaughter[0]->Fill(cent, daughter0->Pt());
196	fHistoGeneratedDaughter[1]->Fill(cent, daughter1->Pt());
197	if (recoParticleArray.Contains(daughter1))
198	fHistoReconstructedDaughter[1]->Fill(cent, daughter1->Pt());
199
200	/* fill particle histos */
201	fHistoGenerated->Fill(cent, particle->Pt());
202	if (recoParticleArray.Contains(daughter0) &&
203	recoParticleArray.Contains(daughter1))
204	fHistoReconstructed->Fill(cent, particle->Pt());
205	break;
206	}
207
208	}
209
210	PostData(1, fHistoList);
211	}
212