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