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