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