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