Update of tracking code for tilted pads
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15                                                       
16 /*
17 $Log$
18 Revision 1.19  2002/10/22 15:53:08  alibrary
19 Introducing Riostream.h
20
21 Revision 1.18  2002/10/14 14:57:44  hristov
22 Merging the VirtualMC branch to the main development branch (HEAD)
23
24 Revision 1.14.6.2  2002/07/24 10:09:31  alibrary
25 Updating VirtualMC
26
27 Revision 1.17  2002/06/13 12:09:58  hristov
28 Minor corrections
29
30 Revision 1.16  2002/06/12 09:54:36  cblume
31 Update of tracking code provided by Sergei
32
33 Revision 1.14  2001/11/14 10:50:46  cblume
34 Changes in digits IO. Add merging of summable digits
35
36 Revision 1.13  2001/05/30 12:17:47  hristov
37 Loop variables declared once
38
39 Revision 1.12  2001/05/28 17:07:58  hristov
40 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
41
42 Revision 1.8  2000/12/20 13:00:44  cblume
43 Modifications for the HP-compiler
44
45 Revision 1.7  2000/12/08 16:07:02  cblume
46 Update of the tracking by Sergei
47
48 Revision 1.6  2000/11/30 17:38:08  cblume
49 Changes to get in line with new STEER and EVGEN
50
51 Revision 1.5  2000/11/14 14:40:27  cblume
52 Correction for the Sun compiler (kTRUE and kFALSE)
53
54 Revision 1.4  2000/11/10 14:57:52  cblume
55 Changes in the geometry constants for the DEC compiler
56
57 Revision 1.3  2000/10/15 23:40:01  cblume
58 Remove AliTRDconst
59
60 Revision 1.2  2000/10/06 16:49:46  cblume
61 Made Getters const
62
63 Revision 1.1.2.2  2000/10/04 16:34:58  cblume
64 Replace include files by forward declarations
65
66 Revision 1.1.2.1  2000/09/22 14:47:52  cblume
67 Add the tracking code
68
69 */   
70
71 #include <Riostream.h>
72
73 #include <TFile.h>
74 #include <TBranch.h>
75 #include <TTree.h>  
76 #include <TObjArray.h> 
77
78 #include "AliTRDgeometry.h"
79 #include "AliTRDparameter.h"
80 #include "AliTRDgeometryDetail.h"
81 #include "AliTRDcluster.h" 
82 #include "AliTRDtrack.h"
83 #include "../TPC/AliTPCtrack.h"
84
85 #include "AliTRDtracker.h"
86
87 ClassImp(AliTRDtracker) 
88
89   const  Float_t     AliTRDtracker::fSeedDepth          = 0.5; 
90   const  Float_t     AliTRDtracker::fSeedStep           = 0.10;   
91   const  Float_t     AliTRDtracker::fSeedGap            = 0.25;  
92
93   const  Float_t     AliTRDtracker::fMaxSeedDeltaZ12    = 40.;  
94   const  Float_t     AliTRDtracker::fMaxSeedDeltaZ      = 25.;  
95   const  Float_t     AliTRDtracker::fMaxSeedC           = 0.0052; 
96   const  Float_t     AliTRDtracker::fMaxSeedTan         = 1.2;  
97   const  Float_t     AliTRDtracker::fMaxSeedVertexZ     = 150.; 
98
99   const  Double_t    AliTRDtracker::fSeedErrorSY        = 0.2;
100   const  Double_t    AliTRDtracker::fSeedErrorSY3       = 2.5;
101   const  Double_t    AliTRDtracker::fSeedErrorSZ        = 0.1;
102
103   const  Float_t     AliTRDtracker::fMinClustersInSeed  = 0.7;  
104
105   const  Float_t     AliTRDtracker::fMinClustersInTrack = 0.5;  
106   const  Float_t     AliTRDtracker::fMinFractionOfFoundClusters = 0.8;  
107
108   const  Float_t     AliTRDtracker::fSkipDepth          = 0.05;
109   const  Float_t     AliTRDtracker::fLabelFraction      = 0.8;  
110   const  Float_t     AliTRDtracker::fWideRoad           = 20.;
111
112   const  Double_t    AliTRDtracker::fMaxChi2            = 24.; 
113
114
115 //____________________________________________________________________
116 AliTRDtracker::AliTRDtracker(const TFile *geomfile)
117 {
118   // 
119   //  Main constructor
120   //  
121
122   Float_t fTzero = 0;
123   Int_t   version = 2;
124   
125   fAddTRDseeds = kFALSE;
126
127   fGeom = NULL;
128   
129   TDirectory *savedir=gDirectory; 
130   TFile *in=(TFile*)geomfile;  
131   if (!in->IsOpen()) {
132     printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
133     printf("    DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
134   }
135   else {
136     in->cd();  
137     in->ls();
138     fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
139     fPar  = (AliTRDparameter*) in->Get("TRDparameter");
140     fGeom->Dump();
141   }
142
143   if(fGeom) {
144     //    fTzero = geo->GetT0();
145     fTzero = 0.;
146     version = fGeom->IsVersion();
147     printf("Found geometry version %d on file \n", version);
148   }
149   else { 
150     printf("AliTRDtracker::AliTRDtracker(): cann't find TRD geometry!\n");
151     printf("    DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
152     fGeom = new AliTRDgeometryDetail(); 
153     fPar = new AliTRDparameter();
154   }
155
156   savedir->cd();  
157
158
159   //  fGeom->SetT0(fTzero);
160
161   fEvent     = 0;
162
163   fNclusters = 0;
164   fClusters  = new TObjArray(2000); 
165   fNseeds    = 0;
166   fSeeds     = new TObjArray(2000);
167   fNtracks   = 0;
168   fTracks    = new TObjArray(1000);
169
170   for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
171     Int_t tr_s = CookSectorIndex(geom_s);
172     fTrSec[tr_s] = new AliTRDtrackingSector(fGeom, geom_s, fPar);
173   }
174
175   fSY2corr = 0.025;
176   fSZ2corr = 12.;      
177
178   // calculate max gap on track
179
180   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
181   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
182
183   Double_t dx = (Double_t) fPar->GetTimeBinSize();   
184   Int_t tbAmp = fPar->GetTimeBefore();
185   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
186   Int_t tbDrift = fPar->GetTimeMax();
187   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
188
189   tbDrift = TMath::Min(tbDrift,maxDrift);
190   tbAmp = TMath::Min(tbAmp,maxAmp);
191
192   fTimeBinsPerPlane = tbAmp + tbDrift;
193   fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fSkipDepth);
194
195   fVocal = kFALSE;
196
197 }   
198
199 //___________________________________________________________________
200 AliTRDtracker::~AliTRDtracker()
201 {
202   delete fClusters;
203   delete fTracks;
204   delete fSeeds;
205   delete fGeom;  
206   delete fPar;  
207
208   for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
209     delete fTrSec[geom_s];
210   }
211 }   
212
213 //_____________________________________________________________________
214 inline Double_t f1trd(Double_t x1,Double_t y1,
215                       Double_t x2,Double_t y2,
216                       Double_t x3,Double_t y3)
217 {
218   //
219   // Initial approximation of the track curvature
220   //
221   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
222   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
223                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
224   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
225                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
226
227   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
228
229   return -xr*yr/sqrt(xr*xr+yr*yr);
230 }          
231
232 //_____________________________________________________________________
233 inline Double_t f2trd(Double_t x1,Double_t y1,
234                       Double_t x2,Double_t y2,
235                       Double_t x3,Double_t y3)
236 {
237   //
238   // Initial approximation of the track curvature times X coordinate
239   // of the center of curvature
240   //
241
242   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
243   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
244                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
245   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
246                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
247
248   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
249
250   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
251 }          
252
253 //_____________________________________________________________________
254 inline Double_t f3trd(Double_t x1,Double_t y1,
255                       Double_t x2,Double_t y2,
256                       Double_t z1,Double_t z2)
257 {
258   //
259   // Initial approximation of the tangent of the track dip angle
260   //
261
262   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
263 }            
264
265 //___________________________________________________________________
266 Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out)
267 {
268   //
269   // Finds tracks within the TRD. File <inp> is expected to contain seeds 
270   // at the outer part of the TRD. If <inp> is NULL, the seeds
271   // are found within the TRD if fAddTRDseeds is TRUE. 
272   // The tracks are propagated to the innermost time bin 
273   // of the TRD and stored in file <out>. 
274   //
275
276   LoadEvent();
277  
278   TDirectory *savedir=gDirectory;
279
280   char   tname[100];
281
282   if (!out->IsOpen()) {
283     cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n";
284     return 1;
285   }    
286
287   sprintf(tname,"seedTRDtoTPC_%d",fEvent); 
288   TTree tpc_tree(tname,"Tree with seeds from TRD at outer TPC pad row");
289   AliTPCtrack *iotrack=0;
290   tpc_tree.Branch("tracks","AliTPCtrack",&iotrack,32000,0); 
291
292   sprintf(tname,"TreeT%d_TRD",fEvent);
293   TTree trd_tree(tname,"TRD tracks at inner TRD time bin");
294   AliTRDtrack *iotrack_trd=0;
295   trd_tree.Branch("tracks","AliTRDtrack",&iotrack_trd,32000,0);  
296
297   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
298   Float_t foundMin = fMinClustersInTrack * timeBins; 
299
300   if (inp) {
301      TFile *in=(TFile*)inp;
302      if (!in->IsOpen()) {
303         cerr<<"AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n";
304         cerr<<" ... going for seeds finding inside the TRD\n";
305      }
306      else {
307        in->cd();
308        sprintf(tname,"TRDb_%d",fEvent);  
309        TTree *seedTree=(TTree*)in->Get(tname);  
310        if (!seedTree) {
311          cerr<<"AliTRDtracker::Clusters2Tracks(): ";
312          cerr<<"can't get a tree with track seeds !\n";
313          return 3;
314        }  
315        AliTRDtrack *seed=new AliTRDtrack;
316        seedTree->SetBranchAddress("tracks",&seed);
317
318        Int_t n=(Int_t)seedTree->GetEntries();
319        for (Int_t i=0; i<n; i++) {
320          seedTree->GetEvent(i);
321          seed->ResetCovariance(); 
322          AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
323          fSeeds->AddLast(tr);
324          fNseeds++;
325        }          
326        delete seed;
327        delete seedTree;
328      }
329   }
330
331   out->cd();
332
333
334
335   // find tracks from loaded seeds
336
337   Int_t nseed=fSeeds->GetEntriesFast();
338   Int_t i, found = 0;
339   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
340
341   for (i=0; i<nseed; i++) {   
342     AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt; 
343     FollowProlongation(t, innerTB); 
344     if (t.GetNumberOfClusters() >= foundMin) {
345       UseClusters(&t);
346       CookLabel(pt, 1-fLabelFraction);
347       //      t.CookdEdx();
348     }
349     iotrack_trd = pt;
350     trd_tree.Fill();
351     cout<<found++<<'\r';     
352
353     if(PropagateToTPC(t)) {
354       AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
355       iotrack = tpc;
356       tpc_tree.Fill();
357       delete tpc;
358     }  
359     delete fSeeds->RemoveAt(i);
360     fNseeds--;
361   }     
362
363   cout<<"Number of loaded seeds: "<<nseed<<endl;  
364   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
365
366   // after tracks from loaded seeds are found and the corresponding 
367   // clusters are used, look for additional seeds from TRD
368
369   if(fAddTRDseeds) { 
370     // Find tracks for the seeds in the TRD
371     Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
372   
373     Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep);
374     Int_t gap = (Int_t) (timeBins * fSeedGap);
375     Int_t step = (Int_t) (timeBins * fSeedStep);
376   
377     // make a first turn with tight cut on initial curvature
378     for(Int_t turn = 1; turn <= 2; turn++) {
379       if(turn == 2) {
380         nSteps = (Int_t) (fSeedDepth / (3*fSeedStep));
381         step = (Int_t) (timeBins * (3*fSeedStep));
382       }
383       for(Int_t i=0; i<nSteps; i++) {
384         Int_t outer=timeBins-1-i*step; 
385         Int_t inner=outer-gap;
386
387         nseed=fSeeds->GetEntriesFast();
388       
389         MakeSeeds(inner, outer, turn);
390       
391         nseed=fSeeds->GetEntriesFast();
392         printf("\n turn %d, step %d: number of seeds for TRD inward %d\n", 
393                turn, i, nseed); 
394               
395         for (Int_t i=0; i<nseed; i++) {   
396           AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt; 
397           FollowProlongation(t,innerTB); 
398           if (t.GetNumberOfClusters() >= foundMin) {
399             UseClusters(&t);
400             CookLabel(pt, 1-fLabelFraction);
401             t.CookdEdx();
402             cout<<found++<<'\r';     
403             iotrack_trd = pt;
404             trd_tree.Fill();
405             if(PropagateToTPC(t)) {
406               AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
407               iotrack = tpc;
408               tpc_tree.Fill();
409               delete tpc;
410             }   
411           }
412           delete fSeeds->RemoveAt(i);
413           fNseeds--;
414         }
415       }
416     }
417   }
418   tpc_tree.Write(); 
419   trd_tree.Write(); 
420   
421   cout<<"Total number of found tracks: "<<found<<endl;
422     
423   UnloadEvent();
424     
425   savedir->cd();  
426   
427   return 0;    
428 }     
429      
430   
431
432 //_____________________________________________________________________________
433 Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
434   //
435   // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and
436   // backpropagated by the TPC tracker. Each seed is first propagated 
437   // to the TRD, and then its prolongation is searched in the TRD.
438   // If sufficiently long continuation of the track is found in the TRD
439   // the track is updated, otherwise it's stored as originaly defined 
440   // by the TPC tracker.   
441   //  
442
443
444   LoadEvent();
445
446   TDirectory *savedir=gDirectory;
447
448   TFile *in=(TFile*)inp;
449
450   if (!in->IsOpen()) {
451      cerr<<"AliTRDtracker::PropagateBack(): ";
452      cerr<<"file with back propagated TPC tracks is not open !\n";
453      return 1;
454   }                   
455
456   if (!out->IsOpen()) {
457      cerr<<"AliTRDtracker::PropagateBack(): ";
458      cerr<<"file for back propagated TRD tracks is not open !\n";
459      return 2;
460   }      
461
462   in->cd();
463   char   tname[100];
464   sprintf(tname,"seedsTPCtoTRD_%d",fEvent);       
465   TTree *seedTree=(TTree*)in->Get(tname);
466   if (!seedTree) {
467      cerr<<"AliTRDtracker::PropagateBack(): ";
468      cerr<<"can't get a tree with seeds from TPC !\n";
469      cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n";
470      return 3;
471   }
472
473   AliTPCtrack *seed=new AliTPCtrack;
474   seedTree->SetBranchAddress("tracks",&seed);
475
476   Int_t n=(Int_t)seedTree->GetEntries();
477   for (Int_t i=0; i<n; i++) {
478      seedTree->GetEvent(i);
479      Int_t lbl = seed->GetLabel();
480      AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
481      tr->SetSeedLabel(lbl);
482      fSeeds->AddLast(tr);
483      fNseeds++;
484   }
485
486   delete seed;
487   delete seedTree;
488
489   out->cd();
490
491   AliTPCtrack *otrack=0;
492
493   sprintf(tname,"seedsTRDtoTOF1_%d",fEvent);  
494   TTree tofTree1(tname,"Tracks back propagated through TPC and TRD");
495   tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
496
497   sprintf(tname,"seedsTRDtoTOF2_%d",fEvent);  
498   TTree tofTree2(tname,"Tracks back propagated through TPC and TRD");
499   tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
500
501   sprintf(tname,"seedsTRDtoPHOS_%d",fEvent);  
502   TTree phosTree(tname,"Tracks back propagated through TPC and TRD");
503   phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
504
505   sprintf(tname,"seedsTRDtoRICH_%d",fEvent);  
506   TTree richTree(tname,"Tracks back propagated through TPC and TRD");
507   richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);  
508
509   sprintf(tname,"TRDb_%d",fEvent);  
510   TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin");
511   AliTRDtrack *otrack_trd=0;
512   trdTree.Branch("tracks","AliTRDtrack",&otrack_trd,32000,0);   
513      
514   Int_t found=0;  
515
516   Int_t nseed=fSeeds->GetEntriesFast();
517
518   Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan(); 
519
520   Int_t outermost_tb  = fTrSec[0]->GetOuterTimeBin();
521
522   for (Int_t i=0; i<nseed; i++) {  
523
524     AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
525     Int_t expectedClr = FollowBackProlongation(s);
526     Int_t foundClr = s.GetNumberOfClusters();
527     Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX());
528
529     //    printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
530     //     i, foundClr, expectedClr, foundMin);
531
532     if (foundClr >= foundMin) {
533       s.CookdEdx(); 
534       CookLabel(ps, 1-fLabelFraction);
535       UseClusters(ps);
536       otrack_trd=ps;
537       trdTree.Fill();
538       cout<<found++<<'\r';
539     }
540
541     if(((expectedClr < 10) && (last_tb == outermost_tb)) ||
542        ((expectedClr >= 10) && 
543         (((Float_t) foundClr) / ((Float_t) expectedClr) >= 
544          fMinFractionOfFoundClusters) && (last_tb == outermost_tb))) {
545
546       Double_t x_tof = 375.5;
547     
548       if(PropagateToOuterPlane(s,x_tof)) {
549         AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
550         otrack = pt;
551         tofTree1.Fill();
552         delete pt;
553
554         x_tof = 381.5;
555     
556         if(PropagateToOuterPlane(s,x_tof)) {
557           AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
558           otrack = pt;
559           tofTree2.Fill();
560           delete pt;
561
562           Double_t x_phos = 460.;
563           
564           if(PropagateToOuterPlane(s,x_phos)) {
565             AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
566             otrack = pt;
567             phosTree.Fill();
568             delete pt;
569             
570             Double_t x_rich = 490+1.267;
571             
572             if(PropagateToOuterPlane(s,x_rich)) {
573               AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
574               otrack = pt;
575               richTree.Fill();
576               delete pt;
577             }   
578           }
579         }
580       }      
581     }
582   }
583   
584   tofTree1.Write(); 
585   tofTree2.Write(); 
586   phosTree.Write(); 
587   richTree.Write(); 
588   trdTree.Write(); 
589
590   savedir->cd();  
591   cerr<<"Number of seeds: "<<nseed<<endl;  
592   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
593
594   UnloadEvent();
595
596   return 0;
597
598 }
599
600
601 //---------------------------------------------------------------------------
602 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
603 {
604   // Starting from current position on track=t this function tries
605   // to extrapolate the track up to timeBin=0 and to confirm prolongation
606   // if a close cluster is found. Returns the number of clusters
607   // expected to be found in sensitive layers
608
609   Float_t  wIndex, wTB, wChi2;
610   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
611   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
612   Float_t  wPx, wPy, wPz, wC;
613   Double_t Px, Py, Pz;
614   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
615
616   Int_t trackIndex = t.GetLabel();  
617
618   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
619
620   Int_t try_again=fMaxGap;
621
622   Double_t alpha=t.GetAlpha();
623
624   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
625   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
626
627   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
628
629   Double_t x0, rho, x, dx, y, ymax, z;
630
631   Int_t expectedNumberOfClusters = 0;
632   Bool_t lookForCluster;
633
634   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
635
636  
637   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) { 
638
639     y = t.GetY(); z = t.GetZ();
640
641     // first propagate to the inner surface of the current time bin 
642     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
643     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
644     if(!t.PropagateTo(x,x0,rho,0.139)) break;
645     y = t.GetY();
646     ymax = x*TMath::Tan(0.5*alpha);
647     if (y > ymax) {
648       s = (s+1) % ns;
649       if (!t.Rotate(alpha)) break;
650       if(!t.PropagateTo(x,x0,rho,0.139)) break;
651     } else if (y <-ymax) {
652       s = (s-1+ns) % ns;                           
653       if (!t.Rotate(-alpha)) break;   
654       if(!t.PropagateTo(x,x0,rho,0.139)) break;
655     } 
656
657     y = t.GetY(); z = t.GetZ();
658
659     // now propagate to the middle plane of the next time bin 
660     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
661     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
662     if(!t.PropagateTo(x,x0,rho,0.139)) break;
663     y = t.GetY();
664     ymax = x*TMath::Tan(0.5*alpha);
665     if (y > ymax) {
666       s = (s+1) % ns;
667       if (!t.Rotate(alpha)) break;
668       if(!t.PropagateTo(x,x0,rho,0.139)) break;
669     } else if (y <-ymax) {
670       s = (s-1+ns) % ns;                           
671       if (!t.Rotate(-alpha)) break;   
672       if(!t.PropagateTo(x,x0,rho,0.139)) break;
673     } 
674
675
676     if(lookForCluster) {
677
678
679       expectedNumberOfClusters++;       
680
681       wIndex = (Float_t) t.GetLabel();
682       wTB = nr;
683
684
685       AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr-1));
686
687
688       Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
689
690
691       Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
692
693       Double_t road;
694       if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
695       else return expectedNumberOfClusters;
696       
697       wYrt = (Float_t) y;
698       wZrt = (Float_t) z;
699       wYwindow = (Float_t) road;
700       t.GetPxPyPz(Px,Py,Pz);
701       wPx = (Float_t) Px;
702       wPy = (Float_t) Py;
703       wPz = (Float_t) Pz;
704       wC  = (Float_t) t.GetC();
705       wSigmaC2 = (Float_t) t.GetSigmaC2();
706       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
707       wSigmaY2 = (Float_t) t.GetSigmaY2();
708       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
709       wChi2 = -1;            
710       
711       if (road>fWideRoad) {
712         if (t.GetNumberOfClusters()>4)
713           cerr<<t.GetNumberOfClusters()
714               <<"FindProlongation warning: Too broad road !\n";
715         return 0;
716       }             
717
718       AliTRDcluster *cl=0;
719       UInt_t index=0;
720
721       Double_t max_chi2=fMaxChi2;
722
723       wYclosest = 12345678;
724       wYcorrect = 12345678;
725       wZclosest = 12345678;
726       wZcorrect = 12345678;
727       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);   
728
729       // Find the closest correct cluster for debugging purposes
730       if (time_bin) {
731         Float_t minDY = 1000000;
732         for (Int_t i=0; i<time_bin; i++) {
733           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
734           if((c->GetLabel(0) != trackIndex) &&
735              (c->GetLabel(1) != trackIndex) &&
736              (c->GetLabel(2) != trackIndex)) continue;
737           if(TMath::Abs(c->GetY() - y) > minDY) continue;
738           minDY = TMath::Abs(c->GetY() - y);
739           wYcorrect = c->GetY();
740           wZcorrect = c->GetZ();
741
742           Double_t h01 = GetTiltFactor(c);
743           wChi2 = t.GetPredictedChi2(c, h01);
744         }
745       }                    
746
747       // Now go for the real cluster search
748
749       if (time_bin) {
750
751         for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
752           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
753           if (c->GetY() > y+road) break;
754           if (c->IsUsed() > 0) continue;
755           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
756
757           Double_t h01 = GetTiltFactor(c);
758           Double_t chi2=t.GetPredictedChi2(c,h01);
759           
760           if (chi2 > max_chi2) continue;
761           max_chi2=chi2;
762           cl=c;
763           index=time_bin.GetIndex(i);
764         }               
765
766         if(!cl) {
767
768           for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
769             AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
770             
771             if (c->GetY() > y+road) break;
772             if (c->IsUsed() > 0) continue;
773             if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
774             
775             Double_t h01 = GetTiltFactor(c);
776             Double_t chi2=t.GetPredictedChi2(c, h01);
777             
778             if (chi2 > max_chi2) continue;
779             max_chi2=chi2;
780             cl=c;
781             index=time_bin.GetIndex(i);
782           }
783         }       
784         
785
786         if (cl) {
787           wYclosest = cl->GetY();
788           wZclosest = cl->GetZ();
789           Double_t h01 = GetTiltFactor(cl);
790
791           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); 
792           if(!t.Update(cl,max_chi2,index,h01)) {
793             if(!try_again--) return 0;
794           }  
795           else try_again=fMaxGap;
796         }
797         else {
798           if (try_again==0) break; 
799           try_again--;
800         }
801
802         /*
803         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
804           
805           printf(" %f", wIndex);       //1
806           printf(" %f", wTB);          //2
807           printf(" %f", wYrt);         //3
808           printf(" %f", wYclosest);    //4
809           printf(" %f", wYcorrect);    //5
810           printf(" %f", wYwindow);     //6
811           printf(" %f", wZrt);         //7
812           printf(" %f", wZclosest);    //8
813           printf(" %f", wZcorrect);    //9
814           printf(" %f", wZwindow);     //10
815           printf(" %f", wPx);          //11
816           printf(" %f", wPy);          //12
817           printf(" %f", wPz);          //13
818           printf(" %f", wSigmaC2*1000000);  //14
819           printf(" %f", wSigmaTgl2*1000);   //15
820           printf(" %f", wSigmaY2);     //16
821           //      printf(" %f", wSigmaZ2);     //17
822           printf(" %f", wChi2);     //17
823           printf(" %f", wC);           //18
824           printf("\n");
825         } 
826         */                        
827       }
828     }  
829   }
830   return expectedNumberOfClusters;
831   
832   
833 }                
834
835 //___________________________________________________________________
836
837 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
838 {
839   // Starting from current radial position of track <t> this function
840   // extrapolates the track up to outer timebin and in the sensitive
841   // layers confirms prolongation if a close cluster is found. 
842   // Returns the number of clusters expected to be found in sensitive layers
843
844   Float_t  wIndex, wTB, wChi2;
845   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
846   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
847   Float_t  wPx, wPy, wPz, wC;
848   Double_t Px, Py, Pz;
849   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
850
851   Int_t trackIndex = t.GetLabel();  
852
853   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
854
855   Int_t try_again=fMaxGap;
856
857   Double_t alpha=t.GetAlpha();
858
859   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
860   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
861
862   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
863
864   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
865
866   Double_t x0, rho, x, dx, y, ymax, z;
867
868   Bool_t lookForCluster;
869
870   Int_t expectedNumberOfClusters = 0;
871   x = t.GetX();
872
873   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
874
875
876   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB; nr++) { 
877
878     y = t.GetY(); z = t.GetZ();
879
880     // first propagate to the outer surface of the current time bin 
881
882     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
883     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
884
885     if(!t.PropagateTo(x,x0,rho,0.139)) break;
886
887     y = t.GetY();
888     ymax = x*TMath::Tan(0.5*alpha);
889
890     if (y > ymax) {
891       s = (s+1) % ns;
892       if (!t.Rotate(alpha)) break;
893       if(!t.PropagateTo(x,x0,rho,0.139)) break;
894     } else if (y <-ymax) {
895       s = (s-1+ns) % ns;                           
896       if (!t.Rotate(-alpha)) break;   
897       if(!t.PropagateTo(x,x0,rho,0.139)) break;
898     } 
899     y = t.GetY(); z = t.GetZ();
900
901     // now propagate to the middle plane of the next time bin 
902     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
903
904     x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
905
906     if(!t.PropagateTo(x,x0,rho,0.139)) break;
907
908     y = t.GetY();
909
910     ymax = x*TMath::Tan(0.5*alpha);
911
912     if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
913
914     if (y > ymax) {
915       s = (s+1) % ns;
916       if (!t.Rotate(alpha)) break;
917       if(!t.PropagateTo(x,x0,rho,0.139)) break;
918     } else if (y <-ymax) {
919       s = (s-1+ns) % ns;              
920       if (!t.Rotate(-alpha)) break;   
921       if(!t.PropagateTo(x,x0,rho,0.139)) break;
922     } 
923
924     //    printf("label %d, pl %d, lookForCluster %d \n",
925     //     trackIndex, nr+1, lookForCluster);
926
927     if(lookForCluster) {
928       expectedNumberOfClusters++;       
929
930       wIndex = (Float_t) t.GetLabel();
931       wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
932
933       AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr+1));
934       Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
935       Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
936       if((t.GetSigmaY2() + sy2) < 0) break;
937       Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2); 
938       Double_t y=t.GetY(), z=t.GetZ();
939
940       wYrt = (Float_t) y;
941       wZrt = (Float_t) z;
942       wYwindow = (Float_t) road;
943       t.GetPxPyPz(Px,Py,Pz);
944       wPx = (Float_t) Px;
945       wPy = (Float_t) Py;
946       wPz = (Float_t) Pz;
947       wC  = (Float_t) t.GetC();
948       wSigmaC2 = (Float_t) t.GetSigmaC2();
949       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
950       wSigmaY2 = (Float_t) t.GetSigmaY2();
951       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
952       wChi2 = -1;            
953       
954       if (road>fWideRoad) {
955         if (t.GetNumberOfClusters()>4)
956           cerr<<t.GetNumberOfClusters()
957               <<"FindProlongation warning: Too broad road !\n";
958         return 0;
959       }      
960
961       AliTRDcluster *cl=0;
962       UInt_t index=0;
963
964       Double_t max_chi2=fMaxChi2;
965
966       wYclosest = 12345678;
967       wYcorrect = 12345678;
968       wZclosest = 12345678;
969       wZcorrect = 12345678;
970       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);   
971
972       // Find the closest correct cluster for debugging purposes
973       if (time_bin) {
974         Float_t minDY = 1000000;
975         for (Int_t i=0; i<time_bin; i++) {
976           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
977           if((c->GetLabel(0) != trackIndex) &&
978              (c->GetLabel(1) != trackIndex) &&
979              (c->GetLabel(2) != trackIndex)) continue;
980           if(TMath::Abs(c->GetY() - y) > minDY) continue;
981           minDY = TMath::Abs(c->GetY() - y);
982           wYcorrect = c->GetY();
983           wZcorrect = c->GetZ();
984
985           Double_t h01 = GetTiltFactor(c);
986           wChi2 = t.GetPredictedChi2(c, h01);
987         }
988       }                    
989
990       // Now go for the real cluster search
991
992       if (time_bin) {
993
994         for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
995           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
996           if (c->GetY() > y+road) break;
997           if (c->IsUsed() > 0) continue;
998           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
999
1000           Double_t h01 = GetTiltFactor(c);
1001           Double_t chi2=t.GetPredictedChi2(c,h01);
1002           
1003           if (chi2 > max_chi2) continue;
1004           max_chi2=chi2;
1005           cl=c;
1006           index=time_bin.GetIndex(i);
1007         }               
1008
1009         if(!cl) {
1010
1011           for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
1012             AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
1013             
1014             if (c->GetY() > y+road) break;
1015             if (c->IsUsed() > 0) continue;
1016             if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
1017             
1018             Double_t h01 = GetTiltFactor(c);
1019             Double_t chi2=t.GetPredictedChi2(c,h01);
1020             
1021             if (chi2 > max_chi2) continue;
1022             max_chi2=chi2;
1023             cl=c;
1024             index=time_bin.GetIndex(i);
1025           }
1026         }       
1027         
1028         if (cl) {
1029           wYclosest = cl->GetY();
1030           wZclosest = cl->GetZ();
1031
1032           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); 
1033           Double_t h01 = GetTiltFactor(cl);
1034           if(!t.Update(cl,max_chi2,index,h01)) {
1035             if(!try_again--) return 0;
1036           }  
1037           else try_again=fMaxGap;
1038         }
1039         else {
1040           if (try_again==0) break; 
1041           try_again--;
1042         }
1043
1044         /*
1045         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1046           
1047           printf(" %f", wIndex);       //1
1048           printf(" %f", wTB);          //2
1049           printf(" %f", wYrt);         //3
1050           printf(" %f", wYclosest);    //4
1051           printf(" %f", wYcorrect);    //5
1052           printf(" %f", wYwindow);     //6
1053           printf(" %f", wZrt);         //7
1054           printf(" %f", wZclosest);    //8
1055           printf(" %f", wZcorrect);    //9
1056           printf(" %f", wZwindow);     //10
1057           printf(" %f", wPx);          //11
1058           printf(" %f", wPy);          //12
1059           printf(" %f", wPz);          //13
1060           printf(" %f", wSigmaC2*1000000);  //14
1061           printf(" %f", wSigmaTgl2*1000);   //15
1062           printf(" %f", wSigmaY2);     //16
1063           //      printf(" %f", wSigmaZ2);     //17
1064           printf(" %f", wChi2);     //17
1065           printf(" %f", wC);           //18
1066           printf("\n");
1067         } 
1068         */                        
1069       }
1070     }  
1071   }
1072   return expectedNumberOfClusters;
1073 }         
1074
1075 //___________________________________________________________________
1076
1077 Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1078 {
1079   // Starting from current radial position of track <t> this function
1080   // extrapolates the track up to radial position <xToGo>. 
1081   // Returns 1 if track reaches the plane, and 0 otherwise 
1082
1083   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
1084
1085   Double_t alpha=t.GetAlpha();
1086
1087   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1088   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1089
1090   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
1091
1092   Bool_t lookForCluster;
1093   Double_t x0, rho, x, dx, y, ymax, z;
1094
1095   x = t.GetX();
1096
1097   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1098
1099   Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1100
1101   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) { 
1102
1103     y = t.GetY(); z = t.GetZ();
1104
1105     // first propagate to the outer surface of the current time bin 
1106     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1107     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
1108     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1109     y = t.GetY();
1110     ymax = x*TMath::Tan(0.5*alpha);
1111     if (y > ymax) {
1112       s = (s+1) % ns;
1113       if (!t.Rotate(alpha)) return 0;
1114     } else if (y <-ymax) {
1115       s = (s-1+ns) % ns;                           
1116       if (!t.Rotate(-alpha)) return 0;   
1117     } 
1118     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1119
1120     y = t.GetY(); z = t.GetZ();
1121
1122     // now propagate to the middle plane of the next time bin 
1123     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1124     x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
1125     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1126     y = t.GetY();
1127     ymax = x*TMath::Tan(0.5*alpha);
1128     if (y > ymax) {
1129       s = (s+1) % ns;
1130       if (!t.Rotate(alpha)) return 0;
1131     } else if (y <-ymax) {
1132       s = (s-1+ns) % ns;                           
1133       if (!t.Rotate(-alpha)) return 0;   
1134     } 
1135     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1136   }
1137   return 1;
1138 }         
1139
1140 //___________________________________________________________________
1141
1142 Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1143 {
1144   // Starting from current radial position of track <t> this function
1145   // extrapolates the track up to radial position of the outermost
1146   // padrow of the TPC. 
1147   // Returns 1 if track reaches the TPC, and 0 otherwise 
1148
1149   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
1150
1151   Double_t alpha=t.GetAlpha();
1152
1153   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1154   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1155
1156   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
1157
1158   Bool_t lookForCluster;
1159   Double_t x0, rho, x, dx, y, ymax, z;
1160
1161   x = t.GetX();
1162
1163   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1164
1165   Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1166
1167   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) { 
1168
1169     y = t.GetY(); z = t.GetZ();
1170
1171     // first propagate to the outer surface of the current time bin 
1172     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1173     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1174     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1175     y = t.GetY();
1176     ymax = x*TMath::Tan(0.5*alpha);
1177     if (y > ymax) {
1178       s = (s+1) % ns;
1179       if (!t.Rotate(alpha)) return 0;
1180     } else if (y <-ymax) {
1181       s = (s-1+ns) % ns;                           
1182       if (!t.Rotate(-alpha)) return 0;   
1183     } 
1184     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1185
1186     y = t.GetY(); z = t.GetZ();
1187
1188     // now propagate to the middle plane of the next time bin 
1189     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,x0,lookForCluster);
1190     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1191     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1192     y = t.GetY();
1193     ymax = x*TMath::Tan(0.5*alpha);
1194     if (y > ymax) {
1195       s = (s+1) % ns;
1196       if (!t.Rotate(alpha)) return 0;
1197     } else if (y <-ymax) {
1198       s = (s-1+ns) % ns;                           
1199       if (!t.Rotate(-alpha)) return 0;   
1200     } 
1201     if(!t.PropagateTo(x,x0,rho,0.139)) return 0;
1202   }
1203   return 1;
1204 }         
1205
1206
1207 //_____________________________________________________________________________
1208 void AliTRDtracker::LoadEvent()
1209 {
1210   // Fills clusters into TRD tracking_sectors 
1211   // Note that the numbering scheme for the TRD tracking_sectors 
1212   // differs from that of TRD sectors
1213
1214   ReadClusters(fClusters);
1215   Int_t ncl=fClusters->GetEntriesFast();
1216   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1217               
1218   UInt_t index;
1219   while (ncl--) {
1220     printf("\r %d left  ",ncl); 
1221     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1222     Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
1223     Int_t sector=fGeom->GetSector(detector);
1224     Int_t plane=fGeom->GetPlane(detector);
1225
1226     Int_t tracking_sector = CookSectorIndex(sector);
1227
1228     Int_t gtb = fTrSec[tracking_sector]->CookTimeBinIndex(plane,local_time_bin);
1229     Int_t layer = fTrSec[tracking_sector]->GetLayerNumber(gtb);
1230
1231     index=ncl;
1232     fTrSec[tracking_sector]->GetLayer(layer)->InsertCluster(c,index);
1233   }    
1234   printf("\r\n");
1235
1236 }
1237
1238 //_____________________________________________________________________________
1239 void AliTRDtracker::UnloadEvent() 
1240
1241   //
1242   // Clears the arrays of clusters and tracks. Resets sectors and timebins 
1243   //
1244
1245   Int_t i, nentr;
1246
1247   nentr = fClusters->GetEntriesFast();
1248   for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1249
1250   nentr = fSeeds->GetEntriesFast();
1251   for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1252
1253   nentr = fTracks->GetEntriesFast();
1254   for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1255
1256   Int_t nsec = AliTRDgeometry::kNsect;
1257
1258   for (i = 0; i < nsec; i++) {    
1259     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1260       fTrSec[i]->GetLayer(pl)->Clear();
1261     }
1262   }
1263
1264 }
1265
1266 //__________________________________________________________________________
1267 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1268 {
1269   // Creates track seeds using clusters in timeBins=i1,i2
1270
1271   if(turn > 2) {
1272     cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
1273     return;
1274   }
1275
1276   Double_t x[5], c[15];
1277   Int_t max_sec=AliTRDgeometry::kNsect;
1278   
1279   Double_t alpha=AliTRDgeometry::GetAlpha();
1280   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1281   Double_t cs=cos(alpha), sn=sin(alpha);
1282   Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1283     
1284       
1285   Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1286   Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1287       
1288   Double_t x1 =fTrSec[0]->GetX(i1);
1289   Double_t xx2=fTrSec[0]->GetX(i2);
1290       
1291   for (Int_t ns=0; ns<max_sec; ns++) {
1292     
1293     Int_t nl2 = *(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1294     Int_t nl=(*fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1295     Int_t nm=(*fTrSec[ns]->GetLayer(i2));
1296     Int_t nu=(*fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1297     Int_t nu2=(*fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1298     
1299     AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1300     
1301     for (Int_t is=0; is < r1; is++) {
1302       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
1303       
1304       for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
1305         
1306         const AliTRDcluster *cl;
1307         Double_t x2,   y2,   z2;
1308         Double_t x3=0., y3=0.;   
1309         
1310         if (js<nl2) {
1311           if(turn != 2) continue;
1312           AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1313           cl=r2[js];
1314           y2=cl->GetY(); z2=cl->GetZ();
1315           
1316           x2= xx2*cs2+y2*sn2;
1317           y2=-xx2*sn2+y2*cs2;
1318         }
1319         else if (js<nl2+nl) {
1320           if(turn != 1) continue;
1321           AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1322           cl=r2[js-nl2];
1323           y2=cl->GetY(); z2=cl->GetZ();
1324           
1325           x2= xx2*cs+y2*sn;
1326           y2=-xx2*sn+y2*cs;
1327         }                                
1328         else if (js<nl2+nl+nm) {
1329           if(turn != 1) continue;
1330           AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1331           cl=r2[js-nl2-nl];
1332           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1333         }
1334         else if (js<nl2+nl+nm+nu) {
1335           if(turn != 1) continue;
1336           AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1337           cl=r2[js-nl2-nl-nm];
1338           y2=cl->GetY(); z2=cl->GetZ();
1339           
1340           x2=xx2*cs-y2*sn;
1341           y2=xx2*sn+y2*cs;
1342         }              
1343         else {
1344           if(turn != 2) continue;
1345           AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1346           cl=r2[js-nl2-nl-nm-nu];
1347           y2=cl->GetY(); z2=cl->GetZ();
1348           
1349           x2=xx2*cs2-y2*sn2;
1350           y2=xx2*sn2+y2*cs2;
1351         }
1352         
1353         if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue;
1354         
1355         Double_t zz=z1 - z1/x1*(x1-x2);
1356         
1357         if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
1358         
1359         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1360         if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
1361         
1362         x[0]=y1;
1363         x[1]=z1;
1364         x[2]=f1trd(x1,y1,x2,y2,x3,y3);
1365         
1366         if (TMath::Abs(x[2]) > fMaxSeedC) continue;      
1367         
1368         x[3]=f2trd(x1,y1,x2,y2,x3,y3);
1369         
1370         if (TMath::Abs(x[2]*x1-x[3]) >= 0.99999) continue;
1371         
1372         x[4]=f3trd(x1,y1,x2,y2,z1,z2);
1373         
1374         if (TMath::Abs(x[4]) > fMaxSeedTan) continue;
1375         
1376         Double_t a=asin(x[3]);
1377         Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
1378         
1379         if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
1380         
1381         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
1382         Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2();
1383         Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;  
1384
1385         // Tilt changes
1386         Double_t h01 = GetTiltFactor(r1[is]);
1387         sy1=sy1+sz1*h01*h01;
1388         Double_t syz=sz1*h01;
1389         // end of tilt changes
1390         
1391         Double_t f20=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1392         Double_t f22=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1393         Double_t f24=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1394         Double_t f30=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[3])/sy;
1395         Double_t f32=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[3])/sy;
1396         Double_t f34=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[3])/sy;
1397         Double_t f40=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[4])/sy;
1398         Double_t f41=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[4])/sz;
1399         Double_t f42=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[4])/sy;
1400         Double_t f43=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[4])/sz;    
1401         
1402         c[0]=sy1;
1403         //      c[1]=0.;       c[2]=sz1;
1404         c[1]=syz;       c[2]=sz1*100;
1405         c[3]=f20*sy1;  c[4]=0.;   c[5]=f20*sy1*f20+f22*sy2*f22+f24*sy3*f24;
1406         c[6]=f30*sy1;  c[7]=0.;   c[8]=f30*sy1*f20+f32*sy2*f22+f34*sy3*f24;
1407         c[9]=f30*sy1*f30+f32*sy2*f32+f34*sy3*f34;
1408         c[10]=f40*sy1; c[11]=f41*sz1;   c[12]=f40*sy1*f20+f42*sy2*f22;
1409         c[13]=f40*sy1*f30+f42*sy2*f32;
1410         c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;
1411         
1412         UInt_t index=r1.GetIndex(is);
1413         
1414         AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
1415
1416         Int_t rc=FollowProlongation(*track, i2);     
1417         
1418         if ((rc < 1) ||
1419             (track->GetNumberOfClusters() < 
1420              (outer-inner)*fMinClustersInSeed)) delete track;
1421         else {
1422           fSeeds->AddLast(track); fNseeds++;
1423           cerr<<"\r found seed "<<fNseeds;
1424         }
1425       }
1426     }
1427   }
1428 }            
1429
1430 //_____________________________________________________________________________
1431 void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp) 
1432 {
1433   //
1434   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) 
1435   // from the file. The names of the cluster tree and branches 
1436   // should match the ones used in AliTRDclusterizer::WriteClusters()
1437   //
1438
1439   TDirectory *savedir=gDirectory; 
1440
1441   if (inp) {
1442      TFile *in=(TFile*)inp;
1443      if (!in->IsOpen()) {
1444         cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n";
1445         return;
1446      }
1447      else{
1448        in->cd();
1449      }
1450   }
1451
1452   Char_t treeName[12];
1453   sprintf(treeName,"TreeR%d_TRD",fEvent);
1454   TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
1455   
1456   TObjArray *ClusterArray = new TObjArray(400); 
1457   
1458   ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); 
1459   
1460   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1461   printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
1462   
1463   // Loop through all entries in the tree
1464   Int_t nbytes;
1465   AliTRDcluster *c = 0;
1466   printf("\n");
1467
1468   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
1469     
1470     // Import the tree
1471     nbytes += ClusterTree->GetEvent(iEntry);  
1472     
1473     // Get the number of points in the detector
1474     Int_t nCluster = ClusterArray->GetEntriesFast();  
1475     printf("\r Read %d clusters from entry %d", nCluster, iEntry);
1476     
1477     // Loop through all TRD digits
1478     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
1479       c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
1480       AliTRDcluster *co = new AliTRDcluster(*c);
1481       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1482       Int_t ltb = co->GetLocalTimeBin();
1483       if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1484       
1485       array->AddLast(co);
1486       delete ClusterArray->RemoveAt(iCluster); 
1487     }
1488   }
1489
1490   delete ClusterArray;
1491   savedir->cd();   
1492
1493 }
1494
1495 //______________________________________________________________________
1496 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
1497 {
1498   //
1499   // Reads AliTRDclusters from file <filename>. The names of the cluster
1500   // tree and branches should match the ones used in
1501   // AliTRDclusterizer::WriteClusters()
1502   // if <array> == 0, clusters are added into AliTRDtracker fCluster array
1503   //
1504
1505   TDirectory *savedir=gDirectory;
1506
1507   TFile *file = TFile::Open(filename);
1508   if (!file->IsOpen()) {
1509     cerr<<"Can't open file with TRD clusters"<<endl;
1510     return;
1511   }
1512
1513   Char_t treeName[12];
1514   sprintf(treeName,"TreeR%d_TRD",fEvent);
1515   TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
1516
1517   if (!ClusterTree) {
1518      cerr<<"AliTRDtracker::ReadClusters(): ";
1519      cerr<<"can't get a tree with clusters !\n";
1520      return;
1521   }
1522
1523   TObjArray *ClusterArray = new TObjArray(400);
1524
1525   ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
1526
1527   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1528   cout<<"found "<<nEntries<<" in ClusterTree"<<endl;   
1529
1530   // Loop through all entries in the tree
1531   Int_t nbytes;
1532   AliTRDcluster *c = 0;
1533
1534   printf("\n");
1535
1536   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1537
1538     // Import the tree
1539     nbytes += ClusterTree->GetEvent(iEntry);
1540
1541     // Get the number of points in the detector
1542     Int_t nCluster = ClusterArray->GetEntriesFast();
1543     printf("\r Read %d clusters from entry %d", nCluster, iEntry);
1544
1545     // Loop through all TRD digits
1546     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
1547       c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
1548       AliTRDcluster *co = new AliTRDcluster(*c);
1549       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1550       Int_t ltb = co->GetLocalTimeBin();
1551       if(ltb != 0) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1552       array->AddLast(co);
1553       delete ClusterArray->RemoveAt(iCluster);
1554     }
1555   }
1556
1557   file->Close();
1558   delete ClusterArray;
1559   savedir->cd();
1560
1561 }                      
1562
1563
1564 //__________________________________________________________________
1565 void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const {
1566
1567   Int_t label=123456789, index, i, j;
1568   Int_t ncl=pt->GetNumberOfClusters();
1569   const Int_t range = fTrSec[0]->GetOuterTimeBin()+1;
1570
1571   Bool_t label_added;
1572
1573   //  Int_t s[range][2];
1574   Int_t **s = new Int_t* [range];
1575   for (i=0; i<range; i++) {
1576     s[i] = new Int_t[2];
1577   }
1578   for (i=0; i<range; i++) {
1579     s[i][0]=-1;
1580     s[i][1]=0;
1581   }
1582
1583   Int_t t0,t1,t2;
1584   for (i=0; i<ncl; i++) {
1585     index=pt->GetClusterIndex(i);
1586     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1587     t0=c->GetLabel(0);
1588     t1=c->GetLabel(1);
1589     t2=c->GetLabel(2);
1590   }
1591
1592   for (i=0; i<ncl; i++) {
1593     index=pt->GetClusterIndex(i);
1594     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1595     for (Int_t k=0; k<3; k++) { 
1596       label=c->GetLabel(k);
1597       label_added=kFALSE; j=0;
1598       if (label >= 0) {
1599         while ( (!label_added) && ( j < range ) ) {
1600           if (s[j][0]==label || s[j][1]==0) {
1601             s[j][0]=label; 
1602             s[j][1]=s[j][1]+1; 
1603             label_added=kTRUE;
1604           }
1605           j++;
1606         }
1607       }
1608     }
1609   }
1610
1611   Int_t max=0;
1612   label = -123456789;
1613
1614   for (i=0; i<range; i++) {
1615     if (s[i][1]>max) {
1616       max=s[i][1]; label=s[i][0];
1617     }
1618   }
1619
1620   for (i=0; i<range; i++) {
1621     delete []s[i];
1622   }        
1623
1624   delete []s;
1625
1626   if ((1.- Float_t(max)/ncl) > wrong) label=-label;   
1627
1628   pt->SetLabel(label); 
1629
1630 }
1631
1632
1633 //__________________________________________________________________
1634 void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const {
1635   Int_t ncl=t->GetNumberOfClusters();
1636   for (Int_t i=from; i<ncl; i++) {
1637     Int_t index = t->GetClusterIndex(i);
1638     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1639     c->Use();
1640   }
1641 }
1642
1643
1644 //_____________________________________________________________________
1645 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
1646 {
1647   // Parametrised "expected" error of the cluster reconstruction in Y 
1648
1649   Double_t s = 0.08 * 0.08;    
1650   return s;
1651 }
1652
1653 //_____________________________________________________________________
1654 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
1655 {
1656   // Parametrised "expected" error of the cluster reconstruction in Z 
1657
1658   Double_t s = 6 * 6 /12.;  
1659   return s;
1660 }                  
1661
1662
1663 //_____________________________________________________________________
1664 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t local_tb) const 
1665 {
1666   //
1667   // Returns radial position which corresponds to time bin <local_tb>
1668   // in tracking sector <sector> and plane <plane>
1669   //
1670
1671   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, local_tb); 
1672   Int_t pl = fTrSec[sector]->GetLayerNumber(index);
1673   return fTrSec[sector]->GetLayer(pl)->GetX();
1674
1675 }
1676
1677
1678 //_______________________________________________________
1679 AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x, 
1680                Double_t dx, Double_t rho, Double_t x0, Int_t tb_index)
1681
1682   //
1683   // AliTRDpropagationLayer constructor
1684   //
1685
1686   fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = x0;
1687   fClusters = NULL; fIndex = NULL; fTimeBinIndex = tb_index;
1688
1689
1690   for(Int_t i=0; i < (Int_t) kZONES; i++) {
1691     fZc[i]=0; fZmax[i] = 0;
1692   }
1693
1694   fYmax = 0;
1695
1696   if(fTimeBinIndex >= 0) { 
1697     fClusters = new AliTRDcluster*[kMAX_CLUSTER_PER_TIME_BIN];
1698     fIndex = new UInt_t[kMAX_CLUSTER_PER_TIME_BIN];
1699   }
1700
1701   fHole = kFALSE;
1702   fHoleZc = 0;
1703   fHoleZmax = 0;
1704   fHoleYc = 0;
1705   fHoleYmax = 0;
1706   fHoleRho = 0;
1707   fHoleX0 = 0;
1708
1709 }
1710
1711 //_______________________________________________________
1712 void AliTRDtracker::AliTRDpropagationLayer::SetHole(
1713           Double_t Zmax, Double_t Ymax, Double_t rho, 
1714           Double_t x0, Double_t Yc, Double_t Zc) 
1715 {
1716   //
1717   // Sets hole in the layer 
1718   //
1719
1720   fHole = kTRUE;
1721   fHoleZc = Zc;
1722   fHoleZmax = Zmax;
1723   fHoleYc = Yc;
1724   fHoleYmax = Ymax;
1725   fHoleRho = rho;
1726   fHoleX0 = x0;
1727 }
1728   
1729
1730 //_______________________________________________________
1731 AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
1732 {
1733   //
1734   // AliTRDtrackingSector Constructor
1735   //
1736
1737   fGeom = geo;
1738   fPar = par;
1739   fGeomSector = gs;
1740   fTzeroShift = 0.13;
1741   fN = 0;
1742
1743   for(UInt_t i=0; i < kMAX_TIME_BIN_INDEX; i++) fTimeBinIndex[i] = -1;
1744
1745
1746   AliTRDpropagationLayer* ppl;
1747
1748   Double_t x, xin, xout, dx, rho, x0;
1749   Int_t    steps;
1750
1751   // set time bins in the gas of the TPC
1752
1753   xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
1754   rho = 0.9e-3;  x0 = 28.94;
1755
1756   for(Int_t i=0; i<steps; i++) {
1757     x = xin + i*dx + dx/2;
1758     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1759     InsertLayer(ppl);
1760   }
1761
1762   // set time bins in the outer field cage vessel
1763
1764   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1765   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1766   InsertLayer(ppl);
1767
1768   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1769   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1770   InsertLayer(ppl);
1771
1772   dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; x0 = 41.28; // Nomex
1773   steps = 5; dx = (xout - xin)/steps;
1774   for(Int_t i=0; i<steps; i++) {
1775     x = xin + i*dx + dx/2;
1776     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1777     InsertLayer(ppl);
1778   }
1779
1780   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1781   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1782   InsertLayer(ppl);
1783
1784   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1785   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1786   InsertLayer(ppl);
1787
1788
1789   // set time bins in CO2
1790
1791   xin = xout; xout = 275.0; 
1792   steps = 50; dx = (xout - xin)/steps;
1793   rho = 1.977e-3;  x0 = 36.2;
1794   
1795   for(Int_t i=0; i<steps; i++) {
1796     x = xin + i*dx + dx/2;
1797     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1798     InsertLayer(ppl);
1799   }
1800
1801   // set time bins in the outer containment vessel
1802
1803   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; x0 = 24.01; // Al
1804   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1805   InsertLayer(ppl);
1806
1807   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1808   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1809   InsertLayer(ppl);
1810
1811   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1812   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1813   InsertLayer(ppl);
1814
1815   dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; x0 = 41.28; // Nomex
1816   steps = 10; dx = (xout - xin)/steps;
1817   for(Int_t i=0; i<steps; i++) {
1818     x = xin + i*dx + dx/2;
1819     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1820     InsertLayer(ppl);
1821   }
1822
1823   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; x0 = 44.86; // prepreg
1824   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1825   InsertLayer(ppl);
1826
1827   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; x0 = 44.77; // Tedlar
1828   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1829   InsertLayer(ppl);
1830   
1831   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; x0 = 24.01; // Al
1832   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1833   InsertLayer(ppl);
1834
1835   Double_t xtrd = (Double_t) fGeom->Rmin();  
1836
1837   // add layers between TPC and TRD (Air temporarily)
1838   xin = xout; xout = xtrd;
1839   steps = 50; dx = (xout - xin)/steps;
1840   rho = 1.2e-3;  x0 = 36.66;
1841   
1842   for(Int_t i=0; i<steps; i++) {
1843     x = xin + i*dx + dx/2;
1844     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1845     InsertLayer(ppl);
1846   }
1847
1848
1849   Double_t alpha=AliTRDgeometry::GetAlpha();
1850
1851   // add layers for each of the planes
1852
1853   Double_t dxRo = (Double_t) fGeom->CroHght();    // Rohacell 
1854   Double_t dxSpace = (Double_t) fGeom->Cspace();  // Spacing between planes
1855   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
1856   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region  
1857   Double_t dxRad = (Double_t) fGeom->CraHght();   // Radiator
1858   Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo; 
1859   Double_t dxPlane = dxTEC + dxSpace; 
1860
1861   Int_t tbBefore = (Int_t) (dxAmp/fPar->GetTimeBinSize());
1862   if(tbBefore > fPar->GetTimeBefore()) tbBefore = fPar->GetTimeBefore();
1863
1864   Int_t tb, tb_index;
1865   const Int_t  nChambers = AliTRDgeometry::Ncham();
1866   Double_t  Ymax = 0, holeYmax = 0;
1867   Double_t *  Zc  = new Double_t[nChambers];
1868   Double_t * Zmax = new Double_t[nChambers];
1869   Double_t  holeZmax = 1000.;   // the whole sector is missing
1870
1871
1872   for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
1873
1874     // Radiator 
1875     xin = xtrd + plane * dxPlane; xout = xin + dxRad;
1876     steps = 12; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6; 
1877     for(Int_t i=0; i<steps; i++) {
1878       x = xin + i*dx + dx/2;
1879       ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
1880       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1881         holeYmax = x*TMath::Tan(0.5*alpha);
1882         ppl->SetHole(holeYmax, holeZmax);
1883       }
1884       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1885         holeYmax = x*TMath::Tan(0.5*alpha);
1886         ppl->SetHole(holeYmax, holeZmax);
1887       }
1888       InsertLayer(ppl);
1889     }
1890
1891     Ymax = fGeom->GetChamberWidth(plane)/2;
1892     for(Int_t ch = 0; ch < nChambers; ch++) {
1893       Zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
1894       Float_t pad = fPar->GetRowPadSize(plane,ch,0);
1895       Float_t row0 = fPar->GetRow0(plane,ch,0);
1896       Int_t nPads = fPar->GetRowMax(plane,ch,0);
1897       Zc[ch] = (pad * nPads)/2 + row0 - pad/2;
1898     }
1899
1900     dx = fPar->GetTimeBinSize(); 
1901     rho = 0.00295 * 0.85; x0 = 11.0;  
1902
1903     Double_t x0 = (Double_t) fPar->GetTime0(plane);
1904     Double_t xbottom = x0 - dxDrift;
1905     Double_t xtop = x0 + dxAmp;
1906
1907     // Amplification region
1908
1909     steps = (Int_t) (dxAmp/dx);
1910
1911     for(tb = 0; tb < steps; tb++) {
1912       x = x0 + tb * dx + dx/2;
1913       tb_index = CookTimeBinIndex(plane, -tb-1);
1914       ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1915       ppl->SetYmax(Ymax);
1916       for(Int_t ch = 0; ch < nChambers; ch++) {
1917         ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1918       }
1919       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1920         holeYmax = x*TMath::Tan(0.5*alpha);
1921         ppl->SetHole(holeYmax, holeZmax);
1922       }
1923       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1924         holeYmax = x*TMath::Tan(0.5*alpha);
1925         ppl->SetHole(holeYmax, holeZmax);
1926       }
1927       InsertLayer(ppl);
1928     }
1929     tb_index = CookTimeBinIndex(plane, -steps);
1930     x = (x + dx/2 + xtop)/2;
1931     dx = 2*(xtop-x);
1932     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1933     ppl->SetYmax(Ymax);
1934     for(Int_t ch = 0; ch < nChambers; ch++) {
1935       ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1936     }
1937     if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1938       holeYmax = x*TMath::Tan(0.5*alpha);
1939       ppl->SetHole(holeYmax, holeZmax);
1940     }
1941     if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1942       holeYmax = x*TMath::Tan(0.5*alpha);
1943       ppl->SetHole(holeYmax, holeZmax);
1944     }
1945     InsertLayer(ppl);
1946
1947     // Drift region
1948     dx = fPar->GetTimeBinSize();
1949     steps = (Int_t) (dxDrift/dx);
1950
1951     for(tb = 0; tb < steps; tb++) {
1952       x = x0 - tb * dx - dx/2;
1953       tb_index = CookTimeBinIndex(plane, tb);
1954
1955       ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1956       ppl->SetYmax(Ymax);
1957       for(Int_t ch = 0; ch < nChambers; ch++) {
1958         ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1959       }
1960       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1961         holeYmax = x*TMath::Tan(0.5*alpha);
1962         ppl->SetHole(holeYmax, holeZmax);
1963       }
1964       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1965         holeYmax = x*TMath::Tan(0.5*alpha);
1966         ppl->SetHole(holeYmax, holeZmax);
1967       }
1968       InsertLayer(ppl);
1969     }
1970     tb_index = CookTimeBinIndex(plane, steps);
1971     x = (x - dx/2 + xbottom)/2;
1972     dx = 2*(x-xbottom);
1973     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,tb_index);
1974     ppl->SetYmax(Ymax);
1975     for(Int_t ch = 0; ch < nChambers; ch++) {
1976       ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
1977     }
1978     if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1979       holeYmax = x*TMath::Tan(0.5*alpha);
1980       ppl->SetHole(holeYmax, holeZmax);
1981     }
1982     if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1983       holeYmax = x*TMath::Tan(0.5*alpha);
1984       ppl->SetHole(holeYmax, holeZmax);
1985     }
1986     InsertLayer(ppl);
1987
1988     // Pad Plane
1989     xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; x0 = 33.0;
1990     ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,x0,-1);
1991     if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
1992       holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
1993       ppl->SetHole(holeYmax, holeZmax);
1994     }
1995     if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
1996       holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
1997       ppl->SetHole(holeYmax, holeZmax);
1998     }
1999     InsertLayer(ppl);
2000
2001     // Rohacell
2002     xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
2003     steps = 5; dx = (xout - xin)/steps; rho = 0.074; x0 = 40.6; 
2004     for(Int_t i=0; i<steps; i++) {
2005       x = xin + i*dx + dx/2;
2006       ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
2007       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2008         holeYmax = x*TMath::Tan(0.5*alpha);
2009         ppl->SetHole(holeYmax, holeZmax);
2010       }
2011       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2012         holeYmax = x*TMath::Tan(0.5*alpha);
2013         ppl->SetHole(holeYmax, holeZmax);
2014       }
2015       InsertLayer(ppl);
2016     }
2017
2018     // Space between the chambers, air
2019     xin = xout; xout = xtrd + (plane + 1) * dxPlane;
2020     steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66; 
2021     for(Int_t i=0; i<steps; i++) {
2022       x = xin + i*dx + dx/2;
2023       ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
2024       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2025         holeYmax = x*TMath::Tan(0.5*alpha);
2026         ppl->SetHole(holeYmax, holeZmax);
2027       }
2028       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2029         holeYmax = x*TMath::Tan(0.5*alpha);
2030         ppl->SetHole(holeYmax, holeZmax);
2031       }
2032       InsertLayer(ppl);
2033     }
2034   }    
2035
2036   // Space between the TRD and RICH
2037   Double_t xRICH = 500.;
2038   xin = xout; xout = xRICH;
2039   steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; x0 = 36.66; 
2040   for(Int_t i=0; i<steps; i++) {
2041     x = xin + i*dx + dx/2;
2042     ppl = new AliTRDpropagationLayer(x,dx,rho,x0,-1);
2043     InsertLayer(ppl);
2044   }
2045
2046   MapTimeBinLayers();
2047   delete [] Zc;
2048   delete [] Zmax;
2049
2050 }
2051
2052 //______________________________________________________
2053
2054 Int_t  AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t local_tb) const
2055 {
2056   //
2057   // depending on the digitization parameters calculates "global"
2058   // time bin index for timebin <local_tb> in plane <plane>
2059   //
2060
2061   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2062   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region  
2063   Double_t dx = (Double_t) fPar->GetTimeBinSize();  
2064
2065   Int_t tbAmp = fPar->GetTimeBefore();
2066   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
2067   Int_t tbDrift = fPar->GetTimeMax();
2068   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
2069
2070   Int_t tb_per_plane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
2071
2072   Int_t gtb = (plane+1) * tb_per_plane - local_tb - 1;
2073
2074   if((local_tb < 0) && 
2075      (TMath::Abs(local_tb) > TMath::Min(tbAmp,maxAmp))) return -1;
2076   if(local_tb >= TMath::Min(tbDrift,maxDrift)) return -1;
2077
2078   return gtb;
2079
2080
2081 }
2082
2083 //______________________________________________________
2084
2085 void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers() 
2086 {
2087   //
2088   // For all sensitive time bins sets corresponding layer index
2089   // in the array fTimeBins 
2090   //
2091
2092   Int_t index;
2093
2094   for(Int_t i = 0; i < fN; i++) {
2095     index = fLayers[i]->GetTimeBinIndex();
2096     
2097     //    printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2098
2099     if(index < 0) continue;
2100     if(index >= (Int_t) kMAX_TIME_BIN_INDEX) {
2101       printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2102       printf("    index %d exceeds allowed maximum of %d!\n",
2103              index, kMAX_TIME_BIN_INDEX-1);
2104       continue;
2105     }
2106     fTimeBinIndex[index] = i;
2107   }
2108
2109   Double_t x1, dx1, x2, dx2, gap;
2110
2111   for(Int_t i = 0; i < fN-1; i++) {
2112     x1 = fLayers[i]->GetX();
2113     dx1 = fLayers[i]->GetdX();
2114     x2 = fLayers[i+1]->GetX();
2115     dx2 = fLayers[i+1]->GetdX();
2116     gap = (x2 - dx2/2) - (x1 + dx1/2);
2117     if(gap < -0.01) {
2118       printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2119       printf("             %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2120     }
2121     if(gap > 0.01) { 
2122       printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2123       printf("             (%f - %f) - (%f + %f) = %f\n", 
2124              x2, dx2/2, x1, dx1, gap);
2125     }
2126   }
2127 }
2128   
2129
2130 //______________________________________________________
2131
2132
2133 Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2134 {
2135   // 
2136   // Returns the number of time bin which in radial position is closest to <x>
2137   //
2138
2139   if(x >= fLayers[fN-1]->GetX()) return fN-1; 
2140   if(x <= fLayers[0]->GetX()) return 0; 
2141
2142   Int_t b=0, e=fN-1, m=(b+e)/2;
2143   for (; b<e; m=(b+e)/2) {
2144     if (x > fLayers[m]->GetX()) b=m+1;
2145     else e=m;
2146   }
2147   if(TMath::Abs(x - fLayers[m]->GetX()) > 
2148      TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2149   else return m;
2150
2151 }
2152
2153 //______________________________________________________
2154
2155 Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const 
2156 {
2157   // 
2158   // Returns number of the innermost SENSITIVE propagation layer
2159   //
2160
2161   return GetLayerNumber(0);
2162 }
2163
2164 //______________________________________________________
2165
2166 Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const 
2167 {
2168   // 
2169   // Returns number of the outermost SENSITIVE time bin
2170   //
2171
2172   return GetLayerNumber(GetNumberOfTimeBins() - 1);
2173 }
2174
2175 //______________________________________________________
2176
2177 Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const 
2178 {
2179   // 
2180   // Returns number of SENSITIVE time bins
2181   //
2182
2183   Int_t tb, layer;
2184   for(tb = kMAX_TIME_BIN_INDEX-1; tb >=0; tb--) {
2185     layer = GetLayerNumber(tb);
2186     if(layer>=0) break;
2187   }
2188   return tb+1;
2189 }
2190
2191 //______________________________________________________
2192
2193 void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2194
2195   //
2196   // Insert layer <pl> in fLayers array.
2197   // Layers are sorted according to X coordinate.
2198
2199   if ( fN == ((Int_t) kMAX_LAYERS_PER_SECTOR)) {
2200     printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2201     return;
2202   }
2203   if (fN==0) {fLayers[fN++] = pl; return;}
2204   Int_t i=Find(pl->GetX());
2205
2206   memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2207   fLayers[i]=pl; fN++;
2208
2209 }              
2210
2211 //______________________________________________________
2212
2213 Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const 
2214 {
2215   //
2216   // Returns index of the propagation layer nearest to X 
2217   //
2218
2219   if (x <= fLayers[0]->GetX()) return 0;
2220   if (x > fLayers[fN-1]->GetX()) return fN;
2221   Int_t b=0, e=fN-1, m=(b+e)/2;
2222   for (; b<e; m=(b+e)/2) {
2223     if (x > fLayers[m]->GetX()) b=m+1;
2224     else e=m;
2225   }
2226   return m;
2227 }             
2228
2229 //______________________________________________________
2230
2231 void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
2232         Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &x0, 
2233         Bool_t &lookForCluster) const
2234 {
2235   //
2236   // Returns radial step <dx>, density <rho>, rad. length <x0>,
2237   // and sensitivity <lookForCluster> in point <y,z>  
2238   //
2239
2240   dx  = fdX;
2241   rho = fRho;
2242   x0  = fX0;
2243   lookForCluster = kFALSE;
2244
2245   // check dead regions
2246   if(fTimeBinIndex >= 0) {
2247     for(Int_t ch = 0; ch < (Int_t) kZONES; ch++) {
2248       if(TMath::Abs(z - fZc[ch]) < fZmax[ch]) 
2249         lookForCluster = kTRUE;
2250     }
2251     if(TMath::Abs(y) > fYmax) lookForCluster = kFALSE;
2252     if(!lookForCluster) { 
2253       //      rho = 1.7; x0 = 33.0; // G10 
2254     }
2255   }
2256
2257   // check hole
2258   if(fHole && (TMath::Abs(y - fHoleYc) < fHoleYmax) && 
2259               (TMath::Abs(z - fHoleZc) < fHoleZmax)) {
2260     lookForCluster = kFALSE;
2261     rho = fHoleRho;
2262     x0  = fHoleX0;
2263   }         
2264
2265   return;
2266 }
2267
2268 //______________________________________________________
2269
2270 void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c, 
2271                                                           UInt_t index) {
2272
2273 // Insert cluster in cluster array.
2274 // Clusters are sorted according to Y coordinate.  
2275
2276   if(fTimeBinIndex < 0) { 
2277     printf("*** attempt to insert cluster into non-sensitive time bin!\n");
2278     return;
2279   }
2280
2281   if (fN== (Int_t) kMAX_CLUSTER_PER_TIME_BIN) {
2282     printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n"); 
2283     return;
2284   }
2285   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2286   Int_t i=Find(c->GetY());
2287   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
2288   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t)); 
2289   fIndex[i]=index; fClusters[i]=c; fN++;
2290 }  
2291
2292 //______________________________________________________
2293
2294 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
2295
2296 // Returns index of the cluster nearest in Y    
2297
2298   if (y <= fClusters[0]->GetY()) return 0;
2299   if (y > fClusters[fN-1]->GetY()) return fN;
2300   Int_t b=0, e=fN-1, m=(b+e)/2;
2301   for (; b<e; m=(b+e)/2) {
2302     if (y > fClusters[m]->GetY()) b=m+1;
2303     else e=m;
2304   }
2305   return m;
2306 }    
2307
2308 //---------------------------------------------------------
2309
2310 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
2311 //
2312 //  Returns correction factor for tilted pads geometry 
2313 //
2314
2315   Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
2316   Int_t det = c->GetDetector();    
2317   Int_t plane = fGeom->GetPlane(det);
2318   if((plane == 0) || (plane == 2) || (plane == 4)) h01=-h01;
2319   
2320   return h01;
2321 }
2322
2323
2324
2325
2326
2327