git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Move pad planes from AliTRDCommomParam to AliTRDgeometry
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackingAnalysis.cxx
CommitLineData
180da4e4 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18////////////////////////////////////////////////////////////////////
19// //
20// Fills a set of QA histograms to check the correctness of //
21// the TRD reconstruction //
22// //
23////////////////////////////////////////////////////////////////////
50f3bbf4 24
25#include "AliTRDtrackingAnalysis.h"
26
50f3bbf4 27#include "TFile.h"
28#include "TTree.h"
29#include "TH1D.h"
30#include "TH2D.h"
31#include "TObjArray.h"
32#include "TCanvas.h"
33#include "TGeoMatrix.h"
34#include "TStyle.h"
35#include "TGraphErrors.h"
36#include "TF1.h"
74cc0efb 37#include "TMath.h"
50f3bbf4 38
39#include "AliRunLoader.h"
40#include "AliTRDgeometry.h"
41#include "AliRun.h"
42#include "AliESD.h"
43#include "AliESDtrack.h"
44#include "AliTrackReference.h"
45
46#include "AliTRDcluster.h"
50f3bbf4 47#include "AliTRDpadPlane.h"
48#include "AliTRDcalibDB.h"
49#include "AliTracker.h"
50#include "AliTRDtracker.h"
51//#include "AliTRDtracklet.h"
52
53#include "TGeoManager.h"
54
55//////////////////////////////////////////////////////////////////////////////////////////
56
57AliTRDtrackingAnalysis::AliTRDtrackingAnalysis():
58 TObject(),
ae6817b3 59 fPath(0),
60 fRefTPC(0),
61 fRefTRD(0),
50f3bbf4 62 fLoader(0),
63 fEsdTree(0),
64 fESD(0),
ae6817b3 65 fTracker(0),
66 fDeltaPt(0),
67 fDeltaZ(0),
68 fDeltaX(0),
69 fDeltaYPos(0),
70 fDeltaYNeg(0),
71 fNPoints(0),
72 fNGood(0),
73 fRefSpace(0),
74 fGeo(0),
75 fClY2(0),
76 fClY3(0),
77 fTgPhi(0),
78 fGrResTgPhi(0),
79 fGrMeanTgPhi(0),
80 fTrklY(0),
81 fTrklZ(0),
82 fClZ(0),
83 fClZZ(0),
84 fClYY(0),
85 fClYX(0),
86 fNLabels(0),
180da4e4 87 fTestBits(0),
ae6817b3 88 fRefDx(0),
89 fClZXref(0),
90 fClZXcl(0),
91 fClPos(0)
50f3bbf4 92{
93
94 fDeltaX = new TH1D("deltaX", ";delta X (cm)", 100, -1, 1);
95 fDeltaZ = new TH1D("deltaZ", ";delta Z (cm)", 100, -2, 2);
96
97 fDeltaYPos = new TH1D("deltaYpos", ";delta Y (mm)", 100, -1, 1);
98 fDeltaYNeg = new TH1D("deltaYneg", ";delta Y (mm)", 100, -1, 1);
99
100 fNPoints = new TH1D("nPoints", ";np", 40, -0.5, 39.5);
101 fNGood = new TH1D("nGood", ";np", 40, -0.5, 39.5);
102
103 fDeltaPt = new TH1D("deltaPt", ";delta Pt/Pt (%)", 100, -10, 10);
104 fRefSpace = new TH2D("refSpace", ";y;x", 120, -60, 60, 200, -4, 1);
105
106 fTrklY = new TH1D("trklY", ";delta Y (mm)", 100, -1, 1);
107 fTrklZ = new TH1D("trklZ", ";delta Z (cm)", 100, -10, 10);
108
109
110 // cluster studies
111 fClY2 = new TH1D("clY2", ";delta Y (mm)", 100, -10, 10);
112 fClY3 = new TH1D("clY3", ";delta Y (mm)", 100, -10, 10);
113
114 for(int i=0; i<12; i++) // bewere hidden constants in the code
115 fClYTgPhi[i] = new TH1D(Form("clYtgPhi%d", i), ";delta Y (mm)", 100, -3, 3);
116
117 fTgPhi = new TH1D("tgPhi", ";Tg(#phi)", 100, -0.3, 0.3);
118 fGrResTgPhi = new TGraphErrors();
119 fGrMeanTgPhi = new TGraphErrors();
120
121 //fPullY2 = new TH1D("pullY2", ";pulls Y", 100, -5, 5);
122 //fPullY3 = new TH1D("pullY3", ";pulls Y", 100, -5, 5);
123
124
125 fClZ = new TH1D("clZ", ";delta Z (cm)", 200, -20, 20);
126 fClZZ = new TH2D("clZZ", ";z ref;z cl", 600, -300, 300, 600, -300, 300);
127
128 fClYY = new TH2D("clYY", ";dY;dY", 100, -3, 3, 100, -3, 3);
129 fClYX = new TH2D("clYX", ";Y;X", 250, -60, 60, 100, -4, 1);
130
131 fNLabels = new TH1D("clLabels", ";n labels", 10, -0.5, 9.5);
180da4e4 132 fTestBits = new TH1D("bits", ";bits", 10, -0.5, 9.5);
50f3bbf4 133
134 fRefDx = new TH1D("refDX", ";delta X", 100, 0, 20);
135 fClPos = new TH2D("clPos", ";z;y", 400, -400, 400, 120, -60, 60);
136
137 fClZXref = new TH2D("clZXref", ";z;x", 36, -54, 54, 300, 280, 380);
138 fClZXcl = new TH2D("clZXcl", ";z;x", 36, -54, 54, 300, 280, 380);
139
140 //fGeo = new AliTRDgeometry();
141}
142
143//////////////////////////////////////////////////////////////////////////////////////////
144
180da4e4 145AliTRDtrackingAnalysis::AliTRDtrackingAnalysis(const AliTRDtrackingAnalysis &t):
146 TObject(t),
147 fPath(0),
148 fRefTPC(0),
149 fRefTRD(0),
150 fLoader(0),
151 fEsdTree(0),
152 fESD(0),
153 fTracker(0),
154 fDeltaPt(0),
155 fDeltaZ(0),
156 fDeltaX(0),
157 fDeltaYPos(0),
158 fDeltaYNeg(0),
159 fNPoints(0),
160 fNGood(0),
161 fRefSpace(0),
162 fGeo(0),
163 fClY2(0),
164 fClY3(0),
165 fTgPhi(0),
166 fGrResTgPhi(0),
167 fGrMeanTgPhi(0),
168 fTrklY(0),
169 fTrklZ(0),
170 fClZ(0),
171 fClZZ(0),
172 fClYY(0),
173 fClYX(0),
174 fNLabels(0),
175 fTestBits(0),
176 fRefDx(0),
177 fClZXref(0),
178 fClZXcl(0),
179 fClPos(0)
180{
181
182 fDeltaX = new TH1D("deltaX", ";delta X (cm)", 100, -1, 1);
183 fDeltaZ = new TH1D("deltaZ", ";delta Z (cm)", 100, -2, 2);
184
185 fDeltaYPos = new TH1D("deltaYpos", ";delta Y (mm)", 100, -1, 1);
186 fDeltaYNeg = new TH1D("deltaYneg", ";delta Y (mm)", 100, -1, 1);
187
188 fNPoints = new TH1D("nPoints", ";np", 40, -0.5, 39.5);
189 fNGood = new TH1D("nGood", ";np", 40, -0.5, 39.5);
190
191 fDeltaPt = new TH1D("deltaPt", ";delta Pt/Pt (%)", 100, -10, 10);
192 fRefSpace = new TH2D("refSpace", ";y;x", 120, -60, 60, 200, -4, 1);
193
194 fTrklY = new TH1D("trklY", ";delta Y (mm)", 100, -1, 1);
195 fTrklZ = new TH1D("trklZ", ";delta Z (cm)", 100, -10, 10);
196
197
198 // cluster studies
199 fClY2 = new TH1D("clY2", ";delta Y (mm)", 100, -10, 10);
200 fClY3 = new TH1D("clY3", ";delta Y (mm)", 100, -10, 10);
201
202 for(int i=0; i<12; i++) // bewere hidden constants in the code
203 fClYTgPhi[i] = new TH1D(Form("clYtgPhi%d", i), ";delta Y (mm)", 100, -3, 3);
204
205 fTgPhi = new TH1D("tgPhi", ";Tg(#phi)", 100, -0.3, 0.3);
206 fGrResTgPhi = new TGraphErrors();
207 fGrMeanTgPhi = new TGraphErrors();
208
209 //fPullY2 = new TH1D("pullY2", ";pulls Y", 100, -5, 5);
210 //fPullY3 = new TH1D("pullY3", ";pulls Y", 100, -5, 5);
211
212
213 fClZ = new TH1D("clZ", ";delta Z (cm)", 200, -20, 20);
214 fClZZ = new TH2D("clZZ", ";z ref;z cl", 600, -300, 300, 600, -300, 300);
215
216 fClYY = new TH2D("clYY", ";dY;dY", 100, -3, 3, 100, -3, 3);
217 fClYX = new TH2D("clYX", ";Y;X", 250, -60, 60, 100, -4, 1);
218
219 fNLabels = new TH1D("clLabels", ";n labels", 10, -0.5, 9.5);
220 fTestBits = new TH1D("bits", ";bits", 10, -0.5, 9.5);
221
222 fRefDx = new TH1D("refDX", ";delta X", 100, 0, 20);
223 fClPos = new TH2D("clPos", ";z;y", 400, -400, 400, 120, -60, 60);
224
225 fClZXref = new TH2D("clZXref", ";z;x", 36, -54, 54, 300, 280, 380);
226 fClZXcl = new TH2D("clZXcl", ";z;x", 36, -54, 54, 300, 280, 380);
227
228 //fGeo = new AliTRDgeometry();
229}
230
231//////////////////////////////////////////////////////////////////////////////////////////
232
233void AliTRDtrackingAnalysis::DrawResolutionPt(int startEvent, int stopEvent)
234{
235 //
236 // Check the pt resolution
237 //
50f3bbf4 238
239 CheckFiles();
240
241 // loop over ESD events
242 int nevents = fEsdTree->GetEntries();
243
244 for(int iEvent=startEvent; iEvent<nevents && iEvent < stopEvent; iEvent++) {
245
246 Info("Draw", "Event = %d", iEvent);
247
248 fEsdTree->GetEvent(iEvent);
249 fLoader->GetEvent(iEvent);
250 LoadRefs();
251
252 int nTracks = fESD->GetNumberOfTracks();
253 for(int iTrack=0; iTrack<nTracks; iTrack++) {
254
255 //Info("Track", "Track = %d", iTrack);
256 AliESDtrack *esdTrack = fESD->GetTrack(iTrack);
257 if (!esdTrack->GetInnerParam()) continue;
258 const AliExternalTrackParam *param = esdTrack->GetOuterParam();
259 int status = esdTrack->GetStatus();
260
261 if (!(status & AliESDtrack::kTRDout)) continue;
262 if (!(status & AliESDtrack::kTRDrefit)) continue;
74cc0efb 263 if (TMath::Abs(esdTrack->GetOuterParam()->GetPt()) < 1.0) continue;
50f3bbf4 264
265 int ch=0;
266 while(param->GetX() > fGeo->GetTime0(ch)+2) ch++;
267 fRefSpace->Fill(2.*ch+0.5, param->GetX() - fGeo->GetTime0(ch));
268 //if (ch < 5) continue;
269
270 double lastX = 0;
271 int label = abs(esdTrack->GetTRDLabel());
272 int ntr = 0;
273 int ngood = 0;
274
275 for(int iPoint=GetReference(label); iPoint<fRefTRD->GetEntries(); iPoint++) {
276
277 AliTrackReference *aRef = (AliTrackReference*)(*fRefTRD)[iPoint];
278 if (aRef->GetTrack() != label) break;
279 ntr++;
280
281 lastX = (aRef->LocalX() < lastX)? lastX : aRef->LocalX();
282 double dx = aRef->LocalX() - param->GetX();
74cc0efb 283 if (TMath::Abs(dx) > 1.) continue;
50f3bbf4 284 ngood++;
285
286 double bz=fESD->GetMagneticField();
287 AliExternalTrackParam out(*param);
288 out.PropagateTo(aRef->LocalX(),bz);
289
290 double dp = aRef->Pt() + out.GetPt();
291 double dy = 10. * (aRef->LocalY() - out.GetY()); // in mm
292
293 fDeltaPt->Fill(100. * dp / aRef->Pt());
294 //fDeltaPt->Fill(out.GetPt() / aRef->Pt());
295 fDeltaX->Fill(dx);
296
297 if (esdTrack->GetSign() > 0) fDeltaYPos->Fill(dy);
298 else fDeltaYNeg->Fill(dy);
299
300 fDeltaZ->Fill(aRef->Z() - out.GetZ());
301 }
302
303 //if (ngood == 0) Info("X", "N = %d, X = %f, DX = %f", ntr, param->GetX(), param->GetX()-lastX);
304
305 fNPoints->Fill(ntr);
306 fNGood->Fill(ngood);
307 }
308 }
309
310 new TCanvas();
311 fDeltaPt->Draw();
312
313 new TCanvas();
314 fDeltaX->Draw();
315
316 /*
317 new TCanvas();
318 fNPoints->Draw();
319
320 new TCanvas();
321 fNGood->Draw();
322 */
323
324 new TCanvas();
325 fDeltaYPos->Draw();
326
327 new TCanvas();
328 fDeltaYNeg->Draw();
329
330 new TCanvas();
331 fDeltaZ->Draw();
332
333 new TCanvas();
334 fRefSpace->Draw();
335
336}
337
338//////////////////////////////////////////////////////////////////////////////////////////
339
340// void AliTRDtrackingAnalysis::DrawTrackletResolution(int startEvent, int stopEvent) {
341
342// LoadRecPointsFile();
343
344// TFile *file = new TFile(Form("%s/TRD.Tracklets.root", fPath), "read");
345// TTree *tree = (TTree*)file->Get("TRDtracklets");
346
347// AliTRDtracklet *tracklet = new AliTRDtracklet();
348// tree->SetBranchAddress("tracklets", &tracklet);
349
350// for(int ev=startEvent; ev<stopEvent; ev++) {
351
352// gAlice->GetEvent(ev);
353// LoadRefs();
354
355// int N = tree->GetEntries();
356// for(int i=0; i<N; i++) {
357
358// tree->GetEntry(i);
359
360// Double_t yref, zref, tgphi;
361// int stat = GetMCPosition(tracklet->GetLabel(), tracklet->GetX(), yref, zref, tgphi);
362// if (stat < 0) continue;
363
364// int plane = tracklet->GetPlane();
365// Double_t h01 = tracklet->GetTilt();
366
367// //printf("Tile = %f\tcorrection = %f um \n", h01, 1e4 * h01 * (tracklet->GetZ()-zref));
368// //double dz = zref - tracklet->GetZ() - cls->GetZ();
369
370// fTrklY->Fill(10 * (tracklet->GetY() - yref));
371// fTrklZ->Fill(tracklet->GetZ() - zref);
372
373// int ch=0;
374// while(tracklet->GetX() > fGeo->GetTime0(ch)+2) ch++;
375// fRefSpace->Fill(tracklet->GetY(), tracklet->GetX() - fGeo->GetTime0(ch));
376// }
377// }
378
379// new TCanvas();
380// fTrklZ->Draw();
381
382// new TCanvas();
383// fRefSpace->Draw();
384
385// gStyle->SetOptFit(1);
386// new TCanvas();
387// fTrklY->Draw();
388// fTrklY->Fit("gaus");
389
390// //new TCanvas();
391// //fClZ->Draw();
392// }
393
394//////////////////////////////////////////////////////////////////////////////////////////
395
180da4e4 396void AliTRDtrackingAnalysis::DrawRecPointResolution(int startEvent, int stopEvent)
397{
398 //
399 // Check the resolution of the reconstructed points
400 //
50f3bbf4 401
402 LoadRecPointsFile();
403 TObjArray *module = new TObjArray();
404
405 int nEvents = gAlice->GetEventsPerRun();
406
407 for(int ev=startEvent; ev<nEvents && ev < stopEvent; ev++) {
408
409 gAlice->GetEvent(ev);
410 LoadRefs();
411
412 TTree *tree = fLoader->GetTreeR("TRD", 0);
413 tree->SetBranchAddress("TRDcluster", &module);
414
415 Info("Res", "Refs Loaded");
416
180da4e4 417 Int_t nn = tree->GetEntries();
418 for(int i=0; i<nn; i++) {
50f3bbf4 419
420 tree->GetEntry(i);
421 int m = module->GetEntries();
422
423 for(int j=0; j<m; j++) {
424
425 AliTRDcluster *cls = (AliTRDcluster*)module->At(j);
426 if (cls->GetQ() < 10) continue;
427 fTracker->Transform(cls);
428 fClPos->Fill(cls->GetZ(), cls->GetY());
429
430 int plane = fGeo->GetPlane(cls->GetDetector());
431
432 int nl = 0;
433 for(int k=0; k<3; k++) if (cls->GetLabel(k) > -1) nl++;
434 fNLabels->Fill(nl);
435
436 Double_t yref, zref, tgphi;
437 int stat = GetMCPosition(cls->GetLabel(0), cls->GetX(), yref, zref, tgphi);
438 if (stat < 0) continue;
439
440 fClZXcl->Fill(cls->GetZ(), cls->GetX());
441 fClZXref->Fill(zref, cls->GetX());
442
f162af62 443 AliTRDpadPlane *padPlane = fGeo->GetPadPlane(plane,0);
50f3bbf4 444 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
445
446 //double dz = zref - padPlane->GetRow0();
447 double dz = zref - cls->GetZ();
448 double dy = dz * h01;
449 double yy = cls->GetY() - dy;
450
451 if (cls->From2pad()) fClY2->Fill(10 * (yy - yref));
452 if (cls->From3pad()) fClY3->Fill(10 * (yy - yref));
453
454 int idx = GetPhiBin(tgphi);
455 if (idx >= 0 && idx < 12) fClYTgPhi[idx]->Fill(10 * (yy - yref));
456
457 fClZZ->Fill(zref, cls->GetZ());
458 fClZ->Fill(dz);
459 fTgPhi->Fill(tgphi);
460 fClYX->Fill(cls->GetY(), cls->GetX() - fGeo->GetTime0(plane));
461 }
462 }
463 }
464
465 new TCanvas();
466 fClYX->Draw();
467
468 //new TCanvas();
469 //fNLabels->Draw();
470 new TCanvas();
471 fClZ->Draw();
472
473 new TCanvas();
474 fTgPhi->Draw();
475
476 gStyle->SetOptFit(1);
477
478 new TCanvas();
479 gPad->SetLogy();
480 fClY2->Draw();
481 fClY2->Fit("gaus", "", "", -2*fClY2->GetRMS(), 2*fClY2->GetRMS());
482
483 new TCanvas();
484 gPad->SetLogy();
485 fClY3->Draw();
486 fClY3->Fit("gaus", "", "", -2*fClY3->GetRMS(), 2*fClY3->GetRMS());
487
488 //new TCanvas();
489 //fRefSpace->Draw();
490
491 //new TCanvas();
492 //fClZXcl->Draw();
493
494 //new TCanvas();
495 //fClZXref->Draw();
496
497 /**/
498 TCanvas *c = new TCanvas();
499 c->Divide(4,3);
500
501 for(int i=0; i<12; i++) {
502
503 c->cd(i+1);
504 fClYTgPhi[i]->Draw();
505 if (fClYTgPhi[i]->GetSum() < 100) continue;
506
507 double mean = fClYTgPhi[i]->GetMean();
508 double rms = fClYTgPhi[i]->GetRMS();
509
510 fClYTgPhi[i]->Fit("gaus", "", "", mean-2*rms, mean+2*rms);
511 TF1 *f = fClYTgPhi[i]->GetFunction("gaus");
512
513 int n = fGrResTgPhi->GetN();
514 fGrResTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(2));
515 fGrResTgPhi->SetPointError(n, 0, f->GetParError(2));
516
517 fGrMeanTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(1));
518 fGrMeanTgPhi->SetPointError(n, 0, f->GetParError(1));
519 }
520
521 //gSystem->Sleep(1000);
522 gStyle->SetOptStat(0);
523
524 c = new TCanvas();
525 TH1D *dummy = new TH1D("dummy", "", 100, -0.3, 0.3);
526
527 dummy->SetTitle(";tg(#phi);resolution (mm)");
528 dummy->SetMinimum(0);
529 dummy->SetMaximum(1);
530
531 //c->cd();
532 (dummy->Clone("dummy1"))->Draw();
533
534 fGrResTgPhi->Draw("PL");
535 fGrResTgPhi->GetHistogram()->SetTitle(";tg(#phi);resolution (mm)");
536 fGrResTgPhi->SetMarkerStyle(20);
537
538 c = new TCanvas();
539 dummy->SetTitle(";tg(#phi);mean value (mm)");
540 dummy->SetMinimum(-0.3);
541 dummy->SetMaximum(0.3);
542
543 //c->cd();
544 (dummy->Clone("dummy2"))->Draw();
545
546 fGrMeanTgPhi->Draw("PL");
547 fGrMeanTgPhi->GetHistogram()->SetTitle(";tg(#phi);mean value (mm)");
548 fGrMeanTgPhi->SetMarkerStyle(20);
549 /**/
550
551
552 //new TCanvas();
553 //fClZZ->Draw("colz");
554
555 //new TCanvas();
556 //fClPos->Draw("colz");
557
558 //new TCanvas();
180da4e4 559 //fTestBits->Draw();
50f3bbf4 560
561 //new TCanvas();
562 //fRefDx->Draw();
563
564}
565
566//////////////////////////////////////////////////////////////////////////////////////////
567
180da4e4 568void AliTRDtrackingAnalysis::LoadRecPointsFile()
569{
570 //
571 // Load the clusters from the input file
572 //
50f3bbf4 573
574 char filename[256];
575 sprintf(filename, "%s/galice.root", fPath);
576
577 fLoader = AliRunLoader::Open(filename);
578 if (!fLoader) {
579 Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
580 return;
581 }
582
583 fLoader->LoadgAlice();
584 gAlice = fLoader->GetAliRun();
585
586 if (!gAlice) {
587 Error("CheckFiles", "no galice object found");
588 return;
589 }
590
591 fLoader->LoadKinematics();
592 fLoader->LoadHeader();
593 fLoader->LoadTrackRefs();
594
595 //TGeoManager::Import("/data/alice_u/radomski/condor/run_0/geometry.root");
596 TGeoManager::Import(Form("%s/geometry.root", fPath));
597
598
599 fLoader->CdGAFile();
600 fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
601 fTracker = new AliTRDtracker(gFile);
602
603 fLoader->LoadRecPoints("TRD");
604
605 AliTracker::SetFieldMap(gAlice->Field(), 1);
606}
607
608//////////////////////////////////////////////////////////////////////////////////////////
609
180da4e4 610void AliTRDtrackingAnalysis::CheckFiles()
611{
612 //
613 // Check the presence of the input files
614 //
50f3bbf4 615
616 // MC info
617
618 char filename[256];
619 sprintf(filename, "%s/galice.root", fPath);
620
621 fLoader = AliRunLoader::Open(filename);
622 if (!fLoader) {
623 Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
624 return;
625 }
626
627 fLoader->LoadgAlice();
628 gAlice = fLoader->GetAliRun();
629
630 if (!gAlice) {
631 Error("CheckFiles", "no galice object found");
632 return;
633 }
634
635 fLoader->LoadKinematics();
636 fLoader->LoadHeader();
637 fLoader->LoadTrackRefs();
638
639 fLoader->CdGAFile();
640 fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
641 //fGeo->ReadGeoMatrices();
642
643 // ESD
644
645 sprintf(filename,"%s/AliESDs.root", fPath);
646 TFile *esdFile = new TFile(filename, "READ");
647
648 if (esdFile->IsZombie()) {
649 Error("CheckFiles", "file not present: AliESDs.root");
650 return;
651 }
652
653 fEsdTree = (TTree*)esdFile->Get("esdTree");
654 fESD = new AliESD();
f162af62 655 fESD->ReadFromTree(fEsdTree);
656 //fEsdTree->SetBranchAddress("ESD", &fESD);
50f3bbf4 657}
658
659//////////////////////////////////////////////////////////////////////////////////////////
660
180da4e4 661void AliTRDtrackingAnalysis::LoadRefs()
662{
663 //
664 // Load the track references
665 //
50f3bbf4 666
667 if (fRefTPC) delete fRefTPC;
668 if (fRefTRD) delete fRefTRD;
669
670 fRefTPC = new TObjArray();
671 fRefTRD = new TObjArray();
672
673 //fLoader->GetEvent(event);
674 //AliStack* stack = gAlice->Stack();
675 TTree *refTree = fLoader->TreeTR();
676
180da4e4 677 const int kNBranch = 2;
50f3bbf4 678 const char *brName[] = {"TPC", "TRD"};
679 TClonesArray *clRefs = new TClonesArray("AliTrackReference");
680
180da4e4 681 for(int b=0; b<kNBranch; b++) {
50f3bbf4 682
683 TBranch *branch = refTree->GetBranch(brName[b]);
684 refTree->SetBranchAddress(brName[b],&clRefs);
685
686 int nEntries = branch->GetEntries();
687 for(int iTrack = 0; iTrack < nEntries; iTrack++) {
688
689 refTree->GetEvent(iTrack);
690 int nPoints = clRefs->GetEntries();
691 for(int iPoint=0; iPoint<nPoints; iPoint++) {
692 AliTrackReference *ref = (AliTrackReference*)clRefs->At(iPoint);
693 if (b == 0) fRefTPC->Add(new AliTrackReference(*ref));
694 if (b == 1) fRefTRD->Add(new AliTrackReference(*ref));
695 }
696 }
697 }
698
699 fRefTPC->Sort();
700 fRefTRD->Sort();
701
702 for(int i=0; i<fRefTRD->GetEntries(); i++) {
703 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
704 fLabels[i] = ref->GetTrack();
705
706 int p=0;
707 while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
708 fRefSpace->Fill(ref->LocalY(), ref->LocalX()-fGeo->GetTime0(p));
709
180da4e4 710 //for(int bit=0; bit<9; bit++) if (ref->TestBit(bit)) fTestBits->Fill(bit);
50f3bbf4 711 }
712
713 delete clRefs;
714 Info("LoadRefs", "TPC = %d\t TRD = %d", fRefTPC->GetEntries(), fRefTRD->GetEntries());
715}
716
717//////////////////////////////////////////////////////////////////////////////////////////
718
180da4e4 719Int_t AliTRDtrackingAnalysis::GetReference(Int_t label)
720{
721 //
722 // Sort the track references
723 //
50f3bbf4 724
725 int start = TMath::BinarySearch(fRefTRD->GetEntries(), fLabels, label);
726
727 while (start >= 0) {
728 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[start];
729 if (ref->GetTrack() != label) return start+1;
730 start--;
731 }
732
733 return 0;
734}
735
736//////////////////////////////////////////////////////////////////////////////////////////
737
180da4e4 738int AliTRDtrackingAnalysis::GetMCPosition(Int_t label, Double_t x, Double_t &Y, Double_t &Z, Double_t &tgphi)
739{
740 //
741 // Determine the MC positions from the track references
742 //
50f3bbf4 743
744 double lowX = 100.;
745 double highX = 100.;
746 int idLow = -1;
747 int idHigh = -1;
748
749 int nref= 0;
750 int idx = GetReference(label);
751 for(int i=idx; i<fRefTRD->GetEntries(); i++) {
752
753 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
754 if (ref->GetTrack() != label) break;
755 nref++;
756
757 //int p=0;
758 //while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
759 //if (p != layer) continue;
760
761 double dX = ref->LocalX()-x;
762 if ( dX > 0 ) {
763 if (dX < highX) {
764 idHigh = i;
765 highX = dX;
766 }
767 } else {
74cc0efb 768 dX = TMath::Abs(dX);
50f3bbf4 769 if (dX < lowX) {
770 idLow = i;
771 lowX = dX;
772 }
773 }
774 }
775
776 if (idLow == -1 || idHigh == -1) return -1;
777
778 AliTrackReference *refI = (AliTrackReference*)(*fRefTRD)[idLow];
779 AliTrackReference *refO = (AliTrackReference*)(*fRefTRD)[idHigh];
780
781
782 double dx = refO->LocalX() - refI->LocalX();
783 double dy = refO->LocalY() - refI->LocalY();
784 double dz = refO->Z() - refI->Z();
785 double ddx = (x - refI->LocalX())/dx;
786
787 fRefDx->Fill(dx);
788
789 Y = refI->LocalY() + ddx * dy;
790 Z = refI->Z() + ddx * dz;
791
792 tgphi = dy/dx;
793
794 return 0;
795}
796
797//////////////////////////////////////////////////////////////////////////////////////////
180da4e4 798Int_t AliTRDtrackingAnalysis::GetPhiBin(Double_t phi) const
799{
800 //
801 // Return the phi bin
802 //
50f3bbf4 803 return (int)((phi+0.3)/0.05);
804}
805
806//////////////////////////////////////////////////////////////////////////////////////////
180da4e4 807Double_t AliTRDtrackingAnalysis::GetPhi(Int_t bin) const
808{
809 //
810 // Return phi for a given bin
811 //
50f3bbf4 812 return bin * 0.05 - 0.3 + 0.025;
813}
814//////////////////////////////////////////////////////////////////////////////////////////