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