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