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