Bug corrected
[u/mrichter/AliRoot.git] / PWG1 / comparison / AliTPCComparisonTask.cxx
CommitLineData
31e081ed 1#include "TChain.h"
2#include "TTree.h"
3#include "TH1F.h"
4#include "TH2F.h"
5#include "TList.h"
6#include "TClonesArray.h"
7#include "TCanvas.h"
8#include "TStyle.h"
9#include "TLine.h"
10#include "TText.h"
11#include "TFile.h"
12#include "TBenchmark.h"
13
14#include "AliAnalysisTask.h"
15#include "AliAnalysisManager.h"
16
17#include "AliLog.h"
18#include "AliVEvent.h"
19#include "AliESDEvent.h"
20#include "AliMCEvent.h"
21#include "AliAODEvent.h"
22
23#include "AliTrackReference.h"
24#include "AliMCComparisonTrack.h"
25
26#include "AliTPCComparisonTask.h"
27
28ClassImp(AliTPCComparisonTask)
29
30AliTPCComparisonTask::AliTPCComparisonTask(const char* name)
31 : AliAnalysisTaskSE(name), fListOfHistos(0), hgood(0), hfound(0),
32 hfake(0), hp(0), hl(0), hpt(0), hmpt(0), he(0), hep(0), hgoodPhi(0),
33 hfoundPhi(0)
34{
35 // Constructor
36 AliInfo("Constructor AliTPCComparisonTask");
37 // Define input and output slots here
38 // Input slot #0 works with a TChain
39 DefineInput(0, TChain::Class());
40 // Output slot #1 TList
41 DefineOutput(1, TList::Class());
42}
43
44
45
46void AliTPCComparisonTask::UserCreateOutputObjects()
47{
48 // Create histograms
49 // Called once
50 AliInfo("AliTPCComparisonTask::UserCreateOutputObjects");
51 // Create output container
52 fListOfHistos = new TList();
53
54 hgood = new TH1F("hgood", "Good tracks", 30, 0.2, 6.1);
55 hfound= new TH1F("hfound", "Found tracks", 30, 0.2, 6.1);
56 hfake = new TH1F("hfake", "Fake tracks", 30, 0.2, 6.1);
57 hp = new TH1F("hp", "PHI resolution", 50, -20., 20.);
58 hl = new TH1F("hl", "LAMBDA resolution", 50, -20., 20.);
59 hpt = new TH1F("hpt", "Relative Pt resolution", 30, -10., 10.);
60 hmpt = new TH1F("hmpt", "Relative Pt resolution (pt>4GeV/c)", 30, -60., 60.);
61 he = new TH1F("he", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 100.);
62 hep = new TH2F("hep", "dE/dX vs momentum", 50, 0., 2., 50, 0., 400.);
63 hgoodPhi = new TH1F("hgoodPhi", "Phi for Good tracks", 90, -TMath::Pi(), TMath::Pi());
64 hfoundPhi = new TH1F("hfoundPhi", "Phi for Found tracks", 90, -TMath::Pi(), TMath::Pi());
65
66 fListOfHistos->Add(hgood);
67 fListOfHistos->Add(hfound);
68 fListOfHistos->Add(hfake);
69 fListOfHistos->Add(hp);
70 fListOfHistos->Add(hl);
71 fListOfHistos->Add(hpt);
72 fListOfHistos->Add(hmpt);
73 fListOfHistos->Add(he);
74 fListOfHistos->Add(hep);
75 fListOfHistos->Add(hgoodPhi);
76 fListOfHistos->Add(hfoundPhi);
77}
78
79
80void AliTPCComparisonTask::UserExec(Option_t *)
81{
82 // Main loop
83 // Called for each event
84
85 // MC information
86 AliMCEvent* mcEvent = MCEvent();
87 if (!mcEvent) {
88 Printf("ERROR: Could not retrieve MC event");
89 return;
90 }
91
92 Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
93
94 Int_t nt = 0;
95
96 TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
97
98 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++)
99 {
100 AliMCParticle* track = mcEvent->GetTrack(iTracks);
101 if (!track)
102 {
103 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
104 continue;
105 }
106
107 // Track selection
108 if (track->Pt()<0.100) continue;
109 if (TMath::Abs(track->Pz()/track->Pt())>0.999) continue;
110
111 Double_t vx=track->Xv(),vy=track->Yv(),vz=track->Zv();
112 if (TMath::Sqrt(vx*vx+vy*vy)>3.5) continue;
113 if (TMath::Abs(vz) > 50.) continue;
114
115 Bool_t LabelTPC = kFALSE;
116 AliTrackReference* trackRef = 0;
117 for (Int_t iTrackRef = 0; iTrackRef < track->GetNumberOfTrackReferences(); iTrackRef++)
118 {
119 trackRef = track->GetTrackReference(iTrackRef);
120 if(trackRef)
121 {
122 Int_t detectorId = trackRef->DetectorId();
123 if (detectorId == AliTrackReference::kTPC)
124 {
125 LabelTPC = kTRUE;
126 break;
127 }
128 }
129 }
130 if (LabelTPC)
131 {
132 AliMCComparisonTrack *ref=new((*refs)[nt]) AliMCComparisonTrack();
133 ref->SetLabel(iTracks);
134 TParticle* particle= track->Particle();
135 Int_t pdg = particle->GetPdgCode();
136 ref->SetPDG(pdg);
137 ref->SetTPCMomentum(trackRef->Px(),trackRef->Py(),trackRef->Pz());
138 Float_t Pt = TMath::Sqrt(trackRef->Px()*trackRef->Px() +
139 trackRef->Py()*trackRef->Py());
140 hgood->Fill(Pt);
141 Float_t phig=TMath::ATan2(trackRef->Py(),trackRef->Px());
142 hgoodPhi->Fill(phig);
143 nt++;
144 }
145 } //track loop
146
147
148 // ESD information
149 AliVEvent *event = InputEvent();
150 if (!event) {
151 Printf("ERROR: Could not retrieve event");
152 return;
153 }
154
155
156 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
157
158 Bool_t iFound;
159 Int_t nfound = 0;
160 Int_t nfake = 0;
161 Int_t nlost = 0;
162
163 Int_t MCgoods = refs->GetEntriesFast();
164
165 for (Int_t k = 0; k < MCgoods; k++)
166 {
167 AliMCComparisonTrack *ref=(AliMCComparisonTrack*)refs->UncheckedAt(k);
168 Int_t MCLabel = ref->GetLabel();
169 Float_t ptg =
170 TMath::Sqrt(ref->GetTPCPx()*ref->GetTPCPx() +
171 ref->GetTPCPy()*ref->GetTPCPy());
172 Float_t PhiG=TMath::ATan2(ref->GetTPCPy(),ref->GetTPCPx());
173
174 iFound = kFALSE;
175 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++)
176 {
177 AliESDtrack* track = esd->GetTrack(iTrack);
178 if (!track)
179 {
180 Printf("ERROR: Could not receive track %d", iTrack);
181 continue;
182 }
183 Int_t TPCLabel = track->GetTPCLabel();
184
185 if (MCLabel == TMath::Abs(TPCLabel))
186 {
187 if (MCLabel == TPCLabel)
188 {
189 nfound++;
190 hfound->Fill(ptg);
191 hfoundPhi->Fill(PhiG);
192 }
193 else
194 {
195 nfake++;
196 hfake->Fill(ptg);
197 }
198 iFound = kTRUE;
199
200
201 Double_t pxpypz[3]; track->GetInnerPxPyPz(pxpypz);
202 Float_t phi=TMath::ATan2(pxpypz[1],pxpypz[0]);
203 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
204 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
205 Double_t pt=TMath::Sqrt(pxpypz[0]*pxpypz[0]+pxpypz[1]*pxpypz[1]);
206 Float_t lam=TMath::ATan2(pxpypz[2],pt);
207 Float_t pt_1=1/pt;
208
209 Int_t pdg=(Int_t)ref->GetPDG();
210 if (TMath::Abs(pdg)==11 && ptg>4.)
211 {
212 //high momentum electrons
213 hmpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
214 }
215 else
216 {
217 Float_t phig=TMath::ATan2(ref->GetTPCPy(),ref->GetTPCPx());
218 hp->Fill((phi - phig)*1000.);
219 Float_t lamg=TMath::ATan2(ref->GetTPCPz(),ptg);
220 hl->Fill((lam - lamg)*1000.);
221 hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
222 }
223
224 Float_t mom=pt/TMath::Cos(lam);
225 Float_t dedx=track->GetTPCsignal();
226 hep->Fill(mom,dedx,1.);
227 if (TMath::Abs(pdg)==211) //pions
228 if (mom>0.4 && mom<0.5)
229 {
230 he->Fill(dedx,1.);
231 }
232
233
234 break;
235 }
236 }
237 if (!iFound)
238 {
239 nlost++;
240 }
241 }
242
243 Printf(" Results: " );
244 Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
245
246 refs->Clear();
247
248 // Post output data.
249 PostData(1, fListOfHistos);
250}
251
252//________________________________________________________________________
253void AliTPCComparisonTask::Terminate(Option_t *)
254{
255 // Draw result to the screen
256 // Called once at the end of the query
257
258 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
259 if (!fListOfHistos) {
260 Printf("ERROR: fListOfHistos not available");
261 return;
262 }
263
264
265 hgood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
266 hfound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
267 hfake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
268 hp = dynamic_cast<TH1F*>(fListOfHistos->At(3));
269 hl = dynamic_cast<TH1F*>(fListOfHistos->At(4));
270 hpt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
271 hmpt = dynamic_cast<TH1F*>(fListOfHistos->At(6));
272 he = dynamic_cast<TH1F*>(fListOfHistos->At(7));
273 hep = dynamic_cast<TH2F*>(fListOfHistos->At(8));
274 hgoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
275 hfoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
276
277
278 gStyle->SetOptStat(111110);
279 gStyle->SetOptFit(1);
280
281 TCanvas* c1=new TCanvas("c1","",0,0,700,850);
282
283 Int_t minc=33;
284
285 TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
286 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
287 hp->SetFillColor(4); hp->SetXTitle("(mrad)");
288 if (hp->GetEntries()<minc) hp->Draw(); else hp->Fit("gaus"); c1->cd();
289
290 TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
291 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
292 hl->SetXTitle("(mrad)");
293 if (hl->GetEntries()<minc) hl->Draw(); else hl->Fit("gaus"); c1->cd();
294
295 TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
296 p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
297 hpt->SetXTitle("(%)");
298 if (hpt->GetEntries()<minc) hpt->Draw(); else hpt->Fit("gaus"); c1->cd();
299
300 TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
301 p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
302 hmpt->SetXTitle("(%)");
303 if (hmpt->GetEntries()<minc) hmpt->Draw(); else hmpt->Fit("gaus"); c1->cd();
304
305
306 TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
307 p5->SetFillColor(41); p5->SetFrameFillColor(10);
308 hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
309 TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 30, 0.2, 6.1);
310 TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 30, 0.2, 6.1);
311 hg->Divide(hfound,hgood,1.,1.,"B");
312 hf->Divide(hfake,hgood,1.,1.,"B");
313 hg->SetMaximum(1.4);
314 hg->SetYTitle("Tracking efficiency");
315 hg->SetXTitle("Pt (GeV/c)");
316 hg->Draw();
317
318 TLine *line1 = new TLine(0.1,1.0,6.1,1.0); line1->SetLineStyle(4);
319 line1->Draw("same");
320 TLine *line2 = new TLine(0.1,0.9,6.1,0.9); line2->SetLineStyle(4);
321 line2->Draw("same");
322
323 hf->SetFillColor(1);
324 hf->SetFillStyle(3013);
325 hf->SetLineColor(2);
326 hf->SetLineWidth(2);
327 hf->Draw("histsame");
328 TText *text = new TText(0.461176,0.248448,"Fake tracks");
329 text->SetTextSize(0.05);
330 text->Draw();
331 text = new TText(0.453919,1.11408,"Good tracks");
332 text->SetTextSize(0.05);
333 text->Draw();
334
335 TCanvas *c2=new TCanvas("c2","",320,32,530,590);
336
337 TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw();
338 p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
339 he->SetFillColor(2); he->SetFillStyle(3005);
340 he->SetXTitle("Arbitrary Units");
341 if (he->GetEntries()<minc) he->Draw(); else he->Fit("gaus"); c2->cd();
342
343 TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw();
344 p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
345 hep->SetFillColor(2); hep->SetFillStyle(3005);
346 hep->SetXTitle("p (GeV/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
347 hep->Draw(); c1->cd();
348
349 TCanvas* c3=new TCanvas("c3","",10,10,510,510);
350 c3->SetFillColor(42); c3->SetFrameFillColor(10);
351 hfoundPhi->Sumw2(); hgoodPhi->Sumw2();
352 TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
353 90, -TMath::Pi(), TMath::Pi());
354 hgphi->Divide(hfoundPhi,hgoodPhi,1.,1.,"B");
355 hgphi->SetYTitle("Tracking efficiency");
356 hgphi->SetXTitle("Phi (rad)");
357 hgphi->Draw();
358
359 TFile fc("AliTPCComparison.root","RECREATE");
360 c1->Write();
361 c2->Write();
362 c3->Write();
363 fc.Close();
364
365}