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