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