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