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