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