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