Made Getters const
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
CommitLineData
46d29e70 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/*
17$Log$
18Revision 1.1.2.2 2000/10/04 16:34:58 cblume
19Replace include files by forward declarations
20
21Revision 1.1.2.1 2000/09/22 14:47:52 cblume
22Add the tracking code
23
24*/
25
26#include <TFile.h>
27#include <TROOT.h>
28#include <TBranch.h>
29#include <TTree.h>
30
31#include "AliRun.h"
32
33#include "AliTRD.h"
34#include "AliTRDconst.h"
35#include "AliTRDgeometry.h"
36#include "AliTRDrecPoint.h"
37#include "AliTRDcluster.h"
38#include "AliTRDtrack.h"
39#include "AliTRDtrackingSector.h"
40#include "AliTRDtimeBin.h"
41
42#include "AliTRDtracker.h"
43
44ClassImp(AliTRDtracker)
45
46//____________________________________________________________________
47AliTRDtracker::AliTRDtracker()
48{
49 //
50 // Default constructor
51 //
52
53 fInputFile = NULL;
54 fEvent = 0;
55
56 fGeom = NULL;
57 fNclusters = 0;
58 fClusters = NULL;
59 fNseeds = 0;
60 fSeeds = NULL;
61 fNtracks = 0;
62 fTracks = NULL;
63
64}
65
66//____________________________________________________________________
67AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
68 :TNamed(name, title)
69{
70 fInputFile = NULL;
71 fEvent = 0;
72
73 fGeom = NULL;
74 fNclusters = 0;
75 fClusters = new TObjArray(2000);
76 fNseeds = 0;
77 fSeeds = new TObjArray(20000);
78 fNtracks = 0;
79 fTracks = new TObjArray(10000);
80
81}
82
83
84//___________________________________________________________________
85AliTRDtracker::~AliTRDtracker()
86{
87 if (fInputFile) {
88 fInputFile->Close();
89 delete fInputFile;
90 }
91 delete fClusters;
92 delete fTracks;
93 delete fSeeds;
94 delete fGeom;
95}
96
97
98//_____________________________________________________________________
99static Double_t SigmaY2trd(Double_t r, Double_t tgl, Double_t pt)
100{
101 // Parametrised "expected" error of the cluster reconstruction in Y
102
103 Double_t s = 1.;
104 return s;
105}
106
107//_____________________________________________________________________
108static Double_t SigmaZ2trd(Double_t r, Double_t tgl)
109{
110 // Parametrised "expected" error of the cluster reconstruction in Z
111
112 AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
113 AliTRDgeometry *fGeom;
114 fGeom = TRD->GetGeometry();
115 Double_t s, pad = fGeom->GetRowPadSize();
116 s = pad * pad /12.;
117 return s;
118}
119
120//_____________________________________________________________________
121inline Double_t f1trd(Double_t x1,Double_t y1,
122 Double_t x2,Double_t y2,
123 Double_t x3,Double_t y3)
124{
125 // Initial approximation of the track curvature
126 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
127
128 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
129 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
130 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
131 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
132 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
133
134 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
135
136 return -xr*yr/sqrt(xr*xr+yr*yr);
137}
138
139//_____________________________________________________________________
140inline Double_t f2trd(Double_t x1,Double_t y1,
141 Double_t x2,Double_t y2,
142 Double_t x3,Double_t y3)
143{
144 // Initial approximation of the track curvature times center of curvature
145 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
146
147 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
148 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
149 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
150 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
151 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
152
153 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
154
155 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
156}
157
158//_____________________________________________________________________
159inline Double_t f3trd(Double_t x1,Double_t y1,
160 Double_t x2,Double_t y2,
161 Double_t z1,Double_t z2)
162{
163 // Initial approximation of the tangent of the track dip angle
164 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
165
166 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
167}
168
169
170//___________________________________________________________________
171
172static Int_t FindProlongation(AliTRDseed& t, AliTRDtrackingSector *sec,
173 Int_t s, Int_t rf=0)
174{
175 // Starting from current position on track=t this function tries
176 // to extrapolate the track up to timeBin=rf and to confirm prolongation
177 // if a close cluster is found. *sec is a pointer to allocated
178 // array of sectors, in which the initial sector has index=s.
179
180 const Int_t TIME_BINS_TO_SKIP=Int_t(0.2*sec->GetNtimeBins());
181 const Double_t MAX_CHI2=12.;
182 Int_t try_again=TIME_BINS_TO_SKIP;
183
184 Double_t alpha=sec->GetAlpha();
185
186 Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
187
188 for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
189 Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
190 if (!t.PropagateTo(x)) {
191 cerr<<"Can't propagate to x = "<<x<<endl;
192 return 0;
193 }
194
195 AliTRDcluster *cl=0;
196 UInt_t index=0;
197
198 Double_t max_chi2=MAX_CHI2;
199 // const AliTRDtimeBin& time_bin=sec[s][nr];
200 AliTRDtimeBin& time_bin=sec[s][nr];
201 Double_t sy2=SigmaY2trd(t.GetX(),t.GetTgl(),t.GetPt());
202 Double_t sz2=SigmaZ2trd(t.GetX(),t.GetTgl());
203 Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
204
205 if (road>30) {
206 if (t.GetNclusters() > 4) {
207 cerr<<t.GetNclusters()<<" FindProlongation: Road is too wide !\n";
208 }
209 return 0;
210 }
211
212 if (time_bin) {
213 for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
214 AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
215 if (c->GetY() > y+road) break;
216 if (c->IsUsed() > 0) continue;
217
218 if((c->GetZ()-z)*(c->GetZ()-z) > 25. + sz2) continue;
219
220 Double_t chi2=t.GetPredictedChi2(c);
221
222 if (chi2 > max_chi2) continue;
223 max_chi2=chi2;
224 cl=c;
225 index=time_bin.GetIndex(i);
226 }
227 }
228 if (cl) {
229
230 // Float_t l=sec->GetPitch();
231 // t.SetSampledEdx(cl->fQ/l,Int_t(t));
232
233 t.Update(cl,max_chi2,index);
234
235 try_again=TIME_BINS_TO_SKIP;
236 } else {
237 if (try_again==0) break;
238 if (y > ymax) {
239 cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
240 s = (s+1) % ns;
241 if (!t.Rotate(alpha)) {
242 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
243 return 0;
244 }
245 } else if (y <-ymax) {
246 cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
247 s = (s-1+ns) % ns;
248 if (!t.Rotate(-alpha)) {
249 cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
250 return 0;
251 }
252 }
253 try_again--;
254 }
255 }
256
257 return 1;
258}
259
260
261
262//_____________________________________________________________________________
263void AliTRDtracker::GetEvent(const Char_t *name, Int_t nEvent)
264{
265 // Opens a ROOT-file with TRD-clusters and reads in the cluster-tree
266
267
268 // Connect the AliRoot file containing Geometry, Kine, and Hits
269 fInputFile = (TFile*) gROOT->GetListOfFiles()->FindObject(name);
270 if (!fInputFile) {
271 printf("AliTRDtracker::Open -- ");
272 printf("Open the ALIROOT-file %s.\n",name);
273 fInputFile = new TFile(name,"UPDATE");
274 }
275 else {
276 printf("AliTRDtracker::Open -- ");
277 printf("%s is already open.\n",name);
278 }
279
280 // Get AliRun object from file or create it if not on file
281 //if (!gAlice) {
282 gAlice = (AliRun*) fInputFile->Get("gAlice");
283 if (gAlice) {
284 printf("AliTRDtracker::GetEvent -- ");
285 printf("AliRun object found on file.\n");
286 }
287 else {
288 printf("AliTRDtracker::GetEvent -- ");
289 printf("Could not find AliRun object.\n");
290 }
291 //}
292
293 fEvent = nEvent;
294
295 // Import the Trees for the event nEvent in the file
296 Int_t nparticles = gAlice->GetEvent(fEvent);
297 cerr<<"nparticles = "<<nparticles<<endl;
298
299 if (nparticles <= 0) {
300 printf("AliTRDtracker::GetEvent -- ");
301 printf("No entries in the trees for event %d.\n",fEvent);
302 }
303
304 AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
305 fGeom = TRD->GetGeometry();
306
307 Char_t treeName[14];
308 sprintf(treeName,"TRDrecPoints%d", nEvent);
309
310 TTree *tree=(TTree*)fInputFile->Get(treeName);
311 TBranch *branch=tree->GetBranch("TRDrecPoints");
312 Int_t nentr = (Int_t) tree->GetEntries();
313 printf("found %d entries in %s tree.\n",nentr,treeName);
314
315 for (Int_t i=0; i<nentr; i++) {
316 TObjArray *ioArray = new TObjArray(400);
317 branch->SetAddress(&ioArray);
318 tree->GetEvent(i);
319 Int_t npoints = ioArray->GetEntriesFast();
320 printf("Read %d rec. points from entry %d \n", npoints, i);
321 for(Int_t j=0; j<npoints; j++) {
322 AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j);
323 AliTRDcluster *c = new AliTRDcluster(p);
324 fClusters->AddLast(c);
325 }
326 }
327
328}
329
330
331//_____________________________________________________________________________
332void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
333{
334 // Fills clusters into TRD tracking_sectors
335 // Note that the numbering scheme for the TRD tracking_sectors
336 // differs from that of TRD sectors
337
338 for (Int_t i=0; i<kNsect; i++) sec[i].SetUp();
339
340 // Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's
341
342 cerr<<"MakeSeeds: sorting clusters"<<endl;
343
344 Int_t ncl=fClusters->GetEntriesFast();
345 UInt_t index;
346 while (ncl--) {
347 printf("\r %d left",ncl);
348 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
349 Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
350 Int_t sector=fGeom->GetSector(detector);
351
352 Int_t tracking_sector=sector;
353 if(sector > 0) tracking_sector = kNsect - sector;
354
355 Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin);
356 index=ncl;
357 sec[tracking_sector][tb].InsertCluster(c,index);
358 }
359 printf("\r\n");
360}
361
362
363//_____________________________________________________________________________
364void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
365{
366 // Creates track seeds using clusters in timeBins=i1,i2
367
368 Int_t i2 = inner, i1=outer;
369
370 if (!fClusters) return;
371
372 AliTRDtrackingSector fTrSec[kNsect];
373 SetUpSectors(fTrSec);
374
375 // find seeds
376
377 Double_t x[5], c[15];
378 Int_t max_sec=kNsect;
379
380 Double_t alpha=fTrSec[0].GetAlpha();
381 Double_t shift=fTrSec[0].GetAlpha()/2.;
382 Double_t cs=cos(alpha), sn=sin(alpha);
383
384 Double_t x1 =fTrSec[0].GetX(i1);
385 Double_t xx2=fTrSec[0].GetX(i2);
386
387 for (Int_t ns=0; ns<max_sec; ns++) {
388
389 Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2];
390 Int_t nm=fTrSec[ns][i2];
391 Int_t nu=fTrSec[(ns+1)%max_sec][i2];
392
393 AliTRDtimeBin& r1=fTrSec[ns][i1];
394
395 for (Int_t is=0; is < r1; is++) {
396 Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
397
398 for (Int_t js=0; js < nl+nm+nu; js++) {
399 const AliTRDcluster *cl;
400 Double_t x2, y2, z2;
401 Double_t x3=0.,y3=0.;
402
403 if (js<nl) {
404 AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2];
405 cl=r2[js];
406 y2=cl->GetY(); z2=cl->GetZ();
407
408 x2= xx2*cs+y2*sn;
409 y2=-xx2*sn+y2*cs;
410
411 } else
412 if (js<nl+nm) {
413 AliTRDtimeBin& r2=fTrSec[ns][i2];
414 cl=r2[js-nl];
415 x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
416 } else {
417 AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
418 cl=r2[js-nl-nm];
419 y2=cl->GetY(); z2=cl->GetZ();
420
421 x2=xx2*cs-y2*sn;
422 y2=xx2*sn+y2*cs;
423
424 }
425
426 Double_t zz=z1 - z1/x1*(x1-x2);
427
428 if (TMath::Abs(zz-z2)>30.) continue;
429
430 Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
431 if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
432
433 x[0]=y1;
434 x[1]=z1;
435 x[2]=f1trd(x1,y1,x2,y2,x3,y3);
436
437 if (TMath::Abs(x[2]) >= 0.01) continue;
438
439 x[3]=f2trd(x1,y1,x2,y2,x3,y3);
440
441 if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
442
443 x[4]=f3trd(x1,y1,x2,y2,z1,z2);
444
445 if (TMath::Abs(x[4]) > 1.2) continue;
446
447 Double_t a=asin(x[3]);
448 Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
449 if (TMath::Abs(zv)>200.) continue;
450
451 Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2()*12;
452 Double_t sy2=cl->GetSigmaY2(), sz2=cl->GetSigmaZ2()*12;
453 Double_t sy3=100*0.025, sy=0.1, sz=0.1;
454
455 Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
456 Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
457 Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
458 Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
459 Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
460 Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
461 Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
462 Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
463 Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
464 Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;
465
466 c[0]=sy1;
467 c[1]=0.; c[2]=sz1;
468 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
469 c[6]=f30*sy1; c[7]=0.; c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
470 c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
471 c[10]=f40*sy1; c[11]=f41*sz1; c[12]=f40*sy1*f20+f42*sy2*f22;
472 c[13]=f40*sy1*f30+f42*sy2*f32;
473 c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
474
475 UInt_t index=r1.GetIndex(is);
476 AliTRDseed *track=new AliTRDseed(index, x, c, x1, ns*alpha+shift);
477
478 // Float_t l=fTrSec->GetPitch();
479 // track->SetSampledEdx(r1[is]->fQ/l,0);
480
481 Int_t rc=FindProlongation(*track,fTrSec,ns,i2);
482
483 if ((rc < 0) || (track->GetNclusters() < (i1-i2)/2) ) delete track;
484 else {
485 fSeeds->AddLast(track); fNseeds++;
486 cerr<<"found seed "<<fNseeds<<endl;
487 }
488 }
489 }
490 }
491
492 // delete[] fTrSec;
493
494 fSeeds->Sort();
495}
496
497//___________________________________________________________________
498void AliTRDtracker::FindTracks()
499{
500 if (!fClusters) return;
501
502 AliTRDtrackingSector fTrSec[kNsect];
503 SetUpSectors(fTrSec);
504
505 // find tracks
506
507 Int_t num_of_time_bins = fTrSec[0].GetNtimeBins();
508 Int_t nseed=fSeeds->GetEntriesFast();
509
510 Int_t nSeedClusters;
511 for (Int_t i=0; i<nseed; i++) {
512 cerr<<"FindTracks: seed "<<i+1<<" out of "<<nseed<<endl;
513
514 AliTRDseed& t=*((AliTRDseed*)fSeeds->UncheckedAt(i));
515
516 nSeedClusters = t.GetNclusters();
517 Double_t alpha=t.GetAlpha();
518
519 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
520 if (alpha < 0. ) alpha += 2.*TMath::Pi();
521 Int_t ns=Int_t(alpha/fTrSec[0].GetAlpha())%kNsect;
522
523 if (FindProlongation(t,fTrSec,ns)) {
524 cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl;
525 if ((t.GetNclusters() >= Int_t(0.3*num_of_time_bins)) &&
526 ((t.GetNclusters()-nSeedClusters)>60)) {
527 t.CookdEdx();
528 Int_t label = GetTrackLabel(t);
529 t.SetLabel(label);
530 AliTRDtrack* pt=&t;
531 fTracks->AddLast(pt); fNtracks++;
532 UseClusters(t);
533 cerr<<"found track "<<fNtracks<<endl;
534 }
535 }
536 }
537}
538
539//__________________________________________________________________
540void AliTRDtracker::UseClusters(AliTRDseed t) {
541 Int_t ncl=t.GetNclusters();
542 for (Int_t i=0; i<ncl; i++) {
543 Int_t index = t.GetClusterIndex(i);
544 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
545 c->Use();
546 }
547}
548
549//__________________________________________________________________
550Int_t AliTRDtracker::GetTrackLabel(AliTRDseed t) {
551
552 Int_t label=123456789, index, i, j;
553 Int_t ncl=t.GetNclusters();
554 const Int_t range=300;
555 Bool_t label_added;
556
557 Int_t s[range][2];
558 for (i=0; i<range; i++) {
559 s[i][0]=-1;
560 s[i][1]=0;
561 }
562
563 Int_t t0,t1,t2;
564 for (i=0; i<ncl; i++) {
565 index=t.GetClusterIndex(i);
566 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
567 t0=c->GetTrackIndex(0);
568 t1=c->GetTrackIndex(1);
569 t2=c->GetTrackIndex(2);
570 }
571
572 for (i=0; i<ncl; i++) {
573 index=t.GetClusterIndex(i);
574 AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
575 for (Int_t k=0; k<3; k++) {
576 label=c->GetTrackIndex(k);
577 label_added=false; j=0;
578 if (label >= 0) {
579 while ( (!label_added) && ( j < range ) ) {
580 if (s[j][0]==label || s[j][1]==0) {
581 s[j][0]=label;
582 s[j][1]=s[j][1]+1;
583 label_added=true;
584 }
585 j++;
586 }
587 }
588 }
589 }
590
591
592 Int_t max=0;
593 label = -123456789;
594
595 for (i=0; i<range; i++) {
596 if (s[i][1]>max) {
597 max=s[i][1]; label=s[i][0];
598 }
599 }
600 if(max > ncl/2) return label;
601 else return -1;
602}
603
604//___________________________________________________________________
605
606Int_t AliTRDtracker::WriteTracks() {
607
608 TDirectory *savedir=gDirectory;
609
610 TFile *out=TFile::Open("AliTRDtracks.root","RECREATE");
611
612 TTree tracktree("TreeT","Tree with TRD tracks");
613
614 AliTRDtrack *iotrack=0;
615 tracktree.Branch("tracks","AliTRDtrack",&iotrack,32000,0);
616
617 Int_t ntracks=fTracks->GetEntriesFast();
618
619 for (Int_t i=0; i<ntracks; i++) {
620 AliTRDtrack *pt=(AliTRDtrack*)fTracks->UncheckedAt(i);
621 iotrack=pt;
622 tracktree.Fill();
623 cerr<<"WriteTracks: put track "<<i<<" in the tree"<<endl;
624 }
625
626 tracktree.Write();
627 out->Close();
628
629 savedir->cd();
630
631 cerr<<"WriteTracks: done"<<endl;
632 return 0;
633}
634
635
636
637//_____________________________________________________________________________
638void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename, Int_t nEvent = 0, Int_t option = 1)
639{
640 //
641 // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
642 // from the file. The names of the cluster tree and branches
643 // should match the ones used in AliTRDclusterizer::WriteClusters()
644 //
645
646 TDirectory *savedir=gDirectory;
647
648 TFile *file = TFile::Open(filename);
649 if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;}
650
651 Char_t treeName[14];
652 sprintf(treeName,"TRDrecPoints%d", nEvent);
653
654 TTree *tree=(TTree*)file->Get(treeName);
655 TBranch *branch=tree->GetBranch("TRDrecPoints");
656 Int_t nentr = (Int_t) tree->GetEntries();
657 printf("found %d entries in %s tree.\n",nentr,treeName);
658
659 for (Int_t i=0; i<nentr; i++) {
660 TObjArray *ioArray = new TObjArray(400);
661 branch->SetAddress(&ioArray);
662 tree->GetEvent(i);
663 Int_t npoints = ioArray->GetEntriesFast();
664 printf("Read %d rec. points from entry %d \n", npoints, i);
665 for(Int_t j=0; j<npoints; j++) {
666 AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j);
667 AliTRDcluster *c = new AliTRDcluster(p);
668 if( option >= 0) array->AddLast(c);
669 else array->AddLast(p);
670 }
671 }
672
673 file->Close();
674 savedir->cd();
675
676}
677
678
679
680