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