First version of the parallel TPC tracking (M.Ivanov)
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.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
18
19
20
21 /*
22   AliTPC parallel tracker - 
23   How to use?  - 
24   run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
25   run AliTPCFindTracksMI.C macro - to find tracks
26   tracks are written to AliTPCtracks.root file
27   for comparison also seeds are written to the same file - to special branch
28 */
29
30 //-------------------------------------------------------
31 //          Implementation of the TPC tracker
32 //
33 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
34 // 
35 //-------------------------------------------------------
36
37 #include <TObjArray.h>
38 #include <TFile.h>
39 #include <TTree.h>
40 #include <iostream.h>
41
42 #include "AliTPCtrackerMI.h"
43 #include "AliTPCclusterMI.h"
44 #include "AliTPCParam.h"
45 #include "AliTPCClustersRow.h"
46 #include "AliComplexCluster.h"
47 #include "TStopwatch.h"
48
49
50 ClassImp(AliTPCseed)
51
52 AliTPCclusterTracks::AliTPCclusterTracks(){
53   // class for storing overlaping info
54   fTrackIndex[0]=-1;
55   fTrackIndex[1]=-1;
56   fTrackIndex[2]=-1;
57   fDistance[0]=1000;
58   fDistance[1]=1000;
59   fDistance[2]=1000;
60 }
61
62
63
64
65
66 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){
67
68   Int_t sec=(i&0xff000000)>>24; 
69   Int_t row = (i&0x00ff0000)>>16;
70   track->fRow=(i&0x00ff0000)>>16;
71   track->fSector = sec;
72   //  Int_t index = i&0xFFFF;
73   if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow(); 
74   track->fClusterIndex[track->fRow] = i;
75   track->fFirstPoint = row;
76   //
77
78   AliTPCTrackPoint   *trpoint =track->GetTrackPoint(track->fRow);
79   Float_t angle2 = track->GetSnp()*track->GetSnp();
80   angle2 = TMath::Sqrt(angle2/(1-angle2)); 
81   //
82   //SET NEW Track Point
83   //
84   if (c!=0){
85     //if we have a cluster
86     trpoint->GetCPoint().SetY(c->GetY());
87     trpoint->GetCPoint().SetZ(c->GetZ());    
88     //
89     trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY));
90     trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ));
91     //
92     trpoint->GetCPoint().SetType(c->GetType());
93     trpoint->GetCPoint().SetQ(c->GetQ());
94     trpoint->GetCPoint().SetMax(c->GetMax());
95     //  
96     trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2));
97     trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2));
98     //
99   }
100   trpoint->GetTPoint().SetX(track->GetX());
101   trpoint->GetTPoint().SetY(track->GetY());
102   trpoint->GetTPoint().SetZ(track->GetZ());
103   //
104   trpoint->GetTPoint().SetAngleY(angle2);
105   trpoint->GetTPoint().SetAngleZ(track->GetTgl());
106   
107
108   
109   if (chi2>10){
110     //    printf("suspicious chi2 %f\n",chi2);
111   }
112   //  if (track->fIsSeeding){ 
113     track->fErrorY2 *= 1.2;
114     track->fErrorY2 += 0.0064;    
115     track->fErrorZ2 *= 1.2;   
116     track->fErrorY2 += 0.005;    
117  
118     //}
119
120   return track->Update(c,chi2,i);
121
122 }
123 //_____________________________________________________________________________
124 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn): 
125 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
126 {
127   //---------------------------------------------------------------------
128   // The main TPC tracker constructor
129   //---------------------------------------------------------------------
130   fInnerSec=new AliTPCSector[fkNIS];         
131   fOuterSec=new AliTPCSector[fkNOS];
132
133   Int_t i;
134   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
135   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
136
137   fN=0;  fSectors=0;
138
139   fClustersArray.Setup(par);
140   fClustersArray.SetClusterType("AliTPCclusterMI");
141
142   char   cname[100];
143   if (eventn==-1) {
144     sprintf(cname,"TreeC_TPC");
145   }
146   else {
147     sprintf(cname,"TreeC_TPC_%d",eventn);
148   }
149
150   fClustersArray.ConnectTree(cname);
151
152   fEventN = eventn;
153   fSeeds=0;
154   fNtracks = 0;
155   fParam = par;
156 }
157
158 //_____________________________________________________________________________
159 AliTPCtrackerMI::~AliTPCtrackerMI() {
160   //------------------------------------------------------------------
161   // TPC tracker destructor
162   //------------------------------------------------------------------
163   delete[] fInnerSec;
164   delete[] fOuterSec;
165   if (fSeeds) {
166     fSeeds->Delete(); 
167     delete fSeeds;
168   }
169 }
170
171
172 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
173   //
174   //
175   Float_t snoise2;
176   Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
177
178   //cluster "quality"
179   Float_t rsigmay = 1;
180   Int_t ctype = 0;
181
182   //standard if we don't have cluster - take MIP
183   const Float_t chmip = 50.; 
184   Float_t amp = chmip/0.3;  
185   Float_t nel;
186   Float_t nprim;
187   if (cl){
188     amp = cl->GetQ();
189     rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY);
190     ctype = cl->GetType();
191   }
192   
193
194   Float_t landau=2 ;    //landau fluctuation part
195   Float_t gg=2;         // gg fluctuation part
196   Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
197
198
199   if (fSectors==fInnerSec){
200     snoise2 = 0.0004;
201     nel     = 0.268*amp;
202     nprim   = 0.155*amp;
203     gg      = (2+0.0002*amp)/nel;
204     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
205     if (landau>1) landau=1;
206   }
207   else {
208     snoise2 = 0.0004;
209     nel     = 0.3*amp;
210     nprim   = 0.133*amp;
211     gg      = (2+0.0002*amp)/nel;
212     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
213     if (landau>1) landau=1;
214   }
215
216
217   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
218   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
219   angle2 = angle2/(1-angle2); 
220   Float_t angular = landau*angle2*padlength*padlength/12.;
221   Float_t res = sdiff + angular;
222   
223   
224   if ((ctype==0) && (fSectors ==fOuterSec))
225     res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2));
226
227   if ((ctype==0) && (fSectors ==fInnerSec))
228     res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2));
229
230
231   if ((ctype>0))
232     res*= TMath::Power((rsigmay+0.5),1.5)+0.0064;
233   
234   if (ctype<0)
235     res*=2.4;  // overestimate error 2 times
236   
237   res+= snoise2;
238  
239   if (res<2*snoise2)
240     res = 2*snoise2;
241
242   seed->SetErrorY2(res);
243   return res;
244 }
245
246
247
248
249 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
250   //
251   //
252   Float_t snoise2;
253   Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
254   //signal quality
255   Float_t rsigmaz=1;
256   Int_t ctype =0;
257
258   const Float_t chmip = 50.;
259   Float_t amp = chmip/0.3;  
260   Float_t nel;
261   Float_t nprim;
262   if (cl){
263     amp = cl->GetQ();
264     rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ);
265     ctype = cl->GetType();
266   }
267
268   //
269   Float_t landau=2 ;    //landau fluctuation part
270   Float_t gg=2;         // gg fluctuation part
271   Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
272
273   if (fSectors==fInnerSec){
274     snoise2 = 0.0004;
275     nel     = 0.268*amp;
276     nprim   = 0.155*amp;
277     gg      = (2+0.0002*amp)/nel;
278     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
279     if (landau>1) landau=1;
280   }
281   else {
282     snoise2 = 0.0004;
283     nel     = 0.3*amp;
284     nprim   = 0.133*amp;
285     gg      = (2+0.0002*amp)/nel;
286     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
287     if (landau>1) landau=1;
288   }
289   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
290   
291   Float_t angle = seed->GetTgl();
292   Float_t angular = landau*angle*angle*padlength*padlength/12.;
293   Float_t res = sdiff + angular;
294
295   if ((ctype==0) && (fSectors ==fOuterSec))
296     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
297
298   if ((ctype==0) && (fSectors ==fInnerSec))
299     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
300   if ((ctype>0))
301     res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064;  //0.31+0.147*ctype;
302   if (ctype<0)
303     res*=1.3;
304   if ((ctype<0) &&amp<70)
305     res*=1.3;  
306
307   res += snoise2;
308   if (res<2*snoise2)
309      res = 2*snoise2;
310
311   seed->SetErrorZ2(res);
312   return res;
313 }
314
315
316
317 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
318 {
319   //-----------------------------------------------------------------
320   // This function find proloncation of a track to a reference plane x=xk.
321   // doesn't change internal state of the track
322   //-----------------------------------------------------------------
323   
324   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
325   //  Double_t y1=fP0, z1=fP1;
326   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
327   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
328   
329   y = fP0;
330   z = fP1;
331   y += dx*(c1+c2)/(r1+r2);
332   z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
333   return 0;  
334 }
335
336
337 //_____________________________________________________________________________
338 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const 
339 {
340   //-----------------------------------------------------------------
341   // This function calculates a predicted chi2 increment.
342   //-----------------------------------------------------------------
343   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
344   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
345   r00+=fC00; r01+=fC10; r11+=fC11;
346
347   Double_t det=r00*r11 - r01*r01;
348   if (TMath::Abs(det) < 1.e-10) {
349     Int_t n=GetNumberOfClusters();
350     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
351     return 1e10;
352   }
353   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
354   
355   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
356   
357   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
358 }
359
360
361 //_________________________________________________________________________________________
362
363
364 Int_t AliTPCseed::Compare(const TObject *o) const {
365   //-----------------------------------------------------------------
366   // This function compares tracks according to the sector - for given sector according z
367   //-----------------------------------------------------------------
368   AliTPCseed *t=(AliTPCseed*)o;
369   if (t->fSector>fSector) return -1;
370   if (t->fSector<fSector) return 1;
371
372   Double_t z2 = t->GetZ();
373   Double_t z1 = GetZ();
374   if (z2>z1) return 1;
375   if (z2<z1) return -1;
376   return 0;
377 }
378
379
380
381
382
383
384 //_____________________________________________________________________________
385 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
386   //-----------------------------------------------------------------
387   // This function associates a cluster with this track.
388   //-----------------------------------------------------------------
389   //  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
390   //Double_t r00=sigmay2, r01=0., r11=sigmaz2;
391   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
392
393   r00+=fC00; r01+=fC10; r11+=fC11;
394   Double_t det=r00*r11 - r01*r01;
395   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
396
397   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
398   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
399   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
400   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
401   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
402
403   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
404   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
405   if (TMath::Abs(cur*fX-eta) >= 0.99999) {
406     Int_t n=GetNumberOfClusters();
407     if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
408     return 0;
409   }
410
411   fP0 += k00*dy + k01*dz;
412   fP1 += k10*dy + k11*dz;
413   fP2  = eta;
414   fP3 += k30*dy + k31*dz;
415   fP4  = cur;
416
417   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
418   Double_t c12=fC21, c13=fC31, c14=fC41;
419
420   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
421   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
422   fC40-=k00*c04+k01*c14; 
423
424   fC11-=k10*c01+k11*fC11;
425   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
426   fC41-=k10*c04+k11*c14; 
427
428   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
429   fC42-=k20*c04+k21*c14; 
430
431   fC33-=k30*c03+k31*c13;
432   fC43-=k40*c03+k41*c13; 
433
434   fC44-=k40*c04+k41*c14; 
435
436   Int_t n=GetNumberOfClusters();
437   fIndex[n]=index;
438   SetNumberOfClusters(n+1);
439   SetChi2(GetChi2()+chisq);
440
441   return 1;
442 }
443
444
445
446 //_____________________________________________________________________________
447 Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
448                    Double_t x2,Double_t y2,
449                    Double_t x3,Double_t y3) 
450 {
451   //-----------------------------------------------------------------
452   // Initial approximation of the track curvature
453   //-----------------------------------------------------------------
454   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
455   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
456                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
457   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
458                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
459
460   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
461
462   return -xr*yr/sqrt(xr*xr+yr*yr); 
463 }
464
465
466 //_____________________________________________________________________________
467 Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
468                    Double_t x2,Double_t y2,
469                    Double_t x3,Double_t y3) 
470 {
471   //-----------------------------------------------------------------
472   // Initial approximation of the track curvature times center of curvature
473   //-----------------------------------------------------------------
474   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
475   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
476                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
477   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
478                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
479
480   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
481   
482   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
483 }
484
485 //_____________________________________________________________________________
486 Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1, 
487                    Double_t x2,Double_t y2,
488                    Double_t z1,Double_t z2) 
489 {
490   //-----------------------------------------------------------------
491   // Initial approximation of the tangent of the track dip angle
492   //-----------------------------------------------------------------
493   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
494 }
495
496
497 void AliTPCtrackerMI::LoadClusters()
498 {
499   //
500   // load clusters to the memory
501   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
502   for (Int_t i=0; i<j; i++) {
503     fClustersArray.LoadEntry(i);
504   }
505 }
506
507 void AliTPCtrackerMI::UnloadClusters()
508 {
509   //
510   // load clusters to the memory
511   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
512   for (Int_t i=0; i<j; i++) {
513     fClustersArray.ClearSegment(i);
514   }
515 }
516
517
518
519 //_____________________________________________________________________________
520 void AliTPCtrackerMI::LoadOuterSectors() {
521   //-----------------------------------------------------------------
522   // This function fills outer TPC sectors with clusters.
523   //-----------------------------------------------------------------
524   UInt_t index;
525   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
526   for (Int_t i=0; i<j; i++) {
527     //  AliSegmentID *s=fClustersArray.LoadEntry(i);
528     AliSegmentID *s= fClustersArray.At(i);
529     Int_t sec,row;
530     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
531     par->AdjustSectorRow(s->GetID(),sec,row);
532     if (sec<fkNIS*2) continue;
533     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
534     Int_t ncl=clrow->GetArray()->GetEntriesFast();
535     while (ncl--) {
536       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
537       index=(((sec<<8)+row)<<16)+ncl;
538       fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
539     }
540   }  
541   fN=fkNOS;
542   fSectors=fOuterSec;
543 }
544
545
546 //_____________________________________________________________________________
547 void AliTPCtrackerMI::LoadInnerSectors() {
548   //-----------------------------------------------------------------
549   // This function fills inner TPC sectors with clusters.
550   //-----------------------------------------------------------------
551   UInt_t index;
552   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
553   for (Int_t i=0; i<j; i++) {
554     //   AliSegmentID *s=fClustersArray.LoadEntry(i);
555     AliSegmentID *s=fClustersArray.At(i);
556     Int_t sec,row;
557     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
558     par->AdjustSectorRow(s->GetID(),sec,row);
559     if (sec>=fkNIS*2) continue;
560     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
561     Int_t ncl=clrow->GetArray()->GetEntriesFast();
562     while (ncl--) {
563       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
564       index=(((sec<<8)+row)<<16)+ncl;
565       fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
566     }
567   }
568
569   fN=fkNIS;
570   fSectors=fInnerSec;
571 }
572
573 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
574   //-----------------------------------------------------------------
575   // This function tries to find a track prolongation to next pad row
576   //-----------------------------------------------------------------
577   //  Double_t xt=t.GetX();
578   //  Int_t row = fSectors->GetRowNumber(xt)-1;
579   //  if (row < nr) return 1; // don't prolongate if not information until now -
580   //
581   Double_t  x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
582   if (!t.PropagateTo(x)) return 0;
583   // update current
584   t.fCurrentSigmaY = GetSigmaY(&t);
585   t.fCurrentSigmaZ = GetSigmaZ(&t);
586   //  
587   AliTPCclusterMI *cl=0;
588   UInt_t index=0;
589   const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
590   Double_t sy2=ErrY2(&t)*2;
591   Double_t sz2=ErrZ2(&t)*2;
592
593
594   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
595   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
596   Double_t  y=t.GetY(), z=t.GetZ();
597
598   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
599     t.fInDead = kTRUE;
600     return 0;
601   } 
602   else
603     {
604       t.fNFoundable++;
605     }   
606   //calculate 
607   Float_t maxdistance = roady*roady + roadz*roadz;
608   if (krow) {
609     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
610       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
611       if (c->GetZ() > z+roadz) break;
612       if ( (c->GetY()-y) >  roady ) continue;
613       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
614       if (maxdistance>distance) {
615         maxdistance = distance;
616         cl=c;
617         index=krow.GetIndex(i);       
618       }
619     }
620   }      
621   if (cl) {
622     //Double_t sy2=
623     ErrY2(&t,cl);
624     //Double_t sz2=
625     ErrZ2(&t,cl);
626     Double_t chi2= t.GetPredictedChi2(cl);    
627     UpdateTrack(&t,cl,chi2,index);   
628   } else {    
629     if (y > ymax) {
630       t.fRelativeSector= (t.fRelativeSector+1) % fN;
631       if (!t.Rotate(fSectors->GetAlpha())) 
632         return 0;
633     } else if (y <-ymax) {
634       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
635       if (!t.Rotate(-fSectors->GetAlpha())) 
636         return 0;
637     }   
638   }
639   return 1;
640 }
641
642
643 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex,  Int_t nr) {
644   //-----------------------------------------------------------------
645   // This function tries to find a track prolongation to next pad row
646   //-----------------------------------------------------------------
647   t.fCurrentCluster  = 0;
648   t.fCurrentClusterIndex1 = 0;   
649   t.fCurrentClusterIndex2 = 0;
650    
651   Double_t xt=t.GetX();
652   Int_t row = fSectors->GetRowNumber(xt)-1;
653   if (row < nr) return 1; // don't prolongate if not information until now -
654   Double_t x=fSectors->GetX(nr);
655   if (!t.PropagateTo(x)) return 0;
656   // update current
657   t.fCurrentSigmaY = GetSigmaY(&t);
658   t.fCurrentSigmaZ = GetSigmaZ(&t);
659     
660   AliTPCclusterMI *cl=0;
661   UInt_t index=0;
662   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
663   //
664   Double_t  y=t.GetY(), z=t.GetZ();
665   Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
666   Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
667   //
668
669   Float_t maxdistance = 1000000;
670   if (krow) {    
671     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
672       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
673       if (c->GetZ() > z+roadz) break;
674       if (TMath::Abs(c->GetY()-y)>roady) continue;
675             
676       //krow.UpdateClusterTrack(i,trindex,&t);
677
678       Float_t dy2 = (c->GetY()- t.GetY());
679       dy2*=dy2;
680       Float_t dz2 = (c->GetZ()- t.GetZ());
681       dz2*=dz2;
682       //
683       Float_t distance = dy2+dz2;
684       //
685       if (distance > maxdistance) continue;
686       maxdistance = distance;
687       cl=c;
688       index=i;       
689     }
690   }
691   t.fCurrentCluster  = cl;
692   t.fCurrentClusterIndex1 = krow.GetIndex(index);   
693   t.fCurrentClusterIndex2 = index;   
694   return 1;
695 }
696
697
698 Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
699   //-----------------------------------------------------------------
700   // This function tries to find a track prolongation to next pad row
701   //-----------------------------------------------------------------
702   AliTPCseed & t  = *((AliTPCseed*)(fSeeds->At(trindex)));
703   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
704   //  Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
705   Double_t y=t.GetY();
706   Double_t ymax=fSectors->GetMaxY(nr);
707
708   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
709     t.fInDead = kTRUE;
710     return 0;
711   } 
712   else
713     {
714       t.fNFoundable++;
715     }
716   
717   if (t.fCurrentCluster) {
718     //    Float_t l=fSectors->GetPadPitchLength();
719     //    AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
720
721     Double_t sy2=ErrY2(&t,t.fCurrentCluster);
722     Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
723
724
725     Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
726     Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
727
728     Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
729     Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
730     
731     Double_t rdistance  = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
732
733
734     //    printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
735     if (rdistance>4) return 0;
736
737     if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)  
738         return 0;  //suspisiouce - will be changed
739
740     if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)  
741         // strict cut on overlaped cluster
742         return 0;  //suspisiouce - will be changed
743
744     if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 ) 
745          && t.fCurrentCluster->GetType()<0)
746       return 0;
747
748     //    t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
749     UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
750    
751   } else {
752     if (y > ymax) {
753       t.fRelativeSector= (t.fRelativeSector+1) % fN;
754       if (!t.Rotate(fSectors->GetAlpha())) 
755         return 0;
756     } else if (y <-ymax) {
757       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
758       if (!t.Rotate(-fSectors->GetAlpha())) 
759         return 0;
760     }
761   }
762   return 1;
763 }
764
765
766
767
768 /*
769 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
770 {
771   //-----------------------------------------------------------------
772   // fast prolongation mathod -
773   // don't update track only after step clusters
774   //-----------------------------------------------------------------
775   Double_t xt=t.GetX();
776   //
777   Double_t alpha=t.GetAlpha(); 
778   alpha =- fSectors->GetAlphaShift();
779   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
780   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
781   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
782   Int_t row0 = fSectors->GetRowNumber(xt); 
783   Double_t x    = fSectors->GetX(row0);
784   Double_t ymax = fSectors->GetMaxY(row0);
785   //
786   Double_t sy2=ErrY2(&t)*2;
787   Double_t sz2=ErrZ2(&t)*2;
788   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
789   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
790   Float_t maxdistance = roady*roady + roadz*roadz; 
791   t.fCurrentSigmaY = GetSigmaY(&t);
792   t.fCurrentSigmaZ = GetSigmaZ(&t);
793   //
794   Int_t nclusters = 0;
795   Double_t y;
796   Double_t z;
797   Double_t yy[200];   //track prolongation
798   Double_t zz[200];
799   Double_t cy[200];  // founded cluster position
800   Double_t cz[200];
801   Double_t sy[200];  // founded cluster error
802   Double_t sz[200];
803   Bool_t   hitted[200]; // indication of cluster presence
804   //
805   
806   //
807   for (Int_t drow = step; drow>=0; drow--) {
808     Int_t row = row0-drow;
809     if (row<0) break;
810     Double_t x    = fSectors->GetX(row);
811     Double_t ymax = fSectors->GetMaxY(row);
812     t.GetProlongation(x,y,z);
813     yy[drow] =y;
814     zz[drow] =z;    
815     const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
816     if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
817       t.fInDead = kTRUE;
818       break;
819     } 
820     else
821       {
822         t.fNFoundable++;
823       }  
824     
825     //find nearest  cluster 
826     AliTPCclusterMI *cl= 0;
827     if (krow) {
828       for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
829         AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
830         if (c->GetZ() > z+roadz) break;
831         if ( (c->GetY()-y) >  roady ) continue;
832         Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
833         if (maxdistance>distance) {
834           maxdistance = distance;
835           cl=c;
836           //      index=krow.GetIndex(i);       
837         }       
838       }              
839     }  // end of seearch
840     //update cluster information
841     if (cl){ 
842       cy[drow] = cl->GetY();
843       cz[drow] = cl->GetZ();
844       sy[drow] = ErrY2(&t,cl);
845       sz[drow] = ErrZ2(&t,cl);
846       hitted[drow] = kTRUE;
847       nclusters++;
848     }
849     else
850       hitted[drow] = kFALSE;
851   }
852   //if we have information - update track
853   if (nclusters>0){
854     Float_t sumyw0   = 0;
855     Float_t sumzw0   = 0;
856     Float_t sumyw   = 0;
857     Float_t sumzw   = 0;
858     Float_t sumrow  = 0;
859     for (Int_t i=0;i<step;i++)
860       if (hitted[i]){
861         sumyw0+= 1/sy[i];
862         sumzw0+= 1/sz[i];
863         //
864         sumyw+= (cy[i]-yy[i])/sy[i];
865         sumzw+= (cz[i]-zz[i])/sz[i];
866         sumrow+= i;
867       }
868     Float_t dy = sumyw/sumyw0;
869     Float_t dz = sumzw/sumzw0;
870     Float_t mrow = sumrow/nclusters+row0;
871     Float_t x = fSectors->GetX(mrow);
872     t.PropagateTo(x);
873     AliTPCclusterMI cvirtual;
874     cvirtual.SetZ(dz+t.GetZ());
875     cvirtual.SetY(dy+t.GetY());
876     t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
877     t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
878     Float_t chi2 = t.GetPredictedChi2(&cvirtual);
879     t.Update(&cvirtual,chi2,0);
880     Int_t ncl = t.GetNumberOfClusters();
881     ncl = ncl-1+nclusters;
882     t.SetN(ncl);
883   }     
884   return  1;
885 }   
886 */
887
888 //_____________________________________________________________________________
889 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
890   //-----------------------------------------------------------------
891   // This function tries to find a track prolongation.
892   //-----------------------------------------------------------------
893   Double_t xt=t.GetX();
894   //
895   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
896   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
897   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
898   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
899     
900   for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
901    
902     if (FollowToNext(t,nr)==0) {
903     }
904   }   
905   return 1;
906 }
907
908
909 //_____________________________________________________________________________
910 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
911   //-----------------------------------------------------------------
912   // This function tries to find a track prolongation.
913   //-----------------------------------------------------------------
914   Double_t xt=t.GetX();
915   //
916   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
917   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
918   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
919   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
920     
921   for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
922     if (t.GetSnp()>0.2)
923       FollowToNext(t,nr);                                                             
924   }   
925   return 1;
926 }
927
928
929
930
931    
932 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
933 {
934   //
935   //
936   sum1=0;
937   sum2=0;
938   Int_t sum=0;
939   if (s1->fSector!=s2->fSector) return 0;
940   //
941   Float_t dz2 =(s1->GetZ() - s2->GetZ());
942   dz2*=dz2;
943   Float_t dy2 =(s1->GetY() - s2->GetY());
944   dy2*=dy2;
945   Float_t distance = TMath::Sqrt(dz2+dy2);
946   if (distance>5.) return 0; // if there are far away  - not overlap - to reduce combinatorics
947  
948   for (Int_t i=0;i<160;i++){
949     if (s1->fClusterIndex[i]>0) sum1++;
950     if (s2->fClusterIndex[i]>0) sum2++;
951     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
952       sum++;
953     }
954   }
955  
956   Float_t summin = TMath::Min(sum1,sum2);
957   Float_t ratio = sum/Float_t(summin);
958   return ratio;
959 }
960
961 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
962 {
963   //
964   //
965   if (s1->fSector!=s2->fSector) return;
966   //
967   Float_t dz2 =(s1->GetZ() - s2->GetZ());
968   dz2*=dz2;
969   Float_t dy2 =(s1->GetY() - s2->GetY());
970   dy2*=dy2;
971   Float_t distance = TMath::Sqrt(dz2+dy2);
972   if (distance>15.) return ; // if there are far away  - not overlap - to reduce combinatorics
973   //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
974   //  TClonesArray &pointarray1 = *(s1->fPoints);
975   //TClonesArray &pointarray2 = *(s2->fPoints);
976   //
977   for (Int_t i=0;i<160;i++){
978     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
979       //  AliTPCTrackPoint *p1  = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
980       //AliTPCTrackPoint *p2  = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i)); 
981       AliTPCTrackPoint *p1  = s1->GetTrackPoint(i);
982       AliTPCTrackPoint *p2  = s2->GetTrackPoint(i);; 
983       p1->fIsShared = kTRUE;
984       p2->fIsShared = kTRUE;
985     }
986   } 
987 }
988
989
990
991
992 void  AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared=kFALSE){
993
994   
995
996   //
997   // remove overlap - used removal factor - removal index stored in the track
998   arr->Sort();  // sorting according z
999   arr->Expand(arr->GetEntries());
1000   Int_t nseed=arr->GetEntriesFast();
1001   //  printf("seeds \t%p \t%d\n",arr, nseed);
1002   //  arr->Expand(arr->GetEntries());  //remove 0 pointers
1003   nseed = arr->GetEntriesFast();
1004   Int_t removed = 0;
1005   for (Int_t i=0; i<nseed; i++) {
1006     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
1007     if (!pt) {
1008       continue;
1009     }
1010     if (!(pt->IsActive())) continue;
1011     for (Int_t j=i+1; j<nseed; j++){
1012       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1013       if ((pt2) && pt2->IsActive())
1014         if (pt->fSector == pt2->fSector)
1015           if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
1016             Int_t sum1,sum2;
1017             Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
1018             if (ratio>factor){
1019               //          if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex);  // arr->RemoveAt(j);             
1020               Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
1021               Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
1022               removed++;
1023               if (TMath::Abs(ratio3)>0.025){  // if much more points  
1024                 if (sum1>sum2) pt2->Desactivate(removalindex);
1025                 else {
1026                   pt->Desactivate(removalindex); // arr->RemoveAt(i); 
1027                   break;
1028                 }
1029               }
1030               else{  //decide on mean chi2
1031                 if (ratio2<1)  
1032                   pt2->Desactivate(removalindex);
1033                 else {
1034                   pt->Desactivate(removalindex); // arr->RemoveAt(i); 
1035                   break;
1036                 }           
1037               }  
1038               
1039             }  // if suspicious ratio
1040           }
1041           else
1042             break;
1043     }
1044   }
1045   //  printf("removed\t%d\n",removed);
1046   Int_t good =0; 
1047   for (Int_t i=0; i<nseed; i++) {
1048     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
1049     if (!pt) break;
1050     if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
1051       //desactivate tracks with small number of points
1052       //      printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
1053       pt->Desactivate(10);  //desactivate  - small muber of points
1054     }
1055     if (!(pt->IsActive())) continue;
1056     good++;
1057   }
1058   
1059   
1060   if (shared)
1061     for (Int_t i=0; i<nseed; i++) {
1062       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
1063       if (!pt) continue;
1064       if (!(pt->IsActive())) continue;
1065       for (Int_t j=i+1; j<nseed; j++){
1066         AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1067         if ((pt2) && pt2->IsActive()) {
1068           if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
1069           SignShared(pt,pt2);
1070         }
1071       }
1072     }
1073   fNtracks = good;
1074   printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
1075
1076
1077 }
1078 void AliTPCtrackerMI::MakeSeedsAll()
1079 {
1080   if (fSeeds == 0) fSeeds = new TObjArray;
1081   TObjArray * arr;
1082   for (Int_t sec=0;sec<fkNOS;sec+=3){
1083      arr = MakeSeedsSectors(sec,sec+3);
1084      Int_t nseed = arr->GetEntriesFast();
1085      for (Int_t i=0;i<nseed;i++)
1086        fSeeds->AddLast(arr->RemoveAt(i));
1087   }
1088    //  fSeeds = MakeSeedsSectors(0,fkNOS);
1089 }
1090
1091 TObjArray *  AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
1092 {
1093   //
1094   // loop over all  sectors and make seed
1095   //find track seeds
1096   TStopwatch timer;
1097   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1098   Int_t nrows=nlow+nup;  
1099   Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
1100   //  if (fSeeds==0) fSeeds = new TObjArray;
1101   TObjArray * arr = new TObjArray;
1102
1103   for (Int_t sec=sec1; sec<sec2;sec++){
1104     MakeSeeds(arr, sec, nup-1, nup-1-gap);
1105     MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
1106   }
1107   Int_t nseed=arr->GetEntriesFast();
1108   gap=Int_t(0.3*nrows);
1109   // continue seeds 
1110   Int_t i;
1111   for (i=0; i<nseed; i++) {
1112     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt; 
1113     if (!pt) continue;
1114     if (FollowProlongation(t,nup-gap)) {
1115       pt->fIsSeeding =kFALSE;
1116       continue;
1117     }
1118     delete arr->RemoveAt(i);
1119   }   
1120   //
1121   //remove seeds which overlaps  
1122   RemoveOverlap(arr,0.6,1);       
1123   //delete seeds - which were sign  
1124   nseed=arr->GetEntriesFast();
1125   for (i=0; i<nseed; i++) {
1126     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1127     //, &t=*pt;
1128     if (!pt) continue;
1129     if ((pt->IsActive())  &&  pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
1130       //FollowBackProlongation(t,nup);
1131       continue;
1132     }
1133     delete arr->RemoveAt(i);    
1134   }  
1135   return arr;
1136 }
1137
1138
1139
1140 //_____________________________________________________________________________
1141 void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1142   //-----------------------------------------------------------------
1143   // This function creates track seeds.
1144   //-----------------------------------------------------------------
1145   //  if (fSeeds==0) fSeeds=new TObjArray(15000);
1146
1147   Double_t x[5], c[15];
1148
1149   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1150   Double_t cs=cos(alpha), sn=sin(alpha);
1151
1152   Double_t x1 =fOuterSec->GetX(i1);
1153   Double_t xx2=fOuterSec->GetX(i2);
1154   //
1155   //  for (Int_t ns=0; ns<fkNOS; ns++) 
1156   Int_t ns =sec;
1157     {
1158     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1159     Int_t nm=fOuterSec[ns][i2];
1160     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
1161     const AliTPCRow& kr1=fOuterSec[ns][i1];
1162     AliTPCRow&  kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1163     AliTPCRow&  kr22 = fOuterSec[(ns)%fkNOS][i2];
1164     AliTPCRow&  kr23 = fOuterSec[(ns+1)%fkNOS][i2];
1165
1166     for (Int_t is=0; is < kr1; is++) {
1167       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1168       Double_t x3=GetX(), y3=GetY(), z3=GetZ();
1169
1170       Float_t anglez = (z1-z3)/(x1-x3); 
1171       Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
1172
1173       for (Int_t js=0; js < nl+nm+nu; js++) {
1174         const AliTPCclusterMI *kcl;
1175         Double_t x2,   y2,   z2;
1176         if (js<nl) {     
1177           if (js==0) {
1178             js = kr21.Find(extraz-15.);
1179             if (js>=nl) continue;
1180           }       
1181           kcl=kr21[js];
1182           z2=kcl->GetZ();
1183           if ((extraz-z2)>10) continue;   
1184           if ((extraz-z2)<-10) {
1185             js = nl;
1186             continue;
1187           }
1188           y2=kcl->GetY(); 
1189           x2= xx2*cs+y2*sn;
1190           y2=-xx2*sn+y2*cs;
1191         } else 
1192           if (js<nl+nm) {
1193             if (js==nl) {
1194               js = nl+kr22.Find(extraz-15.);
1195               if (js>=nl+nm) continue;
1196             }             
1197             kcl=kr22[js-nl];
1198             z2=kcl->GetZ();
1199             if ((extraz-z2)>10) continue;
1200             if ((extraz-z2)<-10) {
1201               js = nl+nm;
1202               continue;
1203             }
1204             x2=xx2; y2=kcl->GetY(); 
1205           } else {
1206             //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];   
1207             if (js==nl+nm) {
1208               js = nl+nm+kr23.Find(extraz-15.);
1209               if (js>=nl+nm+nu) break;
1210             }     
1211             kcl=kr23[js-nl-nm];
1212             z2=kcl->GetZ(); 
1213             if ((extraz-z2)>10) continue;
1214             if ((extraz-z2)<-10) {
1215               break;        
1216             }
1217             y2=kcl->GetY();
1218             x2=xx2*cs-y2*sn;
1219             y2=xx2*sn+y2*cs;
1220           }
1221
1222         Double_t zz=z1 - anglez*(x1-x2); 
1223         if (TMath::Abs(zz-z2)>10.) continue;
1224
1225         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1226         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
1227
1228         x[0]=y1;
1229         x[1]=z1;
1230         x[4]=f1(x1,y1,x2,y2,x3,y3);
1231         if (TMath::Abs(x[4]) >= 0.0066) continue;
1232         x[2]=f2(x1,y1,x2,y2,x3,y3);
1233         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1234         x[3]=f3(x1,y1,x2,y2,z1,z2);
1235         if (TMath::Abs(x[3]) > 1.2) continue;
1236         Double_t a=asin(x[2]);
1237         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1238         if (TMath::Abs(zv-z3)>10.) continue; 
1239
1240         Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
1241         Double_t sy2=kcl->GetSigmaY2()*2,     sz2=kcl->GetSigmaZ2()*4;
1242         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1243         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1244         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1245
1246         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1247         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1248         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1249         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1250         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1251         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1252         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1253         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1254         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1255         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1256
1257         c[0]=sy1;
1258         c[1]=0.;       c[2]=sz1;
1259         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1260         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1261                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1262         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1263         c[13]=f30*sy1*f40+f32*sy2*f42;
1264         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1265
1266         UInt_t index=kr1.GetIndex(is);
1267         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
1268         track->fIsSeeding = kTRUE;
1269         Int_t rc=FollowProlongation(*track, i2);
1270         //FollowProlongationFast(*track, 5);
1271         //FollowProlongationFast(*track, 5);
1272         //FollowProlongationFast(*track, 5);
1273         //FollowProlongationFast(*track, 5);
1274         //Int_t rc = 1;
1275         
1276         track->fLastPoint = i1;  // first cluster in track position
1277         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1278         else arr->AddLast(track); 
1279       }
1280     }
1281   }
1282 }
1283
1284 //_____________________________________________________________________________
1285 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
1286   //-----------------------------------------------------------------
1287   // This function reades track seeds.
1288   //-----------------------------------------------------------------
1289   TDirectory *savedir=gDirectory; 
1290
1291   TFile *in=(TFile*)inp;
1292   if (!in->IsOpen()) {
1293      cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
1294      return 1;
1295   }
1296
1297   in->cd();
1298   TTree *seedTree=(TTree*)in->Get("Seeds");
1299   if (!seedTree) {
1300      cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
1301      cerr<<"can't get a tree with track seeds !\n";
1302      return 2;
1303   }
1304   AliTPCtrack *seed=new AliTPCtrack; 
1305   seedTree->SetBranchAddress("tracks",&seed);
1306   
1307   if (fSeeds==0) fSeeds=new TObjArray(15000);
1308
1309   Int_t n=(Int_t)seedTree->GetEntries();
1310   for (Int_t i=0; i<n; i++) {
1311      seedTree->GetEvent(i);
1312      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
1313   }
1314   
1315   delete seed;
1316   delete seedTree; 
1317   savedir->cd();
1318   return 0;
1319 }
1320
1321 //_____________________________________________________________________________
1322 Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
1323   //-----------------------------------------------------------------
1324   // This is a track finder.
1325   //-----------------------------------------------------------------
1326   TDirectory *savedir=gDirectory; 
1327
1328   if (inp) {
1329      TFile *in=(TFile*)inp;
1330      if (!in->IsOpen()) {
1331         cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
1332         return 1;
1333      }
1334   }
1335
1336   if (!out->IsOpen()) {
1337      cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
1338      return 2;
1339   }
1340
1341   out->cd();
1342
1343   char   tname[100];
1344   sprintf(tname,"TreeT_TPC_%d",fEventN);
1345   TTree tracktree(tname,"Tree with TPC tracks");
1346   TTree seedtree("Seeds","Seeds");
1347   AliTPCtrack *iotrack=0;
1348   AliTPCseed  *ioseed=0;
1349   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
1350   TStopwatch timer;
1351   
1352   printf("Loading clusters \n");
1353   LoadClusters();
1354   printf("Time for loading clusters: \t");timer.Print();timer.Start();
1355
1356   printf("Loading outer sectors\n");
1357   LoadOuterSectors();
1358   printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
1359
1360   printf("Loading inner sectors\n");
1361   LoadInnerSectors();
1362   printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
1363   fSectors = fOuterSec;
1364   fN=fkNOS;
1365   
1366   
1367   //find track seeds
1368   MakeSeedsAll(); 
1369   printf("Time for seeding: \t"); timer.Print();timer.Start();
1370   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1371   Int_t nrows=nlow+nup;
1372   
1373   Int_t gap=Int_t(0.3*nrows);
1374   Int_t i;
1375   Int_t nseed=fSeeds->GetEntriesFast();
1376   // outer sectors parallel tracking
1377   ParallelTracking(fSectors->GetNRows()-gap-1,0); 
1378   printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
1379
1380   RemoveOverlap(fSeeds, 0.6,3);  
1381   printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
1382   //parallel tracking 
1383   fSectors = fInnerSec;
1384   fN=fkNIS;
1385   ParallelTracking(fSectors->GetNRows()-1,0);
1386   printf("Number of tracks after  inner tracking  %d\n",fNtracks); 
1387   printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
1388   //
1389   RemoveOverlap(fSeeds,0.6,5,kTRUE);  // remove overlap -  shared points signed 
1390   printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
1391   //
1392   // 
1393   
1394   ioseed  = (AliTPCseed*)(fSeeds->UncheckedAt(0));
1395   AliTPCseed * vseed = new AliTPCseed;
1396   vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
1397   vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
1398   vseed->fPoints->ExpandCreateFast(2);
1399   
1400   TBranch * seedbranch =   seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
1401   //delete vseed;
1402   nseed=fSeeds->GetEntriesFast();
1403
1404   Int_t found = 0;
1405   for (i=0; i<nseed; i++) {
1406     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;    
1407     if (!pt) continue;    
1408     Int_t nc=t.GetNumberOfClusters();
1409     t.CookdEdx(0.02,0.6);
1410     CookLabel(pt,0.1); //For comparison only
1411     //   if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
1412     if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (nc>Int_t(0.3*nrows)))){
1413       iotrack=pt;
1414       tracktree.Fill();
1415      cerr<<found++<<'\r';      
1416     }   
1417     else 
1418       if ( (pt->IsActive())) fNtracks--;
1419     pt->RebuildSeed();
1420     seedbranch->SetAddress(&pt);
1421
1422     seedtree.Fill();        
1423     for (Int_t j=0;j<160;j++){
1424       delete pt->fPoints->RemoveAt(j);
1425     }
1426     delete pt->fPoints;
1427     pt->fPoints =0;
1428     delete fSeeds->RemoveAt(i);
1429   }
1430   printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
1431
1432   UnloadClusters();
1433   printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1434
1435   tracktree.Write();
1436   seedtree.Write();
1437   cerr<<"Number of found tracks : "<<fNtracks<<"\t"<<found<<endl;
1438   
1439   savedir->cd();
1440   
1441   return 0;
1442 }
1443
1444
1445 void  AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
1446 {
1447   //
1448   // try to track in parralel
1449
1450   Int_t nseed=fSeeds->GetEntriesFast();
1451   //prepare seeds for tracking
1452   for (Int_t i=0; i<nseed; i++) {
1453     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt; 
1454     if (!pt) continue;
1455     if (!t.IsActive()) continue;
1456     // follow prolongation to the first layer
1457     FollowProlongation(t, rfirst+1);
1458   }
1459
1460
1461   //
1462   for (Int_t nr=rfirst; nr>=rlast; nr--){      
1463     // make indexes with the cluster tracks for given       
1464     //    for (Int_t i = 0;i<fN;i++)
1465     //  fSectors[i][nr].MakeClusterTracks();
1466
1467     // find nearest cluster
1468     for (Int_t i=0; i<nseed; i++) {
1469       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt; 
1470       if (!pt) continue;
1471       if (!pt->IsActive()) continue;
1472       UpdateClusters(t,i,nr);
1473     }
1474     // prolonagate to the nearest cluster - if founded
1475     for (Int_t i=0; i<nseed; i++) {
1476       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i); 
1477       if (!pt) continue;
1478       if (!pt->IsActive()) continue;
1479       FollowToNextCluster(i,nr);
1480     }
1481     //    for (Int_t i= 0;i<fN;i++)
1482     //  fSectors[i][nr].ClearClusterTracks();
1483   }  
1484   
1485 }
1486
1487 Float_t  AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
1488 {
1489   //
1490   //  
1491   Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1492   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1493   Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
1494   Float_t angular  = seed->GetSnp();
1495   angular = angular*angular/(1-angular*angular);
1496   //  angular*=angular;
1497   //angular  = TMath::Sqrt(angular/(1-angular));
1498   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
1499   return res;
1500 }
1501 Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
1502 {
1503   //
1504   //
1505   Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1506   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1507   Float_t sres = fParam->GetZSigma();
1508   Float_t angular  = seed->GetTgl();
1509   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
1510   return res;
1511 }
1512
1513
1514
1515 //_________________________________________________________________________
1516 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1517   //--------------------------------------------------------------------
1518   //       Return pointer to a given cluster
1519   //--------------------------------------------------------------------
1520   Int_t sec=(index&0xff000000)>>24; 
1521   Int_t row=(index&0x00ff0000)>>16; 
1522   Int_t ncl=(index&0x0000ffff)>>00;
1523
1524   AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
1525   return (AliTPCclusterMI*)(*clrow)[ncl];      
1526 }
1527
1528 //__________________________________________________________________________
1529 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1530   //--------------------------------------------------------------------
1531   //This function "cooks" a track label. If label<0, this track is fake.
1532   //--------------------------------------------------------------------
1533   Int_t noc=t->GetNumberOfClusters();
1534   Int_t *lb=new Int_t[noc];
1535   Int_t *mx=new Int_t[noc];
1536   AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
1537
1538   Int_t i;
1539   for (i=0; i<noc; i++) {
1540      lb[i]=mx[i]=0;
1541      Int_t index=t->GetClusterIndex(i);
1542      clusters[i]=GetClusterMI(index);
1543   }
1544
1545   Int_t lab=123456789;
1546   for (i=0; i<noc; i++) {
1547     AliTPCclusterMI *c=clusters[i];
1548     lab=TMath::Abs(c->GetLabel(0));
1549     Int_t j;
1550     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1551     lb[j]=lab;
1552     (mx[j])++;
1553   }
1554
1555   Int_t max=0;
1556   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1557     
1558   for (i=0; i<noc; i++) {
1559     AliTPCclusterMI *c=clusters[i];
1560     if (TMath::Abs(c->GetLabel(1)) == lab ||
1561         TMath::Abs(c->GetLabel(2)) == lab ) max++;
1562   }
1563
1564   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1565
1566   else {
1567      Int_t tail=Int_t(0.10*noc);
1568      max=0;
1569      for (i=1; i<=tail; i++) {
1570        AliTPCclusterMI *c=clusters[noc-i];
1571        if (lab == TMath::Abs(c->GetLabel(0)) ||
1572            lab == TMath::Abs(c->GetLabel(1)) ||
1573            lab == TMath::Abs(c->GetLabel(2))) max++;
1574      }
1575      if (max < Int_t(0.5*tail)) lab=-lab;
1576   }
1577
1578   t->SetLabel(lab);
1579
1580   delete[] lb;
1581   delete[] mx;
1582   delete[] clusters;
1583 }
1584
1585 //_________________________________________________________________________
1586 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1587   //-----------------------------------------------------------------------
1588   // Setup inner sector
1589   //-----------------------------------------------------------------------
1590   if (f==0) {
1591      fAlpha=par->GetInnerAngle();
1592      fAlphaShift=par->GetInnerAngleShift();
1593      fPadPitchWidth=par->GetInnerPadPitchWidth();
1594      fPadPitchLength=par->GetInnerPadPitchLength();
1595      fN=par->GetNRowLow();
1596      fRow=new AliTPCRow[fN];
1597      for (Int_t i=0; i<fN; i++) {
1598        fRow[i].SetX(par->GetPadRowRadiiLow(i));
1599        fRow[i].fDeadZone =1.5;  //1.5 cm of dead zone
1600      }
1601   } else {
1602      fAlpha=par->GetOuterAngle();
1603      fAlphaShift=par->GetOuterAngleShift();
1604      fPadPitchWidth  = par->GetOuterPadPitchWidth();
1605      fPadPitchLength = par->GetOuter1PadPitchLength();
1606      f1PadPitchLength = par->GetOuter1PadPitchLength();
1607      f2PadPitchLength = par->GetOuter2PadPitchLength();
1608
1609      fN=par->GetNRowUp();
1610      fRow=new AliTPCRow[fN];
1611      for (Int_t i=0; i<fN; i++) {
1612        fRow[i].SetX(par->GetPadRowRadiiUp(i)); 
1613        fRow[i].fDeadZone =1.5;  // 1.5 cm of dead zone
1614      }
1615   } 
1616 }
1617
1618
1619 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
1620   //
1621   if (fClusterTracks) delete [] fClusterTracks;
1622   fClusterTracks = 0;
1623 }
1624
1625 void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
1626   //create cluster tracks
1627   if (fN>0) 
1628     fClusterTracks = new AliTPCclusterTracks[fN];
1629 }
1630
1631 void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
1632   if (fClusterTracks) delete[] fClusterTracks;
1633   fClusterTracks =0;
1634 }
1635
1636
1637
1638 void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
1639   //
1640   //
1641   // update information of the cluster tracks - if track is nearer then other tracks to the 
1642   // given track
1643   const AliTPCclusterMI * cl = (*this)[clindex];
1644   AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
1645   // find the distance of the cluster to the track
1646   Float_t dy2 = (cl->GetY()- seed->GetY());
1647   dy2*=dy2;
1648   Float_t dz2 = (cl->GetZ()- seed->GetZ());
1649   dz2*=dz2;
1650   //
1651   Float_t distance = TMath::Sqrt(dy2+dz2);
1652   if (distance> 3) 
1653     return;  // MI - to be changed - AliTPCtrackerParam
1654   
1655   if ( distance < cltracks->fDistance[0]){
1656     cltracks->fDistance[2] =cltracks->fDistance[1];
1657     cltracks->fDistance[1] =cltracks->fDistance[0];
1658     cltracks->fDistance[0] =distance;
1659     cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1660     cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
1661     cltracks->fTrackIndex[0] =trindex; 
1662   }
1663   else
1664     if ( distance < cltracks->fDistance[1]){
1665       cltracks->fDistance[2] =cltracks->fDistance[1];  
1666       cltracks->fDistance[1] =distance;
1667       cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1668       cltracks->fTrackIndex[1] =trindex; 
1669     } else
1670       if (distance < cltracks->fDistance[2]){
1671         cltracks->fDistance[2] =distance;
1672         cltracks->fTrackIndex[2] =trindex;
1673       }  
1674 }
1675
1676
1677 //_________________________________________________________________________
1678 void 
1679 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
1680   //-----------------------------------------------------------------------
1681   // Insert a cluster into this pad row in accordence with its y-coordinate
1682   //-----------------------------------------------------------------------
1683   if (fN==kMaxClusterPerRow) {
1684     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1685   }
1686   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
1687   Int_t i=Find(c->GetZ());
1688   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
1689   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
1690   fIndex[i]=index; fClusters[i]=c; fN++;
1691 }
1692
1693 //___________________________________________________________________
1694 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
1695   //-----------------------------------------------------------------------
1696   // Return the index of the nearest cluster 
1697   //-----------------------------------------------------------------------
1698   if (fN==0) return 0;
1699   if (z <= fClusters[0]->GetZ()) return 0;
1700   if (z > fClusters[fN-1]->GetZ()) return fN;
1701   Int_t b=0, e=fN-1, m=(b+e)/2;
1702   for (; b<e; m=(b+e)/2) {
1703     if (z > fClusters[m]->GetZ()) b=m+1;
1704     else e=m; 
1705   }
1706   return m;
1707 }
1708
1709
1710
1711 AliTPCseed::AliTPCseed():AliTPCtrack(){
1712   //
1713   fRow=0; 
1714   fRemoval =0; 
1715   memset(fClusterIndex,0,sizeof(Int_t)*200);
1716   fPoints = 0;
1717   fEPoints = 0;
1718   fNFoundable =0;
1719   fNShared  =0;
1720   fTrackPoints =0;
1721   fRemoval = 0;
1722 }
1723
1724 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
1725   fPoints = 0;
1726   fEPoints = 0;
1727   fNShared  =0; 
1728   fTrackPoints =0;
1729   fRemoval =0;
1730 }
1731
1732 AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
1733   fRow=0;
1734   memset(fClusterIndex,0,sizeof(Int_t)*200); 
1735   fPoints = 0;
1736   fEPoints = 0;
1737   fNFoundable =0; 
1738   fNShared  =0; 
1739   fTrackPoints =0;
1740   fRemoval =0;
1741 }
1742
1743 AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15], 
1744                                         Double_t xr, Double_t alpha):      
1745   AliTPCtrack(index, xx, cc, xr, alpha) {
1746   //
1747   //
1748   fRow =0;
1749   memset(fClusterIndex,0,sizeof(Int_t)*200); 
1750   fPoints = 0;
1751   fEPoints = 0;
1752   fNFoundable =0;
1753   fNShared  = 0;
1754   fTrackPoints =0;
1755   fRemoval =0;
1756 }
1757
1758 AliTPCseed::~AliTPCseed(){
1759   if (fPoints) delete fPoints;
1760   fPoints =0;
1761   fEPoints = 0;
1762   if (fTrackPoints){
1763     for (Int_t i=0;i<8;i++){
1764       delete [] fTrackPoints[i];
1765     }
1766     delete fTrackPoints;
1767     fTrackPoints =0;
1768   }
1769
1770 }
1771
1772 AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
1773 {
1774   //
1775   // 
1776   if (!fTrackPoints) {
1777     fTrackPoints = new AliTPCTrackPoint*[8];
1778     for ( Int_t i=0;i<8;i++)
1779       fTrackPoints[i]=0;
1780   }
1781   Int_t index1 = i/20;
1782   if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
1783   return &(fTrackPoints[index1][i%20]);
1784 }
1785
1786 void AliTPCseed::RebuildSeed()
1787 {
1788   //
1789   // rebuild seed to be ready for storing
1790   fPoints = new TClonesArray("AliTPCTrackPoint",160);
1791   fPoints->ExpandCreateFast(160);
1792   fEPoints = new TClonesArray("AliTPCExactPoint",1);
1793   for (Int_t i=0;i<160;i++){
1794     AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
1795     *trpoint = *(GetTrackPoint(i));
1796   }
1797
1798 }
1799
1800 //_____________________________________________________________________________
1801 void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
1802   //-----------------------------------------------------------------
1803   // This funtion calculates dE/dX within the "low" and "up" cuts.
1804   //-----------------------------------------------------------------
1805
1806   Float_t amp[200];
1807   Float_t angular[200];
1808   Float_t weight[200];
1809   Int_t index[200];
1810   //Int_t nc = 0;
1811   //  TClonesArray & arr = *fPoints; 
1812   Float_t meanlog = 100.;
1813   
1814   Float_t mean[4]  = {0,0,0,0};
1815   Float_t sigma[4] = {1000,1000,1000,1000};
1816   Int_t nc[4]      = {0,0,0,0};
1817   Float_t norm[4]    = {1000,1000,1000,1000};
1818   //
1819   //
1820   fNShared =0;
1821
1822   for (Int_t of =0; of<4; of++){    
1823     for (Int_t i=of;i<160;i+=4)
1824       {
1825         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
1826         AliTPCTrackPoint * point = GetTrackPoint(i);
1827         if (point==0) continue;
1828         if (point->fIsShared){
1829           fNShared++;
1830           continue;
1831         }
1832         if (point->GetCPoint().GetMax()<5) continue;
1833         Float_t angley = point->GetTPoint().GetAngleY();
1834         Float_t anglez = point->GetTPoint().GetAngleZ();
1835         Int_t   type   = point->GetCPoint().GetType();
1836         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
1837         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
1838         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
1839
1840         Float_t ampc   = 0;     // normalization to the number of electrons
1841         if (i>64){
1842           ampc = 1.*point->GetCPoint().GetMax();
1843           //ampc = 1.*point->GetCPoint().GetQ();          
1844           //      AliTPCClusterPoint & p = point->GetCPoint();
1845           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
1846           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
1847           //Float_t dz = 
1848           //  TMath::Abs( Int_t(iz) - iz + 0.5);
1849           //ampc *= 1.15*(1-0.3*dy);
1850           //ampc *= 1.15*(1-0.3*dz);
1851
1852         }
1853         else{ 
1854           ampc = 1.0*point->GetCPoint().GetMax(); 
1855           //ampc = 1.0*point->GetCPoint().GetQ(); 
1856           //AliTPCClusterPoint & p = point->GetCPoint();
1857           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
1858           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
1859           //Float_t dz = 
1860           //  TMath::Abs( Int_t(iz) - iz + 0.5);
1861
1862           //ampc *= 1.15*(1-0.3*dy);
1863           //ampc *= 1.15*(1-0.3*dz);
1864
1865         }
1866         ampc *= 2.0;     // put mean value to channel 50
1867         //ampc *= 0.58;     // put mean value to channel 50
1868         Float_t w      =  1.;
1869         //      if (type>0)  w =  1./(type/2.-0.5); 
1870         if (i<64) {
1871           ampc /= 0.6;
1872         }
1873         if (i>128){
1874           ampc /=1.5;
1875         }
1876         if (type<0) {  //amp at the border - lower weight
1877           // w*= 2.;
1878           continue;
1879         }
1880         if (rsigma>1.5) ampc/=1.3;  // if big backround
1881         amp[nc[of]]        = ampc;
1882         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
1883         weight[nc[of]]     = w;
1884         nc[of]++;
1885       }
1886     
1887     TMath::Sort(nc[of],amp,index,kFALSE);
1888     Float_t sumamp=0;
1889     Float_t sumamp2=0;
1890     Float_t sumw=0;
1891     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
1892     meanlog = 200;
1893     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
1894       Float_t ampl      = amp[index[i]]/angular[index[i]];
1895       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
1896       //
1897       sumw    += weight[index[i]]; 
1898       sumamp  += weight[index[i]]*ampl;
1899       sumamp2 += weight[index[i]]*ampl*ampl;
1900       norm[of]    += angular[index[i]]*weight[index[i]];
1901     }
1902     if (sumw<1){ 
1903       SetdEdx(0);  
1904     }
1905     else {
1906       norm[of] /= sumw;
1907       mean[of]  = sumamp/sumw;
1908       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
1909       if (sigma[of]>0.1) 
1910         sigma[of] = TMath::Sqrt(sigma[of]);
1911       else
1912         sigma[of] = 1000;
1913       
1914     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
1915     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
1916     //mean *=(1-0.1*(norm-1.));
1917     }
1918   }
1919
1920   Float_t dedx =0;
1921   fSdEdx =0;
1922   fMAngular =0;
1923   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
1924   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
1925
1926   
1927   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ 
1928   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
1929
1930   Int_t norm2 = 0;
1931   Int_t norm3 = 0;
1932   for (Int_t i =0;i<4;i++){
1933     if (nc[i]>2&&nc[i]<1000){
1934       dedx      += mean[i] *nc[i];
1935       fSdEdx    += sigma[i]*(nc[i]-2);
1936       fMAngular += norm[i] *nc[i];    
1937       norm2     += nc[i];
1938       norm3     += nc[i]-2;
1939     }
1940     fDEDX[i]  = mean[i];             
1941     fSDEDX[i] = sigma[i];            
1942     fNCDEDX[i]= nc[i]; 
1943   }
1944
1945   if (norm3>0){
1946     dedx   /=norm2;
1947     fSdEdx /=norm3;
1948     fMAngular/=norm2;
1949   }
1950   else{
1951     SetdEdx(0);
1952     return;
1953   }
1954   //  Float_t dedx1 =dedx;
1955   
1956   dedx =0;
1957   for (Int_t i =0;i<4;i++){
1958     if (nc[i]>2&&nc[i]<1000){
1959       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
1960       dedx      += mean[i] *nc[i];
1961     }
1962     fDEDX[i]  = mean[i];                
1963   }
1964   dedx /= norm2;
1965   
1966
1967   
1968   SetdEdx(dedx);
1969     
1970   //mi deDX
1971
1972
1973
1974   //Very rough PID
1975   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
1976
1977   if (p<0.6) {
1978     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1979     if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
1980     SetMass(0.93827); return;
1981   }
1982
1983   if (p<1.2) {
1984     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
1985     SetMass(0.93827); return;
1986   }
1987
1988   SetMass(0.13957); return;
1989
1990 }
1991
1992