]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/comparison/AliTRDComparisonTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
better log validation
[u/mrichter/AliRoot.git] / PWGPP / comparison / AliTRDComparisonTask.cxx
CommitLineData
2d35f5d9 1//-------------------------------------------------------------------------
2//
3// This is the PROOF-enabled version of TRD/Macros/AliTRDComparisonV2.C macro.
4// Origin: Andrei.Zalite@cern.ch
5//
6//-------------------------------------------------------------------------
7
f2247c19 8#include "TChain.h"
9#include "TTree.h"
10#include "TH1F.h"
11#include "TH2F.h"
12#include "TList.h"
13#include "TClonesArray.h"
14#include "TCanvas.h"
15#include "TStyle.h"
16#include "TLine.h"
17#include "TText.h"
18#include "TFile.h"
f2247c19 19
20#include "AliLog.h"
21#include "AliVEvent.h"
22#include "AliESDEvent.h"
23#include "AliMCEvent.h"
24#include "AliESDtrack.h"
25#include "AliTrackReference.h"
26#include "AliMCComparisonTrack.h"
27
28#include "AliTRDComparisonTask.h"
29
30ClassImp(AliTRDComparisonTask)
31
2d35f5d9 32extern TStyle *gStyle;
33
f2247c19 34AliTRDComparisonTask::AliTRDComparisonTask()
35 : AliAnalysisTaskSE("AliTRDComparisonTask"),
36 fListOfHistos(0),
2d35f5d9 37 fGood(0),
38 fFound(0),
39 fFake(0),
40 fP(0),
41 fL(0),
42 fPt(0),
43 fHmpt(0),
44 fE(0),
45 fEp(0),
46 fGoodPhi(0),
47 fFoundPhi(0),
48 fZ(0),
49 fC(0)
f2247c19 50{
51 // Default constructor
52 AliInfo("Default constructor AliTRDComparisonTask");
53 // Define input and output slots here
54 // Input slot #0 works with a TChain
55 DefineInput(0, TChain::Class());
56 // Output slot #1 TList
57 DefineOutput(1, TList::Class());
58}
59
60AliTRDComparisonTask::AliTRDComparisonTask(const char* name)
61 : AliAnalysisTaskSE(name),
62 fListOfHistos(0),
2d35f5d9 63 fGood(0),
64 fFound(0),
65 fFake(0),
66 fP(0),
67 fL(0),
68 fPt(0),
69 fHmpt(0),
70 fE(0),
71 fEp(0),
72 fGoodPhi(0),
73 fFoundPhi(0),
74 fZ(0),
75 fC(0)
f2247c19 76{
77 // Constructor
78 AliInfo("Constructor AliTRDComparisonTask");
79 // Define input and output slots here
80 // Input slot #0 works with a TChain
81 DefineInput(0, TChain::Class());
82 // Output slot #1 TList
83 DefineOutput(1, TList::Class());
84}
85
86
87
88void AliTRDComparisonTask::UserCreateOutputObjects()
89{
90 // Create histograms
91 // Called once
92 AliInfo("AliTRDComparisonTask::UserCreateOutputObjects");
93 // Create output container
94 fListOfHistos = new TList();
95
2d35f5d9 96 fGood = new TH1F("fGood", "Pt for good tracks", 34, 0.2, 7.0);
97 fFound = new TH1F("fFound", "Pt for found tracks", 34, 0.2, 7.0);
98 fFake = new TH1F("fFake", "Pt for fake tracks", 34, 0.2, 7.0);
99 fP = new TH1F("fP", "PHI resolution", 50, -20., 20.);
100 fL = new TH1F("fL", "LAMBDA resolution", 50, -20., 20.);
101 fPt = new TH1F("fPt", "Relative Pt resolution", 30, -10., 10.);
102 fHmpt = new TH1F("fHmpt", "Y and Z resolution", 30, -30., 30.);
103 fE = new TH1F("fE", "dE/dX for pions with 0.4<p<0.5 GeV/c", 50, 0., 1000.);
104 fEp = new TH2F("fEp", "dE/dX vs momentum", 50, 0., 2., 50, 0., 2000.);
105 fGoodPhi = new TH1F("fGoodPhi", "Phi for Good tracks", 90, 0., 2.*TMath::Pi());
106 fFoundPhi = new TH1F("fFoundPhi", "Phi for Found tracks", 90, 0., 2.*TMath::Pi());
107 fZ = new TH1F("fZ", "Z resolution", 30, -30., 30.);
108 fC = new TH1F("fC", "Number of the assigned clusters", 25, 110., 135.);
f2247c19 109
2d35f5d9 110 fListOfHistos->Add(fGood);
111 fListOfHistos->Add(fFound);
112 fListOfHistos->Add(fFake);
113 fListOfHistos->Add(fP);
114 fListOfHistos->Add(fL);
115 fListOfHistos->Add(fPt);
116 fListOfHistos->Add(fHmpt);
117 fListOfHistos->Add(fE);
118 fListOfHistos->Add(fEp);
119 fListOfHistos->Add(fGoodPhi);
120 fListOfHistos->Add(fFoundPhi);
121 fListOfHistos->Add(fZ);
122 fListOfHistos->Add(fC);
f2247c19 123}
124
125
126void AliTRDComparisonTask::UserExec(Option_t *)
127{
128
129 // Main loop
130 // Called for each event
131
132 // MC information
133 AliMCEvent* mcEvent = MCEvent();
134 if (!mcEvent) {
135 Printf("ERROR: Could not retrieve MC event");
136 return;
137 }
138
139 Int_t nt = 0;
140
141 TClonesArray dummy("AliMCComparisonTrack",1000), *refs=&dummy;
142
143 // Loop over all MC tracks and select "good" tracks
144 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
84933864 145 AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
f2247c19 146 if (!track) {
147 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
148 continue;
149 }
150
151 // Track selection
152 if (track->Pt() < 0.2) continue;
153 if (track->Pt() > 7.0) continue;
154 if (TMath::Abs(track->Pz()/track->Pt()) > 0.999) continue;
155
156 Double_t vx = track->Xv(), vy = track->Yv(), vz = track->Zv();
157 if (TMath::Sqrt(vx*vx+vy*vy) > 3.5) continue;
158 if (TMath::Abs(vz) > 50.) continue;
159
160 // Check TRD information
161 Int_t nRefs = track->GetNumberOfTrackReferences();
162 if (nRefs <= 1) continue;
163
164 AliTrackReference* trackRef0 = 0;
165 for (Int_t iTrackRef = 0; iTrackRef < nRefs; iTrackRef++) {
166 trackRef0 = track->GetTrackReference(iTrackRef);
167 if(trackRef0) {
168 Int_t detectorId = trackRef0->DetectorId();
169 if (detectorId == AliTrackReference::kTRD) {
170 break;
171 }
172 trackRef0 = 0;
173 }
174 }
175 if (!trackRef0) continue;
176
177 if (trackRef0->LocalX() > 300.) continue;
178
179 AliTrackReference* trackRef = 0;
180 for (Int_t iTrackRef = nRefs - 1; iTrackRef >= 0; --iTrackRef) {
181 trackRef = track->GetTrackReference(iTrackRef);
182 if (!trackRef) continue;
183 if (trackRef->LocalX() > 363. &&
184 trackRef->DetectorId() == AliTrackReference::kTRD) break;
185 trackRef = 0;
186 }
187 if (!trackRef) continue;
188
189 if (TMath::Abs(trackRef0->Alpha() - trackRef->Alpha()) > 1e-5) continue;
190
191 // Check TPC information
2d35f5d9 192 Bool_t labelTPC = kFALSE;
f2247c19 193 for (Int_t iTrackRef = 0; iTrackRef < nRefs; iTrackRef++) {
194 trackRef = track->GetTrackReference(iTrackRef);
195 if(trackRef) {
196 Int_t detectorId = trackRef->DetectorId();
197 if (detectorId == AliTrackReference::kTPC) {
2d35f5d9 198 labelTPC = kTRUE;
f2247c19 199 break;
200 }
201 }
202 } // track references loop
203
204 // "Good" tracks
2d35f5d9 205 if (labelTPC) {
f2247c19 206 AliMCComparisonTrack* ref = new((*refs)[nt]) AliMCComparisonTrack();
207 ref->SetLabel(iTracks);
208 TParticle* particle = track->Particle();
209 Int_t pdg = particle->GetPdgCode();
210 ref->SetPDG(pdg);
211 ref->SetPz(trackRef->Pz());
2d35f5d9 212 Float_t pt = trackRef->Pt();
213 ref->SetPt(pt);
214 fGood->Fill(pt);
f2247c19 215 Float_t phig = trackRef->Phi();
216 ref->SetPhi(phig);
2d35f5d9 217 fGoodPhi->Fill(phig);
f2247c19 218 ref->SetLocalX(trackRef->LocalX());
219 ref->SetLocalY(trackRef->LocalY());
220 ref->SetZ(trackRef->Z());
221 nt++;
222 }
223 } // track loop
224
225
226 // ESD information
227 AliVEvent* event = InputEvent();
228 if (!event) {
229 Printf("ERROR: Could not retrieve event");
230 return;
231 }
232
233 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
234
235 Bool_t iFound;
236 Int_t nfound = 0;
237 Int_t nfake = 0;
238 Int_t nlost = 0;
239
2d35f5d9 240 Int_t mcGoods = refs->GetEntriesFast();
f2247c19 241
242 // Loop over all "good" MC tracks
2d35f5d9 243 for (Int_t k = 0; k < mcGoods; k++) {
f2247c19 244 AliMCComparisonTrack* ref = (AliMCComparisonTrack*)refs->UncheckedAt(k);
245 if (!ref) continue;
2d35f5d9 246 Int_t mcLabel = ref->GetLabel();
f2247c19 247 Float_t ptg = ref->GetPt();
2d35f5d9 248 Float_t phiG = ref->GetPhi();
f2247c19 249 iFound = kFALSE;
250 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
251 AliESDtrack* track = esd->GetTrack(iTrack);
252 if (!track) {
253 Printf("ERROR: Could not receive track %d", iTrack);
254 continue;
255 }
256
257 if (! track->IsOn(AliESDtrack::kTRDout)) continue;
258 const AliExternalTrackParam * outorig = track->GetOuterParam();
259 if (!outorig) continue;
260
2d35f5d9 261 Int_t lable = track->GetTRDLabel();
f2247c19 262
2d35f5d9 263 if (mcLabel == TMath::Abs(lable)) {
264 if (mcLabel == lable) {
f2247c19 265 nfound++;
2d35f5d9 266 fFound->Fill(ptg);
267 fFoundPhi->Fill(phiG);
f2247c19 268 }
269 else {
270 nfake++;
2d35f5d9 271 fFake->Fill(ptg);
f2247c19 272 }
273 iFound = kTRUE;
274
275 AliExternalTrackParam out(*outorig);
276
277 Double_t bz = esd->GetMagneticField();
278 Double_t xg = ref->GetLocalX();
279 out.PropagateTo(xg,bz);
280
281 Float_t phi = TMath::ASin(out.GetSnp()) + out.GetAlpha();
282 if (phi < 0) phi += 2*TMath::Pi();
283 if (phi >= 2*TMath::Pi()) phi -= 2*TMath::Pi();
284 Float_t phig = ref->GetPhi();
2d35f5d9 285 fP->Fill((phi - phig)*1000.);
f2247c19 286
287 Float_t lam = TMath::ATan(out.GetTgl());
288 Float_t lamg = TMath::ATan2(ref->GetPz(), ptg);
2d35f5d9 289 fL->Fill((lam - lamg)*1000.);
f2247c19 290
2d35f5d9 291 fC->Fill(track->GetTRDclusters(0));
f2247c19 292
2d35f5d9 293 Float_t pt1 = out.OneOverPt();
294 fPt->Fill((pt1 - 1/ptg)/(1/ptg)*100.);
f2247c19 295
296 Float_t y = out.GetY();
297 Float_t yg = ref->GetLocalY();
2d35f5d9 298 fHmpt->Fill(10*(y - yg));
f2247c19 299
300 Float_t z = out.GetZ();
301 Float_t zg = ref->GetZ();
2d35f5d9 302 fZ->Fill(10.*(z - zg));
f2247c19 303
2d35f5d9 304 Float_t mom = 1./(pt1*TMath::Cos(lam));
f2247c19 305 Float_t dedx = track->GetTRDsignal();
306 Printf (" dedx = %f ", dedx);
2d35f5d9 307 fEp->Fill(mom, dedx, 1.);
f2247c19 308
309 Int_t pdg = ref->GetPDG();
310 if (TMath::Abs(pdg)==211) //pions
2d35f5d9 311 if (mom>0.4 && mom<0.5) fE->Fill(dedx,1.);
f2247c19 312
313 break;
314 }
315 }
316 if (!iFound) {
317 nlost++;
318 }
319 }
320
321 Printf(" Results: " );
322 Printf(" Good %d Found %d Fake %d Lost %d ", nt, nfound, nfake, nlost);
323
324 refs->Clear();
325
326 // Post output data.
327 PostData(1, fListOfHistos);
328}
329
330
331void AliTRDComparisonTask::Terminate(Option_t *) {
332 // Draw result to the screen
333 // Called once at the end of the query
334 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
335 if (!fListOfHistos) {
336 Printf("ERROR: fListOfHistos not available");
337 return;
338 }
339
2d35f5d9 340 fGood = dynamic_cast<TH1F*>(fListOfHistos->At(0));
341 fFound = dynamic_cast<TH1F*>(fListOfHistos->At(1));
342 fFake = dynamic_cast<TH1F*>(fListOfHistos->At(2));
343 fP = dynamic_cast<TH1F*>(fListOfHistos->At(3));
344 fL = dynamic_cast<TH1F*>(fListOfHistos->At(4));
345 fPt = dynamic_cast<TH1F*>(fListOfHistos->At(5));
346 fHmpt = dynamic_cast<TH1F*>(fListOfHistos->At(6));
347 fE = dynamic_cast<TH1F*>(fListOfHistos->At(7));
348 fEp = dynamic_cast<TH2F*>(fListOfHistos->At(8));
349 fGoodPhi = dynamic_cast<TH1F*>(fListOfHistos->At(9));
350 fFoundPhi = dynamic_cast<TH1F*>(fListOfHistos->At(10));
351 fZ = dynamic_cast<TH1F*>(fListOfHistos->At(11));
352 fC = dynamic_cast<TH1F*>(fListOfHistos->At(12));
f2247c19 353
354 gStyle->SetOptStat(111110);
355 gStyle->SetOptFit(1);
356
357 TCanvas* c1 = new TCanvas("c1", "", 0, 0, 700, 850);
358
359 Int_t minc = 33;
360
361 TPad* p1 = new TPad("p1", "", 0., 0.3, 0.5, 0.6); p1->Draw();
362 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
2d35f5d9 363 fP->SetFillColor(4); fP->SetXTitle("(mrad)");
364 if (fP->GetEntries() < minc) fP->Draw(); else fP->Fit("gaus"); c1->cd();
f2247c19 365
366 TPad* p2 = new TPad("p2", "", 0.5, 0.3, 1., 0.6); p2->Draw();
367 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
2d35f5d9 368 fL->SetFillColor(4); fL->SetXTitle("(mrad)");
369 if (fL->GetEntries() < minc) fL->Draw(); else fL->Fit("gaus"); c1->cd();
f2247c19 370
371 TPad* p3 = new TPad("p3", "", 0., 0., 0.5, 0.3); p3->Draw();
372 p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
2d35f5d9 373 fPt->SetFillColor(2); fPt->SetXTitle("(%)");
374 if (fPt->GetEntries() < minc) fPt->Draw(); else fPt->Fit("gaus"); c1->cd();
f2247c19 375
376 TPad* p4 = new TPad("p4", "", 0.5, 0., 1., 0.3); p4->Draw();
377 p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
2d35f5d9 378 fHmpt->SetFillColor(6); fHmpt->SetXTitle("(mm)");
379 if (fHmpt->GetEntries() < minc) fHmpt->Draw(); else fHmpt->Fit("gaus");
380 fZ->Draw("same"); c1->cd();
f2247c19 381
382 TPad* p5 = new TPad("p5", "", 0., 0.6, 1., 1.); p5->Draw(); p5->cd();
383 p5->SetFillColor(41); p5->SetFrameFillColor(10);
2d35f5d9 384 fFound->Sumw2(); fGood->Sumw2(); fFake->Sumw2();
f2247c19 385 TH1F* hg = new TH1F("hg", "Efficiency for good tracks", 34, 0.2, 7.0);
386 TH1F* hf = new TH1F("hf", "Efficiency for fake tracks", 34, 0.2, 7.0);
2d35f5d9 387 hg->Divide(fFound,fGood,1.,1.,"B");
388 hf->Divide(fFake,fGood,1.,1.,"B");
f2247c19 389 hg->SetMaximum(1.4);
390 hg->SetYTitle("Tracking efficiency");
391 hg->SetXTitle("Pt (GeV/c)");
392 hg->Draw();
393
394 TLine* line1 = new TLine(0.2, 1.0, 7.0, 1.0); line1->SetLineStyle(4);
395 line1->Draw("same");
396 TLine* line2 = new TLine(0.2, 0.9, 7.0, 0.9); line2->SetLineStyle(4);
397 line2->Draw("same");
398
399 hf->SetFillColor(1);
400 hf->SetFillStyle(3013);
401 hf->SetLineColor(2);
402 hf->SetLineWidth(2);
403 hf->Draw("histsame");
404 TText* text = new TText(0.461176, 0.248448, "Fake tracks");
405 text->SetTextSize(0.05);
406 text->Draw();
407 text = new TText(0.453919, 1.11408, "Good tracks");
408 text->SetTextSize(0.05);
409 text->Draw();
410
411 TCanvas* c2 = new TCanvas("c2", "", 320, 32, 530, 590);
412
413 TPad* p6 = new TPad("p6", "", 0., 0., 1., .5); p6->Draw();
414 p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
2d35f5d9 415 fE->SetFillColor(2); fE->SetFillStyle(3005);
416 fE->SetXTitle("Arbitrary Units");
417 if (fE->GetEntries() < minc) fE->Draw(); else fE->Fit("gaus"); c2->cd();
f2247c19 418
419 TPad* p7 = new TPad("p7", "", 0., 0.5, 1., 1.); p7->Draw();
420 p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
2d35f5d9 421 fEp->SetMarkerStyle(8); fEp->SetMarkerSize(0.4);
422 fEp->SetFillColor(2); fEp->SetFillStyle(3005);
423 fEp->SetXTitle("p (GeV/c)"); fEp->SetYTitle("dE/dX (Arb. Units)");
424 fEp->Draw(); c1->cd();
f2247c19 425
426 TCanvas* c3 = new TCanvas("c3", "", 10, 10, 510, 510);
427 c3->SetFillColor(42); c3->SetFrameFillColor(10);
2d35f5d9 428 fFoundPhi->Sumw2(); fGoodPhi->Sumw2();
f2247c19 429 TH1F* hgphi = new TH1F("hgphi", "Efficiency for good tracks (Phi)",
430 90, 0., 2.*TMath::Pi());
2d35f5d9 431 hgphi->Divide(fFoundPhi, fGoodPhi, 1., 1., "B");
f2247c19 432 hgphi->SetYTitle("Tracking efficiency");
433 hgphi->SetXTitle("Phi (rad)");
434 hgphi->Draw();
435
436 TFile fc("AliTRDComparison.root","RECREATE");
437 c1->Write();
438 c2->Write();
439 c3->Write();
440 fc.Close();
441}