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