fd14328b8bdc74b8d5c566b66107fd77c4eebc0f
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ParticleEfficiency / AliAnalysisTaskParticleEfficiency.cxx
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 *)
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