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