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