770812ca37261d16665fc75a3b553f5d3f890a0c
[u/mrichter/AliRoot.git] / TPC / TPCrec / 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   // Calculate crosstalk estimate
1414   //
1415   const Int_t nROCs   = 72;
1416   const Int_t   nIterations=3;  // 
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     //
1494     // 2.) dump the crosstalk matrices to tree for further investigation
1495     //     a.) to estimate fluctuation of pedestal in indiviula wire segments
1496     //     b.) to check correlation between regions
1497     //     c.) to check relative conribution of signal below threshold to crosstalk
1498     
1499     if (AliTPCReconstructor::StreamLevel()&kStreamCrosstalkMatrix) {
1500       for (Int_t isector=0; isector<nROCs; isector++){  //set all ellemts of crosstalk matrix to 0
1501         TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1502         TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
1503         TMatrixD * crossTalkMatrixCache = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs*2);
1504         TVectorD vecAll(crossTalkMatrix->GetNrows());
1505         TVectorD vecBelow(crossTalkMatrix->GetNrows());
1506         TVectorD vecCache(crossTalkMatrixCache->GetNrows());
1507         //
1508         for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
1509           for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
1510             vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
1511             vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
1512             vecCache[iwire]=(*crossTalkMatrixCache)(iwire,itime);
1513           }
1514           (*fDebugStreamer)<<"crosstalkMatrix"<<
1515             "iter="<<iter<<                      //iteration
1516             "sector="<<isector<<                 // sector
1517             "itime="<<itime<<                    // time bin index
1518             "vecAll.="<<&vecAll<<                // crosstalk charge + charge below threshold
1519             "vecCache.="<<&vecCache<<                // crosstalk charge + charge below threshold         
1520             "vecBelow.="<<&vecBelow<<            // crosstalk contribution from signal below threshold
1521             "\n";
1522         }
1523       }
1524     }
1525     if (iter<nIterations-1) for (Int_t isector=0; isector<nROCs*2; isector++){  //set all ellemts of crosstalk matrix to 0 
1526       TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1527       TMatrixD * crossTalkMatrixCache = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs*2);
1528       if (crossTalkMatrix){
1529         (*crossTalkMatrixCache)*=0;
1530         (*crossTalkMatrixCache)+=(*crossTalkMatrix);
1531         (*crossTalkMatrix)*=0;
1532       }
1533     }      
1534   }
1535
1536
1537 }
1538
1539
1540
1541
1542 void    AliTPCtracker::FilterOutlierClusters(){
1543   //
1544   // filter outlier clusters  
1545   //
1546   /*
1547     1.)..... booking part
1548     nSectors=72;
1549     nTimeBins=fParam->Get....
1550     TH2F hisTime("","", sector,0,sector, nTimeBins,0,nTimeBins);
1551     TH2F hisPadRow("","", sector,0,sector, nPadRows,0,nPadRows);
1552     2.) .... filling part
1553     .... cluster loop { hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin()); }
1554     
1555     3.) ...filtering part 
1556     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... }
1557     
1558     sector loop
1559     { 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 } 
1560     //
1561     4. Disabling clusters
1562     
1563   */
1564   
1565   //
1566   // 1.) booking part 
1567   // 
1568   //  AliTPCcalibDB *db=AliTPCcalibDB::Instance();
1569   Int_t nSectors=AliTPCROC::Instance()->GetNSectors(); 
1570   Int_t nTimeBins= 1100; // *Bug here - we should get NTimeBins from ALTRO - Parameters not relyable
1571   Int_t nPadRows=(AliTPCROC::Instance()->GetNRows(0) + AliTPCROC::Instance()->GetNRows(36));
1572   // parameters for filtering
1573   const Double_t nSigmaCut=9.;           // should be in recoParam ?
1574   const Double_t offsetTime=100;         // should be in RecoParam ?  -
1575   const Double_t offsetPadRow=300;       // should be in RecoParam ?
1576   const Double_t offsetTimeAccept=8;     // should be in RecoParam ?  - obtained as mean +1 rms in high IR pp
1577   TH2F hisTime("hisSectorTime","hisSectorTime", nSectors,0,nSectors, nTimeBins,0,nTimeBins);
1578   TH2F hisPadRow("hisSectorRow","hisSectorRow", nSectors,0,nSectors, nPadRows,0,nPadRows);
1579   //
1580   // 2.) Filling part -- loop over clusters
1581   //
1582   for (Int_t isector=0; isector<36; isector++){      // loop over sectors
1583     for (Int_t iside=0; iside<2; iside++){           // loop over sides A/C
1584       AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1585       Int_t nrows = sector.GetNRows();
1586       for (Int_t row = 0;row<nrows;row++){           // loop over rows       
1587         AliTPCtrackerRow&  tpcrow = sector[row];       
1588         Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
1589         if (iside>0) ncl=tpcrow.GetN2();
1590         for (Int_t i=0;i<ncl;i++) {  // loop over clusters
1591           AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1592           hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin()); 
1593           hisPadRow.Fill(cluster->GetDetector(),cluster->GetRow()); 
1594         }
1595       } 
1596     } 
1597   } 
1598
1599   //
1600   // 3. Filtering part
1601   //
1602   TVectorD vecTime(nTimeBins);
1603   TVectorD vecPadRow(nPadRows);
1604   TVectorD vecMedianSectorTime(nSectors);
1605   TVectorD vecRMSSectorTime(nSectors);
1606   TVectorD vecMedianSectorTimeOut6(nSectors);
1607   TVectorD vecMedianSectorTimeOut9(nSectors);//
1608   TVectorD vecMedianSectorTimeOut(nSectors);//
1609   TVectorD vecMedianSectorPadRow(nSectors);
1610   TVectorD vecRMSSectorPadRow(nSectors);
1611   TVectorD vecMedianSectorPadRowOut6(nSectors);
1612   TVectorD vecMedianSectorPadRowOut9(nSectors);
1613   TVectorD vecMedianSectorPadRowOut(nSectors);
1614   TVectorD vecSectorOut6(nSectors);
1615   TVectorD vecSectorOut9(nSectors);
1616   TMatrixD matSectorCluster(nSectors,2);
1617   //
1618   // 3.a)  median, rms calculations for hisTime 
1619   //
1620   for (Int_t isec=0; isec<nSectors; isec++){
1621     vecMedianSectorTimeOut6[isec]=0;
1622     vecMedianSectorTimeOut9[isec]=0;
1623     for (Int_t itime=0; itime<nTimeBins; itime++){
1624       vecTime[itime]=hisTime.GetBinContent(isec+1, itime+1);      
1625     }
1626     Double_t median= TMath::Mean(nTimeBins,vecTime.GetMatrixArray());
1627     Double_t rms= TMath::RMS(nTimeBins,vecTime.GetMatrixArray()); 
1628     vecMedianSectorTime[isec]=median;
1629     vecRMSSectorTime[isec]=rms;
1630     if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector TimeStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
1631     //
1632     // declare outliers
1633     for (Int_t itime=0; itime<nTimeBins; itime++){
1634       Double_t entries= hisTime.GetBinContent(isec+1, itime+1); 
1635       if (entries>median+6.*rms+offsetTime) {
1636         vecMedianSectorTimeOut6[isec]+=1;
1637       }
1638       if (entries>median+9.*rms+offsetTime) {
1639         vecMedianSectorTimeOut9[isec]+=1;
1640       }
1641     }
1642   }    
1643   //
1644   // 3.b) median, rms calculations for hisPadRow
1645   // 
1646   for (Int_t isec=0; isec<nSectors; isec++){
1647     vecMedianSectorPadRowOut6[isec]=0;
1648     vecMedianSectorPadRowOut9[isec]=0;
1649     for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1650       vecPadRow[ipadrow]=hisPadRow.GetBinContent(isec+1, ipadrow+1);      
1651     }
1652     Int_t nPadRowsSector= AliTPCROC::Instance()->GetNRows(isec);
1653     Double_t median= TMath::Mean(nPadRowsSector,vecPadRow.GetMatrixArray());
1654     Double_t rms= TMath::RMS(nPadRowsSector,vecPadRow.GetMatrixArray());
1655     vecMedianSectorPadRow[isec]=median;
1656     vecRMSSectorPadRow[isec]=rms;
1657     if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector PadRowStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
1658     //
1659     // declare outliers
1660     for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1661       Double_t entries= hisPadRow.GetBinContent(isec+1, ipadrow+1);
1662       if (entries>median+6.*rms+offsetPadRow) {
1663         vecMedianSectorPadRowOut6[isec]+=1;
1664       }
1665       if (entries>median+9.*rms+offsetPadRow) {
1666         vecMedianSectorPadRowOut9[isec]+=1;
1667       }
1668     }
1669   }
1670   //
1671   // 3.c) filter outlier sectors
1672   //
1673   Double_t medianSectorTime = TMath::Median(nSectors, vecTime.GetMatrixArray());
1674   Double_t mean69SectorTime, rms69SectorTime=0;
1675   AliMathBase::EvaluateUni(nSectors,  vecTime.GetMatrixArray(), mean69SectorTime,rms69SectorTime,69);
1676   for (Int_t isec=0; isec<nSectors; isec++){
1677     vecSectorOut6[isec]=0;
1678     vecSectorOut9[isec]=0;
1679     matSectorCluster(isec,0)=0;
1680     matSectorCluster(isec,1)=0;
1681     if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+6.*(rms69SectorTime+ offsetTimeAccept))) {
1682       vecSectorOut6[isec]=1;
1683     }
1684     if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+9.*(rms69SectorTime+ offsetTimeAccept))){
1685       vecSectorOut9[isec]=1;
1686     }
1687   }
1688   // light version of export variable
1689   Int_t filteredSector= vecSectorOut9.Sum();                  // light version of export variable
1690   Int_t filteredSectorTime= vecMedianSectorTimeOut9.Sum();
1691   Int_t filteredSectorPadRow= vecMedianSectorPadRowOut9.Sum();
1692   if (fEvent) if (fEvent->GetHeader()){
1693     fEvent->GetHeader()->SetTPCNoiseFilterCounter(0,TMath::Min(filteredSector,255));
1694     fEvent->GetHeader()->SetTPCNoiseFilterCounter(1,TMath::Min(filteredSectorTime,255));
1695     fEvent->GetHeader()->SetTPCNoiseFilterCounter(2,TMath::Min(filteredSectorPadRow,255));
1696   }
1697  
1698   //
1699   // 4. Disabling clusters in outlier layers
1700   //
1701   Int_t counterAll=0;
1702   Int_t counterOut=0;
1703   for (Int_t isector=0; isector<36; isector++){      // loop over sectors
1704     for (Int_t iside=0; iside<2; iside++){           // loop over sides A/C
1705       AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1706       Int_t nrows = sector.GetNRows();
1707       for (Int_t row = 0;row<nrows;row++){           // loop over rows       
1708         AliTPCtrackerRow&  tpcrow = sector[row];       
1709         Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
1710         if (iside>0) ncl=tpcrow.GetN2();
1711         for (Int_t i=0;i<ncl;i++) {  // loop over clusters
1712           AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1713           Double_t medianTime=vecMedianSectorTime[cluster->GetDetector()];
1714           Double_t medianPadRow=vecMedianSectorPadRow[cluster->GetDetector()];
1715           Double_t rmsTime=vecRMSSectorTime[cluster->GetDetector()];
1716           Double_t rmsPadRow=vecRMSSectorPadRow[cluster->GetDetector()];
1717           Int_t entriesPadRow=hisPadRow.GetBinContent(cluster->GetDetector()+1, cluster->GetRow()+1);
1718           Int_t entriesTime=hisTime.GetBinContent(cluster->GetDetector()+1, cluster->GetTimeBin()+1);
1719           Bool_t isOut=kFALSE;
1720           if (vecSectorOut9[cluster->GetDetector()]>0.5) {
1721             isOut=kTRUE;
1722           }
1723           
1724           if (entriesTime>medianTime+nSigmaCut*rmsTime+offsetTime) {
1725             isOut=kTRUE;
1726             vecMedianSectorTimeOut[cluster->GetDetector()]++;
1727           }
1728           if (entriesPadRow>medianPadRow+nSigmaCut*rmsPadRow+offsetPadRow) {
1729             isOut=kTRUE;
1730             vecMedianSectorPadRowOut[cluster->GetDetector()]++;
1731           }
1732           counterAll++;
1733           matSectorCluster(cluster->GetDetector(),0)+=1;
1734           if (isOut){
1735             cluster->Disable();
1736             counterOut++;
1737             matSectorCluster(cluster->GetDetector(),1)+=1;
1738           }
1739         }
1740       }
1741     }
1742   }
1743   for (Int_t isec=0; isec<nSectors; isec++){
1744     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());
1745   }
1746   //
1747   // dump info to streamer - for later tuning of cuts
1748   //
1749   if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) {  // stream TPC data ouliers filtering infomation
1750     AliLog::Flush();
1751     AliInfo(TString::Format("Cluster counter: (%d/%d) (Filtered/All)",counterOut,counterAll).Data());
1752     for (Int_t iSec=0; iSec<nSectors; iSec++){
1753       if (vecSectorOut9[iSec]>0 ||  matSectorCluster(iSec,1)>0) {
1754         AliInfo(TString::Format("Filtered sector\t%d",iSec).Data());
1755         Double_t vecMedTime =TMath::Median(72,vecMedianSectorTime.GetMatrixArray());
1756         Double_t vecMedPadRow =TMath::Median(72,vecMedianSectorPadRow.GetMatrixArray());
1757         Double_t vecMedCluster=(counterAll-counterOut)/72;
1758         AliInfo(TString::Format("VecMedianSectorTime\t(%4.4f/%4.4f/%4.4f)",       vecMedianSectorTimeOut[iSec],vecMedianSectorTime[iSec],vecMedTime).Data());
1759         AliInfo(TString::Format("VecMedianSectorPadRow\t(%4.4f/%4.4f/%4.4f)",     vecMedianSectorPadRowOut[iSec],vecMedianSectorPadRow[iSec],vecMedPadRow).Data());
1760         AliInfo(TString::Format("MatSectorCluster\t(%4.4f/%4.4f/%4.4f)\n",          matSectorCluster(iSec,1), matSectorCluster(iSec,0),  vecMedCluster).Data());
1761         AliLog::Flush();
1762       }
1763     }
1764     AliLog::Flush();
1765     Int_t eventNr = fEvent->GetEventNumberInFile();
1766     (*fDebugStreamer)<<"filterClusterInfo"<<
1767       // minimal set variables for the ESDevent
1768       "eventNr="<<eventNr<<
1769       "counterAll="<<counterAll<<
1770       "counterOut="<<counterOut<<
1771       "matSectotCluster.="<<&matSectorCluster<<                   // 
1772       //
1773       "filteredSector="<<filteredSector<<                        //  counter filtered sectors                   
1774       "filteredSectorTime="<<filteredSectorTime<<                //  counter filtered time bins
1775       "filteredSectorPadRow="<<filteredSectorPadRow<<            //  counter filtered pad-rows
1776       // per sector outlier information
1777       "medianSectorTime="<<medianSectorTime<<                    // median number of clusters per sector/timebin
1778       "mean69SectorTime="<<mean69SectorTime<<                    // LTM statistic  mean of clusters per sector/timebin
1779       "rms69SectorTime="<<rms69SectorTime<<                      // LTM statistic  RMS of clusters per sector/timebin
1780       "vecSectorOut6.="<<&vecSectorOut6<<                        // flag array sector  - 6 sigma +accept margin outlier
1781       "vecSectorOut9.="<<&vecSectorOut9<<                        // flag array sector  - 9 sigma + accept margin outlier
1782       // per sector/timebin outlier detection
1783       "vecMedianSectorTime.="<<&vecMedianSectorTime<<
1784       "vecRMSSectorTime.="<<&vecRMSSectorTime<<
1785       "vecMedianSectorTimeOut6.="<<&vecMedianSectorTimeOut6<<
1786       "vecMedianSectorTimeOut9.="<<&vecMedianSectorTimeOut9<<
1787       "vecMedianSectorTimeOut0.="<<&vecMedianSectorTimeOut<<
1788       // per sector/pad-row outlier detection
1789       "vecMedianSectorPadRow.="<<&vecMedianSectorPadRow<<
1790       "vecRMSSectorPadRow.="<<&vecRMSSectorPadRow<<
1791       "vecMedianSectorPadRowOut6.="<<&vecMedianSectorPadRowOut6<<
1792       "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut9<<
1793       "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut<<
1794       "\n";
1795     ((*fDebugStreamer)<<"filterClusterInfo").GetTree()->Write();
1796     fDebugStreamer->GetFile()->Flush();
1797   }
1798 }
1799
1800 void AliTPCtracker::UnloadClusters()
1801 {
1802   //
1803   // unload clusters from the memory
1804   //
1805   Int_t nrows = fOuterSec->GetNRows();
1806   for (Int_t sec = 0;sec<fkNOS;sec++)
1807     for (Int_t row = 0;row<nrows;row++){
1808       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1809       //      if (tpcrow){
1810       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1811       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1812       //}
1813       tpcrow->ResetClusters();
1814     }
1815   //
1816   nrows = fInnerSec->GetNRows();
1817   for (Int_t sec = 0;sec<fkNIS;sec++)
1818     for (Int_t row = 0;row<nrows;row++){
1819       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1820       //if (tpcrow){
1821       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1822       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1823       //}
1824       tpcrow->ResetClusters();
1825     }
1826
1827   return ;
1828 }
1829
1830 void AliTPCtracker::FillClusterArray(TObjArray* array) const{
1831   //
1832   // Filling cluster to the array - For visualization purposes
1833   //
1834   Int_t nrows=0;
1835   nrows = fOuterSec->GetNRows();
1836   for (Int_t sec = 0;sec<fkNOS;sec++)
1837     for (Int_t row = 0;row<nrows;row++){
1838       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1839       if (!tpcrow) continue;
1840       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1841         array->AddLast((TObject*)((*tpcrow)[icl]));
1842       }
1843     } 
1844   nrows = fInnerSec->GetNRows();
1845   for (Int_t sec = 0;sec<fkNIS;sec++)
1846     for (Int_t row = 0;row<nrows;row++){
1847       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1848       if (!tpcrow) continue;
1849       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1850         array->AddLast((TObject*)(*tpcrow)[icl]);
1851       }
1852     }
1853 }
1854
1855
1856 void   AliTPCtracker::Transform(AliTPCclusterMI * cluster){
1857   //
1858   // transformation
1859   //
1860   AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1861   AliTPCTransform *transform = calibDB->GetTransform() ;
1862   if (!transform) {
1863     AliFatal("Tranformations not in calibDB");
1864     return;
1865   }
1866   transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1867   Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),static_cast<Double_t>(cluster->GetPad()),static_cast<Double_t>(cluster->GetTimeBin())};
1868   Int_t i[1]={cluster->GetDetector()};
1869   transform->Transform(x,i,0,1);  
1870   //  if (cluster->GetDetector()%36>17){
1871   //  x[1]*=-1;
1872   //}
1873
1874   //
1875   // in debug mode  check the transformation
1876   //
1877   if ((AliTPCReconstructor::StreamLevel()&kStreamTransform)>0) { 
1878     Float_t gx[3];
1879     cluster->GetGlobalXYZ(gx);
1880     Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1881     TTreeSRedirector &cstream = *fDebugStreamer;
1882     cstream<<"Transform"<<  // needed for debugging of the cluster transformation, resp. used for later visualization 
1883       "event="<<event<<
1884       "x0="<<x[0]<<
1885       "x1="<<x[1]<<
1886       "x2="<<x[2]<<
1887       "gx0="<<gx[0]<<
1888       "gx1="<<gx[1]<<
1889       "gx2="<<gx[2]<<
1890       "Cl.="<<cluster<<
1891       "\n"; 
1892   }
1893   cluster->SetX(x[0]);
1894   cluster->SetY(x[1]);
1895   cluster->SetZ(x[2]);
1896   // The old stuff:
1897   //
1898   // 
1899   //
1900   //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1901   if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1902     TGeoHMatrix  *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1903     //TGeoHMatrix  mat;
1904     Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1905     Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1906     if (mat) mat->LocalToMaster(pos,posC);
1907     else{
1908       // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1909     }
1910     cluster->SetX(posC[0]);
1911     cluster->SetY(posC[1]);
1912     cluster->SetZ(posC[2]);
1913   }
1914 }
1915
1916 void  AliTPCtracker::ApplyXtalkCorrection(){
1917   //
1918   // ApplyXtalk correction 
1919   // Loop over all clusters
1920   //      add to each cluster signal corresponding to common Xtalk mode for given time bin at given wire segment
1921   // cluster loop
1922   for (Int_t isector=0; isector<36; isector++){  //loop tracking sectors
1923     for (Int_t iside=0; iside<2; iside++){       // loop over sides A/C
1924       AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1925       Int_t nrows     = sector.GetNRows();       
1926       for (Int_t row = 0;row<nrows;row++){           // loop over rows       
1927         AliTPCtrackerRow&  tpcrow = sector[row];     // row object   
1928         Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
1929         if (iside>0) ncl=tpcrow.GetN2();
1930         Int_t xSector=0;    // sector number in the TPC convention 0-72
1931         if (isector<18){  //if IROC
1932           xSector=isector+(iside>0)*18;
1933         }else{
1934           xSector=isector+18;  // isector -18 +36   
1935           if (iside>0) xSector+=18;
1936         }       
1937         TMatrixD &crossTalkMatrix= *((TMatrixD*)fCrossTalkSignalArray->At(xSector));
1938         Int_t wireSegmentID     = fkParam->GetWireSegment(xSector,row);
1939         for (Int_t i=0;i<ncl;i++) {
1940           AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1941           Int_t iTimeBin=TMath::Nint(cluster->GetTimeBin());
1942           Double_t xTalk= crossTalkMatrix[wireSegmentID][iTimeBin];
1943           cluster->SetMax(cluster->GetMax()+xTalk);
1944           const Double_t kDummy=4;
1945           Double_t sumxTalk=xTalk*kDummy; // should be calculated via time response function
1946           cluster->SetQ(cluster->GetQ()+sumxTalk);
1947
1948
1949           if ((AliTPCReconstructor::StreamLevel()&kStreamXtalk)>0) {  // flag: stream crosstalk correctio as applied to cluster
1950             TTreeSRedirector &cstream = *fDebugStreamer;
1951             if (gRandom->Rndm() > 0.){
1952               cstream<<"Xtalk"<<
1953                 "isector=" << isector <<               // sector [0,36]
1954                 "iside=" << iside <<                   // side A or C
1955                 "row=" << row <<                       // padrow
1956                 "i=" << i <<                           // index of the cluster 
1957                 "xSector=" << xSector <<               // sector [0,72] 
1958                 "wireSegmentID=" << wireSegmentID <<   // anode wire segment id [0,10] 
1959                 "iTimeBin=" << iTimeBin <<             // timebin of the corrected cluster 
1960                 "xTalk=" << xTalk <<                   // Xtalk contribution added to Qmax
1961                 "sumxTalk=" << sumxTalk <<             // Xtalk contribution added to Qtot (roughly 3*Xtalk) 
1962                 "cluster.=" << cluster <<              // corrected cluster object 
1963                 "\n";
1964             }
1965           }// dump the results to the debug streamer if in debug mode
1966
1967
1968
1969
1970
1971         }
1972       }
1973     }
1974   }
1975 }
1976
1977 void  AliTPCtracker::ApplyTailCancellation(){
1978   //
1979   // Correct the cluster charge for the ion tail effect 
1980   // The TimeResponse function accessed via  AliTPCcalibDB (TPC/Calib/IonTail)
1981   //
1982
1983   // Retrieve
1984   TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
1985   if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
1986   TObject *rocFactorIROC  = ionTailArr->FindObject("factorIROC");
1987   TObject *rocFactorOROC  = ionTailArr->FindObject("factorOROC");   
1988   Float_t factorIROC      = (atof(rocFactorIROC->GetTitle()));
1989   Float_t factorOROC      = (atof(rocFactorOROC->GetTitle()));
1990
1991   // find the number of clusters for the whole TPC (nclALL)
1992   Int_t nclALL=0;
1993   for (Int_t isector=0; isector<36; isector++){
1994     AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1995     nclALL += sector.GetNClInSector(0);
1996     nclALL += sector.GetNClInSector(1);
1997   }
1998
1999   // start looping over all clusters 
2000   for (Int_t iside=0; iside<2; iside++){    // loop over sides
2001     //
2002     //
2003     for (Int_t secType=0; secType<2; secType++){  //loop over inner or outer sector
2004       // cache experimantal tuning factor for the different chamber type 
2005       const Float_t ampfactor = (secType==0)?factorIROC:factorOROC;
2006       std::cout << " ampfactor = " << ampfactor << std::endl;
2007       //
2008       for (Int_t sec = 0;sec<fkNOS;sec++){        //loop overs sectors
2009         //
2010         //
2011         // Cache time response functions and their positons to COG of the cluster        
2012         TGraphErrors ** graphRes   = new TGraphErrors *[20];
2013         Float_t * indexAmpGraphs   = new Float_t[20];      
2014         for (Int_t icache=0; icache<20; icache++) 
2015         {
2016           graphRes[icache]       = NULL;
2017           indexAmpGraphs[icache] = 0;
2018         }
2019         /////////////////////////////  --> position fo sie loop
2020         if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(sec+36*secType+18*iside,graphRes,indexAmpGraphs))
2021         {
2022           continue;
2023         }
2024         
2025         AliTPCtrackerSector &sector= (secType==0)?fInnerSec[sec]:fOuterSec[sec];  
2026         Int_t nrows     = sector.GetNRows();                                       // number of rows
2027         Int_t nclSector = sector.GetNClInSector(iside);                            // ncl per sector to be used for debugging
2028
2029         for (Int_t row = 0;row<nrows;row++){           // loop over rows
2030
2031           AliTPCtrackerRow&  tpcrow = sector[row];     // row object   
2032           Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
2033           if (iside>0) ncl=tpcrow.GetN2();
2034         
2035           // Order clusters in time for the proper correction of ion tail
2036           Float_t qTotArray[ncl];                      // arrays to be filled with modified Qtot and Qmax values in order to avoid float->int conversion  
2037           Float_t qMaxArray[ncl];
2038           Int_t sortedClusterIndex[ncl];
2039           Float_t sortedClusterTimeBin[ncl];
2040           TObjArray *rowClusterArray = new TObjArray(ncl);  // cache clusters for each row  
2041           for (Int_t i=0;i<ncl;i++) 
2042           {
2043             qTotArray[i]=0;
2044             qMaxArray[i]=0;
2045             sortedClusterIndex[i]=i;
2046             AliTPCclusterMI *rowcl= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
2047             if (rowcl) {
2048               rowClusterArray->AddAt(rowcl,i);
2049             } else {
2050               rowClusterArray->RemoveAt(i);
2051             }
2052             // Fill the timebin info to the array in order to sort wrt tb
2053             if (!rowcl) {
2054               sortedClusterTimeBin[i]=0.0;
2055             } else {
2056             sortedClusterTimeBin[i] = rowcl->GetTimeBin();
2057             }
2058
2059           } 
2060           TMath::Sort(ncl,sortedClusterTimeBin,sortedClusterIndex,kFALSE);       // sort clusters in time
2061      
2062           // Main cluster correction loops over clusters
2063           for (Int_t icl0=0; icl0<ncl;icl0++){    // first loop over clusters
2064
2065             AliTPCclusterMI *cl0= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl0]));
2066             
2067             if (!cl0) continue;
2068             Int_t nclPad=0;                       
2069             for (Int_t icl1=0; icl1<ncl;icl1++){  // second loop over clusters     
2070               AliTPCclusterMI *cl1= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl1]));
2071               if (!cl1) continue;
2072               if (TMath::Abs(cl0->GetPad()-cl1->GetPad())>4) continue;           // no contribution if far away in pad direction
2073               if (cl0->GetTimeBin()<= cl1->GetTimeBin()) continue;               // no contibution to the tail if later
2074               if (TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin())>600) continue; // out of the range of response function
2075
2076               if (TMath::Abs(cl0->GetPad()-cl1->GetPad())<4) nclPad++;           // count ncl for every pad for debugging
2077             
2078               // Get the correction values for Qmax and Qtot and find total correction for a given cluster
2079               Double_t ionTailMax=0.;  
2080               Double_t ionTailTotal=0.;  
2081               GetTailValue(ampfactor,ionTailMax,ionTailTotal,graphRes,indexAmpGraphs,cl0,cl1);
2082               ionTailMax=TMath::Abs(ionTailMax);
2083               ionTailTotal=TMath::Abs(ionTailTotal);
2084               qTotArray[icl0]+=ionTailTotal;
2085               qMaxArray[icl0]+=ionTailMax;
2086
2087               // Dump some info for debugging while clusters are being corrected
2088               if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {  // flag: stream ion tail correction  as applied to cluster
2089                 TTreeSRedirector &cstream = *fDebugStreamer;
2090                 if (gRandom->Rndm() > 0.999){
2091                   cstream<<"IonTail"<<
2092                       "cl0.="         <<cl0          <<   // cluster 0 (to be corrected)
2093                       "cl1.="         <<cl1          <<   // cluster 1 (previous cluster)
2094                       "ionTailTotal=" <<ionTailTotal <<   // ion Tail from cluster 1 contribution to cluster0
2095                       "ionTailMax="   <<ionTailMax   <<   // ion Tail from cluster 1 contribution to cluster0 
2096                       "\n";
2097                 }
2098               }// dump the results to the debug streamer if in debug mode
2099             
2100             }//end of second loop over clusters
2101             
2102             // Set corrected values of the corrected cluster          
2103             cl0->SetQ(TMath::Nint(Float_t(cl0->GetQ())+Float_t(qTotArray[icl0])));
2104             cl0->SetMax(TMath::Nint(Float_t(cl0->GetMax())+qMaxArray[icl0]));
2105           
2106             // Dump some info for debugging after clusters are corrected 
2107             if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {
2108               TTreeSRedirector &cstream = *fDebugStreamer;
2109               if (gRandom->Rndm() > 0.999){
2110               cstream<<"IonTailCorrected"<<
2111                   "cl0.="                     << cl0              <<   // cluster 0 with huge Qmax
2112                   "ionTailTotalPerCluster="   << qTotArray[icl0]  <<
2113                   "ionTailMaxPerCluster="     << qMaxArray[icl0]  <<
2114                   "nclALL="                   << nclALL           <<
2115                   "nclSector="                << nclSector        <<
2116                   "nclRow="                   << ncl              <<
2117                   "nclPad="                   << nclPad           <<
2118                   "row="                      << row              <<
2119                   "sector="                   << sec              <<
2120                   "icl0="                     << icl0             <<
2121                   "\n";
2122               }
2123             }// dump the results to the debug streamer if in debug mode
2124           
2125           }//end of first loop over cluster
2126           delete rowClusterArray;
2127         }//end of loop over rows
2128         for (int i=0; i<20; i++) delete graphRes[i];
2129         delete [] graphRes;
2130         delete [] indexAmpGraphs;
2131       
2132       }//end of loop over sectors
2133     }//end of loop over IROC/OROC
2134   }// end of side loop
2135 }
2136 //_____________________________________________________________________________
2137 void AliTPCtracker::GetTailValue(Float_t ampfactor,Double_t &ionTailMax, Double_t &ionTailTotal,TGraphErrors **graphRes,Float_t *indexAmpGraphs,AliTPCclusterMI *cl0,AliTPCclusterMI *cl1){
2138
2139   //
2140   // Function in order to calculate the amount of the correction to be added for a given cluster, return values are ionTailTaoltal and ionTailMax
2141   // Parameters:
2142   // cl0 -  cluster to be modified
2143   // cl1 -  source cluster ion tail of this cluster will be added to the cl0 (accroding time and pad response function)
2144   // 
2145   const Double_t kMinPRF    = 0.5;                           // minimal PRF width
2146   ionTailTotal              = 0.;                            // correction value to be added to Qtot of cl0
2147   ionTailMax                = 0.;                            // correction value to be added to Qmax of cl0
2148
2149   Float_t qTot0             =  cl0->GetQ();                  // cl0 Qtot info
2150   Float_t qTot1             =  cl1->GetQ();                  // cl1 Qtot info
2151   Int_t sectorPad           =  cl1->GetDetector();           // sector number
2152   Int_t padcl0              =  TMath::Nint(cl0->GetPad());   // pad0
2153   Int_t padcl1              =  TMath::Nint(cl1->GetPad());   // pad1
2154   Float_t padWidth          = (sectorPad < 36)?0.4:0.6;      // pad width in cm
2155   const Int_t deltaTimebin  =  TMath::Nint(TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin()))+12;  //distance between pads of cl1 and cl0 increased by 12 bins
2156   Double_t rmsPad1          = (cl1->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl1->GetSigmaY2())/padWidth);
2157   Double_t rmsPad0          = (cl0->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl0->GetSigmaY2())/padWidth);
2158   
2159    
2160   
2161   Double_t sumAmp1=0.;
2162   for (Int_t idelta =-2; idelta<=2;idelta++){
2163     sumAmp1+=TMath::Exp(-idelta*idelta/(2*rmsPad1));
2164   }
2165
2166   Double_t sumAmp0=0.;
2167   for (Int_t idelta =-2; idelta<=2;idelta++){
2168     sumAmp0+=TMath::Exp(-idelta*idelta/(2*rmsPad0));
2169   }
2170
2171   // Apply the correction  -->   cl1 corrects cl0 (loop over cl1's pads and find which pads of cl0 are going to be corrected)
2172   Int_t padScan=2;      // +-2 pad-timebin window will be scanned
2173   for (Int_t ipad1=padcl1-padScan; ipad1<=padcl1+padScan; ipad1++) {
2174     //
2175     //
2176     Float_t  deltaPad1  = TMath::Abs(cl1->GetPad()-(Float_t)ipad1);
2177     Double_t amp1       = (TMath::Exp(-(deltaPad1*deltaPad1)/(2*rmsPad1)))/sumAmp1;  // normalized pad response function
2178     Float_t qTotPad1    = amp1*qTot1;                                               // used as a factor to multipliy the response function
2179       
2180     // find closest value of cl1 to COG (among the time response functions' amplitude array --> to select proper t.r.f.)
2181     Int_t ampIndex = 0;
2182     Float_t diffAmp  = TMath::Abs(deltaPad1-indexAmpGraphs[0]);
2183     for (Int_t j=0;j<20;j++) {
2184       if (diffAmp > TMath::Abs(deltaPad1-indexAmpGraphs[j]) && indexAmpGraphs[j]!=0)
2185         {
2186           diffAmp  = TMath::Abs(deltaPad1-indexAmpGraphs[j]);
2187           ampIndex = j;
2188         }
2189     }
2190     if (!graphRes[ampIndex]) continue;
2191     if (deltaTimebin+2 >= graphRes[ampIndex]->GetN()) continue;
2192     if (graphRes[ampIndex]->GetY()[deltaTimebin+2]>=0) continue;
2193      
2194     for (Int_t ipad0=padcl0-padScan; ipad0<=padcl0+padScan; ipad0++) {
2195       //
2196       //
2197       if (ipad1!=ipad0) continue;                                     // check if ipad1 channel sees ipad0 channel, if not no correction to be applied.
2198       
2199       Float_t deltaPad0  = TMath::Abs(cl0->GetPad()-(Float_t)ipad0);
2200       Double_t amp0      = (TMath::Exp(-(deltaPad0*deltaPad0)/(2*rmsPad0)))/sumAmp0;  // normalized pad resp function
2201       Float_t qMaxPad0   = amp0*qTot0;
2202            
2203       // Add 5 timebin range contribution around the max peak (-+2 tb window)
2204       for (Int_t itb=deltaTimebin-2; itb<=deltaTimebin+2; itb++) {
2205
2206         if (itb<0) continue; 
2207         if (itb>=graphRes[ampIndex]->GetN()) continue;
2208        
2209         // calculate contribution to qTot
2210         Float_t tailCorr =  TMath::Abs((qTotPad1*ampfactor)*(graphRes[ampIndex])->GetY()[itb]);
2211         if (ipad1!=padcl0) { 
2212           ionTailTotal += TMath::Min(qMaxPad0,tailCorr);   // for side pad
2213         } else {             
2214           ionTailTotal += tailCorr;                        // for center pad
2215         }
2216         // calculate contribution to qMax
2217         if (itb == deltaTimebin && ipad1 == padcl0) ionTailMax += tailCorr;   
2218         
2219       } // end of tb correction loop which is applied over 5 tb range
2220
2221     } // end of cl0 loop
2222   } // end of cl1 loop
2223   
2224 }
2225
2226 //_____________________________________________________________________________
2227 Int_t AliTPCtracker::LoadOuterSectors() {
2228   //-----------------------------------------------------------------
2229   // This function fills outer TPC sectors with clusters.
2230   //-----------------------------------------------------------------
2231   Int_t nrows = fOuterSec->GetNRows();
2232   UInt_t index=0;
2233   for (Int_t sec = 0;sec<fkNOS;sec++)
2234     for (Int_t row = 0;row<nrows;row++){
2235       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);  
2236       Int_t sec2 = sec+2*fkNIS;
2237       //left
2238       Int_t ncl = tpcrow->GetN1();
2239       while (ncl--) {
2240         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
2241         index=(((sec2<<8)+row)<<16)+ncl;
2242         tpcrow->InsertCluster(c,index);
2243       }
2244       //right
2245       ncl = tpcrow->GetN2();
2246       while (ncl--) {
2247         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
2248         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
2249         tpcrow->InsertCluster(c,index);
2250       }
2251       //
2252       // write indexes for fast acces
2253       //
2254       for (Int_t i=0;i<510;i++)
2255         tpcrow->SetFastCluster(i,-1);
2256       for (Int_t i=0;i<tpcrow->GetN();i++){
2257         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
2258         tpcrow->SetFastCluster(zi,i);  // write index
2259       }
2260       Int_t last = 0;
2261       for (Int_t i=0;i<510;i++){
2262         if (tpcrow->GetFastCluster(i)<0)
2263           tpcrow->SetFastCluster(i,last);
2264         else
2265           last = tpcrow->GetFastCluster(i);
2266       }
2267     }  
2268   fN=fkNOS;
2269   fSectors=fOuterSec;
2270   return 0;
2271 }
2272
2273
2274 //_____________________________________________________________________________
2275 Int_t  AliTPCtracker::LoadInnerSectors() {
2276   //-----------------------------------------------------------------
2277   // This function fills inner TPC sectors with clusters.
2278   //-----------------------------------------------------------------
2279   Int_t nrows = fInnerSec->GetNRows();
2280   UInt_t index=0;
2281   for (Int_t sec = 0;sec<fkNIS;sec++)
2282     for (Int_t row = 0;row<nrows;row++){
2283       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
2284       //
2285       //left
2286       Int_t ncl = tpcrow->GetN1();
2287       while (ncl--) {
2288         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
2289         index=(((sec<<8)+row)<<16)+ncl;
2290         tpcrow->InsertCluster(c,index);
2291       }
2292       //right
2293       ncl = tpcrow->GetN2();
2294       while (ncl--) {
2295         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
2296         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
2297         tpcrow->InsertCluster(c,index);
2298       }
2299       //
2300       // write indexes for fast acces
2301       //
2302       for (Int_t i=0;i<510;i++)
2303         tpcrow->SetFastCluster(i,-1);
2304       for (Int_t i=0;i<tpcrow->GetN();i++){
2305         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
2306         tpcrow->SetFastCluster(zi,i);  // write index
2307       }
2308       Int_t last = 0;
2309       for (Int_t i=0;i<510;i++){
2310         if (tpcrow->GetFastCluster(i)<0)
2311           tpcrow->SetFastCluster(i,last);
2312         else
2313           last = tpcrow->GetFastCluster(i);
2314       }
2315
2316     }  
2317    
2318   fN=fkNIS;
2319   fSectors=fInnerSec;
2320   return 0;
2321 }
2322
2323
2324
2325 //_________________________________________________________________________
2326 AliTPCclusterMI *AliTPCtracker::GetClusterMI(Int_t index) const {
2327   //--------------------------------------------------------------------
2328   //       Return pointer to a given cluster
2329   //--------------------------------------------------------------------
2330   if (index<0) return 0; // no cluster
2331   Int_t sec=(index&0xff000000)>>24; 
2332   Int_t row=(index&0x00ff0000)>>16; 
2333   Int_t ncl=(index&0x00007fff)>>00;
2334
2335   const AliTPCtrackerRow * tpcrow=0;
2336   TClonesArray * clrow =0;
2337
2338   if (sec<0 || sec>=fkNIS*4) {
2339     AliWarning(Form("Wrong sector %d",sec));
2340     return 0x0;
2341   }
2342
2343   if (sec<fkNIS*2){
2344     AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
2345     if (tracksec.GetNRows()<=row) return 0;
2346     tpcrow = &(tracksec[row]);
2347     if (tpcrow==0) return 0;
2348
2349     if (sec<fkNIS) {
2350       if (tpcrow->GetN1()<=ncl) return 0;
2351       clrow = tpcrow->GetClusters1();
2352     }
2353     else {
2354       if (tpcrow->GetN2()<=ncl) return 0;
2355       clrow = tpcrow->GetClusters2();
2356     }
2357   }
2358   else {
2359     AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
2360     if (tracksec.GetNRows()<=row) return 0;
2361     tpcrow = &(tracksec[row]);
2362     if (tpcrow==0) return 0;
2363
2364     if (sec-2*fkNIS<fkNOS) {
2365       if (tpcrow->GetN1()<=ncl) return 0;
2366       clrow = tpcrow->GetClusters1();
2367     }
2368     else {
2369       if (tpcrow->GetN2()<=ncl) return 0;
2370       clrow = tpcrow->GetClusters2();
2371     }
2372   }
2373
2374   return (AliTPCclusterMI*)clrow->At(ncl);
2375   
2376 }
2377
2378
2379
2380 Int_t AliTPCtracker::FollowToNext(AliTPCseed& t, Int_t nr) {
2381   //-----------------------------------------------------------------
2382   // This function tries to find a track prolongation to next pad row
2383   //-----------------------------------------------------------------
2384   //
2385   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
2386   //
2387   //
2388   AliTPCclusterMI *cl=0;
2389   Int_t tpcindex= t.GetClusterIndex2(nr);
2390   //
2391   // update current shape info every 5 pad-row
2392   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
2393     GetShape(&t,nr);    
2394     //}
2395   //  
2396   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters 
2397     //        
2398     if (tpcindex==-1) return 0; //track in dead zone
2399     if (tpcindex >= 0){     //
2400       cl = t.GetClusterPointer(nr);
2401       //if (cl==0) cl = GetClusterMI(tpcindex);
2402       if (!cl) cl = GetClusterMI(tpcindex);
2403       t.SetCurrentClusterIndex1(tpcindex); 
2404     }
2405     if (cl){      
2406       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
2407       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
2408       //
2409       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
2410       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
2411       
2412       if (TMath::Abs(angle-t.GetAlpha())>0.001){
2413         Double_t rotation = angle-t.GetAlpha();
2414         t.SetRelativeSector(relativesector);
2415         if (!t.Rotate(rotation)) {
2416           t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2417           return 0;
2418         }       
2419       }
2420       if (!t.PropagateTo(x)) {
2421         t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2422         return 0;
2423       }
2424       //
2425       t.SetCurrentCluster(cl); 
2426       t.SetRow(nr);
2427       Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
2428       if ((tpcindex&0x8000)==0) accept =0;
2429       if (accept<3) { 
2430         //if founded cluster is acceptible
2431         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
2432           t.SetErrorY2(t.GetErrorY2()+0.03);
2433           t.SetErrorZ2(t.GetErrorZ2()+0.03); 
2434           t.SetErrorY2(t.GetErrorY2()*3);
2435           t.SetErrorZ2(t.GetErrorZ2()*3); 
2436         }
2437         t.SetNFoundable(t.GetNFoundable()+1);
2438         UpdateTrack(&t,accept);
2439         return 1;
2440       }
2441       else { // Remove old cluster from track
2442         t.SetClusterIndex(nr, -3);
2443         t.SetClusterPointer(nr, 0);
2444       }
2445     }
2446   }
2447   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;  // cut on angle
2448   if (fIteration>1 && IsFindable(t)){
2449     // not look for new cluster during refitting    
2450     t.SetNFoundable(t.GetNFoundable()+1);
2451     return 0;
2452   }
2453   //
2454   UInt_t index=0;
2455   //  if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
2456   if (!t.PropagateTo(x)) {
2457     if (fIteration==0) t.SetRemoval(10);
2458     return 0;
2459   }
2460   Double_t y = t.GetY(); 
2461   if (TMath::Abs(y)>ymax){
2462     if (y > ymax) {
2463       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
2464       if (!t.Rotate(fSectors->GetAlpha())) 
2465         return 0;
2466     } else if (y <-ymax) {
2467       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
2468       if (!t.Rotate(-fSectors->GetAlpha())) 
2469         return 0;
2470     }
2471     if (!t.PropagateTo(x)) {
2472       if (fIteration==0) t.SetRemoval(10);
2473       return 0;
2474     }
2475     y = t.GetY();
2476   }
2477   //
2478   Double_t z=t.GetZ();
2479   //
2480
2481   if (!IsActive(t.GetRelativeSector(),nr)) {
2482     t.SetInDead(kTRUE);
2483     t.SetClusterIndex2(nr,-1); 
2484     return 0;
2485   }
2486   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
2487   Bool_t isActive  = IsActive(t.GetRelativeSector(),nr);
2488   Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
2489   
2490   if (!isActive || !isActive2) return 0;
2491
2492   const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
2493   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
2494   Double_t  roady  =1.;
2495   Double_t  roadz = 1.;
2496   //
2497   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2498     t.SetInDead(kTRUE);
2499     t.SetClusterIndex2(nr,-1); 
2500     return 0;
2501   } 
2502   else
2503     {
2504       if (IsFindable(t))
2505           //      if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) 
2506         t.SetNFoundable(t.GetNFoundable()+1);
2507       else
2508         return 0;
2509     }   
2510   //calculate 
2511   if (krow) {
2512     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);    
2513     cl = krow.FindNearest2(y,z,roady,roadz,index);    
2514     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));       
2515   }  
2516   if (cl) {
2517     t.SetCurrentCluster(cl); 
2518     t.SetRow(nr);
2519     if (fIteration==2&&cl->IsUsed(10)) return 0; 
2520     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
2521     if (fIteration==2&&cl->IsUsed(11)) {
2522       t.SetErrorY2(t.GetErrorY2()+0.03);
2523       t.SetErrorZ2(t.GetErrorZ2()+0.03); 
2524       t.SetErrorY2(t.GetErrorY2()*3);
2525       t.SetErrorZ2(t.GetErrorZ2()*3); 
2526     }
2527     /*    
2528     if (t.fCurrentCluster->IsUsed(10)){
2529       //
2530       //     
2531
2532       t.fNShared++;
2533       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
2534         t.fRemoval =10;
2535         return 0;
2536       }
2537     }
2538     */
2539     if (accept<3) UpdateTrack(&t,accept);
2540
2541   } else {  
2542     if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
2543     
2544   }
2545   return 1;
2546 }
2547
2548
2549
2550 //_________________________________________________________________________
2551 Bool_t AliTPCtracker::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
2552 {
2553   // Get track space point by index
2554   // return false in case the cluster doesn't exist
2555   AliTPCclusterMI *cl = GetClusterMI(index);
2556   if (!cl) return kFALSE;
2557   Int_t sector = (index&0xff000000)>>24;
2558   //  Int_t row = (index&0x00ff0000)>>16;
2559   Float_t xyz[3];
2560   //  xyz[0] = fkParam->GetPadRowRadii(sector,row);
2561   xyz[0] = cl->GetX();
2562   xyz[1] = cl->GetY();
2563   xyz[2] = cl->GetZ();
2564   Float_t sin,cos;
2565   fkParam->AdjustCosSin(sector,cos,sin);
2566   Float_t x = cos*xyz[0]-sin*xyz[1];
2567   Float_t y = cos*xyz[1]+sin*xyz[0];
2568   Float_t cov[6];
2569   Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
2570   if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
2571   Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
2572   if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
2573   cov[0] = sin*sin*sigmaY2;
2574   cov[1] = -sin*cos*sigmaY2;
2575   cov[2] = 0.;
2576   cov[3] = cos*cos*sigmaY2;
2577   cov[4] = 0.;
2578   cov[5] = sigmaZ2;
2579   p.SetXYZ(x,y,xyz[2],cov);
2580   AliGeomManager::ELayerID iLayer;
2581   Int_t idet;
2582   if (sector < fkParam->GetNInnerSector()) {
2583     iLayer = AliGeomManager::kTPC1;
2584     idet = sector;
2585   }
2586   else {
2587     iLayer = AliGeomManager::kTPC2;
2588     idet = sector - fkParam->GetNInnerSector();
2589   }
2590   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
2591   p.SetVolumeID(volid);
2592   return kTRUE;
2593 }
2594
2595
2596
2597 Int_t AliTPCtracker::UpdateClusters(AliTPCseed& t,  Int_t nr) {
2598   //-----------------------------------------------------------------
2599   // This function tries to find a track prolongation to next pad row
2600   //-----------------------------------------------------------------
2601   t.SetCurrentCluster(0);
2602   t.SetCurrentClusterIndex1(-3);   
2603    
2604   Double_t xt=t.GetX();
2605   Int_t     row = GetRowNumber(xt)-1; 
2606   Double_t  ymax= GetMaxY(nr);
2607
2608   if (row < nr) return 1; // don't prolongate if not information until now -
2609 //   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
2610 //     t.fRemoval =10;
2611 //     return 0;  // not prolongate strongly inclined tracks
2612 //   } 
2613 //   if (TMath::Abs(t.GetSnp())>0.95) {
2614 //     t.fRemoval =10;
2615 //     return 0;  // not prolongate strongly inclined tracks
2616 //   }// patch 28 fev 06
2617
2618   Double_t x= GetXrow(nr);
2619   Double_t y,z;
2620   //t.PropagateTo(x+0.02);
2621   //t.PropagateTo(x+0.01);
2622   if (!t.PropagateTo(x)){
2623     return 0;
2624   }
2625   //
2626   y=t.GetY();
2627   z=t.GetZ();
2628
2629   if (TMath::Abs(y)>ymax){
2630     if (y > ymax) {
2631       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
2632       if (!t.Rotate(fSectors->GetAlpha())) 
2633         return 0;
2634     } else if (y <-ymax) {
2635       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
2636       if (!t.Rotate(-fSectors->GetAlpha())) 
2637         return 0;
2638     }
2639     //    if (!t.PropagateTo(x)){
2640     //  return 0;
2641     //}
2642     return 1;
2643     //y = t.GetY();    
2644   }
2645   //
2646   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
2647
2648   if (!IsActive(t.GetRelativeSector(),nr)) {
2649     t.SetInDead(kTRUE);
2650     t.SetClusterIndex2(nr,-1); 
2651     return 0;
2652   }
2653   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
2654
2655   AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
2656
2657   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2658     t.SetInDead(kTRUE);
2659     t.SetClusterIndex2(nr,-1); 
2660     return 0;
2661   } 
2662   else
2663     {
2664
2665       //      if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) 
2666       if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
2667       else
2668         return 0;      
2669     }
2670
2671   // update current
2672   if ( (nr%2==0) || t.GetNumberOfClusters()<2){
2673     //    t.fCurrentSigmaY = GetSigmaY(&t);
2674     //t.fCurrentSigmaZ = GetSigmaZ(&t);
2675     GetShape(&t,nr);
2676   }
2677     
2678   AliTPCclusterMI *cl=0;
2679   Int_t index=0;
2680   //
2681   Double_t roady = 1.;
2682   Double_t roadz = 1.;
2683   //
2684
2685   if (!cl){
2686     index = t.GetClusterIndex2(nr);    
2687     if ( (index >= 0) && (index&0x8000)==0){
2688       cl = t.GetClusterPointer(nr);
2689       if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
2690       t.SetCurrentClusterIndex1(index);
2691       if (cl) {
2692         t.SetCurrentCluster(cl);
2693         return 1;
2694       }
2695     }
2696   }
2697
2698   //  if (index<0) return 0;
2699   UInt_t uindex = TMath::Abs(index);
2700
2701   if (krow) {    
2702     //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);      
2703     cl = krow.FindNearest2(y,z,roady,roadz,uindex);      
2704   }
2705
2706   if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));   
2707   t.SetCurrentCluster(cl);
2708
2709   return 1;
2710 }
2711
2712
2713 Int_t AliTPCtracker::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
2714   //-----------------------------------------------------------------
2715   // This function tries to find a track prolongation to next pad row
2716   //-----------------------------------------------------------------
2717
2718   //update error according neighborhoud
2719
2720   if (t.GetCurrentCluster()) {
2721     t.SetRow(nr); 
2722     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
2723     
2724     if (t.GetCurrentCluster()->IsUsed(10)){
2725       //
2726       //
2727       //  t.fErrorZ2*=2;
2728       //  t.fErrorY2*=2;
2729       t.SetNShared(t.GetNShared()+1);
2730       if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
2731         t.SetRemoval(10);
2732         return 0;
2733       }
2734     }   
2735     if (fIteration>0) accept = 0;
2736     if (accept<3)  UpdateTrack(&t,accept);  
2737  
2738   } else {
2739     if (fIteration==0){
2740       if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);      
2741       if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);      
2742
2743       if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
2744     }
2745   }
2746   return 1;
2747 }
2748
2749
2750
2751 //_____________________________________________________________________________
2752 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
2753   //-----------------------------------------------------------------
2754   // This function tries to find a track prolongation.
2755   //-----------------------------------------------------------------
2756   Double_t xt=t.GetX();
2757   //
2758   Double_t alpha=t.GetAlpha();
2759   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
2760   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
2761   //
2762   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2763     
2764   Int_t first = GetRowNumber(xt);
2765   if (!fromSeeds)
2766     first -= step;
2767   if (first < 0)
2768     first = 0;
2769   for (Int_t nr= first; nr>=rf; nr-=step) {
2770     // update kink info
2771     if (t.GetKinkIndexes()[0]>0){
2772       for (Int_t i=0;i<3;i++){
2773         Int_t index = t.GetKinkIndexes()[i];
2774         if (index==0) break;
2775         if (index<0) continue;
2776         //
2777         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2778         if (!kink){
2779           printf("PROBLEM\n");
2780         }
2781         else{
2782           Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2783           if (kinkrow==nr){
2784             AliExternalTrackParam paramd(t);
2785             kink->SetDaughter(paramd);
2786             kink->SetStatus(2,5);
2787             kink->Update();
2788           }
2789         }
2790       }
2791     }
2792
2793     if (nr==80) t.UpdateReference();
2794     if (nr<fInnerSec->GetNRows()) 
2795       fSectors = fInnerSec;
2796     else
2797       fSectors = fOuterSec;
2798     if (FollowToNext(t,nr)==0) 
2799       if (!t.IsActive()) 
2800         return 0;
2801     
2802   }   
2803   return 1;
2804 }
2805
2806
2807
2808
2809
2810
2811 Int_t AliTPCtracker::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2812   //-----------------------------------------------------------------
2813   // This function tries to find a track prolongation.
2814   //-----------------------------------------------------------------
2815   //
2816   Double_t xt=t.GetX();  
2817   Double_t alpha=t.GetAlpha();
2818   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
2819   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
2820   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2821     
2822   Int_t first = t.GetFirstPoint();
2823   Int_t ri = GetRowNumber(xt); 
2824   if (!fromSeeds)
2825     ri += 1;
2826
2827   if (first<ri) first = ri;
2828   //
2829   if (first<0) first=0;
2830   for (Int_t nr=first; nr<=rf; nr++) {
2831     //    if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2832     if (t.GetKinkIndexes()[0]<0){
2833       for (Int_t i=0;i<3;i++){
2834         Int_t index = t.GetKinkIndexes()[i];
2835         if (index==0) break;
2836         if (index>0) continue;
2837         index = TMath::Abs(index);
2838         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2839         if (!kink){
2840           printf("PROBLEM\n");
2841         }
2842         else{
2843           Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2844           if (kinkrow==nr){
2845             AliExternalTrackParam paramm(t);
2846             kink->SetMother(paramm);
2847             kink->SetStatus(2,1);
2848             kink->Update();
2849           }
2850         }
2851       }      
2852     }
2853     //
2854     if (nr<fInnerSec->GetNRows()) 
2855       fSectors = fInnerSec;
2856     else
2857       fSectors = fOuterSec;
2858
2859     FollowToNext(t,nr);                                                             
2860   }   
2861   return 1;
2862 }
2863
2864
2865
2866
2867    
2868 Float_t AliTPCtracker::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2869 {
2870   // overlapping factor
2871   //
2872   sum1=0;
2873   sum2=0;
2874   Int_t sum=0;
2875   //
2876   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2877   dz2*=dz2;  
2878
2879   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2880   dy2*=dy2;
2881   Float_t distance = TMath::Sqrt(dz2+dy2);
2882   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2883  
2884   //  Int_t offset =0;
2885   Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2886   Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2887   if (lastpoint>160) 
2888     lastpoint =160;
2889   if (firstpoint<0) 
2890     firstpoint = 0;
2891   if (firstpoint>lastpoint) {
2892     firstpoint =lastpoint;
2893     //    lastpoint  =160;
2894   }
2895     
2896   
2897   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2898     if (s1->GetClusterIndex2(i)>0) sum1++;
2899     if (s2->GetClusterIndex2(i)>0) sum2++;
2900     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2901       sum++;
2902     }
2903   }
2904   if (sum<5) return 0;
2905
2906   Float_t summin = TMath::Min(sum1+1,sum2+1);
2907   Float_t ratio = (sum+1)/Float_t(summin);
2908   return ratio;
2909 }
2910
2911 void  AliTPCtracker::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2912 {
2913   // shared clusters
2914   //
2915   Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2916   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2917   Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2918   Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2919   
2920   //
2921   Int_t sumshared=0;
2922   //
2923   //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2924   //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2925   Int_t firstpoint = 0;
2926   Int_t lastpoint = 160;
2927   //
2928   //  if (firstpoint>=lastpoint-5) return;;
2929
2930   for (Int_t i=firstpoint;i<lastpoint;i++){
2931     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2932     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2933       sumshared++;
2934     }
2935   }
2936   if (sumshared>cutN0){
2937     // sign clusters
2938     //
2939     for (Int_t i=firstpoint;i<lastpoint;i++){
2940       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2941       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2942         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2943         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);; 
2944         if (s1->IsActive()&&s2->IsActive()){
2945           p1->SetShared(kTRUE);
2946           p2->SetShared(kTRUE);
2947         }       
2948       }
2949     }
2950   }
2951   //  
2952   if (sumshared>cutN0){
2953     for (Int_t i=0;i<4;i++){
2954       if (s1->GetOverlapLabel(3*i)==0){
2955         s1->SetOverlapLabel(3*i,  s2->GetLabel());
2956         s1->SetOverlapLabel(3*i+1,sumshared);
2957         s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2958         break;
2959       } 
2960     }
2961     for (Int_t i=0;i<4;i++){
2962       if (s2->GetOverlapLabel(3*i)==0){
2963         s2->SetOverlapLabel(3*i,  s1->GetLabel());
2964         s2->SetOverlapLabel(3*i+1,sumshared);
2965         s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2966         break;
2967       } 
2968     }    
2969   }
2970 }
2971
2972 void  AliTPCtracker::SignShared(TObjArray * arr)
2973 {
2974   //
2975   //sort trackss according sectors
2976   //  
2977   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2978     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2979     if (!pt) continue;
2980     //if (pt) RotateToLocal(pt);
2981     pt->SetSort(0);
2982   }
2983   arr->UnSort();
2984   arr->Sort();  // sorting according relative sectors
2985   arr->Expand(arr->GetEntries());
2986   //
2987   //
2988   Int_t nseed=arr->GetEntriesFast();
2989   for (Int_t i=0; i<nseed; i++) {
2990     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2991     if (!pt) continue;
2992     for (Int_t j=0;j<12;j++){
2993       pt->SetOverlapLabel(j,0);
2994     }
2995   }
2996   for (Int_t i=0; i<nseed; i++) {
2997     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2998     if (!pt) continue;
2999     if (pt->GetRemoval()>10) continue;
3000     for (Int_t j=i+1; j<nseed; j++){
3001       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
3002       if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
3003       //      if (pt2){
3004       if (pt2->GetRemoval()<=10) {
3005         //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
3006         SignShared(pt,pt2);
3007       }
3008     }  
3009   }
3010 }
3011
3012
3013 void AliTPCtracker::SortTracks(TObjArray * arr, Int_t mode) const
3014 {
3015   //
3016   //sort tracks in array according mode criteria
3017   Int_t nseed = arr->GetEntriesFast();    
3018   for (Int_t i=0; i<nseed; i++) {
3019     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3020     if (!pt) {
3021       continue;
3022     }
3023     pt->SetSort(mode);
3024   }
3025   arr->UnSort();
3026   arr->Sort();
3027 }
3028
3029
3030 void AliTPCtracker::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
3031 {
3032   //
3033   // Loop over all tracks and remove overlaped tracks (with lower quality)
3034   // Algorithm:
3035   //    1. Unsign clusters
3036   //    2. Sort tracks according quality
3037   //       Quality is defined by the number of cluster between first and last points 
3038   //       
3039   //    3. Loop over tracks - decreasing quality order
3040   //       a.) remove - If the fraction of shared cluster less than factor (1- n or 2) 
3041   //       b.) remove - If the minimal number of clusters less than minimal and not ITS
3042   //       c.) if track accepted - sign clusters
3043   //
3044   //Called in - AliTPCtracker::Clusters2Tracks()
3045   //          - AliTPCtracker::PropagateBack()
3046   //          - AliTPCtracker::RefitInward()
3047   //
3048   // Arguments:
3049   //   factor1 - factor for constrained
3050   //   factor2 -        for non constrained tracks 
3051   //            if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
3052   //
3053   UnsignClusters();
3054   //
3055   Int_t nseed = arr->GetEntriesFast();  
3056   Float_t * quality = new Float_t[nseed];
3057   Int_t   * indexes = new Int_t[nseed];
3058   Int_t good =0;
3059   //
3060   //
3061   for (Int_t i=0; i<nseed; i++) {
3062     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3063     if (!pt){
3064       quality[i]=-1;
3065       continue;
3066     }
3067     pt->UpdatePoints();    //select first last max dens points
3068     Float_t * points = pt->GetPoints();
3069     if (points[3]<0.8) quality[i] =-1;
3070     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
3071     //prefer high momenta tracks if overlaps
3072     quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5); 
3073   }
3074   TMath::Sort(nseed,quality,indexes);
3075   //
3076   //
3077   for (Int_t itrack=0; itrack<nseed; itrack++) {
3078     Int_t trackindex = indexes[itrack];
3079     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);   
3080     if (!pt) continue;
3081     //
3082     if (quality[trackindex]<0){
3083       MarkSeedFree( arr->RemoveAt(trackindex) );
3084       continue;
3085     }
3086     //
3087     //
3088     Int_t first = Int_t(pt->GetPoints()[0]);
3089     Int_t last  = Int_t(pt->GetPoints()[2]);
3090     Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
3091     //
3092     Int_t found,foundable,shared;
3093     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
3094     //    pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
3095     Bool_t itsgold =kFALSE;
3096     if (pt->GetESD()){
3097       Int_t dummy[12];
3098       if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
3099     }
3100     if (!itsgold){
3101       //
3102       if (Float_t(shared+1)/Float_t(found+1)>factor){
3103         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
3104         if( (AliTPCReconstructor::StreamLevel()&kStreamRemoveUsed)>0){ // flag:stream  information about TPC tracks which were descarded (double track removal)
3105           TTreeSRedirector &cstream = *fDebugStreamer;
3106           cstream<<"RemoveUsed"<<
3107             "iter="<<fIteration<<
3108             "pt.="<<pt<<
3109             "\n";
3110         }
3111         MarkSeedFree( arr->RemoveAt(trackindex) );
3112         continue;
3113       }      
3114       if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
3115         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
3116         if( (AliTPCReconstructor::StreamLevel()&kStreamRemoveShort)>0){ // flag:stream  information about TPC tracks which were discarded (short track removal)
3117           TTreeSRedirector &cstream = *fDebugStreamer;
3118           cstream<<"RemoveShort"<<
3119             "iter="<<fIteration<<
3120             "pt.="<<pt<<
3121             "\n";
3122         }
3123         MarkSeedFree( arr->RemoveAt(trackindex) );
3124         continue;
3125       }
3126     }
3127
3128     good++;
3129     //if (sharedfactor>0.4) continue;
3130     if (pt->GetKinkIndexes()[0]>0) continue;
3131     //Remove tracks with undefined properties - seems
3132     if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ? 
3133     //
3134     for (Int_t i=first; i<last; i++) {
3135       Int_t index=pt->GetClusterIndex2(i);
3136       // if (index<0 || index&0x8000 ) continue;
3137       if (index<0 || index&0x8000 ) continue;
3138       AliTPCclusterMI *c= pt->GetClusterPointer(i);        
3139       if (!c) continue;
3140       c->Use(10);  
3141     }    
3142   }
3143   fNtracks = good;
3144   if (fDebug>0){
3145     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
3146   }
3147   delete []quality;
3148   delete []indexes;
3149 }
3150
3151 void AliTPCtracker::DumpClusters(Int_t iter, TObjArray *trackArray) 
3152 {
3153   //
3154   // Dump clusters after reco
3155   // signed and unsigned cluster can be visualized   
3156   // 1. Unsign all cluster
3157   // 2. Sign all used clusters
3158   // 3. Dump clusters
3159   UnsignClusters();
3160   Int_t nseed = trackArray->GetEntries();
3161   for (Int_t i=0; i<nseed; i++){
3162     AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);    
3163     if (!pt) {
3164       continue;
3165     }    
3166     Bool_t isKink=pt->GetKinkIndex(0)!=0;
3167     for (Int_t j=0; j<160; ++j) {
3168       Int_t index=pt->GetClusterIndex2(j);
3169       if (index<0) continue;
3170       AliTPCclusterMI *c= pt->GetClusterPointer(j);
3171       if (!c) continue;
3172       if (isKink) c->Use(100);   // kink
3173       c->Use(10);                // by default usage 10
3174     }
3175   }
3176   //
3177   Int_t eventNr = fEvent->GetEventNumberInFile();
3178
3179   for (Int_t sec=0;sec<fkNIS;sec++){
3180     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
3181       TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
3182       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){    
3183         AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
3184         Float_t gx[3];  cl->GetGlobalXYZ(gx);
3185         (*fDebugStreamer)<<"clDump"<< 
3186           "eventNr="<<eventNr<<
3187           "iter="<<iter<<
3188           "cl.="<<cl<<      
3189           "gx0="<<gx[0]<<
3190           "gx1="<<gx[1]<<
3191           "gx2="<<gx[2]<<
3192           "\n";
3193       }
3194       cla = fInnerSec[sec][row].GetClusters2();
3195       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
3196         AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
3197         Float_t gx[3];  cl->GetGlobalXYZ(gx);
3198         (*fDebugStreamer)<<"clDump"<< 
3199           "eventNr="<<eventNr<<
3200           "iter="<<iter<<
3201           "cl.="<<cl<<
3202           "gx0="<<gx[0]<<
3203           "gx1="<<gx[1]<<
3204           "gx2="<<gx[2]<<
3205           "\n";
3206       }
3207     }
3208   }
3209   
3210   for (Int_t sec=0;sec<fkNOS;sec++){
3211     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
3212       TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
3213       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
3214         Float_t gx[3];  
3215         AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
3216         cl->GetGlobalXYZ(gx);
3217         (*fDebugStreamer)<<"clDump"<< 
3218           "eventNr="<<eventNr<<
3219           "iter="<<iter<<
3220           "cl.="<<cl<<
3221           "gx0="<<gx[0]<<
3222           "gx1="<<gx[1]<<
3223           "gx2="<<gx[2]<<
3224           "\n";      
3225       }
3226       cla = fOuterSec[sec][row].GetClusters2();
3227       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
3228         Float_t gx[3];  
3229         AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
3230         cl->GetGlobalXYZ(gx);
3231         (*fDebugStreamer)<<"clDump"<< 
3232           "eventNr="<<eventNr<<
3233           "iter="<<iter<<
3234           "cl.="<<cl<<
3235           "gx0="<<gx[0]<<
3236           "gx1="<<gx[1]<<
3237           "gx2="<<gx[2]<<
3238           "\n";      
3239       }
3240     }
3241   }
3242   
3243 }
3244 void AliTPCtracker::UnsignClusters() 
3245 {
3246   //
3247   // loop over all clusters and unsign them
3248   //
3249   
3250   for (Int_t sec=0;sec<fkNIS;sec++){
3251     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
3252       TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
3253       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
3254         //      if (cl[icl].IsUsed(10))         
3255         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
3256       cla = fInnerSec[sec][row].GetClusters2();
3257       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
3258         //if (cl[icl].IsUsed(10))       
3259         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
3260     }
3261   }
3262   
3263   for (Int_t sec=0;sec<fkNOS;sec++){
3264     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
3265       TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
3266       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
3267         //if (cl[icl].IsUsed(10))       
3268         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
3269       cla = fOuterSec[sec][row].GetClusters2();
3270       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
3271         //if (cl[icl].IsUsed(10))       
3272         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
3273     }
3274   }
3275   
3276 }
3277
3278
3279
3280 void AliTPCtracker::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
3281 {
3282   //
3283   //sign clusters to be "used"
3284   //
3285   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
3286   // loop over "primaries"
3287   
3288   Float_t sumdens=0;
3289   Float_t sumdens2=0;
3290   Float_t sumn   =0;
3291   Float_t sumn2  =0;
3292   Float_t sumchi =0;
3293   Float_t sumchi2 =0;
3294
3295   Float_t sum    =0;
3296
3297   TStopwatch timer;
3298   timer.Start();
3299
3300   Int_t nseed = arr->GetEntriesFast();
3301   for (Int_t i=0; i<nseed; i++) {
3302     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3303     if (!pt) {
3304       continue;
3305     }    
3306     if (!(pt->IsActive())) continue;
3307     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
3308     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
3309       sumdens += dens;
3310       sumdens2+= dens*dens;
3311       sumn    += pt->GetNumberOfClusters();
3312       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
3313       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
3314       if (chi2>5) chi2=5;
3315       sumchi  +=chi2;
3316       sumchi2 +=chi2*chi2;
3317       sum++;
3318     }
3319   }
3320
3321   Float_t mdensity = 0.9;
3322   Float_t meann    = 130;
3323   Float_t meanchi  = 1;
3324   Float_t sdensity = 0.1;
3325   Float_t smeann    = 10;
3326   Float_t smeanchi  =0.4;
3327   
3328
3329   if (sum>20){
3330     mdensity = sumdens/sum;
3331     meann    = sumn/sum;
3332     meanchi  = sumchi/sum;
3333     //
3334     sdensity = sumdens2/sum-mdensity*mdensity;
3335     if (sdensity >= 0)
3336        sdensity = TMath::Sqrt(sdensity);
3337     else
3338        sdensity = 0.1;
3339     //
3340     smeann   = sumn2/sum-meann*meann;
3341     if (smeann >= 0)
3342       smeann   = TMath::Sqrt(smeann);
3343     else 
3344       smeann   = 10;
3345     //
3346     smeanchi = sumchi2/sum - meanchi*meanchi;
3347     if (smeanchi >= 0)
3348       smeanchi = TMath::Sqrt(smeanchi);
3349     else
3350       smeanchi = 0.4;
3351   }
3352
3353
3354   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
3355   //
3356   for (Int_t i=0; i<nseed; i++) {
3357     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3358     if (!pt) {
3359       continue;
3360     }
3361     if (pt->GetBSigned()) continue;
3362     if (pt->GetBConstrain()) continue;    
3363     //if (!(pt->IsActive())) continue;
3364     /*
3365     Int_t found,foundable,shared;    
3366     pt->GetClusterStatistic(0,160,found, foundable,shared);
3367     if (shared/float(found)>0.3) {
3368       if (shared/float(found)>0.9 ){
3369         //MarkSeedFree( arr->RemoveAt(i) );
3370       }
3371       continue;
3372     }
3373     */
3374     Bool_t isok =kFALSE;
3375     if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
3376       isok = kTRUE;
3377     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
3378       isok =kTRUE;
3379     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
3380       isok =kTRUE;
3381     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
3382       isok =kTRUE;
3383     
3384     if (isok)     
3385       for (Int_t j=0; j<160; ++j) {     
3386         Int_t index=pt->GetClusterIndex2(j);
3387         if (index<0) continue;
3388         AliTPCclusterMI *c= pt->GetClusterPointer(j);
3389         if (!c) continue;
3390         //if (!(c->IsUsed(10))) c->Use();  
3391         c->Use(10);  
3392       }
3393   }
3394   
3395   
3396   //
3397   Double_t maxchi  = meanchi+2.*smeanchi;
3398
3399   for (Int_t i=0; i<nseed; i++) {
3400     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3401     if (!pt) {
3402       continue;
3403     }    
3404     //if (!(pt->IsActive())) continue;
3405     if (pt->GetBSigned()) continue;
3406     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
3407     if (chi>maxchi) continue;
3408
3409     Float_t bfactor=1;
3410     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
3411    
3412     //sign only tracks with enoug big density at the beginning
3413     
3414     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue; 
3415     
3416     
3417     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
3418     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
3419    
3420     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
3421     if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
3422       minn=0;
3423
3424     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
3425       //Int_t noc=pt->GetNumberOfClusters();
3426       pt->SetBSigned(kTRUE);
3427       for (Int_t j=0; j<160; ++j) {
3428
3429         Int_t index=pt->GetClusterIndex2(j);
3430         if (index<0) continue;
3431         AliTPCclusterMI *c= pt->GetClusterPointer(j);
3432         if (!c) continue;
3433         //      if (!(c->IsUsed(10))) c->Use();  
3434         c->Use(10);  
3435       }
3436     }
3437   }
3438   //  gLastCheck = nseed;
3439   //  arr->Compress();
3440   if (fDebug>0){
3441     timer.Print();
3442   }
3443 }
3444
3445
3446
3447 Int_t AliTPCtracker::RefitInward(AliESDEvent *event)
3448 {
3449   //
3450   // back propagation of ESD tracks
3451   //
3452   //return 0;
3453   if (!event) return 0;
3454   const Int_t kMaxFriendTracks=2000;
3455   fEvent = event;
3456   fEventHLT = 0;
3457   // extract correction object for multiplicity dependence of dEdx
3458   TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
3459
3460   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
3461   if (!transform) {
3462     AliFatal("Tranformations not in RefitInward");
3463     return 0;
3464   }
3465   transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
3466   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
3467   Int_t nContribut = event->GetNumberOfTracks();
3468   TGraphErrors * graphMultDependenceDeDx = 0x0;
3469   if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
3470     if (recoParam->GetUseTotCharge()) {
3471       graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
3472     } else {
3473       graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
3474     }
3475   }
3476   //
3477   ReadSeeds(event,2);
3478   fIteration=2;
3479   //PrepareForProlongation(fSeeds,1);
3480   PropagateForward2(fSeeds);
3481   RemoveUsed2(fSeeds,0.4,0.4,20);
3482
3483   Int_t entriesSeed=fSeeds->GetEntries();
3484   TObjArray arraySeed(entriesSeed);
3485   for (Int_t i=0;i<entriesSeed;i++) {
3486     arraySeed.AddAt(fSeeds->At(i),i);    
3487   }
3488   SignShared(&arraySeed);
3489   //  FindCurling(fSeeds, event,2); // find multi found tracks
3490   FindSplitted(fSeeds, event,2); // find multi found tracks
3491   if ((AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC)>0)  FindMultiMC(fSeeds, fEvent,2); // flag: stream MC infomation about the multiple find track (ONLY for MC data)
3492
3493   Int_t ntracks=0;
3494   Int_t nseed = fSeeds->GetEntriesFast();
3495   for (Int_t i=0;i<nseed;i++){
3496     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
3497     if (!seed) continue;
3498     if (seed->GetKinkIndex(0)>0)  UpdateKinkQualityD(seed);  // update quality informations for kinks
3499     AliESDtrack *esd=event->GetTrack(i);
3500
3501     if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
3502       AliExternalTrackParam paramIn;
3503       AliExternalTrackParam paramOut;
3504       Int_t ncl = seed->RefitTrack(seed,&paramIn,&paramOut);
3505       if ((AliTPCReconstructor::StreamLevel() & kStreamRecoverIn)>0) {  // flag: stream track information for track  failing in RefitInward function and recovered back
3506         (*fDebugStreamer)<<"RecoverIn"<< 
3507           "seed.="<<seed<<
3508           "esd.="<<esd<<
3509           "pin.="<<&paramIn<<
3510           "pout.="<<&paramOut<<
3511           "ncl="<<ncl<<
3512           "\n";
3513       }
3514       if (ncl>15) {
3515         seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
3516         seed->SetNumberOfClusters(ncl);
3517       }
3518     }
3519
3520     seed->PropagateTo(fkParam->GetInnerRadiusLow());
3521     seed->UpdatePoints();
3522     AddCovariance(seed);
3523     MakeESDBitmaps(seed, esd);
3524     seed->CookdEdx(0.02,0.6);
3525     CookLabel(seed,0.1); //For comparison only
3526     //
3527     if (((AliTPCReconstructor::StreamLevel()&kStreamRefitInward)>0) && seed!=0) {
3528       TTreeSRedirector &cstream = *fDebugStreamer;
3529       cstream<<"RefitInward"<<  // flag: stream track information in RefitInward function (after tracking Iteration 2)
3530         "Esd.="<<esd<<
3531         "Track.="<<seed<<
3532         "\n"; 
3533     }
3534
3535     if (seed->GetNumberOfClusters()>15){
3536       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit); 
3537       esd->SetTPCPoints(seed->GetPoints());
3538       esd->SetTPCPointsF(seed->GetNFoundable());
3539       Int_t   ndedx = seed->GetNCDEDX(0);
3540       Float_t sdedx = seed->GetSDEDX(0);
3541       Float_t dedx  = seed->GetdEdx();
3542       // apply mutliplicity dependent dEdx correction if available
3543       if (graphMultDependenceDeDx) {
3544         Double_t corrGain =  AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
3545         dedx += (1 - corrGain)*50.; // MIP is normalized to 50
3546       }
3547       esd->SetTPCsignal(dedx, sdedx, ndedx);
3548       //
3549       // fill new dEdx information
3550       //
3551       Double32_t signal[4]; 
3552       Double32_t signalMax[4]; 
3553       Char_t ncl[3]; 
3554       Char_t nrows[3];
3555       //
3556       for(Int_t iarr=0;iarr<3;iarr++) {
3557         signal[iarr] = seed->GetDEDXregion(iarr+1);
3558         signalMax[iarr] = seed->GetDEDXregion(iarr+5);
3559         ncl[iarr] = seed->GetNCDEDX(iarr+1);
3560         nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
3561       }
3562       signal[3] = seed->GetDEDXregion(4);
3563       signalMax[3] = seed->GetDEDXregion(8);
3564       
3565       //
3566       AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
3567       infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
3568       infoTpcPid->SetTPCSignalsQmax(signalMax);
3569       esd->SetTPCdEdxInfo(infoTpcPid);
3570       //
3571       // add seed to the esd track in Calib level
3572       //
3573       Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
3574       if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
3575         // RS: this is the only place where the seed is created not in the pool, 
3576         // since it should belong to ESDevent
3577         AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE); 
3578         esd->AddCalibObject(seedCopy);
3579       }
3580       ntracks++;
3581     }
3582     else{
3583       //printf("problem\n");
3584     }
3585   }
3586   //FindKinks(fSeeds,event);
3587   if ((AliTPCReconstructor::StreamLevel()&kStreamClDump)>0)  DumpClusters(2,fSeeds);  // dump clusters at the end of process (signed with useage flags)
3588   Info("RefitInward","Number of refitted tracks %d",ntracks);
3589
3590   AliCosmicTracker::FindCosmic(event, kTRUE);
3591
3592   FillClusterOccupancyInfo();
3593
3594   return 0;
3595 }
3596
3597
3598 Int_t AliTPCtracker::PropagateBack(AliESDEvent *event)
3599 {
3600   //
3601   // back propagation of ESD tracks
3602   //
3603   if (!event) return 0;
3604   fEvent = event;
3605   fEventHLT = 0;
3606   fIteration = 1;
3607   ReadSeeds(event,1);
3608   PropagateBack(fSeeds); 
3609   RemoveUsed2(fSeeds,0.4,0.4,20);
3610   //FindCurling(fSeeds, fEvent,1);  
3611   FindSplitted(fSeeds, event,1); // find multi found tracks
3612   if ((AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC)>0)  FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
3613
3614   //
3615   Int_t nseed = fSeeds->GetEntriesFast();
3616   Int_t ntracks=0;
3617   for (Int_t i=0;i<nseed;i++){
3618     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
3619     if (!seed) continue;
3620     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
3621     seed->UpdatePoints();
3622     AddCovariance(seed);
3623     AliESDtrack *esd=event->GetTrack(i);
3624     if (!esd) continue; //never happen
3625     if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
3626       AliExternalTrackParam paramIn;
3627       AliExternalTrackParam paramOut;
3628       Int_t ncl = seed->RefitTrack(seed,&paramIn,&paramOut);
3629       if ((AliTPCReconstructor::StreamLevel()&kStreamRecoverBack)>0) { // flag: stream track information for track  faling PropagateBack function and recovered back (RefitTrack)
3630         (*fDebugStreamer)<<"RecoverBack"<<
3631           "seed.="<<seed<<
3632           "esd.="<<esd<<
3633           "pin.="<<&paramIn<<
3634           "pout.="<<&paramOut<<
3635           "ncl="<<ncl<<
3636           "\n";
3637       }
3638       if (ncl>15) {
3639         seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
3640         seed->SetNumberOfClusters(ncl);
3641       }
3642     }
3643     seed->CookdEdx(0.02,0.6);
3644     CookLabel(seed,0.1); //For comparison only
3645     if (seed->GetNumberOfClusters()>15){
3646       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
3647       esd->SetTPCPoints(seed->GetPoints());
3648       esd->SetTPCPointsF(seed->GetNFoundable());
3649       Int_t   ndedx = seed->GetNCDEDX(0);
3650       Float_t sdedx = seed->GetSDEDX(0);
3651       Float_t dedx  = seed->GetdEdx();
3652       esd->SetTPCsignal(dedx, sdedx, ndedx);
3653       ntracks++;
3654       Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
3655       // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number      
3656       if (((AliTPCReconstructor::StreamLevel()&kStreamCPropagateBack)>0) && esd) {
3657         (*fDebugStreamer)<<"PropagateBack"<<  // flag: stream track information in PropagateBack function (after tracking Iteration 1)
3658           "Tr0.="<<seed<<
3659           "esd.="<<esd<<
3660           "EventNrInFile="<<eventnumber<<
3661           "\n";       
3662       }
3663     }
3664   }
3665   if (AliTPCReconstructor::StreamLevel()&kStreamClDump)  DumpClusters(1,fSeeds); // flag:stream clusters at the end of process (signed with usage flag)
3666   //FindKinks(fSeeds,event);
3667   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
3668   fEvent =0;
3669   fEventHLT = 0;
3670
3671   return 0;
3672 }
3673
3674
3675 Int_t AliTPCtracker::PostProcess(AliESDEvent *event)
3676 {
3677   //
3678   // Post process events 
3679   //
3680   if (!event) return 0;
3681
3682   //
3683   // Set TPC event status
3684   // 
3685
3686   // event affected by HV dip
3687   // reset TPC status
3688   if(IsTPCHVDipEvent(event)) { 
3689     event->ResetDetectorStatus(AliDAQ::kTPC);
3690   }
3691  
3692   //printf("Status %d \n", event->IsDetectorOn(AliDAQ::kTPC));
3693
3694   return 0;
3695 }
3696
3697
3698  void AliTPCtracker::DeleteSeeds()
3699 {
3700   //
3701   fSeeds->Clear();
3702   ResetSeedsPool();
3703   delete fSeeds;
3704   fSeeds =0;
3705 }
3706
3707 void AliTPCtracker::ReadSeeds(const AliESDEvent *const event, Int_t direction)
3708 {
3709   //
3710   //read seeds from the event
3711   
3712   Int_t nentr=event->GetNumberOfTracks();
3713   if (fDebug>0){
3714     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
3715   }
3716   if (fSeeds) 
3717     DeleteSeeds();
3718   if (!fSeeds){   
3719     fSeeds = new TObjArray(nentr);
3720   }
3721   UnsignClusters();
3722   //  Int_t ntrk=0;  
3723   for (Int_t i=0; i<nentr; i++) {
3724     AliESDtrack *esd=event->GetTrack(i);
3725     ULong_t status=esd->GetStatus();
3726     if (!(status&AliESDtrack::kTPCin)) continue;
3727     AliTPCtrack t(*esd);
3728     t.SetNumberOfClusters(0);
3729     //    AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
3730     AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t/*,t.GetAlpha()*/);
3731     seed->SetPoolID(fLastSeedID);
3732     seed->SetUniqueID(esd->GetID());
3733     AddCovariance(seed);   //add systematic ucertainty
3734     for (Int_t ikink=0;ikink<3;ikink++) {
3735       Int_t index = esd->GetKinkIndex(ikink);
3736       seed->GetKinkIndexes()[ikink] = index;
3737       if (index==0) continue;
3738       index = TMath::Abs(index);
3739       AliESDkink * kink = fEvent->GetKink(index-1);
3740       if (kink&&esd->GetKinkIndex(ikink)<0){
3741         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
3742         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,0);
3743       }
3744       if (kink&&esd->GetKinkIndex(ikink)>0){
3745         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
3746         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,4);
3747       }
3748
3749     }
3750     if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.); 
3751     //RS    if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
3752     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
3753     //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
3754     //  fSeeds->AddAt(0,i);
3755     //  MarkSeedFree( seed );
3756     //  continue;    
3757     //}
3758     
3759     //
3760     //
3761     // rotate to the local coordinate system
3762     //   
3763     fSectors=fInnerSec; fN=fkNIS;    
3764     Double_t alpha=seed->GetAlpha();
3765     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
3766     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
3767     Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
3768     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
3769     alpha-=seed->GetAlpha();  
3770     if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
3771     if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
3772     if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
3773       AliWarning(Form("Rotating track over %f",alpha));
3774       if (!seed->Rotate(alpha)) {
3775         MarkSeedFree( seed );
3776         continue;
3777       }
3778     }
3779     seed->SetESD(esd);
3780     // sign clusters
3781     if (esd->GetKinkIndex(0)<=0){
3782       for (Int_t irow=0;irow<160;irow++){
3783         Int_t index = seed->GetClusterIndex2(irow);    
3784         if (index >= 0){ 
3785           //
3786           AliTPCclusterMI * cl = GetClusterMI(index);
3787           seed->SetClusterPointer(irow,cl);
3788           if (cl){
3789             if ((index & 0x8000)==0){
3790               cl->Use(10);  // accepted cluster   
3791             }else{
3792               cl->Use(6);   // close cluster not accepted
3793             }   
3794           }else{
3795             Info("ReadSeeds","Not found cluster");
3796           }
3797         }
3798       }
3799     }
3800     fSeeds->AddAt(seed,i);
3801   }
3802 }
3803
3804
3805
3806 //_____________________________________________________________________________
3807 void AliTPCtracker::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3808                                  Float_t deltay, Int_t ddsec) {
3809   //-----------------------------------------------------------------
3810   // This function creates track seeds.
3811   // SEEDING WITH VERTEX CONSTRAIN 
3812   //-----------------------------------------------------------------
3813   // cuts[0]   - fP4 cut
3814   // cuts[1]   - tan(phi)  cut
3815   // cuts[2]   - zvertex cut
3816   // cuts[3]   - fP3 cut
3817   Int_t nin0  = 0;
3818   Int_t nin1  = 0;
3819   Int_t nin2  = 0;
3820   Int_t nin   = 0;
3821   Int_t nout1 = 0;
3822   Int_t nout2 = 0;
3823
3824   Double_t x[5], c[15];
3825   //  Int_t di = i1-i2;
3826   //
3827   AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed();
3828   seed->SetPoolID(fLastSeedID);
3829   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3830   Double_t cs=cos(alpha), sn=sin(alpha);
3831   //
3832   //  Double_t x1 =fOuterSec->GetX(i1);
3833   //Double_t xx2=fOuterSec->GetX(i2);
3834   
3835   Double_t x1 =GetXrow(i1);
3836   Double_t xx2=GetXrow(i2);
3837
3838   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3839
3840   Int_t imiddle = (i2+i1)/2;    //middle pad row index
3841   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3842   const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3843   //
3844   Int_t ns =sec;   
3845
3846   const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3847   Double_t ymax  = GetMaxY(i1)-kr1.GetDeadZone()-1.5;  
3848   Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;  
3849
3850   //
3851   // change cut on curvature if it can't reach this layer
3852   // maximal curvature set to reach it
3853   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3854   if (dvertexmax*0.5*cuts[0]>0.85){
3855     cuts[0] = 0.85/(dvertexmax*0.5+1.);
3856   }
3857   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
3858
3859   //  Int_t ddsec = 1;
3860   if (deltay>0) ddsec = 0; 
3861   // loop over clusters  
3862   for (Int_t is=0; is < kr1; is++) {
3863     //
3864     if (kr1[is]->IsUsed(10)) continue;
3865     if (kr1[is]->IsDisabled()) {
3866       continue;
3867     }
3868
3869     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
3870     //if (TMath::Abs(y1)>ymax) continue;
3871
3872     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
3873
3874     // find possible directions    
3875     Float_t anglez = (z1-z3)/(x1-x3); 
3876     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z      
3877     //
3878     //
3879     //find   rotation angles relative to line given by vertex and point 1
3880     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3881     Double_t dvertex  = TMath::Sqrt(dvertex2);
3882     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
3883     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);            
3884     
3885     //
3886     // loop over 2 sectors
3887     Int_t dsec1=-ddsec;
3888     Int_t dsec2= ddsec;
3889     if (y1<0)  dsec2= 0;
3890     if (y1>0)  dsec1= 0;
3891     
3892     Double_t dddz1=0;  // direction of delta inclination in z axis
3893     Double_t dddz2=0;
3894     if ( (z1-z3)>0)
3895       dddz1 =1;    
3896     else
3897       dddz2 =1;
3898     //
3899     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3900       Int_t sec2 = sec + dsec;
3901       // 
3902       //      AliTPCtrackerRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3903       //AliTPCtrackerRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3904       AliTPCtrackerRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
3905       AliTPCtrackerRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3906       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3907       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3908
3909       // rotation angles to p1-p3
3910       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;            
3911       Double_t x2,   y2,   z2; 
3912       //
3913       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3914
3915       //
3916       Double_t dxx0 =  (xx2-x3)*cs13r;
3917       Double_t dyy0 =  (xx2-x3)*sn13r;
3918       for (Int_t js=index1; js < index2; js++) {
3919         const AliTPCclusterMI *kcl = kr2[js];
3920         if (kcl->IsUsed(10)) continue;  
3921         if (kcl->IsDisabled()) {
3922           continue;
3923         }
3924         //
3925         //calcutate parameters
3926         //      
3927         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
3928         // stright track
3929         if (TMath::Abs(yy0)<0.000001) continue;
3930         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
3931         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3932         Double_t r02 = (0.25+y0*y0)*dvertex2;   
3933         //curvature (radius) cut
3934         if (r02<r2min) continue;                
3935        
3936         nin0++; 
3937         //
3938         Double_t c0  = 1/TMath::Sqrt(r02);
3939         if (yy0>0) c0*=-1.;     
3940                
3941        
3942         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
3943         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3944         Double_t dfi0   = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3945         Double_t dfi1   = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);  
3946         //
3947         //
3948         Double_t z0  =  kcl->GetZ();  
3949         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
3950         if (TMath::Abs(zzzz2-z0)>0.5) continue;       
3951         nin1++;              
3952         //      
3953         Double_t dip    = (z1-z0)*c0/dfi1;        
3954         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3955         //
3956         y2 = kcl->GetY(); 
3957         if (dsec==0){
3958           x2 = xx2; 
3959           z2 = kcl->GetZ();       
3960         }
3961         else
3962           {
3963             // rotation 
3964             z2 = kcl->GetZ();  
3965             x2= xx2*cs-y2*sn*dsec;
3966             y2=+xx2*sn*dsec+y2*cs;
3967           }
3968         
3969         x[0] = y1;
3970         x[1] = z1;
3971         x[2] = x0;
3972         x[3] = dip;
3973         x[4] = c0;
3974         //
3975         //
3976         // do we have cluster at the middle ?
3977         Double_t ym,zm;
3978         GetProlongation(x1,xm,x,ym,zm);
3979         UInt_t dummy; 
3980         AliTPCclusterMI * cm=0;
3981         if (TMath::Abs(ym)-ymaxm<0){      
3982           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3983           if ((!cm) || (cm->IsUsed(10))) {        
3984             continue;
3985           }
3986         }
3987         else{     
3988           // rotate y1 to system 0
3989           // get state vector in rotated system 
3990           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
3991           Double_t xr2  =  x0*cs+yr1*sn*dsec;
3992           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3993           //
3994           GetProlongation(xx2,xm,xr,ym,zm);
3995           if (TMath::Abs(ym)-ymaxm<0){
3996             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3997             if ((!cm) || (cm->IsUsed(10))) {      
3998               continue;
3999             }