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