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