]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtrackingAnalysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Reconstruction QA by Sylwester
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackingAnalysis.cxx
CommitLineData
50f3bbf4 1
2#include "AliTRDtrackingAnalysis.h"
3
4
5#include "TFile.h"
6#include "TTree.h"
7#include "TH1D.h"
8#include "TH2D.h"
9#include "TObjArray.h"
10#include "TCanvas.h"
11#include "TGeoMatrix.h"
12#include "TStyle.h"
13#include "TGraphErrors.h"
14#include "TF1.h"
15
16#include "AliRunLoader.h"
17#include "AliTRDgeometry.h"
18#include "AliRun.h"
19#include "AliESD.h"
20#include "AliESDtrack.h"
21#include "AliTrackReference.h"
22
23#include "AliTRDcluster.h"
24#include "AliTRDCommonParam.h"
25#include "AliTRDpadPlane.h"
26#include "AliTRDcalibDB.h"
27#include "AliTracker.h"
28#include "AliTRDtracker.h"
29//#include "AliTRDtracklet.h"
30
31#include "TGeoManager.h"
32
33//////////////////////////////////////////////////////////////////////////////////////////
34
35AliTRDtrackingAnalysis::AliTRDtrackingAnalysis():
36 TObject(),
37 fLoader(0),
38 fEsdTree(0),
39 fESD(0),
40 fRefTPC(0),
41 fRefTRD(0)
42{
43
44 fDeltaX = new TH1D("deltaX", ";delta X (cm)", 100, -1, 1);
45 fDeltaZ = new TH1D("deltaZ", ";delta Z (cm)", 100, -2, 2);
46
47 fDeltaYPos = new TH1D("deltaYpos", ";delta Y (mm)", 100, -1, 1);
48 fDeltaYNeg = new TH1D("deltaYneg", ";delta Y (mm)", 100, -1, 1);
49
50 fNPoints = new TH1D("nPoints", ";np", 40, -0.5, 39.5);
51 fNGood = new TH1D("nGood", ";np", 40, -0.5, 39.5);
52
53 fDeltaPt = new TH1D("deltaPt", ";delta Pt/Pt (%)", 100, -10, 10);
54 fRefSpace = new TH2D("refSpace", ";y;x", 120, -60, 60, 200, -4, 1);
55
56 fTrklY = new TH1D("trklY", ";delta Y (mm)", 100, -1, 1);
57 fTrklZ = new TH1D("trklZ", ";delta Z (cm)", 100, -10, 10);
58
59
60 // cluster studies
61 fClY2 = new TH1D("clY2", ";delta Y (mm)", 100, -10, 10);
62 fClY3 = new TH1D("clY3", ";delta Y (mm)", 100, -10, 10);
63
64 for(int i=0; i<12; i++) // bewere hidden constants in the code
65 fClYTgPhi[i] = new TH1D(Form("clYtgPhi%d", i), ";delta Y (mm)", 100, -3, 3);
66
67 fTgPhi = new TH1D("tgPhi", ";Tg(#phi)", 100, -0.3, 0.3);
68 fGrResTgPhi = new TGraphErrors();
69 fGrMeanTgPhi = new TGraphErrors();
70
71 //fPullY2 = new TH1D("pullY2", ";pulls Y", 100, -5, 5);
72 //fPullY3 = new TH1D("pullY3", ";pulls Y", 100, -5, 5);
73
74
75 fClZ = new TH1D("clZ", ";delta Z (cm)", 200, -20, 20);
76 fClZZ = new TH2D("clZZ", ";z ref;z cl", 600, -300, 300, 600, -300, 300);
77
78 fClYY = new TH2D("clYY", ";dY;dY", 100, -3, 3, 100, -3, 3);
79 fClYX = new TH2D("clYX", ";Y;X", 250, -60, 60, 100, -4, 1);
80
81 fNLabels = new TH1D("clLabels", ";n labels", 10, -0.5, 9.5);
82 fBits = new TH1D("bits", ";bits", 10, -0.5, 9.5);
83
84 fRefDx = new TH1D("refDX", ";delta X", 100, 0, 20);
85 fClPos = new TH2D("clPos", ";z;y", 400, -400, 400, 120, -60, 60);
86
87 fClZXref = new TH2D("clZXref", ";z;x", 36, -54, 54, 300, 280, 380);
88 fClZXcl = new TH2D("clZXcl", ";z;x", 36, -54, 54, 300, 280, 380);
89
90 //fGeo = new AliTRDgeometry();
91}
92
93//////////////////////////////////////////////////////////////////////////////////////////
94
95void AliTRDtrackingAnalysis::DrawResolutionPt(int startEvent, int stopEvent) {
96
97 CheckFiles();
98
99 // loop over ESD events
100 int nevents = fEsdTree->GetEntries();
101
102 for(int iEvent=startEvent; iEvent<nevents && iEvent < stopEvent; iEvent++) {
103
104 Info("Draw", "Event = %d", iEvent);
105
106 fEsdTree->GetEvent(iEvent);
107 fLoader->GetEvent(iEvent);
108 LoadRefs();
109
110 int nTracks = fESD->GetNumberOfTracks();
111 for(int iTrack=0; iTrack<nTracks; iTrack++) {
112
113 //Info("Track", "Track = %d", iTrack);
114 AliESDtrack *esdTrack = fESD->GetTrack(iTrack);
115 if (!esdTrack->GetInnerParam()) continue;
116 const AliExternalTrackParam *param = esdTrack->GetOuterParam();
117 int status = esdTrack->GetStatus();
118
119 if (!(status & AliESDtrack::kTRDout)) continue;
120 if (!(status & AliESDtrack::kTRDrefit)) continue;
121 if (fabs(esdTrack->GetOuterParam()->GetPt()) < 1.0) continue;
122
123 int ch=0;
124 while(param->GetX() > fGeo->GetTime0(ch)+2) ch++;
125 fRefSpace->Fill(2.*ch+0.5, param->GetX() - fGeo->GetTime0(ch));
126 //if (ch < 5) continue;
127
128 double lastX = 0;
129 int label = abs(esdTrack->GetTRDLabel());
130 int ntr = 0;
131 int ngood = 0;
132
133 for(int iPoint=GetReference(label); iPoint<fRefTRD->GetEntries(); iPoint++) {
134
135 AliTrackReference *aRef = (AliTrackReference*)(*fRefTRD)[iPoint];
136 if (aRef->GetTrack() != label) break;
137 ntr++;
138
139 lastX = (aRef->LocalX() < lastX)? lastX : aRef->LocalX();
140 double dx = aRef->LocalX() - param->GetX();
141 if (fabs(dx) > 1.) continue;
142 ngood++;
143
144 double bz=fESD->GetMagneticField();
145 AliExternalTrackParam out(*param);
146 out.PropagateTo(aRef->LocalX(),bz);
147
148 double dp = aRef->Pt() + out.GetPt();
149 double dy = 10. * (aRef->LocalY() - out.GetY()); // in mm
150
151 fDeltaPt->Fill(100. * dp / aRef->Pt());
152 //fDeltaPt->Fill(out.GetPt() / aRef->Pt());
153 fDeltaX->Fill(dx);
154
155 if (esdTrack->GetSign() > 0) fDeltaYPos->Fill(dy);
156 else fDeltaYNeg->Fill(dy);
157
158 fDeltaZ->Fill(aRef->Z() - out.GetZ());
159 }
160
161 //if (ngood == 0) Info("X", "N = %d, X = %f, DX = %f", ntr, param->GetX(), param->GetX()-lastX);
162
163 fNPoints->Fill(ntr);
164 fNGood->Fill(ngood);
165 }
166 }
167
168 new TCanvas();
169 fDeltaPt->Draw();
170
171 new TCanvas();
172 fDeltaX->Draw();
173
174 /*
175 new TCanvas();
176 fNPoints->Draw();
177
178 new TCanvas();
179 fNGood->Draw();
180 */
181
182 new TCanvas();
183 fDeltaYPos->Draw();
184
185 new TCanvas();
186 fDeltaYNeg->Draw();
187
188 new TCanvas();
189 fDeltaZ->Draw();
190
191 new TCanvas();
192 fRefSpace->Draw();
193
194}
195
196//////////////////////////////////////////////////////////////////////////////////////////
197
198// void AliTRDtrackingAnalysis::DrawTrackletResolution(int startEvent, int stopEvent) {
199
200// LoadRecPointsFile();
201
202// TFile *file = new TFile(Form("%s/TRD.Tracklets.root", fPath), "read");
203// TTree *tree = (TTree*)file->Get("TRDtracklets");
204
205// AliTRDtracklet *tracklet = new AliTRDtracklet();
206// tree->SetBranchAddress("tracklets", &tracklet);
207
208// for(int ev=startEvent; ev<stopEvent; ev++) {
209
210// gAlice->GetEvent(ev);
211// LoadRefs();
212
213// int N = tree->GetEntries();
214// for(int i=0; i<N; i++) {
215
216// tree->GetEntry(i);
217
218// Double_t yref, zref, tgphi;
219// int stat = GetMCPosition(tracklet->GetLabel(), tracklet->GetX(), yref, zref, tgphi);
220// if (stat < 0) continue;
221
222// int plane = tracklet->GetPlane();
223// Double_t h01 = tracklet->GetTilt();
224
225// //printf("Tile = %f\tcorrection = %f um \n", h01, 1e4 * h01 * (tracklet->GetZ()-zref));
226// //double dz = zref - tracklet->GetZ() - cls->GetZ();
227
228// fTrklY->Fill(10 * (tracklet->GetY() - yref));
229// fTrklZ->Fill(tracklet->GetZ() - zref);
230
231// int ch=0;
232// while(tracklet->GetX() > fGeo->GetTime0(ch)+2) ch++;
233// fRefSpace->Fill(tracklet->GetY(), tracklet->GetX() - fGeo->GetTime0(ch));
234// }
235// }
236
237// new TCanvas();
238// fTrklZ->Draw();
239
240// new TCanvas();
241// fRefSpace->Draw();
242
243// gStyle->SetOptFit(1);
244// new TCanvas();
245// fTrklY->Draw();
246// fTrklY->Fit("gaus");
247
248// //new TCanvas();
249// //fClZ->Draw();
250// }
251
252//////////////////////////////////////////////////////////////////////////////////////////
253
254void AliTRDtrackingAnalysis::DrawRecPointResolution(int startEvent, int stopEvent) {
255
256 LoadRecPointsFile();
257 TObjArray *module = new TObjArray();
258
259 int nEvents = gAlice->GetEventsPerRun();
260
261 for(int ev=startEvent; ev<nEvents && ev < stopEvent; ev++) {
262
263 gAlice->GetEvent(ev);
264 LoadRefs();
265
266 TTree *tree = fLoader->GetTreeR("TRD", 0);
267 tree->SetBranchAddress("TRDcluster", &module);
268
269 Info("Res", "Refs Loaded");
270
271 int N = tree->GetEntries();
272 for(int i=0; i<N; i++) {
273
274 tree->GetEntry(i);
275 int m = module->GetEntries();
276
277 for(int j=0; j<m; j++) {
278
279 AliTRDcluster *cls = (AliTRDcluster*)module->At(j);
280 if (cls->GetQ() < 10) continue;
281 fTracker->Transform(cls);
282 fClPos->Fill(cls->GetZ(), cls->GetY());
283
284 int plane = fGeo->GetPlane(cls->GetDetector());
285
286 int nl = 0;
287 for(int k=0; k<3; k++) if (cls->GetLabel(k) > -1) nl++;
288 fNLabels->Fill(nl);
289
290 Double_t yref, zref, tgphi;
291 int stat = GetMCPosition(cls->GetLabel(0), cls->GetX(), yref, zref, tgphi);
292 if (stat < 0) continue;
293
294 fClZXcl->Fill(cls->GetZ(), cls->GetX());
295 fClZXref->Fill(zref, cls->GetX());
296
297 AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
298 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
299
300 //double dz = zref - padPlane->GetRow0();
301 double dz = zref - cls->GetZ();
302 double dy = dz * h01;
303 double yy = cls->GetY() - dy;
304
305 if (cls->From2pad()) fClY2->Fill(10 * (yy - yref));
306 if (cls->From3pad()) fClY3->Fill(10 * (yy - yref));
307
308 int idx = GetPhiBin(tgphi);
309 if (idx >= 0 && idx < 12) fClYTgPhi[idx]->Fill(10 * (yy - yref));
310
311 fClZZ->Fill(zref, cls->GetZ());
312 fClZ->Fill(dz);
313 fTgPhi->Fill(tgphi);
314 fClYX->Fill(cls->GetY(), cls->GetX() - fGeo->GetTime0(plane));
315 }
316 }
317 }
318
319 new TCanvas();
320 fClYX->Draw();
321
322 //new TCanvas();
323 //fNLabels->Draw();
324 new TCanvas();
325 fClZ->Draw();
326
327 new TCanvas();
328 fTgPhi->Draw();
329
330 gStyle->SetOptFit(1);
331
332 new TCanvas();
333 gPad->SetLogy();
334 fClY2->Draw();
335 fClY2->Fit("gaus", "", "", -2*fClY2->GetRMS(), 2*fClY2->GetRMS());
336
337 new TCanvas();
338 gPad->SetLogy();
339 fClY3->Draw();
340 fClY3->Fit("gaus", "", "", -2*fClY3->GetRMS(), 2*fClY3->GetRMS());
341
342 //new TCanvas();
343 //fRefSpace->Draw();
344
345 //new TCanvas();
346 //fClZXcl->Draw();
347
348 //new TCanvas();
349 //fClZXref->Draw();
350
351 /**/
352 TCanvas *c = new TCanvas();
353 c->Divide(4,3);
354
355 for(int i=0; i<12; i++) {
356
357 c->cd(i+1);
358 fClYTgPhi[i]->Draw();
359 if (fClYTgPhi[i]->GetSum() < 100) continue;
360
361 double mean = fClYTgPhi[i]->GetMean();
362 double rms = fClYTgPhi[i]->GetRMS();
363
364 fClYTgPhi[i]->Fit("gaus", "", "", mean-2*rms, mean+2*rms);
365 TF1 *f = fClYTgPhi[i]->GetFunction("gaus");
366
367 int n = fGrResTgPhi->GetN();
368 fGrResTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(2));
369 fGrResTgPhi->SetPointError(n, 0, f->GetParError(2));
370
371 fGrMeanTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(1));
372 fGrMeanTgPhi->SetPointError(n, 0, f->GetParError(1));
373 }
374
375 //gSystem->Sleep(1000);
376 gStyle->SetOptStat(0);
377
378 c = new TCanvas();
379 TH1D *dummy = new TH1D("dummy", "", 100, -0.3, 0.3);
380
381 dummy->SetTitle(";tg(#phi);resolution (mm)");
382 dummy->SetMinimum(0);
383 dummy->SetMaximum(1);
384
385 //c->cd();
386 (dummy->Clone("dummy1"))->Draw();
387
388 fGrResTgPhi->Draw("PL");
389 fGrResTgPhi->GetHistogram()->SetTitle(";tg(#phi);resolution (mm)");
390 fGrResTgPhi->SetMarkerStyle(20);
391
392 c = new TCanvas();
393 dummy->SetTitle(";tg(#phi);mean value (mm)");
394 dummy->SetMinimum(-0.3);
395 dummy->SetMaximum(0.3);
396
397 //c->cd();
398 (dummy->Clone("dummy2"))->Draw();
399
400 fGrMeanTgPhi->Draw("PL");
401 fGrMeanTgPhi->GetHistogram()->SetTitle(";tg(#phi);mean value (mm)");
402 fGrMeanTgPhi->SetMarkerStyle(20);
403 /**/
404
405
406 //new TCanvas();
407 //fClZZ->Draw("colz");
408
409 //new TCanvas();
410 //fClPos->Draw("colz");
411
412 //new TCanvas();
413 //fBits->Draw();
414
415 //new TCanvas();
416 //fRefDx->Draw();
417
418}
419
420//////////////////////////////////////////////////////////////////////////////////////////
421
422void AliTRDtrackingAnalysis::LoadRecPointsFile() {
423
424 char filename[256];
425 sprintf(filename, "%s/galice.root", fPath);
426
427 fLoader = AliRunLoader::Open(filename);
428 if (!fLoader) {
429 Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
430 return;
431 }
432
433 fLoader->LoadgAlice();
434 gAlice = fLoader->GetAliRun();
435
436 if (!gAlice) {
437 Error("CheckFiles", "no galice object found");
438 return;
439 }
440
441 fLoader->LoadKinematics();
442 fLoader->LoadHeader();
443 fLoader->LoadTrackRefs();
444
445 //TGeoManager::Import("/data/alice_u/radomski/condor/run_0/geometry.root");
446 TGeoManager::Import(Form("%s/geometry.root", fPath));
447
448
449 fLoader->CdGAFile();
450 fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
451 fTracker = new AliTRDtracker(gFile);
452
453 fLoader->LoadRecPoints("TRD");
454
455 AliTracker::SetFieldMap(gAlice->Field(), 1);
456}
457
458//////////////////////////////////////////////////////////////////////////////////////////
459
460void AliTRDtrackingAnalysis::CheckFiles() {
461
462 // MC info
463
464 char filename[256];
465 sprintf(filename, "%s/galice.root", fPath);
466
467 fLoader = AliRunLoader::Open(filename);
468 if (!fLoader) {
469 Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
470 return;
471 }
472
473 fLoader->LoadgAlice();
474 gAlice = fLoader->GetAliRun();
475
476 if (!gAlice) {
477 Error("CheckFiles", "no galice object found");
478 return;
479 }
480
481 fLoader->LoadKinematics();
482 fLoader->LoadHeader();
483 fLoader->LoadTrackRefs();
484
485 fLoader->CdGAFile();
486 fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
487 //fGeo->ReadGeoMatrices();
488
489 // ESD
490
491 sprintf(filename,"%s/AliESDs.root", fPath);
492 TFile *esdFile = new TFile(filename, "READ");
493
494 if (esdFile->IsZombie()) {
495 Error("CheckFiles", "file not present: AliESDs.root");
496 return;
497 }
498
499 fEsdTree = (TTree*)esdFile->Get("esdTree");
500 fESD = new AliESD();
501
502 fEsdTree->SetBranchAddress("ESD", &fESD);
503}
504
505//////////////////////////////////////////////////////////////////////////////////////////
506
507void AliTRDtrackingAnalysis::LoadRefs() {
508
509 if (fRefTPC) delete fRefTPC;
510 if (fRefTRD) delete fRefTRD;
511
512 fRefTPC = new TObjArray();
513 fRefTRD = new TObjArray();
514
515 //fLoader->GetEvent(event);
516 //AliStack* stack = gAlice->Stack();
517 TTree *refTree = fLoader->TreeTR();
518
519 const int nBranch = 2;
520 const char *brName[] = {"TPC", "TRD"};
521 TClonesArray *clRefs = new TClonesArray("AliTrackReference");
522
523 for(int b=0; b<nBranch; b++) {
524
525 TBranch *branch = refTree->GetBranch(brName[b]);
526 refTree->SetBranchAddress(brName[b],&clRefs);
527
528 int nEntries = branch->GetEntries();
529 for(int iTrack = 0; iTrack < nEntries; iTrack++) {
530
531 refTree->GetEvent(iTrack);
532 int nPoints = clRefs->GetEntries();
533 for(int iPoint=0; iPoint<nPoints; iPoint++) {
534 AliTrackReference *ref = (AliTrackReference*)clRefs->At(iPoint);
535 if (b == 0) fRefTPC->Add(new AliTrackReference(*ref));
536 if (b == 1) fRefTRD->Add(new AliTrackReference(*ref));
537 }
538 }
539 }
540
541 fRefTPC->Sort();
542 fRefTRD->Sort();
543
544 for(int i=0; i<fRefTRD->GetEntries(); i++) {
545 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
546 fLabels[i] = ref->GetTrack();
547
548 int p=0;
549 while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
550 fRefSpace->Fill(ref->LocalY(), ref->LocalX()-fGeo->GetTime0(p));
551
552 //for(int bit=0; bit<9; bit++) if (ref->TestBit(bit)) fBits->Fill(bit);
553 }
554
555 delete clRefs;
556 Info("LoadRefs", "TPC = %d\t TRD = %d", fRefTPC->GetEntries(), fRefTRD->GetEntries());
557}
558
559//////////////////////////////////////////////////////////////////////////////////////////
560
561Int_t AliTRDtrackingAnalysis::GetReference(Int_t label) {
562
563 int start = TMath::BinarySearch(fRefTRD->GetEntries(), fLabels, label);
564
565 while (start >= 0) {
566 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[start];
567 if (ref->GetTrack() != label) return start+1;
568 start--;
569 }
570
571 return 0;
572}
573
574//////////////////////////////////////////////////////////////////////////////////////////
575
576int AliTRDtrackingAnalysis::GetMCPosition(Int_t label, Double_t x, Double_t &Y, Double_t &Z, Double_t &tgphi) {
577
578 double lowX = 100.;
579 double highX = 100.;
580 int idLow = -1;
581 int idHigh = -1;
582
583 int nref= 0;
584 int idx = GetReference(label);
585 for(int i=idx; i<fRefTRD->GetEntries(); i++) {
586
587 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
588 if (ref->GetTrack() != label) break;
589 nref++;
590
591 //int p=0;
592 //while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
593 //if (p != layer) continue;
594
595 double dX = ref->LocalX()-x;
596 if ( dX > 0 ) {
597 if (dX < highX) {
598 idHigh = i;
599 highX = dX;
600 }
601 } else {
602 dX = fabs(dX);
603 if (dX < lowX) {
604 idLow = i;
605 lowX = dX;
606 }
607 }
608 }
609
610 if (idLow == -1 || idHigh == -1) return -1;
611
612 AliTrackReference *refI = (AliTrackReference*)(*fRefTRD)[idLow];
613 AliTrackReference *refO = (AliTrackReference*)(*fRefTRD)[idHigh];
614
615
616 double dx = refO->LocalX() - refI->LocalX();
617 double dy = refO->LocalY() - refI->LocalY();
618 double dz = refO->Z() - refI->Z();
619 double ddx = (x - refI->LocalX())/dx;
620
621 fRefDx->Fill(dx);
622
623 Y = refI->LocalY() + ddx * dy;
624 Z = refI->Z() + ddx * dz;
625
626 tgphi = dy/dx;
627
628 return 0;
629}
630
631//////////////////////////////////////////////////////////////////////////////////////////
632Int_t AliTRDtrackingAnalysis::GetPhiBin(Double_t phi) {
633 return (int)((phi+0.3)/0.05);
634}
635
636//////////////////////////////////////////////////////////////////////////////////////////
637
638Double_t AliTRDtrackingAnalysis::GetPhi(Int_t bin) {
639 return bin * 0.05 - 0.3 + 0.025;
640}
641//////////////////////////////////////////////////////////////////////////////////////////