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