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