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