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