stdlib.h included to declare exit() on HP
[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.23  2002/11/19 11:50:08  hristov
19 Removing CONTAINERS (Yu.Belikov)
20
21 Revision 1.19  2002/07/19 07:31:40  kowal2
22 Improvement in tracking by J. Belikov
23
24 Revision 1.18  2002/05/13 07:33:52  kowal2
25 Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
26 in the case of defined region of interests.
27
28 Revision 1.17  2002/03/18 17:59:13  kowal2
29 Chnges in the pad geometry - 3 pad lengths introduced.
30
31 Revision 1.16  2001/11/08 16:39:03  hristov
32 Additional protection (M.Masera)
33
34 Revision 1.15  2001/11/08 16:36:33  hristov
35 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)
36
37 Revision 1.14  2001/10/21 19:04:55  hristov
38 Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
39
40 Revision 1.13  2001/07/23 12:01:30  hristov
41 Initialisation added
42
43 Revision 1.12  2001/07/20 14:32:44  kowal2
44 Processing of many events possible now
45
46 Revision 1.11  2001/05/23 08:50:10  hristov
47 Weird inline removed
48
49 Revision 1.10  2001/05/16 14:57:25  alibrary
50 New files for folders and Stack
51
52 Revision 1.9  2001/05/11 07:16:56  hristov
53 Fix needed on Sun and Alpha
54
55 Revision 1.8  2001/05/08 15:00:15  hristov
56 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
57
58 Revision 1.5  2000/12/20 07:51:59  kowal2
59 Changes suggested by Alessandra and Paolo to avoid overlapped
60 data fields in encapsulated classes.
61
62 Revision 1.4  2000/11/02 07:27:16  kowal2
63 code corrections
64
65 Revision 1.2  2000/06/30 12:07:50  kowal2
66 Updated from the TPC-PreRelease branch
67
68 Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
69 Splitted from AliTPCtracking
70
71 */
72
73 //-------------------------------------------------------
74 //          Implementation of the TPC tracker
75 //
76 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
77 //-------------------------------------------------------
78 #include <TObjArray.h>
79 #include <TFile.h>
80 #include <TTree.h>
81 #include <Riostream.h>
82
83 #include "AliTPCtracker.h"
84 #include "AliTPCcluster.h"
85 #include "AliTPCParam.h"
86 #include "AliClusters.h"
87
88 #include <stdlib.h>
89 //_____________________________________________________________________________
90 AliTPCtracker::AliTPCtracker(const AliTPCParam *par): 
91 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
92 {
93   //---------------------------------------------------------------------
94   // The main TPC tracker constructor
95   //---------------------------------------------------------------------
96   fInnerSec=new AliTPCSector[fkNIS];         
97   fOuterSec=new AliTPCSector[fkNOS];
98
99   Int_t i;
100   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
101   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
102
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::LoadClusters() {
213   //-----------------------------------------------------------------
214   // This function loads TPC clusters.
215   //-----------------------------------------------------------------
216   if (!gFile->IsOpen()) {
217     cerr<<"AliTPCtracker::LoadClusters : "<<
218       "file with clusters has not been open !\n";
219     return;
220   }
221
222   Char_t name[100];
223   sprintf(name,"TreeC_TPC_%d",GetEventNumber());
224   TTree *cTree=(TTree*)gFile->Get(name);
225   if (!cTree) {
226     cerr<<"AliTPCtracker::LoadClusters : "<<
227       "can't get the tree with TPC clusters !\n";
228     return;
229   }
230
231   TBranch *branch=cTree->GetBranch("Segment");
232   if (!branch) {
233     cerr<<"AliTPCtracker::LoadClusters : "<<
234       "can't get the segment branch !\n";
235     return;
236   }
237   AliClusters carray, *addr=&carray; carray.SetClass("AliTPCcluster");
238   branch->SetAddress(&addr);
239
240   Int_t nentr=(Int_t)cTree->GetEntries();
241
242   for (Int_t i=0; i<nentr; i++) {
243       cTree->GetEvent(i);
244
245       Int_t ncl=carray.GetArray()->GetEntriesFast();
246
247       Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
248       Int_t id=carray.GetID();
249       if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
250          cerr<<"AliTPCtracker::LoadClusters : "<<
251                "wrong index !\n";
252          exit(1);
253       }        
254       Int_t outindex = 2*fkNIS*nir;
255       if (id<outindex) {
256          Int_t sec = id/nir;
257          Int_t row = id - sec*nir;
258          sec %= fkNIS;
259          AliTPCRow &padrow=fInnerSec[sec][row];
260          while (ncl--) {
261            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
262            padrow.InsertCluster(c,sec,row);
263          }           
264       } else {
265          id -= outindex;
266          Int_t sec = id/nor;
267          Int_t row = id - sec*nor;
268          sec %= fkNOS;
269          AliTPCRow &padrow=fOuterSec[sec][row];
270          while (ncl--) {
271            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
272            padrow.InsertCluster(c,sec+fkNIS,row);
273          }
274       }
275
276       carray.GetArray()->Clear();
277   }
278   delete cTree;
279 }
280
281 //_____________________________________________________________________________
282 void AliTPCtracker::UnloadClusters() {
283   //-----------------------------------------------------------------
284   // This function unloads TPC clusters.
285   //-----------------------------------------------------------------
286   Int_t i;
287   for (i=0; i<fkNIS; i++) {
288     Int_t nr=fInnerSec->GetNRows();
289     for (Int_t n=0; n<nr; n++) fInnerSec[i][n].ResetClusters();
290   }
291   for (i=0; i<fkNOS; i++) {
292     Int_t nr=fOuterSec->GetNRows();
293     for (Int_t n=0; n<nr; n++) fOuterSec[i][n].ResetClusters();
294   }
295 }
296
297 //_____________________________________________________________________________
298 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
299   //-----------------------------------------------------------------
300   // This function tries to find a track prolongation.
301   //-----------------------------------------------------------------
302   Double_t xt=t.GetX();
303   const Int_t kSKIP=(t.GetNumberOfClusters()<10) ? kRowsToSkip : 
304                     Int_t(0.5*fSectors->GetNRows());
305   Int_t tryAgain=kSKIP;
306
307   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
308   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
309   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
310   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
311
312   Int_t nrows=fSectors->GetRowNumber(xt)-1;
313   for (Int_t nr=nrows; nr>=rf; nr--) {
314     Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
315     if (!t.PropagateTo(x)) return 0;
316
317     AliTPCcluster *cl=0;
318     UInt_t index=0;
319     Double_t maxchi2=kMaxCHI2;
320     const AliTPCRow &krow=fSectors[s][nr];
321     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
322     Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
323     Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
324     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
325
326     if (road>kMaxROAD) {
327       if (t.GetNumberOfClusters()>4) 
328         cerr<<t.GetNumberOfClusters()
329         <<"FindProlongation warning: Too broad road !\n"; 
330       return 0;
331     }
332
333     if (krow) {
334       for (Int_t i=krow.Find(y-road); i<krow; i++) {
335         AliTPCcluster *c=(AliTPCcluster*)(krow[i]);
336         if (c->GetY() > y+road) break;
337         if (c->IsUsed()) continue;
338         if ((c->GetZ()-z)*(c->GetZ()-z) > 16.*(t.GetSigmaZ2()+sz2)) continue;
339         Double_t chi2=t.GetPredictedChi2(c);
340         if (chi2 > maxchi2) continue;
341         maxchi2=chi2;
342         cl=c;
343         index=krow.GetIndex(i);       
344       }
345     }
346     if (cl) {
347       Float_t l=fSectors->GetPadPitchWidth();
348       Float_t corr=1.; if (nr>63) corr=0.67; // new (third) pad response !
349       t.SetSampledEdx(cl->GetQ()/l*corr,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 >= fkNIS) row=-1; } 
387   else { if (sec <  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 >= fkNIS) row=-1; }
432        else { if (sec < 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       Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
457       seed.SetSampledEdx(cl->GetQ()/l*corr,seed.GetNumberOfClusters());
458       seed.Update(cl,maxchi2,index);
459     }
460
461   }
462
463   seed.SetLabel(nc);
464
465   return 1;
466 }
467
468 //_____________________________________________________________________________
469 void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
470   //-----------------------------------------------------------------
471   // This function creates track seeds.
472   //-----------------------------------------------------------------
473   Double_t x[5], c[15];
474
475   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
476   Double_t cs=cos(alpha), sn=sin(alpha);
477
478   Double_t x1 =fSectors->GetX(i1);
479   Double_t xx2=fSectors->GetX(i2);
480
481   for (Int_t ns=0; ns<fN; ns++) {
482     Int_t nl=fSectors[(ns-1+fN)%fN][i2];
483     Int_t nm=fSectors[ns][i2];
484     Int_t nu=fSectors[(ns+1)%fN][i2];
485     const AliTPCRow& kr1=fSectors[ns][i1];
486     for (Int_t is=0; is < kr1; is++) {
487       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
488       for (Int_t js=0; js < nl+nm+nu; js++) {
489         const AliTPCcluster *kcl;
490         Double_t x2,   y2,   z2;
491         Double_t x3=GetX(), y3=GetY(), z3=GetZ();
492
493         if (js<nl) {
494           const AliTPCRow& kr2=fSectors[(ns-1+fN)%fN][i2];
495           kcl=kr2[js];
496           y2=kcl->GetY(); z2=kcl->GetZ();
497           x2= xx2*cs+y2*sn;
498           y2=-xx2*sn+y2*cs;
499         } else 
500           if (js<nl+nm) {
501             const AliTPCRow& kr2=fSectors[ns][i2];
502             kcl=kr2[js-nl];
503             x2=xx2; y2=kcl->GetY(); z2=kcl->GetZ();
504           } else {
505             const AliTPCRow& kr2=fSectors[(ns+1)%fN][i2];
506             kcl=kr2[js-nl-nm];
507             y2=kcl->GetY(); z2=kcl->GetZ();
508             x2=xx2*cs-y2*sn;
509             y2=xx2*sn+y2*cs;
510           }
511
512         Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2); 
513         if (TMath::Abs(zz-z2)>5.) continue;
514
515         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
516         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
517
518         x[0]=y1;
519         x[1]=z1;
520         x[4]=f1(x1,y1,x2,y2,x3,y3);
521         if (TMath::Abs(x[4]) >= 0.0066) continue;
522         x[2]=f2(x1,y1,x2,y2,x3,y3);
523         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
524         x[3]=f3(x1,y1,x2,y2,z1,z2);
525         if (TMath::Abs(x[3]) > 1.2) continue;
526         Double_t a=asin(x[2]);
527         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
528         if (TMath::Abs(zv-z3)>10.) continue; 
529
530         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
531         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
532         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
533         Double_t sy3=25000*x[4]*x[4]+0.1, 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=fSectors->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   LoadClusters();
629
630   out->cd();
631
632   char   tname[100];
633   sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
634   TTree tracktree(tname,"Tree with TPC tracks");
635   AliTPCtrack *iotrack=0;
636   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
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      fSectors=fOuterSec; fN=fkNOS;
644      fSeeds=new TObjArray(15000);     
645      MakeSeeds(nup-1, nup-1-gap);
646      MakeSeeds(nup-1-shift, nup-1-shift-gap);
647   }
648   fSeeds->Sort();
649
650   Int_t found=0;
651   Int_t nseed=fSeeds->GetEntriesFast();
652   for (Int_t i=0; i<nseed; i++) {
653     //tracking in the outer sectors
654     fSectors=fOuterSec; fN=fkNOS;
655
656     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
657     if (!FollowProlongation(t)) {
658        delete fSeeds->RemoveAt(i);
659        continue;
660     }
661
662     //tracking in the inner sectors
663     fSectors=fInnerSec; fN=fkNIS;
664
665     Double_t alpha=t.GetAlpha() - fInnerSec->GetAlphaShift();
666     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
667     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
668     Int_t ns=Int_t(alpha/fInnerSec->GetAlpha())%fkNIS;
669
670     alpha=ns*fInnerSec->GetAlpha()+fInnerSec->GetAlphaShift()-t.GetAlpha();
671
672     if (t.Rotate(alpha)) {
673        if (FollowProlongation(t)) {
674           if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
675              t.CookdEdx();
676              CookLabel(pt,0.1); //For comparison only
677              iotrack=pt;
678              tracktree.Fill();
679              UseClusters(&t);
680              cerr<<found++<<'\r';
681           }
682        }
683     }
684     delete fSeeds->RemoveAt(i); 
685   }
686
687   tracktree.Write();
688
689   cerr<<"Number of found tracks : "<<found<<endl;
690
691   savedir->cd();
692
693   UnloadClusters();
694   fSeeds->Clear(); delete fSeeds; fSeeds=0;
695
696   return 0;
697 }
698
699 //_____________________________________________________________________________
700 Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
701   //-----------------------------------------------------------------
702   // This function propagates tracks back through the TPC.
703   //-----------------------------------------------------------------
704   fSeeds=new TObjArray(15000);
705
706   TFile *in=(TFile*)inp;
707   TDirectory *savedir=gDirectory; 
708
709   if (!in->IsOpen()) {
710      cerr<<"AliTPCtracker::PropagateBack(): ";
711      cerr<<"file with back propagated ITS tracks is not open !\n";
712      return 1;
713   }
714
715   if (!out->IsOpen()) {
716      cerr<<"AliTPCtracker::PropagateBack(): ";
717      cerr<<"file for back propagated TPC tracks is not open !\n";
718      return 2;
719   }
720
721   LoadClusters();
722
723   in->cd();
724   TTree *bckTree=(TTree*)in->Get("TreeT_ITSb_0");
725   if (!bckTree) {
726      cerr<<"AliTPCtracker::PropagateBack() ";
727      cerr<<"can't get a tree with back propagated ITS tracks !\n";
728      return 3;
729   }
730   AliTPCtrack *bckTrack=new AliTPCtrack; 
731   bckTree->SetBranchAddress("tracks",&bckTrack);
732
733   TTree *tpcTree=(TTree*)in->Get("TreeT_TPC_0");
734   if (!tpcTree) {
735      cerr<<"AliTPCtracker::PropagateBack() ";
736      cerr<<"can't get a tree with TPC tracks !\n";
737      return 4;
738   }
739   AliTPCtrack *tpcTrack=new AliTPCtrack; 
740   tpcTree->SetBranchAddress("tracks",&tpcTrack);
741
742 //*** Prepare an array of tracks to be back propagated
743   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
744   Int_t nrows=nlow+nup;
745
746   TObjArray tracks(15000);
747   Int_t i=0,j=0;
748   Int_t tpcN=(Int_t)tpcTree->GetEntries();
749   Int_t bckN=(Int_t)bckTree->GetEntries();
750   if (j<bckN) bckTree->GetEvent(j++);
751   for (i=0; i<tpcN; i++) {
752      tpcTree->GetEvent(i);
753      Double_t alpha=tpcTrack->GetAlpha();
754      if (j<bckN &&
755      TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
756         if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
757         fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
758         bckTree->GetEvent(j++);
759      } else {
760         tpcTrack->ResetCovariance();
761         fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
762      }
763      tracks.AddLast(new AliTPCtrack(*tpcTrack));
764   }
765
766   out->cd();
767   TTree backTree("TreeT_TPCb_0","Tree with back propagated TPC tracks");
768   AliTPCtrack *otrack=0;
769   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
770
771 //*** Back propagation through inner sectors
772   fSectors=fInnerSec; fN=fkNIS;
773   Int_t nseed=fSeeds->GetEntriesFast();
774   for (i=0; i<nseed; i++) {
775      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
776      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
777
778      Int_t nc=t.GetNumberOfClusters();
779      s.SetLabel(nc-1); //set number of the cluster to start with
780
781      if (FollowBackProlongation(s,t)) {
782         UseClusters(&s);
783         continue;
784      }
785      delete fSeeds->RemoveAt(i);
786   }  
787
788 //*** Back propagation through outer sectors
789   fSectors=fOuterSec; fN=fkNOS;
790   Int_t found=0;
791   for (i=0; i<nseed; i++) {
792      AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
793      if (!ps) continue;
794      Int_t nc=s.GetNumberOfClusters();
795      const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
796
797      Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
798      if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
799      if (alpha < 0.            ) alpha += 2.*TMath::Pi();
800      Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
801
802      alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
803      alpha-=s.GetAlpha();
804
805      if (s.Rotate(alpha)) {
806         if (FollowBackProlongation(s,t)) {
807            if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
808               s.CookdEdx();
809               s.SetLabel(t.GetLabel());
810               UseClusters(&s,nc);
811               otrack=ps;
812               backTree.Fill();
813               cerr<<found++<<'\r';
814               continue;
815            }
816         }
817      }
818      delete fSeeds->RemoveAt(i);
819   }  
820
821   backTree.Write();
822   savedir->cd();
823   cerr<<"Number of seeds: "<<nseed<<endl;
824   cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
825   cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
826
827   delete bckTrack;
828   delete tpcTrack;
829
830   delete bckTree; //Thanks to Mariana Bondila
831   delete tpcTree; //Thanks to Mariana Bondila
832
833   UnloadClusters();  
834
835   return 0;
836 }
837
838 //_________________________________________________________________________
839 AliCluster *AliTPCtracker::GetCluster(Int_t index) const {
840   //--------------------------------------------------------------------
841   //       Return pointer to a given cluster
842   //--------------------------------------------------------------------
843   Int_t sec=(index&0xff000000)>>24; 
844   Int_t row=(index&0x00ff0000)>>16; 
845   Int_t ncl=(index&0x0000ffff)>>00;
846
847   const AliTPCcluster *cl=0;
848
849   if (sec<fkNIS) {
850     cl=fInnerSec[sec][row].GetUnsortedCluster(ncl); 
851   } else {
852     sec-=fkNIS;
853     cl=fOuterSec[sec][row].GetUnsortedCluster(ncl);
854   }
855
856   return (AliCluster*)cl;      
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      f1PadPitchLength=par->GetInnerPadPitchLength();
926      f2PadPitchLength=f1PadPitchLength;
927
928      fN=par->GetNRowLow();
929      fRow=new AliTPCRow[fN];
930      for (Int_t i=0; i<fN; i++) fRow[i].SetX(par->GetPadRowRadiiLow(i));
931   } else {
932      fAlpha=par->GetOuterAngle();
933      fAlphaShift=par->GetOuterAngleShift();
934      fPadPitchWidth=par->GetOuterPadPitchWidth();
935      f1PadPitchLength=par->GetOuter1PadPitchLength();
936      f2PadPitchLength=par->GetOuter2PadPitchLength();
937
938      fN=par->GetNRowUp();
939      fRow=new AliTPCRow[fN];
940      for (Int_t i=0; i<fN; i++){ 
941      fRow[i].SetX(par->GetPadRowRadiiUp(i));
942      }
943   } 
944 }
945
946 //_________________________________________________________________________
947 void AliTPCtracker::
948 AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
949   //-----------------------------------------------------------------------
950   // Insert a cluster into this pad row in accordence with its y-coordinate
951   //-----------------------------------------------------------------------
952   if (fN==kMaxClusterPerRow) {
953     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
954   }
955
956   Int_t index=(((sec<<8)+row)<<16)+fN;
957
958   if (fN==0) {
959      fIndex[0]=index;
960      fClusterArray[0]=*c; fClusters[fN++]=fClusterArray; 
961      return;
962   }
963
964   if (fN==fSize) {
965      Int_t size=fSize*2;
966      AliTPCcluster *buff=new AliTPCcluster[size];
967      memcpy(buff,fClusterArray,fSize*sizeof(AliTPCcluster));
968      for (Int_t i=0; i<fN; i++) 
969         fClusters[i]=buff+(fClusters[i]-fClusterArray);
970      delete[] fClusterArray;
971      fClusterArray=buff;
972      fSize=size;
973   }
974
975   Int_t i=Find(c->GetY());
976   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCcluster*));
977   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
978   fIndex[i]=index; 
979   fClusters[i]=fClusterArray+fN; fClusterArray[fN++]=*c;
980 }
981
982 //___________________________________________________________________
983 Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const {
984   //-----------------------------------------------------------------------
985   // Return the index of the nearest cluster 
986   //-----------------------------------------------------------------------
987   if(fN<=0) return 0;
988   if (y <= fClusters[0]->GetY()) return 0;
989   if (y > fClusters[fN-1]->GetY()) return fN;
990   Int_t b=0, e=fN-1, m=(b+e)/2;
991   for (; b<e; m=(b+e)/2) {
992     if (y > fClusters[m]->GetY()) b=m+1;
993     else e=m; 
994   }
995   return m;
996 }
997
998 //_____________________________________________________________________________
999 void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1000   //-----------------------------------------------------------------
1001   // This funtion calculates dE/dX within the "low" and "up" cuts.
1002   //-----------------------------------------------------------------
1003   Int_t i;
1004   Int_t nc=GetNumberOfClusters();
1005
1006   Int_t swap;//stupid sorting
1007   do {
1008     swap=0;
1009     for (i=0; i<nc-1; i++) {
1010       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
1011       Float_t tmp=fdEdxSample[i];
1012       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
1013       swap++;
1014     }
1015   } while (swap);
1016
1017   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
1018   Float_t dedx=0;
1019   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
1020   dedx /= (nu-nl+1);
1021   SetdEdx(dedx);
1022
1023   //Very rough PID
1024   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1025
1026   Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
1027   if (p<0.6) {
1028     if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
1029     if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
1030     SetMass(0.93827); return;
1031   }
1032
1033   if (p<1.2) {
1034     if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
1035     SetMass(0.93827); return;
1036   }
1037
1038   SetMass(0.13957); return;
1039
1040 }
1041
1042
1043