]> 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
Obsolete.
[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"
af885e0f 42#include "AliESDEvent.h"
50f3bbf4 43#include "AliESDtrack.h"
44#include "AliTrackReference.h"
b17a2d01 45#include "AliTracker.h"
50f3bbf4 46
47#include "AliTRDcluster.h"
50f3bbf4 48#include "AliTRDpadPlane.h"
49#include "AliTRDcalibDB.h"
50f3bbf4 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;
6c23ffed 263 if (esdTrack->GetOuterParam()->Pt() < 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
6c23ffed 290 double dp = aRef->Pt() + out.GetSignedPt();
50f3bbf4 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;
af26ce80 427 //fTracker->Transform(cls);
50f3bbf4 428 fClPos->Fill(cls->GetZ(), cls->GetY());
429
053767a4 430 int layer = fGeo->GetLayer(cls->GetDetector());
50f3bbf4 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
053767a4 443 AliTRDpadPlane *padPlane = fGeo->GetPadPlane(layer,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
34eaaa7e 451 if (cls->GetNPads() == 2) fClY2->Fill(10 * (yy - yref));
452 if (cls->GetNPads() == 3) fClY3->Fill(10 * (yy - yref));
50f3bbf4 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);
053767a4 460 fClYX->Fill(cls->GetY(), cls->GetX() - fGeo->GetTime0(layer));
50f3bbf4 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");
af885e0f 654 fESD = new AliESDEvent();
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
50f3bbf4 677 TClonesArray *clRefs = new TClonesArray("AliTrackReference");
50f3bbf4 678
b17a2d01 679 TBranch *branch = refTree->GetBranch("TrackReferences");
680 refTree->SetBranchAddress("TrackReferences",&clRefs);
50f3bbf4 681
b17a2d01 682 int nEntries = branch->GetEntries();
683 for(int iTrack = 0; iTrack < nEntries; iTrack++) {
50f3bbf4 684
b17a2d01 685 refTree->GetEvent(iTrack);
686 int nPoints = clRefs->GetEntries();
687 for(int iPoint=0; iPoint<nPoints; iPoint++) {
688 AliTrackReference *ref = (AliTrackReference*)clRefs->At(iPoint);
689 if (ref->DetectorId() == AliTrackReference::kTPC) fRefTPC->Add(new AliTrackReference(*ref));
690 if (ref->DetectorId() == AliTrackReference::kTRD) fRefTRD->Add(new AliTrackReference(*ref));
691 }
50f3bbf4 692 }
693
694 fRefTPC->Sort();
695 fRefTRD->Sort();
696
697 for(int i=0; i<fRefTRD->GetEntries(); i++) {
698 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
699 fLabels[i] = ref->GetTrack();
700
701 int p=0;
702 while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
703 fRefSpace->Fill(ref->LocalY(), ref->LocalX()-fGeo->GetTime0(p));
704
180da4e4 705 //for(int bit=0; bit<9; bit++) if (ref->TestBit(bit)) fTestBits->Fill(bit);
50f3bbf4 706 }
707
708 delete clRefs;
709 Info("LoadRefs", "TPC = %d\t TRD = %d", fRefTPC->GetEntries(), fRefTRD->GetEntries());
710}
711
712//////////////////////////////////////////////////////////////////////////////////////////
713
180da4e4 714Int_t AliTRDtrackingAnalysis::GetReference(Int_t label)
715{
716 //
717 // Sort the track references
718 //
50f3bbf4 719
720 int start = TMath::BinarySearch(fRefTRD->GetEntries(), fLabels, label);
721
722 while (start >= 0) {
723 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[start];
724 if (ref->GetTrack() != label) return start+1;
725 start--;
726 }
727
728 return 0;
729}
730
731//////////////////////////////////////////////////////////////////////////////////////////
732
180da4e4 733int AliTRDtrackingAnalysis::GetMCPosition(Int_t label, Double_t x, Double_t &Y, Double_t &Z, Double_t &tgphi)
734{
735 //
736 // Determine the MC positions from the track references
737 //
50f3bbf4 738
739 double lowX = 100.;
740 double highX = 100.;
741 int idLow = -1;
742 int idHigh = -1;
743
744 int nref= 0;
745 int idx = GetReference(label);
746 for(int i=idx; i<fRefTRD->GetEntries(); i++) {
747
748 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
749 if (ref->GetTrack() != label) break;
750 nref++;
751
752 //int p=0;
753 //while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
754 //if (p != layer) continue;
755
756 double dX = ref->LocalX()-x;
757 if ( dX > 0 ) {
758 if (dX < highX) {
759 idHigh = i;
760 highX = dX;
761 }
762 } else {
74cc0efb 763 dX = TMath::Abs(dX);
50f3bbf4 764 if (dX < lowX) {
765 idLow = i;
766 lowX = dX;
767 }
768 }
769 }
770
771 if (idLow == -1 || idHigh == -1) return -1;
772
773 AliTrackReference *refI = (AliTrackReference*)(*fRefTRD)[idLow];
774 AliTrackReference *refO = (AliTrackReference*)(*fRefTRD)[idHigh];
775
776
777 double dx = refO->LocalX() - refI->LocalX();
778 double dy = refO->LocalY() - refI->LocalY();
779 double dz = refO->Z() - refI->Z();
780 double ddx = (x - refI->LocalX())/dx;
781
782 fRefDx->Fill(dx);
783
784 Y = refI->LocalY() + ddx * dy;
785 Z = refI->Z() + ddx * dz;
786
787 tgphi = dy/dx;
788
789 return 0;
790}
791
792//////////////////////////////////////////////////////////////////////////////////////////
180da4e4 793Int_t AliTRDtrackingAnalysis::GetPhiBin(Double_t phi) const
794{
795 //
796 // Return the phi bin
797 //
50f3bbf4 798 return (int)((phi+0.3)/0.05);
799}
800
801//////////////////////////////////////////////////////////////////////////////////////////
180da4e4 802Double_t AliTRDtrackingAnalysis::GetPhi(Int_t bin) const
803{
804 //
805 // Return phi for a given bin
806 //
50f3bbf4 807 return bin * 0.05 - 0.3 + 0.025;
808}
809//////////////////////////////////////////////////////////////////////////////////////////