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