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