Corrected memory leaks
[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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  The standard TRD tracker                                                 //
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include <Riostream.h>
25 #include <TFile.h>
26 #include <TBranch.h>
27 #include <TTree.h>  
28 #include <TObjArray.h> 
29
30 #include "AliTRDgeometry.h"
31 #include "AliTRDparameter.h"
32 #include "AliTRDgeometryDetail.h"
33 #include "AliTRDgeometryHole.h"
34 #include "AliTRDcluster.h" 
35 #include "AliTRDtrack.h"
36 #include "AliTRDPartID.h"
37 #include "../TPC/AliTPCtrack.h"
38
39 #include "AliTRDtracker.h"
40
41 ClassImp(AliTRDtracker) 
42
43   const  Float_t     AliTRDtracker::fgkSeedDepth          = 0.5; 
44   const  Float_t     AliTRDtracker::fgkSeedStep           = 0.10;   
45   const  Float_t     AliTRDtracker::fgkSeedGap            = 0.25;  
46
47   const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ12    = 40.;  
48   const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ      = 25.;  
49   const  Float_t     AliTRDtracker::fgkMaxSeedC           = 0.0052; 
50   const  Float_t     AliTRDtracker::fgkMaxSeedTan         = 1.2;  
51   const  Float_t     AliTRDtracker::fgkMaxSeedVertexZ     = 150.; 
52
53   const  Double_t    AliTRDtracker::fgkSeedErrorSY        = 0.2;
54   const  Double_t    AliTRDtracker::fgkSeedErrorSY3       = 2.5;
55   const  Double_t    AliTRDtracker::fgkSeedErrorSZ        = 0.1;
56
57   const  Float_t     AliTRDtracker::fgkMinClustersInSeed  = 0.7;  
58
59   const  Float_t     AliTRDtracker::fgkMinClustersInTrack = 0.5;  
60   const  Float_t     AliTRDtracker::fgkMinFractionOfFoundClusters = 0.8;  
61
62   const  Float_t     AliTRDtracker::fgkSkipDepth          = 0.3;
63   const  Float_t     AliTRDtracker::fgkLabelFraction      = 0.8;  
64   const  Float_t     AliTRDtracker::fgkWideRoad           = 20.;
65
66   const  Double_t    AliTRDtracker::fgkMaxChi2            = 12.; 
67
68 const Int_t AliTRDtracker::fgkFirstPlane = 5;
69 const Int_t AliTRDtracker::fgkLastPlane = 17;
70
71
72 //____________________________________________________________________
73 AliTRDtracker::AliTRDtracker():AliTracker(),
74                                fGeom(0),
75                                fPar(0),
76                                fNclusters(0),
77                                fClusters(0),
78                                fNseeds(0),
79                                fSeeds(0),
80                                fNtracks(0),
81                                fTracks(0),
82                                fSY2corr(0),
83                                fSZ2corr(0),
84                                fTimeBinsPerPlane(0),
85                                fMaxGap(0),
86                                fVocal(kFALSE),
87                                fAddTRDseeds(kFALSE),
88                                fNoTilt(kFALSE)
89 {
90   for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
91   for(Int_t j=0;j<5;j++)
92     for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
93
94 //____________________________________________________________________
95 AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
96 {
97   // 
98   //  Main constructor
99   //  
100
101   //Float_t fTzero = 0;
102    
103   fAddTRDseeds = kFALSE;
104   fGeom = NULL;
105   fNoTilt = kFALSE;
106   
107   TDirectory *savedir=gDirectory; 
108   TFile *in=(TFile*)geomfile;  
109   if (!in->IsOpen()) {
110     printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
111     printf("    DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
112   }
113   else {
114     in->cd();  
115 //    in->ls();
116     fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
117     fPar  = (AliTRDparameter*) in->Get("TRDparameter");
118 //    fGeom->Dump();
119   }
120
121   if(fGeom) {
122     //    fTzero = geo->GetT0();
123     printf("Found geometry version %d on file \n", fGeom->IsVersion());
124   }
125   else { 
126     printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
127     //printf("The DETAIL TRD geometry will be used\n");
128     //fGeom = new AliTRDgeometryDetail();
129     fGeom = new AliTRDgeometryHole();
130     fGeom->SetPHOShole();
131     fGeom->SetRICHhole();    
132   } 
133
134   if (!fPar) {  
135     printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n");
136     printf("The DEFAULT TRD parameter will be used\n");
137     fPar = new AliTRDparameter();
138   }
139   fPar->ReInit();
140
141   savedir->cd();  
142
143
144   //  fGeom->SetT0(fTzero);
145
146   fNclusters = 0;
147   fClusters  = new TObjArray(2000); 
148   fNseeds    = 0;
149   fSeeds     = new TObjArray(2000);
150   fNtracks   = 0;
151   fTracks    = new TObjArray(1000);
152
153   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
154     Int_t trS = CookSectorIndex(geomS);
155     fTrSec[trS] = new AliTRDtrackingSector(fGeom, geomS, fPar);
156     for (Int_t icham=0;icham<AliTRDgeometry::kNcham; icham++){
157       fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS);
158     }
159   }
160
161   Float_t tiltAngle = TMath::Abs(fPar->GetTiltingAngle()); 
162   if(tiltAngle < 0.1) {
163     fNoTilt = kTRUE;
164   }
165
166   fSY2corr = 0.2;
167   fSZ2corr = 120.;      
168
169   if(fNoTilt && (tiltAngle > 0.1)) fSY2corr = fSY2corr + tiltAngle * 0.05; 
170
171
172   // calculate max gap on track
173
174   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
175   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
176
177   Double_t dx = (Double_t) fPar->GetTimeBinSize();   
178   Int_t tbAmp = fPar->GetTimeBefore();
179   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
180   if(kTRUE) maxAmp = 0;  // intentional until we change the parameter class 
181   Int_t tbDrift = fPar->GetTimeMax();
182   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
183
184   tbDrift = TMath::Min(tbDrift,maxDrift);
185   tbAmp = TMath::Min(tbAmp,maxAmp);
186
187   fTimeBinsPerPlane = tbAmp + tbDrift;
188   fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fgkSkipDepth);
189
190   fVocal = kFALSE;
191
192
193   // Barrel Tracks [SR, 03.04.2003]
194
195   fBarrelFile = 0;
196   fBarrelTree = 0;
197   fBarrelArray = 0;
198   fBarrelTrack = 0;
199
200   savedir->cd();
201 }   
202
203 //___________________________________________________________________
204 AliTRDtracker::~AliTRDtracker()
205 {
206   //
207   // Destructor of AliTRDtracker 
208   //
209
210   if (fClusters) {
211     fClusters->Delete();
212     delete fClusters;
213   }
214   if (fTracks) {
215     fTracks->Delete();
216     delete fTracks;
217   }
218   if (fSeeds) {
219     fSeeds->Delete();
220     delete fSeeds;
221   }
222   delete fGeom;  
223   delete fPar;  
224
225   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
226     delete fTrSec[geomS];
227   }
228 }   
229
230 //_____________________________________________________________________
231
232 void AliTRDtracker::SetBarrelTree(const char *mode) {
233   //
234   //
235   //
236
237   if (!IsStoringBarrel()) return;
238
239   TDirectory *sav = gDirectory;
240   if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
241
242   char buff[40];
243   sprintf(buff,  "BarrelTRD_%d_%s", GetEventNumber(), mode);
244
245   fBarrelFile->cd();
246   fBarrelTree = new TTree(buff, "Barrel TPC tracks");
247   
248   Int_t nRefs = fgkLastPlane - fgkFirstPlane + 1;
249
250   if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", nRefs);
251   for(Int_t i=0; i<nRefs; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
252   
253   fBarrelTree->Branch("tracks", &fBarrelArray);
254   sav->cd();
255 }
256   
257 //_____________________________________________________________________
258
259 void AliTRDtracker::StoreBarrelTrack(AliTRDtrack *ps, Int_t refPlane, Int_t isIn) {
260   //
261   //
262   //
263   
264   if (!IsStoringBarrel()) return;
265   
266   static Int_t nClusters;
267   static Int_t nWrong;
268   static Double_t chi2;
269   static Int_t index;
270   static Bool_t wasLast = kTRUE;
271   
272   Int_t newClusters, newWrong;
273   Double_t newChi2;
274   
275   if (wasLast) {   
276  
277     fBarrelArray->Clear();
278     nClusters = nWrong = 0;
279     chi2 = 0.0;
280     index = 0;
281     wasLast = kFALSE;
282   }
283   
284   fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
285   ps->GetBarrelTrack(fBarrelTrack);
286   
287   newClusters = ps->GetNumberOfClusters() - nClusters; 
288   newWrong = ps->GetNWrong() - nWrong;
289   newChi2 = ps->GetChi2() - chi2;
290   
291   nClusters =  ps->GetNumberOfClusters();
292   nWrong = ps->GetNWrong();
293   chi2 = ps->GetChi2();  
294
295   if (refPlane != fgkLastPlane) {
296     fBarrelTrack->SetNClusters(newClusters, newChi2);
297     fBarrelTrack->SetNWrongClusters(newWrong);
298   } else {
299     wasLast = kTRUE;
300   } 
301
302   fBarrelTrack->SetRefPlane(refPlane, isIn);
303 }
304
305 //_____________________________________________________________________
306
307 Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
308   //
309   // Rotates the track when necessary
310   //
311
312   Double_t alpha = AliTRDgeometry::GetAlpha(); 
313   Double_t y = track->GetY();
314   Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
315
316   //Int_t ns = AliTRDgeometry::kNsect;
317   //Int_t s=Int_t(track->GetAlpha()/alpha)%ns; 
318
319   if (y > ymax) {
320     //s = (s+1) % ns;
321     if (!track->Rotate(alpha)) return kFALSE;
322   } else if (y <-ymax) {
323     //s = (s-1+ns) % ns;                           
324     if (!track->Rotate(-alpha)) return kFALSE;   
325   } 
326
327   return kTRUE;
328 }
329
330 //_____________________________________________________________________
331 inline Double_t f1trd(Double_t x1,Double_t y1,
332                       Double_t x2,Double_t y2,
333                       Double_t x3,Double_t y3)
334 {
335   //
336   // Initial approximation of the track curvature
337   //
338   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
339   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
340                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
341   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
342                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
343
344   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
345
346   return -xr*yr/sqrt(xr*xr+yr*yr);
347 }          
348
349 //_____________________________________________________________________
350 inline Double_t f2trd(Double_t x1,Double_t y1,
351                       Double_t x2,Double_t y2,
352                       Double_t x3,Double_t y3)
353 {
354   //
355   // Initial approximation of the track curvature times X coordinate
356   // of the center of curvature
357   //
358
359   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
360   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
361                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
362   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
363                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
364
365   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
366
367   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
368 }          
369
370 //_____________________________________________________________________
371 inline Double_t f3trd(Double_t x1,Double_t y1,
372                       Double_t x2,Double_t y2,
373                       Double_t z1,Double_t z2)
374 {
375   //
376   // Initial approximation of the tangent of the track dip angle
377   //
378
379   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
380 }            
381
382
383 AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin){
384   //
385   //try to find cluster in the backup list
386   //
387   AliTRDcluster * cl =0;
388   UInt_t *indexes = track->GetBackupIndexes();
389   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
390     if (indexes[i]==0) break;  
391     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
392     if (!cli) break;
393     if (cli->GetLocalTimeBin()!=timebin) continue;
394     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
395     if (iplane==plane) {
396       cl = cli;
397       break;
398     }
399   }
400   return cl;
401 }
402
403
404 Int_t  AliTRDtracker::GetLastPlane(AliTRDtrack * track){
405   //
406   //return last updated plane
407   Int_t lastplane=0;
408   UInt_t *indexes = track->GetBackupIndexes();
409   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
410     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
411     if (!cli) break;
412     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
413     if (iplane>lastplane) {
414       lastplane = iplane;
415     }
416   }
417   return lastplane;
418 }
419 //___________________________________________________________________
420 Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out)
421 {
422   //
423   // Finds tracks within the TRD. File <inp> is expected to contain seeds 
424   // at the outer part of the TRD. If <inp> is NULL, the seeds
425   // are found within the TRD if fAddTRDseeds is TRUE. 
426   // The tracks are propagated to the innermost time bin 
427   // of the TRD and stored in file <out>. 
428   //
429
430   LoadEvent();
431  
432   TDirectory *savedir=gDirectory;
433
434   char   tname[100];
435
436   if (!out->IsOpen()) {
437     cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n";
438     return 1;
439   }    
440
441   sprintf(tname,"seedTRDtoTPC_%d",GetEventNumber()); 
442   TTree tpcTree(tname,"Tree with seeds from TRD at outer TPC pad row");
443   AliTPCtrack *iotrack=0;
444   tpcTree.Branch("tracks","AliTPCtrack",&iotrack,32000,0); 
445
446   sprintf(tname,"TreeT%d_TRD",GetEventNumber());
447   TTree trdTree(tname,"TRD tracks at inner TRD time bin");
448   AliTRDtrack *iotrackTRD=0;
449   trdTree.Branch("tracks","AliTRDtrack",&iotrackTRD,32000,0);  
450
451   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
452   Float_t foundMin = fgkMinClustersInTrack * timeBins; 
453
454   if (inp) {
455      TFile *in=(TFile*)inp;
456      if (!in->IsOpen()) {
457         cerr<<
458           "AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n";
459         cerr<<" ... going for seeds finding inside the TRD\n";
460      }
461      else {
462        in->cd();
463        sprintf(tname,"TRDb_%d",GetEventNumber());  
464        TTree *seedTree=(TTree*)in->Get(tname);  
465        if (!seedTree) {
466          cerr<<"AliTRDtracker::Clusters2Tracks(): ";
467          cerr<<"can't get a tree with track seeds !\n";
468          return 3;
469        }  
470        AliTRDtrack *seed=new AliTRDtrack;
471        seedTree->SetBranchAddress("tracks",&seed);
472
473        Int_t n=(Int_t)seedTree->GetEntries();
474        for (Int_t i=0; i<n; i++) {
475          seedTree->GetEvent(i);
476          seed->ResetCovariance(); 
477          AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
478          fSeeds->AddLast(tr);
479          fNseeds++;
480        }          
481        delete seed;
482        delete seedTree;
483      }
484   }
485
486   out->cd();
487
488
489   // find tracks from loaded seeds
490
491   Int_t nseed=fSeeds->GetEntriesFast();
492   Int_t i, found = 0;
493   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
494
495   for (i=0; i<nseed; i++) {   
496     AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt; 
497     FollowProlongation(t, innerTB); 
498     if (t.GetNumberOfClusters() >= foundMin) {
499       UseClusters(&t);
500       CookLabel(pt, 1-fgkLabelFraction);
501       //      t.CookdEdx();
502     }
503     iotrackTRD = pt;
504     trdTree.Fill();
505     found++;
506 //    cout<<found<<'\r';     
507
508     if(PropagateToTPC(t)) {
509       AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
510       iotrack = tpc;
511       tpcTree.Fill();
512       delete tpc;
513     }  
514     delete fSeeds->RemoveAt(i);
515     fNseeds--;
516   }     
517
518   cout<<"Number of loaded seeds: "<<nseed<<endl;  
519   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
520
521   // after tracks from loaded seeds are found and the corresponding 
522   // clusters are used, look for additional seeds from TRD
523
524   if(fAddTRDseeds) { 
525     // Find tracks for the seeds in the TRD
526     Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
527   
528     Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
529     Int_t gap = (Int_t) (timeBins * fgkSeedGap);
530     Int_t step = (Int_t) (timeBins * fgkSeedStep);
531   
532     // make a first turn with tight cut on initial curvature
533     for(Int_t turn = 1; turn <= 2; turn++) {
534       if(turn == 2) {
535         nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep));
536         step = (Int_t) (timeBins * (3*fgkSeedStep));
537       }
538       for(Int_t i=0; i<nSteps; i++) {
539         Int_t outer=timeBins-1-i*step; 
540         Int_t inner=outer-gap;
541
542         nseed=fSeeds->GetEntriesFast();
543       
544         MakeSeeds(inner, outer, turn);
545       
546         nseed=fSeeds->GetEntriesFast();
547         printf("\n turn %d, step %d: number of seeds for TRD inward %d\n", 
548                turn, i, nseed); 
549               
550         for (Int_t i=0; i<nseed; i++) {   
551           AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt; 
552           FollowProlongation(t,innerTB); 
553           if (t.GetNumberOfClusters() >= foundMin) {
554             UseClusters(&t);
555             CookLabel(pt, 1-fgkLabelFraction);
556             t.CookdEdx();
557             found++;
558 //            cout<<found<<'\r';     
559             iotrackTRD = pt;
560             trdTree.Fill();
561             if(PropagateToTPC(t)) {
562               AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
563               iotrack = tpc;
564               tpcTree.Fill();
565               delete tpc;
566             }        
567           }
568           delete fSeeds->RemoveAt(i);
569           fNseeds--;
570         }
571       }
572     }
573   }
574   tpcTree.Write(); 
575   trdTree.Write(); 
576   
577   cout<<"Total number of found tracks: "<<found<<endl;
578     
579   UnloadEvent();
580     
581   savedir->cd();  
582   
583   return 0;    
584 }     
585           
586 //___________________________________________________________________
587 Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
588 {
589   //
590   // Finds tracks within the TRD. The ESD event is expected to contain seeds 
591   // at the outer part of the TRD. The seeds
592   // are found within the TRD if fAddTRDseeds is TRUE. 
593   // The tracks are propagated to the innermost time bin 
594   // of the TRD and the ESD event is updated
595   //
596
597   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
598   Float_t foundMin = fgkMinClustersInTrack * timeBins; 
599   Int_t nseed = 0;
600   Int_t found = 0;
601   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
602
603   Int_t n = event->GetNumberOfTracks();
604   for (Int_t i=0; i<n; i++) {
605     AliESDtrack* seed=event->GetTrack(i);
606     ULong_t status=seed->GetStatus();
607     if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
608     if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
609     nseed++;
610
611     AliTRDtrack* seed2 = new AliTRDtrack(*seed);
612     //seed2->ResetCovariance(); 
613     AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
614     AliTRDtrack &t=*pt; 
615     FollowProlongation(t, innerTB); 
616     if (t.GetNumberOfClusters() >= foundMin) {
617       UseClusters(&t);
618       CookLabel(pt, 1-fgkLabelFraction);
619       //      t.CookdEdx();
620     }
621     found++;
622 //    cout<<found<<'\r';     
623
624     if(PropagateToTPC(t)) {
625       seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
626     }  
627     delete seed2;
628     delete pt;
629   }     
630
631   cout<<"Number of loaded seeds: "<<nseed<<endl;  
632   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
633
634   // after tracks from loaded seeds are found and the corresponding 
635   // clusters are used, look for additional seeds from TRD
636
637   if(fAddTRDseeds) { 
638     // Find tracks for the seeds in the TRD
639     Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
640   
641     Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
642     Int_t gap = (Int_t) (timeBins * fgkSeedGap);
643     Int_t step = (Int_t) (timeBins * fgkSeedStep);
644   
645     // make a first turn with tight cut on initial curvature
646     for(Int_t turn = 1; turn <= 2; turn++) {
647       if(turn == 2) {
648         nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep));
649         step = (Int_t) (timeBins * (3*fgkSeedStep));
650       }
651       for(Int_t i=0; i<nSteps; i++) {
652         Int_t outer=timeBins-1-i*step; 
653         Int_t inner=outer-gap;
654
655         nseed=fSeeds->GetEntriesFast();
656       
657         MakeSeeds(inner, outer, turn);
658       
659         nseed=fSeeds->GetEntriesFast();
660         printf("\n turn %d, step %d: number of seeds for TRD inward %d\n", 
661                turn, i, nseed); 
662               
663         for (Int_t i=0; i<nseed; i++) {   
664           AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt; 
665           FollowProlongation(t,innerTB); 
666           if (t.GetNumberOfClusters() >= foundMin) {
667             UseClusters(&t);
668             CookLabel(pt, 1-fgkLabelFraction);
669             t.CookdEdx();
670             found++;
671 //            cout<<found<<'\r';     
672             if(PropagateToTPC(t)) {
673               AliESDtrack track;
674               track.UpdateTrackParams(pt,AliESDtrack::kTRDin);
675               event->AddTrack(&track);
676             }        
677           }
678           delete fSeeds->RemoveAt(i);
679           fNseeds--;
680         }
681       }
682     }
683   }
684   
685   cout<<"Total number of found tracks: "<<found<<endl;
686     
687   return 0;    
688 }     
689      
690   
691
692 //_____________________________________________________________________________
693 Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
694   //
695   // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and
696   // backpropagated by the TPC tracker. Each seed is first propagated 
697   // to the TRD, and then its prolongation is searched in the TRD.
698   // If sufficiently long continuation of the track is found in the TRD
699   // the track is updated, otherwise it's stored as originaly defined 
700   // by the TPC tracker.   
701   //  
702
703   LoadEvent();
704
705   TDirectory *savedir=gDirectory;
706
707   TFile *in=(TFile*)inp;
708
709   if (!in->IsOpen()) {
710      cerr<<"AliTRDtracker::PropagateBack(): ";
711      cerr<<"file with back propagated TPC tracks is not open !\n";
712      return 1;
713   }                   
714
715   if (!out->IsOpen()) {
716      cerr<<"AliTRDtracker::PropagateBack(): ";
717      cerr<<"file for back propagated TRD tracks is not open !\n";
718      return 2;
719   }      
720
721   in->cd();
722   char   tname[100];
723   sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber());       
724   TTree *seedTree=(TTree*)in->Get(tname);
725   if (!seedTree) {
726      cerr<<"AliTRDtracker::PropagateBack(): ";
727      cerr<<"can't get a tree with seeds from TPC !\n";
728      cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n";
729      return 3;
730   }
731
732   AliTPCtrack *seed=new AliTPCtrack;
733   seedTree->SetBranchAddress("tracks",&seed);
734
735   Int_t n=(Int_t)seedTree->GetEntries();
736   for (Int_t i=0; i<n; i++) {
737      seedTree->GetEvent(i);
738      Int_t lbl = seed->GetLabel();
739      AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
740      tr->SetSeedLabel(lbl);
741      fSeeds->AddLast(tr);
742      fNseeds++;
743   }
744
745   delete seed;
746   delete seedTree;
747
748   out->cd();
749
750   AliTPCtrack *otrack=0;
751
752   sprintf(tname,"seedsTRDtoTOF1_%d",GetEventNumber());  
753   TTree tofTree1(tname,"Tracks back propagated through TPC and TRD");
754   tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
755
756   sprintf(tname,"seedsTRDtoTOF2_%d",GetEventNumber());  
757   TTree tofTree2(tname,"Tracks back propagated through TPC and TRD");
758   tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
759
760   sprintf(tname,"seedsTRDtoPHOS_%d",GetEventNumber());  
761   TTree phosTree(tname,"Tracks back propagated through TPC and TRD");
762   phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
763
764   sprintf(tname,"seedsTRDtoRICH_%d",GetEventNumber());  
765   TTree richTree(tname,"Tracks back propagated through TPC and TRD");
766   richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
767
768   sprintf(tname,"TRDb_%d",GetEventNumber());  
769   TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin");
770   AliTRDtrack *otrackTRD=0;
771   trdTree.Branch("tracks","AliTRDtrack",&otrackTRD,32000,0);   
772      
773   if (IsStoringBarrel()) SetBarrelTree("back");
774   out->cd();
775
776   Int_t found=0;  
777   Int_t nseed=fSeeds->GetEntriesFast();
778
779   //  Float_t foundMin = fgkMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan(); 
780   Float_t foundMin = 40;
781
782   Int_t outermostTB  = fTrSec[0]->GetOuterTimeBin();
783
784   for (Int_t i=0; i<nseed; i++) {  
785
786     AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
787     Int_t expectedClr = FollowBackProlongation(s);
788
789     if (IsStoringBarrel()) {
790       StoreBarrelTrack(ps, fgkLastPlane, kTrackBack);
791       fBarrelTree->Fill();        
792     }
793
794     Int_t foundClr = s.GetNumberOfClusters();
795     Int_t lastTB = fTrSec[0]->GetLayerNumber(s.GetX());
796
797     //    printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
798     //     i, foundClr, expectedClr, foundMin);
799
800     if (foundClr >= foundMin) {
801       if(foundClr >= 2) {
802         s.CookdEdx(); 
803         CookLabel(ps, 1-fgkLabelFraction);
804         UseClusters(ps);
805       }
806       
807       // Propagate to outer reference plane [SR, GSI, 18.02.2003]
808       ps->PropagateTo(364.8);
809       otrackTRD=ps;
810       trdTree.Fill();
811       found++;
812 //      cout<<found<<'\r';
813     }
814
815     if(((expectedClr < 10) && (lastTB == outermostTB)) ||
816        ((expectedClr >= 10) && 
817         (((Float_t) foundClr) / ((Float_t) expectedClr) >= 
818          fgkMinFractionOfFoundClusters) && (lastTB == outermostTB))) {
819
820       Double_t xTOF = 375.5;
821     
822       if(PropagateToOuterPlane(s,xTOF)) {
823         AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
824         otrack = pt;
825         tofTree1.Fill();
826         delete pt;
827
828         xTOF = 381.5;
829     
830         if(PropagateToOuterPlane(s,xTOF)) {
831           AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
832           otrack = pt;
833           tofTree2.Fill();
834           delete pt;
835
836           Double_t xPHOS = 460.;
837           
838           if(PropagateToOuterPlane(s,xPHOS)) {
839             AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
840             otrack = pt;
841             phosTree.Fill();
842             delete pt;
843             
844             Double_t xRICH = 490+1.267;
845             
846             if(PropagateToOuterPlane(s,xRICH)) {
847               AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
848               otrack = pt;
849               richTree.Fill();
850               delete pt;
851             }   
852           }
853         }
854       }      
855     }
856   }
857   
858   
859   out->cd();
860   tofTree1.Write(); 
861   tofTree2.Write(); 
862   phosTree.Write(); 
863   richTree.Write(); 
864   trdTree.Write(); 
865
866   if (IsStoringBarrel()) { // [SR, 03.04.2003]
867     fBarrelFile->cd();
868     fBarrelTree->Write();
869     fBarrelFile->Flush();
870   }
871
872   savedir->cd();  
873   cerr<<"Number of seeds: "<<nseed<<endl;  
874   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
875
876   UnloadEvent();
877
878   return 0;
879
880 }
881
882 //_____________________________________________________________________________
883 Int_t AliTRDtracker::PropagateBack(AliESD* event) {
884   //
885   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
886   // backpropagated by the TPC tracker. Each seed is first propagated 
887   // to the TRD, and then its prolongation is searched in the TRD.
888   // If sufficiently long continuation of the track is found in the TRD
889   // the track is updated, otherwise it's stored as originaly defined 
890   // by the TPC tracker.   
891   //  
892
893   Int_t found=0;  
894   Float_t foundMin = 40;
895
896   Int_t n = event->GetNumberOfTracks();
897   for (Int_t i=0; i<n; i++) {
898     AliESDtrack* seed=event->GetTrack(i);
899     ULong_t status=seed->GetStatus();
900     if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
901     if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
902
903     Int_t lbl = seed->GetLabel();
904     AliTRDtrack *track = new AliTRDtrack(*seed);
905     track->SetSeedLabel(lbl);
906     fNseeds++;
907
908     Int_t expectedClr = FollowBackProlongation(*track);
909     if (track->GetNumberOfClusters()<expectedClr/3){
910       AliTRDtrack *track1 = new AliTRDtrack(*seed);
911       track1->SetSeedLabel(lbl);
912       FollowBackProlongation(*track1);
913       AliTRDtrack *track2= new AliTRDtrack(*seed);
914       track->SetSeedLabel(lbl);
915       FollowBackProlongation(*track2);      
916       delete track1;
917       delete track2;
918     }
919     
920     Int_t foundClr = track->GetNumberOfClusters();
921     if (foundClr >= foundMin) {
922       if(foundClr >= 2) {
923         track->CookdEdx(); 
924 //      CookLabel(track, 1-fgkLabelFraction);
925         UseClusters(track);
926       }
927       
928       // Propagate to outer reference plane [SR, GSI, 18.02.2003]
929 //      track->PropagateTo(364.8);  why?
930       
931       //seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
932       //found++;
933     }
934
935     //Propagation to the TOF (I.Belikov)
936     
937     if (track->GetStop()==kFALSE){
938       Double_t xTOF = 375.5;    
939       PropagateToOuterPlane(*track,xTOF); 
940       
941       Double_t xtof=378.;
942       Double_t c2=track->GetC()*xtof - track->GetEta();
943       if (TMath::Abs(c2)>=0.9999999) continue;
944       
945       Double_t ymax=xtof*TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
946       Double_t y=track->GetYat(xtof);
947       if (y > ymax) {
948         if (!track->Rotate(AliTRDgeometry::GetAlpha())) return 1;
949       } else if (y <-ymax) {
950         if (!track->Rotate(-AliTRDgeometry::GetAlpha())) return 1;
951       }
952       
953       if (track->PropagateTo(xtof)) {
954         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);    
955         if (track->GetNumberOfClusters()>foundMin) found++;
956       }
957     }else{
958       if (track->GetNumberOfClusters()>15&&track->GetNumberOfClusters()>0.5*expectedClr){
959         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
960         seed->UpdateTrackParams(track, AliESDtrack::kTRDStop);    
961         found++;
962       }
963     }
964
965     delete track;
966     
967     //End of propagation to the TOF
968     //if (foundClr>foundMin)
969     //  seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
970     
971
972   }
973   
974   cerr<<"Number of seeds: "<<fNseeds<<endl;  
975   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
976
977   fSeeds->Clear(); fNseeds=0;
978
979   return 0;
980
981 }
982
983 //_____________________________________________________________________________
984 Int_t AliTRDtracker::RefitInward(AliESD* event)
985 {
986   //
987   // Refits tracks within the TRD. The ESD event is expected to contain seeds 
988   // at the outer part of the TRD. 
989   // The tracks are propagated to the innermost time bin 
990   // of the TRD and the ESD event is updated
991   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
992   //
993
994   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
995   Float_t foundMin = fgkMinClustersInTrack * timeBins; 
996   Int_t nseed = 0;
997   Int_t found = 0;
998   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
999
1000   Int_t n = event->GetNumberOfTracks();
1001   for (Int_t i=0; i<n; i++) {
1002     AliESDtrack* seed=event->GetTrack(i);
1003     ULong_t status=seed->GetStatus();
1004     if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
1005     if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
1006     nseed++;
1007
1008     AliTRDtrack* seed2 = new AliTRDtrack(*seed);
1009     seed2->ResetCovariance(5.); 
1010     AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
1011     UInt_t * indexes2 = seed2->GetIndexes();
1012     UInt_t * indexes3 = pt->GetBackupIndexes();
1013     for (Int_t i=0;i<200;i++) {
1014       if (indexes2[i]==0) break;
1015       indexes3[i] = indexes2[i];
1016     }          
1017     //AliTRDtrack *pt = seed2;
1018     AliTRDtrack &t=*pt; 
1019     FollowProlongation(t, innerTB); 
1020
1021     if (t.GetNumberOfClusters()<seed->GetTRDclusters(indexes3)*0.5){
1022       // debug  - why we dont go back?
1023       AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
1024       UInt_t * indexes2 = seed2->GetIndexes();
1025       UInt_t * indexes3 = pt2->GetBackupIndexes();
1026       for (Int_t i=0;i<200;i++) {
1027         if (indexes2[i]==0) break;
1028         indexes3[i] = indexes2[i];
1029       }  
1030       FollowProlongation(*pt2, innerTB);
1031       delete pt2;
1032     }
1033
1034     if (t.GetNumberOfClusters() >= foundMin) {
1035       //      UseClusters(&t);
1036       //CookLabel(pt, 1-fgkLabelFraction);
1037       //      t.CookdEdx();
1038     }
1039     found++;
1040 //    cout<<found<<'\r';     
1041
1042     if(PropagateToTPC(t)) {
1043       seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
1044     }  
1045     delete seed2;
1046     delete pt;
1047   }     
1048
1049   cout<<"Number of loaded seeds: "<<nseed<<endl;  
1050   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
1051
1052   return 0;
1053
1054 }
1055
1056
1057 //---------------------------------------------------------------------------
1058 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
1059 {
1060   // Starting from current position on track=t this function tries
1061   // to extrapolate the track up to timeBin=0 and to confirm prolongation
1062   // if a close cluster is found. Returns the number of clusters
1063   // expected to be found in sensitive layers
1064
1065   Float_t  wIndex, wTB, wChi2;
1066   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
1067   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
1068   Float_t  wPx, wPy, wPz, wC;
1069   Double_t px, py, pz;
1070   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
1071   Int_t lastplane = GetLastPlane(&t);
1072
1073   Int_t trackIndex = t.GetLabel();  
1074
1075   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
1076
1077   Int_t tryAgain=fMaxGap;
1078
1079   Double_t alpha=t.GetAlpha();
1080   alpha = TVector2::Phi_0_2pi(alpha);
1081
1082   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
1083   Double_t radLength, rho, x, dx, y, ymax, z;
1084
1085   Int_t expectedNumberOfClusters = 0;
1086   Bool_t lookForCluster;
1087
1088   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1089
1090  
1091   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) { 
1092
1093     y = t.GetY(); z = t.GetZ();
1094
1095     // first propagate to the inner surface of the current time bin 
1096     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1097     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1098     if(!t.PropagateTo(x,radLength,rho)) break;
1099     y = t.GetY();
1100     ymax = x*TMath::Tan(0.5*alpha);
1101     if (y > ymax) {
1102       s = (s+1) % ns;
1103       if (!t.Rotate(alpha)) break;
1104       if(!t.PropagateTo(x,radLength,rho)) break;
1105     } else if (y <-ymax) {
1106       s = (s-1+ns) % ns;                           
1107       if (!t.Rotate(-alpha)) break;   
1108       if(!t.PropagateTo(x,radLength,rho)) break;
1109     } 
1110
1111     y = t.GetY(); z = t.GetZ();
1112
1113     // now propagate to the middle plane of the next time bin 
1114     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1115     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1116     if(!t.PropagateTo(x,radLength,rho)) break;
1117     y = t.GetY();
1118     ymax = x*TMath::Tan(0.5*alpha);
1119     if (y > ymax) {
1120       s = (s+1) % ns;
1121       if (!t.Rotate(alpha)) break;
1122       if(!t.PropagateTo(x,radLength,rho)) break;
1123     } else if (y <-ymax) {
1124       s = (s-1+ns) % ns;                           
1125       if (!t.Rotate(-alpha)) break;   
1126       if(!t.PropagateTo(x,radLength,rho)) break;
1127     } 
1128
1129
1130     if(lookForCluster) {
1131
1132       expectedNumberOfClusters++;       
1133       wIndex = (Float_t) t.GetLabel();
1134       wTB = nr;
1135
1136       AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr-1));
1137
1138       Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
1139       Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
1140
1141       Double_t road;
1142       if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
1143       else return expectedNumberOfClusters;
1144       
1145       wYrt = (Float_t) y;
1146       wZrt = (Float_t) z;
1147       wYwindow = (Float_t) road;
1148       t.GetPxPyPz(px,py,pz);
1149       wPx = (Float_t) px;
1150       wPy = (Float_t) py;
1151       wPz = (Float_t) pz;
1152       wC  = (Float_t) t.GetC();
1153       wSigmaC2 = (Float_t) t.GetSigmaC2();
1154       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
1155       wSigmaY2 = (Float_t) t.GetSigmaY2();
1156       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
1157       wChi2 = -1;            
1158       
1159
1160       AliTRDcluster *cl=0;
1161       UInt_t index=0;
1162
1163       Double_t maxChi2=fgkMaxChi2;
1164
1165       wYclosest = 12345678;
1166       wYcorrect = 12345678;
1167       wZclosest = 12345678;
1168       wZcorrect = 12345678;
1169       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);   
1170
1171       // Find the closest correct cluster for debugging purposes
1172       if (timeBin) {
1173         Float_t minDY = 1000000;
1174         for (Int_t i=0; i<timeBin; i++) {
1175           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1176           if((c->GetLabel(0) != trackIndex) &&
1177              (c->GetLabel(1) != trackIndex) &&
1178              (c->GetLabel(2) != trackIndex)) continue;
1179           if(TMath::Abs(c->GetY() - y) > minDY) continue;
1180           minDY = TMath::Abs(c->GetY() - y);
1181           wYcorrect = c->GetY();
1182           wZcorrect = c->GetZ();
1183
1184           Double_t h01 = GetTiltFactor(c);
1185           wChi2 = t.GetPredictedChi2(c, h01);
1186         }
1187       }                    
1188
1189       // Now go for the real cluster search
1190
1191       if (timeBin) {
1192         //
1193         //find cluster in history
1194         cl =0;
1195         
1196         AliTRDcluster * cl0 = timeBin[0];
1197         if (!cl0) {
1198           continue;
1199         }
1200         Int_t plane = fGeom->GetPlane(cl0->GetDetector());
1201         if (plane>lastplane) continue;
1202         Int_t timebin = cl0->GetLocalTimeBin();
1203         AliTRDcluster * cl2= GetCluster(&t,plane, timebin);
1204         if (cl2) {
1205           cl =cl2;      
1206           Double_t h01 = GetTiltFactor(cl);
1207           maxChi2=t.GetPredictedChi2(cl,h01);
1208         }
1209         if ((!cl) && road>fgkWideRoad) {
1210           //if (t.GetNumberOfClusters()>4)
1211           //  cerr<<t.GetNumberOfClusters()
1212           //    <<"FindProlongation warning: Too broad road !\n";
1213           continue;
1214         }             
1215
1216
1217         if(!cl){
1218
1219           for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1220             AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1221             if (c->GetY() > y+road) break;
1222             if (c->IsUsed() > 0) continue;
1223             if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1224             
1225             Double_t h01 = GetTiltFactor(c);
1226             Double_t chi2=t.GetPredictedChi2(c,h01);
1227             
1228             if (chi2 > maxChi2) continue;
1229             maxChi2=chi2;
1230             cl=c;
1231             index=timeBin.GetIndex(i);
1232           }               
1233         }
1234
1235         if(!cl) {
1236
1237           for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1238             AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1239             
1240             if (c->GetY() > y+road) break;
1241             if (c->IsUsed() > 0) continue;
1242             if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue;
1243             
1244             Double_t h01 = GetTiltFactor(c);
1245             Double_t chi2=t.GetPredictedChi2(c, h01);
1246             
1247             if (chi2 > maxChi2) continue;
1248             maxChi2=chi2;
1249             cl=c;
1250             index=timeBin.GetIndex(i);
1251           }
1252         }        
1253         if (cl) {
1254           wYclosest = cl->GetY();
1255           wZclosest = cl->GetZ();
1256           Double_t h01 = GetTiltFactor(cl);
1257
1258           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); 
1259           //printf("Track   position\t%f\t%f\t%f\n",t.GetX(),t.GetY(),t.GetZ());
1260           //printf("Cluster position\t%d\t%f\t%f\n",cl->GetLocalTimeBin(),cl->GetY(),cl->GetZ());
1261           Int_t det = cl->GetDetector();    
1262           Int_t plane = fGeom->GetPlane(det);
1263
1264           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1265             //if(!t.Update(cl,maxChi2,index,h01)) {
1266             //if(!tryAgain--) return 0;
1267           }  
1268           else tryAgain=fMaxGap;
1269         }
1270         else {
1271           //if (tryAgain==0) break; 
1272           tryAgain--;
1273         }
1274
1275         /*
1276         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1277           
1278           printf(" %f", wIndex);       //1
1279           printf(" %f", wTB);          //2
1280           printf(" %f", wYrt);         //3
1281           printf(" %f", wYclosest);    //4
1282           printf(" %f", wYcorrect);    //5
1283           printf(" %f", wYwindow);     //6
1284           printf(" %f", wZrt);         //7
1285           printf(" %f", wZclosest);    //8
1286           printf(" %f", wZcorrect);    //9
1287           printf(" %f", wZwindow);     //10
1288           printf(" %f", wPx);          //11
1289           printf(" %f", wPy);          //12
1290           printf(" %f", wPz);          //13
1291           printf(" %f", wSigmaC2*1000000);  //14
1292           printf(" %f", wSigmaTgl2*1000);   //15
1293           printf(" %f", wSigmaY2);     //16
1294           //      printf(" %f", wSigmaZ2);     //17
1295           printf(" %f", wChi2);     //17
1296           printf(" %f", wC);           //18
1297           printf("\n");
1298         } 
1299         */                        
1300       }
1301     }  
1302   }
1303   return expectedNumberOfClusters;
1304   
1305   
1306 }                
1307
1308 //___________________________________________________________________
1309
1310 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
1311 {
1312   // Starting from current radial position of track <t> this function
1313   // extrapolates the track up to outer timebin and in the sensitive
1314   // layers confirms prolongation if a close cluster is found. 
1315   // Returns the number of clusters expected to be found in sensitive layers
1316
1317   Float_t  wIndex, wTB, wChi2;
1318   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
1319   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
1320   Float_t  wPx, wPy, wPz, wC;
1321   Double_t px, py, pz;
1322   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
1323
1324   Int_t trackIndex = t.GetLabel();  
1325   Int_t tryAgain=fMaxGap;
1326
1327   Double_t alpha=t.GetAlpha();
1328   TVector2::Phi_0_2pi(alpha);
1329
1330   Int_t s;
1331
1332   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
1333   Double_t radLength, rho, x, dx, y, ymax = 0, z;
1334   Bool_t lookForCluster;
1335
1336   Int_t expectedNumberOfClusters = 0;
1337   x = t.GetX();
1338
1339   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1340
1341   Int_t nRefPlane = fgkFirstPlane;
1342   Bool_t isNewLayer = kFALSE; 
1343
1344   Double_t chi2;
1345   Double_t minDY;
1346   Int_t zone =-10;
1347   Int_t nr;
1348   for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) { 
1349     
1350     y = t.GetY(); 
1351     z = t.GetZ();
1352
1353     // first propagate to the outer surface of the current time bin 
1354
1355     s = t.GetSector();
1356     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1357     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; 
1358     y = t.GetY(); 
1359     z = t.GetZ();
1360
1361     if(!t.PropagateTo(x,radLength,rho)) break;
1362     //    if (!AdjustSector(&t)) break;
1363     //
1364     // MI -fix untill correct material desription will be implemented
1365     //
1366     Float_t angle =  t.GetAlpha();  // MI - if rotation - we go through the material - 
1367     if (!AdjustSector(&t)) break;
1368     if (TMath::Abs(angle -  t.GetAlpha())>0.000001) break; //better to stop track
1369     Int_t currentzone = fTrSec[s]->GetLayer(nr)->GetZone(z);
1370     if (currentzone==-10) break;  // we are in the frame
1371     if (currentzone>-10){   // layer knows where we are
1372       if (zone==-10) zone = currentzone;
1373       if (zone!=currentzone) break;  
1374     }
1375     //
1376     //
1377     s = t.GetSector();
1378     if (!t.PropagateTo(x,radLength,rho)) break;
1379
1380     y = t.GetY();
1381     z = t.GetZ();
1382
1383     // Barrel Tracks [SR, 04.04.2003]
1384
1385     s = t.GetSector();
1386     if (fTrSec[s]->GetLayer(nr)->IsSensitive() != 
1387         fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) {
1388
1389 //      if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
1390     }
1391
1392     if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() && 
1393           ! fTrSec[s]->GetLayer(nr)->IsSensitive()) {
1394       isNewLayer = kTRUE;
1395     } else {isNewLayer = kFALSE;}
1396
1397     y = t.GetY();
1398     z = t.GetZ();
1399
1400     // now propagate to the middle plane of the next time bin 
1401     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1402
1403     x = fTrSec[s]->GetLayer(nr+1)->GetX(); 
1404       if(!t.PropagateTo(x,radLength,rho)) break;
1405     if (!AdjustSector(&t)) break;
1406     s = t.GetSector();
1407       if(!t.PropagateTo(x,radLength,rho)) break;
1408
1409     y = t.GetY();
1410     z = t.GetZ();
1411
1412     if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
1413     //    printf("label %d, pl %d, lookForCluster %d \n",
1414     //     trackIndex, nr+1, lookForCluster);
1415
1416     if(lookForCluster) {
1417       expectedNumberOfClusters++;       
1418
1419       wIndex = (Float_t) t.GetLabel();
1420       wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
1421
1422       AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr+1));
1423       Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
1424       Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
1425       if((t.GetSigmaY2() + sy2) < 0) break;
1426       Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2); 
1427       Double_t y=t.GetY(), z=t.GetZ();
1428
1429       wYrt = (Float_t) y;
1430       wZrt = (Float_t) z;
1431       wYwindow = (Float_t) road;
1432       t.GetPxPyPz(px,py,pz);
1433       wPx = (Float_t) px;
1434       wPy = (Float_t) py;
1435       wPz = (Float_t) pz;
1436       wC  = (Float_t) t.GetC();
1437       wSigmaC2 = (Float_t) t.GetSigmaC2();
1438       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
1439       wSigmaY2 = (Float_t) t.GetSigmaY2();
1440       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
1441       wChi2 = -1;            
1442       
1443       if (road>fgkWideRoad) {
1444         if (t.GetNumberOfClusters()>4)
1445           cerr<<t.GetNumberOfClusters()
1446               <<"FindProlongation warning: Too broad road !\n";
1447         return 0;
1448       }      
1449
1450       AliTRDcluster *cl=0;
1451       UInt_t index=0;
1452
1453       Double_t maxChi2=fgkMaxChi2;
1454
1455       if (isNewLayer) { 
1456         road = 3 * road;
1457         //sz2 = 3 * sz2;
1458         maxChi2 = 10 * fgkMaxChi2;
1459       }
1460       
1461       if (nRefPlane == fgkFirstPlane) maxChi2 = 20 * fgkMaxChi2; 
1462       if (nRefPlane == fgkFirstPlane+2) maxChi2 = 15 * fgkMaxChi2;
1463       if (t.GetNRotate() > 0) maxChi2 = 3 * maxChi2;
1464       
1465
1466       wYclosest = 12345678;
1467       wYcorrect = 12345678;
1468       wZclosest = 12345678;
1469       wZcorrect = 12345678;
1470       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);   
1471
1472       // Find the closest correct cluster for debugging purposes
1473       if (timeBin) {
1474         minDY = 1000000;
1475         for (Int_t i=0; i<timeBin; i++) {
1476           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1477           if((c->GetLabel(0) != trackIndex) &&
1478              (c->GetLabel(1) != trackIndex) &&
1479              (c->GetLabel(2) != trackIndex)) continue;
1480           if(TMath::Abs(c->GetY() - y) > minDY) continue;
1481           //minDY = TMath::Abs(c->GetY() - y);
1482           minDY = c->GetY() - y;
1483           wYcorrect = c->GetY();
1484           wZcorrect = c->GetZ();
1485
1486           Double_t h01 = GetTiltFactor(c);
1487           wChi2 = t.GetPredictedChi2(c, h01);
1488         }
1489       }                    
1490
1491       // Now go for the real cluster search
1492
1493       if (timeBin) {
1494
1495         for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1496           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1497           if (c->GetY() > y+road) break;
1498           if (c->IsUsed() > 0) continue;
1499           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1500
1501           Double_t h01 = GetTiltFactor(c);
1502           chi2=t.GetPredictedChi2(c,h01);
1503           
1504           if (chi2 > maxChi2) continue;
1505           maxChi2=chi2;
1506           cl=c;
1507           index=timeBin.GetIndex(i);
1508
1509           //check is correct
1510           if((c->GetLabel(0) != trackIndex) &&
1511              (c->GetLabel(1) != trackIndex) &&
1512              (c->GetLabel(2) != trackIndex)) t.AddNWrong();
1513         }               
1514         
1515         if(!cl) {
1516
1517           for (Int_t i=timeBin.Find(y-road); i<timeBin; i++) {
1518             AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
1519             
1520             if (c->GetY() > y+road) break;
1521             if (c->IsUsed() > 0) continue;
1522             if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
1523             
1524             Double_t h01 = GetTiltFactor(c);
1525             chi2=t.GetPredictedChi2(c,h01);
1526             
1527             if (chi2 > maxChi2) continue;
1528             maxChi2=chi2;
1529             cl=c;
1530             index=timeBin.GetIndex(i);
1531           }
1532         }        
1533         
1534         if (cl) {
1535           wYclosest = cl->GetY();
1536           wZclosest = cl->GetZ();
1537
1538           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); 
1539           Double_t h01 = GetTiltFactor(cl);
1540           Int_t det = cl->GetDetector();    
1541           Int_t plane = fGeom->GetPlane(det);
1542
1543           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1544           //if(!t.Update(cl,maxChi2,index,h01)) {
1545             if(!tryAgain--) return 0;
1546           }  
1547           else tryAgain=fMaxGap;
1548         }
1549         else {
1550           if (tryAgain==0) break; 
1551           tryAgain--;
1552           
1553           //if (minDY < 1000000 && isNewLayer) 
1554             //cout << "\t" << nRefPlane << "\t" << "\t" << t.GetNRotate() <<  "\t" << 
1555             //  road << "\t" << minDY << "\t" << chi2 << "\t" << wChi2 << "\t" << maxChi2 << endl;
1556                                                                      
1557         }
1558
1559         isNewLayer = kFALSE;
1560
1561         /*
1562         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1563           
1564           printf(" %f", wIndex);       //1
1565           printf(" %f", wTB);          //2
1566           printf(" %f", wYrt);         //3
1567           printf(" %f", wYclosest);    //4
1568           printf(" %f", wYcorrect);    //5
1569           printf(" %f", wYwindow);     //6
1570           printf(" %f", wZrt);         //7
1571           printf(" %f", wZclosest);    //8
1572           printf(" %f", wZcorrect);    //9
1573           printf(" %f", wZwindow);     //10
1574           printf(" %f", wPx);          //11
1575           printf(" %f", wPy);          //12
1576           printf(" %f", wPz);          //13
1577           printf(" %f", wSigmaC2*1000000);  //14
1578           printf(" %f", wSigmaTgl2*1000);   //15
1579           printf(" %f", wSigmaY2);     //16
1580           //      printf(" %f", wSigmaZ2);     //17
1581           printf(" %f", wChi2);     //17
1582           printf(" %f", wC);           //18
1583           printf("\n");
1584         } 
1585         */                        
1586       }
1587     }  
1588   }
1589   if (nr<outerTB) 
1590     t.SetStop(kTRUE);
1591   else
1592     t.SetStop(kFALSE);
1593   return expectedNumberOfClusters;
1594
1595
1596 }         
1597
1598 //---------------------------------------------------------------------------
1599 Int_t AliTRDtracker::Refit(AliTRDtrack& t, Int_t rf)
1600 {
1601   // Starting from current position on track=t this function tries
1602   // to extrapolate the track up to timeBin=0 and to reuse already
1603   // assigned clusters. Returns the number of clusters
1604   // expected to be found in sensitive layers
1605   // get indices of assigned clusters for each layer
1606   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
1607
1608   Int_t iCluster[90];
1609   for (Int_t i = 0; i < 90; i++) iCluster[i] = 0;
1610   for (Int_t i = 0; i < t.GetNumberOfClusters(); i++) {
1611     Int_t index = t.GetClusterIndex(i);
1612     AliTRDcluster *cl=(AliTRDcluster*) GetCluster(index);
1613     if (!cl) continue;
1614     Int_t detector=cl->GetDetector();
1615     Int_t localTimeBin=cl->GetLocalTimeBin();
1616     Int_t sector=fGeom->GetSector(detector);
1617     Int_t plane=fGeom->GetPlane(detector);
1618
1619     Int_t trackingSector = CookSectorIndex(sector);
1620
1621     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1622     if(gtb < 0) continue; 
1623     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1624     iCluster[layer] = index;
1625   }
1626   t.ResetClusters();
1627
1628   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
1629
1630   Double_t alpha=t.GetAlpha();
1631   alpha = TVector2::Phi_0_2pi(alpha);
1632
1633   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
1634   Double_t radLength, rho, x, dx, y, ymax, z;
1635
1636   Int_t expectedNumberOfClusters = 0;
1637   Bool_t lookForCluster;
1638
1639   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1640
1641  
1642   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) { 
1643
1644     y = t.GetY(); z = t.GetZ();
1645
1646     // first propagate to the inner surface of the current time bin 
1647     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1648     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1649     if(!t.PropagateTo(x,radLength,rho)) break;
1650     y = t.GetY();
1651     ymax = x*TMath::Tan(0.5*alpha);
1652     if (y > ymax) {
1653       s = (s+1) % ns;
1654       if (!t.Rotate(alpha)) break;
1655       if(!t.PropagateTo(x,radLength,rho)) break;
1656     } else if (y <-ymax) {
1657       s = (s-1+ns) % ns;                           
1658       if (!t.Rotate(-alpha)) break;   
1659       if(!t.PropagateTo(x,radLength,rho)) break;
1660     } 
1661
1662     y = t.GetY(); z = t.GetZ();
1663
1664     // now propagate to the middle plane of the next time bin 
1665     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1666     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1667     if(!t.PropagateTo(x,radLength,rho)) break;
1668     y = t.GetY();
1669     ymax = x*TMath::Tan(0.5*alpha);
1670     if (y > ymax) {
1671       s = (s+1) % ns;
1672       if (!t.Rotate(alpha)) break;
1673       if(!t.PropagateTo(x,radLength,rho)) break;
1674     } else if (y <-ymax) {
1675       s = (s-1+ns) % ns;                           
1676       if (!t.Rotate(-alpha)) break;   
1677       if(!t.PropagateTo(x,radLength,rho)) break;
1678     } 
1679
1680     if(lookForCluster) expectedNumberOfClusters++;       
1681
1682     // use assigned cluster
1683     if (!iCluster[nr-1]) continue;
1684     AliTRDcluster *cl=(AliTRDcluster*)GetCluster(iCluster[nr-1]);
1685     Double_t h01 = GetTiltFactor(cl);
1686     Double_t chi2=t.GetPredictedChi2(cl, h01);
1687     t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); 
1688     t.Update(cl,chi2,iCluster[nr-1],h01);
1689   }
1690
1691   return expectedNumberOfClusters;
1692 }                
1693
1694 //___________________________________________________________________
1695
1696 Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1697 {
1698   // Starting from current radial position of track <t> this function
1699   // extrapolates the track up to radial position <xToGo>. 
1700   // Returns 1 if track reaches the plane, and 0 otherwise 
1701
1702   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
1703
1704   Double_t alpha=t.GetAlpha();
1705
1706   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1707   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1708
1709   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
1710
1711   Bool_t lookForCluster;
1712   Double_t radLength, rho, x, dx, y, ymax, z;
1713
1714   x = t.GetX();
1715
1716   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1717
1718   Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1719
1720   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) { 
1721
1722     y = t.GetY(); z = t.GetZ();
1723
1724     // first propagate to the outer surface of the current time bin 
1725     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1726     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
1727     if(!t.PropagateTo(x,radLength,rho)) return 0;
1728     y = t.GetY();
1729     ymax = x*TMath::Tan(0.5*alpha);
1730     if (y > ymax) {
1731       s = (s+1) % ns;
1732       if (!t.Rotate(alpha)) return 0;
1733     } else if (y <-ymax) {
1734       s = (s-1+ns) % ns;                           
1735       if (!t.Rotate(-alpha)) return 0;   
1736     } 
1737     if(!t.PropagateTo(x,radLength,rho)) return 0;
1738
1739     y = t.GetY(); z = t.GetZ();
1740
1741     // now propagate to the middle plane of the next time bin 
1742     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1743     x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
1744     if(!t.PropagateTo(x,radLength,rho)) return 0;
1745     y = t.GetY();
1746     ymax = x*TMath::Tan(0.5*alpha);
1747     if (y > ymax) {
1748       s = (s+1) % ns;
1749       if (!t.Rotate(alpha)) return 0;
1750     } else if (y <-ymax) {
1751       s = (s-1+ns) % ns;                           
1752       if (!t.Rotate(-alpha)) return 0;   
1753     } 
1754     if(!t.PropagateTo(x,radLength,rho)) return 0;
1755   }
1756   return 1;
1757 }         
1758
1759 //___________________________________________________________________
1760
1761 Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1762 {
1763   // Starting from current radial position of track <t> this function
1764   // extrapolates the track up to radial position of the outermost
1765   // padrow of the TPC. 
1766   // Returns 1 if track reaches the TPC, and 0 otherwise 
1767
1768   //Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
1769
1770   Double_t alpha=t.GetAlpha();
1771   alpha = TVector2::Phi_0_2pi(alpha);
1772
1773   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
1774
1775   Bool_t lookForCluster;
1776   Double_t radLength, rho, x, dx, y, /*ymax,*/ z;
1777
1778   x = t.GetX();
1779
1780   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1781   Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1782
1783   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) { 
1784
1785     y = t.GetY(); 
1786     z = t.GetZ();
1787
1788     // first propagate to the outer surface of the current time bin 
1789     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1790     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; 
1791     
1792     if(!t.PropagateTo(x,radLength,rho)) return 0;
1793     AdjustSector(&t);
1794     if(!t.PropagateTo(x,radLength,rho)) return 0;
1795
1796     y = t.GetY(); 
1797     z = t.GetZ();
1798
1799     // now propagate to the middle plane of the next time bin 
1800     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1801     x = fTrSec[s]->GetLayer(nr-1)->GetX(); 
1802     
1803     if(!t.PropagateTo(x,radLength,rho)) return 0;
1804     AdjustSector(&t);
1805     if(!t.PropagateTo(x,radLength,rho)) return 0;
1806   }
1807   return 1;
1808 }         
1809
1810 void AliTRDtracker::LoadEvent()
1811 {
1812   // Fills clusters into TRD tracking_sectors 
1813   // Note that the numbering scheme for the TRD tracking_sectors 
1814   // differs from that of TRD sectors
1815
1816   ReadClusters(fClusters);
1817   Int_t ncl=fClusters->GetEntriesFast();
1818   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1819               
1820   UInt_t index;
1821   while (ncl--) {
1822 //    printf("\r %d left  ",ncl); 
1823     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1824     Int_t detector=c->GetDetector();
1825     Int_t localTimeBin=c->GetLocalTimeBin();
1826     Int_t sector=fGeom->GetSector(detector);
1827     Int_t plane=fGeom->GetPlane(detector);
1828
1829     Int_t trackingSector = CookSectorIndex(sector);
1830
1831     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1832     if(gtb < 0) continue; 
1833     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1834
1835     index=ncl;
1836     fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
1837   }    
1838   printf("\r\n");
1839
1840 }
1841
1842 //_____________________________________________________________________________
1843 Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1844 {
1845   // Fills clusters into TRD tracking_sectors 
1846   // Note that the numbering scheme for the TRD tracking_sectors 
1847   // differs from that of TRD sectors
1848
1849   if (ReadClusters(fClusters,cTree)) {
1850      Error("LoadClusters","Problem with reading the clusters !");
1851      return 1;
1852   }
1853   Int_t ncl=fClusters->GetEntriesFast();
1854   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1855               
1856   UInt_t index;
1857   for (Int_t ichamber=0;ichamber<5;ichamber++)
1858     for (Int_t isector=0;isector<18;isector++){
1859       fHoles[ichamber][isector]=kTRUE;
1860     }
1861
1862
1863   while (ncl--) {
1864 //    printf("\r %d left  ",ncl); 
1865     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1866     Int_t detector=c->GetDetector();
1867     Int_t localTimeBin=c->GetLocalTimeBin();
1868     Int_t sector=fGeom->GetSector(detector);
1869     Int_t plane=fGeom->GetPlane(detector);
1870       
1871     Int_t trackingSector = CookSectorIndex(sector);
1872     if (c->GetLabel(0)>0){
1873       Int_t chamber = fGeom->GetChamber(detector);
1874       fHoles[chamber][trackingSector]=kFALSE;
1875     }
1876
1877     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1878     if(gtb < 0) continue; 
1879     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1880
1881     index=ncl;
1882     fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
1883   }    
1884   printf("\r\n");
1885   //
1886   //
1887   /*
1888   for (Int_t isector=0;isector<18;isector++){
1889     for (Int_t ichamber=0;ichamber<5;ichamber++)      
1890       if (fHoles[ichamber][isector]!=fGeom->IsHole(0,ichamber,17-isector)) 
1891         printf("Problem \t%d\t%d\t%d\t%d\n",isector,ichamber,fHoles[ichamber][isector],
1892              fGeom->IsHole(0,ichamber,17-isector));
1893   }
1894   */
1895   return 0;
1896 }
1897
1898 //_____________________________________________________________________________
1899 void AliTRDtracker::UnloadEvent() 
1900
1901   //
1902   // Clears the arrays of clusters and tracks. Resets sectors and timebins 
1903   //
1904
1905   Int_t i, nentr;
1906
1907   nentr = fClusters->GetEntriesFast();
1908   for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1909
1910   nentr = fSeeds->GetEntriesFast();
1911   for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1912
1913   nentr = fTracks->GetEntriesFast();
1914   for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1915
1916   Int_t nsec = AliTRDgeometry::kNsect;
1917
1918   for (i = 0; i < nsec; i++) {    
1919     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1920       fTrSec[i]->GetLayer(pl)->Clear();
1921     }
1922   }
1923
1924 }
1925
1926 //__________________________________________________________________________
1927 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1928 {
1929   // Creates track seeds using clusters in timeBins=i1,i2
1930
1931   if(turn > 2) {
1932     cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
1933     return;
1934   }
1935
1936   Double_t x[5], c[15];
1937   Int_t maxSec=AliTRDgeometry::kNsect;
1938   
1939   Double_t alpha=AliTRDgeometry::GetAlpha();
1940   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1941   Double_t cs=cos(alpha), sn=sin(alpha);
1942   Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1943     
1944       
1945   Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1946   Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1947       
1948   Double_t x1 =fTrSec[0]->GetX(i1);
1949   Double_t xx2=fTrSec[0]->GetX(i2);
1950       
1951   for (Int_t ns=0; ns<maxSec; ns++) {
1952     
1953     Int_t nl2 = *(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
1954     Int_t nl=(*fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
1955     Int_t nm=(*fTrSec[ns]->GetLayer(i2));
1956     Int_t nu=(*fTrSec[(ns+1)%maxSec]->GetLayer(i2));
1957     Int_t nu2=(*fTrSec[(ns+2)%maxSec]->GetLayer(i2));
1958     
1959     AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1960     
1961     for (Int_t is=0; is < r1; is++) {
1962       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
1963       
1964       for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
1965         
1966         const AliTRDcluster *cl;
1967         Double_t x2,   y2,   z2;
1968         Double_t x3=0., y3=0.;   
1969         
1970         if (js<nl2) {
1971           if(turn != 2) continue;
1972           AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
1973           cl=r2[js];
1974           y2=cl->GetY(); z2=cl->GetZ();
1975           
1976           x2= xx2*cs2+y2*sn2;
1977           y2=-xx2*sn2+y2*cs2;
1978         }
1979         else if (js<nl2+nl) {
1980           if(turn != 1) continue;
1981           AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
1982           cl=r2[js-nl2];
1983           y2=cl->GetY(); z2=cl->GetZ();
1984           
1985           x2= xx2*cs+y2*sn;
1986           y2=-xx2*sn+y2*cs;
1987         }                                
1988         else if (js<nl2+nl+nm) {
1989           if(turn != 1) continue;
1990           AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1991           cl=r2[js-nl2-nl];
1992           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1993         }
1994         else if (js<nl2+nl+nm+nu) {
1995           if(turn != 1) continue;
1996           AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%maxSec]->GetLayer(i2));
1997           cl=r2[js-nl2-nl-nm];
1998           y2=cl->GetY(); z2=cl->GetZ();
1999           
2000           x2=xx2*cs-y2*sn;
2001           y2=xx2*sn+y2*cs;
2002         }              
2003         else {
2004           if(turn != 2) continue;
2005           AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%maxSec]->GetLayer(i2));
2006           cl=r2[js-nl2-nl-nm-nu];
2007           y2=cl->GetY(); z2=cl->GetZ();
2008           
2009           x2=xx2*cs2-y2*sn2;
2010           y2=xx2*sn2+y2*cs2;
2011         }
2012         
2013         if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue;
2014         
2015         Double_t zz=z1 - z1/x1*(x1-x2);
2016         
2017         if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue;
2018         
2019         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
2020         if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
2021         
2022         x[0]=y1;
2023         x[1]=z1;
2024         x[4]=f1trd(x1,y1,x2,y2,x3,y3);
2025         
2026         if (TMath::Abs(x[4]) > fgkMaxSeedC) continue;      
2027         
2028         x[2]=f2trd(x1,y1,x2,y2,x3,y3);
2029         
2030         if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
2031         
2032         x[3]=f3trd(x1,y1,x2,y2,z1,z2);
2033         
2034         if (TMath::Abs(x[3]) > fgkMaxSeedTan) continue;
2035         
2036         Double_t a=asin(x[2]);
2037         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
2038         
2039         if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue;
2040         
2041         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
2042         Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2();
2043         Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ;  
2044
2045         // Tilt changes
2046         Double_t h01 = GetTiltFactor(r1[is]);
2047         Double_t xuFactor = 100.;
2048         if(fNoTilt) { 
2049           h01 = 0;
2050           xuFactor = 1;
2051         }
2052
2053         sy1=sy1+sz1*h01*h01;
2054         Double_t syz=sz1*(-h01);
2055         // end of tilt changes
2056         
2057         Double_t f40=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2058         Double_t f42=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2059         Double_t f43=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2060         Double_t f20=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2061         Double_t f22=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2062         Double_t f23=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2063         Double_t f30=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2064         Double_t f31=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2065         Double_t f32=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2066         Double_t f34=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;    
2067
2068         
2069         c[0]=sy1;
2070         //        c[1]=0.;       c[2]=sz1;
2071         c[1]=syz;       c[2]=sz1*xuFactor;
2072         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2073         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
2074                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2075         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2076         c[13]=f30*sy1*f40+f32*sy2*f42;
2077         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;      
2078         
2079         UInt_t index=r1.GetIndex(is);
2080         
2081         AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
2082
2083         Int_t rc=FollowProlongation(*track, i2);     
2084         
2085         if ((rc < 1) ||
2086             (track->GetNumberOfClusters() < 
2087              (outer-inner)*fgkMinClustersInSeed)) delete track;
2088         else {
2089           fSeeds->AddLast(track); fNseeds++;
2090 //          cerr<<"\r found seed "<<fNseeds;
2091         }
2092       }
2093     }
2094   }
2095 }            
2096
2097 //_____________________________________________________________________________
2098 Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) 
2099 {
2100   //
2101   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) 
2102   // from the file. The names of the cluster tree and branches 
2103   // should match the ones used in AliTRDclusterizer::WriteClusters()
2104   //
2105   TObjArray *clusterArray = new TObjArray(400); 
2106   
2107   TBranch *branch=ClusterTree->GetBranch("TRDcluster");
2108   if (!branch) {
2109     Error("ReadClusters","Can't get the branch !");
2110     return 1;
2111   }
2112   branch->SetAddress(&clusterArray); 
2113   
2114   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
2115   printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
2116   
2117   // Loop through all entries in the tree
2118   Int_t nbytes;
2119   AliTRDcluster *c = 0;
2120   printf("\n");
2121
2122   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
2123     
2124     // Import the tree
2125     nbytes += ClusterTree->GetEvent(iEntry);  
2126     
2127     // Get the number of points in the detector
2128     Int_t nCluster = clusterArray->GetEntriesFast();  
2129 //    printf("\r Read %d clusters from entry %d", nCluster, iEntry);
2130     
2131     // Loop through all TRD digits
2132     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
2133       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2134       AliTRDcluster *co = new AliTRDcluster(*c);
2135       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
2136       Int_t ltb = co->GetLocalTimeBin();
2137       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
2138       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
2139       array->AddLast(co);
2140       delete clusterArray->RemoveAt(iCluster); 
2141     }
2142   }
2143
2144   delete clusterArray;
2145
2146   return 0;
2147 }
2148
2149 //______________________________________________________________________
2150 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
2151 {
2152   //
2153   // Reads AliTRDclusters from file <filename>. The names of the cluster
2154   // tree and branches should match the ones used in
2155   // AliTRDclusterizer::WriteClusters()
2156   // if <array> == 0, clusters are added into AliTRDtracker fCluster array
2157   //
2158
2159   TDirectory *savedir=gDirectory;
2160
2161   TFile *file = TFile::Open(filename);
2162   if (!file->IsOpen()) {
2163     cerr<<"Can't open file with TRD clusters"<<endl;
2164     return;
2165   }
2166
2167   Char_t treeName[12];
2168   sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
2169   TTree *clusterTree = (TTree*) gDirectory->Get(treeName);
2170
2171   if (!clusterTree) {
2172      cerr<<"AliTRDtracker::ReadClusters(): ";
2173      cerr<<"can't get a tree with clusters !\n";
2174      return;
2175   }
2176
2177   TObjArray *clusterArray = new TObjArray(400);
2178
2179   clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
2180
2181   Int_t nEntries = (Int_t) clusterTree->GetEntries();
2182   cout<<"found "<<nEntries<<" in clusterTree"<<endl;   
2183
2184   // Loop through all entries in the tree
2185   Int_t nbytes;
2186   AliTRDcluster *c = 0;
2187
2188   printf("\n");
2189
2190   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2191
2192     // Import the tree
2193     nbytes += clusterTree->GetEvent(iEntry);
2194
2195     // Get the number of points in the detector
2196     Int_t nCluster = clusterArray->GetEntriesFast();
2197     printf("\n Read %d clusters from entry %d", nCluster, iEntry);
2198
2199     // Loop through all TRD digits
2200     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2201       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2202       AliTRDcluster *co = new AliTRDcluster(*c);
2203       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
2204       Int_t ltb = co->GetLocalTimeBin();
2205       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
2206       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
2207       array->AddLast(co);
2208       delete clusterArray->RemoveAt(iCluster);
2209     }
2210   }
2211
2212   file->Close();
2213   delete clusterArray;
2214   savedir->cd();
2215
2216 }                      
2217
2218 void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp) 
2219 {
2220   //
2221   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) 
2222   // from the file. The names of the cluster tree and branches 
2223   // should match the ones used in AliTRDclusterizer::WriteClusters()
2224   //
2225
2226   TDirectory *savedir=gDirectory; 
2227
2228   if (inp) {
2229      TFile *in=(TFile*)inp;
2230      if (!in->IsOpen()) {
2231         cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n";
2232         return;
2233      }
2234      else{
2235        in->cd();
2236      }
2237   }
2238
2239   Char_t treeName[12];
2240   sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
2241   TTree *clusterTree = (TTree*) gDirectory->Get(treeName);
2242   
2243   TObjArray *clusterArray = new TObjArray(400); 
2244   
2245   clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray); 
2246   
2247   Int_t nEntries = (Int_t) clusterTree->GetEntries();
2248   printf("found %d entries in %s.\n",nEntries,clusterTree->GetName());
2249   
2250   // Loop through all entries in the tree
2251   Int_t nbytes;
2252   AliTRDcluster *c = 0;
2253   printf("\n");
2254
2255   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
2256     
2257     // Import the tree
2258     nbytes += clusterTree->GetEvent(iEntry);  
2259     
2260     // Get the number of points in the detector
2261     Int_t nCluster = clusterArray->GetEntriesFast();  
2262 //    printf("\r Read %d clusters from entry %d", nCluster, iEntry);
2263     
2264     // Loop through all TRD digits
2265     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
2266       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2267       AliTRDcluster *co = new AliTRDcluster(*c);
2268       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
2269       Int_t ltb = co->GetLocalTimeBin();
2270       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
2271       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
2272       array->AddLast(co);
2273       delete clusterArray->RemoveAt(iCluster); 
2274     }
2275   }
2276
2277   delete clusterArray;
2278   savedir->cd();   
2279
2280 }
2281
2282 //__________________________________________________________________
2283 void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const 
2284 {
2285   //
2286   // This cooks a label. Mmmmh, smells good...
2287   //
2288
2289   Int_t label=123456789, index, i, j;
2290   Int_t ncl=pt->GetNumberOfClusters();
2291   const Int_t kRange = fTrSec[0]->GetOuterTimeBin()+1;
2292
2293   Bool_t labelAdded;
2294
2295   //  Int_t s[kRange][2];
2296   Int_t **s = new Int_t* [kRange];
2297   for (i=0; i<kRange; i++) {
2298     s[i] = new Int_t[2];
2299   }
2300   for (i=0; i<kRange; i++) {
2301     s[i][0]=-1;
2302     s[i][1]=0;
2303   }
2304
2305   Int_t t0,t1,t2;
2306   for (i=0; i<ncl; i++) {
2307     index=pt->GetClusterIndex(i);
2308     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2309     t0=c->GetLabel(0);
2310     t1=c->GetLabel(1);
2311     t2=c->GetLabel(2);
2312   }
2313
2314   for (i=0; i<ncl; i++) {
2315     index=pt->GetClusterIndex(i);
2316     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2317     for (Int_t k=0; k<3; k++) { 
2318       label=c->GetLabel(k);
2319       labelAdded=kFALSE; j=0;
2320       if (label >= 0) {
2321         while ( (!labelAdded) && ( j < kRange ) ) {
2322           if (s[j][0]==label || s[j][1]==0) {
2323             s[j][0]=label; 
2324             s[j][1]=s[j][1]+1; 
2325             labelAdded=kTRUE;
2326           }
2327           j++;
2328         }
2329       }
2330     }
2331   }
2332
2333   Int_t max=0;
2334   label = -123456789;
2335
2336   for (i=0; i<kRange; i++) {
2337     if (s[i][1]>max) {
2338       max=s[i][1]; label=s[i][0];
2339     }
2340   }
2341
2342   for (i=0; i<kRange; i++) {
2343     delete []s[i];
2344   }        
2345
2346   delete []s;
2347
2348   if ((1.- Float_t(max)/ncl) > wrong) label=-label;   
2349
2350   pt->SetLabel(label); 
2351
2352 }
2353
2354
2355 //__________________________________________________________________
2356 void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const 
2357 {
2358   //
2359   // Use clusters, but don't abuse them!
2360   //
2361
2362   Int_t ncl=t->GetNumberOfClusters();
2363   for (Int_t i=from; i<ncl; i++) {
2364     Int_t index = t->GetClusterIndex(i);
2365     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2366     c->Use();
2367   }
2368 }
2369
2370
2371 //_____________________________________________________________________
2372 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
2373 {
2374   // Parametrised "expected" error of the cluster reconstruction in Y 
2375
2376   Double_t s = 0.08 * 0.08;    
2377   return s;
2378 }
2379
2380 //_____________________________________________________________________
2381 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
2382 {
2383   // Parametrised "expected" error of the cluster reconstruction in Z 
2384
2385   Double_t s = 9 * 9 /12.;  
2386   return s;
2387 }                  
2388
2389 //_____________________________________________________________________
2390 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const 
2391 {
2392   //
2393   // Returns radial position which corresponds to time bin <localTB>
2394   // in tracking sector <sector> and plane <plane>
2395   //
2396
2397   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB); 
2398   Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2399   return fTrSec[sector]->GetLayer(pl)->GetX();
2400
2401 }
2402
2403
2404 //_______________________________________________________
2405 AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x, 
2406                Double_t dx, Double_t rho, Double_t radLength, Int_t tbIndex)
2407
2408   //
2409   // AliTRDpropagationLayer constructor
2410   //
2411
2412   fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = radLength;
2413   fClusters = NULL; fIndex = NULL; fTimeBinIndex = tbIndex;
2414
2415
2416   for(Int_t i=0; i < (Int_t) kZones; i++) {
2417     fZc[i]=0; fZmax[i] = 0;
2418   }
2419
2420   fYmax = 0;
2421
2422   if(fTimeBinIndex >= 0) { 
2423     fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
2424     fIndex = new UInt_t[kMaxClusterPerTimeBin];
2425   }
2426
2427   for (Int_t i=0;i<5;i++) fIsHole[i] = kFALSE;
2428   fHole = kFALSE;
2429   fHoleZc = 0;
2430   fHoleZmax = 0;
2431   fHoleYc = 0;
2432   fHoleYmax = 0;
2433   fHoleRho = 0;
2434   fHoleX0 = 0;
2435
2436 }
2437
2438 //_______________________________________________________
2439 void AliTRDtracker::AliTRDpropagationLayer::SetHole(
2440           Double_t Zmax, Double_t Ymax, Double_t rho, 
2441           Double_t radLength, Double_t Yc, Double_t Zc) 
2442 {
2443   //
2444   // Sets hole in the layer 
2445   //
2446   fHole = kTRUE;
2447   fHoleZc = Zc;
2448   fHoleZmax = Zmax;
2449   fHoleYc = Yc;
2450   fHoleYmax = Ymax;
2451   fHoleRho = rho;
2452   fHoleX0 = radLength;
2453 }
2454   
2455
2456 //_______________________________________________________
2457 AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
2458 {
2459   //
2460   // AliTRDtrackingSector Constructor
2461   //
2462
2463   fGeom = geo;
2464   fPar = par;
2465   fGeomSector = gs;
2466   fTzeroShift = 0.13;
2467   fN = 0;
2468   //
2469   // get holes description from geometry
2470   Bool_t holes[AliTRDgeometry::kNcham];
2471   //printf("sector\t%d\t",gs);
2472   for (Int_t icham=0; icham<AliTRDgeometry::kNcham;icham++){
2473     holes[icham] = fGeom->IsHole(0,icham,gs);
2474     //printf("%d",holes[icham]);
2475   } 
2476   //printf("\n");
2477   
2478   for(UInt_t i=0; i < kMaxTimeBinIndex; i++) fTimeBinIndex[i] = -1;
2479
2480
2481   AliTRDpropagationLayer* ppl;
2482
2483   Double_t x, xin, xout, dx, rho, radLength;
2484   Int_t    steps;
2485
2486   // set time bins in the gas of the TPC
2487
2488   xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
2489   rho = 0.9e-3;  radLength = 28.94;
2490
2491   for(Int_t i=0; i<steps; i++) {
2492     x = xin + i*dx + dx/2;
2493     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2494     InsertLayer(ppl);
2495   }
2496
2497   // set time bins in the outer field cage vessel
2498
2499   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2500   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2501   InsertLayer(ppl);
2502
2503   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2504   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2505   InsertLayer(ppl);
2506
2507   dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
2508   steps = 5; dx = (xout - xin)/steps;
2509   for(Int_t i=0; i<steps; i++) {
2510     x = xin + i*dx + dx/2;
2511     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2512     InsertLayer(ppl);
2513   }
2514
2515   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2516   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2517   InsertLayer(ppl);
2518
2519   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2520   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2521   InsertLayer(ppl);
2522
2523
2524   // set time bins in CO2
2525
2526   xin = xout; xout = 275.0; 
2527   steps = 50; dx = (xout - xin)/steps;
2528   rho = 1.977e-3;  radLength = 36.2;
2529   
2530   for(Int_t i=0; i<steps; i++) {
2531     x = xin + i*dx + dx/2;
2532     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2533     InsertLayer(ppl);
2534   }
2535
2536   // set time bins in the outer containment vessel
2537
2538   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2539   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2540   InsertLayer(ppl);
2541
2542   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2543   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2544   InsertLayer(ppl);
2545
2546   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2547   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2548   InsertLayer(ppl);
2549
2550   dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
2551   steps = 10; dx = (xout - xin)/steps;
2552   for(Int_t i=0; i<steps; i++) {
2553     x = xin + i*dx + dx/2;
2554     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2555     InsertLayer(ppl);
2556   }
2557
2558   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
2559   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2560   InsertLayer(ppl);
2561
2562   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
2563   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2564   InsertLayer(ppl);
2565   
2566   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
2567   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2568   InsertLayer(ppl);
2569
2570   Double_t xtrd = (Double_t) fGeom->Rmin();  
2571
2572   // add layers between TPC and TRD (Air temporarily)
2573   xin = xout; xout = xtrd;
2574   steps = 50; dx = (xout - xin)/steps;
2575   rho = 1.2e-3;  radLength = 36.66;
2576   
2577   for(Int_t i=0; i<steps; i++) {
2578     x = xin + i*dx + dx/2;
2579     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2580     InsertLayer(ppl);
2581   }
2582
2583
2584   //  Double_t alpha=AliTRDgeometry::GetAlpha();
2585
2586   // add layers for each of the planes
2587
2588   Double_t dxRo = (Double_t) fGeom->CroHght();    // Rohacell 
2589   Double_t dxSpace = (Double_t) fGeom->Cspace();  // Spacing between planes
2590   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2591   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region  
2592   Double_t dxRad = (Double_t) fGeom->CraHght();   // Radiator
2593   Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo; 
2594   Double_t dxPlane = dxTEC + dxSpace; 
2595
2596   Int_t tb, tbIndex;
2597   const Int_t  kNchambers = AliTRDgeometry::Ncham();
2598   Double_t  ymax = 0;
2599   //, holeYmax = 0;
2600   Double_t ymaxsensitive=0;
2601   Double_t *zc = new Double_t[kNchambers];
2602   Double_t *zmax = new Double_t[kNchambers];
2603   Double_t *zmaxsensitive = new Double_t[kNchambers];  
2604   //  Double_t  holeZmax = 1000.;   // the whole sector is missing
2605
2606   for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2607     //
2608     // Radiator 
2609     xin = xtrd + plane * dxPlane; xout = xin + dxRad;
2610     steps = 12; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6; 
2611     for(Int_t i=0; i<steps; i++) {
2612       x = xin + i*dx + dx/2;
2613       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);      
2614       InsertLayer(ppl);
2615     }
2616
2617     ymax          = fGeom->GetChamberWidth(plane)/2.;
2618     ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
2619     
2620     for(Int_t ch = 0; ch < kNchambers; ch++) {
2621       zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
2622       Float_t pad = fPar->GetRowPadSize(plane,ch,0);
2623       Float_t row0 = fPar->GetRow0(plane,ch,0);
2624       Int_t nPads = fPar->GetRowMax(plane,ch,0);
2625       zmaxsensitive[ch] = Float_t(nPads)*pad/2.;      
2626       //      zc[ch] = (pad * nPads)/2 + row0 - pad/2;
2627       zc[ch] = (pad * nPads)/2 + row0;
2628       //zc[ch] = row0+zmax[ch]-AliTRDgeometry::RpadW();
2629
2630     }
2631
2632     dx = fPar->GetTimeBinSize(); 
2633     rho = 0.00295 * 0.85; radLength = 11.0;  
2634
2635     Double_t x0 = (Double_t) fPar->GetTime0(plane);
2636     Double_t xbottom = x0 - dxDrift;
2637     Double_t xtop = x0 + dxAmp;
2638     //
2639     // Amplification region
2640     steps = (Int_t) (dxAmp/dx);
2641
2642     for(tb = 0; tb < steps; tb++) {
2643       x = x0 + tb * dx + dx/2;
2644       tbIndex = CookTimeBinIndex(plane, -tb-1);
2645       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2646       ppl->SetYmax(ymax,ymaxsensitive);
2647       ppl->SetZ(zc, zmax, zmaxsensitive);
2648       ppl->SetHoles(holes);
2649       InsertLayer(ppl);
2650     }
2651     tbIndex = CookTimeBinIndex(plane, -steps);
2652     x = (x + dx/2 + xtop)/2;
2653     dx = 2*(xtop-x);
2654     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2655     ppl->SetYmax(ymax,ymaxsensitive);
2656     ppl->SetZ(zc, zmax,zmaxsensitive);
2657     ppl->SetHoles(holes);
2658     InsertLayer(ppl);
2659
2660     // Drift region
2661     dx = fPar->GetTimeBinSize();
2662     steps = (Int_t) (dxDrift/dx);
2663
2664     for(tb = 0; tb < steps; tb++) {
2665       x = x0 - tb * dx - dx/2;
2666       tbIndex = CookTimeBinIndex(plane, tb);
2667
2668       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2669       ppl->SetYmax(ymax,ymaxsensitive);
2670       ppl->SetZ(zc, zmax, zmaxsensitive);
2671       ppl->SetHoles(holes);
2672       InsertLayer(ppl);
2673     }
2674     tbIndex = CookTimeBinIndex(plane, steps);
2675     x = (x - dx/2 + xbottom)/2;
2676     dx = 2*(x-xbottom);
2677     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
2678     ppl->SetYmax(ymax,ymaxsensitive);
2679     ppl->SetZ(zc, zmax, zmaxsensitive);
2680     ppl->SetHoles(holes);    
2681     InsertLayer(ppl);
2682
2683     // Pad Plane
2684     xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; radLength = 33.0;
2685     ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
2686     ppl->SetYmax(ymax,ymaxsensitive);
2687     ppl->SetZ(zc, zmax,zmax);
2688     ppl->SetHoles(holes);         
2689     InsertLayer(ppl);
2690
2691     // Rohacell
2692     xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
2693     steps = 5; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6; 
2694     for(Int_t i=0; i<steps; i++) {
2695       x = xin + i*dx + dx/2;
2696       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2697       ppl->SetYmax(ymax,ymaxsensitive);
2698       ppl->SetZ(zc, zmax,zmax);
2699       ppl->SetHoles(holes);
2700       InsertLayer(ppl);
2701     }
2702
2703     // Space between the chambers, air
2704     xin = xout; xout = xtrd + (plane + 1) * dxPlane;
2705     steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66; 
2706     for(Int_t i=0; i<steps; i++) {
2707       x = xin + i*dx + dx/2;
2708       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2709       InsertLayer(ppl);
2710     }
2711   }    
2712
2713   // Space between the TRD and RICH
2714   Double_t xRICH = 500.;
2715   xin = xout; xout = xRICH;
2716   steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66; 
2717   for(Int_t i=0; i<steps; i++) {
2718     x = xin + i*dx + dx/2;
2719     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
2720     InsertLayer(ppl);
2721   }
2722
2723   MapTimeBinLayers();
2724   delete [] zc;
2725   delete [] zmax;
2726
2727 }
2728
2729 //______________________________________________________
2730
2731 Int_t  AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t localTB) const
2732 {
2733   //
2734   // depending on the digitization parameters calculates "global"
2735   // time bin index for timebin <localTB> in plane <plane>
2736   //
2737
2738   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2739   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region  
2740   Double_t dx = (Double_t) fPar->GetTimeBinSize();  
2741
2742   Int_t tbAmp = fPar->GetTimeBefore();
2743   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
2744   if(kTRUE) maxAmp = 0;   // intentional until we change parameter class 
2745   Int_t tbDrift = fPar->GetTimeMax();
2746   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
2747
2748   Int_t tbPerPlane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
2749
2750   Int_t gtb = (plane+1) * tbPerPlane - localTB - 1 - TMath::Min(tbAmp,maxAmp);
2751
2752   if((localTB < 0) && 
2753      (TMath::Abs(localTB) > TMath::Min(tbAmp,maxAmp))) return -1;
2754   if(localTB >= TMath::Min(tbDrift,maxDrift)) return -1;
2755
2756   return gtb;
2757
2758
2759 }
2760
2761 //______________________________________________________
2762
2763 void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers() 
2764 {
2765   //
2766   // For all sensitive time bins sets corresponding layer index
2767   // in the array fTimeBins 
2768   //
2769
2770   Int_t index;
2771
2772   for(Int_t i = 0; i < fN; i++) {
2773     index = fLayers[i]->GetTimeBinIndex();
2774     
2775     //    printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2776
2777     if(index < 0) continue;
2778     if(index >= (Int_t) kMaxTimeBinIndex) {
2779       printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2780       printf("    index %d exceeds allowed maximum of %d!\n",
2781              index, kMaxTimeBinIndex-1);
2782       continue;
2783     }
2784     fTimeBinIndex[index] = i;
2785   }
2786
2787   Double_t x1, dx1, x2, dx2, gap;
2788
2789   for(Int_t i = 0; i < fN-1; i++) {
2790     x1 = fLayers[i]->GetX();
2791     dx1 = fLayers[i]->GetdX();
2792     x2 = fLayers[i+1]->GetX();
2793     dx2 = fLayers[i+1]->GetdX();
2794     gap = (x2 - dx2/2) - (x1 + dx1/2);
2795     if(gap < -0.01) {
2796       printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2797       printf("             %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2798     }
2799     if(gap > 0.01) { 
2800       printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2801       printf("             (%f - %f) - (%f + %f) = %f\n", 
2802              x2, dx2/2, x1, dx1, gap);
2803     }
2804   }
2805 }
2806   
2807
2808 //______________________________________________________
2809
2810
2811 Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2812 {
2813   // 
2814   // Returns the number of time bin which in radial position is closest to <x>
2815   //
2816
2817   if(x >= fLayers[fN-1]->GetX()) return fN-1; 
2818   if(x <= fLayers[0]->GetX()) return 0; 
2819
2820   Int_t b=0, e=fN-1, m=(b+e)/2;
2821   for (; b<e; m=(b+e)/2) {
2822     if (x > fLayers[m]->GetX()) b=m+1;
2823     else e=m;
2824   }
2825   if(TMath::Abs(x - fLayers[m]->GetX()) > 
2826      TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2827   else return m;
2828
2829 }
2830
2831 //______________________________________________________
2832
2833 Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const 
2834 {
2835   // 
2836   // Returns number of the innermost SENSITIVE propagation layer
2837   //
2838
2839   return GetLayerNumber(0);
2840 }
2841
2842 //______________________________________________________
2843
2844 Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const 
2845 {
2846   // 
2847   // Returns number of the outermost SENSITIVE time bin
2848   //
2849
2850   return GetLayerNumber(GetNumberOfTimeBins() - 1);
2851 }
2852
2853 //______________________________________________________
2854
2855 Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const 
2856 {
2857   // 
2858   // Returns number of SENSITIVE time bins
2859   //
2860
2861   Int_t tb, layer;
2862   for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
2863     layer = GetLayerNumber(tb);
2864     if(layer>=0) break;
2865   }
2866   return tb+1;
2867 }
2868
2869 //______________________________________________________
2870
2871 void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2872
2873   //
2874   // Insert layer <pl> in fLayers array.
2875   // Layers are sorted according to X coordinate.
2876
2877   if ( fN == ((Int_t) kMaxLayersPerSector)) {
2878     printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2879     return;
2880   }
2881   if (fN==0) {fLayers[fN++] = pl; return;}
2882   Int_t i=Find(pl->GetX());
2883
2884   memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2885   fLayers[i]=pl; fN++;
2886
2887 }              
2888
2889 //______________________________________________________
2890
2891 Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const 
2892 {
2893   //
2894   // Returns index of the propagation layer nearest to X 
2895   //
2896
2897   if (x <= fLayers[0]->GetX()) return 0;
2898   if (x > fLayers[fN-1]->GetX()) return fN;
2899   Int_t b=0, e=fN-1, m=(b+e)/2;
2900   for (; b<e; m=(b+e)/2) {
2901     if (x > fLayers[m]->GetX()) b=m+1;
2902     else e=m;
2903   }
2904   return m;
2905 }             
2906
2907 //______________________________________________________
2908 void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
2909 {
2910   //
2911   // set centers and the width of sectors
2912   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
2913     fZc[icham] = center[icham];  
2914     fZmax[icham] = w[icham];
2915     fZmaxSensitive[icham] = wsensitive[icham];
2916     //   printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
2917   }  
2918 }
2919 //______________________________________________________
2920
2921 void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
2922 {
2923   //
2924   // set centers and the width of sectors
2925   fHole = kFALSE;
2926   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
2927     fIsHole[icham] = holes[icham]; 
2928     if (holes[icham]) fHole = kTRUE;
2929   }  
2930 }
2931
2932
2933
2934 void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
2935         Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength, 
2936         Bool_t &lookForCluster) const
2937 {
2938   //
2939   // Returns radial step <dx>, density <rho>, rad. length <radLength>,
2940   // and sensitivity <lookForCluster> in point <y,z>  
2941   //
2942
2943   dx  = fdX;
2944   rho = fRho;
2945   radLength  = fX0;
2946   lookForCluster = kFALSE;
2947   //
2948   // check dead regions in sensitive volume 
2949   if(fTimeBinIndex >= 0) {
2950     //
2951     Int_t zone=-1;
2952     for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
2953       if  (TMath::Abs(z - fZc[ch]) < fZmaxSensitive[ch]){ 
2954         zone = ch;
2955         lookForCluster = !(fIsHole[zone]);
2956         if(TMath::Abs(y) > fYmaxSensitive){  
2957           lookForCluster = kFALSE;
2958         }
2959         if (fIsHole[zone]) {
2960           //if hole
2961           rho = 1.29e-3;
2962           radLength = 36.66;
2963         }
2964       }    
2965     }
2966     return;
2967   }
2968   //
2969   //
2970   // check hole
2971   if (fHole==kFALSE) return;
2972   //
2973   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
2974     if  (TMath::Abs(z - fZc[ch]) < fZmax[ch]){ 
2975       if (fIsHole[ch]) {
2976         //if hole
2977         rho = 1.29e-3;
2978         radLength = 36.66;
2979       }
2980     }
2981   }
2982   return;
2983 }
2984
2985 Int_t  AliTRDtracker::AliTRDpropagationLayer::GetZone( Double_t z) const
2986 {
2987   //
2988   //
2989   if (fTimeBinIndex < 0) return -20;  //unknown 
2990   Int_t zone=-10;   // dead zone
2991   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
2992     if(TMath::Abs(z - fZc[ch]) < fZmax[ch]) 
2993       zone = ch;
2994   }
2995   return zone;
2996 }
2997
2998
2999 //______________________________________________________
3000
3001 void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c, 
3002                                                           UInt_t index) {
3003
3004 // Insert cluster in cluster array.
3005 // Clusters are sorted according to Y coordinate.  
3006
3007   if(fTimeBinIndex < 0) { 
3008     printf("*** attempt to insert cluster into non-sensitive time bin!\n");
3009     return;
3010   }
3011
3012   if (fN== (Int_t) kMaxClusterPerTimeBin) {
3013     printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n"); 
3014     return;
3015   }
3016   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
3017   Int_t i=Find(c->GetY());
3018   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3019   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t)); 
3020   fIndex[i]=index; fClusters[i]=c; fN++;
3021 }  
3022
3023 //______________________________________________________
3024
3025 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
3026
3027 // Returns index of the cluster nearest in Y    
3028
3029   if (y <= fClusters[0]->GetY()) return 0;
3030   if (y > fClusters[fN-1]->GetY()) return fN;
3031   Int_t b=0, e=fN-1, m=(b+e)/2;
3032   for (; b<e; m=(b+e)/2) {
3033     if (y > fClusters[m]->GetY()) b=m+1;
3034     else e=m;
3035   }
3036   return m;
3037 }    
3038
3039 //---------------------------------------------------------
3040
3041 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
3042 //
3043 //  Returns correction factor for tilted pads geometry 
3044 //
3045
3046   Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
3047   Int_t det = c->GetDetector();    
3048   Int_t plane = fGeom->GetPlane(det);
3049
3050   if((plane == 1) || (plane == 3) || (plane == 5)) h01=-h01;
3051
3052   if(fNoTilt) h01 = 0;
3053   
3054   return h01;
3055 }
3056
3057
3058
3059
3060
3061