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