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