from Marco :
[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 //          Implementation of the TPC tracker
19 //
20 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
21 // 
22 //  AliTPC parallel tracker
23 //
24 //  The track fitting is based on Kalaman filtering approach
25
26 //  The track finding steps:
27 //      1. Seeding - with and without vertex constraint
28 //                 - seeding with vertex constain done at first n^2 proble
29 //                 - seeding without vertex constraint n^3 problem
30 //      2. Tracking - follow prolongation road - find cluster - update kalman track
31
32 //  The seeding and tracking is repeated several times, in different seeding region.
33 //  This approach enables to find the track which cannot be seeded in some region of TPC
34 //  This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
35
36 //  With this approach we reach almost 100 % efficiency also for high occupancy events.
37 //  (If the seeding efficiency in a region is about 90 % than with logical or of several 
38 //  regions we will reach 100% (in theory - supposing independence) 
39
40 //  Repeating several seeding - tracking procedures some of the tracks can be find 
41 //  several times. 
42
43 //  The procedures to remove multi find tacks are impremented:
44 //  RemoveUsed2 - fast procedure n problem - 
45 //                 Algorithm - Sorting tracks according quality
46 //                             remove tracks with some shared fraction 
47 //                             Sharing in respect to all tacks 
48 //                             Signing clusters in gold region
49 //  FindSplitted - slower algorithm n^2
50 //                 Sort the tracks according quality
51 //                 Loop over pair of tracks
52 //                 If overlap with other track bigger than threshold - remove track
53 //  
54 //  FindCurling  - Finds the pair of tracks which are curling
55 //               - About 10% of tracks can be find with this procedure
56 //                 The combinatorial background is too big to be used in High 
57 //                  multiplicity environment 
58 //               - n^2 problem - Slow procedure - currently it is disabled because of 
59 //                  low efficiency
60 //                 
61 //  The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 //  tpcRecoParam-> SetClusterSharing(kFALSE);
63 //  IT IS HIGHLY non recomended to use it in high flux enviroonment
64 //  Even using this switch some tracks can be found more than once 
65 //  (because of multiple seeding and low quality tracks which will not cross full chamber)
66 //                          
67 //
68 // The tracker itself can be debugged  - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); 
71 //
72 // The debug level -  different procedure produce tree for numerical debugging
73 //                    To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
74 //
75
76 //
77 // Adding systematic errors to the covariance:
78 // 
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0. 
84 //
85 // The sytematic errors are added to the covariance matrix in following places:
86 //
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration - 
89 //      2.a ITS->TPC   - AliTPCtrackerMI::ReadSeeds 
90 //      2.b TPC->TRD   - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration  -
92 //      3.a TRD->TPC   - AliTPCtrackerMI::ReadSeeds
93 //      3.b TPC->ITS   - AliTPCtrackerMI::RefitInward
94 //
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
97 //
98 //      1. ErrParam stream  - dump information about 
99 //         1.a) cluster
100 //         2.a) cluster error estimate
101 //         3.a) cluster shape estimate
102 //
103 //
104 //  Debug streamer levels:
105 //  
106 //-------------------------------------------------------
107
108
109 /* $Id$ */
110
111 #include "Riostream.h"
112 #include <TClonesArray.h>
113 #include <TFile.h>
114 #include <TObjArray.h>
115 #include <TTree.h>
116 #include <TGraphErrors.h>
117 #include "AliLog.h"
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
122 #include "AliKink.h"
123 #include "AliV0.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
132
133 #include "AliTPCtrackerSector.h" 
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
139 #include "TRandom.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
144
145 //
146
147 ClassImp(AliTPCtrackerMI)
148
149
150
151 class AliTPCFastMath {
152 public:
153   AliTPCFastMath();  
154   static Double_t FastAsin(Double_t x);   
155  private: 
156   static Double_t fgFastAsin[20000];  //lookup table for fast asin computation
157 };
158
159 Double_t AliTPCFastMath::fgFastAsin[20000];
160 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
161
162 AliTPCFastMath::AliTPCFastMath(){
163   //
164   // initialized lookup table;
165   for (Int_t i=0;i<10000;i++){
166     fgFastAsin[2*i] = TMath::ASin(i/10000.);
167     fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
168   }
169 }
170
171 Double_t AliTPCFastMath::FastAsin(Double_t x){
172   //
173   // return asin using lookup table
174   if (x>0){
175     Int_t index = int(x*10000);
176     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177   }
178   x*=-1;
179   Int_t index = int(x*10000);
180   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
181 }
182 //__________________________________________________________________
183 AliTPCtrackerMI::AliTPCtrackerMI()
184                 :AliTracker(),
185                  fkNIS(0),
186                  fInnerSec(0),
187                  fkNOS(0),
188                  fOuterSec(0),
189                  fN(0),
190                  fSectors(0),
191                  fInput(0),
192                  fOutput(0),
193                  fSeedTree(0),
194                  fTreeDebug(0),
195                  fEvent(0),
196                  fDebug(0),
197                  fNewIO(kFALSE),
198                  fNtracks(0),
199                  fSeeds(0),
200                  fIteration(0),
201                  fkParam(0),
202                  fDebugStreamer(0)
203 {
204   //
205   // default constructor
206   //
207   for (Int_t irow=0; irow<200; irow++){
208     fXRow[irow]=0;
209     fYMax[irow]=0;
210     fPadLength[irow]=0;
211   }
212
213 }
214 //_____________________________________________________________________
215
216
217
218 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
219   //
220   //update track information using current cluster - track->fCurrentCluster
221
222
223   AliTPCclusterMI* c =track->GetCurrentCluster();
224   if (accept > 0) //sign not accepted clusters
225     track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);  
226   else // unsign accpeted clusters
227     track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);  
228   UInt_t i = track->GetCurrentClusterIndex1();
229
230   Int_t sec=(i&0xff000000)>>24; 
231   //Int_t row = (i&0x00ff0000)>>16; 
232   track->SetRow((i&0x00ff0000)>>16);
233   track->SetSector(sec);
234   //  Int_t index = i&0xFFFF;
235   if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow()); 
236   track->SetClusterIndex2(track->GetRow(), i);  
237   //track->fFirstPoint = row;
238   //if ( track->fLastPoint<row) track->fLastPoint =row;
239   //  if (track->fRow<0 || track->fRow>160) {
240   //  printf("problem\n");
241   //}
242   if (track->GetFirstPoint()>track->GetRow()) 
243     track->SetFirstPoint(track->GetRow());
244   if (track->GetLastPoint()<track->GetRow()) 
245     track->SetLastPoint(track->GetRow());
246   
247
248   track->SetClusterPointer(track->GetRow(),c);  
249   //
250
251   Double_t angle2 = track->GetSnp()*track->GetSnp();
252   //
253   //SET NEW Track Point
254   //
255   if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
256   {
257     angle2 = TMath::Sqrt(angle2/(1-angle2)); 
258     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->GetRow()));
259     //
260     point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
261     point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
262     point.SetErrY(sqrt(track->GetErrorY2()));
263     point.SetErrZ(sqrt(track->GetErrorZ2()));
264     //
265     point.SetX(track->GetX());
266     point.SetY(track->GetY());
267     point.SetZ(track->GetZ());
268     point.SetAngleY(angle2);
269     point.SetAngleZ(track->GetTgl());
270     if (point.IsShared()){
271       track->SetErrorY2(track->GetErrorY2()*4);
272       track->SetErrorZ2(track->GetErrorZ2()*4);
273     }
274   }  
275
276   Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
277   //
278 //   track->SetErrorY2(track->GetErrorY2()*1.3);
279 //   track->SetErrorY2(track->GetErrorY2()+0.01);    
280 //   track->SetErrorZ2(track->GetErrorZ2()*1.3);   
281 //   track->SetErrorZ2(track->GetErrorZ2()+0.005);      
282     //}
283   if (accept>0) return 0;
284   if (track->GetNumberOfClusters()%20==0){
285     //    if (track->fHelixIn){
286     //  TClonesArray & larr = *(track->fHelixIn);    
287     //  Int_t ihelix = larr.GetEntriesFast();
288     //  new(larr[ihelix]) AliHelix(*track) ;    
289     //}
290   }
291   track->SetNoCluster(0);
292   return track->Update(c,chi2,i);
293 }
294
295
296
297 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
298 {
299   //
300   // decide according desired precision to accept given 
301   // cluster for tracking
302   Double_t  yt=0,zt=0;
303   seed->GetProlongation(cluster->GetX(),yt,zt);
304   Double_t sy2=ErrY2(seed,cluster);
305   Double_t sz2=ErrZ2(seed,cluster);
306   
307   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
308   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
309   Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
310   Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
311   Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
312     (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
313   Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
314     (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
315   
316   Double_t rdistance2  = rdistancey2+rdistancez2;
317   //Int_t  accept =0;
318   
319   if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
320     Float_t rmsy2 = seed->GetCurrentSigmaY2();
321     Float_t rmsz2 = seed->GetCurrentSigmaZ2();
322     Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
323     Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
324     Float_t rmsy2p30R  = seed->GetCMeanSigmaY2p30R();
325     Float_t rmsz2p30R  = seed->GetCMeanSigmaZ2p30R();
326     AliExternalTrackParam param(*seed);
327     static TVectorD gcl(3),gtr(3);
328     Float_t gclf[3];
329     param.GetXYZ(gcl.GetMatrixArray());
330     cluster->GetGlobalXYZ(gclf);
331     gcl[0]=gclf[0];    gcl[1]=gclf[1];    gcl[2]=gclf[2];
332
333     
334     if (AliTPCReconstructor::StreamLevel()>2) {
335     (*fDebugStreamer)<<"ErrParam"<<
336       "Cl.="<<cluster<<
337       "T.="<<&param<<
338       "dy="<<dy<<
339       "dz="<<dz<<
340       "yt="<<yt<<
341       "zt="<<zt<<
342       "gcl.="<<&gcl<<
343       "gtr.="<<&gtr<<
344       "erry2="<<sy2<<
345       "errz2="<<sz2<<
346       "rmsy2="<<rmsy2<<
347       "rmsz2="<<rmsz2<< 
348       "rmsy2p30="<<rmsy2p30<<
349       "rmsz2p30="<<rmsz2p30<<   
350       "rmsy2p30R="<<rmsy2p30R<<
351       "rmsz2p30R="<<rmsz2p30R<< 
352       // normalize distance - 
353       "rdisty="<<rdistancey2<<
354       "rdistz="<<rdistancez2<<
355       "rdist="<<rdistance2<< //       
356       "\n";
357     }
358   }
359   //return 0;  // temporary
360   if (rdistance2>32) return 3;
361   
362   
363   if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)  
364     return 2;  //suspisiouce - will be changed
365   
366   if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)  
367     // strict cut on overlaped cluster
368     return  2;  //suspisiouce - will be changed
369   
370   if ( (rdistancey2>1. || rdistancez2>6.25 ) 
371        && cluster->GetType()<0){
372     seed->SetNFoundable(seed->GetNFoundable()-1);
373     return 2;    
374   }
375
376   if (AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 3 ||
377       AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 4 ) {
378     if(!AliTPCReconstructor::GetRecoParam()->GetUseOnePadCluster())
379       if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
380         return 2;
381   }
382   return 0;
383 }
384
385
386
387
388
389 //_____________________________________________________________________________
390 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par): 
391 AliTracker(), 
392                  fkNIS(par->GetNInnerSector()/2),
393                  fInnerSec(0),
394                  fkNOS(par->GetNOuterSector()/2),
395                  fOuterSec(0),
396                  fN(0),
397                  fSectors(0),
398                  fInput(0),
399                  fOutput(0),
400                  fSeedTree(0),
401                  fTreeDebug(0),
402                  fEvent(0),
403                  fDebug(0),
404                  fNewIO(0),
405                  fNtracks(0),
406                  fSeeds(0),
407                  fIteration(0),
408                  fkParam(0),
409                  fDebugStreamer(0)
410 {
411   //---------------------------------------------------------------------
412   // The main TPC tracker constructor
413   //---------------------------------------------------------------------
414   fInnerSec=new AliTPCtrackerSector[fkNIS];         
415   fOuterSec=new AliTPCtrackerSector[fkNOS];
416  
417   Int_t i;
418   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
419   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
420
421   fkParam = par;  
422   Int_t nrowlow = par->GetNRowLow();
423   Int_t nrowup = par->GetNRowUp();
424
425   
426   for (i=0;i<nrowlow;i++){
427     fXRow[i]     = par->GetPadRowRadiiLow(i);
428     fPadLength[i]= par->GetPadPitchLength(0,i);
429     fYMax[i]     = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
430   }
431
432   
433   for (i=0;i<nrowup;i++){
434     fXRow[i+nrowlow]      = par->GetPadRowRadiiUp(i);
435     fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
436     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
437   }
438
439   if (AliTPCReconstructor::StreamLevel()>0) {
440     fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
441   }
442 }
443 //________________________________________________________________________
444 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
445   AliTracker(t),
446                  fkNIS(t.fkNIS),
447                  fInnerSec(0),
448                  fkNOS(t.fkNOS),
449                  fOuterSec(0),
450                  fN(0),
451                  fSectors(0),
452                  fInput(0),
453                  fOutput(0),
454                  fSeedTree(0),
455                  fTreeDebug(0),
456                  fEvent(0),
457                  fDebug(0),
458                  fNewIO(kFALSE),
459                  fNtracks(0),
460                  fSeeds(0),
461                  fIteration(0),
462                  fkParam(0),
463                  fDebugStreamer(0)
464 {
465   //------------------------------------
466   // dummy copy constructor
467   //------------------------------------------------------------------
468   fOutput=t.fOutput;
469   for (Int_t irow=0; irow<200; irow++){
470     fXRow[irow]=0;
471     fYMax[irow]=0;
472     fPadLength[irow]=0;
473   }
474
475 }
476 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
477 {
478   //------------------------------
479   // dummy 
480   //--------------------------------------------------------------
481   return *this;
482 }
483 //_____________________________________________________________________________
484 AliTPCtrackerMI::~AliTPCtrackerMI() {
485   //------------------------------------------------------------------
486   // TPC tracker destructor
487   //------------------------------------------------------------------
488   delete[] fInnerSec;
489   delete[] fOuterSec;
490   if (fSeeds) {
491     fSeeds->Delete(); 
492     delete fSeeds;
493   }
494   if (fDebugStreamer) delete fDebugStreamer;
495 }
496
497
498 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
499 {
500   //
501   //
502   //fill esds using updated tracks
503   if (!fEvent) return;
504   
505     // write tracks to the event
506     // store index of the track
507     Int_t nseed=arr->GetEntriesFast();
508     //FindKinks(arr,fEvent);
509     for (Int_t i=0; i<nseed; i++) {
510       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
511       if (!pt) continue; 
512       pt->UpdatePoints();
513       AddCovariance(pt);
514       if (AliTPCReconstructor::StreamLevel()>1) {
515         (*fDebugStreamer)<<"Track0"<<
516           "Tr0.="<<pt<<
517           "\n";       
518       }
519       //      pt->PropagateTo(fkParam->GetInnerRadiusLow());
520       if (pt->GetKinkIndex(0)<=0){  //don't propagate daughter tracks 
521         pt->PropagateTo(fkParam->GetInnerRadiusLow());
522       }
523  
524       if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
525         AliESDtrack iotrack;
526         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
527         iotrack.SetTPCPoints(pt->GetPoints());
528         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
529         iotrack.SetV0Indexes(pt->GetV0Indexes());
530         //      iotrack.SetTPCpid(pt->fTPCr);
531         //iotrack.SetTPCindex(i);
532         fEvent->AddTrack(&iotrack);
533         continue;
534       }
535        
536       if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
537         AliESDtrack iotrack;
538         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
539         iotrack.SetTPCPoints(pt->GetPoints());
540         //iotrack.SetTPCindex(i);
541         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
542         iotrack.SetV0Indexes(pt->GetV0Indexes());
543         //      iotrack.SetTPCpid(pt->fTPCr);
544         fEvent->AddTrack(&iotrack);
545         continue;
546       } 
547       //
548       // short tracks  - maybe decays
549
550       if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
551         Int_t found,foundable,shared;
552         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
553         if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
554           AliESDtrack iotrack;
555           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);    
556           //iotrack.SetTPCindex(i);
557           iotrack.SetTPCPoints(pt->GetPoints());
558           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
559           iotrack.SetV0Indexes(pt->GetV0Indexes());
560           //iotrack.SetTPCpid(pt->fTPCr);
561           fEvent->AddTrack(&iotrack);
562           continue;
563         }
564       }       
565       
566       if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
567         Int_t found,foundable,shared;
568         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
569         if (found<20) continue;
570         if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
571         //
572         AliESDtrack iotrack;
573         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
574         iotrack.SetTPCPoints(pt->GetPoints());
575         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
576         iotrack.SetV0Indexes(pt->GetV0Indexes());
577         //iotrack.SetTPCpid(pt->fTPCr);
578         //iotrack.SetTPCindex(i);
579         fEvent->AddTrack(&iotrack);
580         continue;
581       }   
582       // short tracks  - secondaties
583       //
584       if ( (pt->GetNumberOfClusters()>30) ) {
585         Int_t found,foundable,shared;
586         pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
587         if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
588           AliESDtrack iotrack;
589           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);    
590           iotrack.SetTPCPoints(pt->GetPoints());
591           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
592           iotrack.SetV0Indexes(pt->GetV0Indexes());
593           //iotrack.SetTPCpid(pt->fTPCr);       
594           //iotrack.SetTPCindex(i);
595           fEvent->AddTrack(&iotrack);
596           continue;
597         }
598       }       
599       
600       if ( (pt->GetNumberOfClusters()>15)) {
601         Int_t found,foundable,shared;
602         pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
603         if (found<15) continue;
604         if (foundable<=0) continue;
605         if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
606         if (float(found)/float(foundable)<0.8) continue;
607         //
608         AliESDtrack iotrack;
609         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
610         iotrack.SetTPCPoints(pt->GetPoints());
611         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
612         iotrack.SetV0Indexes(pt->GetV0Indexes());
613         //      iotrack.SetTPCpid(pt->fTPCr);
614         //iotrack.SetTPCindex(i);
615         fEvent->AddTrack(&iotrack);
616         continue;
617       }   
618     }
619     // >> account for suppressed tracks in the kink indices (RS)
620     int nESDtracks = fEvent->GetNumberOfTracks();
621     for (int it=nESDtracks;it--;) {
622       AliESDtrack* esdTr = fEvent->GetTrack(it);
623       if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
624       for (int ik=0;ik<3;ik++) {
625         int knkId=0;
626         if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
627         AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
628         if (!kink) {
629           AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
630           continue;
631         }
632         kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
633       }
634     }
635     // << account for suppressed tracks in the kink indices (RS)  
636     AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
637   
638 }
639
640
641
642
643
644 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
645   //
646   //
647   // Use calibrated cluster error from OCDB
648   //
649   AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
650   //
651   Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
652   Int_t ctype = cl->GetType();  
653   Int_t    type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
654   Double_t angle = seed->GetSnp()*seed->GetSnp();
655   angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
656   Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
657   if (ctype<0) {
658     erry2+=0.5;  // edge cluster
659   }
660   erry2*=erry2;
661   seed->SetErrorY2(erry2);
662   //
663   return erry2;
664
665 //calculate look-up table at the beginning
666 //   static Bool_t  ginit = kFALSE;
667 //   static Float_t gnoise1,gnoise2,gnoise3;
668 //   static Float_t ggg1[10000];
669 //   static Float_t ggg2[10000];
670 //   static Float_t ggg3[10000];
671 //   static Float_t glandau1[10000];
672 //   static Float_t glandau2[10000];
673 //   static Float_t glandau3[10000];
674 //   //
675 //   static Float_t gcor01[500];
676 //   static Float_t gcor02[500];
677 //   static Float_t gcorp[500];
678 //   //
679
680 //   //
681 //   if (ginit==kFALSE){
682 //     for (Int_t i=1;i<500;i++){
683 //       Float_t rsigma = float(i)/100.;
684 //       gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
685 //       gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
686 //       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
687 //     }
688
689 //     //
690 //     for (Int_t i=3;i<10000;i++){
691 //       //
692 //       //
693 //       // inner sector
694 //       Float_t amp = float(i);
695 //       Float_t padlength =0.75;
696 //       gnoise1 = 0.0004/padlength;
697 //       Float_t nel     = 0.268*amp;
698 //       Float_t nprim   = 0.155*amp;
699 //       ggg1[i]          = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
700 //       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
701 //       if (glandau1[i]>1) glandau1[i]=1;
702 //       glandau1[i]*=padlength*padlength/12.;      
703 //       //
704 //       // outer short
705 //       padlength =1.;
706 //       gnoise2   = 0.0004/padlength;
707 //       nel       = 0.3*amp;
708 //       nprim     = 0.133*amp;
709 //       ggg2[i]      = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
710 //       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
711 //       if (glandau2[i]>1) glandau2[i]=1;
712 //       glandau2[i]*=padlength*padlength/12.;
713 //       //
714 //       //
715 //       // outer long
716 //       padlength =1.5;
717 //       gnoise3   = 0.0004/padlength;
718 //       nel       = 0.3*amp;
719 //       nprim     = 0.133*amp;
720 //       ggg3[i]      = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
721 //       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
722 //       if (glandau3[i]>1) glandau3[i]=1;
723 //       glandau3[i]*=padlength*padlength/12.;
724 //       //
725 //     }
726 //     ginit = kTRUE;
727 //   }
728 //   //
729 //   //
730 //   //
731 //   Int_t amp = int(TMath::Abs(cl->GetQ()));  
732 //   if (amp>9999) {
733 //     seed->SetErrorY2(1.);
734 //     return 1.;
735 //   }
736 //   Float_t snoise2;
737 //   Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
738 //   Int_t ctype = cl->GetType();  
739 //   Float_t padlength= GetPadPitchLength(seed->GetRow());
740 //   Double_t angle2 = seed->GetSnp()*seed->GetSnp();
741 //   angle2 = angle2/(1-angle2); 
742 //   //
743 //   //cluster "quality"
744 //   Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
745 //   Float_t res;
746 //   //
747 //   if (fSectors==fInnerSec){
748 //     snoise2 = gnoise1;
749 //     res     = ggg1[amp]*z+glandau1[amp]*angle2;     
750 //     if (ctype==0) res *= gcor01[rsigmay];
751 //     if ((ctype>0)){
752 //       res+=0.002;
753 //       res*= gcorp[rsigmay];
754 //     }
755 //   }
756 //   else {
757 //     if (padlength<1.1){
758 //       snoise2 = gnoise2;
759 //       res     = ggg2[amp]*z+glandau2[amp]*angle2; 
760 //       if (ctype==0) res *= gcor02[rsigmay];      
761 //       if ((ctype>0)){
762 //      res+=0.002;
763 //      res*= gcorp[rsigmay];
764 //       }
765 //     }
766 //     else{
767 //       snoise2 = gnoise3;      
768 //       res     = ggg3[amp]*z+glandau3[amp]*angle2; 
769 //       if (ctype==0) res *= gcor02[rsigmay];
770 //       if ((ctype>0)){
771 //      res+=0.002;
772 //      res*= gcorp[rsigmay];
773 //       }
774 //     }
775 //   }  
776
777 //   if (ctype<0){
778 //     res+=0.005;
779 //     res*=2.4;  // overestimate error 2 times
780 //   }
781 //   res+= snoise2;
782  
783 //   if (res<2*snoise2)
784 //     res = 2*snoise2;
785   
786 //   seed->SetErrorY2(res);
787 //   return res;
788
789
790 }
791
792
793
794 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
795   //
796   //
797   // Use calibrated cluster error from OCDB
798   //
799   AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
800   //
801   Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
802   Int_t ctype = cl->GetType();  
803   Int_t    type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
804   //
805   Double_t angle2 = seed->GetSnp()*seed->GetSnp();
806   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2)); 
807   Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
808   Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
809   if (ctype<0) {
810     errz2+=0.5;  // edge cluster
811   }
812   errz2*=errz2;
813   seed->SetErrorZ2(errz2);
814   //
815   return errz2;
816
817
818
819 //   //seed->SetErrorY2(0.1);
820 //   //return 0.1;
821 //   //calculate look-up table at the beginning
822 //   static Bool_t  ginit = kFALSE;
823 //   static Float_t gnoise1,gnoise2,gnoise3;
824 //   static Float_t ggg1[10000];
825 //   static Float_t ggg2[10000];
826 //   static Float_t ggg3[10000];
827 //   static Float_t glandau1[10000];
828 //   static Float_t glandau2[10000];
829 //   static Float_t glandau3[10000];
830 //   //
831 //   static Float_t gcor01[1000];
832 //   static Float_t gcor02[1000];
833 //   static Float_t gcorp[1000];
834 //   //
835
836 //   //
837 //   if (ginit==kFALSE){
838 //     for (Int_t i=1;i<1000;i++){
839 //       Float_t rsigma = float(i)/100.;
840 //       gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
841 //       gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
842 //       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
843 //     }
844
845 //     //
846 //     for (Int_t i=3;i<10000;i++){
847 //       //
848 //       //
849 //       // inner sector
850 //       Float_t amp = float(i);
851 //       Float_t padlength =0.75;
852 //       gnoise1 = 0.0004/padlength;
853 //       Float_t nel     = 0.268*amp;
854 //       Float_t nprim   = 0.155*amp;
855 //       ggg1[i]          = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
856 //       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
857 //       if (glandau1[i]>1) glandau1[i]=1;
858 //       glandau1[i]*=padlength*padlength/12.;      
859 //       //
860 //       // outer short
861 //       padlength =1.;
862 //       gnoise2   = 0.0004/padlength;
863 //       nel       = 0.3*amp;
864 //       nprim     = 0.133*amp;
865 //       ggg2[i]      = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
866 //       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
867 //       if (glandau2[i]>1) glandau2[i]=1;
868 //       glandau2[i]*=padlength*padlength/12.;
869 //       //
870 //       //
871 //       // outer long
872 //       padlength =1.5;
873 //       gnoise3   = 0.0004/padlength;
874 //       nel       = 0.3*amp;
875 //       nprim     = 0.133*amp;
876 //       ggg3[i]      = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
877 //       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
878 //       if (glandau3[i]>1) glandau3[i]=1;
879 //       glandau3[i]*=padlength*padlength/12.;
880 //       //
881 //     }
882 //     ginit = kTRUE;
883 //   }
884 //   //
885 //   //
886 //   //
887 //   Int_t amp = int(TMath::Abs(cl->GetQ()));  
888 //   if (amp>9999) {
889 //     seed->SetErrorY2(1.);
890 //     return 1.;
891 //   }
892 //   Float_t snoise2;
893 //   Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
894 //   Int_t ctype = cl->GetType();  
895 //   Float_t padlength= GetPadPitchLength(seed->GetRow());
896 //   //
897 //   Double_t angle2 = seed->GetSnp()*seed->GetSnp();
898 //   //  if (angle2<0.6) angle2 = 0.6;
899 //   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2)); 
900 //   //
901 //   //cluster "quality"
902 //   Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
903 //   Float_t res;
904 //   //
905 //   if (fSectors==fInnerSec){
906 //     snoise2 = gnoise1;
907 //     res     = ggg1[amp]*z+glandau1[amp]*angle2;     
908 //     if (ctype==0) res *= gcor01[rsigmaz];
909 //     if ((ctype>0)){
910 //       res+=0.002;
911 //       res*= gcorp[rsigmaz];
912 //     }
913 //   }
914 //   else {
915 //     if (padlength<1.1){
916 //       snoise2 = gnoise2;
917 //       res     = ggg2[amp]*z+glandau2[amp]*angle2; 
918 //       if (ctype==0) res *= gcor02[rsigmaz];      
919 //       if ((ctype>0)){
920 //      res+=0.002;
921 //      res*= gcorp[rsigmaz];
922 //       }
923 //     }
924 //     else{
925 //       snoise2 = gnoise3;      
926 //       res     = ggg3[amp]*z+glandau3[amp]*angle2; 
927 //       if (ctype==0) res *= gcor02[rsigmaz];
928 //       if ((ctype>0)){
929 //      res+=0.002;
930 //      res*= gcorp[rsigmaz];
931 //       }
932 //     }
933 //   }  
934
935 //   if (ctype<0){
936 //     res+=0.002;
937 //     res*=1.3;
938 //   }
939 //   if ((ctype<0) &&amp<70){
940 //     res+=0.002;
941 //     res*=1.3;  
942 //   }
943 //   res += snoise2;
944 //   if (res<2*snoise2)
945 //      res = 2*snoise2;
946 //   if (res>3) res =3;
947 //   seed->SetErrorZ2(res);
948 //   return res;
949 }
950
951
952
953
954
955 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
956 {
957   //rotate to track "local coordinata
958   Float_t x = seed->GetX();
959   Float_t y = seed->GetY();
960   Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
961   
962   if (y > ymax) {
963     seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
964     if (!seed->Rotate(fSectors->GetAlpha())) 
965       return;
966   } else if (y <-ymax) {
967     seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
968     if (!seed->Rotate(-fSectors->GetAlpha())) 
969       return;
970   }   
971
972 }
973
974
975
976 //_____________________________________________________________________________
977 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
978                    Double_t x2,Double_t y2,
979                    Double_t x3,Double_t y3) const
980 {
981   //-----------------------------------------------------------------
982   // Initial approximation of the track curvature
983   //-----------------------------------------------------------------
984   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
985   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
986                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
987   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
988                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
989
990   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
991   if ( xr*xr+yr*yr<=0.00000000000001) return 100;
992   return -xr*yr/sqrt(xr*xr+yr*yr); 
993 }
994
995
996
997 //_____________________________________________________________________________
998 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
999                    Double_t x2,Double_t y2,
1000                    Double_t x3,Double_t y3) const
1001 {
1002   //-----------------------------------------------------------------
1003   // Initial approximation of the track curvature
1004   //-----------------------------------------------------------------
1005   x3 -=x1;
1006   x2 -=x1;
1007   y3 -=y1;
1008   y2 -=y1;
1009   //  
1010   Double_t det = x3*y2-x2*y3;
1011   if (TMath::Abs(det)<1e-10){
1012     return 100;
1013   }
1014   //
1015   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1016   Double_t x0 = x3*0.5-y3*u;
1017   Double_t y0 = y3*0.5+x3*u;
1018   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1019   if (det<0) c2*=-1;
1020   return c2;
1021 }
1022
1023
1024 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1025                    Double_t x2,Double_t y2,
1026                    Double_t x3,Double_t y3) const 
1027 {
1028   //-----------------------------------------------------------------
1029   // Initial approximation of the track curvature
1030   //-----------------------------------------------------------------
1031   x3 -=x1;
1032   x2 -=x1;
1033   y3 -=y1;
1034   y2 -=y1;
1035   //  
1036   Double_t det = x3*y2-x2*y3;
1037   if (TMath::Abs(det)<1e-10) {
1038     return 100;
1039   }
1040   //
1041   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1042   Double_t x0 = x3*0.5-y3*u; 
1043   Double_t y0 = y3*0.5+x3*u;
1044   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1045   if (det<0) c2*=-1;
1046   x0+=x1;
1047   x0*=c2;  
1048   return x0;
1049 }
1050
1051
1052
1053 //_____________________________________________________________________________
1054 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1055                    Double_t x2,Double_t y2,
1056                    Double_t x3,Double_t y3) const
1057 {
1058   //-----------------------------------------------------------------
1059   // Initial approximation of the track curvature times center of curvature
1060   //-----------------------------------------------------------------
1061   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1062   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1063                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1064   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1065                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1066
1067   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1068   
1069   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1070 }
1071
1072 //_____________________________________________________________________________
1073 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1, 
1074                    Double_t x2,Double_t y2,
1075                    Double_t z1,Double_t z2) const
1076 {
1077   //-----------------------------------------------------------------
1078   // Initial approximation of the tangent of the track dip angle
1079   //-----------------------------------------------------------------
1080   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1081 }
1082
1083
1084 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1, 
1085                    Double_t x2,Double_t y2,
1086                    Double_t z1,Double_t z2, Double_t c) const
1087 {
1088   //-----------------------------------------------------------------
1089   // Initial approximation of the tangent of the track dip angle
1090   //-----------------------------------------------------------------
1091
1092   //  Double_t angle1;
1093   
1094   //angle1    =  (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1095   //
1096   Double_t d  =  TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1097   if (TMath::Abs(d*c*0.5)>1) return 0;
1098   //  Double_t   angle2    =  TMath::ASin(d*c*0.5);
1099   //  Double_t   angle2    =  AliTPCFastMath::FastAsin(d*c*0.5);
1100   Double_t   angle2    = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1101
1102   angle2  = (z1-z2)*c/(angle2*2.);
1103   return angle2;
1104 }
1105
1106 Bool_t   AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1107 {//-----------------------------------------------------------------
1108   // This function find proloncation of a track to a reference plane x=x2.
1109   //-----------------------------------------------------------------
1110   
1111   Double_t dx=x2-x1;
1112
1113   if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {   
1114     return kFALSE;
1115   }
1116
1117   Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1118   Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));  
1119   y = x[0];
1120   z = x[1];
1121   
1122   Double_t dy = dx*(c1+c2)/(r1+r2);
1123   Double_t dz = 0;
1124   //
1125   Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1126   
1127   if (TMath::Abs(delta)>0.01){
1128     dz = x[3]*TMath::ASin(delta)/x[4];
1129   }else{
1130     dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1131   }
1132   
1133   //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1134
1135   y+=dy;
1136   z+=dz;
1137   
1138   return kTRUE;  
1139 }
1140
1141 Int_t  AliTPCtrackerMI::LoadClusters (TTree *const tree)
1142 {
1143   // load clusters
1144   //
1145   fInput = tree;
1146   return LoadClusters();
1147 }
1148
1149
1150 Int_t  AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1151 {
1152   //
1153   // load clusters to the memory
1154   AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1155   Int_t lower   = arr->LowerBound();
1156   Int_t entries = arr->GetEntriesFast();
1157   
1158   for (Int_t i=lower; i<entries; i++) {
1159     clrow = (AliTPCClustersRow*) arr->At(i);
1160     if(!clrow) continue;
1161     if(!clrow->GetArray()) continue;
1162     
1163     //  
1164     Int_t sec,row;
1165     fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1166     
1167     for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1168       Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1169     }
1170     //
1171     if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1172     AliTPCtrackerRow * tpcrow=0;
1173     Int_t left=0;
1174     if (sec<fkNIS*2){
1175       tpcrow = &(fInnerSec[sec%fkNIS][row]);    
1176       left = sec/fkNIS;
1177     }
1178     else{
1179       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1180       left = (sec-fkNIS*2)/fkNOS;
1181     }
1182     if (left ==0){
1183       tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1184       tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1185       for (Int_t j=0;j<tpcrow->GetN1();++j) 
1186         tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1187     }
1188     if (left ==1){
1189       tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1190       tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1191       for (Int_t j=0;j<tpcrow->GetN2();++j) 
1192         tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1193     }
1194     clrow->GetArray()->Clear("C");
1195   }
1196   //
1197   delete clrow;
1198   LoadOuterSectors();
1199   LoadInnerSectors();
1200   return 0;
1201 }
1202
1203 Int_t  AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1204 {
1205   //
1206   // load clusters to the memory from one 
1207   // TClonesArray
1208   //
1209   AliTPCclusterMI *clust=0;
1210   Int_t count[72][96] = { {0} , {0} }; 
1211
1212   // loop over clusters
1213   for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1214     clust = (AliTPCclusterMI*)arr->At(icl);
1215     if(!clust) continue;
1216     //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1217
1218     // transform clusters
1219     Transform(clust);
1220
1221     // count clusters per pad row
1222     count[clust->GetDetector()][clust->GetRow()]++;
1223   }
1224
1225   // insert clusters to sectors
1226   for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1227     clust = (AliTPCclusterMI*)arr->At(icl);
1228     if(!clust) continue;
1229
1230     Int_t sec = clust->GetDetector();
1231     Int_t row = clust->GetRow();
1232
1233     // filter overlapping pad rows needed by HLT
1234     if(sec<fkNIS*2) { //IROCs
1235      if(row == 30) continue;
1236     }
1237     else { // OROCs
1238       if(row == 27 || row == 76) continue;
1239     }
1240
1241     Int_t left=0;
1242     if (sec<fkNIS*2){
1243       left = sec/fkNIS;
1244       fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);    
1245     }
1246     else{
1247       left = (sec-fkNIS*2)/fkNOS;
1248       fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1249     }
1250   }
1251
1252   // Load functions must be called behind LoadCluster(TClonesArray*)
1253   // needed by HLT
1254   //LoadOuterSectors();
1255   //LoadInnerSectors();
1256
1257   return 0;
1258 }
1259
1260
1261 Int_t  AliTPCtrackerMI::LoadClusters()
1262 {
1263   //
1264   // load clusters to the memory
1265   AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1266   //
1267   //  TTree * tree = fClustersArray.GetTree();
1268
1269   TTree * tree = fInput;
1270   TBranch * br = tree->GetBranch("Segment");
1271   br->SetAddress(&clrow);
1272   //
1273   Int_t j=Int_t(tree->GetEntries());
1274   for (Int_t i=0; i<j; i++) {
1275     br->GetEntry(i);
1276     //  
1277     Int_t sec,row;
1278     fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1279     for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1280       Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1281     }
1282     //
1283     AliTPCtrackerRow * tpcrow=0;
1284     Int_t left=0;
1285     if (sec<fkNIS*2){
1286       tpcrow = &(fInnerSec[sec%fkNIS][row]);    
1287       left = sec/fkNIS;
1288     }
1289     else{
1290       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1291       left = (sec-fkNIS*2)/fkNOS;
1292     }
1293     if (left ==0){
1294       tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1295       tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1296       for (Int_t k=0;k<tpcrow->GetN1();++k) 
1297         tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1298     }
1299     if (left ==1){
1300       tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1301       tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1302       for (Int_t k=0;k<tpcrow->GetN2();++k) 
1303         tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1304     }
1305   }
1306   //
1307   delete clrow;
1308   LoadOuterSectors();
1309   LoadInnerSectors();
1310   return 0;
1311 }
1312
1313
1314 void AliTPCtrackerMI::UnloadClusters()
1315 {
1316   //
1317   // unload clusters from the memory
1318   //
1319   Int_t nrows = fOuterSec->GetNRows();
1320   for (Int_t sec = 0;sec<fkNOS;sec++)
1321     for (Int_t row = 0;row<nrows;row++){
1322       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1323       //      if (tpcrow){
1324       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1325       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1326       //}
1327       tpcrow->ResetClusters();
1328     }
1329   //
1330   nrows = fInnerSec->GetNRows();
1331   for (Int_t sec = 0;sec<fkNIS;sec++)
1332     for (Int_t row = 0;row<nrows;row++){
1333       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1334       //if (tpcrow){
1335       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1336       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1337       //}
1338       tpcrow->ResetClusters();
1339     }
1340
1341   return ;
1342 }
1343
1344 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1345   //
1346   // Filling cluster to the array - For visualization purposes
1347   //
1348   Int_t nrows=0;
1349   nrows = fOuterSec->GetNRows();
1350   for (Int_t sec = 0;sec<fkNOS;sec++)
1351     for (Int_t row = 0;row<nrows;row++){
1352       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1353       if (!tpcrow) continue;
1354       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1355         array->AddLast((TObject*)((*tpcrow)[icl]));
1356       }
1357     } 
1358   nrows = fInnerSec->GetNRows();
1359   for (Int_t sec = 0;sec<fkNIS;sec++)
1360     for (Int_t row = 0;row<nrows;row++){
1361       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1362       if (!tpcrow) continue;
1363       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1364         array->AddLast((TObject*)(*tpcrow)[icl]);
1365       }
1366     }
1367 }
1368
1369
1370 void   AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1371   //
1372   // transformation
1373   //
1374   AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1375   AliTPCTransform *transform = calibDB->GetTransform() ;
1376   if (!transform) {
1377     AliFatal("Tranformations not in calibDB");
1378     return;
1379   }
1380   transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1381   Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1382   Int_t i[1]={cluster->GetDetector()};
1383   transform->Transform(x,i,0,1);  
1384   //  if (cluster->GetDetector()%36>17){
1385   //  x[1]*=-1;
1386   //}
1387
1388   //
1389   // in debug mode  check the transformation
1390   //
1391   if (AliTPCReconstructor::StreamLevel()>2) {
1392     Float_t gx[3];
1393     cluster->GetGlobalXYZ(gx);
1394     Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1395     TTreeSRedirector &cstream = *fDebugStreamer;
1396     cstream<<"Transform"<<
1397       "event="<<event<<
1398       "x0="<<x[0]<<
1399       "x1="<<x[1]<<
1400       "x2="<<x[2]<<
1401       "gx0="<<gx[0]<<
1402       "gx1="<<gx[1]<<
1403       "gx2="<<gx[2]<<
1404       "Cl.="<<cluster<<
1405       "\n"; 
1406   }
1407   cluster->SetX(x[0]);
1408   cluster->SetY(x[1]);
1409   cluster->SetZ(x[2]);
1410   // The old stuff:
1411   //
1412   // 
1413   //
1414   //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1415   if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1416     TGeoHMatrix  *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1417     //TGeoHMatrix  mat;
1418     Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1419     Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1420     if (mat) mat->LocalToMaster(pos,posC);
1421     else{
1422       // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1423     }
1424     cluster->SetX(posC[0]);
1425     cluster->SetY(posC[1]);
1426     cluster->SetZ(posC[2]);
1427   }
1428 }
1429
1430 //_____________________________________________________________________________
1431 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1432   //-----------------------------------------------------------------
1433   // This function fills outer TPC sectors with clusters.
1434   //-----------------------------------------------------------------
1435   Int_t nrows = fOuterSec->GetNRows();
1436   UInt_t index=0;
1437   for (Int_t sec = 0;sec<fkNOS;sec++)
1438     for (Int_t row = 0;row<nrows;row++){
1439       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);  
1440       Int_t sec2 = sec+2*fkNIS;
1441       //left
1442       Int_t ncl = tpcrow->GetN1();
1443       while (ncl--) {
1444         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1445         index=(((sec2<<8)+row)<<16)+ncl;
1446         tpcrow->InsertCluster(c,index);
1447       }
1448       //right
1449       ncl = tpcrow->GetN2();
1450       while (ncl--) {
1451         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1452         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1453         tpcrow->InsertCluster(c,index);
1454       }
1455       //
1456       // write indexes for fast acces
1457       //
1458       for (Int_t i=0;i<510;i++)
1459         tpcrow->SetFastCluster(i,-1);
1460       for (Int_t i=0;i<tpcrow->GetN();i++){
1461         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1462         tpcrow->SetFastCluster(zi,i);  // write index
1463       }
1464       Int_t last = 0;
1465       for (Int_t i=0;i<510;i++){
1466         if (tpcrow->GetFastCluster(i)<0)
1467           tpcrow->SetFastCluster(i,last);
1468         else
1469           last = tpcrow->GetFastCluster(i);
1470       }
1471     }  
1472   fN=fkNOS;
1473   fSectors=fOuterSec;
1474   return 0;
1475 }
1476
1477
1478 //_____________________________________________________________________________
1479 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1480   //-----------------------------------------------------------------
1481   // This function fills inner TPC sectors with clusters.
1482   //-----------------------------------------------------------------
1483   Int_t nrows = fInnerSec->GetNRows();
1484   UInt_t index=0;
1485   for (Int_t sec = 0;sec<fkNIS;sec++)
1486     for (Int_t row = 0;row<nrows;row++){
1487       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1488       //
1489       //left
1490       Int_t ncl = tpcrow->GetN1();
1491       while (ncl--) {
1492         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1493         index=(((sec<<8)+row)<<16)+ncl;
1494         tpcrow->InsertCluster(c,index);
1495       }
1496       //right
1497       ncl = tpcrow->GetN2();
1498       while (ncl--) {
1499         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1500         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1501         tpcrow->InsertCluster(c,index);
1502       }
1503       //
1504       // write indexes for fast acces
1505       //
1506       for (Int_t i=0;i<510;i++)
1507         tpcrow->SetFastCluster(i,-1);
1508       for (Int_t i=0;i<tpcrow->GetN();i++){
1509         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1510         tpcrow->SetFastCluster(zi,i);  // write index
1511       }
1512       Int_t last = 0;
1513       for (Int_t i=0;i<510;i++){
1514         if (tpcrow->GetFastCluster(i)<0)
1515           tpcrow->SetFastCluster(i,last);
1516         else
1517           last = tpcrow->GetFastCluster(i);
1518       }
1519
1520     }  
1521    
1522   fN=fkNIS;
1523   fSectors=fInnerSec;
1524   return 0;
1525 }
1526
1527
1528
1529 //_________________________________________________________________________
1530 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1531   //--------------------------------------------------------------------
1532   //       Return pointer to a given cluster
1533   //--------------------------------------------------------------------
1534   if (index<0) return 0; // no cluster
1535   Int_t sec=(index&0xff000000)>>24; 
1536   Int_t row=(index&0x00ff0000)>>16; 
1537   Int_t ncl=(index&0x00007fff)>>00;
1538
1539   const AliTPCtrackerRow * tpcrow=0;
1540   AliTPCclusterMI * clrow =0;
1541
1542   if (sec<0 || sec>=fkNIS*4) {
1543     AliWarning(Form("Wrong sector %d",sec));
1544     return 0x0;
1545   }
1546
1547   if (sec<fkNIS*2){
1548     AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1549     if (tracksec.GetNRows()<=row) return 0;
1550     tpcrow = &(tracksec[row]);
1551     if (tpcrow==0) return 0;
1552
1553     if (sec<fkNIS) {
1554       if (tpcrow->GetN1()<=ncl) return 0;
1555       clrow = tpcrow->GetClusters1();
1556     }
1557     else {
1558       if (tpcrow->GetN2()<=ncl) return 0;
1559       clrow = tpcrow->GetClusters2();
1560     }
1561   }
1562   else {
1563     AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1564     if (tracksec.GetNRows()<=row) return 0;
1565     tpcrow = &(tracksec[row]);
1566     if (tpcrow==0) return 0;
1567
1568     if (sec-2*fkNIS<fkNOS) {
1569       if (tpcrow->GetN1()<=ncl) return 0;
1570       clrow = tpcrow->GetClusters1();
1571     }
1572     else {
1573       if (tpcrow->GetN2()<=ncl) return 0;
1574       clrow = tpcrow->GetClusters2();
1575     }
1576   }
1577
1578   return &(clrow[ncl]);      
1579   
1580 }
1581
1582
1583
1584 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1585   //-----------------------------------------------------------------
1586   // This function tries to find a track prolongation to next pad row
1587   //-----------------------------------------------------------------
1588   //
1589   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1590   //
1591   //
1592   AliTPCclusterMI *cl=0;
1593   Int_t tpcindex= t.GetClusterIndex2(nr);
1594   //
1595   // update current shape info every 5 pad-row
1596   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1597     GetShape(&t,nr);    
1598     //}
1599   //  
1600   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters 
1601     //        
1602     if (tpcindex==-1) return 0; //track in dead zone
1603     if (tpcindex >= 0){     //
1604       cl = t.GetClusterPointer(nr);
1605       if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1606       t.SetCurrentClusterIndex1(tpcindex); 
1607     }
1608     if (cl){      
1609       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
1610       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1611       //
1612       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1613       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1614       
1615       if (TMath::Abs(angle-t.GetAlpha())>0.001){
1616         Double_t rotation = angle-t.GetAlpha();
1617         t.SetRelativeSector(relativesector);
1618         if (!t.Rotate(rotation)) return 0;      
1619       }
1620       if (!t.PropagateTo(x)) return 0;
1621       //
1622       t.SetCurrentCluster(cl); 
1623       t.SetRow(nr);
1624       Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1625       if ((tpcindex&0x8000)==0) accept =0;
1626       if (accept<3) { 
1627         //if founded cluster is acceptible
1628         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
1629           t.SetErrorY2(t.GetErrorY2()+0.03);
1630           t.SetErrorZ2(t.GetErrorZ2()+0.03); 
1631           t.SetErrorY2(t.GetErrorY2()*3);
1632           t.SetErrorZ2(t.GetErrorZ2()*3); 
1633         }
1634         t.SetNFoundable(t.GetNFoundable()+1);
1635         UpdateTrack(&t,accept);
1636         return 1;
1637       }
1638       else { // Remove old cluster from track
1639         t.SetClusterIndex(nr, -3);
1640         t.SetClusterPointer(nr, 0);
1641       }
1642     }
1643   }
1644   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;  // cut on angle
1645   if (fIteration>1 && IsFindable(t)){
1646     // not look for new cluster during refitting    
1647     t.SetNFoundable(t.GetNFoundable()+1);
1648     return 0;
1649   }
1650   //
1651   UInt_t index=0;
1652   //  if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1653   if (!t.PropagateTo(x)) {
1654     if (fIteration==0) t.SetRemoval(10);
1655     return 0;
1656   }
1657   Double_t y = t.GetY(); 
1658   if (TMath::Abs(y)>ymax){
1659     if (y > ymax) {
1660       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1661       if (!t.Rotate(fSectors->GetAlpha())) 
1662         return 0;
1663     } else if (y <-ymax) {
1664       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1665       if (!t.Rotate(-fSectors->GetAlpha())) 
1666         return 0;
1667     }
1668     if (!t.PropagateTo(x)) {
1669       if (fIteration==0) t.SetRemoval(10);
1670       return 0;
1671     }
1672     y = t.GetY();
1673   }
1674   //
1675   Double_t z=t.GetZ();
1676   //
1677
1678   if (!IsActive(t.GetRelativeSector(),nr)) {
1679     t.SetInDead(kTRUE);
1680     t.SetClusterIndex2(nr,-1); 
1681     return 0;
1682   }
1683   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1684   Bool_t isActive  = IsActive(t.GetRelativeSector(),nr);
1685   Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1686   
1687   if (!isActive || !isActive2) return 0;
1688
1689   const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1690   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1691   Double_t  roady  =1.;
1692   Double_t  roadz = 1.;
1693   //
1694   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1695     t.SetInDead(kTRUE);
1696     t.SetClusterIndex2(nr,-1); 
1697     return 0;
1698   } 
1699   else
1700     {
1701       if (IsFindable(t))
1702           //      if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) 
1703         t.SetNFoundable(t.GetNFoundable()+1);
1704       else
1705         return 0;
1706     }   
1707   //calculate 
1708   if (krow) {
1709     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);    
1710     cl = krow.FindNearest2(y,z,roady,roadz,index);    
1711     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));       
1712   }  
1713   if (cl) {
1714     t.SetCurrentCluster(cl); 
1715     t.SetRow(nr);
1716     if (fIteration==2&&cl->IsUsed(10)) return 0; 
1717     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1718     if (fIteration==2&&cl->IsUsed(11)) {
1719       t.SetErrorY2(t.GetErrorY2()+0.03);
1720       t.SetErrorZ2(t.GetErrorZ2()+0.03); 
1721       t.SetErrorY2(t.GetErrorY2()*3);
1722       t.SetErrorZ2(t.GetErrorZ2()*3); 
1723     }
1724     /*    
1725     if (t.fCurrentCluster->IsUsed(10)){
1726       //
1727       //     
1728
1729       t.fNShared++;
1730       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1731         t.fRemoval =10;
1732         return 0;
1733       }
1734     }
1735     */
1736     if (accept<3) UpdateTrack(&t,accept);
1737
1738   } else {  
1739     if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1740     
1741   }
1742   return 1;
1743 }
1744
1745
1746
1747 //_________________________________________________________________________
1748 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1749 {
1750   // Get track space point by index
1751   // return false in case the cluster doesn't exist
1752   AliTPCclusterMI *cl = GetClusterMI(index);
1753   if (!cl) return kFALSE;
1754   Int_t sector = (index&0xff000000)>>24;
1755   //  Int_t row = (index&0x00ff0000)>>16;
1756   Float_t xyz[3];
1757   //  xyz[0] = fkParam->GetPadRowRadii(sector,row);
1758   xyz[0] = cl->GetX();
1759   xyz[1] = cl->GetY();
1760   xyz[2] = cl->GetZ();
1761   Float_t sin,cos;
1762   fkParam->AdjustCosSin(sector,cos,sin);
1763   Float_t x = cos*xyz[0]-sin*xyz[1];
1764   Float_t y = cos*xyz[1]+sin*xyz[0];
1765   Float_t cov[6];
1766   Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1767   if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1768   Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1769   if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1770   cov[0] = sin*sin*sigmaY2;
1771   cov[1] = -sin*cos*sigmaY2;
1772   cov[2] = 0.;
1773   cov[3] = cos*cos*sigmaY2;
1774   cov[4] = 0.;
1775   cov[5] = sigmaZ2;
1776   p.SetXYZ(x,y,xyz[2],cov);
1777   AliGeomManager::ELayerID iLayer;
1778   Int_t idet;
1779   if (sector < fkParam->GetNInnerSector()) {
1780     iLayer = AliGeomManager::kTPC1;
1781     idet = sector;
1782   }
1783   else {
1784     iLayer = AliGeomManager::kTPC2;
1785     idet = sector - fkParam->GetNInnerSector();
1786   }
1787   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1788   p.SetVolumeID(volid);
1789   return kTRUE;
1790 }
1791
1792
1793
1794 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
1795   //-----------------------------------------------------------------
1796   // This function tries to find a track prolongation to next pad row
1797   //-----------------------------------------------------------------
1798   t.SetCurrentCluster(0);
1799   t.SetCurrentClusterIndex1(-3);   
1800    
1801   Double_t xt=t.GetX();
1802   Int_t     row = GetRowNumber(xt)-1; 
1803   Double_t  ymax= GetMaxY(nr);
1804
1805   if (row < nr) return 1; // don't prolongate if not information until now -
1806 //   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1807 //     t.fRemoval =10;
1808 //     return 0;  // not prolongate strongly inclined tracks
1809 //   } 
1810 //   if (TMath::Abs(t.GetSnp())>0.95) {
1811 //     t.fRemoval =10;
1812 //     return 0;  // not prolongate strongly inclined tracks
1813 //   }// patch 28 fev 06
1814
1815   Double_t x= GetXrow(nr);
1816   Double_t y,z;
1817   //t.PropagateTo(x+0.02);
1818   //t.PropagateTo(x+0.01);
1819   if (!t.PropagateTo(x)){
1820     return 0;
1821   }
1822   //
1823   y=t.GetY();
1824   z=t.GetZ();
1825
1826   if (TMath::Abs(y)>ymax){
1827     if (y > ymax) {
1828       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1829       if (!t.Rotate(fSectors->GetAlpha())) 
1830         return 0;
1831     } else if (y <-ymax) {
1832       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1833       if (!t.Rotate(-fSectors->GetAlpha())) 
1834         return 0;
1835     }
1836     //    if (!t.PropagateTo(x)){
1837     //  return 0;
1838     //}
1839     return 1;
1840     //y = t.GetY();    
1841   }
1842   //
1843   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1844
1845   if (!IsActive(t.GetRelativeSector(),nr)) {
1846     t.SetInDead(kTRUE);
1847     t.SetClusterIndex2(nr,-1); 
1848     return 0;
1849   }
1850   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1851
1852   AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1853
1854   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1855     t.SetInDead(kTRUE);
1856     t.SetClusterIndex2(nr,-1); 
1857     return 0;
1858   } 
1859   else
1860     {
1861
1862       //      if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) 
1863       if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1864       else
1865         return 0;      
1866     }
1867
1868   // update current
1869   if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1870     //    t.fCurrentSigmaY = GetSigmaY(&t);
1871     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1872     GetShape(&t,nr);
1873   }
1874     
1875   AliTPCclusterMI *cl=0;
1876   Int_t index=0;
1877   //
1878   Double_t roady = 1.;
1879   Double_t roadz = 1.;
1880   //
1881
1882   if (!cl){
1883     index = t.GetClusterIndex2(nr);    
1884     if ( (index>0) && (index&0x8000)==0){
1885       cl = t.GetClusterPointer(nr);
1886       if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1887       t.SetCurrentClusterIndex1(index);
1888       if (cl) {
1889         t.SetCurrentCluster(cl);
1890         return 1;
1891       }
1892     }
1893   }
1894
1895   //  if (index<0) return 0;
1896   UInt_t uindex = TMath::Abs(index);
1897
1898   if (krow) {    
1899     //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);      
1900     cl = krow.FindNearest2(y,z,roady,roadz,uindex);      
1901   }
1902
1903   if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));   
1904   t.SetCurrentCluster(cl);
1905
1906   return 1;
1907 }
1908
1909
1910 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1911   //-----------------------------------------------------------------
1912   // This function tries to find a track prolongation to next pad row
1913   //-----------------------------------------------------------------
1914
1915   //update error according neighborhoud
1916
1917   if (t.GetCurrentCluster()) {
1918     t.SetRow(nr); 
1919     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1920     
1921     if (t.GetCurrentCluster()->IsUsed(10)){
1922       //
1923       //
1924       //  t.fErrorZ2*=2;
1925       //  t.fErrorY2*=2;
1926       t.SetNShared(t.GetNShared()+1);
1927       if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1928         t.SetRemoval(10);
1929         return 0;
1930       }
1931     }   
1932     if (fIteration>0) accept = 0;
1933     if (accept<3)  UpdateTrack(&t,accept);  
1934  
1935   } else {
1936     if (fIteration==0){
1937       if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);      
1938       if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);      
1939
1940       if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1941     }
1942   }
1943   return 1;
1944 }
1945
1946
1947
1948 //_____________________________________________________________________________
1949 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1950   //-----------------------------------------------------------------
1951   // This function tries to find a track prolongation.
1952   //-----------------------------------------------------------------
1953   Double_t xt=t.GetX();
1954   //
1955   Double_t alpha=t.GetAlpha();
1956   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
1957   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
1958   //
1959   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1960     
1961   Int_t first = GetRowNumber(xt);
1962   if (!fromSeeds)
1963     first -= step;
1964   if (first < 0)
1965     first = 0;
1966   for (Int_t nr= first; nr>=rf; nr-=step) {
1967     // update kink info
1968     if (t.GetKinkIndexes()[0]>0){
1969       for (Int_t i=0;i<3;i++){
1970         Int_t index = t.GetKinkIndexes()[i];
1971         if (index==0) break;
1972         if (index<0) continue;
1973         //
1974         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1975         if (!kink){
1976           printf("PROBLEM\n");
1977         }
1978         else{
1979           Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1980           if (kinkrow==nr){
1981             AliExternalTrackParam paramd(t);
1982             kink->SetDaughter(paramd);
1983             kink->SetStatus(2,5);
1984             kink->Update();
1985           }
1986         }
1987       }
1988     }
1989
1990     if (nr==80) t.UpdateReference();
1991     if (nr<fInnerSec->GetNRows()) 
1992       fSectors = fInnerSec;
1993     else
1994       fSectors = fOuterSec;
1995     if (FollowToNext(t,nr)==0) 
1996       if (!t.IsActive()) 
1997         return 0;
1998     
1999   }   
2000   return 1;
2001 }
2002
2003
2004
2005
2006
2007
2008 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2009   //-----------------------------------------------------------------
2010   // This function tries to find a track prolongation.
2011   //-----------------------------------------------------------------
2012   //
2013   Double_t xt=t.GetX();  
2014   Double_t alpha=t.GetAlpha();
2015   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
2016   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
2017   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha())%fN);
2018     
2019   Int_t first = t.GetFirstPoint();
2020   Int_t ri = GetRowNumber(xt); 
2021   if (!fromSeeds)
2022     ri += 1;
2023
2024   if (first<ri) first = ri;
2025   //
2026   if (first<0) first=0;
2027   for (Int_t nr=first; nr<=rf; nr++) {
2028     //    if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2029     if (t.GetKinkIndexes()[0]<0){
2030       for (Int_t i=0;i<3;i++){
2031         Int_t index = t.GetKinkIndexes()[i];
2032         if (index==0) break;
2033         if (index>0) continue;
2034         index = TMath::Abs(index);
2035         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2036         if (!kink){
2037           printf("PROBLEM\n");
2038         }
2039         else{
2040           Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2041           if (kinkrow==nr){
2042             AliExternalTrackParam paramm(t);
2043             kink->SetMother(paramm);
2044             kink->SetStatus(2,1);
2045             kink->Update();
2046           }
2047         }
2048       }      
2049     }
2050     //
2051     if (nr<fInnerSec->GetNRows()) 
2052       fSectors = fInnerSec;
2053     else
2054       fSectors = fOuterSec;
2055
2056     FollowToNext(t,nr);                                                             
2057   }   
2058   return 1;
2059 }
2060
2061
2062
2063
2064    
2065 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2066 {
2067   // overlapping factor
2068   //
2069   sum1=0;
2070   sum2=0;
2071   Int_t sum=0;
2072   //
2073   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2074   dz2*=dz2;  
2075
2076   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2077   dy2*=dy2;
2078   Float_t distance = TMath::Sqrt(dz2+dy2);
2079   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2080  
2081   //  Int_t offset =0;
2082   Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2083   Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2084   if (lastpoint>160) 
2085     lastpoint =160;
2086   if (firstpoint<0) 
2087     firstpoint = 0;
2088   if (firstpoint>lastpoint) {
2089     firstpoint =lastpoint;
2090     //    lastpoint  =160;
2091   }
2092     
2093   
2094   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2095     if (s1->GetClusterIndex2(i)>0) sum1++;
2096     if (s2->GetClusterIndex2(i)>0) sum2++;
2097     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2098       sum++;
2099     }
2100   }
2101   if (sum<5) return 0;
2102
2103   Float_t summin = TMath::Min(sum1+1,sum2+1);
2104   Float_t ratio = (sum+1)/Float_t(summin);
2105   return ratio;
2106 }
2107
2108 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2109 {
2110   // shared clusters
2111   //
2112   Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2113   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2114   Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2115   Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2116   
2117   //
2118   Int_t sumshared=0;
2119   //
2120   //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2121   //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2122   Int_t firstpoint = 0;
2123   Int_t lastpoint = 160;
2124   //
2125   //  if (firstpoint>=lastpoint-5) return;;
2126
2127   for (Int_t i=firstpoint;i<lastpoint;i++){
2128     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2129     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2130       sumshared++;
2131       s1->SetSharedMapBit(i, kTRUE);
2132       s2->SetSharedMapBit(i, kTRUE);
2133     }
2134     if (s1->GetClusterIndex2(i)>0)
2135       s1->SetClusterMapBit(i, kTRUE);
2136   }
2137   if (sumshared>cutN0){
2138     // sign clusters
2139     //
2140     for (Int_t i=firstpoint;i<lastpoint;i++){
2141       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2142       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2143         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2144         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);; 
2145         if (s1->IsActive()&&s2->IsActive()){
2146           p1->SetShared(kTRUE);
2147           p2->SetShared(kTRUE);
2148         }       
2149       }
2150     }
2151   }
2152   //  
2153   if (sumshared>cutN0){
2154     for (Int_t i=0;i<4;i++){
2155       if (s1->GetOverlapLabel(3*i)==0){
2156         s1->SetOverlapLabel(3*i,  s2->GetLabel());
2157         s1->SetOverlapLabel(3*i+1,sumshared);
2158         s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2159         break;
2160       } 
2161     }
2162     for (Int_t i=0;i<4;i++){
2163       if (s2->GetOverlapLabel(3*i)==0){
2164         s2->SetOverlapLabel(3*i,  s1->GetLabel());
2165         s2->SetOverlapLabel(3*i+1,sumshared);
2166         s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2167         break;
2168       } 
2169     }    
2170   }
2171 }
2172
2173 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
2174 {
2175   //
2176   //sort trackss according sectors
2177   //  
2178   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2179     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2180     if (!pt) continue;
2181     //if (pt) RotateToLocal(pt);
2182     pt->SetSort(0);
2183   }
2184   arr->UnSort();
2185   arr->Sort();  // sorting according relative sectors
2186   arr->Expand(arr->GetEntries());
2187   //
2188   //
2189   Int_t nseed=arr->GetEntriesFast();
2190   for (Int_t i=0; i<nseed; i++) {
2191     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2192     if (!pt) continue;
2193     for (Int_t j=0;j<12;j++){
2194       pt->SetOverlapLabel(j,0);
2195     }
2196   }
2197   for (Int_t i=0; i<nseed; i++) {
2198     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2199     if (!pt) continue;
2200     if (pt->GetRemoval()>10) continue;
2201     for (Int_t j=i+1; j<nseed; j++){
2202       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2203       if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2204       //      if (pt2){
2205       if (pt2->GetRemoval()<=10) {
2206         //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2207         SignShared(pt,pt2);
2208       }
2209     }  
2210   }
2211 }
2212
2213
2214 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2215 {
2216   //
2217   //sort tracks in array according mode criteria
2218   Int_t nseed = arr->GetEntriesFast();    
2219   for (Int_t i=0; i<nseed; i++) {
2220     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2221     if (!pt) {
2222       continue;
2223     }
2224     pt->SetSort(mode);
2225   }
2226   arr->UnSort();
2227   arr->Sort();
2228 }
2229
2230
2231 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2232 {
2233   //
2234   // Loop over all tracks and remove overlaped tracks (with lower quality)
2235   // Algorithm:
2236   //    1. Unsign clusters
2237   //    2. Sort tracks according quality
2238   //       Quality is defined by the number of cluster between first and last points 
2239   //       
2240   //    3. Loop over tracks - decreasing quality order
2241   //       a.) remove - If the fraction of shared cluster less than factor (1- n or 2) 
2242   //       b.) remove - If the minimal number of clusters less than minimal and not ITS
2243   //       c.) if track accepted - sign clusters
2244   //
2245   //Called in - AliTPCtrackerMI::Clusters2Tracks()
2246   //          - AliTPCtrackerMI::PropagateBack()
2247   //          - AliTPCtrackerMI::RefitInward()
2248   //
2249   // Arguments:
2250   //   factor1 - factor for constrained
2251   //   factor2 -        for non constrained tracks 
2252   //            if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2253   //
2254   UnsignClusters();
2255   //
2256   Int_t nseed = arr->GetEntriesFast();  
2257   Float_t * quality = new Float_t[nseed];
2258   Int_t   * indexes = new Int_t[nseed];
2259   Int_t good =0;
2260   //
2261   //
2262   for (Int_t i=0; i<nseed; i++) {
2263     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2264     if (!pt){
2265       quality[i]=-1;
2266       continue;
2267     }
2268     pt->UpdatePoints();    //select first last max dens points
2269     Float_t * points = pt->GetPoints();
2270     if (points[3]<0.8) quality[i] =-1;
2271     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2272     //prefer high momenta tracks if overlaps
2273     quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5); 
2274   }
2275   TMath::Sort(nseed,quality,indexes);
2276   //
2277   //
2278   for (Int_t itrack=0; itrack<nseed; itrack++) {
2279     Int_t trackindex = indexes[itrack];
2280     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);   
2281     if (!pt) continue;
2282     //
2283     if (quality[trackindex]<0){
2284         delete arr->RemoveAt(trackindex);
2285         continue;
2286     }
2287     //
2288     //
2289     Int_t first = Int_t(pt->GetPoints()[0]);
2290     Int_t last  = Int_t(pt->GetPoints()[2]);
2291     Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2292     //
2293     Int_t found,foundable,shared;
2294     pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE); // better to get statistic in "high-dens"  region do't use full track as in line bellow
2295     //    pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2296     Bool_t itsgold =kFALSE;
2297     if (pt->GetESD()){
2298       Int_t dummy[12];
2299       if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2300     }
2301     if (!itsgold){
2302       //
2303       if (Float_t(shared+1)/Float_t(found+1)>factor){
2304         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2305         if( AliTPCReconstructor::StreamLevel()>3){
2306           TTreeSRedirector &cstream = *fDebugStreamer;
2307           cstream<<"RemoveUsed"<<
2308             "iter="<<fIteration<<
2309             "pt.="<<pt<<
2310             "\n";
2311         }
2312         delete arr->RemoveAt(trackindex);
2313         continue;
2314       }      
2315       if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
2316         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2317         if( AliTPCReconstructor::StreamLevel()>3){
2318           TTreeSRedirector &cstream = *fDebugStreamer;
2319           cstream<<"RemoveShort"<<
2320             "iter="<<fIteration<<
2321             "pt.="<<pt<<
2322             "\n";
2323         }
2324         delete arr->RemoveAt(trackindex);
2325         continue;
2326       }
2327     }
2328
2329     good++;
2330     //if (sharedfactor>0.4) continue;
2331     if (pt->GetKinkIndexes()[0]>0) continue;
2332     //Remove tracks with undefined properties - seems
2333     if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ? 
2334     //
2335     for (Int_t i=first; i<last; i++) {
2336       Int_t index=pt->GetClusterIndex2(i);
2337       // if (index<0 || index&0x8000 ) continue;
2338       if (index<0 || index&0x8000 ) continue;
2339       AliTPCclusterMI *c= pt->GetClusterPointer(i);        
2340       if (!c) continue;
2341       c->Use(10);  
2342     }    
2343   }
2344   fNtracks = good;
2345   if (fDebug>0){
2346     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2347   }
2348   delete []quality;
2349   delete []indexes;
2350 }
2351
2352 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray) 
2353 {
2354   //
2355   // Dump clusters after reco
2356   // signed and unsigned cluster can be visualized   
2357   // 1. Unsign all cluster
2358   // 2. Sign all used clusters
2359   // 3. Dump clusters
2360   UnsignClusters();
2361   Int_t nseed = trackArray->GetEntries();
2362   for (Int_t i=0; i<nseed; i++){
2363     AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);    
2364     if (!pt) {
2365       continue;
2366     }    
2367     Bool_t isKink=pt->GetKinkIndex(0)!=0;
2368     for (Int_t j=0; j<160; ++j) {
2369       Int_t index=pt->GetClusterIndex2(j);
2370       if (index<0) continue;
2371       AliTPCclusterMI *c= pt->GetClusterPointer(j);
2372       if (!c) continue;
2373       if (isKink) c->Use(100);   // kink
2374       c->Use(10);                // by default usage 10
2375     }
2376   }
2377   //
2378
2379   for (Int_t sec=0;sec<fkNIS;sec++){
2380     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2381       AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2382       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){    
2383         Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
2384         (*fDebugStreamer)<<"clDump"<< 
2385           "iter="<<iter<<
2386           "cl.="<<&cl[icl]<<      
2387           "gx0="<<gx[0]<<
2388           "gx1="<<gx[1]<<
2389           "gx2="<<gx[2]<<
2390           "\n";
2391       }
2392       cl = fInnerSec[sec][row].GetClusters2();
2393       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2394         Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
2395         (*fDebugStreamer)<<"clDump"<< 
2396           "iter="<<iter<<
2397           "cl.="<<&cl[icl]<<
2398           "gx0="<<gx[0]<<
2399           "gx1="<<gx[1]<<
2400           "gx2="<<gx[2]<<
2401           "\n";
2402       }
2403     }
2404   }
2405   
2406   for (Int_t sec=0;sec<fkNOS;sec++){
2407     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2408       AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2409       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2410         Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
2411         (*fDebugStreamer)<<"clDump"<< 
2412           "iter="<<iter<<
2413           "cl.="<<&cl[icl]<<
2414           "gx0="<<gx[0]<<
2415           "gx1="<<gx[1]<<
2416           "gx2="<<gx[2]<<
2417           "\n";      
2418       }
2419       cl = fOuterSec[sec][row].GetClusters2();
2420       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2421         Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
2422         (*fDebugStreamer)<<"clDump"<< 
2423           "iter="<<iter<<
2424           "cl.="<<&cl[icl]<<
2425           "gx0="<<gx[0]<<
2426           "gx1="<<gx[1]<<
2427           "gx2="<<gx[2]<<
2428           "\n";      
2429       }
2430     }
2431   }
2432   
2433 }
2434 void AliTPCtrackerMI::UnsignClusters() 
2435 {
2436   //
2437   // loop over all clusters and unsign them
2438   //
2439   
2440   for (Int_t sec=0;sec<fkNIS;sec++){
2441     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2442       AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2443       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2444         //      if (cl[icl].IsUsed(10))         
2445         cl[icl].Use(-1);
2446       cl = fInnerSec[sec][row].GetClusters2();
2447       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2448         //if (cl[icl].IsUsed(10))       
2449           cl[icl].Use(-1);      
2450     }
2451   }
2452   
2453   for (Int_t sec=0;sec<fkNOS;sec++){
2454     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2455       AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2456       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2457         //if (cl[icl].IsUsed(10))       
2458           cl[icl].Use(-1);
2459       cl = fOuterSec[sec][row].GetClusters2();
2460       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2461         //if (cl[icl].IsUsed(10))       
2462         cl[icl].Use(-1);      
2463     }
2464   }
2465   
2466 }
2467
2468
2469
2470 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2471 {
2472   //
2473   //sign clusters to be "used"
2474   //
2475   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2476   // loop over "primaries"
2477   
2478   Float_t sumdens=0;
2479   Float_t sumdens2=0;
2480   Float_t sumn   =0;
2481   Float_t sumn2  =0;
2482   Float_t sumchi =0;
2483   Float_t sumchi2 =0;
2484
2485   Float_t sum    =0;
2486
2487   TStopwatch timer;
2488   timer.Start();
2489
2490   Int_t nseed = arr->GetEntriesFast();
2491   for (Int_t i=0; i<nseed; i++) {
2492     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2493     if (!pt) {
2494       continue;
2495     }    
2496     if (!(pt->IsActive())) continue;
2497     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2498     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2499       sumdens += dens;
2500       sumdens2+= dens*dens;
2501       sumn    += pt->GetNumberOfClusters();
2502       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2503       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2504       if (chi2>5) chi2=5;
2505       sumchi  +=chi2;
2506       sumchi2 +=chi2*chi2;
2507       sum++;
2508     }
2509   }
2510
2511   Float_t mdensity = 0.9;
2512   Float_t meann    = 130;
2513   Float_t meanchi  = 1;
2514   Float_t sdensity = 0.1;
2515   Float_t smeann    = 10;
2516   Float_t smeanchi  =0.4;
2517   
2518
2519   if (sum>20){
2520     mdensity = sumdens/sum;
2521     meann    = sumn/sum;
2522     meanchi  = sumchi/sum;
2523     //
2524     sdensity = sumdens2/sum-mdensity*mdensity;
2525     if (sdensity >= 0)
2526        sdensity = TMath::Sqrt(sdensity);
2527     else
2528        sdensity = 0.1;
2529     //
2530     smeann   = sumn2/sum-meann*meann;
2531     if (smeann >= 0)
2532       smeann   = TMath::Sqrt(smeann);
2533     else 
2534       smeann   = 10;
2535     //
2536     smeanchi = sumchi2/sum - meanchi*meanchi;
2537     if (smeanchi >= 0)
2538       smeanchi = TMath::Sqrt(smeanchi);
2539     else
2540       smeanchi = 0.4;
2541   }
2542
2543
2544   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2545   //
2546   for (Int_t i=0; i<nseed; i++) {
2547     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2548     if (!pt) {
2549       continue;
2550     }
2551     if (pt->GetBSigned()) continue;
2552     if (pt->GetBConstrain()) continue;    
2553     //if (!(pt->IsActive())) continue;
2554     /*
2555     Int_t found,foundable,shared;    
2556     pt->GetClusterStatistic(0,160,found, foundable,shared);
2557     if (shared/float(found)>0.3) {
2558       if (shared/float(found)>0.9 ){
2559         //delete arr->RemoveAt(i);
2560       }
2561       continue;
2562     }
2563     */
2564     Bool_t isok =kFALSE;
2565     if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2566       isok = kTRUE;
2567     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2568       isok =kTRUE;
2569     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2570       isok =kTRUE;
2571     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2572       isok =kTRUE;
2573     
2574     if (isok)     
2575       for (Int_t j=0; j<160; ++j) {     
2576         Int_t index=pt->GetClusterIndex2(j);
2577         if (index<0) continue;
2578         AliTPCclusterMI *c= pt->GetClusterPointer(j);
2579         if (!c) continue;
2580         //if (!(c->IsUsed(10))) c->Use();  
2581         c->Use(10);  
2582       }
2583   }
2584   
2585   
2586   //
2587   Double_t maxchi  = meanchi+2.*smeanchi;
2588
2589   for (Int_t i=0; i<nseed; i++) {
2590     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2591     if (!pt) {
2592       continue;
2593     }    
2594     //if (!(pt->IsActive())) continue;
2595     if (pt->GetBSigned()) continue;
2596     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2597     if (chi>maxchi) continue;
2598
2599     Float_t bfactor=1;
2600     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2601    
2602     //sign only tracks with enoug big density at the beginning
2603     
2604     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue; 
2605     
2606     
2607     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2608     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2609    
2610     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2611     if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2612       minn=0;
2613
2614     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2615       //Int_t noc=pt->GetNumberOfClusters();
2616       pt->SetBSigned(kTRUE);
2617       for (Int_t j=0; j<160; ++j) {
2618
2619         Int_t index=pt->GetClusterIndex2(j);
2620         if (index<0) continue;
2621         AliTPCclusterMI *c= pt->GetClusterPointer(j);
2622         if (!c) continue;
2623         //      if (!(c->IsUsed(10))) c->Use();  
2624         c->Use(10);  
2625       }
2626     }
2627   }
2628   //  gLastCheck = nseed;
2629   //  arr->Compress();
2630   if (fDebug>0){
2631     timer.Print();
2632   }
2633 }
2634
2635
2636 void  AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2637 {
2638   // stop not active tracks
2639   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2640   // take th2 as threshold for number of founded to number of foundable on last 20 active rows 
2641   Int_t nseed = arr->GetEntriesFast();  
2642   //
2643   for (Int_t i=0; i<nseed; i++) {
2644     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2645     if (!pt) {
2646       continue;
2647     }
2648     if (!(pt->IsActive())) continue;
2649     StopNotActive(pt,row0,th0, th1,th2);
2650   }
2651 }
2652
2653
2654
2655 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2656  Float_t th2) const
2657 {
2658   // stop not active tracks
2659   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2660   // take th2 as threshold for number of founded to number of foundable on last 20 active rows 
2661   Int_t sumgood1  = 0;
2662   Int_t sumgood2  = 0;
2663   Int_t foundable = 0;
2664   Int_t maxindex = seed->GetLastPoint();  //last foundable row
2665   if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2666     seed->Desactivate(10) ;
2667     return;
2668   }
2669
2670   for (Int_t i=row0; i<maxindex; i++){
2671     Int_t index = seed->GetClusterIndex2(i);
2672     if (index!=-1) foundable++;
2673     //if (!c) continue;
2674     if (foundable<=30) sumgood1++;
2675     if (foundable<=50) {
2676       sumgood2++;
2677     }
2678     else{ 
2679       break;
2680     }        
2681   }
2682   if (foundable>=30.){ 
2683      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2684   }
2685   if (foundable>=50)
2686     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2687 }
2688
2689
2690 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2691 {
2692   //
2693   // back propagation of ESD tracks
2694   //
2695   //return 0;
2696   if (!event) return 0;
2697   const Int_t kMaxFriendTracks=2000;
2698   fEvent = event;
2699   // extract correction object for multiplicity dependence of dEdx
2700   TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2701
2702   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2703   if (!transform) {
2704     AliFatal("Tranformations not in RefitInward");
2705     return 0;
2706   }
2707   transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2708   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2709   Int_t nContribut = event->GetNumberOfTracks();
2710   TGraphErrors * graphMultDependenceDeDx = 0x0;
2711   if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2712     if (recoParam->GetUseTotCharge()) {
2713       graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2714     } else {
2715       graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2716     }
2717   }
2718   //
2719   ReadSeeds(event,2);
2720   fIteration=2;
2721   //PrepareForProlongation(fSeeds,1);
2722   PropagateForward2(fSeeds);
2723   RemoveUsed2(fSeeds,0.4,0.4,20);
2724
2725   TObjArray arraySeed(fSeeds->GetEntries());
2726   for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2727     arraySeed.AddAt(fSeeds->At(i),i);    
2728   }
2729   SignShared(&arraySeed);
2730   //  FindCurling(fSeeds, event,2); // find multi found tracks
2731   FindSplitted(fSeeds, event,2); // find multi found tracks
2732   if (AliTPCReconstructor::StreamLevel()>5)  FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2733
2734   Int_t ntracks=0;
2735   Int_t nseed = fSeeds->GetEntriesFast();
2736   for (Int_t i=0;i<nseed;i++){
2737     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2738     if (!seed) continue;
2739     if (seed->GetKinkIndex(0)>0)  UpdateKinkQualityD(seed);  // update quality informations for kinks
2740     AliESDtrack *esd=event->GetTrack(i);
2741
2742     if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2743       AliExternalTrackParam paramIn;
2744       AliExternalTrackParam paramOut;
2745       Int_t ncl = seed->RefitTrack(seed,&paramIn,&paramOut);
2746       if (AliTPCReconstructor::StreamLevel()>2) {
2747         (*fDebugStreamer)<<"RecoverIn"<<
2748           "seed.="<<seed<<
2749           "esd.="<<esd<<
2750           "pin.="<<&paramIn<<
2751           "pout.="<<&paramOut<<
2752           "ncl="<<ncl<<
2753           "\n";
2754       }
2755       if (ncl>15) {
2756         seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2757         seed->SetNumberOfClusters(ncl);
2758       }
2759     }
2760
2761     seed->PropagateTo(fkParam->GetInnerRadiusLow());
2762     seed->UpdatePoints();
2763     AddCovariance(seed);
2764     MakeBitmaps(seed);
2765     seed->CookdEdx(0.02,0.6);
2766     CookLabel(seed,0.1); //For comparison only
2767     esd->SetTPCClusterMap(seed->GetClusterMap());
2768     esd->SetTPCSharedMap(seed->GetSharedMap());
2769     //
2770     if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2771       TTreeSRedirector &cstream = *fDebugStreamer;
2772       cstream<<"Crefit"<<
2773         "Esd.="<<esd<<
2774         "Track.="<<seed<<
2775         "\n"; 
2776     }
2777
2778     if (seed->GetNumberOfClusters()>15){
2779       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit); 
2780       esd->SetTPCPoints(seed->GetPoints());
2781       esd->SetTPCPointsF(seed->GetNFoundable());
2782       Int_t   ndedx = seed->GetNCDEDX(0);
2783       Float_t sdedx = seed->GetSDEDX(0);
2784       Float_t dedx  = seed->GetdEdx();
2785       // apply mutliplicity dependent dEdx correction if available
2786       if (graphMultDependenceDeDx) {
2787         Double_t corrGain =  AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2788         dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2789       }
2790       esd->SetTPCsignal(dedx, sdedx, ndedx);
2791       //
2792       // add seed to the esd track in Calib level
2793       //
2794       Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2795       if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2796         AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE); 
2797         esd->AddCalibObject(seedCopy);
2798       }
2799       ntracks++;
2800     }
2801     else{
2802       //printf("problem\n");
2803     }
2804   }
2805   //FindKinks(fSeeds,event);
2806   if (AliTPCReconstructor::StreamLevel()>3)  DumpClusters(2,fSeeds);
2807   Info("RefitInward","Number of refitted tracks %d",ntracks);
2808   return 0;
2809 }
2810
2811
2812 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2813 {
2814   //
2815   // back propagation of ESD tracks
2816   //
2817   if (!event) return 0;
2818   fEvent = event;
2819   fIteration = 1;
2820   ReadSeeds(event,1);
2821   PropagateBack(fSeeds); 
2822   RemoveUsed2(fSeeds,0.4,0.4,20);
2823   //FindCurling(fSeeds, fEvent,1);  
2824   FindSplitted(fSeeds, event,1); // find multi found tracks
2825   if (AliTPCReconstructor::StreamLevel()>5)  FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2826
2827   //
2828   Int_t nseed = fSeeds->GetEntriesFast();
2829   Int_t ntracks=0;
2830   for (Int_t i=0;i<nseed;i++){
2831     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2832     if (!seed) continue;
2833     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
2834     seed->UpdatePoints();
2835     AddCovariance(seed);
2836     AliESDtrack *esd=event->GetTrack(i);
2837     if (!esd) continue; //never happen
2838     if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2839       AliExternalTrackParam paramIn;
2840       AliExternalTrackParam paramOut;
2841       Int_t ncl = seed->RefitTrack(seed,&paramIn,&paramOut);
2842       if (AliTPCReconstructor::StreamLevel()>2) {
2843         (*fDebugStreamer)<<"RecoverBack"<<
2844           "seed.="<<seed<<
2845           "esd.="<<esd<<
2846           "pin.="<<&paramIn<<
2847           "pout.="<<&paramOut<<
2848           "ncl="<<ncl<<
2849           "\n";
2850       }
2851       if (ncl>15) {
2852         seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2853         seed->SetNumberOfClusters(ncl);
2854       }
2855     }
2856     seed->CookdEdx(0.02,0.6);
2857     CookLabel(seed,0.1); //For comparison only
2858     if (seed->GetNumberOfClusters()>15){
2859       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2860       esd->SetTPCPoints(seed->GetPoints());
2861       esd->SetTPCPointsF(seed->GetNFoundable());
2862       Int_t   ndedx = seed->GetNCDEDX(0);
2863       Float_t sdedx = seed->GetSDEDX(0);
2864       Float_t dedx  = seed->GetdEdx();
2865       esd->SetTPCsignal(dedx, sdedx, ndedx);
2866       ntracks++;
2867       Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2868       // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number      
2869       if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2870         (*fDebugStreamer)<<"Cback"<<
2871           "Tr0.="<<seed<<
2872           "esd.="<<esd<<
2873           "EventNrInFile="<<eventnumber<<
2874           "\n";       
2875       }
2876     }
2877   }
2878   if (AliTPCReconstructor::StreamLevel()>3)  DumpClusters(1,fSeeds);
2879   //FindKinks(fSeeds,event);
2880   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2881   fEvent =0;
2882   
2883   return 0;
2884 }
2885
2886
2887 void AliTPCtrackerMI::DeleteSeeds()
2888 {
2889   //
2890   //delete Seeds
2891
2892   Int_t nseed = fSeeds->GetEntriesFast();
2893   for (Int_t i=0;i<nseed;i++){
2894     AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2895     if (seed) delete fSeeds->RemoveAt(i);
2896   }
2897   delete fSeeds;
2898
2899   fSeeds =0;
2900 }
2901
2902 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2903 {
2904   //
2905   //read seeds from the event
2906   
2907   Int_t nentr=event->GetNumberOfTracks();
2908   if (fDebug>0){
2909     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2910   }
2911   if (fSeeds) 
2912     DeleteSeeds();
2913   if (!fSeeds){   
2914     fSeeds = new TObjArray(nentr);
2915   }
2916   UnsignClusters();
2917   //  Int_t ntrk=0;  
2918   for (Int_t i=0; i<nentr; i++) {
2919     AliESDtrack *esd=event->GetTrack(i);
2920     ULong_t status=esd->GetStatus();
2921     if (!(status&AliESDtrack::kTPCin)) continue;
2922     AliTPCtrack t(*esd);
2923     t.SetNumberOfClusters(0);
2924     //    AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2925     AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2926     seed->SetUniqueID(esd->GetID());
2927     AddCovariance(seed);   //add systematic ucertainty
2928     for (Int_t ikink=0;ikink<3;ikink++) {
2929       Int_t index = esd->GetKinkIndex(ikink);
2930       seed->GetKinkIndexes()[ikink] = index;
2931       if (index==0) continue;
2932       index = TMath::Abs(index);
2933       AliESDkink * kink = fEvent->GetKink(index-1);
2934       if (kink&&esd->GetKinkIndex(ikink)<0){
2935         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2936         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,0);
2937       }
2938       if (kink&&esd->GetKinkIndex(ikink)>0){
2939         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2940         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,4);
2941       }
2942
2943     }
2944     if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.); 
2945     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2946     //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2947     //  fSeeds->AddAt(0,i);
2948     //  delete seed;
2949     //  continue;    
2950     //}
2951     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 )  {
2952       Double_t par0[5],par1[5],alpha,x;
2953       esd->GetInnerExternalParameters(alpha,x,par0);
2954       esd->GetExternalParameters(x,par1);
2955       Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2956       Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2957       Double_t trdchi2=0;
2958       if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2959       //reset covariance if suspicious 
2960       if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2961         seed->ResetCovariance(10.);
2962     }
2963
2964     //
2965     //
2966     // rotate to the local coordinate system
2967     //   
2968     fSectors=fInnerSec; fN=fkNIS;    
2969     Double_t alpha=seed->GetAlpha();
2970     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2971     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2972     Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2973     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2974     alpha-=seed->GetAlpha();  
2975     if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2976     if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2977     if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2978       AliWarning(Form("Rotating track over %f",alpha));
2979       if (!seed->Rotate(alpha)) {
2980         delete seed;
2981         continue;
2982       }
2983     }
2984     seed->SetESD(esd);
2985     // sign clusters
2986     if (esd->GetKinkIndex(0)<=0){
2987       for (Int_t irow=0;irow<160;irow++){
2988         Int_t index = seed->GetClusterIndex2(irow);    
2989         if (index>0){ 
2990           //
2991           AliTPCclusterMI * cl = GetClusterMI(index);
2992           seed->SetClusterPointer(irow,cl);
2993           if (cl){
2994             if ((index & 0x8000)==0){
2995               cl->Use(10);  // accepted cluster   
2996             }else{
2997               cl->Use(6);   // close cluster not accepted
2998             }   
2999           }else{
3000             Info("ReadSeeds","Not found cluster");
3001           }
3002         }
3003       }
3004     }
3005     fSeeds->AddAt(seed,i);
3006   }
3007 }
3008
3009
3010
3011 //_____________________________________________________________________________
3012 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3013                                  Float_t deltay, Int_t ddsec) {
3014   //-----------------------------------------------------------------
3015   // This function creates track seeds.
3016   // SEEDING WITH VERTEX CONSTRAIN 
3017   //-----------------------------------------------------------------
3018   // cuts[0]   - fP4 cut
3019   // cuts[1]   - tan(phi)  cut
3020   // cuts[2]   - zvertex cut
3021   // cuts[3]   - fP3 cut
3022   Int_t nin0  = 0;
3023   Int_t nin1  = 0;
3024   Int_t nin2  = 0;
3025   Int_t nin   = 0;
3026   Int_t nout1 = 0;
3027   Int_t nout2 = 0;
3028
3029   Double_t x[5], c[15];
3030   //  Int_t di = i1-i2;
3031   //
3032   AliTPCseed * seed = new AliTPCseed();
3033   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3034   Double_t cs=cos(alpha), sn=sin(alpha);
3035   //
3036   //  Double_t x1 =fOuterSec->GetX(i1);
3037   //Double_t xx2=fOuterSec->GetX(i2);
3038   
3039   Double_t x1 =GetXrow(i1);
3040   Double_t xx2=GetXrow(i2);
3041
3042   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3043
3044   Int_t imiddle = (i2+i1)/2;    //middle pad row index
3045   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3046   const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3047   //
3048   Int_t ns =sec;   
3049
3050   const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3051   Double_t ymax  = GetMaxY(i1)-kr1.GetDeadZone()-1.5;  
3052   Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;  
3053
3054   //
3055   // change cut on curvature if it can't reach this layer
3056   // maximal curvature set to reach it
3057   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3058   if (dvertexmax*0.5*cuts[0]>0.85){
3059     cuts[0] = 0.85/(dvertexmax*0.5+1.);
3060   }
3061   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
3062
3063   //  Int_t ddsec = 1;
3064   if (deltay>0) ddsec = 0; 
3065   // loop over clusters  
3066   for (Int_t is=0; is < kr1; is++) {
3067     //
3068     if (kr1[is]->IsUsed(10)) continue;
3069     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
3070     //if (TMath::Abs(y1)>ymax) continue;
3071
3072     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
3073
3074     // find possible directions    
3075     Float_t anglez = (z1-z3)/(x1-x3); 
3076     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z      
3077     //
3078     //
3079     //find   rotation angles relative to line given by vertex and point 1
3080     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3081     Double_t dvertex  = TMath::Sqrt(dvertex2);
3082     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
3083     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);            
3084     
3085     //
3086     // loop over 2 sectors
3087     Int_t dsec1=-ddsec;
3088     Int_t dsec2= ddsec;
3089     if (y1<0)  dsec2= 0;
3090     if (y1>0)  dsec1= 0;
3091     
3092     Double_t dddz1=0;  // direction of delta inclination in z axis
3093     Double_t dddz2=0;
3094     if ( (z1-z3)>0)
3095       dddz1 =1;    
3096     else
3097       dddz2 =1;
3098     //
3099     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3100       Int_t sec2 = sec + dsec;
3101       // 
3102       //      AliTPCtrackerRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3103       //AliTPCtrackerRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3104       AliTPCtrackerRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
3105       AliTPCtrackerRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3106       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3107       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3108
3109       // rotation angles to p1-p3
3110       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;            
3111       Double_t x2,   y2,   z2; 
3112       //
3113       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3114
3115       //
3116       Double_t dxx0 =  (xx2-x3)*cs13r;
3117       Double_t dyy0 =  (xx2-x3)*sn13r;
3118       for (Int_t js=index1; js < index2; js++) {
3119         const AliTPCclusterMI *kcl = kr2[js];
3120         if (kcl->IsUsed(10)) continue;  
3121         //
3122         //calcutate parameters
3123         //      
3124         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
3125         // stright track
3126         if (TMath::Abs(yy0)<0.000001) continue;
3127         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
3128         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3129         Double_t r02 = (0.25+y0*y0)*dvertex2;   
3130         //curvature (radius) cut
3131         if (r02<r2min) continue;                
3132        
3133         nin0++; 
3134         //
3135         Double_t c0  = 1/TMath::Sqrt(r02);
3136         if (yy0>0) c0*=-1.;     
3137                
3138        
3139         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
3140         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3141         Double_t dfi0   = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3142         Double_t dfi1   = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);  
3143         //
3144         //
3145         Double_t z0  =  kcl->GetZ();  
3146         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
3147         if (TMath::Abs(zzzz2-z0)>0.5) continue;       
3148         nin1++;              
3149         //      
3150         Double_t dip    = (z1-z0)*c0/dfi1;        
3151         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3152         //
3153         y2 = kcl->GetY(); 
3154         if (dsec==0){
3155           x2 = xx2; 
3156           z2 = kcl->GetZ();       
3157         }
3158         else
3159           {
3160             // rotation 
3161             z2 = kcl->GetZ();  
3162             x2= xx2*cs-y2*sn*dsec;
3163             y2=+xx2*sn*dsec+y2*cs;
3164           }
3165         
3166         x[0] = y1;
3167         x[1] = z1;
3168         x[2] = x0;
3169         x[3] = dip;
3170         x[4] = c0;
3171         //
3172         //
3173         // do we have cluster at the middle ?
3174         Double_t ym,zm;
3175         GetProlongation(x1,xm,x,ym,zm);
3176         UInt_t dummy; 
3177         AliTPCclusterMI * cm=0;
3178         if (TMath::Abs(ym)-ymaxm<0){      
3179           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3180           if ((!cm) || (cm->IsUsed(10))) {        
3181             continue;
3182           }
3183         }
3184         else{     
3185           // rotate y1 to system 0
3186           // get state vector in rotated system 
3187           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
3188           Double_t xr2  =  x0*cs+yr1*sn*dsec;
3189           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3190           //
3191           GetProlongation(xx2,xm,xr,ym,zm);
3192           if (TMath::Abs(ym)-ymaxm<0){
3193             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3194             if ((!cm) || (cm->IsUsed(10))) {      
3195               continue;
3196             }
3197           }
3198         }
3199        
3200
3201         Double_t dym = 0;
3202         Double_t dzm = 0;
3203         if (cm){
3204           dym = ym - cm->GetY();
3205           dzm = zm - cm->GetZ();
3206         }
3207         nin2++;
3208
3209
3210         //
3211         //
3212         Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3213         Double_t sy2=kcl->GetSigmaY2()*2.,     sz2=kcl->GetSigmaZ2()*2.;
3214         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3215         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3216         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3217
3218         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3219         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3220         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3221         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3222         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3223         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3224         
3225         Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3226         Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3227         Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3228         Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3229         
3230         c[0]=sy1;
3231         c[1]=0.;       c[2]=sz1;
3232         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3233         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3234                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3235         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3236         c[13]=f30*sy1*f40+f32*sy2*f42;
3237         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3238         
3239         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3240         
3241         UInt_t index=kr1.GetIndex(is);
3242         seed->~AliTPCseed(); // this does not set the pointer to 0...
3243         AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3244         
3245         track->SetIsSeeding(kTRUE);
3246         track->SetSeed1(i1);
3247         track->SetSeed2(i2);
3248         track->SetSeedType(3);
3249
3250        
3251         //if (dsec==0) {
3252           FollowProlongation(*track, (i1+i2)/2,1);
3253           Int_t foundable,found,shared;
3254           track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3255           if ((found<0.55*foundable)  || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3256             seed->Reset();
3257             seed->~AliTPCseed();
3258             continue;
3259           }
3260           //}
3261         
3262         nin++;
3263         FollowProlongation(*track, i2,1);
3264         
3265         
3266         //Int_t rc = 1;
3267         track->SetBConstrain(1);
3268         //      track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3269         track->SetLastPoint(i1);  // first cluster in track position
3270         track->SetFirstPoint(track->GetLastPoint());
3271         
3272         if (track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3273             track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || 
3274             track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3275           seed->Reset();
3276           seed->~AliTPCseed();
3277           continue;
3278         }
3279         nout1++;
3280         // Z VERTEX CONDITION
3281         Double_t zv, bz=GetBz();
3282         if ( !track->GetZAt(0.,bz,zv) ) continue;
3283         if (TMath::Abs(zv-z3)>cuts[2]) {
3284           FollowProlongation(*track, TMath::Max(i2-20,0));
3285           if ( !track->GetZAt(0.,bz,zv) ) continue;
3286           if (TMath::Abs(zv-z3)>cuts[2]){
3287             FollowProlongation(*track, TMath::Max(i2-40,0));
3288             if ( !track->GetZAt(0.,bz,zv) ) continue;
3289             if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3290               // make seed without constrain
3291               AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3292               FollowProlongation(*track2, i2,1);
3293               track2->SetBConstrain(kFALSE);
3294               track2->SetSeedType(1);
3295               arr->AddLast(track2); 
3296               seed->Reset();
3297               seed->~AliTPCseed();
3298               continue;         
3299             }
3300             else{
3301               seed->Reset();
3302               seed->~AliTPCseed();
3303               continue;
3304             
3305             }
3306           }
3307         }
3308       
3309         track->SetSeedType(0);
3310         arr->AddLast(track); 
3311         seed = new AliTPCseed;  
3312         nout2++;
3313         // don't consider other combinations
3314         if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3315           break;
3316       }
3317     }
3318   }
3319   if (fDebug>3){
3320     Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3321   }
3322   delete seed;
3323 }
3324
3325
3326 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3327                                  Float_t deltay) {
3328   
3329
3330
3331   //-----------------------------------------------------------------
3332   // This function creates track seeds.
3333   //-----------------------------------------------------------------
3334   // cuts[0]   - fP4 cut
3335   // cuts[1]   - tan(phi)  cut
3336   // cuts[2]   - zvertex cut
3337   // cuts[3]   - fP3 cut
3338
3339
3340   Int_t nin0  = 0;
3341   Int_t nin1  = 0;
3342   Int_t nin2  = 0;
3343   Int_t nin   = 0;
3344   Int_t nout1 = 0;
3345   Int_t nout2 = 0;
3346   Int_t nout3 =0;
3347   Double_t x[5], c[15];
3348   //
3349   // make temporary seed
3350   AliTPCseed * seed = new AliTPCseed;
3351   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3352   //  Double_t cs=cos(alpha), sn=sin(alpha);
3353   //
3354   //
3355
3356   // first 3 padrows
3357   Double_t x1 = GetXrow(i1-1);
3358   const    AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3359   Double_t y1max  = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;  
3360   //
3361   Double_t x1p = GetXrow(i1);
3362   const    AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3363   //
3364   Double_t x1m = GetXrow(i1-2);
3365   const    AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3366
3367   //
3368   //last 3 padrow for seeding
3369   AliTPCtrackerRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
3370   Double_t    x3   =  GetXrow(i1-7);
3371   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;  
3372   //
3373   AliTPCtrackerRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
3374   Double_t    x3p   = GetXrow(i1-6);
3375   //
3376   AliTPCtrackerRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
3377   Double_t    x3m   = GetXrow(i1-8);
3378
3379   //
3380   //
3381   // middle padrow
3382   Int_t im = i1-4;                           //middle pad row index
3383   Double_t xm         = GetXrow(im);         // radius of middle pad-row
3384   const AliTPCtrackerRow& krm=GetRow(sec,im);   //middle pad -row
3385   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;  
3386   //
3387   //
3388   Double_t deltax  = x1-x3;
3389   Double_t dymax   = deltax*cuts[1];
3390   Double_t dzmax   = deltax*cuts[3];
3391   //
3392   // loop over clusters  
3393   for (Int_t is=0; is < kr1; is++) {
3394     //
3395     if (kr1[is]->IsUsed(10)) continue;
3396     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
3397     //
3398     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge    
3399     // 
3400     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3401     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3402     //    
3403     Double_t y3,   z3;
3404     //
3405     //
3406     UInt_t index;
3407     for (Int_t js=index1; js < index2; js++) {
3408       const AliTPCclusterMI *kcl = kr3[js];
3409       if (kcl->IsUsed(10)) continue;
3410       y3 = kcl->GetY(); 
3411       // apply angular cuts
3412       if (TMath::Abs(y1-y3)>dymax) continue;
3413       x3 = x3; 
3414       z3 = kcl->GetZ(); 
3415       if (TMath::Abs(z1-z3)>dzmax) continue;
3416       //
3417       Double_t angley = (y1-y3)/(x1-x3);
3418       Double_t anglez = (z1-z3)/(x1-x3);
3419       //
3420       Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3421       Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3422       //
3423       Double_t yyym = angley*(xm-x1)+y1;
3424       Double_t zzzm = anglez*(xm-x1)+z1;
3425
3426       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3427       if (!kcm) continue;
3428       if (kcm->IsUsed(10)) continue;
3429       
3430       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3431       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3432       //
3433       //
3434       //
3435       Int_t used  =0;
3436       Int_t found =0;
3437       //
3438       // look around first
3439       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3440                                                       anglez*(x1m-x1)+z1,
3441                                                       erry,errz,index);
3442       //
3443       if (kc1m){
3444         found++;
3445         if (kc1m->IsUsed(10)) used++;
3446       }
3447       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3448                                                       anglez*(x1p-x1)+z1,
3449                                                       erry,errz,index);
3450       //
3451       if (kc1p){
3452         found++;
3453         if (kc1p->IsUsed(10)) used++;
3454       }
3455       if (used>1)  continue;
3456       if (found<1) continue; 
3457
3458       //
3459       // look around last
3460       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3461                                                       anglez*(x3m-x3)+z3,
3462                                                       erry,errz,index);
3463       //
3464       if (kc3m){
3465         found++;
3466         if (kc3m->IsUsed(10)) used++;
3467       }
3468       else 
3469         continue;
3470       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3471                                                       anglez*(x3p-x3)+z3,
3472                                                       erry,errz,index);
3473       //
3474       if (kc3p){
3475         found++;
3476         if (kc3p->IsUsed(10)) used++;
3477       }
3478       else 
3479         continue;
3480       if (used>1)  continue;
3481       if (found<3) continue;       
3482       //
3483       Double_t x2,y2,z2;
3484       x2 = xm;
3485       y2 = kcm->GetY();
3486       z2 = kcm->GetZ();
3487       //
3488                         
3489       x[0]=y1;
3490       x[1]=z1;
3491       x[4]=F1(x1,y1,x2,y2,x3,y3);
3492       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3493       nin0++;
3494       //
3495       x[2]=F2(x1,y1,x2,y2,x3,y3);
3496       nin1++;
3497       //
3498       x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3499       //if (TMath::Abs(x[3]) > cuts[3]) continue;
3500       nin2++;
3501       //
3502       //
3503       Double_t sy1=0.1,  sz1=0.1;
3504       Double_t sy2=0.1,  sz2=0.1;
3505       Double_t sy3=0.1,  sy=0.1, sz=0.1;
3506       
3507       Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3508       Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3509       Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3510       Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3511       Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3512       Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3513       
3514       Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3515       Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3516       Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3517       Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3518       
3519       c[0]=sy1;
3520       c[1]=0.;       c[2]=sz1; 
3521       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3522       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3523       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3524       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3525       c[13]=f30*sy1*f40+f32*sy2*f42;
3526       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3527       
3528       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3529       
3530       index=kr1.GetIndex(is);
3531       seed->~AliTPCseed();
3532       AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3533       
3534       track->SetIsSeeding(kTRUE);
3535
3536       nin++;      
3537       FollowProlongation(*track, i1-7,1);
3538       if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 || 
3539           track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3540         seed->Reset();
3541         seed->~AliTPCseed();
3542         continue;
3543       }
3544       nout1++;
3545       nout2++;  
3546       //Int_t rc = 1;
3547       FollowProlongation(*track, i2,1);
3548       track->SetBConstrain(0);
3549       track->SetLastPoint(i1+fInnerSec->GetNRows());  // first cluster in track position
3550       track->SetFirstPoint(track->GetLastPoint());
3551       
3552       if (track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3553           track->GetNumberOfClusters()<track->GetNFoundable()*0.7 || 
3554           track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3555         seed->Reset();
3556         seed->~AliTPCseed();
3557         continue;
3558       }
3559    
3560       {
3561         FollowProlongation(*track, TMath::Max(i2-10,0),1);
3562         AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3563         FollowProlongation(*track2, i2,1);
3564         track2->SetBConstrain(kFALSE);
3565         track2->SetSeedType(4);
3566         arr->AddLast(track2); 
3567         seed->Reset();
3568         seed->~AliTPCseed();
3569       }
3570       
3571    
3572       //arr->AddLast(track); 
3573       //seed = new AliTPCseed;  
3574       nout3++;
3575     }
3576   }
3577   
3578   if (fDebug>3){
3579     Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3580   }
3581   delete seed;
3582 }
3583
3584
3585 //_____________________________________________________________________________
3586 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3587                                  Float_t deltay, Bool_t /*bconstrain*/) {
3588   //-----------------------------------------------------------------
3589   // This function creates track seeds - without vertex constraint
3590   //-----------------------------------------------------------------
3591   // cuts[0]   - fP4 cut        - not applied
3592   // cuts[1]   - tan(phi)  cut
3593   // cuts[2]   - zvertex cut    - not applied 
3594   // cuts[3]   - fP3 cut
3595   Int_t nin0=0;
3596   Int_t nin1=0;
3597   Int_t nin2=0;
3598   Int_t nin3=0;
3599   //  Int_t nin4=0;
3600   //Int_t nin5=0;
3601
3602   
3603
3604   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3605   //  Double_t cs=cos(alpha), sn=sin(alpha);
3606   Int_t row0 = (i1+i2)/2;
3607   Int_t drow = (i1-i2)/2;
3608   const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3609   AliTPCtrackerRow * kr=0;
3610
3611   AliTPCpolyTrack polytrack;
3612   Int_t nclusters=fSectors[sec][row0];
3613   AliTPCseed * seed = new AliTPCseed;
3614
3615   Int_t sumused=0;
3616   Int_t cused=0;
3617   Int_t cnused=0;
3618   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3619     Int_t nfound =0;
3620     Int_t nfoundable =0;
3621     for (Int_t iter =1; iter<2; iter++){   //iterations
3622       const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3623       const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];      
3624       const AliTPCclusterMI * cl= kr0[is];
3625       
3626       if (cl->IsUsed(10)) {
3627         cused++;
3628       }
3629       else{
3630         cnused++;
3631       }
3632       Double_t x = kr0.GetX();
3633       // Initialization of the polytrack
3634       nfound =0;
3635       nfoundable =0;
3636       polytrack.Reset();
3637       //
3638       Double_t y0= cl->GetY();
3639       Double_t z0= cl->GetZ();
3640       Float_t erry = 0;
3641       Float_t errz = 0;
3642       
3643       Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3644       if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;  // seed only at the edge
3645       
3646       erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;      
3647       errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;      
3648       polytrack.AddPoint(x,y0,z0,erry, errz);
3649
3650       sumused=0;
3651       if (cl->IsUsed(10)) sumused++;
3652
3653
3654       Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3655       Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3656       //
3657       x = krm.GetX();
3658       AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3659       if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3660         erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;          
3661         errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3662         if (cl1->IsUsed(10))  sumused++;
3663         polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3664       }
3665       //
3666       x = krp.GetX();
3667       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3668       if (cl2) {
3669         erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;          
3670         errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3671         if (cl2->IsUsed(10)) sumused++;  
3672         polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3673       }
3674       //
3675       if (sumused>0) continue;
3676       nin0++;
3677       polytrack.UpdateParameters();
3678       // follow polytrack
3679       roadz = 1.2;
3680       roady = 1.2;
3681       //
3682       Double_t yn,zn;
3683       nfoundable = polytrack.GetN();
3684       nfound     = nfoundable; 
3685       //
3686       for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3687         Float_t maxdist = 0.8*(1.+3./(ddrow));
3688         for (Int_t delta = -1;delta<=1;delta+=2){
3689           Int_t row = row0+ddrow*delta;
3690           kr = &(fSectors[sec][row]);
3691           Double_t xn = kr->GetX();
3692           Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3693           polytrack.GetFitPoint(xn,yn,zn);
3694           if (TMath::Abs(yn)>ymax1) continue;
3695           nfoundable++;
3696           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3697           if (cln) {
3698             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3699             if (dist<maxdist){
3700               /*
3701               erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));         
3702               errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3703               if (cln->IsUsed(10)) {
3704                 //      printf("used\n");
3705                 sumused++;
3706                 erry*=2;
3707                 errz*=2;
3708               }
3709               */
3710               erry=0.1;
3711               errz=0.1;
3712               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3713               nfound++;
3714             }
3715           }
3716         }
3717         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;     
3718         polytrack.UpdateParameters();
3719       }           
3720     }
3721     if ( (sumused>3) || (sumused>0.5*nfound))  {
3722       //printf("sumused   %d\n",sumused);
3723       continue;
3724     }
3725     nin1++;
3726     Double_t dy,dz;
3727     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3728     AliTPCpolyTrack track2;
3729     
3730     polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3731     if (track2.GetN()<0.5*nfoundable) continue;
3732     nin2++;
3733
3734     if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3735       //
3736       // test seed with and without constrain
3737       for (Int_t constrain=0; constrain<=0;constrain++){
3738         // add polytrack candidate
3739
3740         Double_t x[5], c[15];
3741         Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3742         track2.GetBoundaries(x3,x1);    
3743         x2 = (x1+x3)/2.;
3744         track2.GetFitPoint(x1,y1,z1);
3745         track2.GetFitPoint(x2,y2,z2);
3746         track2.GetFitPoint(x3,y3,z3);
3747         //
3748         //is track pointing to the vertex ?
3749         Double_t x0,y0,z0;
3750         x0=0;
3751         polytrack.GetFitPoint(x0,y0,z0);
3752
3753         if (constrain) {
3754           x2 = x3;
3755           y2 = y3;
3756           z2 = z3;
3757           
3758           x3 = 0;
3759           y3 = 0;
3760           z3 = 0;
3761         }
3762         x[0]=y1;
3763         x[1]=z1;
3764         x[4]=F1(x1,y1,x2,y2,x3,y3);
3765                 
3766         //      if (TMath::Abs(x[4]) >= cuts[0]) continue;  //
3767         x[2]=F2(x1,y1,x2,y2,x3,y3);
3768         
3769         //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3770         //x[3]=F3(x1,y1,x2,y2,z1,z2);
3771         x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3772         //if (TMath::Abs(x[3]) > cuts[3]) continue;
3773
3774         
3775         Double_t sy =0.1, sz =0.1;
3776         Double_t sy1=0.02, sz1=0.02;
3777         Double_t sy2=0.02, sz2=0.02;
3778         Double_t sy3=0.02;
3779
3780         if (constrain){
3781           sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3782         }
3783         
3784         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3785         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3786         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3787         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3788         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3789         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3790
3791         Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3792         Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3793         Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3794         Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3795
3796         
3797         c[0]=sy1;
3798         c[1]=0.;       c[2]=sz1;
3799         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3800         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3801         c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3802         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3803         c[13]=f30*sy1*f40+f32*sy2*f42;
3804         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3805         
3806         //Int_t row1 = fSectors->GetRowNumber(x1);
3807         Int_t row1 = GetRowNumber(x1);
3808
3809         UInt_t index=0;
3810         //kr0.GetIndex(is);
3811         seed->~AliTPCseed();
3812         AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3813         track->SetIsSeeding(kTRUE);
3814         Int_t rc=FollowProlongation(*track, i2);        
3815         if (constrain) track->SetBConstrain(1);
3816         else
3817           track->SetBConstrain(0);
3818         track->SetLastPoint(row1+fInnerSec->GetNRows());  // first cluster in track position
3819         track->SetFirstPoint(track->GetLastPoint());
3820
3821         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3822             track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || 
3823             track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3824           //delete track;
3825           seed->Reset();
3826           seed->~AliTPCseed();
3827         }
3828         else {
3829           arr->AddLast(track);
3830           seed = new AliTPCseed;
3831         }
3832         nin3++;
3833       }
3834     }  // if accepted seed
3835   }
3836   if (fDebug>3){
3837     Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3838   }
3839   delete seed;
3840 }
3841
3842
3843 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3844 {
3845   //
3846   //
3847   //reseed using track points
3848   Int_t p0 = int(r0*track->GetNumberOfClusters());     // point 0 
3849   Int_t p1 = int(r1*track->GetNumberOfClusters());
3850   Int_t p2 = int(r2*track->GetNumberOfClusters());   // last point
3851   Int_t pp2=0;
3852   Double_t  x0[3],x1[3],x2[3];
3853   for (Int_t i=0;i<3;i++){
3854     x0[i]=-1;
3855     x1[i]=-1;
3856     x2[i]=-1;
3857   }
3858
3859   // find track position at given ratio of the length
3860   Int_t  sec0=0, sec1=0, sec2=0;
3861   Int_t index=-1;
3862   Int_t clindex;
3863   for (Int_t i=0;i<160;i++){
3864     if (track->GetClusterPointer(i)){
3865       index++;
3866       AliTPCTrackerPoint   *trpoint =track->GetTrackPoint(i);
3867       if ( (index<p0) || x0[0]<0 ){
3868         if (trpoint->GetX()>1){
3869           clindex = track->GetClusterIndex2(i);
3870           if (clindex>0){       
3871             x0[0] = trpoint->GetX();
3872             x0[1] = trpoint->GetY();
3873             x0[2] = trpoint->GetZ();
3874             sec0  = ((clindex&0xff000000)>>24)%18;
3875           }
3876         }
3877       }
3878
3879       if ( (index<p1) &&(trpoint->GetX()>1)){
3880         clindex = track->GetClusterIndex2(i);
3881         if (clindex>0){
3882           x1[0] = trpoint->GetX();
3883           x1[1] = trpoint->GetY();
3884           x1[2] = trpoint->GetZ();
3885           sec1  = ((clindex&0xff000000)>>24)%18;
3886         }
3887       }
3888       if ( (index<p2) &&(trpoint->GetX()>1)){
3889         clindex = track->GetClusterIndex2(i);
3890         if (clindex>0){
3891           x2[0] = trpoint->GetX();
3892           x2[1] = trpoint->GetY();
3893           x2[2] = trpoint->GetZ(); 
3894           sec2  = ((clindex&0xff000000)>>24)%18;
3895           pp2 = i;
3896         }
3897       }
3898     }
3899   }
3900   
3901   Double_t alpha, cs,sn, xx2,yy2;
3902   //
3903   alpha = (sec1-sec2)*fSectors->GetAlpha();
3904   cs = TMath::Cos(alpha);
3905   sn = TMath::Sin(alpha); 
3906   xx2= x1[0]*cs-x1[1]*sn;
3907   yy2= x1[0]*sn+x1[1]*cs;
3908   x1[0] = xx2;
3909   x1[1] = yy2;
3910   //
3911   alpha = (sec0-sec2)*fSectors->GetAlpha();
3912   cs = TMath::Cos(alpha);
3913   sn = TMath::Sin(alpha); 
3914   xx2= x0[0]*cs-x0[1]*sn;
3915   yy2= x0[0]*sn+x0[1]*cs;
3916   x0[0] = xx2;
3917   x0[1] = yy2;
3918   //
3919   //
3920   //
3921   Double_t x[5],c[15];
3922   //
3923   x[0]=x2[1];
3924   x[1]=x2[2];
3925   x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3926   //  if (x[4]>1) return 0;
3927   x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3928   x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3929   //if (TMath::Abs(x[3]) > 2.2)  return 0;
3930   //if (TMath::Abs(x[2]) > 1.99) return 0;
3931   //  
3932   Double_t sy =0.1,  sz =0.1;
3933   //
3934   Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3935   Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3936   Double_t sy3=0.01+track->GetSigmaY2();
3937   //
3938   Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3939   Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3940   Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3941   Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3942   Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3943   Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3944   //
3945   Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3946   Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3947   Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3948   Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3949   
3950   
3951   c[0]=sy1;
3952   c[1]=0.;       c[2]=sz1;
3953   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3954   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3955   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3956   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3957   c[13]=f30*sy1*f40+f32*sy2*f42;
3958   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3959   
3960   //  Int_t row1 = fSectors->GetRowNumber(x2[0]);
3961   AliTPCseed *seed=new  AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3962   //  Double_t y0,z0,y1,z1, y2,z2;
3963   //seed->GetProlongation(x0[0],y0,z0);
3964   // seed->GetProlongation(x1[0],y1,z1);
3965   //seed->GetProlongation(x2[0],y2,z2);
3966   //  seed =0;
3967   seed->SetLastPoint(pp2);
3968   seed->SetFirstPoint(pp2);
3969   
3970
3971   return seed;
3972 }
3973
3974
3975 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3976 {
3977   //
3978   //
3979   //reseed using founded clusters 
3980   //
3981   // Find the number of clusters
3982   Int_t nclusters = 0;
3983   for (Int_t irow=0;irow<160;irow++){
3984     if (track->GetClusterIndex(irow)>0) nclusters++;
3985   }
3986   //
3987   Int_t ipos[3];
3988   ipos[0] = TMath::Max(int(r0*nclusters),0);             // point 0 cluster
3989   ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1);   // 
3990   ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1);   // last point
3991   //
3992   //
3993   Double_t  xyz[3][3]={{0}};
3994   Int_t     row[3]={0},sec[3]={0,0,0};
3995   //
3996   // find track row position at given ratio of the length
3997   Int_t index=-1;
3998   for (Int_t irow=0;irow<160;irow++){    
3999     if (track->GetClusterIndex2(irow)<0) continue;
4000     index++;
4001     for (Int_t ipoint=0;ipoint<3;ipoint++){
4002       if (index<=ipos[ipoint]) row[ipoint] = irow;
4003     }        
4004   }
4005   //
4006   //Get cluster and sector position
4007   for (Int_t ipoint=0;ipoint<3;ipoint++){
4008     Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4009     AliTPCclusterMI * cl = GetClusterMI(clindex);
4010     if (cl==0) {
4011       //Error("Bug\n");
4012       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
4013       return 0;
4014     }
4015     sec[ipoint]     = ((clindex&0xff000000)>>24)%18;
4016     xyz[ipoint][0]  = GetXrow(row[ipoint]);
4017     xyz[ipoint][1]  = cl->GetY();
4018     xyz[ipoint][2]  = cl->GetZ();
4019   }
4020   //
4021   //
4022   // Calculate seed state vector and covariance matrix
4023
4024   Double_t alpha, cs,sn, xx2,yy2;
4025   //
4026   alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4027   cs = TMath::Cos(alpha);
4028   sn = TMath::Sin(alpha); 
4029   xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4030   yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4031   xyz[1][0] = xx2;
4032   xyz[1][1] = yy2;
4033   //
4034   alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4035   cs = TMath::Cos(alpha);
4036   sn = TMath::Sin(alpha); 
4037   xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4038   yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4039   xyz[0][0] = xx2;
4040   xyz[0][1] = yy2;
4041   //
4042   //
4043   //
4044   Double_t x[5],c[15];
4045   //
4046   x[0]=xyz[2][1];
4047   x[1]=xyz[2][2];
4048   x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4049   x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4050   x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4051   //  
4052   Double_t sy =0.1,  sz =0.1;
4053   //
4054   Double_t sy1=0.2, sz1=0.2;
4055   Double_t sy2=0.2, sz2=0.2;
4056   Double_t sy3=0.2;
4057   //
4058   Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
4059   Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
4060   Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
4061   Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
4062   Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
4063   Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
4064   //
4065   Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
4066   Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
4067   Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
4068   Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
4069   
4070   
4071   c[0]=sy1;
4072   c[1]=0.;       c[2]=sz1;
4073   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4074   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
4075   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4076   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4077   c[13]=f30*sy1*f40+f32*sy2*f42;
4078   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4079   
4080   //  Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4081   AliTPCseed *seed=new  AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4082   seed->SetLastPoint(row[2]);
4083   seed->SetFirstPoint(row[2]);  
4084   return seed;
4085 }
4086
4087
4088 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4089 {
4090   //
4091   //
4092   //reseed using founded clusters 
4093   //
4094   Double_t  xyz[3][3];
4095   Int_t     row[3]={0,0,0};
4096   Int_t     sec[3]={0,0,0};
4097   //
4098   // forward direction
4099   if (forward){
4100     for (Int_t irow=r0;irow<160;irow++){
4101       if (track->GetClusterIndex(irow)>0){
4102         row[0] = irow;
4103         break;
4104       }
4105     }
4106     for (Int_t irow=160;irow>r0;irow--){
4107       if (track->GetClusterIndex(irow)>0){
4108         row[2] = irow;
4109         break;
4110       }
4111     }
4112     for (Int_t irow=row[2]-15;irow>row[0];irow--){
4113       if (track->GetClusterIndex(irow)>0){
4114         row[1] = irow;
4115         break;
4116       }
4117     }
4118     //
4119   }
4120   if (!forward){
4121     for (Int_t irow=0;irow<r0;irow++){
4122       if (track->GetClusterIndex(irow)>0){
4123         row[0] = irow;
4124         break;
4125       }
4126     }
4127     for (Int_t irow=r0;irow>0;irow--){
4128       if (track->GetClusterIndex(irow)>0){
4129         row[2] = irow;
4130         break;
4131       }
4132     }    
4133     for (Int_t irow=row[2]-15;irow>row[0];irow--){
4134       if (track->GetClusterIndex(irow)>0){
4135         row[1] = irow;
4136         break;
4137       }
4138     } 
4139   }
4140   //
4141   if ((row[2]-row[0])<20) return 0;
4142   if (row[1]==0) return 0;
4143   //
4144   //
4145   //Get cluster and sector position
4146   for (Int_t ipoint=0;ipoint<3;ipoint++){
4147     Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4148     AliTPCclusterMI * cl = GetClusterMI(clindex);
4149     if (cl==0) {
4150       //Error("Bug\n");
4151       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
4152       return 0;
4153     }
4154     sec[ipoint]     = ((clindex&0xff000000)>>24)%18;
4155     xyz[ipoint][0]  = GetXrow(row[ipoint]);
4156     AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);    
4157     if (point&&ipoint<2){
4158       //
4159        xyz[ipoint][1]  = point->GetY();
4160        xyz[ipoint][2]  = point->GetZ();
4161     }
4162     else{
4163       xyz[ipoint][1]  = cl->GetY();
4164       xyz[ipoint][2]  = cl->GetZ();
4165     }
4166   }
4167   //
4168   //
4169   //
4170   //
4171   // Calculate seed state vector and covariance matrix
4172
4173   Double_t alpha, cs,sn, xx2,yy2;
4174   //
4175   alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4176   cs = TMath::Cos(alpha);
4177   sn = TMath::Sin(alpha); 
4178   xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4179   yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4180   xyz[1][0] = xx2;
4181   xyz[1][1] = yy2;
4182   //
4183   alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4184   cs = TMath::Cos(alpha);
4185   sn = TMath::Sin(alpha); 
4186   xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4187   yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4188   xyz[0][0] = xx2;
4189   xyz[0][1] = yy2;
4190   //
4191   //
4192   //
4193   Double_t x[5],c[15];
4194   //
4195   x[0]=xyz[2][1];
4196   x[1]=xyz[2][2];
4197   x[4]=F1(xyz[2][0],xyz[2