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