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