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