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