5e63f8cd7e504c2bf1fc261bafcb90ba6f031c11
[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 //  How to use?  - 
24 //  run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
25 //  run AliTPCFindTracksMI.C macro - to find tracks
26 //  tracks are written to AliTPCtracks.root file
27 //  for comparison also seeds are written to the same file - to special branch
28 //-------------------------------------------------------
29
30
31 /* $Id$ */
32
33
34
35 #include "Riostream.h"
36 #include <TClonesArray.h>
37 #include <TFile.h>
38 #include <TObjArray.h>
39 #include <TTree.h>
40
41 #include "AliComplexCluster.h"
42 #include "AliESD.h"
43 #include "AliHelix.h"
44 #include "AliRunLoader.h"
45 #include "AliTPCClustersRow.h"
46 #include "AliTPCParam.h"
47 #include "AliTPCclusterMI.h"
48 #include "AliTPCpolyTrack.h"
49 #include "AliTPCreco.h" 
50 #include "AliTPCtrackerMI.h"
51 #include "TStopwatch.h"
52
53 #include "AliTPCReconstructor.h"
54 #include "AliESDkink.h"
55 //
56
57 ClassImp(AliTPCseed)
58 ClassImp(AliTPCtrackerMI)
59
60
61 class AliTPCFastMath {
62 public:
63   AliTPCFastMath();  
64   static Double_t FastAsin(Double_t x);   
65  private: 
66   static Double_t fgFastAsin[20000];  //lookup table for fast asin computation
67 };
68
69 Double_t AliTPCFastMath::fgFastAsin[20000];
70 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
71
72 AliTPCFastMath::AliTPCFastMath(){
73   //
74   // initialized lookup table;
75   for (Int_t i=0;i<10000;i++){
76     fgFastAsin[2*i] = TMath::ASin(i/10000.);
77     fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
78   }
79 }
80
81 Double_t AliTPCFastMath::FastAsin(Double_t x){
82   //
83   // return asin using lookup table
84   if (x>0){
85     Int_t index = int(x*10000);
86     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
87   }
88   x*=-1;
89   Int_t index = int(x*10000);
90   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
91 }
92
93
94
95
96 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
97   //
98   //update track information using current cluster - track->fCurrentCluster
99
100
101   AliTPCclusterMI* c =track->fCurrentCluster;
102   if (accept>0) track->fCurrentClusterIndex1 |=0x8000;  //sign not accepted clusters
103
104   UInt_t i = track->fCurrentClusterIndex1;
105
106   Int_t sec=(i&0xff000000)>>24; 
107   //Int_t row = (i&0x00ff0000)>>16; 
108   track->fRow=(i&0x00ff0000)>>16;
109   track->fSector = sec;
110   //  Int_t index = i&0xFFFF;
111   if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow(); 
112   track->SetClusterIndex2(track->fRow, i);  
113   //track->fFirstPoint = row;
114   //if ( track->fLastPoint<row) track->fLastPoint =row;
115   //  if (track->fRow<0 || track->fRow>160) {
116   //  printf("problem\n");
117   //}
118   if (track->fFirstPoint>track->fRow) 
119     track->fFirstPoint = track->fRow;
120   if (track->fLastPoint<track->fRow) 
121     track->fLastPoint  = track->fRow;
122   
123
124   track->fClusterPointer[track->fRow] = c;  
125   //
126
127   Float_t angle2 = track->GetSnp()*track->GetSnp();
128   angle2 = TMath::Sqrt(angle2/(1-angle2)); 
129   //
130   //SET NEW Track Point
131   //
132   //  if (debug)
133   {
134     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->fRow));
135     //
136     point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
137     point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
138     point.SetErrY(sqrt(track->fErrorY2));
139     point.SetErrZ(sqrt(track->fErrorZ2));
140     //
141     point.SetX(track->GetX());
142     point.SetY(track->GetY());
143     point.SetZ(track->GetZ());
144     point.SetAngleY(angle2);
145     point.SetAngleZ(track->GetTgl());
146     if (point.fIsShared){
147       track->fErrorY2 *= 4;
148       track->fErrorZ2 *= 4;
149     }
150   }  
151
152   Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
153   //
154   track->fErrorY2 *= 1.3;
155   track->fErrorY2 += 0.01;    
156   track->fErrorZ2 *= 1.3;   
157   track->fErrorZ2 += 0.005;      
158     //}
159   if (accept>0) return 0;
160   if (track->GetNumberOfClusters()%20==0){
161     //    if (track->fHelixIn){
162     //  TClonesArray & larr = *(track->fHelixIn);    
163     //  Int_t ihelix = larr.GetEntriesFast();
164     //  new(larr[ihelix]) AliHelix(*track) ;    
165     //}
166   }
167   track->fNoCluster =0;
168   return track->Update(c,chi2,i);
169 }
170
171
172
173 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor, 
174                                       Float_t cory, Float_t corz)
175 {
176   //
177   // decide according desired precision to accept given 
178   // cluster for tracking
179   Double_t sy2=ErrY2(seed,cluster)*cory;
180   Double_t sz2=ErrZ2(seed,cluster)*corz;
181   //sy2=ErrY2(seed,cluster)*cory;
182   //sz2=ErrZ2(seed,cluster)*cory;
183   
184   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
185   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
186   
187   Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
188     (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
189   Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
190     (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
191   
192   Double_t rdistance2  = rdistancey2+rdistancez2;
193   //Int_t  accept =0;
194   
195   if (rdistance2>16) return 3;
196   
197   
198   if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)  
199     return 2;  //suspisiouce - will be changed
200   
201   if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)  
202     // strict cut on overlaped cluster
203     return  2;  //suspisiouce - will be changed
204   
205   if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor ) 
206        && cluster->GetType()<0){
207     seed->fNFoundable--;
208     return 2;    
209   }
210   return 0;
211 }
212
213
214
215
216 //_____________________________________________________________________________
217 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par): 
218 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
219 {
220   //---------------------------------------------------------------------
221   // The main TPC tracker constructor
222   //---------------------------------------------------------------------
223   fInnerSec=new AliTPCSector[fkNIS];         
224   fOuterSec=new AliTPCSector[fkNOS];
225  
226   Int_t i;
227   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
228   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
229
230   fN=0;  fSectors=0;
231
232   fSeeds=0;
233   fNtracks = 0;
234   fParam = par;  
235   Int_t nrowlow = par->GetNRowLow();
236   Int_t nrowup = par->GetNRowUp();
237
238   
239   for (Int_t i=0;i<nrowlow;i++){
240     fXRow[i]     = par->GetPadRowRadiiLow(i);
241     fPadLength[i]= par->GetPadPitchLength(0,i);
242     fYMax[i]     = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
243   }
244
245   
246   for (Int_t i=0;i<nrowup;i++){
247     fXRow[i+nrowlow]      = par->GetPadRowRadiiUp(i);
248     fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
249     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
250   }
251   fSeeds=0;
252   //
253   fInput    = 0;
254   fOutput   = 0;
255   fSeedTree = 0;
256   fTreeDebug =0;
257   fNewIO     =0;
258   fDebug     =0;
259   fEvent     =0;
260 }
261 //________________________________________________________________________
262 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
263   AliTracker(t),
264   fkNIS(t.fkNIS),
265   fkNOS(t.fkNOS)
266 {
267   //------------------------------------
268   // dummy copy constructor
269   //------------------------------------------------------------------
270 }
271 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
272   //------------------------------
273   // dummy 
274   //--------------------------------------------------------------
275   return *this;
276 }
277 //_____________________________________________________________________________
278 AliTPCtrackerMI::~AliTPCtrackerMI() {
279   //------------------------------------------------------------------
280   // TPC tracker destructor
281   //------------------------------------------------------------------
282   delete[] fInnerSec;
283   delete[] fOuterSec;
284   if (fSeeds) {
285     fSeeds->Delete(); 
286     delete fSeeds;
287   }
288 }
289
290 void AliTPCtrackerMI::SetIO()
291 {
292   //
293   fNewIO   =  kTRUE;
294   fInput   =  AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
295   
296   fOutput  =  AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
297   if (fOutput){
298     AliTPCtrack *iotrack= new AliTPCtrack;
299     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
300     delete iotrack;
301   }
302 }
303
304
305 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
306 {
307
308   // set input
309   fNewIO = kFALSE;
310   fInput    = 0;
311   fOutput   = 0;
312   fSeedTree = 0;
313   fTreeDebug =0;
314   fInput = input;
315   if (input==0){
316     return;
317   }  
318   //set output
319   fOutput = output;
320   if (output){
321     AliTPCtrack *iotrack= new AliTPCtrack;
322     //    iotrack->fHelixIn   = new TClonesArray("AliHelix");
323     //iotrack->fHelixOut  = new TClonesArray("AliHelix");    
324     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
325     delete iotrack;
326   }
327   if (output && (fDebug&2)){
328     //write the full seed information if specified in debug mode
329     //
330     fSeedTree =  new TTree("Seeds","Seeds");
331     AliTPCseed * vseed = new AliTPCseed;
332     //
333     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
334     arrtr->ExpandCreateFast(160);
335     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
336     //
337     vseed->fPoints = arrtr;
338     vseed->fEPoints = arre;
339     //    vseed->fClusterPoints = arrcl;
340     fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
341     delete arrtr;
342     delete arre;    
343     fTreeDebug = new TTree("trackDebug","trackDebug");
344     TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
345     fTreeDebug->Branch("debug",&arrd,32000,99);
346   }
347
348
349   //set ESD event  
350   fEvent  = event;  
351 }
352
353 void AliTPCtrackerMI::FillESD(TObjArray* arr)
354 {
355   //
356   //
357   //fill esds using updated tracks
358   if (fEvent){
359     // write tracks to the event
360     // store index of the track
361     Int_t nseed=arr->GetEntriesFast();
362     //FindKinks(arr,fEvent);
363     for (Int_t i=0; i<nseed; i++) {
364       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
365       if (!pt) continue; 
366       pt->UpdatePoints();
367       pt->PropagateTo(fParam->GetInnerRadiusLow());
368  
369       if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
370         AliESDtrack iotrack;
371         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
372         iotrack.SetTPCPoints(pt->GetPoints());
373         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
374         //iotrack.SetTPCindex(i);
375         fEvent->AddTrack(&iotrack);
376         continue;
377       }
378        
379       if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) {
380         AliESDtrack iotrack;
381         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
382         iotrack.SetTPCPoints(pt->GetPoints());
383         //iotrack.SetTPCindex(i);
384         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
385         fEvent->AddTrack(&iotrack);
386         continue;
387       } 
388       //
389       // short tracks  - maybe decays
390
391       if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.70) {
392         Int_t found,foundable,shared;
393         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
394         if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2)){
395           AliESDtrack iotrack;
396           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);    
397           //iotrack.SetTPCindex(i);
398           iotrack.SetTPCPoints(pt->GetPoints());
399           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
400           fEvent->AddTrack(&iotrack);
401           continue;
402         }
403       }       
404       
405       if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.8) {
406         Int_t found,foundable,shared;
407         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
408         if (found<20) continue;
409         if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
410         //
411         AliESDtrack iotrack;
412         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
413         iotrack.SetTPCPoints(pt->GetPoints());
414         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
415         //iotrack.SetTPCindex(i);
416         fEvent->AddTrack(&iotrack);
417         continue;
418       }   
419       // short tracks  - secondaties
420       //
421       if ( (pt->GetNumberOfClusters()>30) ) {
422         Int_t found,foundable,shared;
423         pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
424         if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
425           AliESDtrack iotrack;
426           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);    
427           iotrack.SetTPCPoints(pt->GetPoints());
428           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
429           //iotrack.SetTPCindex(i);
430           fEvent->AddTrack(&iotrack);
431           continue;
432         }
433       }       
434       
435       if ( (pt->GetNumberOfClusters()>15)) {
436         Int_t found,foundable,shared;
437         pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
438         if (found<15) continue;
439         if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
440         if (float(found)/float(foundable)<0.8) continue;
441         //
442         AliESDtrack iotrack;
443         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
444         iotrack.SetTPCPoints(pt->GetPoints());
445         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
446         //iotrack.SetTPCindex(i);
447         fEvent->AddTrack(&iotrack);
448         continue;
449       }   
450     }
451   }
452   printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
453 }
454
455 void AliTPCtrackerMI::WriteTracks(TTree * tree)
456 {
457   //
458   // write tracks from seed array to selected tree
459   //
460   fOutput  = tree;
461   if (fOutput){
462     AliTPCtrack *iotrack= new AliTPCtrack;
463     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
464   }
465   WriteTracks();
466 }
467
468 void AliTPCtrackerMI::WriteTracks()
469 {
470   //
471   // write tracks to the given output tree -
472   // output specified with SetIO routine
473   if (!fSeeds)  return;
474   if (!fOutput){
475     SetIO();
476   }
477
478   if (fOutput){
479     AliTPCtrack *iotrack= 0;
480     Int_t nseed=fSeeds->GetEntriesFast();
481     //for (Int_t i=0; i<nseed; i++) {
482     //  iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
483     //  if (iotrack) break;      
484     //}    
485     //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
486     TBranch * br = fOutput->GetBranch("tracks");
487     br->SetAddress(&iotrack);
488     //
489     for (Int_t i=0; i<nseed; i++) {
490       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);    
491       if (!pt) continue;    
492       AliTPCtrack * track = new AliTPCtrack(*pt);
493       iotrack = track;
494       pt->fLab2 =i; 
495       //      br->SetAddress(&iotrack);
496       fOutput->Fill();
497       delete track;
498       iotrack =0;
499     }
500     //fOutput->GetDirectory()->cd();
501     //fOutput->Write();
502   }
503   // delete iotrack;
504   //
505   if (fSeedTree){
506     //write the full seed information if specified in debug mode
507       
508     AliTPCseed * vseed = new AliTPCseed;
509     //
510     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
511     arrtr->ExpandCreateFast(160);
512     //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
513     //arrcl->ExpandCreateFast(160);
514     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
515     //
516     vseed->fPoints = arrtr;
517     vseed->fEPoints = arre;
518     //    vseed->fClusterPoints = arrcl;
519     //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
520     TBranch * brseed = fSeedTree->GetBranch("seeds");
521     
522     Int_t nseed=fSeeds->GetEntriesFast();
523     
524     for (Int_t i=0; i<nseed; i++) {
525       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);    
526       if (!pt) continue;     
527       pt->fPoints = arrtr;
528       //      pt->fClusterPoints = arrcl;
529       pt->fEPoints       = arre;
530       pt->RebuildSeed();
531       vseed = pt;
532       brseed->SetAddress(&vseed);
533       fSeedTree->Fill();
534       pt->fPoints  = 0;
535       pt->fEPoints = 0;
536       //      pt->fClusterPoints = 0;
537     }
538     fSeedTree->Write();
539     if (fTreeDebug) fTreeDebug->Write();
540   }
541
542 }
543   
544
545
546
547 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
548   //
549   //
550   //seed->SetErrorY2(0.1);
551   //return 0.1;
552   //calculate look-up table at the beginning
553   static Bool_t  ginit = kFALSE;
554   static Float_t gnoise1,gnoise2,gnoise3;
555   static Float_t ggg1[10000];
556   static Float_t ggg2[10000];
557   static Float_t ggg3[10000];
558   static Float_t glandau1[10000];
559   static Float_t glandau2[10000];
560   static Float_t glandau3[10000];
561   //
562   static Float_t gcor01[500];
563   static Float_t gcor02[500];
564   static Float_t gcorp[500];
565   //
566
567   //
568   if (ginit==kFALSE){
569     for (Int_t i=1;i<500;i++){
570       Float_t rsigma = float(i)/100.;
571       gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
572       gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
573       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
574     }
575
576     //
577     for (Int_t i=3;i<10000;i++){
578       //
579       //
580       // inner sector
581       Float_t amp = float(i);
582       Float_t padlength =0.75;
583       gnoise1 = 0.0004/padlength;
584       Float_t nel     = 0.268*amp;
585       Float_t nprim   = 0.155*amp;
586       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
587       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
588       if (glandau1[i]>1) glandau1[i]=1;
589       glandau1[i]*=padlength*padlength/12.;      
590       //
591       // outer short
592       padlength =1.;
593       gnoise2   = 0.0004/padlength;
594       nel       = 0.3*amp;
595       nprim     = 0.133*amp;
596       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
597       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
598       if (glandau2[i]>1) glandau2[i]=1;
599       glandau2[i]*=padlength*padlength/12.;
600       //
601       //
602       // outer long
603       padlength =1.5;
604       gnoise3   = 0.0004/padlength;
605       nel       = 0.3*amp;
606       nprim     = 0.133*amp;
607       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
608       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
609       if (glandau3[i]>1) glandau3[i]=1;
610       glandau3[i]*=padlength*padlength/12.;
611       //
612     }
613     ginit = kTRUE;
614   }
615   //
616   //
617   //
618   Int_t amp = int(TMath::Abs(cl->GetQ()));  
619   if (amp>9999) {
620     seed->SetErrorY2(1.);
621     return 1.;
622   }
623   Float_t snoise2;
624   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
625   Int_t ctype = cl->GetType();  
626   Float_t padlength= GetPadPitchLength(seed->fRow);
627   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
628   angle2 = angle2/(1-angle2); 
629   //
630   //cluster "quality"
631   Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
632   Float_t res;
633   //
634   if (fSectors==fInnerSec){
635     snoise2 = gnoise1;
636     res     = ggg1[amp]*z+glandau1[amp]*angle2;     
637     if (ctype==0) res *= gcor01[rsigmay];
638     if ((ctype>0)){
639       res+=0.002;
640       res*= gcorp[rsigmay];
641     }
642   }
643   else {
644     if (padlength<1.1){
645       snoise2 = gnoise2;
646       res     = ggg2[amp]*z+glandau2[amp]*angle2; 
647       if (ctype==0) res *= gcor02[rsigmay];      
648       if ((ctype>0)){
649         res+=0.002;
650         res*= gcorp[rsigmay];
651       }
652     }
653     else{
654       snoise2 = gnoise3;      
655       res     = ggg3[amp]*z+glandau3[amp]*angle2; 
656       if (ctype==0) res *= gcor02[rsigmay];
657       if ((ctype>0)){
658         res+=0.002;
659         res*= gcorp[rsigmay];
660       }
661     }
662   }  
663
664   if (ctype<0){
665     res+=0.005;
666     res*=2.4;  // overestimate error 2 times
667   }
668   res+= snoise2;
669  
670   if (res<2*snoise2)
671     res = 2*snoise2;
672   
673   seed->SetErrorY2(res);
674   return res;
675
676
677 }
678
679
680
681 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
682   //
683   //
684   //seed->SetErrorY2(0.1);
685   //return 0.1;
686   //calculate look-up table at the beginning
687   static Bool_t  ginit = kFALSE;
688   static Float_t gnoise1,gnoise2,gnoise3;
689   static Float_t ggg1[10000];
690   static Float_t ggg2[10000];
691   static Float_t ggg3[10000];
692   static Float_t glandau1[10000];
693   static Float_t glandau2[10000];
694   static Float_t glandau3[10000];
695   //
696   static Float_t gcor01[1000];
697   static Float_t gcor02[1000];
698   static Float_t gcorp[1000];
699   //
700
701   //
702   if (ginit==kFALSE){
703     for (Int_t i=1;i<1000;i++){
704       Float_t rsigma = float(i)/100.;
705       gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
706       gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
707       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
708     }
709
710     //
711     for (Int_t i=3;i<10000;i++){
712       //
713       //
714       // inner sector
715       Float_t amp = float(i);
716       Float_t padlength =0.75;
717       gnoise1 = 0.0004/padlength;
718       Float_t nel     = 0.268*amp;
719       Float_t nprim   = 0.155*amp;
720       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
721       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
722       if (glandau1[i]>1) glandau1[i]=1;
723       glandau1[i]*=padlength*padlength/12.;      
724       //
725       // outer short
726       padlength =1.;
727       gnoise2   = 0.0004/padlength;
728       nel       = 0.3*amp;
729       nprim     = 0.133*amp;
730       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
731       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
732       if (glandau2[i]>1) glandau2[i]=1;
733       glandau2[i]*=padlength*padlength/12.;
734       //
735       //
736       // outer long
737       padlength =1.5;
738       gnoise3   = 0.0004/padlength;
739       nel       = 0.3*amp;
740       nprim     = 0.133*amp;
741       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
742       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
743       if (glandau3[i]>1) glandau3[i]=1;
744       glandau3[i]*=padlength*padlength/12.;
745       //
746     }
747     ginit = kTRUE;
748   }
749   //
750   //
751   //
752   Int_t amp = int(TMath::Abs(cl->GetQ()));  
753   if (amp>9999) {
754     seed->SetErrorY2(1.);
755     return 1.;
756   }
757   Float_t snoise2;
758   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
759   Int_t ctype = cl->GetType();  
760   Float_t padlength= GetPadPitchLength(seed->fRow);
761   //
762   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
763   //  if (angle2<0.6) angle2 = 0.6;
764   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2)); 
765   //
766   //cluster "quality"
767   Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
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[rsigmaz];
774     if ((ctype>0)){
775       res+=0.002;
776       res*= gcorp[rsigmaz];
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[rsigmaz];      
784       if ((ctype>0)){
785         res+=0.002;
786         res*= gcorp[rsigmaz];
787       }
788     }
789     else{
790       snoise2 = gnoise3;      
791       res     = ggg3[amp]*z+glandau3[amp]*angle2; 
792       if (ctype==0) res *= gcor02[rsigmaz];
793       if ((ctype>0)){
794         res+=0.002;
795         res*= gcorp[rsigmaz];
796       }
797     }
798   }  
799
800   if (ctype<0){
801     res+=0.002;
802     res*=1.3;
803   }
804   if ((ctype<0) &&amp<70){
805     res+=0.002;
806     res*=1.3;  
807   }
808   res += snoise2;
809   if (res<2*snoise2)
810      res = 2*snoise2;
811   if (res>3) res =3;
812   seed->SetErrorZ2(res);
813   return res;
814 }
815
816
817
818 /*
819 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
820   //
821   //
822   //seed->SetErrorZ2(0.1);
823   //return 0.1;
824
825   Float_t snoise2;
826   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
827   //
828   Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
829   Int_t ctype = cl->GetType();
830   Float_t amp = TMath::Abs(cl->GetQ());
831   
832   Float_t nel;
833   Float_t nprim;
834   //
835   Float_t landau=2 ;    //landau fluctuation part
836   Float_t gg=2;         // gg fluctuation part
837   Float_t padlength= GetPadPitchLength(seed->GetX());
838  
839   if (fSectors==fInnerSec){
840     snoise2 = 0.0004/padlength;
841     nel     = 0.268*amp;
842     nprim   = 0.155*amp;
843     gg      = (2+0.001*nel/(padlength*padlength))/nel;
844     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
845     if (landau>1) landau=1;
846   }
847   else {
848     snoise2 = 0.0004/padlength;
849     nel     = 0.3*amp;
850     nprim   = 0.133*amp;
851     gg      = (2+0.0008*nel/(padlength*padlength))/nel;
852     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
853     if (landau>1) landau=1;
854   }
855   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
856
857   //
858   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
859   angle2 = TMath::Sqrt((1-angle2));
860   if (angle2<0.6) angle2 = 0.6;
861   //angle2 = 1;
862
863   Float_t angle = seed->GetTgl()/angle2;
864   Float_t angular = landau*angle*angle*padlength*padlength/12.;
865   Float_t res = sdiff + angular;
866
867   
868   if ((ctype==0) && (fSectors ==fOuterSec))
869     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
870
871   if ((ctype==0) && (fSectors ==fInnerSec))
872     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
873   
874   if ((ctype>0)){
875     res+=0.005;
876     res*= TMath::Power(rsigmaz+0.5,1.5);  //0.31+0.147*ctype;
877   }
878   if (ctype<0){
879     res+=0.002;
880     res*=1.3;
881   }
882   if ((ctype<0) &&amp<70){
883     res+=0.002;
884     res*=1.3;  
885   }
886   res += snoise2;
887   if (res<2*snoise2)
888      res = 2*snoise2;
889
890   seed->SetErrorZ2(res);
891   return res;
892 }
893 */
894
895
896
897 void AliTPCseed::Reset(Bool_t all)
898 {
899   //
900   //
901   SetNumberOfClusters(0);
902   fNFoundable = 0;
903   SetChi2(0);
904   ResetCovariance();
905   /*
906   if (fTrackPoints){
907     for (Int_t i=0;i<8;i++){
908       delete [] fTrackPoints[i];
909     }
910     delete fTrackPoints;
911     fTrackPoints =0;
912   }
913   */
914
915   if (all){   
916     for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
917     for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
918   }
919
920 }
921
922
923 void AliTPCseed::Modify(Double_t factor)
924 {
925
926   //------------------------------------------------------------------
927   //This function makes a track forget its history :)  
928   //------------------------------------------------------------------
929   if (factor<=0) {
930     ResetCovariance();
931     return;
932   }
933   fC00*=factor;
934   fC10*=0;  fC11*=factor;
935   fC20*=0;  fC21*=0;  fC22*=factor;
936   fC30*=0;  fC31*=0;  fC32*=0;  fC33*=factor;
937   fC40*=0;  fC41*=0;  fC42*=0;  fC43*=0;  fC44*=factor;
938   SetNumberOfClusters(0);
939   fNFoundable =0;
940   SetChi2(0);
941   fRemoval = 0;
942   fCurrentSigmaY2 = 0.000005;
943   fCurrentSigmaZ2 = 0.000005;
944   fNoCluster     = 0;
945   //fFirstPoint = 160;
946   //fLastPoint  = 0;
947 }
948
949
950
951
952 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
953 {
954   //-----------------------------------------------------------------
955   // This function find proloncation of a track to a reference plane x=xk.
956   // doesn't change internal state of the track
957   //-----------------------------------------------------------------
958   
959   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
960
961   if (TMath::Abs(fP4*xk - fP2) >= 0.999) {   
962     return 0;
963   }
964
965   //  Double_t y1=fP0, z1=fP1;
966   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
967   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
968   
969   y = fP0;
970   z = fP1;
971   //y += dx*(c1+c2)/(r1+r2);
972   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
973   
974   Double_t dy = dx*(c1+c2)/(r1+r2);
975   Double_t dz = 0;
976   //
977   Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
978   /*
979   if (TMath::Abs(delta)>0.0001){
980     dz = fP3*TMath::ASin(delta)/fP4;
981   }else{
982     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
983   }
984   */
985   dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
986   //
987   y+=dy;
988   z+=dz;
989   
990
991   return 1;  
992 }
993
994
995 //_____________________________________________________________________________
996 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const 
997 {
998   //-----------------------------------------------------------------
999   // This function calculates a predicted chi2 increment.
1000   //-----------------------------------------------------------------
1001   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
1002   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1003   r00+=fC00; r01+=fC10; r11+=fC11;
1004
1005   Double_t det=r00*r11 - r01*r01;
1006   if (TMath::Abs(det) < 1.e-10) {
1007     Int_t n=GetNumberOfClusters();
1008     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
1009     return 1e10;
1010   }
1011   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1012   
1013   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1014   
1015   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
1016 }
1017
1018
1019 //_________________________________________________________________________________________
1020
1021
1022 Int_t AliTPCseed::Compare(const TObject *o) const {
1023   //-----------------------------------------------------------------
1024   // This function compares tracks according to the sector - for given sector according z
1025   //-----------------------------------------------------------------
1026   AliTPCseed *t=(AliTPCseed*)o;
1027
1028   if (fSort == 0){
1029     if (t->fRelativeSector>fRelativeSector) return -1;
1030     if (t->fRelativeSector<fRelativeSector) return 1;
1031     Double_t z2 = t->GetZ();
1032     Double_t z1 = GetZ();
1033     if (z2>z1) return 1;
1034     if (z2<z1) return -1;
1035     return 0;
1036   }
1037   else {
1038     Float_t f2 =1;
1039     f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
1040     if (t->fBConstrain) f2=1.2;
1041
1042     Float_t f1 =1;
1043     f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
1044
1045     if (fBConstrain)   f1=1.2;
1046  
1047     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
1048     else return +1;
1049   }
1050 }
1051
1052 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
1053 {
1054   //rotate to track "local coordinata
1055   Float_t x = seed->GetX();
1056   Float_t y = seed->GetY();
1057   Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
1058   
1059   if (y > ymax) {
1060     seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
1061     if (!seed->Rotate(fSectors->GetAlpha())) 
1062       return;
1063   } else if (y <-ymax) {
1064     seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
1065     if (!seed->Rotate(-fSectors->GetAlpha())) 
1066       return;
1067   }   
1068
1069 }
1070
1071
1072
1073
1074 //_____________________________________________________________________________
1075 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
1076   //-----------------------------------------------------------------
1077   // This function associates a cluster with this track.
1078   //-----------------------------------------------------------------
1079   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1080
1081   r00+=fC00; r01+=fC10; r11+=fC11;
1082   Double_t det=r00*r11 - r01*r01;
1083   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
1084
1085   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
1086   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
1087   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
1088   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
1089   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
1090
1091   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1092   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1093   if (TMath::Abs(cur*fX-eta) >= 0.9) {
1094     return 0;
1095   }
1096
1097   fP0 += k00*dy + k01*dz;
1098   fP1 += k10*dy + k11*dz;
1099   fP2  = eta;
1100   fP3 += k30*dy + k31*dz;
1101   fP4  = cur;
1102
1103   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1104   Double_t c12=fC21, c13=fC31, c14=fC41;
1105
1106   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1107   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
1108   fC40-=k00*c04+k01*c14; 
1109
1110   fC11-=k10*c01+k11*fC11;
1111   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
1112   fC41-=k10*c04+k11*c14; 
1113
1114   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
1115   fC42-=k20*c04+k21*c14; 
1116
1117   fC33-=k30*c03+k31*c13;
1118   fC43-=k40*c03+k41*c13; 
1119
1120   fC44-=k40*c04+k41*c14; 
1121
1122   Int_t n=GetNumberOfClusters();
1123   //  fIndex[n]=index;
1124   SetNumberOfClusters(n+1);
1125   SetChi2(GetChi2()+chisq);
1126
1127   return 1;
1128 }
1129
1130
1131
1132 //_____________________________________________________________________________
1133 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1134                    Double_t x2,Double_t y2,
1135                    Double_t x3,Double_t y3) 
1136 {
1137   //-----------------------------------------------------------------
1138   // Initial approximation of the track curvature
1139   //-----------------------------------------------------------------
1140   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1141   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1142                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1143   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1144                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1145
1146   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1147   if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1148   return -xr*yr/sqrt(xr*xr+yr*yr); 
1149 }
1150
1151
1152
1153 //_____________________________________________________________________________
1154 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1155                    Double_t x2,Double_t y2,
1156                    Double_t x3,Double_t y3) 
1157 {
1158   //-----------------------------------------------------------------
1159   // Initial approximation of the track curvature
1160   //-----------------------------------------------------------------
1161   x3 -=x1;
1162   x2 -=x1;
1163   y3 -=y1;
1164   y2 -=y1;
1165   //  
1166   Double_t det = x3*y2-x2*y3;
1167   if (det==0) {
1168     return 100;
1169   }
1170   //
1171   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1172   Double_t x0 = x3*0.5-y3*u;
1173   Double_t y0 = y3*0.5+x3*u;
1174   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1175   if (det<0) c2*=-1;
1176   return c2;
1177 }
1178
1179
1180 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1181                    Double_t x2,Double_t y2,
1182                    Double_t x3,Double_t y3) 
1183 {
1184   //-----------------------------------------------------------------
1185   // Initial approximation of the track curvature
1186   //-----------------------------------------------------------------
1187   x3 -=x1;
1188   x2 -=x1;
1189   y3 -=y1;
1190   y2 -=y1;
1191   //  
1192   Double_t det = x3*y2-x2*y3;
1193   if (det==0) {
1194     return 100;
1195   }
1196   //
1197   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1198   Double_t x0 = x3*0.5-y3*u; 
1199   Double_t y0 = y3*0.5+x3*u;
1200   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1201   if (det<0) c2*=-1;
1202   x0+=x1;
1203   x0*=c2;  
1204   return x0;
1205 }
1206
1207
1208
1209 //_____________________________________________________________________________
1210 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1211                    Double_t x2,Double_t y2,
1212                    Double_t x3,Double_t y3) 
1213 {
1214   //-----------------------------------------------------------------
1215   // Initial approximation of the track curvature times center of curvature
1216   //-----------------------------------------------------------------
1217   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1218   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1219                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1220   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1221                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1222
1223   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1224   
1225   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1226 }
1227
1228 //_____________________________________________________________________________
1229 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1, 
1230                    Double_t x2,Double_t y2,
1231                    Double_t z1,Double_t z2) 
1232 {
1233   //-----------------------------------------------------------------
1234   // Initial approximation of the tangent of the track dip angle
1235   //-----------------------------------------------------------------
1236   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1237 }
1238
1239
1240 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1, 
1241                    Double_t x2,Double_t y2,
1242                    Double_t z1,Double_t z2, Double_t c) 
1243 {
1244   //-----------------------------------------------------------------
1245   // Initial approximation of the tangent of the track dip angle
1246   //-----------------------------------------------------------------
1247
1248   //  Double_t angle1;
1249   
1250   //angle1    =  (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1251   //
1252   Double_t d  =  TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1253   if (TMath::Abs(d*c*0.5)>1) return 0;
1254   //  Double_t   angle2    =  TMath::ASin(d*c*0.5);
1255   //  Double_t   angle2    =  AliTPCFastMath::FastAsin(d*c*0.5);
1256   Double_t   angle2    = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1257
1258   angle2  = (z1-z2)*c/(angle2*2.);
1259   return angle2;
1260 }
1261
1262 Bool_t   AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1263 {//-----------------------------------------------------------------
1264   // This function find proloncation of a track to a reference plane x=x2.
1265   //-----------------------------------------------------------------
1266   
1267   Double_t dx=x2-x1;
1268
1269   if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {   
1270     return kFALSE;
1271   }
1272
1273   Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1274   Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);  
1275   y = x[0];
1276   z = x[1];
1277   
1278   Double_t dy = dx*(c1+c2)/(r1+r2);
1279   Double_t dz = 0;
1280   //
1281   Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1282   
1283   if (TMath::Abs(delta)>0.01){
1284     dz = x[3]*TMath::ASin(delta)/x[4];
1285   }else{
1286     dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1287   }
1288   
1289   //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1290
1291   y+=dy;
1292   z+=dz;
1293   
1294   return kTRUE;  
1295 }
1296
1297 Int_t  AliTPCtrackerMI::LoadClusters (TTree *tree)
1298 {
1299   //
1300   //
1301   fInput = tree;
1302   return LoadClusters();
1303 }
1304
1305 Int_t  AliTPCtrackerMI::LoadClusters()
1306 {
1307   //
1308   // load clusters to the memory
1309   AliTPCClustersRow *clrow= new AliTPCClustersRow;
1310   clrow->SetClass("AliTPCclusterMI");
1311   clrow->SetArray(0);
1312   clrow->GetArray()->ExpandCreateFast(10000);
1313   //
1314   //  TTree * tree = fClustersArray.GetTree();
1315
1316   TTree * tree = fInput;
1317   TBranch * br = tree->GetBranch("Segment");
1318   br->SetAddress(&clrow);
1319   //
1320   Int_t j=Int_t(tree->GetEntries());
1321   for (Int_t i=0; i<j; i++) {
1322     br->GetEntry(i);
1323     //  
1324     Int_t sec,row;
1325     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1326     //
1327     AliTPCRow * tpcrow=0;
1328     Int_t left=0;
1329     if (sec<fkNIS*2){
1330       tpcrow = &(fInnerSec[sec%fkNIS][row]);    
1331       left = sec/fkNIS;
1332     }
1333     else{
1334       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1335       left = (sec-fkNIS*2)/fkNOS;
1336     }
1337     if (left ==0){
1338       tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1339       tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1340       for (Int_t i=0;i<tpcrow->fN1;i++) 
1341         tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1342     }
1343     if (left ==1){
1344       tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1345       tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1346       for (Int_t i=0;i<tpcrow->fN2;i++) 
1347         tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1348     }
1349   }
1350   //
1351   delete clrow;
1352   LoadOuterSectors();
1353   LoadInnerSectors();
1354   return 0;
1355 }
1356
1357
1358 void AliTPCtrackerMI::UnloadClusters()
1359 {
1360   //
1361   // unload clusters from the memory
1362   //
1363   Int_t nrows = fOuterSec->GetNRows();
1364   for (Int_t sec = 0;sec<fkNOS;sec++)
1365     for (Int_t row = 0;row<nrows;row++){
1366       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1367       //      if (tpcrow){
1368       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1369       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1370       //}
1371       tpcrow->ResetClusters();
1372     }
1373   //
1374   nrows = fInnerSec->GetNRows();
1375   for (Int_t sec = 0;sec<fkNIS;sec++)
1376     for (Int_t row = 0;row<nrows;row++){
1377       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1378       //if (tpcrow){
1379       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1380       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1381       //}
1382       tpcrow->ResetClusters();
1383     }
1384
1385   return ;
1386 }
1387
1388
1389 //_____________________________________________________________________________
1390 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1391   //-----------------------------------------------------------------
1392   // This function fills outer TPC sectors with clusters.
1393   //-----------------------------------------------------------------
1394   Int_t nrows = fOuterSec->GetNRows();
1395   UInt_t index=0;
1396   for (Int_t sec = 0;sec<fkNOS;sec++)
1397     for (Int_t row = 0;row<nrows;row++){
1398       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);  
1399       Int_t sec2 = sec+2*fkNIS;
1400       //left
1401       Int_t ncl = tpcrow->fN1;
1402       while (ncl--) {
1403         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1404         index=(((sec2<<8)+row)<<16)+ncl;
1405         tpcrow->InsertCluster(c,index);
1406       }
1407       //right
1408       ncl = tpcrow->fN2;
1409       while (ncl--) {
1410         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1411         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1412         tpcrow->InsertCluster(c,index);
1413       }
1414       //
1415       // write indexes for fast acces
1416       //
1417       for (Int_t i=0;i<510;i++)
1418         tpcrow->fFastCluster[i]=-1;
1419       for (Int_t i=0;i<tpcrow->GetN();i++){
1420         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1421         tpcrow->fFastCluster[zi]=i;  // write index
1422       }
1423       Int_t last = 0;
1424       for (Int_t i=0;i<510;i++){
1425         if (tpcrow->fFastCluster[i]<0)
1426           tpcrow->fFastCluster[i] = last;
1427         else
1428           last = tpcrow->fFastCluster[i];
1429       }
1430     }  
1431   fN=fkNOS;
1432   fSectors=fOuterSec;
1433   return 0;
1434 }
1435
1436
1437 //_____________________________________________________________________________
1438 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1439   //-----------------------------------------------------------------
1440   // This function fills inner TPC sectors with clusters.
1441   //-----------------------------------------------------------------
1442   Int_t nrows = fInnerSec->GetNRows();
1443   UInt_t index=0;
1444   for (Int_t sec = 0;sec<fkNIS;sec++)
1445     for (Int_t row = 0;row<nrows;row++){
1446       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1447       //
1448       //left
1449       Int_t ncl = tpcrow->fN1;
1450       while (ncl--) {
1451         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1452         index=(((sec<<8)+row)<<16)+ncl;
1453         tpcrow->InsertCluster(c,index);
1454       }
1455       //right
1456       ncl = tpcrow->fN2;
1457       while (ncl--) {
1458         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1459         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1460         tpcrow->InsertCluster(c,index);
1461       }
1462       //
1463       // write indexes for fast acces
1464       //
1465       for (Int_t i=0;i<510;i++)
1466         tpcrow->fFastCluster[i]=-1;
1467       for (Int_t i=0;i<tpcrow->GetN();i++){
1468         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1469         tpcrow->fFastCluster[zi]=i;  // write index
1470       }
1471       Int_t last = 0;
1472       for (Int_t i=0;i<510;i++){
1473         if (tpcrow->fFastCluster[i]<0)
1474           tpcrow->fFastCluster[i] = last;
1475         else
1476           last = tpcrow->fFastCluster[i];
1477       }
1478
1479     }  
1480    
1481   fN=fkNIS;
1482   fSectors=fInnerSec;
1483   return 0;
1484 }
1485
1486
1487
1488 //_________________________________________________________________________
1489 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1490   //--------------------------------------------------------------------
1491   //       Return pointer to a given cluster
1492   //--------------------------------------------------------------------
1493   Int_t sec=(index&0xff000000)>>24; 
1494   Int_t row=(index&0x00ff0000)>>16; 
1495   Int_t ncl=(index&0x00007fff)>>00;
1496
1497   const AliTPCRow * tpcrow=0;
1498   AliTPCclusterMI * clrow =0;
1499
1500   if (sec<fkNIS*2){
1501     tpcrow = &(fInnerSec[sec%fkNIS][row]);
1502     if (tpcrow==0) return 0;
1503
1504     if (sec<fkNIS) {
1505       if (tpcrow->fN1<=ncl) return 0;
1506       clrow = tpcrow->fClusters1;
1507     }
1508     else {
1509       if (tpcrow->fN2<=ncl) return 0;
1510       clrow = tpcrow->fClusters2;
1511     }
1512   }
1513   else {
1514     tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1515     if (tpcrow==0) return 0;
1516
1517     if (sec-2*fkNIS<fkNOS) {
1518       if (tpcrow->fN1<=ncl) return 0;
1519       clrow = tpcrow->fClusters1;
1520     }
1521     else {
1522       if (tpcrow->fN2<=ncl) return 0;
1523       clrow = tpcrow->fClusters2;
1524     }
1525   }
1526
1527   return &(clrow[ncl]);      
1528   
1529 }
1530
1531
1532
1533 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1534   //-----------------------------------------------------------------
1535   // This function tries to find a track prolongation to next pad row
1536   //-----------------------------------------------------------------
1537   //
1538   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1539   AliTPCclusterMI *cl=0;
1540   Int_t tpcindex= t.GetClusterIndex2(nr);
1541   //
1542   // update current shape info every 5 pad-row
1543   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1544     GetShape(&t,nr);    
1545     //}
1546   //  
1547   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters 
1548     //        
1549     if (tpcindex==-1) return 0; //track in dead zone
1550     if (tpcindex>0){     //
1551       cl = t.fClusterPointer[nr];
1552       if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1553       t.fCurrentClusterIndex1 = tpcindex; 
1554     }
1555     if (cl){      
1556       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
1557       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1558       //
1559       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1560       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1561       
1562       if (TMath::Abs(angle-t.GetAlpha())>0.001){
1563         Double_t rotation = angle-t.GetAlpha();
1564         t.fRelativeSector= relativesector;
1565         t.Rotate(rotation);     
1566       }
1567       t.PropagateTo(x);
1568       //
1569       t.fCurrentCluster = cl; 
1570       t.fRow = nr;
1571       Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1572       if ((tpcindex&0x8000)==0) accept =0;
1573       if (accept<3) { 
1574         //if founded cluster is acceptible
1575         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
1576           t.fErrorY2 += 0.03;
1577           t.fErrorZ2 += 0.03; 
1578           t.fErrorY2 *= 3;
1579           t.fErrorZ2 *= 3; 
1580         }
1581         t.fNFoundable++;
1582         UpdateTrack(&t,accept);
1583         return 1;
1584       }    
1585     }
1586   }
1587   if (fIteration>1) return 0;  // not look for new cluster during refitting
1588   //
1589   UInt_t index=0;
1590   if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;
1591   Double_t  y=t.GetYat(x);
1592   if (TMath::Abs(y)>ymax){
1593     if (y > ymax) {
1594       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1595       if (!t.Rotate(fSectors->GetAlpha())) 
1596         return 0;
1597     } else if (y <-ymax) {
1598       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1599       if (!t.Rotate(-fSectors->GetAlpha())) 
1600         return 0;
1601     }
1602     //return 1;
1603   }
1604   //
1605   if (!t.PropagateTo(x)) {
1606     if (fIteration==0) t.fRemoval = 10;
1607     return 0;
1608   }
1609   y=t.GetY(); 
1610   Double_t z=t.GetZ();
1611   //
1612   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1613   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1614   Double_t  roady  =1.;
1615   Double_t  roadz = 1.;
1616   //
1617   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1618     t.fInDead = kTRUE;
1619     t.SetClusterIndex2(nr,-1); 
1620     return 0;
1621   } 
1622   else
1623     {
1624       if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10)) t.fNFoundable++;
1625       else
1626         return 0;
1627     }   
1628   //calculate 
1629   if (krow) {
1630     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);    
1631     cl = krow.FindNearest2(y,z,roady,roadz,index);    
1632     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);       
1633   }  
1634   if (cl) {
1635     t.fCurrentCluster = cl; 
1636     t.fRow = nr;
1637     if (fIteration==2&&cl->IsUsed(10)) return 0; 
1638     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1639     if (fIteration==2&&cl->IsUsed(11)) {
1640       t.fErrorY2 += 0.03;
1641       t.fErrorZ2 += 0.03; 
1642       t.fErrorY2 *= 3;
1643       t.fErrorZ2 *= 3; 
1644     }
1645     /*    
1646     if (t.fCurrentCluster->IsUsed(10)){
1647       //
1648       //     
1649
1650       t.fNShared++;
1651       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1652         t.fRemoval =10;
1653         return 0;
1654       }
1655     }
1656     */
1657     if (accept<3) UpdateTrack(&t,accept);
1658
1659   } else {  
1660     if ( fIteration==0 && t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1661     
1662   }
1663   return 1;
1664 }
1665
1666 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1667   //-----------------------------------------------------------------
1668   // This function tries to find a track prolongation to next pad row
1669   //-----------------------------------------------------------------
1670   //
1671   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1672   Double_t y,z; 
1673   if (!t.GetProlongation(x,y,z)) {
1674     t.fRemoval = 10;
1675     return 0;
1676   }
1677   //
1678   //
1679   if (TMath::Abs(y)>ymax){
1680     
1681     if (y > ymax) {
1682       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1683       if (!t.Rotate(fSectors->GetAlpha())) 
1684         return 0;
1685     } else if (y <-ymax) {
1686       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1687       if (!t.Rotate(-fSectors->GetAlpha())) 
1688         return 0;
1689     }
1690     if (!t.PropagateTo(x)) {
1691       return 0;
1692     } 
1693     t.GetProlongation(x,y,z);
1694   }
1695   //
1696   // update current shape info every 3 pad-row
1697   if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1698     //    t.fCurrentSigmaY = GetSigmaY(&t);
1699     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1700     GetShape(&t,nr);
1701   }
1702   //  
1703   AliTPCclusterMI *cl=0;
1704   UInt_t index=0;
1705   
1706   
1707   //Int_t nr2 = nr;
1708   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1709   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1710   Double_t  roady  =1.;
1711   Double_t  roadz = 1.;
1712   //
1713   Int_t row = nr;
1714   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1715     t.fInDead = kTRUE;
1716     t.SetClusterIndex2(row,-1); 
1717     return 0;
1718   } 
1719   else
1720     {
1721       if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1722     }   
1723   //calculate 
1724   
1725   if ((cl==0)&&(krow)) {
1726     //    cl = krow.FindNearest2(y+10,z,roady,roadz,index);    
1727     cl = krow.FindNearest2(y,z,roady,roadz,index);    
1728
1729     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);       
1730   }  
1731
1732   if (cl) {
1733     t.fCurrentCluster = cl; 
1734     //    Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);        
1735     //if (accept<3){
1736       t.SetClusterIndex2(row,index);
1737       t.fClusterPointer[row] = cl;
1738       //}
1739   }
1740   return 1;
1741 }
1742
1743
1744
1745 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
1746   //-----------------------------------------------------------------
1747   // This function tries to find a track prolongation to next pad row
1748   //-----------------------------------------------------------------
1749   t.fCurrentCluster  = 0;
1750   t.fCurrentClusterIndex1 = 0;   
1751    
1752   Double_t xt=t.GetX();
1753   Int_t     row = GetRowNumber(xt)-1; 
1754   Double_t  ymax= GetMaxY(nr);
1755
1756   if (row < nr) return 1; // don't prolongate if not information until now -
1757   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1758     t.fRemoval =10;
1759     return 0;  // not prolongate strongly inclined tracks
1760   } 
1761   if (TMath::Abs(t.GetSnp())>0.95) {
1762     t.fRemoval =10;
1763     return 0;  // not prolongate strongly inclined tracks
1764   }
1765
1766   Double_t x= GetXrow(nr);
1767   Double_t y,z;
1768   //t.PropagateTo(x+0.02);
1769   //t.PropagateTo(x+0.01);
1770   if (!t.PropagateTo(x)){
1771     return 0;
1772   }
1773   //
1774   y=t.GetY();
1775   z=t.GetZ();
1776
1777   if (TMath::Abs(y)>ymax){
1778     if (y > ymax) {
1779       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1780       if (!t.Rotate(fSectors->GetAlpha())) 
1781         return 0;
1782     } else if (y <-ymax) {
1783       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1784       if (!t.Rotate(-fSectors->GetAlpha())) 
1785         return 0;
1786     }
1787     //    if (!t.PropagateTo(x)){
1788     //  return 0;
1789     //}
1790     return 1;
1791     //y = t.GetY();    
1792   }
1793   //
1794
1795   AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1796
1797   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1798     t.fInDead = kTRUE;
1799     t.SetClusterIndex2(nr,-1); 
1800     return 0;
1801   } 
1802   else
1803     {
1804       if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10)) t.fNFoundable++;
1805       else
1806         return 0;      
1807     }
1808
1809   // update current
1810   if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1811     //    t.fCurrentSigmaY = GetSigmaY(&t);
1812     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1813     GetShape(&t,nr);
1814   }
1815     
1816   AliTPCclusterMI *cl=0;
1817   UInt_t index=0;
1818   //
1819   Double_t roady = 1.;
1820   Double_t roadz = 1.;
1821   //
1822
1823   if (!cl){
1824     index = t.GetClusterIndex2(nr);    
1825     if ( (index>0) && (index&0x8000)==0){
1826       cl = t.fClusterPointer[nr];
1827       if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1828       t.fCurrentClusterIndex1 = index;
1829       if (cl) {
1830         t.fCurrentCluster  = cl;
1831         return 1;
1832       }
1833     }
1834   }
1835
1836   if (krow) {    
1837     //cl = krow.FindNearest2(y+10,z,roady,roadz,index);      
1838     cl = krow.FindNearest2(y,z,roady,roadz,index);      
1839   }
1840
1841   if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);   
1842   t.fCurrentCluster  = cl;
1843
1844   return 1;
1845 }
1846
1847
1848 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1849   //-----------------------------------------------------------------
1850   // This function tries to find a track prolongation to next pad row
1851   //-----------------------------------------------------------------
1852
1853   //update error according neighborhoud
1854
1855   if (t.fCurrentCluster) {
1856     t.fRow = nr; 
1857     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1858     
1859     if (t.fCurrentCluster->IsUsed(10)){
1860       //
1861       //
1862       //  t.fErrorZ2*=2;
1863       //  t.fErrorY2*=2;
1864       t.fNShared++;
1865       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1866         t.fRemoval =10;
1867         return 0;
1868       }
1869     }   
1870     if (fIteration>0) accept = 0;
1871     if (accept<3)  UpdateTrack(&t,accept);  
1872  
1873   } else {
1874     if (fIteration==0){
1875       if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;      
1876       if (  t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;      
1877
1878       if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1879     }
1880   }
1881   return 1;
1882 }
1883
1884
1885
1886 //_____________________________________________________________________________
1887 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1888   //-----------------------------------------------------------------
1889   // This function tries to find a track prolongation.
1890   //-----------------------------------------------------------------
1891   Double_t xt=t.GetX();
1892   //
1893   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1894   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
1895   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
1896   //
1897   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1898     
1899   Int_t first = GetRowNumber(xt)-1;
1900   for (Int_t nr= first; nr>=rf; nr-=step) {
1901     // update kink info
1902     if (t.GetKinkIndexes()[0]>0){
1903       for (Int_t i=0;i<3;i++){
1904         Int_t index = t.GetKinkIndexes()[i];
1905         if (index==0) break;
1906         if (index<0) continue;
1907         //
1908         AliESDkink * kink = fEvent->GetKink(index-1);
1909         if (!kink){
1910           printf("PROBLEM\n");
1911         }
1912         else{
1913           Int_t kinkrow = kink->fRow0+Int_t(1./(0.1+kink->fAngle[2]));
1914           if (kinkrow==nr){
1915             AliExternalTrackParam paramd(t);
1916             kink->SetDaughter(paramd);
1917             kink->Update();
1918             kink->fStatus+=100;
1919           }
1920         }
1921       }
1922     }
1923
1924     if (nr==80) t.UpdateReference();
1925     if (nr<fInnerSec->GetNRows()) 
1926       fSectors = fInnerSec;
1927     else
1928       fSectors = fOuterSec;
1929     if (FollowToNext(t,nr)==0) 
1930       if (!t.IsActive()) 
1931         return 0;
1932     
1933   }   
1934   return 1;
1935 }
1936
1937
1938 //_____________________________________________________________________________
1939 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1940   //-----------------------------------------------------------------
1941   // This function tries to find a track prolongation.
1942   //-----------------------------------------------------------------
1943   Double_t xt=t.GetX();
1944   //
1945   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1946   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
1947   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
1948   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1949     
1950   for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1951     
1952     if (FollowToNextFast(t,nr)==0) 
1953       if (!t.IsActive()) return 0;
1954     
1955   }   
1956   return 1;
1957 }
1958
1959
1960
1961
1962
1963 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1964   //-----------------------------------------------------------------
1965   // This function tries to find a track prolongation.
1966   //-----------------------------------------------------------------
1967   //  Double_t xt=t.GetX();  
1968   //
1969   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1970   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
1971   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
1972   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1973     
1974   Int_t first = t.fFirstPoint;
1975   //
1976   if (first<0) first=0;
1977   for (Int_t nr=first; nr<=rf; nr++) {
1978     //if ( (t.GetSnp()<0.9))
1979     if (t.GetKinkIndexes()[0]<0){
1980       for (Int_t i=0;i<3;i++){
1981         Int_t index = t.GetKinkIndexes()[i];
1982         if (index==0) break;
1983         if (index>0) continue;
1984         index = TMath::Abs(index);
1985         AliESDkink * kink = fEvent->GetKink(index-1);
1986         if (!kink){
1987           printf("PROBLEM\n");
1988         }
1989         else{
1990           Int_t kinkrow = kink->fRow0-Int_t(1./(0.1+kink->fAngle[2]));
1991           if (kinkrow==nr){
1992             AliExternalTrackParam paramm(t);
1993             kink->SetMother(paramm);
1994             kink->Update();
1995             kink->fStatus+=10;
1996           }
1997         }
1998       }
1999     }
2000     if (nr<fInnerSec->GetNRows()) 
2001       fSectors = fInnerSec;
2002     else
2003       fSectors = fOuterSec;
2004     FollowToNext(t,nr);                                                             
2005   }   
2006   return 1;
2007 }
2008
2009
2010
2011
2012    
2013 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2014 {
2015   //
2016   //
2017   sum1=0;
2018   sum2=0;
2019   Int_t sum=0;
2020   //
2021   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2022   dz2*=dz2;  
2023
2024   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2025   dy2*=dy2;
2026   Float_t distance = TMath::Sqrt(dz2+dy2);
2027   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2028  
2029   //  Int_t offset =0;
2030   Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
2031   Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
2032   if (lastpoint>160) 
2033     lastpoint =160;
2034   if (firstpoint<0) 
2035     firstpoint = 0;
2036   if (firstpoint>lastpoint) {
2037     firstpoint =lastpoint;
2038     //    lastpoint  =160;
2039   }
2040     
2041   
2042   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2043     if (s1->GetClusterIndex2(i)>0) sum1++;
2044     if (s2->GetClusterIndex2(i)>0) sum2++;
2045     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2046       sum++;
2047     }
2048   }
2049   if (sum<5) return 0;
2050
2051   Float_t summin = TMath::Min(sum1+1,sum2+1);
2052   Float_t ratio = (sum+1)/Float_t(summin);
2053   return ratio;
2054 }
2055
2056 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2057 {
2058   //
2059   //
2060   if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
2061   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
2062
2063   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2064   dz2*=dz2;
2065   Float_t dy2 =(s1->GetY() - s2->GetY());
2066   dy2*=dy2;
2067   Float_t distance = dz2+dy2;
2068   if (distance>325.) return ; // if there are far away  - not overlap - to reduce combinatorics
2069   
2070   //
2071   Int_t sumshared=0;
2072   //
2073   Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
2074   Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
2075   //
2076   if (firstpoint>=lastpoint-5) return;;
2077
2078   for (Int_t i=firstpoint;i<lastpoint;i++){
2079     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2080     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2081       sumshared++;
2082     }
2083   }
2084   if (sumshared>4){
2085     // sign clusters
2086     //
2087     for (Int_t i=firstpoint;i<lastpoint;i++){
2088       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2089       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2090         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2091         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);; 
2092         if (s1->IsActive()&&s2->IsActive()){
2093           p1->fIsShared = kTRUE;
2094           p2->fIsShared = kTRUE;
2095         }       
2096       }
2097     }
2098   }
2099   //  
2100   if (sumshared>10){
2101     for (Int_t i=0;i<4;i++){
2102       if (s1->fOverlapLabels[3*i]==0){
2103         s1->fOverlapLabels[3*i] = s2->GetLabel();
2104         s1->fOverlapLabels[3*i+1] = sumshared;
2105         s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
2106         break;
2107       } 
2108     }
2109     for (Int_t i=0;i<4;i++){
2110       if (s2->fOverlapLabels[3*i]==0){
2111         s2->fOverlapLabels[3*i] = s1->GetLabel();
2112         s2->fOverlapLabels[3*i+1] = sumshared;
2113         s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
2114         break;
2115       } 
2116     }    
2117   }
2118   
2119 }
2120
2121 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
2122 {
2123   //
2124   //sort trackss according sectors
2125   //  
2126   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2127     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2128     if (!pt) continue;
2129     //if (pt) RotateToLocal(pt);
2130     pt->fSort = 0;
2131   }
2132   arr->UnSort();
2133   arr->Sort();  // sorting according z
2134   arr->Expand(arr->GetEntries());
2135   //
2136   //
2137   Int_t nseed=arr->GetEntriesFast();
2138   for (Int_t i=0; i<nseed; i++) {
2139     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2140     if (!pt) continue;
2141     for (Int_t j=0;j<=12;j++){
2142       pt->fOverlapLabels[j] =0;
2143     }
2144   }
2145   for (Int_t i=0; i<nseed; i++) {
2146     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2147     if (!pt) continue;
2148     if (pt->fRemoval>10) continue;
2149     for (Int_t j=i+1; j<nseed; j++){
2150       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2151       //      if (pt2){
2152       if (pt2->fRemoval<=10) {
2153         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2154         SignShared(pt,pt2);
2155       }
2156     }  
2157   }
2158 }
2159
2160 void  AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2,  Int_t removalindex)
2161 {
2162   //
2163   //sort trackss according sectors
2164   //
2165   if (fDebug&1) {
2166     Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries());
2167   }
2168   //
2169   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2170     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2171     if (!pt) continue;
2172     pt->fSort = 0;
2173   }
2174   arr->UnSort();
2175   arr->Sort();  // sorting according z
2176   arr->Expand(arr->GetEntries());
2177   //
2178   //reset overlap labels
2179   //
2180   Int_t nseed=arr->GetEntriesFast();
2181   for (Int_t i=0; i<nseed; i++) {
2182     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2183     if (!pt) continue;
2184     pt->SetUniqueID(i);
2185     for (Int_t j=0;j<=12;j++){
2186       pt->fOverlapLabels[j] =0;
2187     }
2188   }
2189   //
2190   //sign shared tracks
2191   for (Int_t i=0; i<nseed; i++) {
2192     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2193     if (!pt) continue;
2194     if (pt->fRemoval>10) continue;
2195     Float_t deltac = pt->GetC()*0.1;
2196     for (Int_t j=i+1; j<nseed; j++){
2197       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2198       //      if (pt2){
2199       if (pt2->fRemoval<=10) {
2200         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2201         if (TMath::Abs(pt->GetC()  -pt2->GetC())>deltac) continue;
2202         if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2203         //
2204         SignShared(pt,pt2);
2205       }
2206     }
2207   }
2208   //
2209   // remove highly shared tracks
2210   for (Int_t i=0; i<nseed; i++) {
2211     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2212     if (!pt) continue;
2213     if (pt->fRemoval>10) continue;
2214     //
2215     Int_t sumshared =0;
2216     for (Int_t j=0;j<4;j++){
2217       sumshared = pt->fOverlapLabels[j*3+1];      
2218     }
2219     Float_t factor = factor1;
2220     if (pt->fRemoval>0) factor = factor2;
2221     if (sumshared/pt->GetNumberOfClusters()>factor){
2222       for (Int_t j=0;j<4;j++){
2223         if (pt->fOverlapLabels[3*j]==0) continue;
2224         if (pt->fOverlapLabels[3*j+1]<5) continue; 
2225         if (pt->fRemoval==removalindex) continue;      
2226         AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2227         if (!pt2) continue;
2228         if (pt2->GetSigma2C()<pt->GetSigma2C()){
2229           //      pt->fRemoval = removalindex;
2230           delete arr->RemoveAt(i);        
2231           break;
2232         }
2233       }      
2234     }
2235   }
2236   arr->Compress();
2237   if (fDebug&1) {
2238     Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries());
2239   }
2240 }
2241
2242
2243
2244
2245
2246
2247 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2248 {
2249   //
2250   //sort tracks in array according mode criteria
2251   Int_t nseed = arr->GetEntriesFast();    
2252   for (Int_t i=0; i<nseed; i++) {
2253     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2254     if (!pt) {
2255       continue;
2256     }
2257     pt->fSort = mode;
2258   }
2259   arr->UnSort();
2260   arr->Sort();
2261 }
2262
2263 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t removalindex)
2264 {
2265
2266   //Loop over all tracks and remove "overlaps"
2267   //
2268   //
2269   Int_t nseed = arr->GetEntriesFast();  
2270   Int_t good =0;
2271
2272   for (Int_t i=0; i<nseed; i++) {
2273     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2274     if (!pt) {
2275       delete arr->RemoveAt(i);
2276     }
2277     else{
2278       pt->fSort =1;
2279       pt->fBSigned = kFALSE;
2280     }
2281   }
2282   arr->Compress();
2283   nseed = arr->GetEntriesFast();
2284   arr->UnSort();
2285   arr->Sort();
2286   //
2287   //unsign used
2288   UnsignClusters();
2289   //
2290   for (Int_t i=0; i<nseed; i++) {
2291     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2292     if (!pt) {
2293       continue;
2294     }    
2295     Int_t found,foundable,shared;
2296     if (pt->IsActive()) 
2297       pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2298     else
2299       pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE); 
2300     //
2301     Double_t factor = factor2;
2302     if (pt->fBConstrain) factor = factor1;
2303
2304     if ((Float_t(shared)/Float_t(found))>factor){
2305       pt->Desactivate(removalindex);
2306       continue;
2307     }
2308
2309     good++;
2310     for (Int_t i=0; i<160; i++) {
2311       Int_t index=pt->GetClusterIndex2(i);
2312       if (index<0 || index&0x8000 ) continue;
2313       AliTPCclusterMI *c= pt->fClusterPointer[i];        
2314       if (!c) continue;
2315       //      if (!c->IsUsed(10)) c->Use(10);
2316       //if (pt->IsActive()) 
2317       c->Use(10);  
2318       //else
2319       //        c->Use(5);
2320     }
2321     
2322   }
2323   fNtracks = good;
2324   if (fDebug>0){
2325     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2326   }
2327 }
2328
2329
2330 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2331 {
2332
2333   //Loop over all tracks and remove "overlaps"
2334   //
2335   //
2336   UnsignClusters();
2337   //
2338   Int_t nseed = arr->GetEntriesFast();  
2339   Float_t * quality = new Float_t[nseed];
2340   Int_t   * indexes = new Int_t[nseed];
2341   Int_t good =0;
2342   //
2343   //
2344   for (Int_t i=0; i<nseed; i++) {
2345     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2346     if (!pt){
2347       quality[i]=-1;
2348       continue;
2349     }
2350     pt->UpdatePoints();    //select first last max dens points
2351     Float_t * points = pt->GetPoints();
2352     if (points[3]<0.8) quality[i] =-1;
2353     //
2354     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2355   }
2356   TMath::Sort(nseed,quality,indexes);
2357   //
2358   //
2359   for (Int_t itrack=0; itrack<nseed; itrack++) {
2360     Int_t trackindex = indexes[itrack];
2361     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);    
2362     if (quality[trackindex]<0){
2363       if (pt) {
2364         delete arr->RemoveAt(trackindex);
2365       }
2366       else{
2367         arr->RemoveAt(trackindex);
2368       }
2369       continue;
2370     }
2371     //
2372     Int_t first = Int_t(pt->GetPoints()[0]);
2373     Int_t last  = Int_t(pt->GetPoints()[2]);
2374     Double_t factor = (pt->fBConstrain) ? factor1: factor2;
2375     //
2376     Int_t found,foundable,shared;
2377     pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
2378     Float_t sharedfactor = Float_t(shared)/Float_t(found);
2379     //
2380     if (Float_t(shared)/Float_t(found)>factor){
2381       if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2382       delete arr->RemoveAt(trackindex);
2383       continue;
2384     }
2385
2386     if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
2387       if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2388       delete arr->RemoveAt(trackindex);
2389       continue;
2390     }
2391
2392     good++;
2393     if (sharedfactor>0.4) continue;
2394     for (Int_t i=first; i<last; i++) {
2395       Int_t index=pt->GetClusterIndex2(i);
2396       // if (index<0 || index&0x8000 ) continue;
2397       if (index<0 || index&0x8000 ) continue;
2398       AliTPCclusterMI *c= pt->fClusterPointer[i];        
2399       if (!c) continue;
2400       c->Use(10);  
2401     }    
2402   }
2403   fNtracks = good;
2404   if (fDebug>0){
2405     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2406   }
2407   delete []quality;
2408   delete []indexes;
2409 }
2410
2411 void AliTPCtrackerMI::UnsignClusters() 
2412 {
2413   //
2414   // loop over all clusters and unsign them
2415   //
2416   
2417   for (Int_t sec=0;sec<fkNIS;sec++){
2418     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2419       AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2420       for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2421         //      if (cl[icl].IsUsed(10))         
2422         cl[icl].Use(-1);
2423       cl = fInnerSec[sec][row].fClusters2;
2424       for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2425         //if (cl[icl].IsUsed(10))       
2426           cl[icl].Use(-1);      
2427     }
2428   }
2429   
2430   for (Int_t sec=0;sec<fkNOS;sec++){
2431     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2432       AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2433       for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2434         //if (cl[icl].IsUsed(10))       
2435           cl[icl].Use(-1);
2436       cl = fOuterSec[sec][row].fClusters2;
2437       for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2438         //if (cl[icl].IsUsed(10))       
2439         cl[icl].Use(-1);      
2440     }
2441   }
2442   
2443 }
2444
2445
2446
2447 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2448 {
2449   //
2450   //sign clusters to be "used"
2451   //
2452   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2453   // loop over "primaries"
2454   
2455   Float_t sumdens=0;
2456   Float_t sumdens2=0;
2457   Float_t sumn   =0;
2458   Float_t sumn2  =0;
2459   Float_t sumchi =0;
2460   Float_t sumchi2 =0;
2461
2462   Float_t sum    =0;
2463
2464   TStopwatch timer;
2465   timer.Start();
2466
2467   Int_t nseed = arr->GetEntriesFast();
2468   for (Int_t i=0; i<nseed; i++) {
2469     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2470     if (!pt) {
2471       continue;
2472     }    
2473     if (!(pt->IsActive())) continue;
2474     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2475     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2476       sumdens += dens;
2477       sumdens2+= dens*dens;
2478       sumn    += pt->GetNumberOfClusters();
2479       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2480       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2481       if (chi2>5) chi2=5;
2482       sumchi  +=chi2;
2483       sumchi2 +=chi2*chi2;
2484       sum++;
2485     }
2486   }
2487
2488   Float_t mdensity = 0.9;
2489   Float_t meann    = 130;
2490   Float_t meanchi  = 1;
2491   Float_t sdensity = 0.1;
2492   Float_t smeann    = 10;
2493   Float_t smeanchi  =0.4;
2494   
2495
2496   if (sum>20){
2497     mdensity = sumdens/sum;
2498     meann    = sumn/sum;
2499     meanchi  = sumchi/sum;
2500     //
2501     sdensity = sumdens2/sum-mdensity*mdensity;
2502     sdensity = TMath::Sqrt(sdensity);
2503     //
2504     smeann   = sumn2/sum-meann*meann;
2505     smeann   = TMath::Sqrt(smeann);
2506     //
2507     smeanchi = sumchi2/sum - meanchi*meanchi;
2508     smeanchi = TMath::Sqrt(smeanchi);
2509   }
2510
2511
2512   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2513   //
2514   for (Int_t i=0; i<nseed; i++) {
2515     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2516     if (!pt) {
2517       continue;
2518     }
2519     if (pt->fBSigned) continue;
2520     if (pt->fBConstrain) continue;    
2521     //if (!(pt->IsActive())) continue;
2522     /*
2523     Int_t found,foundable,shared;    
2524     pt->GetClusterStatistic(0,160,found, foundable,shared);
2525     if (shared/float(found)>0.3) {
2526       if (shared/float(found)>0.9 ){
2527         //delete arr->RemoveAt(i);
2528       }
2529       continue;
2530     }
2531     */
2532     Bool_t isok =kFALSE;
2533     if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2534       isok = kTRUE;
2535     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2536       isok =kTRUE;
2537     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2538       isok =kTRUE;
2539     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2540       isok =kTRUE;
2541     
2542     if (isok)     
2543       for (Int_t i=0; i<160; i++) {     
2544         Int_t index=pt->GetClusterIndex2(i);
2545         if (index<0) continue;
2546         AliTPCclusterMI *c= pt->fClusterPointer[i];
2547         if (!c) continue;
2548         //if (!(c->IsUsed(10))) c->Use();  
2549         c->Use(10);  
2550       }
2551   }
2552   
2553   
2554   //
2555   Double_t maxchi  = meanchi+2.*smeanchi;
2556
2557   for (Int_t i=0; i<nseed; i++) {
2558     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2559     if (!pt) {
2560       continue;
2561     }    
2562     //if (!(pt->IsActive())) continue;
2563     if (pt->fBSigned) continue;
2564     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2565     if (chi>maxchi) continue;
2566
2567     Float_t bfactor=1;
2568     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2569    
2570     //sign only tracks with enoug big density at the beginning
2571     
2572     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue; 
2573     
2574     
2575     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2576     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2577    
2578     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2579     if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2580       minn=0;
2581
2582     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2583       //Int_t noc=pt->GetNumberOfClusters();
2584       pt->fBSigned = kTRUE;
2585       for (Int_t i=0; i<160; i++) {
2586
2587         Int_t index=pt->GetClusterIndex2(i);
2588         if (index<0) continue;
2589         AliTPCclusterMI *c= pt->fClusterPointer[i];
2590         if (!c) continue;
2591         //      if (!(c->IsUsed(10))) c->Use();  
2592         c->Use(10);  
2593       }
2594     }
2595   }
2596   //  gLastCheck = nseed;
2597   //  arr->Compress();
2598   if (fDebug>0){
2599     timer.Print();
2600   }
2601 }
2602
2603
2604 void  AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2605 {
2606   // stop not active tracks
2607   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2608   // take th2 as threshold for number of founded to number of foundable on last 20 active rows 
2609   Int_t nseed = arr->GetEntriesFast();  
2610   //
2611   for (Int_t i=0; i<nseed; i++) {
2612     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2613     if (!pt) {
2614       continue;
2615     }
2616     if (!(pt->IsActive())) continue;
2617     StopNotActive(pt,row0,th0, th1,th2);
2618   }
2619 }
2620
2621
2622
2623 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2624  Float_t th2) const
2625 {
2626   // stop not active tracks
2627   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2628   // take th2 as threshold for number of founded to number of foundable on last 20 active rows 
2629   Int_t sumgood1  = 0;
2630   Int_t sumgood2  = 0;
2631   Int_t foundable = 0;
2632   Int_t maxindex = seed->fLastPoint;  //last foundable row
2633   if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2634     seed->Desactivate(10) ;
2635     return;
2636   }
2637
2638   for (Int_t i=row0; i<maxindex; i++){
2639     Int_t index = seed->GetClusterIndex2(i);
2640     if (index!=-1) foundable++;
2641     //if (!c) continue;
2642     if (foundable<=30) sumgood1++;
2643     if (foundable<=50) {
2644       sumgood2++;
2645     }
2646     else{ 
2647       break;
2648     }        
2649   }
2650   if (foundable>=30.){ 
2651      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2652   }
2653   if (foundable>=50)
2654     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2655 }
2656
2657
2658 Int_t AliTPCtrackerMI::RefitInward(AliESD *event)
2659 {
2660   //
2661   // back propagation of ESD tracks
2662   //
2663   //return 0;
2664   fEvent = event;
2665   ReadSeeds(event,2);
2666   fIteration=2;
2667   //PrepareForProlongation(fSeeds,1);
2668   PropagateForward2(fSeeds);
2669   Int_t ntracks=0;
2670   Int_t nseed = fSeeds->GetEntriesFast();
2671   for (Int_t i=0;i<nseed;i++){
2672     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2673     if (!seed) continue;
2674     seed->PropagateTo(fParam->GetInnerRadiusLow());
2675     seed->UpdatePoints();
2676     AliESDtrack *esd=event->GetTrack(i);
2677     seed->CookdEdx(0.02,0.6);
2678     CookLabel(seed,0.1); //For comparison only
2679     if (seed->GetNumberOfClusters()>15){
2680       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit); 
2681       esd->SetTPCPoints(seed->GetPoints());
2682       ntracks++;
2683     }
2684     else{
2685       //printf("problem\n");
2686     }
2687   }
2688   //FindKinks(fSeeds,event);
2689   Info("RefitInward","Number of refitted tracks %d",ntracks);
2690   fEvent =0;
2691   //WriteTracks();
2692   return 0;
2693 }
2694
2695
2696 Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2697 {
2698   //
2699   // back propagation of ESD tracks
2700   //
2701
2702   fEvent = event;
2703   fIteration = 1;
2704   ReadSeeds(event,0);
2705   PropagateBack(fSeeds);
2706   Int_t nseed = fSeeds->GetEntriesFast();
2707   Int_t ntracks=0;
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)  UpdateKinkQualityM(seed);  // update quality informations for kinks
2712     seed->UpdatePoints();
2713     AliESDtrack *esd=event->GetTrack(i);
2714     seed->CookdEdx(0.02,0.6);
2715     CookLabel(seed,0.1); //For comparison only
2716     if (seed->GetNumberOfClusters()>15){
2717       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2718       esd->SetTPCPoints(seed->GetPoints());
2719       ntracks++;
2720     }
2721   }
2722   //FindKinks(fSeeds,event);
2723   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2724   fEvent =0;
2725   //WriteTracks();
2726   return 0;
2727 }
2728
2729
2730 void AliTPCtrackerMI::DeleteSeeds()
2731 {
2732   //
2733   //delete Seeds
2734   Int_t nseed = fSeeds->GetEntriesFast();
2735   for (Int_t i=0;i<nseed;i++){
2736     AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2737     if (seed) delete fSeeds->RemoveAt(i);
2738   }
2739   delete fSeeds;
2740   fSeeds =0;
2741 }
2742
2743 void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
2744 {
2745   //
2746   //read seeds from the event
2747   
2748   Int_t nentr=event->GetNumberOfTracks();
2749   if (fDebug>0){
2750     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2751   }
2752   if (fSeeds) 
2753     DeleteSeeds();
2754   if (!fSeeds){   
2755     fSeeds = new TObjArray(nentr);
2756   }
2757   UnsignClusters();
2758   //  Int_t ntrk=0;  
2759   for (Int_t i=0; i<nentr; i++) {
2760     AliESDtrack *esd=event->GetTrack(i);
2761     ULong_t status=esd->GetStatus();
2762     if (!(status&AliESDtrack::kTPCin)) continue;
2763     AliTPCtrack t(*esd);
2764     AliTPCseed *seed=new AliTPCseed(t/*,t.GetAlpha()*/);seed->ResetClusters();
2765     for (Int_t ikink=0;ikink<3;ikink++) seed->GetKinkIndexes()[ikink] = esd->GetKinkIndex(ikink);
2766     if ((status==AliESDtrack::kTPCin)&&(direction==1)) seed->ResetCovariance(); 
2767     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
2768     if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2769       fSeeds->AddAt(0,i);
2770       delete seed;
2771       continue;    
2772     }
2773     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 )  {
2774       Double_t par0[5],par1[5],x;
2775       esd->GetInnerExternalParameters(x,par0);
2776       esd->GetExternalParameters(x,par1);
2777       Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2778       Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2779       Double_t trdchi2=0;
2780       if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2781       //reset covariance if suspicious 
2782       if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2783         seed->ResetCovariance();
2784     }
2785
2786     //
2787     //
2788     // rotate to the local coordinate system
2789    
2790     fSectors=fInnerSec; fN=fkNIS;
2791     
2792     Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2793     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2794     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2795     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2796     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2797     if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2798     if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2799     alpha-=seed->GetAlpha();  
2800     if (!seed->Rotate(alpha)) {
2801       delete seed;
2802       continue;
2803     }
2804     seed->fEsd = esd;
2805     //
2806     //seed->PropagateTo(fSectors->GetX(0));
2807     //
2808     //    Int_t index = esd->GetTPCindex();
2809     //AliTPCseed * seed2= (AliTPCseed*)fSeeds->At(index);
2810     //if (direction==2){
2811     //  AliTPCseed * seed2  = ReSeed(seed,0.,0.5,1.);
2812     //  if (seed2) {
2813     //  delete seed;
2814     //  seed = seed2;
2815     //  }
2816     //}
2817     //
2818     // sign clusters
2819     for (Int_t irow=0;irow<160;irow++){
2820       Int_t index = seed->GetClusterIndex2(irow);    
2821       if (index>0){ 
2822         //
2823         AliTPCclusterMI * cl = GetClusterMI(index);
2824         seed->fClusterPointer[irow] = cl;
2825         if (cl){
2826           if ((index & 0x8000)==0){
2827             cl->Use(10);  // accepted cluster     
2828           }else{
2829             cl->Use(6);   // close cluster not accepted
2830           }     
2831         }else{
2832            Info("ReadSeeds","Not found cluster");
2833         }
2834       }
2835     }
2836     fSeeds->AddAt(seed,i);
2837   }
2838 }
2839
2840
2841
2842 //_____________________________________________________________________________
2843 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
2844                                  Float_t deltay, Int_t ddsec) {
2845   //-----------------------------------------------------------------
2846   // This function creates track seeds.
2847   // SEEDING WITH VERTEX CONSTRAIN 
2848   //-----------------------------------------------------------------
2849   // cuts[0]   - fP4 cut
2850   // cuts[1]   - tan(phi)  cut
2851   // cuts[2]   - zvertex cut
2852   // cuts[3]   - fP3 cut
2853   Int_t nin0  = 0;
2854   Int_t nin1  = 0;
2855   Int_t nin2  = 0;
2856   Int_t nin   = 0;
2857   Int_t nout1 = 0;
2858   Int_t nout2 = 0;
2859
2860   Double_t x[5], c[15];
2861   //  Int_t di = i1-i2;
2862   //
2863   AliTPCseed * seed = new AliTPCseed;
2864   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2865   Double_t cs=cos(alpha), sn=sin(alpha);
2866   //
2867   //  Double_t x1 =fOuterSec->GetX(i1);
2868   //Double_t xx2=fOuterSec->GetX(i2);
2869   
2870   Double_t x1 =GetXrow(i1);
2871   Double_t xx2=GetXrow(i2);
2872
2873   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2874
2875   Int_t imiddle = (i2+i1)/2;    //middle pad row index
2876   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2877   const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2878   //
2879   Int_t ns =sec;   
2880
2881   const AliTPCRow& kr1=GetRow(ns,i1);
2882   Double_t ymax  = GetMaxY(i1)-kr1.fDeadZone-1.5;  
2883   Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;  
2884
2885   //
2886   // change cut on curvature if it can't reach this layer
2887   // maximal curvature set to reach it
2888   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2889   if (dvertexmax*0.5*cuts[0]>0.85){
2890     cuts[0] = 0.85/(dvertexmax*0.5+1.);
2891   }
2892   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
2893
2894   //  Int_t ddsec = 1;
2895   if (deltay>0) ddsec = 0; 
2896   // loop over clusters  
2897   for (Int_t is=0; is < kr1; is++) {
2898     //
2899     if (kr1[is]->IsUsed(10)) continue;
2900     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
2901     //if (TMath::Abs(y1)>ymax) continue;
2902
2903     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2904
2905     // find possible directions    
2906     Float_t anglez = (z1-z3)/(x1-x3); 
2907     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z      
2908     //
2909     //
2910     //find   rotation angles relative to line given by vertex and point 1
2911     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2912     Double_t dvertex  = TMath::Sqrt(dvertex2);
2913     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
2914     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);            
2915     
2916     //
2917     // loop over 2 sectors
2918     Int_t dsec1=-ddsec;
2919     Int_t dsec2= ddsec;
2920     if (y1<0)  dsec2= 0;
2921     if (y1>0)  dsec1= 0;
2922     
2923     Double_t dddz1=0;  // direction of delta inclination in z axis
2924     Double_t dddz2=0;
2925     if ( (z1-z3)>0)
2926       dddz1 =1;    
2927     else
2928       dddz2 =1;
2929     //
2930     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2931       Int_t sec2 = sec + dsec;
2932       // 
2933       //      AliTPCRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2934       //AliTPCRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2935       AliTPCRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
2936       AliTPCRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2937       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2938       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2939
2940       // rotation angles to p1-p3
2941       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;            
2942       Double_t x2,   y2,   z2; 
2943       //
2944       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2945
2946       //
2947       Double_t dxx0 =  (xx2-x3)*cs13r;
2948       Double_t dyy0 =  (xx2-x3)*sn13r;
2949       for (Int_t js=index1; js < index2; js++) {
2950         const AliTPCclusterMI *kcl = kr2[js];
2951         if (kcl->IsUsed(10)) continue;  
2952         //
2953         //calcutate parameters
2954         //      
2955         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
2956         // stright track
2957         if (TMath::Abs(yy0)<0.000001) continue;
2958         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
2959         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2960         Double_t r02 = (0.25+y0*y0)*dvertex2;   
2961         //curvature (radius) cut
2962         if (r02<r2min) continue;                
2963        
2964         nin0++; 
2965         //
2966         Double_t c0  = 1/TMath::Sqrt(r02);
2967         if (yy0>0) c0*=-1.;     
2968                
2969        
2970         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
2971         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2972         Double_t dfi0   = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2973         Double_t dfi1   = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);  
2974         //
2975         //
2976         Double_t z0  =  kcl->GetZ();  
2977         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
2978         if (TMath::Abs(zzzz2-z0)>0.5) continue;       
2979         nin1++;              
2980         //      
2981         Double_t dip    = (z1-z0)*c0/dfi1;        
2982         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2983         //
2984         y2 = kcl->GetY(); 
2985         if (dsec==0){
2986           x2 = xx2; 
2987           z2 = kcl->GetZ();       
2988         }
2989         else
2990           {
2991             // rotation 
2992             z2 = kcl->GetZ();  
2993             x2= xx2*cs-y2*sn*dsec;
2994             y2=+xx2*sn*dsec+y2*cs;
2995           }
2996         
2997         x[0] = y1;
2998         x[1] = z1;
2999         x[2] = x0;
3000         x[3] = dip;
3001         x[4] = c0;
3002         //
3003         //
3004         // do we have cluster at the middle ?
3005         Double_t ym,zm;
3006         GetProlongation(x1,xm,x,ym,zm);
3007         UInt_t dummy; 
3008         AliTPCclusterMI * cm=0;
3009         if (TMath::Abs(ym)-ymaxm<0){      
3010           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3011           if ((!cm) || (cm->IsUsed(10))) {        
3012             continue;
3013           }
3014         }
3015         else{     
3016           // rotate y1 to system 0
3017           // get state vector in rotated system 
3018           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
3019           Double_t xr2  =  x0*cs+yr1*sn*dsec;
3020           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3021           //
3022           GetProlongation(xx2,xm,xr,ym,zm);
3023           if (TMath::Abs(ym)-ymaxm<0){
3024             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3025             if ((!cm) || (cm->IsUsed(10))) {      
3026               continue;
3027             }
3028           }
3029         }
3030        
3031
3032         Double_t dym = 0;
3033         Double_t dzm = 0;
3034         if (cm){
3035           dym = ym - cm->GetY();
3036           dzm = zm - cm->GetZ();
3037         }
3038         nin2++;
3039
3040
3041         //
3042         //
3043         Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3044         Double_t sy2=kcl->GetSigmaY2()*2.,     sz2=kcl->GetSigmaZ2()*2.;
3045         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3046         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3047         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3048
3049         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3050         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3051         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3052         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3053         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3054         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3055         
3056         Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3057         Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3058         Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3059         Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3060         
3061         c[0]=sy1;
3062         c[1]=0.;       c[2]=sz1;
3063         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3064         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3065                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3066         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3067         c[13]=f30*sy1*f40+f32*sy2*f42;
3068         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3069         
3070         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3071         
3072         UInt_t index=kr1.GetIndex(is);
3073         AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
3074         
3075         track->fIsSeeding = kTRUE;
3076         track->fSeed1 = i1;
3077         track->fSeed2 = i2;
3078         track->fSeedType=3;
3079
3080        
3081         //if (dsec==0) {
3082           FollowProlongation(*track, (i1+i2)/2,1);
3083           Int_t foundable,found,shared;
3084           track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3085           if ((found<0.55*foundable)  || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3086             seed->Reset();
3087             seed->~AliTPCseed();
3088             continue;
3089           }
3090           //}
3091         
3092         nin++;
3093         FollowProlongation(*track, i2,1);
3094         
3095         
3096         //Int_t rc = 1;
3097         track->fBConstrain =1;
3098         //      track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3099         track->fLastPoint = i1;  // first cluster in track position
3100         track->fFirstPoint = track->fLastPoint;
3101         
3102         if (track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3103             track->GetNumberOfClusters() < track->fNFoundable*0.6 || 
3104             track->fNShared>0.4*track->GetNumberOfClusters() ) {
3105           seed->Reset();
3106           seed->~AliTPCseed();
3107           continue;
3108         }
3109         nout1++;
3110         // Z VERTEX CONDITION
3111         Double_t zv;
3112         zv = track->GetZ()+track->GetTgl()/track->GetC()*
3113           ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3114         if (TMath::Abs(zv-z3)>cuts[2]) {
3115           FollowProlongation(*track, TMath::Max(i2-20,0));
3116           zv = track->GetZ()+track->GetTgl()/track->GetC()*
3117             ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3118           if (TMath::Abs(zv-z3)>cuts[2]){
3119             FollowProlongation(*track, TMath::Max(i2-40,0));
3120             zv = track->GetZ()+track->GetTgl()/track->GetC()*
3121               ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3122             if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
3123               // make seed without constrain
3124               AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3125               FollowProlongation(*track2, i2,1);
3126               track2->fBConstrain = kFALSE;
3127               track2->fSeedType = 1;
3128               arr->AddLast(track2); 
3129               seed->Reset();
3130               seed->~AliTPCseed();
3131               continue;         
3132             }
3133             else{
3134               seed->Reset();
3135               seed->~AliTPCseed();
3136               continue;
3137             
3138             }
3139           }
3140         }
3141         
3142         track->fSeedType =0;
3143         arr->AddLast(track); 
3144         seed = new AliTPCseed;  
3145         nout2++;
3146         // don't consider other combinations
3147         if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
3148           break;
3149       }
3150     }
3151   }
3152   if (fDebug>3){
3153     Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3154   }
3155   delete seed;
3156 }
3157
3158
3159 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3160                                  Float_t deltay) {
3161   
3162
3163
3164   //-----------------------------------------------------------------
3165   // This function creates track seeds.
3166   //-----------------------------------------------------------------
3167   // cuts[0]   - fP4 cut
3168   // cuts[1]   - tan(phi)  cut
3169   // cuts[2]   - zvertex cut
3170   // cuts[3]   - fP3 cut
3171
3172
3173   Int_t nin0  = 0;
3174   Int_t nin1  = 0;
3175   Int_t nin2  = 0;
3176   Int_t nin   = 0;
3177   Int_t nout1 = 0;
3178   Int_t nout2 = 0;
3179   Int_t nout3 =0;
3180   Double_t x[5], c[15];
3181   //
3182   // make temporary seed
3183   AliTPCseed * seed = new AliTPCseed;
3184   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3185   //  Double_t cs=cos(alpha), sn=sin(alpha);
3186   //
3187   //
3188
3189   // first 3 padrows
3190   Double_t x1 = GetXrow(i1-1);
3191   const    AliTPCRow& kr1=GetRow(sec,i1-1);
3192   Double_t y1max  = GetMaxY(i1-1)-kr1.fDeadZone-1.5;  
3193   //
3194   Double_t x1p = GetXrow(i1);
3195   const    AliTPCRow& kr1p=GetRow(sec,i1);
3196   //
3197   Double_t x1m = GetXrow(i1-2);
3198   const    AliTPCRow& kr1m=GetRow(sec,i1-2);
3199
3200   //
3201   //last 3 padrow for seeding
3202   AliTPCRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
3203   Double_t    x3   =  GetXrow(i1-7);
3204   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;  
3205   //
3206   AliTPCRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
3207   Double_t    x3p   = GetXrow(i1-6);
3208   //
3209   AliTPCRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
3210   Double_t    x3m   = GetXrow(i1-8);
3211
3212   //
3213   //
3214   // middle padrow
3215   Int_t im = i1-4;                           //middle pad row index
3216   Double_t xm         = GetXrow(im);         // radius of middle pad-row
3217   const AliTPCRow& krm=GetRow(sec,im);   //middle pad -row
3218   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;  
3219   //
3220   //
3221   Double_t deltax  = x1-x3;
3222   Double_t dymax   = deltax*cuts[1];
3223   Double_t dzmax   = deltax*cuts[3];
3224   //
3225   // loop over clusters  
3226   for (Int_t is=0; is < kr1; is++) {
3227     //
3228     if (kr1[is]->IsUsed(10)) continue;
3229     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
3230     //
3231     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge    
3232     // 
3233     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3234     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3235     //    
3236     Double_t y3,   z3;
3237     //
3238     //
3239     UInt_t index;
3240     for (Int_t js=index1; js < index2; js++) {
3241       const AliTPCclusterMI *kcl = kr3[js];
3242       if (kcl->IsUsed(10)) continue;
3243       y3 = kcl->GetY(); 
3244       // apply angular cuts
3245       if (TMath::Abs(y1-y3)>dymax) continue;
3246       x3 = x3; 
3247       z3 = kcl->GetZ(); 
3248       if (TMath::Abs(z1-z3)>dzmax) continue;
3249       //
3250       Double_t angley = (y1-y3)/(x1-x3);
3251       Double_t anglez = (z1-z3)/(x1-x3);
3252       //
3253       Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3254       Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3255       //
3256       Double_t yyym = angley*(xm-x1)+y1;
3257       Double_t zzzm = anglez*(xm-x1)+z1;
3258
3259       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3260       if (!kcm) continue;
3261       if (kcm->IsUsed(10)) continue;
3262       
3263       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3264       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3265       //
3266       //
3267       //
3268       Int_t used  =0;
3269       Int_t found =0;
3270       //
3271       // look around first
3272       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3273                                                       anglez*(x1m-x1)+z1,
3274                                                       erry,errz,index);
3275       //
3276       if (kc1m){
3277         found++;
3278         if (kc1m->IsUsed(10)) used++;
3279       }
3280       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3281                                                       anglez*(x1p-x1)+z1,
3282                                                       erry,errz,index);
3283       //
3284       if (kc1p){
3285         found++;
3286         if (kc1p->IsUsed(10)) used++;
3287       }
3288       if (used>1)  continue;
3289       if (found<1) continue; 
3290
3291       //
3292       // look around last
3293       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3294                                                       anglez*(x3m-x3)+z3,
3295                                                       erry,errz,index);
3296       //
3297       if (kc3m){
3298         found++;
3299         if (kc3m->IsUsed(10)) used++;
3300       }
3301       else 
3302         continue;
3303       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3304                                                       anglez*(x3p-x3)+z3,
3305                                                       erry,errz,index);
3306       //
3307       if (kc3p){
3308         found++;
3309         if (kc3p->IsUsed(10)) used++;
3310       }
3311       else 
3312         continue;
3313       if (used>1)  continue;
3314       if (found<3) continue;       
3315       //
3316       Double_t x2,y2,z2;
3317       x2 = xm;
3318       y2 = kcm->GetY();
3319       z2 = kcm->GetZ();
3320       //
3321                         
3322       x[0]=y1;
3323       x[1]=z1;
3324       x[4]=F1(x1,y1,x2,y2,x3,y3);
3325       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3326       nin0++;
3327       //
3328       x[2]=F2(x1,y1,x2,y2,x3,y3);
3329       nin1++;
3330       //
3331       x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3332       //if (TMath::Abs(x[3]) > cuts[3]) continue;
3333       nin2++;
3334       //
3335       //
3336       Double_t sy1=0.1,  sz1=0.1;
3337       Double_t sy2=0.1,  sz2=0.1;
3338       Double_t sy3=0.1,  sy=0.1, sz=0.1;
3339       
3340       Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3341       Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3342       Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3343       Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3344       Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3345       Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3346       
3347       Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3348       Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3349       Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3350       Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3351       
3352       c[0]=sy1;
3353       c[1]=0.;       c[2]=sz1; 
3354       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3355       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3356       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3357       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3358       c[13]=f30*sy1*f40+f32*sy2*f42;
3359       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3360       
3361       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3362       
3363       UInt_t index=kr1.GetIndex(is);
3364       AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3365       
3366       track->fIsSeeding = kTRUE;
3367
3368       nin++;      
3369       FollowProlongation(*track, i1-7,1);
3370       if (track->GetNumberOfClusters() < track->fNFoundable*0.75 || 
3371           track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3372         seed->Reset();
3373         seed->~AliTPCseed();
3374         continue;
3375       }
3376       nout1++;
3377       nout2++;  
3378       //Int_t rc = 1;
3379       FollowProlongation(*track, i2,1);
3380       track->fBConstrain =0;
3381       track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3382       track->fFirstPoint = track->fLastPoint;
3383       
3384       if (track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3385           track->GetNumberOfClusters()<track->fNFoundable*0.7 || 
3386           track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3387         seed->Reset();
3388         seed->~AliTPCseed();
3389         continue;
3390       }
3391    
3392       {
3393         FollowProlongation(*track, TMath::Max(i2-10,0),1);
3394         AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3395         FollowProlongation(*track2, i2,1);
3396         track2->fBConstrain = kFALSE;
3397         track2->fSeedType = 4;
3398         arr->AddLast(track2); 
3399         seed->Reset();
3400         seed->~AliTPCseed();
3401       }
3402       
3403    
3404       //arr->AddLast(track); 
3405       //seed = new AliTPCseed;  
3406       nout3++;
3407     }
3408   }
3409   
3410   if (fDebug>3){
3411     Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3412   }
3413   delete seed;
3414 }
3415
3416
3417 //_____________________________________________________________________________
3418 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3419                                  Float_t deltay, Bool_t /*bconstrain*/) {
3420   //-----------------------------------------------------------------
3421   // This function creates track seeds - without vertex constraint
3422   //-----------------------------------------------------------------
3423   // cuts[0]   - fP4 cut        - not applied
3424   // cuts[1]   - tan(phi)  cut
3425   // cuts[2]   - zvertex cut    - not applied 
3426   // cuts[3]   - fP3 cut
3427   Int_t nin0=0;
3428   Int_t nin1=0;
3429   Int_t nin2=0;
3430   Int_t nin3=0;
3431   //  Int_t nin4=0;
3432   //Int_t nin5=0;
3433
3434   
3435
3436   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3437   //  Double_t cs=cos(alpha), sn=sin(alpha);
3438   Int_t row0 = (i1+i2)/2;
3439   Int_t drow = (i1-i2)/2;
3440   const AliTPCRow& kr0=fSectors[sec][row0];
3441   AliTPCRow * kr=0;
3442
3443   AliTPCpolyTrack polytrack;
3444   Int_t nclusters=fSectors[sec][row0];
3445   AliTPCseed * seed = new AliTPCseed;
3446
3447   Int_t sumused=0;
3448   Int_t cused=0;
3449   Int_t cnused=0;
3450   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3451     Int_t nfound =0;
3452     Int_t nfoundable =0;
3453     for (Int_t iter =1; iter<2; iter++){   //iterations
3454       const AliTPCRow& krm=fSectors[sec][row0-iter];
3455       const AliTPCRow& krp=fSectors[sec][row0+iter];      
3456       const AliTPCclusterMI * cl= kr0[is];
3457       
3458       if (cl->IsUsed(10)) {
3459         cused++;
3460       }
3461       else{
3462         cnused++;
3463       }
3464       Double_t x = kr0.GetX();
3465       // Initialization of the polytrack
3466       nfound =0;
3467       nfoundable =0;
3468       polytrack.Reset();
3469       //
3470       Double_t y0= cl->GetY();
3471       Double_t z0= cl->GetZ();
3472       Float_t erry = 0;
3473       Float_t errz = 0;
3474       
3475       Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3476       if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;  // seed only at the edge
3477       
3478       erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;      
3479       errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;      
3480       polytrack.AddPoint(x,y0,z0,erry, errz);
3481
3482       sumused=0;
3483       if (cl->IsUsed(10)) sumused++;
3484
3485
3486       Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3487       Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3488       //
3489       x = krm.GetX();
3490       AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3491       if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3492         erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;          
3493         errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3494         if (cl1->IsUsed(10))  sumused++;
3495         polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3496       }
3497       //
3498       x = krp.GetX();
3499       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3500       if (cl2) {
3501         erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;          
3502         errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3503         if (cl2->IsUsed(10)) sumused++;  
3504         polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3505       }
3506       //
3507       if (sumused>0) continue;
3508       nin0++;
3509       polytrack.UpdateParameters();
3510       // follow polytrack
3511       roadz = 1.2;
3512       roady = 1.2;
3513       //
3514       Double_t yn,zn;
3515       nfoundable = polytrack.GetN();
3516       nfound     = nfoundable; 
3517       //
3518       for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3519         Float_t maxdist = 0.8*(1.+3./(ddrow));
3520         for (Int_t delta = -1;delta<=1;delta+=2){
3521           Int_t row = row0+ddrow*delta;
3522           kr = &(fSectors[sec][row]);
3523           Double_t xn = kr->GetX();
3524           Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3525           polytrack.GetFitPoint(xn,yn,zn);
3526           if (TMath::Abs(yn)>ymax) continue;
3527           nfoundable++;
3528           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3529           if (cln) {
3530             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3531             if (dist<maxdist){
3532               /*
3533               erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));         
3534               errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3535               if (cln->IsUsed(10)) {
3536                 //      printf("used\n");
3537                 sumused++;
3538                 erry*=2;
3539                 errz*=2;
3540               }
3541               */
3542               erry=0.1;
3543               errz=0.1;
3544               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3545               nfound++;
3546             }
3547           }
3548         }
3549         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;     
3550         polytrack.UpdateParameters();
3551       }           
3552     }
3553     if ( (sumused>3) || (sumused>0.5*nfound))  {
3554       //printf("sumused   %d\n",sumused);
3555       continue;
3556     }
3557     nin1++;
3558     Double_t dy,dz;
3559     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3560     AliTPCpolyTrack track2;
3561     
3562     polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3563     if (track2.GetN()<0.5*nfoundable) continue;
3564     nin2++;
3565
3566     if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3567       //
3568       // test seed with and without constrain
3569       for (Int_t constrain=0; constrain<=0;constrain++){
3570         // add polytrack candidate
3571
3572         Double_t x[5], c[15];
3573         Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3574         track2.GetBoundaries(x3,x1);    
3575         x2 = (x1+x3)/2.;
3576         track2.GetFitPoint(x1,y1,z1);
3577         track2.GetFitPoint(x2,y2,z2);
3578         track2.GetFitPoint(x3,y3,z3);
3579         //
3580         //is track pointing to the vertex ?
3581         Double_t x0,y0,z0;
3582         x0=0;
3583         polytrack.GetFitPoint(x0,y0,z0);
3584
3585         if (constrain) {
3586           x2 = x3;
3587           y2 = y3;
3588           z2 = z3;
3589           
3590           x3 = 0;
3591           y3 = 0;
3592           z3 = 0;
3593         }
3594         x[0]=y1;
3595         x[1]=z1;
3596         x[4]=F1(x1,y1,x2,y2,x3,y3);
3597                 
3598         //      if (TMath::Abs(x[4]) >= cuts[0]) continue;  //
3599         x[2]=F2(x1,y1,x2,y2,x3,y3);
3600         
3601         //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3602         //x[3]=F3(x1,y1,x2,y2,z1,z2);
3603         x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3604         //if (TMath::Abs(x[3]) > cuts[3]) continue;
3605
3606         
3607         Double_t sy =0.1, sz =0.1;
3608         Double_t sy1=0.02, sz1=0.02;
3609         Double_t sy2=0.02, sz2=0.02;
3610         Double_t sy3=0.02;
3611
3612         if (constrain){
3613           sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3614         }
3615         
3616         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3617         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3618         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3619         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3620         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3621         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3622
3623         Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3624         Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3625         Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3626         Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3627
3628         
3629         c[0]=sy1;
3630         c[1]=0.;       c[2]=sz1;
3631         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3632         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3633         c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3634         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3635         c[13]=f30*sy1*f40+f32*sy2*f42;
3636         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3637         
3638         //Int_t row1 = fSectors->GetRowNumber(x1);
3639         Int_t row1 = GetRowNumber(x1);
3640
3641         UInt_t index=0;
3642         //kr0.GetIndex(is);
3643         AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);