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