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