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