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