]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/comparison/AliITSComparisonTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
The present commit corresponds to an important change in the way the
[u/mrichter/AliRoot.git] / PWG1 / comparison / AliITSComparisonTask.cxx
CommitLineData
f2247c19 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 "AliESDtrack.h"
22#include "AliTrackReference.h"
23#include "AliMCComparisonTrack.h"
24
25#include "AliITSComparisonTask.h"
26
27ClassImp(AliITSComparisonTask)
28
29AliITSComparisonTask::AliITSComparisonTask()
30 : AliAnalysisTaskSE("AliITSComaprisonTask"),
31 fListOfHistos(0),
32 hgood(0),
33 hfound(0),
34 hfake(0),
35 hp(0),
36 hl(0),
37 hpt(0),
38 htip(0),
39 he(0),
40 hep(0),
41 hgoodPhi(0),
42 hfoundPhi(0),
43 hlip(0)
44{
45 // Default constructor
46 AliInfo("Default constructor AliITSComparisonTask");
47 // Define input and output slots here
48 // Input slot #0 works with a TChain
49 DefineInput(0, TChain::Class());
50 // Output slot #1 TList
51 DefineOutput(1, TList::Class());
52}
53
54
55
56AliITSComparisonTask::AliITSComparisonTask(const char* name)
57 : AliAnalysisTaskSE(name),
58 fListOfHistos(0),
59 hgood(0),
60 hfound(0),
61 hfake(0),
62 hp(0),
63 hl(0),
64 hpt(0),
65 htip(0),
66 he(0),
67 hep(0),
68 hgoodPhi(0),
69 hfoundPhi(0),
70 hlip(0)
71{
72 // Constructor
73 AliInfo("Constructor AliITSComparisonTask");
74 // Define input and output slots here
75 // Input slot #0 works with a TChain
76 DefineInput(0, TChain::Class());
77 // Output slot #1 TList
78 DefineOutput(1, TList::Class());
79}
80
81
82
83void AliITSComparisonTask::UserCreateOutputObjects()
84{
85 // Create histograms
86 // Called once
87 AliInfo("AliITSComparisonTask::UserCreateOutputObjects");
88 // Create output container
89 fListOfHistos = new TList();
90
91 hgood = new TH1F("hgood", "Pt for good tracks", 34, 0.2, 7.0);
92 hfound = new TH1F("hfound", "Pt for found tracks", 34, 0.2, 7.0);
93 hfake = new TH1F("hfake", "Pt for fake tracks", 34, 0.2, 7.0);
94 hp = new TH1F("hp", "PHI resolution", 50, -20., 20.);
95 hl = new TH1F("hl", "LAMBDA resolution", 50, -20., 20.);
96 hpt = new TH1F("hpt", "Relative Pt resolution", 30, -10., 10.);
97 htip = new TH1F("htip", "Transverse impact parameter", 40, -800., 800.);
98 he = new TH1F("he", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 100.);
99 hep = new TH2F("hep", "dE/dX vs momentum", 50, 0., 2., 50, 0., 400.);
100 hgoodPhi = new TH1F("hgoodPhi", "Phi for good tracks", 90, 0., 2.*TMath::Pi());
101 hfoundPhi = new TH1F("hfoundPhi", "Phi for found tracks", 90, 0., 2.*TMath::Pi());
102 hlip = new TH1F("hlip","Longitudinal impact parameter", 40, -800., 800.);
103
104
105 fListOfHistos->Add(hgood);
106 fListOfHistos->Add(hfound);
107 fListOfHistos->Add(hfake);
108 fListOfHistos->Add(hp);
109 fListOfHistos->Add(hl);
110 fListOfHistos->Add(hpt);
111 fListOfHistos->Add(htip);
112 fListOfHistos->Add(he);
113 fListOfHistos->Add(hep);
114 fListOfHistos->Add(hgoodPhi);
115 fListOfHistos->Add(hfoundPhi);
116 fListOfHistos->Add(hlip);
117}
118
119
120void AliITSComparisonTask::UserExec(Option_t *)
121{
122 // Main loop
123 // Called for each event
124
125 // MC information
126 AliMCEvent* mcEvent = MCEvent();
127 if (!mcEvent) {
128 Printf("ERROR: Could not retrieve MC event");
129 return;
130 }
131
132 Int_t nt = 0;
133
134 TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
135
136 // Loop over all MC tracks and select "good" tracks
137 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
138 AliMCParticle* track = mcEvent->GetTrack(iTracks);
139 if (!track) {
140 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
141 continue;
142 }
143
144 // Track selection
145 if (track->Pt() < 0.2) continue;
146 if (track->Pt() > 7.0) continue;
147 if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
148
149 Double_t vx = track->Xv(), vy = track->Yv(), vz = track->Zv();
150 if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
151 if (TMath::Abs(vz) > 50.) continue;
152
153 // Loop over Track References
154 Bool_t LabelTPC = kFALSE;
155 AliTrackReference* trackRef = 0;
156 UInt_t ITSLayerMap = 0;
157 for (Int_t iTrackRef = 0; iTrackRef < track->GetNumberOfTrackReferences(); iTrackRef++) {
158 trackRef = track->GetTrackReference(iTrackRef);
159 if(trackRef) {
160 Int_t detectorId = trackRef->DetectorId();
161 if (detectorId == AliTrackReference::kITS) {
162 Float_t Radius = trackRef->R();
163 if (Radius > 2.5 && Radius < 5.5)
164 ITSLayerMap |= 1;
165 else if (Radius > 5.5 && Radius < 8.5)
166 ITSLayerMap |= 1 << 1;
167 else if (Radius > 13.5 && Radius < 16.5)
168 ITSLayerMap |= 1 << 2;
169 else if (Radius > 22. && Radius < 26.)
170 ITSLayerMap |= 1 << 3;
171 else if (Radius > 36. && Radius < 41.)
172 ITSLayerMap |= 1 << 4;
173 else if (Radius > 42. && Radius < 46.)
174 ITSLayerMap |= 1 << 5;
175 else {
176 Printf("Wrong radius %f ", Radius);
177 return;
178 }
179
180 }
181 if (detectorId == AliTrackReference::kTPC) {
182 LabelTPC = kTRUE;
183 break;
184 }
185 }
186 } // track references loop
187
188 // Skip tracks that passed not all ITS layers
189 if (ITSLayerMap != 0x3F) continue;
190
191 // "Good" tracks
192 if (LabelTPC) {
193 AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
194 ref->SetLabel(iTracks);
195 TParticle* particle = track->Particle();
196 Int_t pdg = particle->GetPdgCode();
197 ref->SetPDG(pdg);
198 ref->SetPz(track->Pz());
199 Float_t Pt = track->Pt();
200 ref->SetPt(Pt);
201 hgood->Fill(Pt);
202 Float_t phig = track->Phi();
203 ref->SetPhi(phig);
204 hgoodPhi->Fill(phig);
205 nt++;
206 }
207 } // track loop
208
209
210 // ESD information
211 AliVEvent* event = InputEvent();
212 if (!event) {
213 Printf("ERROR: Could not retrieve event");
214 return;
215 }
216
217 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
218
219 Bool_t iFound;
220 Int_t nfound = 0;
221 Int_t nfake = 0;
222 Int_t nlost = 0;
223
224 Int_t MCgoods = refs->GetEntriesFast();
225
226 // Loop over all "good" MC tracks
227 for (Int_t k = 0; k < MCgoods; k++) {
228 AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
229 Int_t MCLabel = ref->GetLabel();
230 Float_t ptg = ref->GetPt();
231 Float_t PhiG = ref->GetPhi();
232 iFound = kFALSE;
233 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
234 AliESDtrack* track = esd->GetTrack(iTrack);
235 if (!track) {
236 Printf("ERROR: Could not receive track %d", iTrack);
237 continue;
238 }
239
240 // if (! track->IsOn(AliESDtrack::kTPCrefit)) continue;
241 if (! track->IsOn(AliESDtrack::kITSrefit)) continue;
242
243 Int_t Label = track->GetLabel();
244
245 if (MCLabel == TMath::Abs(Label)) {
246 if (MCLabel == Label) {
247 nfound++;
248 hfound->Fill(ptg);
249 hfoundPhi->Fill(PhiG);
250 }
251 else {
252 nfake++;
253 hfake->Fill(ptg);
254 }
255 iFound = kTRUE;
256
257 Float_t phi = track->Phi();
258 Double_t pt = track->Pt();
259
260 Float_t lam = TMath::ATan2(track->Pz(),pt);
261 Float_t pt_1 = 1/pt;
262
263 Float_t phig = ref->GetPhi();
264 hp->Fill((phi - phig)*1000.);
265
266 Float_t lamg = TMath::ATan2(ref->GetPz(),ptg);
267 hl->Fill((lam - lamg)*1000.);
268
269 hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
270
271 Float_t d,z; track->GetImpactParameters(d,z);
272 htip->Fill(10000.*d);
273 hlip->Fill(10000.*z);
274
275 Float_t mom = pt/TMath::Cos(lam);
276 Float_t dedx = track->GetTPCsignal();
277 hep->Fill(mom, dedx, 1.);
278
279 Int_t pdg = ref->GetPDG();
280 if (TMath::Abs(pdg )== 211) { // pions
281 if (mom > 0.4 && mom < 0.5) {
282 he->Fill(dedx, 1.);
283 }
284 }
285
286 break;
287 }
288 }
289 if (!iFound) {
290 nlost++;
291 }
292 }
293
294 Printf(" Results: " );
295 Printf(" Found %d Fake %d Lost %d ", nfound, nfake, nlost);
296
297 refs->Clear();
298
299 // Post output data.
300 PostData(1, fListOfHistos);
301}
302
303//________________________________________________________________________
304void AliITSComparisonTask::Terminate(Option_t *)
305{
306
307 // Draw result to the screen
308 // Called once at the end of the query
309 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
310 if (!fListOfHistos) {
311 Printf("ERROR: fListOfHistos not available");
312 return;
313 }
314
315 hgood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
316 hfound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
317 hfake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
318 hp = dynamic_cast<TH1F*>(fListOfHistos->At(3));
319 hl = dynamic_cast<TH1F*>(fListOfHistos->At(4));
320 hpt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
321 htip = dynamic_cast<TH1F*>(fListOfHistos->At(6));
322 he = dynamic_cast<TH1F*>(fListOfHistos->At(7));
323 hep = dynamic_cast<TH2F*>(fListOfHistos->At(8));
324 hgoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
325 hfoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
326 hlip = dynamic_cast<TH1F*>(fListOfHistos->At(11));
327
328 gStyle->SetOptStat(111110);
329 gStyle->SetOptFit(1);
330
331 TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
332
333 Int_t minc = 33;
334
335 TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
336 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
337 hp->SetFillColor(4); hp->SetXTitle("(mrad)");
338 if (hp->GetEntries() < minc) hp->Draw(); else hp->Fit("gaus"); c1->cd();
339
340 TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1., 0.6); p2->Draw();
341 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
342 hl->SetFillColor(4); hl->SetXTitle("(mrad)");
343 if (hl->GetEntries() < minc) hl->Draw(); else hl->Fit("gaus"); c1->cd();
344
345 TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
346 p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
347 hpt->SetFillColor(2); hpt->SetXTitle("(%)");
348 if (hpt->GetEntries() < minc) hpt->Draw(); else hpt->Fit("gaus"); c1->cd();
349
350 TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
351 p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
352 htip->SetFillColor(6); htip->SetXTitle("(micron)");
353 if (htip->GetEntries() < minc) htip->Draw(); else htip->Fit("gaus");
354 hlip->Draw("same"); c1->cd();
355
356 TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
357 p5->SetFillColor(41); p5->SetFrameFillColor(10);
358 hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
359 TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
360 TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
361 hg->Divide(hfound,hgood,1.,1.,"B");
362 hf->Divide(hfake,hgood,1.,1.,"B");
363 hg->SetLineColor(4); hg->SetLineWidth(2);
364 hg->SetMaximum(1.4);
365 hg->SetYTitle("Tracking efficiency");
366 hg->SetXTitle("Pt (GeV/c)");
367 hg->Draw();
368
369 TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
370 line1->Draw("same");
371 TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
372 line2->Draw("same");
373
374 hf->SetFillColor(1);
375 hf->SetFillStyle(3013);
376 hf->SetLineColor(2);
377 hf->SetLineWidth(2);
378 hf->Draw("histsame");
379 TText* text = new TText(0.461176, 0.248448, "Fake tracks");
380 text->SetTextSize(0.05);
381 text->Draw();
382 text = new TText(0.453919, 1.11408, "Good tracks");
383 text->SetTextSize(0.05);
384 text->Draw();
385
386 TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
387
388 TPad* p6 = new TPad("p6", "", 0., 0., 1., .5); p6->Draw();
389 p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
390 he->SetFillColor(2); he->SetFillStyle(3005);
391 he->SetXTitle("Arbitrary Units");
392 if (he->GetEntries() < minc) he->Draw(); else he->Fit("gaus"); c2->cd();
393
394 TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
395 p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
396 hep->SetFillColor(2); hep->SetFillStyle(3005);
397 hep->SetMarkerStyle(8);
398 hep->SetMarkerSize(0.4);
399 hep->SetXTitle("p (GeV/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
400 hep->Draw(); c1->cd();
401
402 TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
403 c3->SetFillColor(42); c3->SetFrameFillColor(10);
404 hfoundPhi->Sumw2(); hgoodPhi->Sumw2();
405 TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
406 90, 0., 2.*TMath::Pi());
407 hgphi->Divide(hfoundPhi, hgoodPhi, 1., 1., "B");
408 hgphi->SetYTitle("Tracking efficiency");
409 hgphi->SetXTitle("Phi (rad)");
410 hgphi->Draw();
411
412 TFile fc("AliITSComparison.root","RECREATE");
413 c1->Write();
414 c2->Write();
415 c3->Write();
416 fc.Close();
417}