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