fd39eb1876d1ae50fbb634acfb6c34973c02aaeb
[u/mrichter/AliRoot.git] / TPC / AliTPCtracker.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.10  2001/05/16 14:57:25  alibrary
19 New files for folders and Stack
20
21 Revision 1.9  2001/05/11 07:16:56  hristov
22 Fix needed on Sun and Alpha
23
24 Revision 1.8  2001/05/08 15:00:15  hristov
25 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
26
27 Revision 1.5  2000/12/20 07:51:59  kowal2
28 Changes suggested by Alessandra and Paolo to avoid overlapped
29 data fields in encapsulated classes.
30
31 Revision 1.4  2000/11/02 07:27:16  kowal2
32 code corrections
33
34 Revision 1.2  2000/06/30 12:07:50  kowal2
35 Updated from the TPC-PreRelease branch
36
37 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
38 Splitted from AliTPCtracking
39
40 */
41
42 //-------------------------------------------------------
43 //          Implementation of the TPC tracker
44 //
45 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
46 //-------------------------------------------------------
47
48 #include <TObjArray.h>
49 #include <TFile.h>
50 #include <TTree.h>
51 #include <iostream.h>
52
53 #include "AliTPCtracker.h"
54 #include "AliTPCcluster.h"
55 #include "AliTPCParam.h"
56 #include "AliTPCClustersRow.h"
57
58 //_____________________________________________________________________________
59 AliTPCtracker::AliTPCtracker(const AliTPCParam *par): 
60 fkNIS(par->GetNInnerSector()/2), 
61 fkNOS(par->GetNOuterSector()/2)
62 {
63   //---------------------------------------------------------------------
64   // The main TPC tracker constructor
65   //---------------------------------------------------------------------
66   fInnerSec=new AliTPCSector[fkNIS];         
67   fOuterSec=new AliTPCSector[fkNOS];
68
69   Int_t i;
70   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
71   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
72
73   fN=0;  fSectors=0;
74
75   fClustersArray.Setup(par);
76   fClustersArray.SetClusterType("AliTPCcluster");
77   fClustersArray.ConnectTree("Segment Tree");
78
79   fSeeds=0;
80 }
81
82 //_____________________________________________________________________________
83 AliTPCtracker::~AliTPCtracker() {
84   //------------------------------------------------------------------
85   // TPC tracker destructor
86   //------------------------------------------------------------------
87   delete[] fInnerSec;
88   delete[] fOuterSec;
89   delete fSeeds;
90 }
91
92 //_____________________________________________________________________________
93 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
94 {
95   //
96   // Parametrised error of the cluster reconstruction (pad direction)   
97   //
98   // Sigma rphi
99   const Float_t kArphi=0.41818e-2;
100   const Float_t kBrphi=0.17460e-4;
101   const Float_t kCrphi=0.30993e-2;
102   const Float_t kDrphi=0.41061e-3;
103   
104   pt=TMath::Abs(pt)*1000.;
105   Double_t x=r/pt;
106   tgl=TMath::Abs(tgl);
107   Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
108   if (s<0.4e-3) s=0.4e-3;
109   s*=1.3; //Iouri Belikov
110
111   return s;
112 }
113
114 //_____________________________________________________________________________
115 Double_t SigmaZ2(Double_t r, Double_t tgl) 
116 {
117   //
118   // Parametrised error of the cluster reconstruction (drift direction)
119   //
120   // Sigma z
121   const Float_t kAz=0.39614e-2;
122   const Float_t kBz=0.22443e-4;
123   const Float_t kCz=0.51504e-1;
124   
125
126   tgl=TMath::Abs(tgl);
127   Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
128   if (s<0.4e-3) s=0.4e-3;
129   s*=1.3; //Iouri Belikov
130
131   return s;
132 }
133
134 //_____________________________________________________________________________
135 Double_t f1(Double_t x1,Double_t y1,
136                    Double_t x2,Double_t y2,
137                    Double_t x3,Double_t y3) 
138 {
139   //-----------------------------------------------------------------
140   // Initial approximation of the track curvature
141   //-----------------------------------------------------------------
142   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
143   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
144                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
145   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
146                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
147
148   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
149
150   return -xr*yr/sqrt(xr*xr+yr*yr); 
151 }
152
153
154 //_____________________________________________________________________________
155 Double_t f2(Double_t x1,Double_t y1,
156                    Double_t x2,Double_t y2,
157                    Double_t x3,Double_t y3) 
158 {
159   //-----------------------------------------------------------------
160   // Initial approximation of the track curvature times center of curvature
161   //-----------------------------------------------------------------
162   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
163   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
164                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
165   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
166                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
167
168   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
169   
170   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
171 }
172
173 //_____________________________________________________________________________
174 Double_t f3(Double_t x1,Double_t y1, 
175                    Double_t x2,Double_t y2,
176                    Double_t z1,Double_t z2) 
177 {
178   //-----------------------------------------------------------------
179   // Initial approximation of the tangent of the track dip angle
180   //-----------------------------------------------------------------
181   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
182 }
183
184 //_____________________________________________________________________________
185 void AliTPCtracker::LoadOuterSectors() {
186   //-----------------------------------------------------------------
187   // This function fills outer TPC sectors with clusters.
188   //-----------------------------------------------------------------
189   UInt_t index;
190   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
191   for (Int_t i=0; i<j; i++) {
192       AliSegmentID *s=fClustersArray.LoadEntry(i);
193       Int_t sec,row;
194       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
195       par->AdjustSectorRow(s->GetID(),sec,row);
196       if (sec<fkNIS*2) continue;
197       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
198       Int_t ncl=clrow->GetArray()->GetEntriesFast();
199       while (ncl--) {
200          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
201          index=(((sec<<8)+row)<<16)+ncl;
202          fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
203       }
204   }
205
206   fN=fkNOS;
207   fSectors=fOuterSec;
208 }
209
210 //_____________________________________________________________________________
211 void AliTPCtracker::UnloadOuterSectors() {
212   //-----------------------------------------------------------------
213   // This function clears outer TPC sectors.
214   //-----------------------------------------------------------------
215   Int_t nup=fOuterSec->GetNRows();
216   for (Int_t i=0; i<fkNOS; i++) {
217     for (Int_t j=0; j<nup; j++) {
218       if (fClustersArray.GetRow(i+fkNIS*2,j)) 
219           fClustersArray.ClearRow(i+fkNIS*2,j);
220       if (fClustersArray.GetRow(i+fkNIS*2+fkNOS,j)) 
221           fClustersArray.ClearRow(i+fkNIS*2+fkNOS,j);
222     }
223   }
224
225   fN=0;
226   fSectors=0;
227 }
228
229 //_____________________________________________________________________________
230 void AliTPCtracker::LoadInnerSectors() {
231   //-----------------------------------------------------------------
232   // This function fills inner TPC sectors with clusters.
233   //-----------------------------------------------------------------
234   UInt_t index;
235   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
236   for (Int_t i=0; i<j; i++) {
237       AliSegmentID *s=fClustersArray.LoadEntry(i);
238       Int_t sec,row;
239       AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
240       par->AdjustSectorRow(s->GetID(),sec,row);
241       if (sec>=fkNIS*2) continue;
242       AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
243       Int_t ncl=clrow->GetArray()->GetEntriesFast();
244       while (ncl--) {
245          AliTPCcluster *c=(AliTPCcluster*)(*clrow)[ncl];
246          index=(((sec<<8)+row)<<16)+ncl;
247          fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
248       }
249   }
250
251   fN=fkNIS;
252   fSectors=fInnerSec;
253 }
254
255 //_____________________________________________________________________________
256 void AliTPCtracker::UnloadInnerSectors() {
257   //-----------------------------------------------------------------
258   // This function clears inner TPC sectors.
259   //-----------------------------------------------------------------
260   Int_t nlow=fInnerSec->GetNRows();
261   for (Int_t i=0; i<fkNIS; i++) {
262     for (Int_t j=0; j<nlow; j++) {
263       if (fClustersArray.GetRow(i,j)) fClustersArray.ClearRow(i,j);
264       if (fClustersArray.GetRow(i+fkNIS,j)) fClustersArray.ClearRow(i+fkNIS,j);
265     }
266   }
267
268   fN=0;
269   fSectors=0;
270 }
271
272 //_____________________________________________________________________________
273 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
274   //-----------------------------------------------------------------
275   // This function tries to find a track prolongation.
276   //-----------------------------------------------------------------
277   Double_t xt=t.GetX();
278   const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip : 
279                     Int_t(0.5*fSectors->GetNRows());
280   Int_t tryAgain=kSKIP;
281
282   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
283   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
284   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
285   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
286
287   for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
288     Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
289     if (!t.PropagateTo(x)) return 0;
290
291     AliTPCcluster *cl=0;
292     UInt_t index=0;
293     Double_t maxchi2=kMaxCHI2;
294     const AliTPCRow &krow=fSectors[s][nr];
295     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
296     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
297     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
298     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
299
300     if (road>kMaxROAD) {
301       if (t.GetNumberOfClusters()>4) 
302         cerr<<t.GetNumberOfClusters()
303         <<"FindProlongation warning: Too broad road !\n"; 
304       return 0;
305     }
306
307     if (krow) {
308       for (Int_t i=krow.Find(y-road); i<krow; i++) {
309         AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
310         if (c->GetY() > y+road) break;
311         if (c->IsUsed()) continue;
312         if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
313         Double_t chi2=t.GetPredictedChi2(c);
314         if (chi2 > maxchi2) continue;
315         maxchi2=chi2;
316         cl=c;
317         index=krow.GetIndex(i);       
318       }
319     }
320     if (cl) {
321       Float_t l=fSectors->GetPadPitchWidth();
322       t.SetSampledEdx(cl->GetQ()/l,t.GetNumberOfClusters());
323       if (!t.Update(cl,maxchi2,index)) {
324          if (!tryAgain--) return 0;
325       } else tryAgain=kSKIP;
326     } else {
327       if (tryAgain==0) break;
328       if (y > ymax) {
329          s = (s+1) % fN;
330          if (!t.Rotate(fSectors->GetAlpha())) return 0;
331       } else if (y <-ymax) {
332          s = (s-1+fN) % fN;
333          if (!t.Rotate(-fSectors->GetAlpha())) return 0;
334       }
335       tryAgain--;
336     }
337   }
338
339   return 1;
340 }
341
342 //_____________________________________________________________________________
343 Int_t AliTPCtracker::FollowBackProlongation
344 (AliTPCseed& seed, const AliTPCtrack &track) {
345   //-----------------------------------------------------------------
346   // This function propagates tracks back through the TPC
347   //-----------------------------------------------------------------
348   Double_t alpha=seed.GetAlpha() - fSectors->GetAlphaShift();
349   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
350   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
351   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
352
353   Int_t idx=-1, sec=-1, row=-1;
354   Int_t nc=seed.GetLabel(); //index of the cluster to start with
355   if (nc--) {
356      idx=track.GetClusterIndex(nc);
357      sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
358   }
359   if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; } 
360   else { if (sec <  2*fkNIS) row=-1; }   
361
362   Int_t nr=fSectors->GetNRows();
363   for (Int_t i=0; i<nr; i++) {
364     Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
365
366     if (!seed.PropagateTo(x)) return 0;
367
368     Double_t y=seed.GetY();
369     if (y > ymax) {
370        s = (s+1) % fN;
371        if (!seed.Rotate(fSectors->GetAlpha())) return 0;
372     } else if (y <-ymax) {
373        s = (s-1+fN) % fN;
374        if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
375     }
376
377     AliTPCcluster *cl=0;
378     Int_t index=0;
379     Double_t maxchi2=kMaxCHI2;
380     Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
381     Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
382     Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
383     Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
384     if (road>kMaxROAD) {
385       cerr<<seed.GetNumberOfClusters()
386           <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n"; 
387       return 0;
388     }
389
390
391     Int_t accepted=seed.GetNumberOfClusters();
392     if (row==i) {
393        //try to accept already found cluster
394        AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
395        Double_t chi2;
396        if ((chi2=seed.GetPredictedChi2(c))<maxchi2 || accepted<27) { 
397           index=idx; cl=c; maxchi2=chi2;
398        } //else cerr<<"AliTPCtracker::FollowBackProlongation: oulier !\n";
399           
400        if (nc--) {
401           idx=track.GetClusterIndex(nc); 
402           sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
403        } 
404        if (fSectors==fInnerSec) { if (sec >= 2*fkNIS) row=-1; }
405        else { if (sec <  2*fkNIS) row=-1; }   
406
407     }
408     if (!cl) {
409        //try to fill the gap
410        const AliTPCRow &krow=fSectors[s][i];
411        if (accepted>27)
412        if (krow) {
413           for (Int_t i=krow.Find(y-road); i<krow; i++) {
414             AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
415             if (c->GetY() > y+road) break;
416             if (c->IsUsed()) continue;
417          if ((c->GetZ()-z)*(c->GetZ()-z)>16.*(seed.GetSigmaZ2()+sz2)) continue;
418             Double_t chi2=seed.GetPredictedChi2(c);
419             if (chi2 > maxchi2) continue;
420             maxchi2=chi2;
421             cl=c;
422             index=krow.GetIndex(i);
423           }
424        }
425     }
426
427     if (cl) {
428       Float_t l=fSectors->GetPadPitchWidth();
429       seed.SetSampledEdx(cl->GetQ()/l,seed.GetNumberOfClusters());
430       seed.Update(cl,maxchi2,index);
431     }
432
433   }
434
435   seed.SetLabel(nc);
436
437   return 1;
438 }
439
440 //_____________________________________________________________________________
441 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
442   //-----------------------------------------------------------------
443   // This function creates track seeds.
444   //-----------------------------------------------------------------
445   if (fSeeds==0) fSeeds=new TObjArray(15000);
446
447   Double_t x[5], c[15];
448
449   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
450   Double_t cs=cos(alpha), sn=sin(alpha);
451
452   Double_t x1 =fOuterSec->GetX(i1);
453   Double_t xx2=fOuterSec->GetX(i2);
454
455   for (Int_t ns=0; ns<fkNOS; ns++) {
456     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
457     Int_t nm=fOuterSec[ns][i2];
458     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
459     const AliTPCRow& kr1=fOuterSec[ns][i1];
460     for (Int_t is=0; is < kr1; is++) {
461       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
462       for (Int_t js=0; js < nl+nm+nu; js++) {
463         const AliTPCcluster *kcl;
464         Double_t x2,   y2,   z2;
465         Double_t x3=0.,y3=0.;
466
467         if (js<nl) {
468           const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
469           kcl=kr2[js];
470           y2=kcl->GetY(); z2=kcl->GetZ();
471           x2= xx2*cs+y2*sn;
472           y2=-xx2*sn+y2*cs;
473         } else 
474           if (js<nl+nm) {
475             const AliTPCRow& kr2=fOuterSec[ns][i2];
476             kcl=kr2[js-nl];
477             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
478           } else {
479             const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
480             kcl=kr2[js-nl-nm];
481             y2=kcl->GetY(); z2=kcl->GetZ();
482             x2=xx2*cs-y2*sn;
483             y2=xx2*sn+y2*cs;
484           }
485
486         Double_t zz=z1 - z1/x1*(x1-x2); 
487         if (TMath::Abs(zz-z2)>5.) continue;
488
489         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
490         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
491
492         x[0]=y1;
493         x[1]=z1;
494         x[4]=f1(x1,y1,x2,y2,x3,y3);
495         if (TMath::Abs(x[4]) >= 0.0066) continue;
496         x[2]=f2(x1,y1,x2,y2,x3,y3);
497         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
498         x[3]=f3(x1,y1,x2,y2,z1,z2);
499         if (TMath::Abs(x[3]) > 1.2) continue;
500         Double_t a=asin(x[2]);
501         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
502         if (TMath::Abs(zv)>10.) continue; 
503
504         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
505         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
506         Double_t sy3=100*0.025, sy=0.1, sz=0.1;
507
508         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
509         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
510         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
511         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
512         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
513         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
514         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
515         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
516         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
517         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
518
519         c[0]=sy1;
520         c[1]=0.;       c[2]=sz1;
521         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
522         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
523                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
524         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
525         c[13]=f30*sy1*f40+f32*sy2*f42;
526         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
527
528         UInt_t index=kr1.GetIndex(is);
529         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
530         Float_t l=fOuterSec->GetPadPitchWidth();
531         track->SetSampledEdx(kr1[is]->GetQ()/l,0);
532
533         Int_t rc=FollowProlongation(*track, i2);
534         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/2) delete track;
535         else fSeeds->AddLast(track); 
536       }
537     }
538   }
539 }
540
541 //_____________________________________________________________________________
542 Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
543   //-----------------------------------------------------------------
544   // This function reades track seeds.
545   //-----------------------------------------------------------------
546   TDirectory *savedir=gDirectory; 
547
548   TFile *in=(TFile*)inp;
549   if (!in->IsOpen()) {
550      cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
551      return 1;
552   }
553
554   in->cd();
555   TTree *seedTree=(TTree*)in->Get("Seeds");
556   if (!seedTree) {
557      cerr<<"AliTPCtracker::ReadSeeds(): ";
558      cerr<<"can't get a tree with track seeds !\n";
559      return 2;
560   }
561   AliTPCtrack *seed=new AliTPCtrack; 
562   seedTree->SetBranchAddress("tracks",&seed);
563   
564   if (fSeeds==0) fSeeds=new TObjArray(15000);
565
566   Int_t n=(Int_t)seedTree->GetEntries();
567   for (Int_t i=0; i<n; i++) {
568      seedTree->GetEvent(i);
569      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
570   }
571   
572   delete seed;
573   savedir->cd();
574
575   return 0;
576 }
577
578 //_____________________________________________________________________________
579 Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
580   //-----------------------------------------------------------------
581   // This is a track finder.
582   //-----------------------------------------------------------------
583   TDirectory *savedir=gDirectory; 
584
585   if (inp) {
586      TFile *in=(TFile*)inp;
587      if (!in->IsOpen()) {
588         cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
589         return 1;
590      }
591   }
592
593   if (!out->IsOpen()) {
594      cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
595      return 2;
596   }
597
598   out->cd();
599   TTree tracktree("TPCf","Tree with TPC tracks");
600   AliTPCtrack *iotrack=0;
601   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
602
603   LoadOuterSectors();
604
605   //find track seeds
606   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
607   Int_t nrows=nlow+nup;
608   if (fSeeds==0) {
609      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
610      MakeSeeds(nup-1, nup-1-gap);
611      MakeSeeds(nup-1-shift, nup-1-shift-gap);
612   }    
613   fSeeds->Sort();
614
615   //tracking in outer sectors
616   Int_t nseed=fSeeds->GetEntriesFast();
617   Int_t i;
618   for (i=0; i<nseed; i++) {
619     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
620     if (FollowProlongation(t)) {
621        UseClusters(&t);
622        continue;
623     }
624     delete fSeeds->RemoveAt(i);
625   }  
626   UnloadOuterSectors();
627
628   //tracking in inner sectors
629   LoadInnerSectors();
630   Int_t found=0;
631   for (i=0; i<nseed; i++) {
632     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
633     if (!pt) continue;
634     Int_t nc=t.GetNumberOfClusters();
635
636     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
637     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
638     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
639     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
640
641     alpha=ns*fInnerSec->GetAlpha() + fInnerSec->GetAlphaShift() - t.GetAlpha();
642
643     if (t.Rotate(alpha)) {
644        if (FollowProlongation(t)) {
645           if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
646              t.CookdEdx();
647              iotrack=pt;
648              tracktree.Fill();
649              UseClusters(&t,nc);
650              cerr<<found++<<'\r';
651           }
652        }
653     }
654     delete fSeeds->RemoveAt(i); 
655   }  
656   UnloadInnerSectors();
657
658   tracktree.Write();
659
660   cerr<<"Number of found tracks : "<<found<<endl;
661
662   savedir->cd();
663
664   return 0;
665 }
666
667 //_____________________________________________________________________________
668 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
669   //-----------------------------------------------------------------
670   // This function propagates tracks back through the TPC.
671   //-----------------------------------------------------------------
672   fSeeds=new TObjArray(15000);
673
674   TFile *in=(TFile*)inp;
675   TDirectory *savedir=gDirectory; 
676
677   if (!in->IsOpen()) {
678      cerr<<"AliTPCtracker::PropagateBack(): ";
679      cerr<<"file with back propagated ITS tracks is not open !\n";
680      return 1;
681   }
682
683   if (!out->IsOpen()) {
684      cerr<<"AliTPCtracker::PropagateBack(): ";
685      cerr<<"file for back propagated TPC tracks is not open !\n";
686      return 2;
687   }
688
689   in->cd();
690   TTree *bckTree=(TTree*)in->Get("ITSb");
691   if (!bckTree) {
692      cerr<<"AliTPCtracker::PropagateBack() ";
693      cerr<<"can't get a tree with back propagated ITS tracks !\n";
694      return 3;
695   }
696   AliTPCtrack *bckTrack=new AliTPCtrack; 
697   bckTree->SetBranchAddress("tracks",&bckTrack);
698
699   TTree *tpcTree=(TTree*)in->Get("TPCf");
700   if (!tpcTree) {
701      cerr<<"AliTPCtracker::PropagateBack() ";
702      cerr<<"can't get a tree with TPC tracks !\n";
703      return 4;
704   }
705   AliTPCtrack *tpcTrack=new AliTPCtrack; 
706   tpcTree->SetBranchAddress("tracks",&tpcTrack);
707
708 //*** Prepare an array of tracks to be back propagated
709   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
710   Int_t nrows=nlow+nup;
711
712   TObjArray tracks(15000);
713   Int_t i=0,j=0;
714   Int_t tpcN=(Int_t)tpcTree->GetEntries();
715   Int_t bckN=(Int_t)bckTree->GetEntries();
716   if (j<bckN) bckTree->GetEvent(j++);
717   for (i=0; i<tpcN; i++) {
718      tpcTree->GetEvent(i);
719      Double_t alpha=tpcTrack->GetAlpha();
720      if (j<bckN &&
721      TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
722         if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
723         fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
724         bckTree->GetEvent(j++);
725      } else {
726         tpcTrack->ResetCovariance();
727         fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
728      }
729      tracks.AddLast(new AliTPCtrack(*tpcTrack));
730   }
731
732   out->cd();
733   TTree backTree("TPCb","Tree with back propagated TPC tracks");
734   AliTPCtrack *otrack=0;
735   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
736
737 //*** Back propagation through inner sectors
738   LoadInnerSectors();
739   Int_t nseed=fSeeds->GetEntriesFast();
740   for (i=0; i<nseed; i++) {
741      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
742      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
743
744      Int_t nc=t.GetNumberOfClusters();
745      s.SetLabel(nc-1); //set number of the cluster to start with
746
747      if (FollowBackProlongation(s,t)) {
748         UseClusters(&s);
749         continue;
750      }
751      delete fSeeds->RemoveAt(i);
752   }  
753   UnloadInnerSectors();
754
755 //*** Back propagation through outer sectors
756   LoadOuterSectors();
757   Int_t found=0;
758   for (i=0; i<nseed; i++) {
759      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
760      if (!ps) continue;
761      Int_t nc=s.GetNumberOfClusters();
762      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
763
764      Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
765      if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
766      if (alpha < 0.            ) alpha += 2.*TMath::Pi();
767      Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
768
769      alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
770      alpha-=s.GetAlpha();
771
772      if (s.Rotate(alpha)) {
773         if (FollowBackProlongation(s,t)) {
774            if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
775               s.CookdEdx();
776               s.SetLabel(t.GetLabel());
777               UseClusters(&s,nc);
778               otrack=ps;
779               backTree.Fill();
780               cerr<<found++<<'\r';
781               continue;
782            }
783         }
784      }
785      delete fSeeds->RemoveAt(i);
786   }  
787   UnloadOuterSectors();
788
789   backTree.Write();
790   savedir->cd();
791   cerr<<"Number of seeds: "<<nseed<<endl;
792   cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
793   cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
794
795   delete bckTrack;
796   delete tpcTrack;
797
798   return 0;
799 }
800
801 //_________________________________________________________________________
802 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
803   //--------------------------------------------------------------------
804   //       Return pointer to a given cluster
805   //--------------------------------------------------------------------
806   Int_t sec=(index&0xff000000)>>24; 
807   Int_t row=(index&0x00ff0000)>>16; 
808   Int_t ncl=(index&0x0000ffff)>>00;
809
810   AliTPCClustersRow *clrow=((AliTPCtracker *) this)->fClustersArray.GetRow(sec,row);
811   return (AliCluster*)(*clrow)[ncl];      
812 }
813
814 //__________________________________________________________________________
815 void AliTPCtracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
816   //--------------------------------------------------------------------
817   //This function "cooks" a track label. If label<0, this track is fake.
818   //--------------------------------------------------------------------
819   Int_t noc=t->GetNumberOfClusters();
820   Int_t *lb=new Int_t[noc];
821   Int_t *mx=new Int_t[noc];
822   AliCluster **clusters=new AliCluster*[noc];
823
824   Int_t i;
825   for (i=0; i<noc; i++) {
826      lb[i]=mx[i]=0;
827      Int_t index=t->GetClusterIndex(i);
828      clusters[i]=GetCluster(index);
829   }
830
831   Int_t lab=123456789;
832   for (i=0; i<noc; i++) {
833     AliCluster *c=clusters[i];
834     lab=TMath::Abs(c->GetLabel(0));
835     Int_t j;
836     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
837     lb[j]=lab;
838     (mx[j])++;
839   }
840
841   Int_t max=0;
842   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
843     
844   for (i=0; i<noc; i++) {
845     AliCluster *c=clusters[i];
846     if (TMath::Abs(c->GetLabel(1)) == lab ||
847         TMath::Abs(c->GetLabel(2)) == lab ) max++;
848   }
849
850   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
851
852   else {
853      Int_t tail=Int_t(0.10*noc);
854      max=0;
855      for (i=1; i<=tail; i++) {
856        AliCluster *c=clusters[noc-i];
857        if (lab == TMath::Abs(c->GetLabel(0)) ||
858            lab == TMath::Abs(c->GetLabel(1)) ||
859            lab == TMath::Abs(c->GetLabel(2))) max++;
860      }
861      if (max < Int_t(0.5*tail)) lab=-lab;
862   }
863
864   t->SetLabel(lab);
865
866   delete[] lb;
867   delete[] mx;
868   delete[] clusters;
869 }
870
871 //_________________________________________________________________________
872 void AliTPCtracker::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
873   //-----------------------------------------------------------------------
874   // Setup inner sector
875   //-----------------------------------------------------------------------
876   if (f==0) {
877      fAlpha=par->GetInnerAngle();
878      fAlphaShift=par->GetInnerAngleShift();
879      fPadPitchWidth=par->GetInnerPadPitchWidth();
880      fPadPitchLength=par->GetInnerPadPitchLength();
881
882      fN=par->GetNRowLow();
883      fRow=new AliTPCRow[fN];
884      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
885   } else {
886      fAlpha=par->GetOuterAngle();
887      fAlphaShift=par->GetOuterAngleShift();
888      fPadPitchWidth=par->GetOuterPadPitchWidth();
889      fPadPitchLength=par->GetOuterPadPitchLength();
890
891      fN=par->GetNRowUp();
892      fRow=new AliTPCRow[fN];
893      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiUp(i));
894   } 
895 }
896
897 //_________________________________________________________________________
898 void 
899 AliTPCtracker::AliTPCRow::InsertCluster(const AliTPCcluster* c, UInt_t index) {
900   //-----------------------------------------------------------------------
901   // Insert a cluster into this pad row in accordence with its y-coordinate
902   //-----------------------------------------------------------------------
903   if (fN==kMaxClusterPerRow) {
904     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
905   }
906   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
907   Int_t i=Find(c->GetY());
908   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
909   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
910   fIndex[i]=index; fClusters[i]=c; fN++;
911 }
912
913 //___________________________________________________________________
914 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
915   //-----------------------------------------------------------------------
916   // Return the index of the nearest cluster 
917   //-----------------------------------------------------------------------
918   if (y <= fClusters[0]->GetY()) return 0;
919   if (y > fClusters[fN-1]->GetY()) return fN;
920   Int_t b=0, e=fN-1, m=(b+e)/2;
921   for (; b<e; m=(b+e)/2) {
922     if (y > fClusters[m]->GetY()) b=m+1;
923     else e=m; 
924   }
925   return m;
926 }
927
928 //_____________________________________________________________________________
929 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
930   //-----------------------------------------------------------------
931   // This funtion calculates dE/dX within the "low" and "up" cuts.
932   //-----------------------------------------------------------------
933   Int_t i;
934   Int_t nc=GetNumberOfClusters();
935
936   Int_t swap;//stupid sorting
937   do {
938     swap=0;
939     for (i=0; i<nc-1; i++) {
940       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
941       Float_t tmp=fdEdxSample[i];
942       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
943       swap++;
944     }
945   } while (swap);
946
947   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
948   Float_t dedx=0;
949   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
950   dedx /= (nu-nl+1);
951   SetdEdx(dedx);
952 }
953
954
955