Adding fClusterMap and fSharedMap of TPC clusters for the tracks
[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   AliTPCTransform trafo;
1249
1250   Double_t x[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1251   Int_t i[1]={cluster->GetDetector()};
1252
1253   trafo.Transform(x,i,0,1);
1254
1255   cluster->SetX(x[0]);
1256   cluster->SetY(x[1]);
1257   cluster->SetZ(x[2]);
1258
1259   return;
1260
1261   // The old stuff:
1262
1263   //
1264   // 
1265   //
1266   //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1267   TGeoHMatrix  *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1268   //TGeoHMatrix  mat;
1269   Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1270   Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1271   if (mat) mat->LocalToMaster(pos,posC);
1272   else{
1273     // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1274   }
1275   //cluster->SetX(posC[0]);
1276   //cluster->SetY(posC[1]);
1277   //cluster->SetZ(posC[2]);
1278 }
1279
1280 //_____________________________________________________________________________
1281 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1282   //-----------------------------------------------------------------
1283   // This function fills outer TPC sectors with clusters.
1284   //-----------------------------------------------------------------
1285   Int_t nrows = fOuterSec->GetNRows();
1286   UInt_t index=0;
1287   for (Int_t sec = 0;sec<fkNOS;sec++)
1288     for (Int_t row = 0;row<nrows;row++){
1289       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);  
1290       Int_t sec2 = sec+2*fkNIS;
1291       //left
1292       Int_t ncl = tpcrow->GetN1();
1293       while (ncl--) {
1294         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1295         index=(((sec2<<8)+row)<<16)+ncl;
1296         tpcrow->InsertCluster(c,index);
1297       }
1298       //right
1299       ncl = tpcrow->GetN2();
1300       while (ncl--) {
1301         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1302         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1303         tpcrow->InsertCluster(c,index);
1304       }
1305       //
1306       // write indexes for fast acces
1307       //
1308       for (Int_t i=0;i<510;i++)
1309         tpcrow->SetFastCluster(i,-1);
1310       for (Int_t i=0;i<tpcrow->GetN();i++){
1311         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1312         tpcrow->SetFastCluster(zi,i);  // write index
1313       }
1314       Int_t last = 0;
1315       for (Int_t i=0;i<510;i++){
1316         if (tpcrow->GetFastCluster(i)<0)
1317           tpcrow->SetFastCluster(i,last);
1318         else
1319           last = tpcrow->GetFastCluster(i);
1320       }
1321     }  
1322   fN=fkNOS;
1323   fSectors=fOuterSec;
1324   return 0;
1325 }
1326
1327
1328 //_____________________________________________________________________________
1329 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1330   //-----------------------------------------------------------------
1331   // This function fills inner TPC sectors with clusters.
1332   //-----------------------------------------------------------------
1333   Int_t nrows = fInnerSec->GetNRows();
1334   UInt_t index=0;
1335   for (Int_t sec = 0;sec<fkNIS;sec++)
1336     for (Int_t row = 0;row<nrows;row++){
1337       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1338       //
1339       //left
1340       Int_t ncl = tpcrow->GetN1();
1341       while (ncl--) {
1342         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1343         index=(((sec<<8)+row)<<16)+ncl;
1344         tpcrow->InsertCluster(c,index);
1345       }
1346       //right
1347       ncl = tpcrow->GetN2();
1348       while (ncl--) {
1349         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1350         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1351         tpcrow->InsertCluster(c,index);
1352       }
1353       //
1354       // write indexes for fast acces
1355       //
1356       for (Int_t i=0;i<510;i++)
1357         tpcrow->SetFastCluster(i,-1);
1358       for (Int_t i=0;i<tpcrow->GetN();i++){
1359         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1360         tpcrow->SetFastCluster(zi,i);  // write index
1361       }
1362       Int_t last = 0;
1363       for (Int_t i=0;i<510;i++){
1364         if (tpcrow->GetFastCluster(i)<0)
1365           tpcrow->SetFastCluster(i,last);
1366         else
1367           last = tpcrow->GetFastCluster(i);
1368       }
1369
1370     }  
1371    
1372   fN=fkNIS;
1373   fSectors=fInnerSec;
1374   return 0;
1375 }
1376
1377
1378
1379 //_________________________________________________________________________
1380 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1381   //--------------------------------------------------------------------
1382   //       Return pointer to a given cluster
1383   //--------------------------------------------------------------------
1384   if (index<0) return 0; // no cluster
1385   Int_t sec=(index&0xff000000)>>24; 
1386   Int_t row=(index&0x00ff0000)>>16; 
1387   Int_t ncl=(index&0x00007fff)>>00;
1388
1389   const AliTPCRow * tpcrow=0;
1390   AliTPCclusterMI * clrow =0;
1391
1392   if (sec<0 || sec>=fkNIS*4) {
1393     AliWarning(Form("Wrong sector %d",sec));
1394     return 0x0;
1395   }
1396
1397   if (sec<fkNIS*2){
1398     tpcrow = &(fInnerSec[sec%fkNIS][row]);
1399     if (tpcrow==0) return 0;
1400
1401     if (sec<fkNIS) {
1402       if (tpcrow->GetN1()<=ncl) return 0;
1403       clrow = tpcrow->GetClusters1();
1404     }
1405     else {
1406       if (tpcrow->GetN2()<=ncl) return 0;
1407       clrow = tpcrow->GetClusters2();
1408     }
1409   }
1410   else {
1411     tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1412     if (tpcrow==0) return 0;
1413
1414     if (sec-2*fkNIS<fkNOS) {
1415       if (tpcrow->GetN1()<=ncl) return 0;
1416       clrow = tpcrow->GetClusters1();
1417     }
1418     else {
1419       if (tpcrow->GetN2()<=ncl) return 0;
1420       clrow = tpcrow->GetClusters2();
1421     }
1422   }
1423
1424   return &(clrow[ncl]);      
1425   
1426 }
1427
1428
1429
1430 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1431   //-----------------------------------------------------------------
1432   // This function tries to find a track prolongation to next pad row
1433   //-----------------------------------------------------------------
1434   //
1435   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1436   AliTPCclusterMI *cl=0;
1437   Int_t tpcindex= t.GetClusterIndex2(nr);
1438   //
1439   // update current shape info every 5 pad-row
1440   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1441     GetShape(&t,nr);    
1442     //}
1443   //  
1444   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters 
1445     //        
1446     if (tpcindex==-1) return 0; //track in dead zone
1447     if (tpcindex>0){     //
1448       cl = t.GetClusterPointer(nr);
1449       if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1450       t.SetCurrentClusterIndex1(tpcindex); 
1451     }
1452     if (cl){      
1453       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
1454       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1455       //
1456       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1457       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1458       
1459       if (TMath::Abs(angle-t.GetAlpha())>0.001){
1460         Double_t rotation = angle-t.GetAlpha();
1461         t.SetRelativeSector(relativesector);
1462         if (!t.Rotate(rotation)) return 0;      
1463       }
1464       if (!t.PropagateTo(x)) return 0;
1465       //
1466       t.SetCurrentCluster(cl); 
1467       t.SetRow(nr);
1468       Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1469       if ((tpcindex&0x8000)==0) accept =0;
1470       if (accept<3) { 
1471         //if founded cluster is acceptible
1472         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
1473           t.SetErrorY2(t.GetErrorY2()+0.03);
1474           t.SetErrorZ2(t.GetErrorZ2()+0.03); 
1475           t.SetErrorY2(t.GetErrorY2()*3);
1476           t.SetErrorZ2(t.GetErrorZ2()*3); 
1477         }
1478         t.SetNFoundable(t.GetNFoundable()+1);
1479         UpdateTrack(&t,accept);
1480         return 1;
1481       }    
1482     }
1483   }
1484   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;  // cut on angle
1485   if (fIteration>1){
1486     // not look for new cluster during refitting    
1487     t.SetNFoundable(t.GetNFoundable()+1);
1488     return 0;
1489   }
1490   //
1491   UInt_t index=0;
1492   //  if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1493   Double_t  y=t.GetYat(x);
1494   if (TMath::Abs(y)>ymax){
1495     if (y > ymax) {
1496       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1497       if (!t.Rotate(fSectors->GetAlpha())) 
1498         return 0;
1499     } else if (y <-ymax) {
1500       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1501       if (!t.Rotate(-fSectors->GetAlpha())) 
1502         return 0;
1503     }
1504     //return 1;
1505   }
1506   //
1507   if (!t.PropagateTo(x)) {
1508     if (fIteration==0) t.SetRemoval(10);
1509     return 0;
1510   }
1511   y=t.GetY(); 
1512   Double_t z=t.GetZ();
1513   //
1514
1515   if (!IsActive(t.GetRelativeSector(),nr)) {
1516     t.SetInDead(kTRUE);
1517     t.SetClusterIndex2(nr,-1); 
1518     return 0;
1519   }
1520   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1521   Bool_t isActive  = IsActive(t.GetRelativeSector(),nr);
1522   Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1523   
1524   if (!isActive || !isActive2) return 0;
1525
1526   const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1527   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1528   Double_t  roady  =1.;
1529   Double_t  roadz = 1.;
1530   //
1531   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1532     t.SetInDead(kTRUE);
1533     t.SetClusterIndex2(nr,-1); 
1534     return 0;
1535   } 
1536   else
1537     {
1538       if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) 
1539         t.SetNFoundable(t.GetNFoundable()+1);
1540       else
1541         return 0;
1542     }   
1543   //calculate 
1544   if (krow) {
1545     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);    
1546     cl = krow.FindNearest2(y,z,roady,roadz,index);    
1547     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));       
1548   }  
1549   if (cl) {
1550     t.SetCurrentCluster(cl); 
1551     t.SetRow(nr);
1552     if (fIteration==2&&cl->IsUsed(10)) return 0; 
1553     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1554     if (fIteration==2&&cl->IsUsed(11)) {
1555       t.SetErrorY2(t.GetErrorY2()+0.03);
1556       t.SetErrorZ2(t.GetErrorZ2()+0.03); 
1557       t.SetErrorY2(t.GetErrorY2()*3);
1558       t.SetErrorZ2(t.GetErrorZ2()*3); 
1559     }
1560     /*    
1561     if (t.fCurrentCluster->IsUsed(10)){
1562       //
1563       //     
1564
1565       t.fNShared++;
1566       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1567         t.fRemoval =10;
1568         return 0;
1569       }
1570     }
1571     */
1572     if (accept<3) UpdateTrack(&t,accept);
1573
1574   } else {  
1575     if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1576     
1577   }
1578   return 1;
1579 }
1580
1581 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1582   //-----------------------------------------------------------------
1583   // This function tries to find a track prolongation to next pad row
1584   //-----------------------------------------------------------------
1585   //
1586   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1587   Double_t y,z; 
1588   if (!t.GetProlongation(x,y,z)) {
1589     t.SetRemoval(10);
1590     return 0;
1591   }
1592   //
1593   //
1594   if (TMath::Abs(y)>ymax){
1595     
1596     if (y > ymax) {
1597       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1598       if (!t.Rotate(fSectors->GetAlpha())) 
1599         return 0;
1600     } else if (y <-ymax) {
1601       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1602       if (!t.Rotate(-fSectors->GetAlpha())) 
1603         return 0;
1604     }
1605     if (!t.PropagateTo(x)) {
1606       return 0;
1607     } 
1608     t.GetProlongation(x,y,z);
1609   }
1610   //
1611   // update current shape info every 3 pad-row
1612   if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1613     //    t.fCurrentSigmaY = GetSigmaY(&t);
1614     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1615     GetShape(&t,nr);
1616   }
1617   //  
1618   AliTPCclusterMI *cl=0;
1619   UInt_t index=0;
1620   
1621   
1622   //Int_t nr2 = nr;
1623   const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1624   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1625   Double_t  roady  =1.;
1626   Double_t  roadz = 1.;
1627   //
1628   Int_t row = nr;
1629   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1630     t.SetInDead(kTRUE);
1631     t.SetClusterIndex2(row,-1); 
1632     return 0;
1633   } 
1634   else
1635     {
1636       if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1637     }   
1638   //calculate 
1639   
1640   if ((cl==0)&&(krow)) {
1641     //    cl = krow.FindNearest2(y+10,z,roady,roadz,index);    
1642     cl = krow.FindNearest2(y,z,roady,roadz,index);    
1643
1644     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));       
1645   }  
1646
1647   if (cl) {
1648     t.SetCurrentCluster(cl); 
1649     //    Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);        
1650     //if (accept<3){
1651       t.SetClusterIndex2(row,index);
1652       t.SetClusterPointer(row, cl);
1653       //}
1654   }
1655   return 1;
1656 }
1657
1658
1659
1660 //_________________________________________________________________________
1661 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1662 {
1663   // Get track space point by index
1664   // return false in case the cluster doesn't exist
1665   AliTPCclusterMI *cl = GetClusterMI(index);
1666   if (!cl) return kFALSE;
1667   Int_t sector = (index&0xff000000)>>24;
1668   //  Int_t row = (index&0x00ff0000)>>16;
1669   Float_t xyz[3];
1670   //  xyz[0] = fParam->GetPadRowRadii(sector,row);
1671   xyz[0] = cl->GetX();
1672   xyz[1] = cl->GetY();
1673   xyz[2] = cl->GetZ();
1674   Float_t sin,cos;
1675   fParam->AdjustCosSin(sector,cos,sin);
1676   Float_t x = cos*xyz[0]-sin*xyz[1];
1677   Float_t y = cos*xyz[1]+sin*xyz[0];
1678   Float_t cov[6];
1679   Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1680   if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1681   Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1682   if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1683   cov[0] = sin*sin*sigmaY2;
1684   cov[1] = -sin*cos*sigmaY2;
1685   cov[2] = 0.;
1686   cov[3] = cos*cos*sigmaY2;
1687   cov[4] = 0.;
1688   cov[5] = sigmaZ2;
1689   p.SetXYZ(x,y,xyz[2],cov);
1690   AliGeomManager::ELayerID iLayer;
1691   Int_t idet;
1692   if (sector < fParam->GetNInnerSector()) {
1693     iLayer = AliGeomManager::kTPC1;
1694     idet = sector;
1695   }
1696   else {
1697     iLayer = AliGeomManager::kTPC2;
1698     idet = sector - fParam->GetNInnerSector();
1699   }
1700   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1701   p.SetVolumeID(volid);
1702   return kTRUE;
1703 }
1704
1705
1706
1707 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
1708   //-----------------------------------------------------------------
1709   // This function tries to find a track prolongation to next pad row
1710   //-----------------------------------------------------------------
1711   t.SetCurrentCluster(0);
1712   t.SetCurrentClusterIndex1(0);   
1713    
1714   Double_t xt=t.GetX();
1715   Int_t     row = GetRowNumber(xt)-1; 
1716   Double_t  ymax= GetMaxY(nr);
1717
1718   if (row < nr) return 1; // don't prolongate if not information until now -
1719 //   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1720 //     t.fRemoval =10;
1721 //     return 0;  // not prolongate strongly inclined tracks
1722 //   } 
1723 //   if (TMath::Abs(t.GetSnp())>0.95) {
1724 //     t.fRemoval =10;
1725 //     return 0;  // not prolongate strongly inclined tracks
1726 //   }// patch 28 fev 06
1727
1728   Double_t x= GetXrow(nr);
1729   Double_t y,z;
1730   //t.PropagateTo(x+0.02);
1731   //t.PropagateTo(x+0.01);
1732   if (!t.PropagateTo(x)){
1733     return 0;
1734   }
1735   //
1736   y=t.GetY();
1737   z=t.GetZ();
1738
1739   if (TMath::Abs(y)>ymax){
1740     if (y > ymax) {
1741       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1742       if (!t.Rotate(fSectors->GetAlpha())) 
1743         return 0;
1744     } else if (y <-ymax) {
1745       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1746       if (!t.Rotate(-fSectors->GetAlpha())) 
1747         return 0;
1748     }
1749     //    if (!t.PropagateTo(x)){
1750     //  return 0;
1751     //}
1752     return 1;
1753     //y = t.GetY();    
1754   }
1755   //
1756   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1757
1758   if (!IsActive(t.GetRelativeSector(),nr)) {
1759     t.SetInDead(kTRUE);
1760     t.SetClusterIndex2(nr,-1); 
1761     return 0;
1762   }
1763   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1764
1765   AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1766
1767   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1768     t.SetInDead(kTRUE);
1769     t.SetClusterIndex2(nr,-1); 
1770     return 0;
1771   } 
1772   else
1773     {
1774       if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1775       else
1776         return 0;      
1777     }
1778
1779   // update current
1780   if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1781     //    t.fCurrentSigmaY = GetSigmaY(&t);
1782     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1783     GetShape(&t,nr);
1784   }
1785     
1786   AliTPCclusterMI *cl=0;
1787   Int_t index=0;
1788   //
1789   Double_t roady = 1.;
1790   Double_t roadz = 1.;
1791   //
1792
1793   if (!cl){
1794     index = t.GetClusterIndex2(nr);    
1795     if ( (index>0) && (index&0x8000)==0){
1796       cl = t.GetClusterPointer(nr);
1797       if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1798       t.SetCurrentClusterIndex1(index);
1799       if (cl) {
1800         t.SetCurrentCluster(cl);
1801         return 1;
1802       }
1803     }
1804   }
1805
1806   //  if (index<0) return 0;
1807   UInt_t uindex = TMath::Abs(index);
1808
1809   if (krow) {    
1810     //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);      
1811     cl = krow.FindNearest2(y,z,roady,roadz,uindex);      
1812   }
1813
1814   if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));   
1815   t.SetCurrentCluster(cl);
1816
1817   return 1;
1818 }
1819
1820
1821 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1822   //-----------------------------------------------------------------
1823   // This function tries to find a track prolongation to next pad row
1824   //-----------------------------------------------------------------
1825
1826   //update error according neighborhoud
1827
1828   if (t.GetCurrentCluster()) {
1829     t.SetRow(nr); 
1830     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1831     
1832     if (t.GetCurrentCluster()->IsUsed(10)){
1833       //
1834       //
1835       //  t.fErrorZ2*=2;
1836       //  t.fErrorY2*=2;
1837       t.SetNShared(t.GetNShared()+1);
1838       if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1839         t.SetRemoval(10);
1840         return 0;
1841       }
1842     }   
1843     if (fIteration>0) accept = 0;
1844     if (accept<3)  UpdateTrack(&t,accept);  
1845  
1846   } else {
1847     if (fIteration==0){
1848       if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);      
1849       if (  t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);      
1850
1851       if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1852     }
1853   }
1854   return 1;
1855 }
1856
1857
1858
1859 //_____________________________________________________________________________
1860 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1861   //-----------------------------------------------------------------
1862   // This function tries to find a track prolongation.
1863   //-----------------------------------------------------------------
1864   Double_t xt=t.GetX();
1865   //
1866   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1867   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
1868   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
1869   //
1870   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1871     
1872   Int_t first = GetRowNumber(xt)-1;
1873   for (Int_t nr= first; nr>=rf; nr-=step) {
1874     // update kink info
1875     if (t.GetKinkIndexes()[0]>0){
1876       for (Int_t i=0;i<3;i++){
1877         Int_t index = t.GetKinkIndexes()[i];
1878         if (index==0) break;
1879         if (index<0) continue;
1880         //
1881         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1882         if (!kink){
1883           printf("PROBLEM\n");
1884         }
1885         else{
1886           Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1887           if (kinkrow==nr){
1888             AliExternalTrackParam paramd(t);
1889             kink->SetDaughter(paramd);
1890             kink->SetStatus(2,5);
1891             kink->Update();
1892           }
1893         }
1894       }
1895     }
1896
1897     if (nr==80) t.UpdateReference();
1898     if (nr<fInnerSec->GetNRows()) 
1899       fSectors = fInnerSec;
1900     else
1901       fSectors = fOuterSec;
1902     if (FollowToNext(t,nr)==0) 
1903       if (!t.IsActive()) 
1904         return 0;
1905     
1906   }   
1907   return 1;
1908 }
1909
1910
1911 //_____________________________________________________________________________
1912 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1913   //-----------------------------------------------------------------
1914   // This function tries to find a track prolongation.
1915   //-----------------------------------------------------------------
1916   Double_t xt=t.GetX();
1917   //
1918   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1919   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
1920   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
1921   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1922     
1923   for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1924     
1925     if (FollowToNextFast(t,nr)==0) 
1926       if (!t.IsActive()) return 0;
1927     
1928   }   
1929   return 1;
1930 }
1931
1932
1933
1934
1935
1936 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1937   //-----------------------------------------------------------------
1938   // This function tries to find a track prolongation.
1939   //-----------------------------------------------------------------
1940   //
1941   Double_t xt=t.GetX();  
1942   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1943   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
1944   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
1945   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1946     
1947   Int_t first = t.GetFirstPoint();
1948   if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1949   //
1950   if (first<0) first=0;
1951   for (Int_t nr=first; nr<=rf; nr++) {
1952     //    if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1953     if (t.GetKinkIndexes()[0]<0){
1954       for (Int_t i=0;i<3;i++){
1955         Int_t index = t.GetKinkIndexes()[i];
1956         if (index==0) break;
1957         if (index>0) continue;
1958         index = TMath::Abs(index);
1959         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1960         if (!kink){
1961           printf("PROBLEM\n");
1962         }
1963         else{
1964           Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1965           if (kinkrow==nr){
1966             AliExternalTrackParam paramm(t);
1967             kink->SetMother(paramm);
1968             kink->SetStatus(2,1);
1969             kink->Update();
1970           }
1971         }
1972       }      
1973     }
1974     //
1975     if (nr<fInnerSec->GetNRows()) 
1976       fSectors = fInnerSec;
1977     else
1978       fSectors = fOuterSec;
1979
1980     FollowToNext(t,nr);                                                             
1981   }   
1982   return 1;
1983 }
1984
1985
1986
1987
1988    
1989 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1990 {
1991   //
1992   //
1993   sum1=0;
1994   sum2=0;
1995   Int_t sum=0;
1996   //
1997   Float_t dz2 =(s1->GetZ() - s2->GetZ());
1998   dz2*=dz2;  
1999
2000   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2001   dy2*=dy2;
2002   Float_t distance = TMath::Sqrt(dz2+dy2);
2003   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2004  
2005   //  Int_t offset =0;
2006   Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2007   Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2008   if (lastpoint>160) 
2009     lastpoint =160;
2010   if (firstpoint<0) 
2011     firstpoint = 0;
2012   if (firstpoint>lastpoint) {
2013     firstpoint =lastpoint;
2014     //    lastpoint  =160;
2015   }
2016     
2017   
2018   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2019     if (s1->GetClusterIndex2(i)>0) sum1++;
2020     if (s2->GetClusterIndex2(i)>0) sum2++;
2021     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2022       sum++;
2023     }
2024   }
2025   if (sum<5) return 0;
2026
2027   Float_t summin = TMath::Min(sum1+1,sum2+1);
2028   Float_t ratio = (sum+1)/Float_t(summin);
2029   return ratio;
2030 }
2031
2032 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2033 {
2034   //
2035   //
2036   Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2037   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2038   Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2039   Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2040   
2041   //
2042   Int_t sumshared=0;
2043   //
2044   //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2045   //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2046   Int_t firstpoint = 0;
2047   Int_t lastpoint = 160;
2048   //
2049   //  if (firstpoint>=lastpoint-5) return;;
2050
2051   for (Int_t i=firstpoint;i<lastpoint;i++){
2052     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2053     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2054       sumshared++;
2055       s1->SetSharedMapBit(i, kTRUE);
2056       s2->SetSharedMapBit(i, kTRUE);
2057     }
2058     if (s1->GetClusterIndex2(i)>0)
2059       s1->SetClusterMapBit(i, kTRUE);
2060   }
2061   if (sumshared>cutN0){
2062     // sign clusters
2063     //
2064     for (Int_t i=firstpoint;i<lastpoint;i++){
2065       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2066       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2067         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2068         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);; 
2069         if (s1->IsActive()&&s2->IsActive()){
2070           p1->SetShared(kTRUE);
2071           p2->SetShared(kTRUE);
2072         }       
2073       }
2074     }
2075   }
2076   //  
2077   if (sumshared>cutN0){
2078     for (Int_t i=0;i<4;i++){
2079       if (s1->GetOverlapLabel(3*i)==0){
2080         s1->SetOverlapLabel(3*i,  s2->GetLabel());
2081         s1->SetOverlapLabel(3*i+1,sumshared);
2082         s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2083         break;
2084       } 
2085     }
2086     for (Int_t i=0;i<4;i++){
2087       if (s2->GetOverlapLabel(3*i)==0){
2088         s2->SetOverlapLabel(3*i,  s1->GetLabel());
2089         s2->SetOverlapLabel(3*i+1,sumshared);
2090         s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2091         break;
2092       } 
2093     }    
2094   }
2095 }
2096
2097 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
2098 {
2099   //
2100   //sort trackss according sectors
2101   //  
2102   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2103     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2104     if (!pt) continue;
2105     //if (pt) RotateToLocal(pt);
2106     pt->SetSort(0);
2107   }
2108   arr->UnSort();
2109   arr->Sort();  // sorting according z
2110   arr->Expand(arr->GetEntries());
2111   //
2112   //
2113   Int_t nseed=arr->GetEntriesFast();
2114   for (Int_t i=0; i<nseed; i++) {
2115     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2116     if (!pt) continue;
2117     for (Int_t j=0;j<=12;j++){
2118       pt->SetOverlapLabel(j,0);
2119     }
2120   }
2121   for (Int_t i=0; i<nseed; i++) {
2122     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2123     if (!pt) continue;
2124     if (pt->GetRemoval()>10) continue;
2125     for (Int_t j=i+1; j<nseed; j++){
2126       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2127       //      if (pt2){
2128       if (pt2->GetRemoval()<=10) {
2129         if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2130         SignShared(pt,pt2);
2131       }
2132     }  
2133   }
2134 }
2135
2136 void  AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2,  Int_t removalindex)
2137 {
2138   //
2139   //sort trackss according sectors
2140   //
2141   if (fDebug&1) {
2142     Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries());
2143   }
2144   //
2145   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2146     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2147     if (!pt) continue;
2148     pt->SetSort(0);
2149   }
2150   arr->UnSort();
2151   arr->Sort();  // sorting according z
2152   arr->Expand(arr->GetEntries());
2153   //
2154   //reset overlap labels
2155   //
2156   Int_t nseed=arr->GetEntriesFast();
2157   for (Int_t i=0; i<nseed; i++) {
2158     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2159     if (!pt) continue;
2160     pt->SetUniqueID(i);
2161     for (Int_t j=0;j<=12;j++){
2162       pt->SetOverlapLabel(j,0);
2163     }
2164   }
2165   //
2166   //sign shared tracks
2167   for (Int_t i=0; i<nseed; i++) {
2168     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2169     if (!pt) continue;
2170     if (pt->GetRemoval()>10) continue;
2171     Float_t deltac = pt->GetC()*0.1;
2172     for (Int_t j=i+1; j<nseed; j++){
2173       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2174       //      if (pt2){
2175       if (pt2->GetRemoval()<=10) {
2176         if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2177         if (TMath::Abs(pt->GetC()  -pt2->GetC())>deltac) continue;
2178         if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2179         //
2180         SignShared(pt,pt2);
2181       }
2182     }
2183   }
2184   //
2185   // remove highly shared tracks
2186   for (Int_t i=0; i<nseed; i++) {
2187     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2188     if (!pt) continue;
2189     if (pt->GetRemoval()>10) continue;
2190     //
2191     Int_t sumshared =0;
2192     for (Int_t j=0;j<4;j++){
2193       sumshared = pt->GetOverlapLabel(j*3+1);      
2194     }
2195     Float_t factor = factor1;
2196     if (pt->GetRemoval()>0) factor = factor2;
2197     if (sumshared/pt->GetNumberOfClusters()>factor){
2198       for (Int_t j=0;j<4;j++){
2199         if (pt->GetOverlapLabel(3*j)==0) continue;
2200         if (pt->GetOverlapLabel(3*j+1)<5) continue; 
2201         if (pt->GetRemoval()==removalindex) continue;      
2202         AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->GetOverlapLabel(3*j+2));
2203         if (!pt2) continue;
2204         if (pt2->GetSigma2C()<pt->GetSigma2C()){
2205           //      pt->fRemoval = removalindex;
2206           delete arr->RemoveAt(i);        
2207           break;
2208         }
2209       }      
2210     }
2211   }
2212   arr->Compress();
2213   if (fDebug&1) {
2214     Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries());
2215   }
2216 }
2217
2218
2219
2220
2221
2222
2223 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2224 {
2225   //
2226   //sort tracks in array according mode criteria
2227   Int_t nseed = arr->GetEntriesFast();    
2228   for (Int_t i=0; i<nseed; i++) {
2229     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2230     if (!pt) {
2231       continue;
2232     }
2233     pt->SetSort(mode);
2234   }
2235   arr->UnSort();
2236   arr->Sort();
2237 }
2238
2239 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t removalindex)
2240 {
2241
2242   //Loop over all tracks and remove "overlaps"
2243   //
2244   //
2245   Int_t nseed = arr->GetEntriesFast();  
2246   Int_t good =0;
2247
2248   for (Int_t i=0; i<nseed; i++) {
2249     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2250     if (!pt) {
2251       delete arr->RemoveAt(i);
2252     }
2253     else{
2254       pt->SetSort(1);
2255       pt->SetBSigned(kFALSE);
2256     }
2257   }
2258   arr->Compress();
2259   nseed = arr->GetEntriesFast();
2260   arr->UnSort();
2261   arr->Sort();
2262   //
2263   //unsign used
2264   UnsignClusters();
2265   //
2266   for (Int_t i=0; i<nseed; i++) {
2267     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2268     if (!pt) {
2269       continue;
2270     }    
2271     Int_t found,foundable,shared;
2272     if (pt->IsActive()) 
2273       pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2274     else
2275       pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE); 
2276     //
2277     Double_t factor = factor2;
2278     if (pt->GetBConstrain()) factor = factor1;
2279
2280     if ((Float_t(shared)/Float_t(found))>factor){
2281       pt->Desactivate(removalindex);
2282       continue;
2283     }
2284
2285     good++;
2286     for (Int_t i=0; i<160; i++) {
2287       Int_t index=pt->GetClusterIndex2(i);
2288       if (index<0 || index&0x8000 ) continue;
2289       AliTPCclusterMI *c= pt->GetClusterPointer(i);        
2290       if (!c) continue;
2291       //      if (!c->IsUsed(10)) c->Use(10);
2292       //if (pt->IsActive()) 
2293       c->Use(10);  
2294       //else
2295       //        c->Use(5);
2296     }
2297     
2298   }
2299   fNtracks = good;
2300   if (fDebug>0){
2301     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2302   }
2303 }
2304
2305
2306 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2307 {
2308
2309   //Loop over all tracks and remove "overlaps"
2310   //
2311   //
2312   UnsignClusters();
2313   //
2314   Int_t nseed = arr->GetEntriesFast();  
2315   Float_t * quality = new Float_t[nseed];
2316   Int_t   * indexes = new Int_t[nseed];
2317   Int_t good =0;
2318   //
2319   //
2320   for (Int_t i=0; i<nseed; i++) {
2321     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2322     if (!pt){
2323       quality[i]=-1;
2324       continue;
2325     }
2326     pt->UpdatePoints();    //select first last max dens points
2327     Float_t * points = pt->GetPoints();
2328     if (points[3]<0.8) quality[i] =-1;
2329     //
2330     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2331   }
2332   TMath::Sort(nseed,quality,indexes);
2333   //
2334   //
2335   for (Int_t itrack=0; itrack<nseed; itrack++) {
2336     Int_t trackindex = indexes[itrack];
2337     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);    
2338     if (quality[trackindex]<0){
2339       if (pt) {
2340         delete arr->RemoveAt(trackindex);
2341       }
2342       else{
2343         arr->RemoveAt(trackindex);
2344       }
2345       continue;
2346     }
2347     //
2348     Int_t first = Int_t(pt->GetPoints()[0]);
2349     Int_t last  = Int_t(pt->GetPoints()[2]);
2350     Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2351     //
2352     Int_t found,foundable,shared;
2353     pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
2354     Float_t sharedfactor = Float_t(shared+1)/Float_t(found+1);
2355     Bool_t itsgold =kFALSE;
2356     if (pt->GetESD()){
2357       Int_t dummy[12];
2358       if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2359     }
2360     if (!itsgold){
2361       //
2362       if (Float_t(shared+1)/Float_t(found+1)>factor){
2363         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2364         delete arr->RemoveAt(trackindex);
2365         continue;
2366       }      
2367       if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
2368         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2369         delete arr->RemoveAt(trackindex);
2370         continue;
2371       }
2372     }
2373
2374     good++;
2375     if (sharedfactor>0.4) continue;
2376     if (pt->GetKinkIndexes()[0]>0) continue;
2377     for (Int_t i=first; i<last; i++) {
2378       Int_t index=pt->GetClusterIndex2(i);
2379       // if (index<0 || index&0x8000 ) continue;
2380       if (index<0 || index&0x8000 ) continue;
2381       AliTPCclusterMI *c= pt->GetClusterPointer(i);        
2382       if (!c) continue;
2383       c->Use(10);  
2384     }    
2385   }
2386   fNtracks = good;
2387   if (fDebug>0){
2388     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2389   }
2390   delete []quality;
2391   delete []indexes;
2392 }
2393
2394 void AliTPCtrackerMI::UnsignClusters() 
2395 {
2396   //
2397   // loop over all clusters and unsign them
2398   //
2399   
2400   for (Int_t sec=0;sec<fkNIS;sec++){
2401     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2402       AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2403       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2404         //      if (cl[icl].IsUsed(10))         
2405         cl[icl].Use(-1);
2406       cl = fInnerSec[sec][row].GetClusters2();
2407       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2408         //if (cl[icl].IsUsed(10))       
2409           cl[icl].Use(-1);      
2410     }
2411   }
2412   
2413   for (Int_t sec=0;sec<fkNOS;sec++){
2414     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2415       AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2416       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2417         //if (cl[icl].IsUsed(10))       
2418           cl[icl].Use(-1);
2419       cl = fOuterSec[sec][row].GetClusters2();
2420       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2421         //if (cl[icl].IsUsed(10))       
2422         cl[icl].Use(-1);      
2423     }
2424   }
2425   
2426 }
2427
2428
2429
2430 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2431 {
2432   //
2433   //sign clusters to be "used"
2434   //
2435   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2436   // loop over "primaries"
2437   
2438   Float_t sumdens=0;
2439   Float_t sumdens2=0;
2440   Float_t sumn   =0;
2441   Float_t sumn2  =0;
2442   Float_t sumchi =0;
2443   Float_t sumchi2 =0;
2444
2445   Float_t sum    =0;
2446
2447   TStopwatch timer;
2448   timer.Start();
2449
2450   Int_t nseed = arr->GetEntriesFast();
2451   for (Int_t i=0; i<nseed; i++) {
2452     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2453     if (!pt) {
2454       continue;
2455     }    
2456     if (!(pt->IsActive())) continue;
2457     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2458     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2459       sumdens += dens;
2460       sumdens2+= dens*dens;
2461       sumn    += pt->GetNumberOfClusters();
2462       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2463       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2464       if (chi2>5) chi2=5;
2465       sumchi  +=chi2;
2466       sumchi2 +=chi2*chi2;
2467       sum++;
2468     }
2469   }
2470
2471   Float_t mdensity = 0.9;
2472   Float_t meann    = 130;
2473   Float_t meanchi  = 1;
2474   Float_t sdensity = 0.1;
2475   Float_t smeann    = 10;
2476   Float_t smeanchi  =0.4;
2477   
2478
2479   if (sum>20){
2480     mdensity = sumdens/sum;
2481     meann    = sumn/sum;
2482     meanchi  = sumchi/sum;
2483     //
2484     sdensity = sumdens2/sum-mdensity*mdensity;
2485     sdensity = TMath::Sqrt(sdensity);
2486     //
2487     smeann   = sumn2/sum-meann*meann;
2488     smeann   = TMath::Sqrt(smeann);
2489     //
2490     smeanchi = sumchi2/sum - meanchi*meanchi;
2491     smeanchi = TMath::Sqrt(smeanchi);
2492   }
2493
2494
2495   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2496   //
2497   for (Int_t i=0; i<nseed; i++) {
2498     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2499     if (!pt) {
2500       continue;
2501     }
2502     if (pt->GetBSigned()) continue;
2503     if (pt->GetBConstrain()) continue;    
2504     //if (!(pt->IsActive())) continue;
2505     /*
2506     Int_t found,foundable,shared;    
2507     pt->GetClusterStatistic(0,160,found, foundable,shared);
2508     if (shared/float(found)>0.3) {
2509       if (shared/float(found)>0.9 ){
2510         //delete arr->RemoveAt(i);
2511       }
2512       continue;
2513     }
2514     */
2515     Bool_t isok =kFALSE;
2516     if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2517       isok = kTRUE;
2518     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2519       isok =kTRUE;
2520     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2521       isok =kTRUE;
2522     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2523       isok =kTRUE;
2524     
2525     if (isok)     
2526       for (Int_t i=0; i<160; i++) {     
2527         Int_t index=pt->GetClusterIndex2(i);
2528         if (index<0) continue;
2529         AliTPCclusterMI *c= pt->GetClusterPointer(i);
2530         if (!c) continue;
2531         //if (!(c->IsUsed(10))) c->Use();  
2532         c->Use(10);  
2533       }
2534   }
2535   
2536   
2537   //
2538   Double_t maxchi  = meanchi+2.*smeanchi;
2539
2540   for (Int_t i=0; i<nseed; i++) {
2541     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2542     if (!pt) {
2543       continue;
2544     }    
2545     //if (!(pt->IsActive())) continue;
2546     if (pt->GetBSigned()) continue;
2547     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2548     if (chi>maxchi) continue;
2549
2550     Float_t bfactor=1;
2551     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2552    
2553     //sign only tracks with enoug big density at the beginning
2554     
2555     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue; 
2556     
2557     
2558     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2559     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2560    
2561     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2562     if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2563       minn=0;
2564
2565     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2566       //Int_t noc=pt->GetNumberOfClusters();
2567       pt->SetBSigned(kTRUE);
2568       for (Int_t i=0; i<160; i++) {
2569
2570         Int_t index=pt->GetClusterIndex2(i);
2571         if (index<0) continue;
2572         AliTPCclusterMI *c= pt->GetClusterPointer(i);
2573         if (!c) continue;
2574         //      if (!(c->IsUsed(10))) c->Use();  
2575         c->Use(10);  
2576       }
2577     }
2578   }
2579   //  gLastCheck = nseed;
2580   //  arr->Compress();
2581   if (fDebug>0){
2582     timer.Print();
2583   }
2584 }
2585
2586
2587 void  AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2588 {
2589   // stop not active tracks
2590   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2591   // take th2 as threshold for number of founded to number of foundable on last 20 active rows 
2592   Int_t nseed = arr->GetEntriesFast();  
2593   //
2594   for (Int_t i=0; i<nseed; i++) {
2595     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2596     if (!pt) {
2597       continue;
2598     }
2599     if (!(pt->IsActive())) continue;
2600     StopNotActive(pt,row0,th0, th1,th2);
2601   }
2602 }
2603
2604
2605
2606 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2607  Float_t th2) const
2608 {
2609   // stop not active tracks
2610   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2611   // take th2 as threshold for number of founded to number of foundable on last 20 active rows 
2612   Int_t sumgood1  = 0;
2613   Int_t sumgood2  = 0;
2614   Int_t foundable = 0;
2615   Int_t maxindex = seed->GetLastPoint();  //last foundable row
2616   if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2617     seed->Desactivate(10) ;
2618     return;
2619   }
2620
2621   for (Int_t i=row0; i<maxindex; i++){
2622     Int_t index = seed->GetClusterIndex2(i);
2623     if (index!=-1) foundable++;
2624     //if (!c) continue;
2625     if (foundable<=30) sumgood1++;
2626     if (foundable<=50) {
2627       sumgood2++;
2628     }
2629     else{ 
2630       break;
2631     }        
2632   }
2633   if (foundable>=30.){ 
2634      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2635   }
2636   if (foundable>=50)
2637     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2638 }
2639
2640
2641 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2642 {
2643   //
2644   // back propagation of ESD tracks
2645   //
2646   //return 0;
2647   fEvent = event;
2648   ReadSeeds(event,2);
2649   fIteration=2;
2650   //PrepareForProlongation(fSeeds,1);
2651   PropagateForward2(fSeeds);
2652   TObjArray arraySeed(fSeeds->GetEntries());
2653   for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2654     arraySeed.AddAt(fSeeds->At(i),i);    
2655   }
2656   SignShared(&arraySeed);
2657
2658   Int_t ntracks=0;
2659   Int_t nseed = fSeeds->GetEntriesFast();
2660   for (Int_t i=0;i<nseed;i++){
2661     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2662     if (!seed) continue;
2663     if (seed->GetKinkIndex(0)>0)  UpdateKinkQualityD(seed);  // update quality informations for kinks
2664
2665     seed->PropagateTo(fParam->GetInnerRadiusLow());
2666     seed->UpdatePoints();
2667     MakeBitmaps(seed);
2668     AliESDtrack *esd=event->GetTrack(i);
2669     seed->CookdEdx(0.02,0.6);
2670     CookLabel(seed,0.1); //For comparison only
2671     esd->SetTPCClusterMap(seed->GetClusterMap());
2672     esd->SetTPCSharedMap(seed->GetSharedMap());
2673     //
2674     if (AliTPCReconstructor::StreamLevel()>0 && seed!=0&&esd!=0) {
2675       TTreeSRedirector &cstream = *fDebugStreamer;
2676       cstream<<"Crefit"<<
2677         "Esd.="<<esd<<
2678         "Track.="<<seed<<
2679         "\n"; 
2680     }
2681
2682     if (seed->GetNumberOfClusters()>15){
2683       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit); 
2684       esd->SetTPCPoints(seed->GetPoints());
2685       esd->SetTPCPointsF(seed->GetNFoundable());
2686       Int_t ndedx   = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2687       Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2688       Float_t dedx  = seed->GetdEdx();
2689       esd->SetTPCsignal(dedx, sdedx, ndedx);
2690       //
2691       // add seed to the esd track in Calib level
2692       //
2693       if (AliTPCReconstructor::StreamLevel()>0){
2694         AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE); 
2695         esd->AddCalibObject(seedCopy);
2696       }
2697       ntracks++;
2698     }
2699     else{
2700       //printf("problem\n");
2701     }
2702   }
2703   //FindKinks(fSeeds,event);
2704   Info("RefitInward","Number of refitted tracks %d",ntracks);
2705   fEvent =0;
2706   //WriteTracks();
2707   return 0;
2708 }
2709
2710
2711 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2712 {
2713   //
2714   // back propagation of ESD tracks
2715   //
2716
2717   fEvent = event;
2718   fIteration = 1;
2719   ReadSeeds(event,1);
2720   PropagateBack(fSeeds); 
2721   RemoveUsed2(fSeeds,0.4,0.4,20);
2722   //
2723   Int_t nseed = fSeeds->GetEntriesFast();
2724   Int_t ntracks=0;
2725   for (Int_t i=0;i<nseed;i++){
2726     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2727     if (!seed) continue;
2728     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
2729     seed->UpdatePoints();
2730     AliESDtrack *esd=event->GetTrack(i);
2731     seed->CookdEdx(0.02,0.6);
2732     CookLabel(seed,0.1); //For comparison only
2733     if (seed->GetNumberOfClusters()>15){
2734       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2735       esd->SetTPCPoints(seed->GetPoints());
2736       esd->SetTPCPointsF(seed->GetNFoundable());
2737       Int_t ndedx   = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2738       Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2739       Float_t dedx  = seed->GetdEdx();
2740       esd->SetTPCsignal(dedx, sdedx, ndedx);
2741       ntracks++;
2742       Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2743       // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number      
2744       (*fDebugStreamer)<<"Cback"<<
2745         "Tr0.="<<seed<<
2746         "EventNrInFile="<<eventnumber<<
2747         "\n"; // patch 28 fev 06   
2748     }
2749   }
2750   //FindKinks(fSeeds,event);
2751   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2752   fEvent =0;
2753   //WriteTracks();
2754
2755   return 0;
2756 }
2757
2758
2759 void AliTPCtrackerMI::DeleteSeeds()
2760 {
2761   //
2762   //delete Seeds
2763
2764   Int_t nseed = fSeeds->GetEntriesFast();
2765   for (Int_t i=0;i<nseed;i++){
2766     AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2767     if (seed) delete fSeeds->RemoveAt(i);
2768   }
2769   delete fSeeds;
2770
2771   fSeeds =0;
2772 }
2773
2774 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2775 {
2776   //
2777   //read seeds from the event
2778   
2779   Int_t nentr=event->GetNumberOfTracks();
2780   if (fDebug>0){
2781     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2782   }
2783   if (fSeeds) 
2784     DeleteSeeds();
2785   if (!fSeeds){   
2786     fSeeds = new TObjArray(nentr);
2787   }
2788   UnsignClusters();
2789   //  Int_t ntrk=0;  
2790   for (Int_t i=0; i<nentr; i++) {
2791     AliESDtrack *esd=event->GetTrack(i);
2792     ULong_t status=esd->GetStatus();
2793     if (!(status&AliESDtrack::kTPCin)) continue;
2794     AliTPCtrack t(*esd);
2795     t.SetNumberOfClusters(0);
2796     //    AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2797     AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2798     seed->SetUniqueID(esd->GetID());
2799     for (Int_t ikink=0;ikink<3;ikink++) {
2800       Int_t index = esd->GetKinkIndex(ikink);
2801       seed->GetKinkIndexes()[ikink] = index;
2802       if (index==0) continue;
2803       index = TMath::Abs(index);
2804       AliESDkink * kink = fEvent->GetKink(index-1);
2805       if (kink&&esd->GetKinkIndex(ikink)<0){
2806         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2807         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,0);
2808       }
2809       if (kink&&esd->GetKinkIndex(ikink)>0){
2810         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2811         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,4);
2812       }
2813
2814     }
2815     if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.); 
2816     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2817     if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2818       fSeeds->AddAt(0,i);
2819       delete seed;
2820       continue;    
2821     }
2822     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 )  {
2823       Double_t par0[5],par1[5],alpha,x;
2824       esd->GetInnerExternalParameters(alpha,x,par0);
2825       esd->GetExternalParameters(x,par1);
2826       Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2827       Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2828       Double_t trdchi2=0;
2829       if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2830       //reset covariance if suspicious 
2831       if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2832         seed->ResetCovariance(10.);
2833     }
2834
2835     //
2836     //
2837     // rotate to the local coordinate system
2838     //   
2839     fSectors=fInnerSec; fN=fkNIS;    
2840     Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2841     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2842     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2843     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2844     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2845     if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2846     if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2847     alpha-=seed->GetAlpha();  
2848     if (!seed->Rotate(alpha)) {
2849       delete seed;
2850       continue;
2851     }
2852     seed->SetESD(esd);
2853     // sign clusters
2854     if (esd->GetKinkIndex(0)<=0){
2855       for (Int_t irow=0;irow<160;irow++){
2856         Int_t index = seed->GetClusterIndex2(irow);    
2857         if (index>0){ 
2858           //
2859           AliTPCclusterMI * cl = GetClusterMI(index);
2860           seed->SetClusterPointer(irow,cl);
2861           if (cl){
2862             if ((index & 0x8000)==0){
2863               cl->Use(10);  // accepted cluster   
2864             }else{
2865               cl->Use(6);   // close cluster not accepted
2866             }   
2867           }else{
2868             Info("ReadSeeds","Not found cluster");
2869           }
2870         }
2871       }
2872     }
2873     fSeeds->AddAt(seed,i);
2874   }
2875 }
2876
2877
2878
2879 //_____________________________________________________________________________
2880 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
2881                                  Float_t deltay, Int_t ddsec) {
2882   //-----------------------------------------------------------------
2883   // This function creates track seeds.
2884   // SEEDING WITH VERTEX CONSTRAIN 
2885   //-----------------------------------------------------------------
2886   // cuts[0]   - fP4 cut
2887   // cuts[1]   - tan(phi)  cut
2888   // cuts[2]   - zvertex cut
2889   // cuts[3]   - fP3 cut
2890   Int_t nin0  = 0;
2891   Int_t nin1  = 0;
2892   Int_t nin2  = 0;
2893   Int_t nin   = 0;
2894   Int_t nout1 = 0;
2895   Int_t nout2 = 0;
2896
2897   Double_t x[5], c[15];
2898   //  Int_t di = i1-i2;
2899   //
2900   AliTPCseed * seed = new AliTPCseed;
2901   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2902   Double_t cs=cos(alpha), sn=sin(alpha);
2903   //
2904   //  Double_t x1 =fOuterSec->GetX(i1);
2905   //Double_t xx2=fOuterSec->GetX(i2);
2906   
2907   Double_t x1 =GetXrow(i1);
2908   Double_t xx2=GetXrow(i2);
2909
2910   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2911
2912   Int_t imiddle = (i2+i1)/2;    //middle pad row index
2913   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2914   const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2915   //
2916   Int_t ns =sec;   
2917
2918   const AliTPCRow& kr1=GetRow(ns,i1);
2919   Double_t ymax  = GetMaxY(i1)-kr1.GetDeadZone()-1.5;  
2920   Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;  
2921
2922   //
2923   // change cut on curvature if it can't reach this layer
2924   // maximal curvature set to reach it
2925   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2926   if (dvertexmax*0.5*cuts[0]>0.85){
2927     cuts[0] = 0.85/(dvertexmax*0.5+1.);
2928   }
2929   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
2930
2931   //  Int_t ddsec = 1;
2932   if (deltay>0) ddsec = 0; 
2933   // loop over clusters  
2934   for (Int_t is=0; is < kr1; is++) {
2935     //
2936     if (kr1[is]->IsUsed(10)) continue;
2937     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
2938     //if (TMath::Abs(y1)>ymax) continue;
2939
2940     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2941
2942     // find possible directions    
2943     Float_t anglez = (z1-z3)/(x1-x3); 
2944     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z      
2945     //
2946     //
2947     //find   rotation angles relative to line given by vertex and point 1
2948     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2949     Double_t dvertex  = TMath::Sqrt(dvertex2);
2950     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
2951     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);            
2952     
2953     //
2954     // loop over 2 sectors
2955     Int_t dsec1=-ddsec;
2956     Int_t dsec2= ddsec;
2957     if (y1<0)  dsec2= 0;
2958     if (y1>0)  dsec1= 0;
2959     
2960     Double_t dddz1=0;  // direction of delta inclination in z axis
2961     Double_t dddz2=0;
2962     if ( (z1-z3)>0)
2963       dddz1 =1;    
2964     else
2965       dddz2 =1;
2966     //
2967     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2968       Int_t sec2 = sec + dsec;
2969       // 
2970       //      AliTPCRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2971       //AliTPCRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2972       AliTPCRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
2973       AliTPCRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2974       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2975       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2976
2977       // rotation angles to p1-p3
2978       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;            
2979       Double_t x2,   y2,   z2; 
2980       //
2981       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2982
2983       //
2984       Double_t dxx0 =  (xx2-x3)*cs13r;
2985       Double_t dyy0 =  (xx2-x3)*sn13r;
2986       for (Int_t js=index1; js < index2; js++) {
2987         const AliTPCclusterMI *kcl = kr2[js];
2988         if (kcl->IsUsed(10)) continue;  
2989         //
2990         //calcutate parameters
2991         //      
2992         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
2993         // stright track
2994         if (TMath::Abs(yy0)<0.000001) continue;
2995         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
2996         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2997         Double_t r02 = (0.25+y0*y0)*dvertex2;   
2998         //curvature (radius) cut
2999         if (r02<r2min) continue;                
3000        
3001         nin0++; 
3002         //
3003         Double_t c0  = 1/TMath::Sqrt(r02);
3004         if (yy0>0) c0*=-1.;     
3005                
3006        
3007         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
3008         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3009         Double_t dfi0   = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3010         Double_t dfi1   = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);  
3011         //
3012         //
3013         Double_t z0  =  kcl->GetZ();  
3014         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
3015         if (TMath::Abs(zzzz2-z0)>0.5) continue;       
3016         nin1++;              
3017         //      
3018         Double_t dip    = (z1-z0)*c0/dfi1;        
3019         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3020         //
3021         y2 = kcl->GetY(); 
3022         if (dsec==0){
3023           x2 = xx2; 
3024           z2 = kcl->GetZ();       
3025         }
3026         else
3027           {
3028             // rotation 
3029             z2 = kcl->GetZ();  
3030             x2= xx2*cs-y2*sn*dsec;
3031             y2=+xx2*sn*dsec+y2*cs;
3032           }
3033         
3034         x[0] = y1;
3035         x[1] = z1;
3036         x[2] = x0;
3037         x[3] = dip;
3038         x[4] = c0;
3039         //
3040         //
3041         // do we have cluster at the middle ?
3042         Double_t ym,zm;
3043         GetProlongation(x1,xm,x,ym,zm);
3044         UInt_t dummy; 
3045         AliTPCclusterMI * cm=0;
3046         if (TMath::Abs(ym)-ymaxm<0){      
3047           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3048           if ((!cm) || (cm->IsUsed(10))) {        
3049             continue;
3050           }
3051         }
3052         else{     
3053           // rotate y1 to system 0
3054           // get state vector in rotated system 
3055           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
3056           Double_t xr2  =  x0*cs+yr1*sn*dsec;
3057           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3058           //
3059           GetProlongation(xx2,xm,xr,ym,zm);
3060           if (TMath::Abs(ym)-ymaxm<0){
3061             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3062             if ((!cm) || (cm->IsUsed(10))) {      
3063               continue;
3064             }
3065           }
3066         }
3067        
3068
3069         Double_t dym = 0;
3070         Double_t dzm = 0;
3071         if (cm){
3072           dym = ym - cm->GetY();
3073           dzm = zm - cm->GetZ();
3074         }
3075         nin2++;
3076
3077
3078         //
3079         //
3080         Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3081         Double_t sy2=kcl->GetSigmaY2()*2.,     sz2=kcl->GetSigmaZ2()*2.;
3082         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3083         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3084         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3085
3086         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3087         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3088         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3089         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3090         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3091         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3092         
3093         Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3094         Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3095         Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3096         Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3097         
3098         c[0]=sy1;
3099         c[1]=0.;       c[2]=sz1;
3100         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3101         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3102                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3103         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3104         c[13]=f30*sy1*f40+f32*sy2*f42;
3105         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3106         
3107         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3108         
3109         UInt_t index=kr1.GetIndex(is);
3110         AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3111         
3112         track->SetIsSeeding(kTRUE);
3113         track->SetSeed1(i1);
3114         track->SetSeed2(i2);
3115         track->SetSeedType(3);
3116
3117        
3118         //if (dsec==0) {
3119           FollowProlongation(*track, (i1+i2)/2,1);
3120           Int_t foundable,found,shared;
3121           track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3122           if ((found<0.55*foundable)  || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3123             seed->Reset();
3124             seed->~AliTPCseed();
3125             continue;
3126           }
3127           //}
3128         
3129         nin++;
3130         FollowProlongation(*track, i2,1);
3131         
3132         
3133         //Int_t rc = 1;
3134         track->SetBConstrain(1);
3135         //      track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3136         track->SetLastPoint(i1);  // first cluster in track position
3137         track->SetFirstPoint(track->GetLastPoint());
3138         
3139         if (track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3140             track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || 
3141             track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3142           seed->Reset();
3143           seed->~AliTPCseed();
3144           continue;
3145         }
3146         nout1++;
3147         // Z VERTEX CONDITION
3148         Double_t zv, bz=GetBz();
3149         if ( !track->GetZAt(0.,bz,zv) ) continue;
3150         if (TMath::Abs(zv-z3)>cuts[2]) {
3151           FollowProlongation(*track, TMath::Max(i2-20,0));
3152           if ( !track->GetZAt(0.,bz,zv) ) continue;
3153           if (TMath::Abs(zv-z3)>cuts[2]){
3154             FollowProlongation(*track, TMath::Max(i2-40,0));
3155             if ( !track->GetZAt(0.,bz,zv) ) continue;
3156             if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3157               // make seed without constrain
3158               AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3159               FollowProlongation(*track2, i2,1);
3160               track2->SetBConstrain(kFALSE);
3161               track2->SetSeedType(1);
3162               arr->AddLast(track2); 
3163               seed->Reset();
3164               seed->~AliTPCseed();
3165               continue;         
3166             }
3167             else{
3168               seed->Reset();
3169               seed->~AliTPCseed();
3170               continue;
3171             
3172             }
3173           }
3174         }
3175         
3176         track->SetSeedType(0);
3177         arr->AddLast(track); 
3178         seed = new AliTPCseed;  
3179         nout2++;
3180         // don't consider other combinations
3181         if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3182           break;
3183       }
3184     }
3185   }
3186   if (fDebug>3){
3187     Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3188   }
3189   delete seed;
3190 }
3191
3192
3193 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3194                                  Float_t deltay) {
3195   
3196
3197
3198   //-----------------------------------------------------------------
3199   // This function creates track seeds.
3200   //-----------------------------------------------------------------
3201   // cuts[0]   - fP4 cut
3202   // cuts[1]   - tan(phi)  cut
3203   // cuts[2]   - zvertex cut
3204   // cuts[3]   - fP3 cut
3205
3206
3207   Int_t nin0  = 0;
3208   Int_t nin1  = 0;
3209   Int_t nin2  = 0;
3210   Int_t nin   = 0;
3211   Int_t nout1 = 0;
3212   Int_t nout2 = 0;
3213   Int_t nout3 =0;
3214   Double_t x[5], c[15];
3215   //
3216   // make temporary seed
3217   AliTPCseed * seed = new AliTPCseed;
3218   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3219   //  Double_t cs=cos(alpha), sn=sin(alpha);
3220   //
3221   //
3222
3223   // first 3 padrows
3224   Double_t x1 = GetXrow(i1-1);
3225   const    AliTPCRow& kr1=GetRow(sec,i1-1);
3226   Double_t y1max  = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;  
3227   //
3228   Double_t x1p = GetXrow(i1);
3229   const    AliTPCRow& kr1p=GetRow(sec,i1);
3230   //
3231   Double_t x1m = GetXrow(i1-2);
3232   const    AliTPCRow& kr1m=GetRow(sec,i1-2);
3233
3234   //
3235   //last 3 padrow for seeding
3236   AliTPCRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
3237   Double_t    x3   =  GetXrow(i1-7);
3238   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;  
3239   //
3240   AliTPCRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
3241   Double_t    x3p   = GetXrow(i1-6);
3242   //
3243   AliTPCRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
3244   Double_t    x3m   = GetXrow(i1-8);
3245
3246   //
3247   //
3248   // middle padrow
3249   Int_t im = i1-4;                           //middle pad row index
3250   Double_t xm         = GetXrow(im);         // radius of middle pad-row
3251   const AliTPCRow& krm=GetRow(sec,im);   //middle pad -row
3252   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;  
3253   //
3254   //
3255   Double_t deltax  = x1-x3;
3256   Double_t dymax   = deltax*cuts[1];
3257   Double_t dzmax   = deltax*cuts[3];
3258   //
3259   // loop over clusters  
3260   for (Int_t is=0; is < kr1; is++) {
3261     //
3262     if (kr1[is]->IsUsed(10)) continue;
3263     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
3264     //
3265     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge    
3266     // 
3267     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3268     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3269     //    
3270     Double_t y3,   z3;
3271     //
3272     //
3273     UInt_t index;
3274     for (Int_t js=index1; js < index2; js++) {
3275       const AliTPCclusterMI *kcl = kr3[js];
3276       if (kcl->IsUsed(10)) continue;
3277       y3 = kcl->GetY(); 
3278       // apply angular cuts
3279       if (TMath::Abs(y1-y3)>dymax) continue;
3280       x3 = x3; 
3281       z3 = kcl->GetZ(); 
3282       if (TMath::Abs(z1-z3)>dzmax) continue;
3283       //
3284       Double_t angley = (y1-y3)/(x1-x3);
3285       Double_t anglez = (z1-z3)/(x1-x3);
3286       //
3287       Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3288       Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3289       //
3290       Double_t yyym = angley*(xm-x1)+y1;
3291       Double_t zzzm = anglez*(xm-x1)+z1;
3292
3293       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3294       if (!kcm) continue;
3295       if (kcm->IsUsed(10)) continue;
3296       
3297       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3298       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3299       //
3300       //
3301       //
3302       Int_t used  =0;
3303       Int_t found =0;
3304       //
3305       // look around first
3306       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3307                                                       anglez*(x1m-x1)+z1,
3308                                                       erry,errz,index);
3309       //
3310       if (kc1m){
3311         found++;
3312         if (kc1m->IsUsed(10)) used++;
3313       }
3314       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3315                                                       anglez*(x1p-x1)+z1,
3316                                                       erry,errz,index);
3317       //
3318       if (kc1p){
3319         found++;
3320         if (kc1p->IsUsed(10)) used++;
3321       }
3322       if (used>1)  continue;
3323       if (found<1) continue; 
3324
3325       //
3326       // look around last
3327       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3328                                                       anglez*(x3m-x3)+z3,
3329                                                       erry,errz,index);
3330       //
3331       if (kc3m){
3332         found++;
3333         if (kc3m->IsUsed(10)) used++;
3334       }
3335       else 
3336         continue;
3337       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3338                                                       anglez*(x3p-x3)+z3,
3339                                                       erry,errz,index);
3340       //
3341       if (kc3p){
3342         found++;
3343         if (kc3p->IsUsed(10)) used++;
3344       }
3345       else 
3346         continue;
3347       if (used>1)  continue;
3348       if (found<3) continue;       
3349       //
3350       Double_t x2,y2,z2;
3351       x2 = xm;
3352       y2 = kcm->GetY();
3353       z2 = kcm->GetZ();
3354       //
3355                         
3356       x[0]=y1;
3357       x[1]=z1;
3358       x[4]=F1(x1,y1,x2,y2,x3,y3);
3359       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3360       nin0++;
3361       //
3362       x[2]=F2(x1,y1,x2,y2,x3,y3);
3363       nin1++;
3364       //
3365       x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3366       //if (TMath::Abs(x[3]) > cuts[3]) continue;
3367       nin2++;
3368       //
3369       //
3370       Double_t sy1=0.1,  sz1=0.1;
3371       Double_t sy2=0.1,  sz2=0.1;
3372       Double_t sy3=0.1,  sy=0.1, sz=0.1;
3373       
3374       Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3375       Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3376       Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3377       Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3378       Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3379       Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3380       
3381       Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3382       Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3383       Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3384       Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3385       
3386       c[0]=sy1;
3387       c[1]=0.;       c[2]=sz1; 
3388       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3389       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3390       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3391       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3392       c[13]=f30*sy1*f40+f32*sy2*f42;
3393       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3394       
3395       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3396       
3397       UInt_t index=kr1.GetIndex(is);
3398       AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3399       
3400       track->SetIsSeeding(kTRUE);
3401
3402       nin++;      
3403       FollowProlongation(*track, i1-7,1);
3404       if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 || 
3405           track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3406         seed->Reset();
3407         seed->~AliTPCseed();
3408         continue;
3409       }
3410       nout1++;
3411       nout2++;  
3412       //Int_t rc = 1;
3413       FollowProlongation(*track, i2,1);
3414       track->SetBConstrain(0);
3415       track->SetLastPoint(i1+fInnerSec->GetNRows());  // first cluster in track position
3416       track->SetFirstPoint(track->GetLastPoint());
3417       
3418       if (track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3419           track->GetNumberOfClusters()<track->GetNFoundable()*0.7 || 
3420           track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3421         seed->Reset();
3422         seed->~AliTPCseed();
3423         continue;
3424       }
3425    
3426       {
3427         FollowProlongation(*track, TMath::Max(i2-10,0),1);
3428         AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3429         FollowProlongation(*track2, i2,1);
3430         track2->SetBConstrain(kFALSE);
3431         track2->SetSeedType(4);
3432         arr->AddLast(track2); 
3433         seed->Reset();
3434         seed->~AliTPCseed();
3435       }
3436       
3437    
3438       //arr->AddLast(track); 
3439       //seed = new AliTPCseed;  
3440       nout3++;
3441     }
3442   }
3443   
3444   if (fDebug>3){
3445     Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3446   }
3447   delete seed;
3448 }
3449
3450
3451 //_____________________________________________________________________________
3452 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3453                                  Float_t deltay, Bool_t /*bconstrain*/) {
3454   //-----------------------------------------------------------------
3455   // This function creates track seeds - without vertex constraint
3456   //-----------------------------------------------------------------
3457   // cuts[0]   - fP4 cut        - not applied
3458   // cuts[1]   - tan(phi)  cut
3459   // cuts[2]   - zvertex cut    - not applied 
3460   // cuts[3]   - fP3 cut
3461   Int_t nin0=0;
3462   Int_t nin1=0;
3463   Int_t nin2=0;
3464   Int_t nin3=0;
3465   //  Int_t nin4=0;
3466   //Int_t nin5=0;
3467
3468   
3469
3470   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3471   //  Double_t cs=cos(alpha), sn=sin(alpha);
3472   Int_t row0 = (i1+i2)/2;
3473   Int_t drow = (i1-i2)/2;
3474   const AliTPCRow& kr0=fSectors[sec][row0];
3475   AliTPCRow * kr=0;
3476
3477   AliTPCpolyTrack polytrack;
3478   Int_t nclusters=fSectors[sec][row0];
3479   AliTPCseed * seed = new AliTPCseed;
3480
3481   Int_t sumused=0;
3482   Int_t cused=0;
3483   Int_t cnused=0;
3484   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3485     Int_t nfound =0;
3486     Int_t nfoundable =0;
3487     for (Int_t iter =1; iter<2; iter++){   //iterations
3488       const AliTPCRow& krm=fSectors[sec][row0-iter];
3489       const AliTPCRow& krp=fSectors[sec][row0+iter];      
3490       const AliTPCclusterMI * cl= kr0[is];
3491       
3492       if (cl->IsUsed(10)) {
3493         cused++;
3494       }
3495       else{
3496         cnused++;
3497       }
3498       Double_t x = kr0.GetX();
3499       // Initialization of the polytrack
3500       nfound =0;
3501       nfoundable =0;
3502       polytrack.Reset();
3503       //
3504       Double_t y0= cl->GetY();
3505       Double_t z0= cl->GetZ();
3506       Float_t erry = 0;
3507       Float_t errz = 0;
3508       
3509       Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3510       if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;  // seed only at the edge
3511       
3512       erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;      
3513       errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;      
3514       polytrack.AddPoint(x,y0,z0,erry, errz);
3515
3516       sumused=0;
3517       if (cl->IsUsed(10)) sumused++;
3518
3519
3520       Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3521       Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3522       //
3523       x = krm.GetX();
3524       AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3525       if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3526         erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;          
3527         errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3528         if (cl1->IsUsed(10))  sumused++;
3529         polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3530       }
3531       //
3532       x = krp.GetX();
3533       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3534       if (cl2) {
3535         erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;          
3536         errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3537         if (cl2->IsUsed(10)) sumused++;  
3538         polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3539       }
3540       //
3541       if (sumused>0) continue;
3542       nin0++;
3543       polytrack.UpdateParameters();
3544       // follow polytrack
3545       roadz = 1.2;
3546       roady = 1.2;
3547       //
3548       Double_t yn,zn;
3549       nfoundable = polytrack.GetN();
3550       nfound     = nfoundable; 
3551       //
3552       for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3553         Float_t maxdist = 0.8*(1.+3./(ddrow));
3554         for (Int_t delta = -1;delta<=1;delta+=2){
3555           Int_t row = row0+ddrow*delta;
3556           kr = &(fSectors[sec][row]);
3557           Double_t xn = kr->GetX();
3558           Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3559           polytrack.GetFitPoint(xn,yn,zn);
3560           if (TMath::Abs(yn)>ymax) continue;
3561           nfoundable++;
3562           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3563           if (cln) {
3564             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3565             if (dist<maxdist){
3566               /*
3567               erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));         
3568               errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3569               if (cln->IsUsed(10)) {
3570                 //      printf("used\n");
3571                 sumused++;
3572                 erry*=2;
3573                 errz*=2;
3574               }
3575               */
3576               erry=0.1;
3577               errz=0.1;
3578               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3579               nfound++;
3580             }
3581           }
3582         }
3583         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;     
3584         polytrack.UpdateParameters();
3585       }           
3586     }
3587     if ( (sumused>3) || (sumused>0.5*nfound))  {
3588       //printf("sumused   %d\n",sumused);
3589       continue;
3590     }
3591     nin1++;
3592     Double_t dy,dz;
3593     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3594     AliTPCpolyTrack track2;
3595     
3596     polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3597     if (track2.GetN()<0.5*nfoundable) continue;
3598     nin2++;
3599
3600     if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3601       //
3602       // test seed with and without constrain
3603       for (Int_t constrain=0; constrain<=0;constrain++){
3604         // add polytrack candidate
3605
3606         Double_t x[5], c[15];
3607         Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3608         track2.GetBoundaries(x3,x1);    
3609         x2 = (x1+x3)/2.;
3610         track2.GetFitPoint(x1,y1,z1);
3611         track2.GetFitPoint(x2,y2,z2);
3612         track2.GetFitPoint(x3,y3,z3);
3613         //
3614         //is track pointing to the vertex ?
3615         Double_t x0,y0,z0;
3616         x0=0;
3617         polytrack.GetFitPoint(x0,y0,z0);
3618
3619         if (constrain) {
3620           x2 = x3;
3621           y2 = y3;
3622           z2 = z3;
3623           
3624           x3 = 0;
3625           y3 = 0;
3626           z3 = 0;
3627         }
3628         x[0]=y1;
3629         x[1]=z1;
3630         x[4]=F1(x1,y1,x2,y2,x3,y3);
3631                 
3632         //      if (TMath::Abs(x[4]) >= cuts[0]) continue;  //
3633         x[2]=F2(x1,y1,x2,y2,x3,y3);
3634         
3635         //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3636         //x[3]=F3(x1,y1,x2,y2,z1,z2);
3637         x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3638         //if (TMath::Abs(x[3]) > cuts[3]) continue;
3639
3640         
3641         Double_t sy =0.1, sz =0.1;
3642         Double_t sy1=0.02, sz1=0.02;
3643         Double_t sy2=0.02, sz2=0.02;
3644         Double_t sy3=0.02;
3645
3646         if (constrain){
3647           sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3648         }
3649         
3650         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3651         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3652         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3653         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3654         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3655         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3656
3657         Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3658         Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3659         Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3660         Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3661
3662         
3663         c[0]=sy1;
3664         c[1]=0.;       c[2]=sz1;
3665         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3666         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3667         c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3668         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3669         c[13]=f30*sy1*f40+f32*sy2*f42;
3670         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3671         
3672         //Int_t row1 = fSectors->GetRowNumber(x1);
3673         Int_t row1 = GetRowNumber(x1);
3674
3675         UInt_t index=0;
3676         //kr0.GetIndex(is);
3677         AliTPCseed *track=new (seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3678         track->SetIsSeeding(kTRUE);
3679         Int_t rc=FollowProlongation(*track, i2);        
3680         if (constrain) track->SetBConstrain(1);
3681         else
3682           track->SetBConstrain(0);
3683         track->SetLastPoint(row1+fInnerSec->GetNRows());  // first cluster in track position
3684         track->SetFirstPoint(track->GetLastPoint());
3685
3686         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3687             track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || 
3688             track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3689           //delete track;
3690           seed->Reset();
3691           seed->~AliTPCseed();
3692         }
3693         else {
3694           arr->AddLast(track);
3695           seed = new AliTPCseed;
3696         }
3697         nin3++;
3698       }
3699     }  // if accepted seed
3700   }
3701   if (fDebug>3){
3702     Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3703   }
3704   delete seed;
3705 }
3706
3707
3708 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3709 {
3710   //
3711   //
3712   //reseed using track points
3713   Int_t p0 = int(r0*track->GetNumberOfClusters());     // point 0 
3714   Int_t p1 = int(r1*track->GetNumberOfClusters());
3715   Int_t p2 = int(r2*track->GetNumberOfClusters());   // last point
3716   Int_t pp2=0;
3717   Double_t  x0[3],x1[3],x2[3];
3718   for (Int_t i=0;i<3;i++){
3719     x0[i]=-1;
3720     x1[i]=-1;
3721     x2[i]=-1;
3722   }
3723
3724   // find track position at given ratio of the length
3725   Int_t  sec0=0, sec1=0, sec2=0;
3726   Int_t index=-1;
3727   Int_t clindex;
3728   for (Int_t i=0;i<160;i++){
3729     if (track->GetClusterPointer(i)){
3730       index++;
3731       AliTPCTrackerPoint   *trpoint =track->GetTrackPoint(i);
3732       if ( (index<p0) || x0[0]<0 ){
3733         if (trpoint->GetX()>1){
3734           clindex = track->GetClusterIndex2(i);
3735           if (clindex>0){       
3736             x0[0] = trpoint->GetX();
3737             x0[1] = trpoint->GetY();
3738             x0[2] = trpoint->GetZ();
3739             sec0  = ((clindex&0xff000000)>>24)%18;
3740           }
3741         }
3742       }
3743
3744       if ( (index<p1) &&(trpoint->GetX()>1)){
3745         clindex = track->GetClusterIndex2(i);
3746         if (clindex>0){
3747           x1[0] = trpoint->GetX();
3748           x1[1] = trpoint->GetY();
3749           x1[2] = trpoint->GetZ();
3750           sec1  = ((clindex&0xff000000)>>24)%18;
3751         }
3752       }
3753       if ( (index<p2) &&(trpoint->GetX()>1)){
3754         clindex = track->GetClusterIndex2(i);
3755         if (clindex>0){
3756           x2[0] = trpoint->GetX();
3757           x2[1] = trpoint->GetY();
3758           x2[2] = trpoint->GetZ(); 
3759           sec2  = ((clindex&0xff000000)>>24)%18;
3760           pp2 = i;
3761         }
3762       }
3763     }
3764   }
3765   
3766   Double_t alpha, cs,sn, xx2,yy2;
3767   //
3768   alpha = (sec1-sec2)*fSectors->GetAlpha();
3769   cs = TMath::Cos(alpha);
3770   sn = TMath::Sin(alpha); 
3771   xx2= x1[0]*cs-x1[1]*sn;
3772   yy2= x1[0]*sn+x1[1]*cs;
3773   x1[0] = xx2;
3774   x1[1] = yy2;
3775   //
3776   alpha = (sec0-sec2)*fSectors->GetAlpha();
3777   cs = TMath::Cos(alpha);
3778   sn = TMath::Sin(alpha); 
3779   xx2= x0[0]*cs-x0[1]*sn;
3780   yy2= x0[0]*sn+x0[1]*cs;
3781   x0[0] = xx2;
3782   x0[1] = yy2;
3783   //
3784   //
3785   //
3786   Double_t x[5],c[15];
3787   //
3788   x[0]=x2[1];
3789   x[1]=x2[2];
3790   x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3791   //  if (x[4]>1) return 0;
3792   x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3793   x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3794   //if (TMath::Abs(x[3]) > 2.2)  return 0;
3795   //if (TMath::Abs(x[2]) > 1.99) return 0;
3796   //  
3797   Double_t sy =0.1,  sz =0.1;
3798   //
3799   Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3800   Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3801   Double_t sy3=0.01+track->GetSigmaY2();
3802   //
3803   Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3804   Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3805   Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3806   Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3807   Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3808   Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3809   //
3810   Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3811   Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3812   Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3813   Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3814   
3815   
3816   c[0]=sy1;
3817   c[1]=0.;       c[2]=sz1;
3818   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3819   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3820   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3821   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3822   c[13]=f30*sy1*f40+f32*sy2*f42;
3823   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3824   
3825   //  Int_t row1 = fSectors->GetRowNumber(x2[0]);
3826   AliTPCseed *seed=new  AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3827   //  Double_t y0,z0,y1,z1, y2,z2;
3828   //seed->GetProlongation(x0[0],y0,z0);
3829   // seed->GetProlongation(x1[0],y1,z1);
3830   //seed->GetProlongation(x2[0],y2,z2);
3831   //  seed =0;
3832   seed->SetLastPoint(pp2);
3833   seed->SetFirstPoint(pp2);
3834   
3835
3836   return seed;
3837 }
3838
3839
3840 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3841 {
3842   //
3843   //
3844   //reseed using founded clusters 
3845   //
3846   // Find the number of clusters
3847   Int_t nclusters = 0;
3848   for (Int_t irow=0;irow<160;irow++){
3849     if (track->GetClusterIndex(irow)>0) nclusters++;
3850   }
3851   //
3852   Int_t ipos[3];
3853   ipos[0] = TMath::Max(int(r0*nclusters),0);             // point 0 cluster
3854   ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1);   // 
3855   ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1);   // last point
3856   //
3857   //
3858   Double_t  xyz[3][3];
3859   Int_t     row[3],sec[3]={0,0,0};
3860   //
3861   // find track row position at given ratio of the length
3862   Int_t index=-1;
3863   for (Int_t irow=0;irow<160;irow++){    
3864     if (track->GetClusterIndex2(irow)<0) continue;
3865     index++;
3866     for (Int_t ipoint=0;ipoint<3;ipoint++){
3867       if (index<=ipos[ipoint]) row[ipoint] = irow;
3868     }        
3869   }
3870   //
3871   //Get cluster and sector position
3872   for (Int_t ipoint=0;ipoint<3;ipoint++){
3873     Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3874     AliTPCclusterMI * cl = GetClusterMI(clindex);
3875     if (cl==0) {
3876       //Error("Bug\n");
3877       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
3878       return 0;
3879     }
3880     sec[ipoint]     = ((clindex&0xff000000)>>24)%18;
3881     xyz[ipoint][0]  = GetXrow(row[ipoint]);
3882     xyz[ipoint][1]  = cl->GetY();
3883     xyz[ipoint][2]  = cl->GetZ();
3884   }
3885   //
3886   //
3887   // Calculate seed state vector and covariance matrix
3888
3889   Double_t alpha, cs,sn, xx2,yy2;
3890   //
3891   alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3892   cs = TMath::Cos(alpha);
3893   sn = TMath::Sin(alpha); 
3894   xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3895   yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3896   xyz[1][0] = xx2;
3897   xyz[1][1] = yy2;
3898   //
3899   alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3900   cs = TMath::Cos(alpha);
3901   sn = TMath::Sin(alpha); 
3902   xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3903   yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3904   xyz[0][0] = xx2;
3905   xyz[0][1] = yy2;
3906   //
3907   //
3908   //
3909   Double_t x[5],c[15];
3910   //
3911   x[0]=xyz[2][1];
3912   x[1]=xyz[2][2];
3913   x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3914   x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3915   x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3916   //  
3917   Double_t sy =0.1,  sz =0.1;
3918   //
3919   Double_t sy1=0.2, sz1=0.2;
3920   Double_t sy2=0.2, sz2=0.2;
3921   Double_t sy3=0.2;
3922   //
3923   Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
3924   Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
3925   Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
3926   Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
3927   Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
3928   Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
3929   //
3930   Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
3931   Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
3932   Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
3933   Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
3934   
3935   
3936   c[0]=sy1;
3937   c[1]=0.;       c[2]=sz1;
3938   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3939   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3940   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3941   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3942   c[13]=f30*sy1*f40+f32*sy2*f42;
3943   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3944   
3945   //  Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3946   AliTPCseed *seed=new  AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3947   seed->SetLastPoint(row[2]);
3948   seed->SetFirstPoint(row[2]);  
3949   return seed;
3950 }
3951
3952
3953 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3954 {
3955   //
3956   //
3957   //reseed using founded clusters 
3958   //
3959   Double_t  xyz[3][3];
3960   Int_t     row[3]={0,0,0};
3961   Int_t     sec[3]={0,0,0};
3962   //
3963   // forward direction
3964   if (forward){
3965     for (Int_t irow=r0;irow<160;irow++){
3966       if (track->GetClusterIndex(irow)>0){
3967         row[0] = irow;
3968         break;
3969       }
3970     }
3971     for (Int_t irow=160;irow>r0;irow--){
3972       if (track->GetClusterIndex(irow)>0){
3973         row[2] = irow;
3974         break;
3975       }
3976     }
3977     for (Int_t irow=row[2]-15;irow>row[0];irow--){
3978       if (track->GetClusterIndex(irow)>0){
3979         row[1] = irow;
3980         break;
3981       }
3982     }
3983     //
3984   }
3985   if (!forward){
3986     for (Int_t irow=0;irow<r0;irow++){
3987       if (track->GetClusterIndex(irow)>0){
3988         row[0] = irow;
3989         break;
3990       }
3991     }
3992     for (Int_t irow=r0;irow>0;irow--){
3993       if (track->GetClusterIndex(irow)>0){
3994         row[2] = irow;
3995         break;
3996       }
3997     }    
3998     for (Int_t irow=row[2]-15;irow>row[0];irow--){
3999       if (track->GetClusterIndex(irow)>0){
4000         row[1] = irow;
4001         break;
4002       }
4003     } 
4004   }
4005   //
4006   if ((row[2]-row[0])<20) return 0;
4007   if (row[1]==0) return 0;
4008   //
4009   //
4010   //Get cluster and sector position
4011   for (Int_t ipoint=0;ipoint<3;ipoint++){
4012     Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4013     AliTPCclusterMI * cl = GetClusterMI(clindex);
4014     if (cl==0) {
4015       //Error("Bug\n");
4016       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
4017       return 0;
4018     }
4019     sec[ipoint]     = ((clindex&0xff000000)>>24)%18;
4020     xyz[ipoint][0]  = GetXrow(row[ipoint]);
4021     AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);    
4022     if (point&&ipoint<2){
4023       //
4024        xyz[ipoint][1]  = point->GetY();
4025        xyz[ipoint][2]  = point->GetZ();
4026     }
4027     else{
4028       xyz[ipoint][1]  = cl->GetY();
4029       xyz[ipoint][2]  = cl->GetZ();
4030     }
4031   }
4032   //
4033   //
4034   //
4035   //
4036   // Calculate seed state vector and covariance matrix
4037
4038   Double_t alpha, cs,sn, xx2,yy2;
4039   //
4040   alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4041   cs = TMath::Cos(alpha);
4042   sn = TMath::Sin(alpha); 
4043   xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4044   yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4045   xyz[1][0] = xx2;
4046   xyz[1][1] = yy2;
4047   //
4048   alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4049   cs = TMath::Cos(alpha);
4050   sn = TMath::Sin(alpha); 
4051   xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4052   yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4053   xyz[0][0] = xx2;
4054   xyz[0][1] = yy2;
4055   //
4056   //
4057   //
4058   Double_t x[5],c[15];
4059   //
4060   x[0]=xyz[2][1];
4061   x[1]=xyz[2][2];
4062   x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4063   x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4064   x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4065   //  
4066   Double_t sy =0.1,  sz =0.1;
4067   //
4068   Double_t sy1=0.2, sz1=0.2;
4069   Double_t sy2=0.2, sz2=0.2;
4070   Double_t sy3=0.2;
4071   //
4072   Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
4073   Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
4074   Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
4075   Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
4076   Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
4077   Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
4078   //
4079   Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
4080   Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
4081   Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
4082   Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
4083   
4084   
4085   c[0]=sy1;
4086   c[1]=0.;       c[2]=sz1;
4087   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4088   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
4089   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4090   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4091   c[13]=f30*sy1*f40+f32*sy2*f42;
4092   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4093   
4094   //  Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4095   AliTPCseed *seed=new  AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4096   seed->SetLastPoint(row[2]);
4097   seed->SetFirstPoint(row[2]);  
4098   for (Int_t i=row[0];i<row[2];i++){
4099     seed->SetClusterIndex(i, track->GetClusterIndex(i));
4100   }
4101
4102   return seed;
4103 }
4104
4105 void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4106 {
4107   //
4108   //  find kinks
4109   //
4110   //
4111
4112   TObjArray *kinks= new TObjArray(10000);
4113   //  TObjArray *v0s= new TObjArray(10000);
4114   Int_t nentries = array->GetEntriesFast();
4115   AliHelix *helixes      = new AliHelix[nentries];
4116   Int_t    *sign         = new Int_t[nentries];
4117   Int_t    *nclusters    = new Int_t[nentries];
4118   Float_t  *alpha        = new Float_t[nentries];
4119   AliKink  *kink         = new AliKink();
4120   Int_t      * usage     = new Int_t[nentries];
4121   Float_t  *zm           = new Float_t[nentries];
4122   Float_t  *z0           = new Float_t[nentries]; 
4123   Float_t  *fim          = new Float_t[nentries];
4124   Float_t  *shared       = new Float_t[nentries];
4125   Bool_t   *circular     = new Bool_t[nentries];
4126   Float_t *dca          = new Float_t[nentries];
4127   //const AliESDVertex * primvertex = esd->GetVertex();
4128   //
4129   //  nentries = array->GetEntriesFast();
4130   //
4131   
4132   //
4133   //
4134   for (Int_t i=0;i<nentries;i++){
4135     sign[i]=0;
4136     usage[i]=0;
4137     AliTPCseed* track = (AliTPCseed*)array->At(i);    
4138     if (!track) continue;
4139     track->SetCircular(0);
4140     shared[i] = kFALSE;
4141     track->UpdatePoints();
4142     if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4143     }
4144     nclusters[i]=track->GetNumberOfClusters();
4145     alpha[i] = track->GetAlpha();
4146     new (&helixes[i]) AliHelix(*track);
4147     Double_t xyz[3];
4148     helixes[i].Evaluate(0,xyz);
4149     sign[i] = (track->GetC()>0) ? -1:1;
4150     Double_t x,y,z;
4151     x=160;
4152     if (track->GetProlongation(x,y,z)){
4153       zm[i]  = z;
4154       fim[i] = alpha[i]+TMath::ATan2(y,x);
4155     }
4156     else{
4157       zm[i]  = track->GetZ();
4158       fim[i] = alpha[i];
4159     }   
4160     z0[i]=1000;
4161     circular[i]= kFALSE;
4162     if (track->GetProlongation(0,y,z))  z0[i] = z;
4163     dca[i] = track->GetD(0,0);    
4164   }
4165   //
4166   //
4167   TStopwatch timer;
4168   timer.Start();
4169   Int_t ncandidates =0;
4170   Int_t nall =0;
4171   Int_t ntracks=0; 
4172   Double_t phase[2][2],radius[2];
4173
4174   //
4175   // Find circling track
4176   TTreeSRedirector &cstream = *fDebugStreamer;
4177   //
4178   for (Int_t i0=0;i0<nentries;i0++){
4179     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4180     if (!track0) continue;    
4181     if (track0->GetNumberOfClusters()<40) continue;
4182     if (TMath::Abs(1./track0->GetC())>200) continue;
4183     for (Int_t i1=i0+1;i1<nentries;i1++){
4184       AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4185       if (!track1) continue;
4186       if (track1->GetNumberOfClusters()<40)                  continue;
4187       if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4188       if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4189       if (TMath::Abs(1./track1->GetC())>200) continue;
4190       if (track1->Get1Pt()*track0->Get1Pt()>0)      continue;
4191       if (track1->GetTgl()*track0->GetTgl()>0)      continue;
4192       if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4193       if (track0->GetBConstrain()&&TMath::Abs(track1->Get1Pt())<TMath::Abs(track0->Get1Pt())) continue; //returning - lower momenta
4194       if (track1->GetBConstrain()&&TMath::Abs(track0->Get1Pt())<TMath::Abs(track1->Get1Pt())) continue; //returning - lower momenta
4195       //
4196       Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4197       if (mindcar<5)   continue;
4198       Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4199       if (mindcaz<5) continue;
4200       if (mindcar+mindcaz<20) continue;
4201       //
4202       //
4203       Float_t xc0 = helixes[i0].GetHelix(6);
4204       Float_t yc0 = helixes[i0].GetHelix(7);
4205       Float_t r0  = helixes[i0].GetHelix(8);
4206       Float_t xc1 = helixes[i1].GetHelix(6);
4207       Float_t yc1 = helixes[i1].GetHelix(7);
4208       Float_t r1  = helixes[i1].GetHelix(8);
4209         
4210       Float_t rmean = (r0+r1)*0.5;
4211       Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4212       //if (delta>30) continue;
4213       if (delta>rmean*0.25) continue;
4214       if (TMath::Abs(r0-r1)/rmean>0.3) continue; 
4215       //
4216       Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4217       if (npoints==0) continue;
4218       helixes[i0].GetClosestPhases(helixes[i1], phase);
4219       //
4220       Double_t xyz0[3];
4221       Double_t xyz1[3];
4222       Double_t hangles[3];
4223       helixes[i0].Evaluate(phase[0][0],xyz0);
4224       helixes[i1].Evaluate(phase[0][1],xyz1);
4225
4226       helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4227       Double_t deltah[2],deltabest;
4228       if (hangles[2]<2.8) continue;
4229       /*
4230       cstream<<"C"<<track0->fLab<<track1->fLab<<
4231         track0->fP3<<track1->fP3<<
4232         track0->fP4<<track1->fP4<<
4233         delta<<rmean<<npoints<<
4234         hangles[0]<<hangles[2]<<
4235         xyz0[2]<<xyz1[2]<<radius[0]<<"\n"; 
4236       */
4237       if (npoints>0){
4238         Int_t ibest=0;
4239         helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4240         if (npoints==2){
4241           helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4242           if (deltah[1]<deltah[0]) ibest=1;
4243         }
4244         deltabest  = TMath::Sqrt(deltah[ibest]);
4245         helixes[i0].Evaluate(phase[ibest][0],xyz0);
4246         helixes[i1].Evaluate(phase[ibest][1],xyz1);