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