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