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