particle-efficiency task added to spectra library
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ParticleEfficiency / AliAnalysisTaskParticleEfficiency.cxx
CommitLineData
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
18ClassImp(AliAnalysisTaskParticleEfficiency)
19
20//_______________________________________________________
21
22AliAnalysisTaskParticleEfficiency::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
49AliAnalysisTaskParticleEfficiency::~AliAnalysisTaskParticleEfficiency()
50{
51
52 /*
53 * default destructor
54 */
55
56}
57
58//_______________________________________________________
59
60void
61AliAnalysisTaskParticleEfficiency::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
92void
16e017a0 93AliAnalysisTaskParticleEfficiency::UserExec(Option_t *)
52961962 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: /* K*0 */
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