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