PWGPP-71 Move TPC noise filtering information form the AliESDevent to the AliESDHeader
[u/mrichter/AliRoot.git] / TPC / Rec / AliTPCtracker.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 //-------------------------------------------------------
18 //          Implementation of the TPC tracker
19 //
20 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
21 // 
22 //  AliTPC parallel tracker
23 //
24 //  The track fitting is based on Kalman filtering approach
25 //  The track finding steps:
26 //      1. Seeding - with and without vertex constraint
27 //                 - seeding with vertex constain done at first n^2 proble
28 //                 - seeding without vertex constraint n^3 problem
29 //      2. Tracking - follow prolongation road - find cluster - update kalman track
30 //  The seeding and tracking is repeated several times, in different seeding region.
31 //  This approach enables to find the track which cannot be seeded in some region of TPC
32 //  This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
33
34 //  With this approach we reach almost 100 % efficiency also for high occupancy events.
35 //  (If the seeding efficiency in a region is about 90 % than with logical or of several 
36 //  regions we will reach 100% (in theory - supposing independence) 
37
38 //  Repeating several seeding - tracking procedures some of the tracks can be find 
39 //  several times. 
40
41 //  The procedures to remove multi find tacks are impremented:
42 //  RemoveUsed2 - fast procedure n problem - 
43 //                 Algorithm - Sorting tracks according quality
44 //                             remove tracks with some shared fraction 
45 //                             Sharing in respect to all tacks 
46 //                             Signing clusters in gold region
47 //  FindSplitted - slower algorithm n^2
48 //                 Sort the tracks according quality
49 //                 Loop over pair of tracks
50 //                 If overlap with other track bigger than threshold - remove track
51 //  
52 //  FindCurling  - Finds the pair of tracks which are curling
53 //               - About 10% of tracks can be find with this procedure
54 //                 The combinatorial background is too big to be used in High 
55 //                  multiplicity environment 
56 //               - n^2 problem - Slow procedure - currently it is disabled because of 
57 //                  low efficiency
58 //                 
59 //  The number of splitted tracks can be reduced disabling the sharing of the cluster.
60 //  tpcRecoParam-> SetClusterSharing(kFALSE);
61 //  IT IS HIGHLY non recomended to use it in high flux enviroonment
62 //  Even using this switch some tracks can be found more than once 
63 //  (because of multiple seeding and low quality tracks which will not cross full chamber)
64 //                          
65 //
66 // The tracker itself can be debugged  - the information about tracks can be stored in several // phases of the reconstruction
67 // To enable storage of the TPC tracks in the ESD friend track
68 // use AliTPCReconstructor::SetStreamLevel(n); 
69 //
70 // The debug level -  different procedure produce tree for numerical debugging of code and data (see comments foEStreamFlags in AliTPCtracker.h  )
71 //
72
73 //
74 // Adding systematic errors to the covariance:
75 // 
76 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
77 // of the tracks (not to the clusters as they are dependent):
78 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
79 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/GeV)
80 // The default values are 0. 
81 //
82 // The systematic errors are added to the covariance matrix in following places:
83 //
84 // 1. During fisrt itteration - AliTPCtracker::FillESD
85 // 2. Second iteration - 
86 //      2.a ITS->TPC   - AliTPCtracker::ReadSeeds 
87 //      2.b TPC->TRD   - AliTPCtracker::PropagateBack
88 // 3. Third iteration  -
89 //      3.a TRD->TPC   - AliTPCtracker::ReadSeeds
90 //      3.b TPC->ITS   - AliTPCtracker::RefitInward
91 //
92
93 /* $Id$ */
94
95 #include "Riostream.h"
96 #include <TClonesArray.h>
97 #include <TFile.h>
98 #include <TObjArray.h>
99 #include <TTree.h>
100 #include <TMatrixD.h>
101 #include <TGraphErrors.h>
102 #include <TTimeStamp.h>
103 #include "AliLog.h"
104 #include "AliComplexCluster.h"
105 #include "AliESDEvent.h"
106 #include "AliESDtrack.h"
107 #include "AliESDVertex.h"
108 #include "AliKink.h"
109 #include "AliV0.h"
110 #include "AliHelix.h"
111 #include "AliRunLoader.h"
112 #include "AliTPCClustersRow.h"
113 #include "AliTPCParam.h"
114 #include "AliTPCReconstructor.h"
115 #include "AliTPCpolyTrack.h"
116 #include "AliTPCreco.h"
117 #include "AliTPCseed.h"
118
119 #include "AliTPCtrackerSector.h" 
120 #include "AliTPCtracker.h"
121 #include "TStopwatch.h"
122 #include "AliTPCReconstructor.h"
123 #include "AliAlignObj.h"
124 #include "AliTrackPointArray.h"
125 #include "TRandom.h"
126 #include "AliTPCcalibDB.h"
127 #include "AliTPCcalibDButil.h"
128 #include "AliTPCTransform.h"
129 #include "AliTPCClusterParam.h"
130 #include "AliTPCdEdxInfo.h"
131 #include "AliDCSSensorArray.h"
132 #include "AliDCSSensor.h"
133 #include "AliDAQ.h"
134 #include "AliCosmicTracker.h"
135 #include "AliTPCROC.h"
136 #include "AliMathBase.h"
137 //
138
139 using std::cerr;
140 using std::endl;
141 ClassImp(AliTPCtracker)
142
143
144
145 class AliTPCFastMath {
146 public:
147   AliTPCFastMath();  
148   static Double_t FastAsin(Double_t x);   
149  private: 
150   static Double_t fgFastAsin[20000];  //lookup table for fast asin computation
151 };
152
153 Double_t AliTPCFastMath::fgFastAsin[20000];
154 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
155
156 AliTPCFastMath::AliTPCFastMath(){
157   //
158   // initialized lookup table;
159   for (Int_t i=0;i<10000;i++){
160     fgFastAsin[2*i] = TMath::ASin(i/10000.);
161     fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
162   }
163 }
164
165 Double_t AliTPCFastMath::FastAsin(Double_t x){
166   //
167   // return asin using lookup table
168   if (x>0){
169     Int_t index = int(x*10000);
170     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
171   }
172   x*=-1;
173   Int_t index = int(x*10000);
174   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
175 }
176 //__________________________________________________________________
177 AliTPCtracker::AliTPCtracker()
178                 :AliTracker(),
179                  fkNIS(0),
180                  fInnerSec(0),
181                  fkNOS(0),
182                  fOuterSec(0),
183                  fN(0),
184                  fSectors(0),
185                  fInput(0),
186                  fOutput(0),
187                  fSeedTree(0),
188                  fTreeDebug(0),
189                  fEvent(0),
190                  fEventHLT(0),
191                  fDebug(0),
192                  fNewIO(kFALSE),
193                  fNtracks(0),
194                  fSeeds(0),
195                  fIteration(0),
196                  fkParam(0),
197                  fDebugStreamer(0),
198                  fUseHLTClusters(4),
199          fCrossTalkSignalArray(0),
200                  fSeedsPool(0),
201                  fFreeSeedsID(500),
202                  fNFreeSeeds(0),
203                  fLastSeedID(-1)
204 {
205   //
206   // default constructor
207   //
208   for (Int_t irow=0; irow<200; irow++){
209     fXRow[irow]=0;
210     fYMax[irow]=0;
211     fPadLength[irow]=0;
212   }
213 }
214 //_____________________________________________________________________
215
216
217
218 Int_t AliTPCtracker::UpdateTrack(AliTPCseed * track, Int_t accept){
219   //
220   //update track information using current cluster - track->fCurrentCluster
221
222
223   AliTPCclusterMI* c =track->GetCurrentCluster();
224   if (accept > 0) //sign not accepted clusters
225     track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);  
226   else // unsign accpeted clusters
227     track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);  
228   UInt_t i = track->GetCurrentClusterIndex1();
229
230   Int_t sec=(i&0xff000000)>>24; 
231   //Int_t row = (i&0x00ff0000)>>16; 
232   track->SetRow((i&0x00ff0000)>>16);
233   track->SetSector(sec);
234   //  Int_t index = i&0xFFFF;
235   if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow()); 
236   track->SetClusterIndex2(track->GetRow(), i);  
237   //track->fFirstPoint = row;
238   //if ( track->fLastPoint<row) track->fLastPoint =row;
239   //  if (track->fRow<0 || track->fRow>160) {
240   //  printf("problem\n");
241   //}
242   if (track->GetFirstPoint()>track->GetRow()) 
243     track->SetFirstPoint(track->GetRow());
244   if (track->GetLastPoint()<track->GetRow()) 
245     track->SetLastPoint(track->GetRow());
246   
247
248   track->SetClusterPointer(track->GetRow(),c);  
249   //
250
251   Double_t angle2 = track->GetSnp()*track->GetSnp();
252   //
253   //SET NEW Track Point
254   //
255   if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
256   {
257     angle2 = TMath::Sqrt(angle2/(1-angle2)); 
258     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->GetRow()));
259     //
260     point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
261     point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
262     point.SetErrY(sqrt(track->GetErrorY2()));
263     point.SetErrZ(sqrt(track->GetErrorZ2()));
264     //
265     point.SetX(track->GetX());
266     point.SetY(track->GetY());
267     point.SetZ(track->GetZ());
268     point.SetAngleY(angle2);
269     point.SetAngleZ(track->GetTgl());
270     if (point.IsShared()){
271       track->SetErrorY2(track->GetErrorY2()*4);
272       track->SetErrorZ2(track->GetErrorZ2()*4);
273     }
274   }  
275
276   Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
277   //
278 //   track->SetErrorY2(track->GetErrorY2()*1.3);
279 //   track->SetErrorY2(track->GetErrorY2()+0.01);    
280 //   track->SetErrorZ2(track->GetErrorZ2()*1.3);   
281 //   track->SetErrorZ2(track->GetErrorZ2()+0.005);      
282     //}
283   if (accept>0) return 0;
284   if (track->GetNumberOfClusters()%20==0){
285     //    if (track->fHelixIn){
286     //  TClonesArray & larr = *(track->fHelixIn);    
287     //  Int_t ihelix = larr.GetEntriesFast();
288     //  new(larr[ihelix]) AliHelix(*track) ;    
289     //}
290   }
291   if (AliTPCReconstructor::StreamLevel()&kStreamUpdateTrack) {
292     Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
293     AliExternalTrackParam param(*track);
294     TTreeSRedirector &cstream = *fDebugStreamer;
295     cstream<<"UpdateTrack"<<
296       "cl.="<<c<<
297       "event="<<event<<
298       "track.="<<&param<<
299       "\n";
300   }
301   track->SetNoCluster(0);
302   return track->Update(c,chi2,i);
303 }
304
305
306
307 Int_t AliTPCtracker::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
308 {
309   //
310   // decide according desired precision to accept given 
311   // cluster for tracking
312   Double_t  yt=0,zt=0;
313   seed->GetProlongation(cluster->GetX(),yt,zt);
314   Double_t sy2=ErrY2(seed,cluster);
315   Double_t sz2=ErrZ2(seed,cluster);
316   
317   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
318   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
319   Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
320   Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
321   Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
322     (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
323   Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
324     (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
325   
326   Double_t rdistance2  = rdistancey2+rdistancez2;
327   //Int_t  accept =0;
328   
329   if (AliTPCReconstructor::StreamLevel()>2 && ( (fIteration>0)|| (seed->GetNumberOfClusters()>20))) {
330     //  if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
331     Float_t rmsy2 = seed->GetCurrentSigmaY2();
332     Float_t rmsz2 = seed->GetCurrentSigmaZ2();
333     Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
334     Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
335     Float_t rmsy2p30R  = seed->GetCMeanSigmaY2p30R();
336     Float_t rmsz2p30R  = seed->GetCMeanSigmaZ2p30R();
337     AliExternalTrackParam param(*seed);
338     static TVectorD gcl(3),gtr(3);
339     Float_t gclf[3];
340     param.GetXYZ(gcl.GetMatrixArray());
341     cluster->GetGlobalXYZ(gclf);
342     gcl[0]=gclf[0];    gcl[1]=gclf[1];    gcl[2]=gclf[2];
343     Int_t nclSeed=seed->GetNumberOfClusters();
344     
345     if (AliTPCReconstructor::StreamLevel()&kStreamErrParam) { // flag:stream in debug mode cluster and track extrapolation at given row together with error nad shape estimate
346       Int_t eventNr = fEvent->GetEventNumberInFile();
347         
348     (*fDebugStreamer)<<"ErrParam"<<
349       "iter="<<fIteration<<
350       "eventNr="<<eventNr<<
351       "Cl.="<<cluster<<
352       "nclSeed="<<nclSeed<<
353       "T.="<<&param<<
354       "dy="<<dy<<
355       "dz="<<dz<<
356       "yt="<<yt<<
357       "zt="<<zt<<
358       "gcl.="<<&gcl<<
359       "gtr.="<<&gtr<<
360       "erry2="<<sy2<<
361       "errz2="<<sz2<<
362       "rmsy2="<<rmsy2<<
363       "rmsz2="<<rmsz2<< 
364       "rmsy2p30="<<rmsy2p30<<
365       "rmsz2p30="<<rmsz2p30<<   
366       "rmsy2p30R="<<rmsy2p30R<<
367       "rmsz2p30R="<<rmsz2p30R<< 
368       // normalize distance - 
369       "rdisty="<<rdistancey2<<
370       "rdistz="<<rdistancez2<<
371       "rdist="<<rdistance2<< //       
372       "\n";
373     }
374   }
375   //return 0;  // temporary
376   if (rdistance2>32) return 3;
377   
378   
379   if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)  
380     return 2;  //suspisiouce - will be changed
381   
382   if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)  
383     // strict cut on overlaped cluster
384     return  2;  //suspisiouce - will be changed
385   
386   if ( (rdistancey2>1. || rdistancez2>6.25 ) 
387        && cluster->GetType()<0){
388     seed->SetNFoundable(seed->GetNFoundable()-1);
389     return 2;    
390   }
391
392   if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
393     if (fIteration==2){
394       if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
395         if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
396           return 2;
397       }
398     }
399   }
400
401   return 0;
402 }
403
404
405
406
407
408 //_____________________________________________________________________________
409 AliTPCtracker::AliTPCtracker(const AliTPCParam *par): 
410 AliTracker(), 
411                  fkNIS(par->GetNInnerSector()/2),
412                  fInnerSec(0),
413                  fkNOS(par->GetNOuterSector()/2),
414                  fOuterSec(0),
415                  fN(0),
416                  fSectors(0),
417                  fInput(0),
418                  fOutput(0),
419                  fSeedTree(0),
420                  fTreeDebug(0),
421                  fEvent(0),
422                  fEventHLT(0),
423                  fDebug(0),
424                  fNewIO(0),
425                  fNtracks(0),
426                  fSeeds(0),
427                  fIteration(0),
428                  fkParam(0),
429          fDebugStreamer(0),
430          fUseHLTClusters(4),
431          fCrossTalkSignalArray(0),
432          fSeedsPool(0),
433                  fFreeSeedsID(500),
434                  fNFreeSeeds(0),
435                  fLastSeedID(-1)
436 {
437   //---------------------------------------------------------------------
438   // The main TPC tracker constructor
439   //---------------------------------------------------------------------
440   fInnerSec=new AliTPCtrackerSector[fkNIS];         
441   fOuterSec=new AliTPCtrackerSector[fkNOS];
442  
443   Int_t i;
444   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
445   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
446
447   fkParam = par;  
448   Int_t nrowlow = par->GetNRowLow();
449   Int_t nrowup = par->GetNRowUp();
450
451   
452   for (i=0;i<nrowlow;i++){
453     fXRow[i]     = par->GetPadRowRadiiLow(i);
454     fPadLength[i]= par->GetPadPitchLength(0,i);
455     fYMax[i]     = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
456   }
457
458   
459   for (i=0;i<nrowup;i++){
460     fXRow[i+nrowlow]      = par->GetPadRowRadiiUp(i);
461     fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
462     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
463   }
464
465   if (AliTPCReconstructor::StreamLevel()>0) {
466     fDebugStreamer = new TTreeSRedirector("TPCdebug.root","recreate");
467   }
468   //
469   fSeedsPool = new TClonesArray("AliTPCseed",1000);
470
471   // crosstalk array and matrix initialization
472   Int_t nROCs   = 72;
473   Int_t nTimeBinsAll  = par->GetMaxTBin();
474   Int_t nWireSegments = 11;
475   fCrossTalkSignalArray = new TObjArray(nROCs*2);  //  
476   fCrossTalkSignalArray->SetOwner(kTRUE);
477   for (Int_t isector=0; isector<2*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   Int_t nROCs   = 72;
1351   static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1352   //
1353   //  TTree * tree = fClustersArray.GetTree();
1354   AliInfo("LoadClusters()\n");
1355   // reset crosstalk matrix
1356   //
1357   for (Int_t isector=0; isector<nROCs*2; isector++){  //set all ellemts of crosstalk matrix to 0
1358     TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1359     if (crossTalkMatrix)(*crossTalkMatrix)*=0;
1360   }
1361
1362   TTree * tree = fInput;
1363   TBranch * br = tree->GetBranch("Segment");
1364   br->SetAddress(&clrow);
1365
1366   // Conversion of pad, row coordinates in local tracking coords.
1367   // Could be skipped here; is already done in clusterfinder
1368
1369
1370   Int_t j=Int_t(tree->GetEntries());
1371   for (Int_t i=0; i<j; i++) {
1372     br->GetEntry(i);
1373     //  
1374     Int_t sec,row;
1375     fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1376     
1377     // wire segmentID and nPadsPerSegment to be used for Xtalk calculation
1378     Int_t wireSegmentID     = fkParam->GetWireSegment(sec,row);
1379     Float_t nPadsPerSegment = (Float_t)(fkParam->GetNPadsPerSegment(wireSegmentID));
1380     TMatrixD &crossTalkSignal =  *((TMatrixD*)fCrossTalkSignalArray->At(sec));
1381     TMatrixD &crossTalkSignalBelow =  *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs));
1382     Int_t nCols=crossTalkSignal.GetNcols();
1383     Double_t missingChargeFactor= AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrectionMissingCharge();
1384     for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1385       AliTPCclusterMI *clXtalk= static_cast<AliTPCclusterMI*>(clrow->GetArray()->At(icl));
1386       Transform((AliTPCclusterMI*)(clXtalk));
1387       Int_t timeBinXtalk = clXtalk->GetTimeBin();      
1388       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 
1389       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
1390       Double_t rmsTime2   = clXtalk->GetSigmaZ2()/(fkParam->GetZWidth()*fkParam->GetZWidth()); 
1391       Double_t rmsPad2    = clXtalk->GetSigmaY2()/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec)); 
1392       if (rmsPadMin2>rmsPad2){
1393         rmsPad2=rmsPadMin2;
1394       }
1395       if (rmsTimeMin2>rmsTime2){
1396         rmsTime2=rmsTimeMin2;
1397       }
1398
1399       Double_t norm= 2.*TMath::Exp(-1.0/(2.*rmsTime2))+2.*TMath::Exp(-4.0/(2.*rmsTime2))+1.;
1400       Double_t qTotXtalk = 0.;   
1401       Double_t qTotXtalkMissing = 0.;   
1402       for (Int_t itb=timeBinXtalk-2, idelta=-2; itb<=timeBinXtalk+2; itb++,idelta++) {        
1403         if (itb<0 || itb>=nCols) continue;
1404         Double_t missingCharge=0;
1405         Double_t trf= TMath::Exp(-idelta*idelta/(2.*rmsTime2));
1406         if (missingChargeFactor>0) {
1407           for (Int_t dpad=-2; dpad<=2; dpad++){
1408             Double_t qPad =   clXtalk->GetMax()*TMath::Exp(-dpad*dpad/(2.*rmsPad2))*trf;
1409             qPad-=2.*crossTalkSignal[wireSegmentID][itb];  // missing part due crosttalk - feedback part - factor 2 used
1410             if (TMath::Nint(qPad)<=fkParam->GetZeroSup()){
1411               missingCharge+=qPad+2.*crossTalkSignal[wireSegmentID][itb];
1412             }else{
1413               missingCharge+=2.*crossTalkSignal[wireSegmentID][itb];
1414             }
1415           }
1416         }
1417         qTotXtalk = clXtalk->GetQ()*trf/norm+missingCharge*missingChargeFactor;
1418         qTotXtalkMissing = missingCharge;
1419         crossTalkSignal[wireSegmentID][itb]+= qTotXtalk/nPadsPerSegment; 
1420         crossTalkSignalBelow[wireSegmentID][itb]+= qTotXtalkMissing/nPadsPerSegment; 
1421       }
1422     }
1423     
1424     //
1425     AliTPCtrackerRow * tpcrow=0;
1426     Int_t left=0;
1427     if (sec<fkNIS*2){
1428       tpcrow = &(fInnerSec[sec%fkNIS][row]);    
1429       left = sec/fkNIS;
1430     }
1431     else{
1432       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1433       left = (sec-fkNIS*2)/fkNOS;
1434     }
1435     if (left ==0){
1436       tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1437       for (Int_t k=0;k<tpcrow->GetN1();++k) 
1438         tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1439     }
1440     if (left ==1){
1441       tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1442       for (Int_t k=0;k<tpcrow->GetN2();++k) 
1443         tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1444     }
1445   }
1446   if (AliTPCReconstructor::StreamLevel()&kStreamCrosstalkMatrix) {
1447     //
1448     // dump the crosstalk matrices to tree for further investigation
1449     //     a.) to estimate fluctuation of pedestal in indiviula wire segments
1450     //     b.) to check correlation between regions
1451     //     c.) to check relative conribution of signal below threshold to crosstalk
1452     for (Int_t isector=0; isector<nROCs; isector++){  //set all ellemts of crosstalk matrix to 0
1453       TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1454       TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
1455       TVectorD vecAll(crossTalkMatrix->GetNrows());
1456       TVectorD vecBelow(crossTalkMatrix->GetNrows());
1457       //
1458       for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
1459         for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
1460           vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
1461           vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
1462         }
1463         (*fDebugStreamer)<<"crosstalkMatrix"<<
1464           "sector="<<isector<<
1465           "itime="<<itime<<
1466           "vecAll.="<<&vecAll<<                // crosstalk charge + charge below threshold
1467           "vecBelow.="<<&vecBelow<<            // crosstalk contribution from signal below threshold
1468           "\n";
1469       }
1470     }
1471   }
1472
1473
1474   //
1475   clrow->Clear("C");
1476   LoadOuterSectors();
1477   LoadInnerSectors();
1478
1479   cout << " =================================================================================================================================== " << endl;
1480   cout << " AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() =  " << AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() << endl;
1481   cout << " AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()  =  " << AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()  << endl;
1482   cout << " =================================================================================================================================== " << endl;
1483
1484   if (AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection()) ApplyTailCancellation();
1485   if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) ApplyXtalkCorrection();
1486   //if (AliTPCReconstructor::GetRecoParam()->GetUseOulierClusterFilter()) FilterOutlierClusters();  
1487   return 0;
1488 }
1489
1490 void    AliTPCtracker::FilterOutlierClusters(){
1491   //
1492   // filter outlier clusters  
1493   //
1494   /*
1495     1.)..... booking part
1496     nSectors=72;
1497     nTimeBins=fParam->Get....
1498     TH2F hisTime("","", sector,0,sector, nTimeBins,0,nTimeBins);
1499     TH2F hisPadRow("","", sector,0,sector, nPadRows,0,nPadRows);
1500     2.) .... filling part
1501     .... cluster loop { hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin()); }
1502     
1503     3.) ...filtering part 
1504     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... }
1505     
1506     sector loop
1507     { 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 } 
1508     //
1509     4. Disabling clusters
1510     
1511   */
1512   
1513   //
1514   // 1.) booking part 
1515   // 
1516   //  AliTPCcalibDB *db=AliTPCcalibDB::Instance();
1517   Int_t nSectors=AliTPCROC::Instance()->GetNSectors(); 
1518   Int_t nTimeBins= 1100; // *Bug here - we should get NTimeBins from ALTRO - Parameters not relyable
1519   Int_t nPadRows=(AliTPCROC::Instance()->GetNRows(0) + AliTPCROC::Instance()->GetNRows(36));
1520   // parameters for filtering
1521   const Double_t nSigmaCut=9.;           // should be in recoParam ?
1522   const Double_t offsetTime=100;         // should be in RecoParam ?  -
1523   const Double_t offsetPadRow=300;       // should be in RecoParam ?
1524   const Double_t offsetTimeAccept=8;     // should be in RecoParam ?  - obtained as mean +1 rms in high IR pp
1525   TH2F hisTime("hisSectorTime","hisSectorTime", nSectors,0,nSectors, nTimeBins,0,nTimeBins);
1526   TH2F hisPadRow("hisSectorRow","hisSectorRow", nSectors,0,nSectors, nPadRows,0,nPadRows);
1527   //
1528   // 2.) Filling part -- loop over clusters
1529   //
1530   for (Int_t isector=0; isector<36; isector++){      // loop over sectors
1531     for (Int_t iside=0; iside<2; iside++){           // loop over sides A/C
1532       AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1533       Int_t nrows = sector.GetNRows();
1534       for (Int_t row = 0;row<nrows;row++){           // loop over rows       
1535         AliTPCtrackerRow&  tpcrow = sector[row];       
1536         Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
1537         if (iside>0) ncl=tpcrow.GetN2();
1538         for (Int_t i=0;i<ncl;i++) {  // loop over clusters
1539           AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1540           hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin()); 
1541           hisPadRow.Fill(cluster->GetDetector(),cluster->GetRow()); 
1542         }
1543       } 
1544     } 
1545   } 
1546
1547   //
1548   // 3. Filtering part
1549   //
1550   TVectorD vecTime(nTimeBins);
1551   TVectorD vecPadRow(nPadRows);
1552   TVectorD vecMedianSectorTime(nSectors);
1553   TVectorD vecRMSSectorTime(nSectors);
1554   TVectorD vecMedianSectorTimeOut6(nSectors);
1555   TVectorD vecMedianSectorTimeOut9(nSectors);//
1556   TVectorD vecMedianSectorTimeOut(nSectors);//
1557   TVectorD vecMedianSectorPadRow(nSectors);
1558   TVectorD vecRMSSectorPadRow(nSectors);
1559   TVectorD vecMedianSectorPadRowOut6(nSectors);
1560   TVectorD vecMedianSectorPadRowOut9(nSectors);
1561   TVectorD vecMedianSectorPadRowOut(nSectors);
1562   TVectorD vecSectorOut6(nSectors);
1563   TVectorD vecSectorOut9(nSectors);
1564   TMatrixD matSectorCluster(nSectors,2);
1565   //
1566   // 3.a)  median, rms calculations for hisTime 
1567   //
1568   for (Int_t isec=0; isec<nSectors; isec++){
1569     vecMedianSectorTimeOut6[isec]=0;
1570     vecMedianSectorTimeOut9[isec]=0;
1571     for (Int_t itime=0; itime<nTimeBins; itime++){
1572       vecTime[itime]=hisTime.GetBinContent(isec+1, itime+1);      
1573     }
1574     Double_t median= TMath::Mean(nTimeBins,vecTime.GetMatrixArray());
1575     Double_t rms= TMath::RMS(nTimeBins,vecTime.GetMatrixArray()); 
1576     vecMedianSectorTime[isec]=median;
1577     vecRMSSectorTime[isec]=rms;
1578     if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector TimeStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
1579     //
1580     // declare outliers
1581     for (Int_t itime=0; itime<nTimeBins; itime++){
1582       Double_t entries= hisTime.GetBinContent(isec+1, itime+1); 
1583       if (entries>median+6.*rms+offsetTime) {
1584         vecMedianSectorTimeOut6[isec]+=1;
1585       }
1586       if (entries>median+9.*rms+offsetTime) {
1587         vecMedianSectorTimeOut9[isec]+=1;
1588       }
1589     }
1590   }    
1591   //
1592   // 3.b) median, rms calculations for hisPadRow
1593   // 
1594   for (Int_t isec=0; isec<nSectors; isec++){
1595     vecMedianSectorPadRowOut6[isec]=0;
1596     vecMedianSectorPadRowOut9[isec]=0;
1597     for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1598       vecPadRow[ipadrow]=hisPadRow.GetBinContent(isec+1, ipadrow+1);      
1599     }
1600     Int_t nPadRowsSector= AliTPCROC::Instance()->GetNRows(isec);
1601     Double_t median= TMath::Mean(nPadRowsSector,vecPadRow.GetMatrixArray());
1602     Double_t rms= TMath::RMS(nPadRowsSector,vecPadRow.GetMatrixArray());
1603     vecMedianSectorPadRow[isec]=median;
1604     vecRMSSectorPadRow[isec]=rms;
1605     if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector PadRowStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
1606     //
1607     // declare outliers
1608     for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1609       Double_t entries= hisPadRow.GetBinContent(isec+1, ipadrow+1);
1610       if (entries>median+6.*rms+offsetPadRow) {
1611         vecMedianSectorPadRowOut6[isec]+=1;
1612       }
1613       if (entries>median+9.*rms+offsetPadRow) {
1614         vecMedianSectorPadRowOut9[isec]+=1;
1615       }
1616     }
1617   }
1618   //
1619   // 3.c) filter outlier sectors
1620   //
1621   Double_t medianSectorTime = TMath::Median(nSectors, vecTime.GetMatrixArray());
1622   Double_t mean69SectorTime, rms69SectorTime=0;
1623   AliMathBase::EvaluateUni(nSectors,  vecTime.GetMatrixArray(), mean69SectorTime,rms69SectorTime,69);
1624   for (Int_t isec=0; isec<nSectors; isec++){
1625     vecSectorOut6[isec]=0;
1626     vecSectorOut9[isec]=0;
1627     matSectorCluster(isec,0)=0;
1628     matSectorCluster(isec,1)=0;
1629     if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+6.*(rms69SectorTime+ offsetTimeAccept))) {
1630       vecSectorOut6[isec]=1;
1631     }
1632     if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+9.*(rms69SectorTime+ offsetTimeAccept))){
1633       vecSectorOut9[isec]=1;
1634     }
1635   }
1636   // light version of export variable
1637   Int_t filteredSector= vecSectorOut9.Sum();                  // light version of export variable
1638   Int_t filteredSectorTime= vecMedianSectorTimeOut9.Sum();
1639   Int_t filteredSectorPadRow= vecMedianSectorPadRowOut9.Sum();
1640   if (fEvent) if (fEvent->GetHeader()){
1641     fEvent->GetHeader()->SetTPCNoiseFilterCounter(0,TMath::Min(filteredSector,255));
1642     fEvent->GetHeader()->SetTPCNoiseFilterCounter(1,TMath::Min(filteredSectorTime,255));
1643     fEvent->GetHeader()->SetTPCNoiseFilterCounter(2,TMath::Min(filteredSectorPadRow,255));
1644   }
1645  
1646   //
1647   // 4. Disabling clusters in outlier layers
1648   //
1649   Int_t counterAll=0;
1650   Int_t counterOut=0;
1651   for (Int_t isector=0; isector<36; isector++){      // loop over sectors
1652     for (Int_t iside=0; iside<2; iside++){           // loop over sides A/C
1653       AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1654       Int_t nrows = sector.GetNRows();
1655       for (Int_t row = 0;row<nrows;row++){           // loop over rows       
1656         AliTPCtrackerRow&  tpcrow = sector[row];       
1657         Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
1658         if (iside>0) ncl=tpcrow.GetN2();
1659         for (Int_t i=0;i<ncl;i++) {  // loop over clusters
1660           AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1661           Double_t medianTime=vecMedianSectorTime[cluster->GetDetector()];
1662           Double_t medianPadRow=vecMedianSectorPadRow[cluster->GetDetector()];
1663           Double_t rmsTime=vecRMSSectorTime[cluster->GetDetector()];
1664           Double_t rmsPadRow=vecRMSSectorPadRow[cluster->GetDetector()];
1665           Int_t entriesPadRow=hisPadRow.GetBinContent(cluster->GetDetector()+1, cluster->GetRow()+1);
1666           Int_t entriesTime=hisTime.GetBinContent(cluster->GetDetector()+1, cluster->GetTimeBin()+1);
1667           Bool_t isOut=kFALSE;
1668           if (vecSectorOut9[cluster->GetDetector()]>0.5) {
1669             isOut=kTRUE;
1670           }
1671           
1672           if (entriesTime>medianTime+nSigmaCut*rmsTime+offsetTime) {
1673             isOut=kTRUE;
1674             vecMedianSectorTimeOut[cluster->GetDetector()]++;
1675           }
1676           if (entriesPadRow>medianPadRow+nSigmaCut*rmsPadRow+offsetPadRow) {
1677             isOut=kTRUE;
1678             vecMedianSectorPadRowOut[cluster->GetDetector()]++;
1679           }
1680           counterAll++;
1681           matSectorCluster(cluster->GetDetector(),0)+=1;
1682           if (isOut){
1683             cluster->Disable();
1684             counterOut++;
1685             matSectorCluster(cluster->GetDetector(),1)+=1;
1686           }
1687         }
1688       }
1689     }
1690   }
1691   for (Int_t isec=0; isec<nSectors; isec++){
1692     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());
1693   }
1694   //
1695   // dump info to streamer - for later tuning of cuts
1696   //
1697   if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) {  // stream TPC data ouliers filtering infomation
1698     AliLog::Flush();
1699     AliInfo(TString::Format("Cluster counter: (%d/%d) (Filtered/All)",counterOut,counterAll).Data());
1700     for (Int_t iSec=0; iSec<nSectors; iSec++){
1701       if (vecSectorOut9[iSec]>0 ||  matSectorCluster(iSec,1)>0) {
1702         AliInfo(TString::Format("Filtered sector\t%d",iSec).Data());
1703         Double_t vecMedTime =TMath::Median(72,vecMedianSectorTime.GetMatrixArray());
1704         Double_t vecMedPadRow =TMath::Median(72,vecMedianSectorPadRow.GetMatrixArray());
1705         Double_t vecMedCluster=(counterAll-counterOut)/72;
1706         AliInfo(TString::Format("VecMedianSectorTime\t(%4.4f/%4.4f/%4.4f)",       vecMedianSectorTimeOut[iSec],vecMedianSectorTime[iSec],vecMedTime).Data());
1707         AliInfo(TString::Format("VecMedianSectorPadRow\t(%4.4f/%4.4f/%4.4f)",     vecMedianSectorPadRowOut[iSec],vecMedianSectorPadRow[iSec],vecMedPadRow).Data());
1708         AliInfo(TString::Format("MatSectorCluster\t(%4.4f/%4.4f/%4.4f)\n",          matSectorCluster(iSec,1), matSectorCluster(iSec,0),  vecMedCluster).Data());
1709         AliLog::Flush();
1710       }
1711     }
1712     AliLog::Flush();
1713     Int_t eventNr = fEvent->GetEventNumberInFile();
1714     (*fDebugStreamer)<<"filterClusterInfo"<<
1715       // minimal set variables for the ESDevent
1716       "eventNr="<<eventNr<<
1717       "counterAll="<<counterAll<<
1718       "counterOut="<<counterOut<<
1719       //
1720       "filteredSector="<<filteredSector<<                        //  counter filtered sectors                   
1721       "filteredSectorTime="<<filteredSectorTime<<                //  counter filtered time bins
1722       "filteredSectorPadRow="<<filteredSectorPadRow<<            //  counter filtered pad-rows
1723       // per sector outlier information
1724       "medianSectorTime="<<medianSectorTime<<                    // median number of clusters per sector/timebin
1725       "mean69SectorTime="<<mean69SectorTime<<                    // LTM statistic  mean of clusters per sector/timebin
1726       "rms69SectorTime="<<rms69SectorTime<<                      // LTM statistic  RMS of clusters per sector/timebin
1727       "vecSectorOut6.="<<&vecSectorOut6<<                        // flag array sector  - 6 sigma +accept margin outlier
1728       "vecSectorOut9.="<<&vecSectorOut9<<                        // flag array sector  - 9 sigma + accept margin outlier
1729       // per sector/timebin outlier detection
1730       "vecMedianSectorTime.="<<&vecMedianSectorTime<<
1731       "vecRMSSectorTime.="<<&vecRMSSectorTime<<
1732       "vecMedianSectorTimeOut6.="<<&vecMedianSectorTimeOut6<<
1733       "vecMedianSectorTimeOut9.="<<&vecMedianSectorTimeOut9<<
1734       "vecMedianSectorTimeOut0.="<<&vecMedianSectorTimeOut<<
1735       // per sector/pad-row outlier detection
1736       "vecMedianSectorPadRow.="<<&vecMedianSectorPadRow<<
1737       "vecRMSSectorPadRow.="<<&vecRMSSectorPadRow<<
1738       "vecMedianSectorPadRowOut6.="<<&vecMedianSectorPadRowOut6<<
1739       "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut9<<
1740       "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut<<
1741       "\n";
1742     ((*fDebugStreamer)<<"filterClusterInfo").GetTree()->Write();
1743     fDebugStreamer->GetFile()->Flush();
1744   }
1745 }
1746
1747 void AliTPCtracker::UnloadClusters()
1748 {
1749   //
1750   // unload clusters from the memory
1751   //
1752   Int_t nrows = fOuterSec->GetNRows();
1753   for (Int_t sec = 0;sec<fkNOS;sec++)
1754     for (Int_t row = 0;row<nrows;row++){
1755       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1756       //      if (tpcrow){
1757       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1758       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1759       //}
1760       tpcrow->ResetClusters();
1761     }
1762   //
1763   nrows = fInnerSec->GetNRows();
1764   for (Int_t sec = 0;sec<fkNIS;sec++)
1765     for (Int_t row = 0;row<nrows;row++){
1766       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1767       //if (tpcrow){
1768       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1769       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1770       //}
1771       tpcrow->ResetClusters();
1772     }
1773
1774   return ;
1775 }
1776
1777 void AliTPCtracker::FillClusterArray(TObjArray* array) const{
1778   //
1779   // Filling cluster to the array - For visualization purposes
1780   //
1781   Int_t nrows=0;
1782   nrows = fOuterSec->GetNRows();
1783   for (Int_t sec = 0;sec<fkNOS;sec++)
1784     for (Int_t row = 0;row<nrows;row++){
1785       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1786       if (!tpcrow) continue;
1787       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1788         array->AddLast((TObject*)((*tpcrow)[icl]));
1789       }
1790     } 
1791   nrows = fInnerSec->GetNRows();
1792   for (Int_t sec = 0;sec<fkNIS;sec++)
1793     for (Int_t row = 0;row<nrows;row++){
1794       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1795       if (!tpcrow) continue;
1796       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1797         array->AddLast((TObject*)(*tpcrow)[icl]);
1798       }
1799     }
1800 }
1801
1802
1803 void   AliTPCtracker::Transform(AliTPCclusterMI * cluster){
1804   //
1805   // transformation
1806   //
1807   AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1808   AliTPCTransform *transform = calibDB->GetTransform() ;
1809   if (!transform) {
1810     AliFatal("Tranformations not in calibDB");
1811     return;
1812   }
1813   transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1814   Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),static_cast<Double_t>(cluster->GetPad()),static_cast<Double_t>(cluster->GetTimeBin())};
1815   Int_t i[1]={cluster->GetDetector()};
1816   transform->Transform(x,i,0,1);  
1817   //  if (cluster->GetDetector()%36>17){
1818   //  x[1]*=-1;
1819   //}
1820
1821   //
1822   // in debug mode  check the transformation
1823   //
1824   if ((AliTPCReconstructor::StreamLevel()&kStreamTransform)>0) { 
1825     Float_t gx[3];
1826     cluster->GetGlobalXYZ(gx);
1827     Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1828     TTreeSRedirector &cstream = *fDebugStreamer;
1829     cstream<<"Transform"<<  // needed for debugging of the cluster transformation, resp. used for later visualization 
1830       "event="<<event<<
1831       "x0="<<x[0]<<
1832       "x1="<<x[1]<<
1833       "x2="<<x[2]<<
1834       "gx0="<<gx[0]<<
1835       "gx1="<<gx[1]<<
1836       "gx2="<<gx[2]<<
1837       "Cl.="<<cluster<<
1838       "\n"; 
1839   }
1840   cluster->SetX(x[0]);
1841   cluster->SetY(x[1]);
1842   cluster->SetZ(x[2]);
1843   // The old stuff:
1844   //
1845   // 
1846   //
1847   //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1848   if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1849     TGeoHMatrix  *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1850     //TGeoHMatrix  mat;
1851     Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1852     Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1853     if (mat) mat->LocalToMaster(pos,posC);
1854     else{
1855       // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1856     }
1857     cluster->SetX(posC[0]);
1858     cluster->SetY(posC[1]);
1859     cluster->SetZ(posC[2]);
1860   }
1861 }
1862
1863 void  AliTPCtracker::ApplyXtalkCorrection(){
1864   //
1865   // ApplyXtalk correction 
1866   // Loop over all clusters
1867   //      add to each cluster signal corresponding to common Xtalk mode for given time bin at given wire segment
1868   // cluster loop
1869   for (Int_t isector=0; isector<36; isector++){  //loop tracking sectors
1870     for (Int_t iside=0; iside<2; iside++){       // loop over sides A/C
1871       AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1872       Int_t nrows     = sector.GetNRows();       
1873       for (Int_t row = 0;row<nrows;row++){           // loop over rows       
1874         AliTPCtrackerRow&  tpcrow = sector[row];     // row object   
1875         Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
1876         if (iside>0) ncl=tpcrow.GetN2();
1877         Int_t xSector=0;    // sector number in the TPC convention 0-72
1878         if (isector<18){  //if IROC
1879           xSector=isector+(iside>0)*18;
1880         }else{
1881           xSector=isector+18;  // isector -18 +36   
1882           if (iside>0) xSector+=18;
1883         }       
1884         TMatrixD &crossTalkMatrix= *((TMatrixD*)fCrossTalkSignalArray->At(xSector));
1885         Int_t wireSegmentID     = fkParam->GetWireSegment(xSector,row);
1886         for (Int_t i=0;i<ncl;i++) {
1887           AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1888           Int_t iTimeBin=TMath::Nint(cluster->GetTimeBin());
1889           Double_t xTalk= crossTalkMatrix[wireSegmentID][iTimeBin];
1890           cluster->SetMax(cluster->GetMax()+xTalk);
1891           const Double_t kDummy=4;
1892           Double_t sumxTalk=xTalk*kDummy; // should be calculated via time response function
1893           cluster->SetQ(cluster->GetQ()+sumxTalk);
1894
1895
1896           if ((AliTPCReconstructor::StreamLevel()&kStreamXtalk)>0) {  // flag: stream crosstalk correctio as applied to cluster
1897             TTreeSRedirector &cstream = *fDebugStreamer;
1898             if (gRandom->Rndm() > 0.){
1899               cstream<<"Xtalk"<<
1900                 "isector=" << isector <<               // sector [0,36]
1901                 "iside=" << iside <<                   // side A or C
1902                 "row=" << row <<                       // padrow
1903                 "i=" << i <<                           // index of the cluster 
1904                 "xSector=" << xSector <<               // sector [0,72] 
1905                 "wireSegmentID=" << wireSegmentID <<   // anode wire segment id [0,10] 
1906                 "iTimeBin=" << iTimeBin <<             // timebin of the corrected cluster 
1907                 "xTalk=" << xTalk <<                   // Xtalk contribution added to Qmax
1908                 "sumxTalk=" << sumxTalk <<             // Xtalk contribution added to Qtot (roughly 3*Xtalk) 
1909                 "cluster.=" << cluster <<              // corrected cluster object 
1910                 "\n";
1911             }
1912           }// dump the results to the debug streamer if in debug mode
1913
1914
1915
1916
1917
1918         }
1919       }
1920     }
1921   }
1922 }
1923
1924 void  AliTPCtracker::ApplyTailCancellation(){
1925   //
1926   // Correct the cluster charge for the ion tail effect 
1927   // The TimeResponse function accessed via  AliTPCcalibDB (TPC/Calib/IonTail)
1928   //
1929
1930   // Retrieve
1931   TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
1932   if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
1933   TObject *rocFactorIROC  = ionTailArr->FindObject("factorIROC");
1934   TObject *rocFactorOROC  = ionTailArr->FindObject("factorOROC");   
1935   Float_t factorIROC      = (atof(rocFactorIROC->GetTitle()));
1936   Float_t factorOROC      = (atof(rocFactorOROC->GetTitle()));
1937
1938   // find the number of clusters for the whole TPC (nclALL)
1939   Int_t nclALL=0;
1940   for (Int_t isector=0; isector<36; isector++){
1941     AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1942     nclALL += sector.GetNClInSector(0);
1943     nclALL += sector.GetNClInSector(1);
1944   }
1945
1946   // start looping over all clusters 
1947   for (Int_t iside=0; iside<2; iside++){    // loop over sides
1948     //
1949     //
1950     for (Int_t secType=0; secType<2; secType++){  //loop over inner or outer sector
1951       // cache experimantal tuning factor for the different chamber type 
1952       const Float_t ampfactor = (secType==0)?factorIROC:factorOROC;
1953       std::cout << " ampfactor = " << ampfactor << std::endl;
1954       //
1955       for (Int_t sec = 0;sec<fkNOS;sec++){        //loop overs sectors
1956         //
1957         //
1958         // Cache time response functions and their positons to COG of the cluster        
1959         TGraphErrors ** graphRes   = new TGraphErrors *[20];
1960         Float_t * indexAmpGraphs   = new Float_t[20];      
1961         for (Int_t icache=0; icache<20; icache++) 
1962         {
1963           graphRes[icache]       = NULL;
1964           indexAmpGraphs[icache] = 0;
1965         }
1966         /////////////////////////////  --> position fo sie loop
1967         if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(sec+36*secType+18*iside,graphRes,indexAmpGraphs))
1968         {
1969           continue;
1970         }
1971         
1972         AliTPCtrackerSector &sector= (secType==0)?fInnerSec[sec]:fOuterSec[sec];  
1973         Int_t nrows     = sector.GetNRows();                                       // number of rows
1974         Int_t nclSector = sector.GetNClInSector(iside);                            // ncl per sector to be used for debugging
1975
1976         for (Int_t row = 0;row<nrows;row++){           // loop over rows
1977
1978           AliTPCtrackerRow&  tpcrow = sector[row];     // row object   
1979           Int_t ncl = tpcrow.GetN1();                  // number of clusters in the row
1980           if (iside>0) ncl=tpcrow.GetN2();
1981         
1982           // Order clusters in time for the proper correction of ion tail
1983           Float_t qTotArray[ncl];                      // arrays to be filled with modified Qtot and Qmax values in order to avoid float->int conversion  
1984           Float_t qMaxArray[ncl];
1985           Int_t sortedClusterIndex[ncl];
1986           Float_t sortedClusterTimeBin[ncl];
1987           TObjArray *rowClusterArray = new TObjArray(ncl);  // cache clusters for each row  
1988           for (Int_t i=0;i<ncl;i++) 
1989           {
1990             qTotArray[i]=0;
1991             qMaxArray[i]=0;
1992             sortedClusterIndex[i]=i;
1993             AliTPCclusterMI *rowcl= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1994             if (rowcl) {
1995               rowClusterArray->AddAt(rowcl,i);
1996             } else {
1997               rowClusterArray->RemoveAt(i);
1998             }
1999             // Fill the timebin info to the array in order to sort wrt tb
2000             if (!rowcl) {
2001               sortedClusterTimeBin[i]=0.0;
2002             } else {
2003             sortedClusterTimeBin[i] = rowcl->GetTimeBin();
2004             }
2005
2006           } 
2007           TMath::Sort(ncl,sortedClusterTimeBin,sortedClusterIndex,kFALSE);       // sort clusters in time
2008      
2009           // Main cluster correction loops over clusters
2010           for (Int_t icl0=0; icl0<ncl;icl0++){    // first loop over clusters
2011
2012             AliTPCclusterMI *cl0= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl0]));
2013             
2014             if (!cl0) continue;
2015             Int_t nclPad=0;                       
2016             for (Int_t icl1=0; icl1<ncl;icl1++){  // second loop over clusters     
2017               AliTPCclusterMI *cl1= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl1]));
2018               if (!cl1) continue;
2019               if (TMath::Abs(cl0->GetPad()-cl1->GetPad())>4) continue;           // no contribution if far away in pad direction
2020               if (cl0->GetTimeBin()<= cl1->GetTimeBin()) continue;               // no contibution to the tail if later
2021               if (TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin())>600) continue; // out of the range of response function
2022
2023               if (TMath::Abs(cl0->GetPad()-cl1->GetPad())<4) nclPad++;           // count ncl for every pad for debugging
2024             
2025               // Get the correction values for Qmax and Qtot and find total correction for a given cluster
2026               Double_t ionTailMax=0.;  
2027               Double_t ionTailTotal=0.;  
2028               GetTailValue(ampfactor,ionTailMax,ionTailTotal,graphRes,indexAmpGraphs,cl0,cl1);
2029               ionTailMax=TMath::Abs(ionTailMax);
2030               ionTailTotal=TMath::Abs(ionTailTotal);
2031               qTotArray[icl0]+=ionTailTotal;
2032               qMaxArray[icl0]+=ionTailMax;
2033
2034               // Dump some info for debugging while clusters are being corrected
2035               if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {  // flag: stream ion tail correction  as applied to cluster
2036                 TTreeSRedirector &cstream = *fDebugStreamer;
2037                 if (gRandom->Rndm() > 0.999){
2038                   cstream<<"IonTail"<<
2039                       "cl0.="         <<cl0          <<   // cluster 0 (to be corrected)
2040                       "cl1.="         <<cl1          <<   // cluster 1 (previous cluster)
2041                       "ionTailTotal=" <<ionTailTotal <<   // ion Tail from cluster 1 contribution to cluster0
2042                       "ionTailMax="   <<ionTailMax   <<   // ion Tail from cluster 1 contribution to cluster0 
2043                       "\n";
2044                 }
2045               }// dump the results to the debug streamer if in debug mode
2046             
2047             }//end of second loop over clusters
2048             
2049             // Set corrected values of the corrected cluster          
2050             cl0->SetQ(TMath::Nint(Float_t(cl0->GetQ())+Float_t(qTotArray[icl0])));
2051             cl0->SetMax(TMath::Nint(Float_t(cl0->GetMax())+qMaxArray[icl0]));
2052           
2053             // Dump some info for debugging after clusters are corrected 
2054             if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {
2055               TTreeSRedirector &cstream = *fDebugStreamer;
2056               if (gRandom->Rndm() > 0.999){
2057               cstream<<"IonTailCorrected"<<
2058                   "cl0.="                     << cl0              <<   // cluster 0 with huge Qmax
2059                   "ionTailTotalPerCluster="   << qTotArray[icl0]  <<
2060                   "ionTailMaxPerCluster="     << qMaxArray[icl0]  <<
2061                   "nclALL="                   << nclALL           <<
2062                   "nclSector="                << nclSector        <<
2063                   "nclRow="                   << ncl              <<
2064                   "nclPad="                   << nclPad           <<
2065                   "row="                      << row              <<
2066                   "sector="                   << sec              <<
2067                   "icl0="                     << icl0             <<
2068                   "\n";
2069               }
2070             }// dump the results to the debug streamer if in debug mode
2071           
2072           }//end of first loop over cluster
2073           delete rowClusterArray;
2074         }//end of loop over rows
2075         for (int i=0; i<20; i++) delete graphRes[i];
2076         delete [] graphRes;
2077         delete [] indexAmpGraphs;
2078       
2079       }//end of loop over sectors
2080     }//end of loop over IROC/OROC
2081   }// end of side loop
2082 }
2083 //_____________________________________________________________________________
2084 void AliTPCtracker::GetTailValue(Float_t ampfactor,Double_t &ionTailMax, Double_t &ionTailTotal,TGraphErrors **graphRes,Float_t *indexAmpGraphs,AliTPCclusterMI *cl0,AliTPCclusterMI *cl1){
2085
2086   //
2087   // Function in order to calculate the amount of the correction to be added for a given cluster, return values are ionTailTaoltal and ionTailMax
2088   // Parameters:
2089   // cl0 -  cluster to be modified
2090   // cl1 -  source cluster ion tail of this cluster will be added to the cl0 (accroding time and pad response function)
2091   // 
2092   const Double_t kMinPRF    = 0.5;                           // minimal PRF width
2093   ionTailTotal              = 0.;                            // correction value to be added to Qtot of cl0
2094   ionTailMax                = 0.;                            // correction value to be added to Qmax of cl0
2095
2096   Float_t qTot0             =  cl0->GetQ();                  // cl0 Qtot info
2097   Float_t qTot1             =  cl1->GetQ();                  // cl1 Qtot info
2098   Int_t sectorPad           =  cl1->GetDetector();           // sector number
2099   Int_t padcl0              =  TMath::Nint(cl0->GetPad());   // pad0
2100   Int_t padcl1              =  TMath::Nint(cl1->GetPad());   // pad1
2101   Float_t padWidth          = (sectorPad < 36)?0.4:0.6;      // pad width in cm
2102   const Int_t deltaTimebin  =  TMath::Nint(TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin()))+12;  //distance between pads of cl1 and cl0 increased by 12 bins
2103   Double_t rmsPad1          = (cl1->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl1->GetSigmaY2())/padWidth);
2104   Double_t rmsPad0          = (cl0->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl0->GetSigmaY2())/padWidth);
2105   
2106    
2107   
2108   Double_t sumAmp1=0.;
2109   for (Int_t idelta =-2; idelta<=2;idelta++){
2110     sumAmp1+=TMath::Exp(-idelta*idelta/(2*rmsPad1));
2111   }
2112
2113   Double_t sumAmp0=0.;
2114   for (Int_t idelta =-2; idelta<=2;idelta++){
2115     sumAmp0+=TMath::Exp(-idelta*idelta/(2*rmsPad0));
2116   }
2117
2118   // Apply the correction  -->   cl1 corrects cl0 (loop over cl1's pads and find which pads of cl0 are going to be corrected)
2119   Int_t padScan=2;      // +-2 pad-timebin window will be scanned
2120   for (Int_t ipad1=padcl1-padScan; ipad1<=padcl1+padScan; ipad1++) {
2121     //
2122     //
2123     Float_t  deltaPad1  = TMath::Abs(cl1->GetPad()-(Float_t)ipad1);
2124     Double_t amp1       = (TMath::Exp(-(deltaPad1*deltaPad1)/(2*rmsPad1)))/sumAmp1;  // normalized pad response function
2125     Float_t qTotPad1    = amp1*qTot1;                                               // used as a factor to multipliy the response function
2126       
2127     // find closest value of cl1 to COG (among the time response functions' amplitude array --> to select proper t.r.f.)
2128     Int_t ampIndex = 0;
2129     Float_t diffAmp  = TMath::Abs(deltaPad1-indexAmpGraphs[0]);
2130     for (Int_t j=0;j<20;j++) {
2131       if (diffAmp > TMath::Abs(deltaPad1-indexAmpGraphs[j]) && indexAmpGraphs[j]!=0)
2132         {
2133           diffAmp  = TMath::Abs(deltaPad1-indexAmpGraphs[j]);
2134           ampIndex = j;
2135         }
2136     }
2137     if (!graphRes[ampIndex]) continue;
2138     if (deltaTimebin+2 >= graphRes[ampIndex]->GetN()) continue;
2139     if (graphRes[ampIndex]->GetY()[deltaTimebin+2]>=0) continue;
2140      
2141     for (Int_t ipad0=padcl0-padScan; ipad0<=padcl0+padScan; ipad0++) {
2142       //
2143       //
2144       if (ipad1!=ipad0) continue;                                     // check if ipad1 channel sees ipad0 channel, if not no correction to be applied.
2145       
2146       Float_t deltaPad0  = TMath::Abs(cl0->GetPad()-(Float_t)ipad0);
2147       Double_t amp0      = (TMath::Exp(-(deltaPad0*deltaPad0)/(2*rmsPad0)))/sumAmp0;  // normalized pad resp function
2148       Float_t qMaxPad0   = amp0*qTot0;
2149            
2150       // Add 5 timebin range contribution around the max peak (-+2 tb window)
2151       for (Int_t itb=deltaTimebin-2; itb<=deltaTimebin+2; itb++) {
2152
2153         if (itb<0) continue; 
2154         if (itb>=graphRes[ampIndex]->GetN()) continue;
2155        
2156         // calculate contribution to qTot
2157         Float_t tailCorr =  TMath::Abs((qTotPad1*ampfactor)*(graphRes[ampIndex])->GetY()[itb]);
2158         if (ipad1!=padcl0) { 
2159           ionTailTotal += TMath::Min(qMaxPad0,tailCorr);   // for side pad
2160         } else {             
2161           ionTailTotal += tailCorr;                        // for center pad
2162         }
2163         // calculate contribution to qMax
2164         if (itb == deltaTimebin && ipad1 == padcl0) ionTailMax += tailCorr;   
2165         
2166       } // end of tb correction loop which is applied over 5 tb range
2167
2168     } // end of cl0 loop
2169   } // end of cl1 loop
2170   
2171 }
2172
2173 //_____________________________________________________________________________
2174 Int_t AliTPCtracker::LoadOuterSectors() {
2175   //-----------------------------------------------------------------
2176   // This function fills outer TPC sectors with clusters.
2177   //-----------------------------------------------------------------
2178   Int_t nrows = fOuterSec->GetNRows();
2179   UInt_t index=0;
2180   for (Int_t sec = 0;sec<fkNOS;sec++)
2181     for (Int_t row = 0;row<nrows;row++){
2182       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);  
2183       Int_t sec2 = sec+2*fkNIS;
2184       //left
2185       Int_t ncl = tpcrow->GetN1();
2186       while (ncl--) {
2187         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
2188         index=(((sec2<<8)+row)<<16)+ncl;
2189         tpcrow->InsertCluster(c,index);
2190       }
2191       //right
2192       ncl = tpcrow->GetN2();
2193       while (ncl--) {
2194         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
2195         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
2196         tpcrow->InsertCluster(c,index);
2197       }
2198       //
2199       // write indexes for fast acces
2200       //
2201       for (Int_t i=0;i<510;i++)
2202         tpcrow->SetFastCluster(i,-1);
2203       for (Int_t i=0;i<tpcrow->GetN();i++){
2204         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
2205         tpcrow->SetFastCluster(zi,i);  // write index
2206       }
2207       Int_t last = 0;
2208       for (Int_t i=0;i<510;i++){
2209         if (tpcrow->GetFastCluster(i)<0)
2210           tpcrow->SetFastCluster(i,last);
2211         else
2212           last = tpcrow->GetFastCluster(i);
2213       }
2214     }  
2215   fN=fkNOS;
2216   fSectors=fOuterSec;
2217   return 0;
2218 }
2219
2220
2221 //_____________________________________________________________________________
2222 Int_t  AliTPCtracker::LoadInnerSectors() {
2223   //-----------------------------------------------------------------
2224   // This function fills inner TPC sectors with clusters.
2225   //-----------------------------------------------------------------
2226   Int_t nrows = fInnerSec->GetNRows();
2227   UInt_t index=0;
2228   for (Int_t sec = 0;sec<fkNIS;sec++)
2229     for (Int_t row = 0;row<nrows;row++){
2230       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
2231       //
2232       //left
2233       Int_t ncl = tpcrow->GetN1();
2234       while (ncl--) {
2235         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
2236         index=(((sec<<8)+row)<<16)+ncl;
2237         tpcrow->InsertCluster(c,index);
2238       }
2239       //right
2240       ncl = tpcrow->GetN2();
2241       while (ncl--) {
2242         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
2243         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
2244         tpcrow->InsertCluster(c,index);
2245       }
2246       //
2247       // write indexes for fast acces
2248       //
2249       for (Int_t i=0;i<510;i++)
2250         tpcrow->SetFastCluster(i,-1);
2251       for (Int_t i=0;i<tpcrow->GetN();i++){
2252         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
2253         tpcrow->SetFastCluster(zi,i);  // write index
2254       }
2255       Int_t last = 0;
2256       for (Int_t i=0;i<510;i++){
2257         if (tpcrow->GetFastCluster(i)<0)
2258           tpcrow->SetFastCluster(i,last);
2259         else
2260           last = tpcrow->GetFastCluster(i);
2261       }
2262
2263     }  
2264    
2265   fN=fkNIS;
2266   fSectors=fInnerSec;
2267   return 0;
2268 }
2269
2270
2271
2272 //_________________________________________________________________________
2273 AliTPCclusterMI *AliTPCtracker::GetClusterMI(Int_t index) const {
2274   //--------------------------------------------------------------------
2275   //       Return pointer to a given cluster
2276   //--------------------------------------------------------------------
2277   if (index<0) return 0; // no cluster
2278   Int_t sec=(index&0xff000000)>>24; 
2279   Int_t row=(index&0x00ff0000)>>16; 
2280   Int_t ncl=(index&0x00007fff)>>00;
2281
2282   const AliTPCtrackerRow * tpcrow=0;
2283   TClonesArray * clrow =0;
2284
2285   if (sec<0 || sec>=fkNIS*4) {
2286     AliWarning(Form("Wrong sector %d",sec));
2287     return 0x0;
2288   }
2289
2290   if (sec<fkNIS*2){
2291     AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
2292     if (tracksec.GetNRows()<=row) return 0;
2293     tpcrow = &(tracksec[row]);
2294     if (tpcrow==0) return 0;
2295
2296     if (sec<fkNIS) {
2297       if (tpcrow->GetN1()<=ncl) return 0;
2298       clrow = tpcrow->GetClusters1();
2299     }
2300     else {
2301       if (tpcrow->GetN2()<=ncl) return 0;
2302       clrow = tpcrow->GetClusters2();
2303     }
2304   }
2305   else {
2306     AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
2307     if (tracksec.GetNRows()<=row) return 0;
2308     tpcrow = &(tracksec[row]);
2309     if (tpcrow==0) return 0;
2310
2311     if (sec-2*fkNIS<fkNOS) {
2312       if (tpcrow->GetN1()<=ncl) return 0;
2313       clrow = tpcrow->GetClusters1();
2314     }
2315     else {
2316       if (tpcrow->GetN2()<=ncl) return 0;
2317       clrow = tpcrow->GetClusters2();
2318     }
2319   }
2320
2321   return (AliTPCclusterMI*)clrow->At(ncl);
2322   
2323 }
2324
2325
2326
2327 Int_t AliTPCtracker::FollowToNext(AliTPCseed& t, Int_t nr) {
2328   //-----------------------------------------------------------------
2329   // This function tries to find a track prolongation to next pad row
2330   //-----------------------------------------------------------------
2331   //
2332   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
2333   //
2334   //
2335   AliTPCclusterMI *cl=0;
2336   Int_t tpcindex= t.GetClusterIndex2(nr);
2337   //
2338   // update current shape info every 5 pad-row
2339   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
2340     GetShape(&t,nr);    
2341     //}
2342   //  
2343   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters 
2344     //        
2345     if (tpcindex==-1) return 0; //track in dead zone
2346     if (tpcindex >= 0){     //
2347       cl = t.GetClusterPointer(nr);
2348       //if (cl==0) cl = GetClusterMI(tpcindex);
2349       if (!cl) cl = GetClusterMI(tpcindex);
2350       t.SetCurrentClusterIndex1(tpcindex); 
2351     }
2352     if (cl){      
2353       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
2354       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
2355       //
2356       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
2357       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
2358       
2359       if (TMath::Abs(angle-t.GetAlpha())>0.001){
2360         Double_t rotation = angle-t.GetAlpha();
2361         t.SetRelativeSector(relativesector);
2362         if (!t.Rotate(rotation)) {
2363           t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2364           return 0;
2365         }       
2366       }
2367       if (!t.PropagateTo(x)) {
2368         t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2369         return 0;
2370       }
2371       //
2372       t.SetCurrentCluster(cl); 
2373       t.SetRow(nr);
2374       Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
2375       if ((tpcindex&0x8000)==0) accept =0;
2376       if (accept<3) { 
2377         //if founded cluster is acceptible
2378         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
2379           t.SetErrorY2(t.GetErrorY2()+0.03);
2380           t.SetErrorZ2(t.GetErrorZ2()+0.03); 
2381           t.SetErrorY2(t.GetErrorY2()*3);
2382           t.SetErrorZ2(t.GetErrorZ2()*3); 
2383         }
2384         t.SetNFoundable(t.GetNFoundable()+1);
2385         UpdateTrack(&t,accept);
2386         return 1;
2387       }
2388       else { // Remove old cluster from track
2389         t.SetClusterIndex(nr, -3);
2390         t.SetClusterPointer(nr, 0);
2391       }
2392     }
2393   }
2394   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;  // cut on angle
2395   if (fIteration>1 && IsFindable(t)){
2396     // not look for new cluster during refitting    
2397     t.SetNFoundable(t.GetNFoundable()+1);
2398     return 0;
2399   }
2400   //
2401   UInt_t index=0;
2402   //  if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
2403   if (!t.PropagateTo(x)) {
2404     if (fIteration==0) t.SetRemoval(10);
2405     return 0;
2406   }
2407   Double_t y = t.GetY(); 
2408   if (TMath::Abs(y)>ymax){
2409     if (y > ymax) {
2410       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
2411       if (!t.Rotate(fSectors->GetAlpha())) 
2412         return 0;
2413     } else if (y <-ymax) {
2414       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
2415       if (!t.Rotate(-fSectors->GetAlpha())) 
2416         return 0;
2417     }
2418     if (!t.PropagateTo(x)) {
2419       if (fIteration==0) t.SetRemoval(10);
2420       return 0;
2421     }
2422     y = t.GetY();
2423   }
2424   //
2425   Double_t z=t.GetZ();
2426   //
2427
2428   if (!IsActive(t.GetRelativeSector(),nr)) {
2429     t.SetInDead(kTRUE);
2430     t.SetClusterIndex2(nr,-1); 
2431     return 0;
2432   }
2433   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
2434   Bool_t isActive  = IsActive(t.GetRelativeSector(),nr);
2435   Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
2436   
2437   if (!isActive || !isActive2) return 0;
2438
2439   const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
2440   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
2441   Double_t  roady  =1.;
2442   Double_t  roadz = 1.;
2443   //
2444   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2445     t.SetInDead(kTRUE);
2446     t.SetClusterIndex2(nr,-1); 
2447     return 0;
2448   } 
2449   else
2450     {
2451       if (IsFindable(t))
2452           //      if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) 
2453         t.SetNFoundable(t.GetNFoundable()+1);
2454       else
2455         return 0;
2456     }   
2457   //calculate 
2458   if (krow) {
2459     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);    
2460     cl = krow.FindNearest2(y,z,roady,roadz,index);    
2461     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));       
2462   }  
2463   if (cl) {
2464     t.SetCurrentCluster(cl); 
2465     t.SetRow(nr);
2466     if (fIteration==2&&cl->IsUsed(10)) return 0; 
2467     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
2468     if (fIteration==2&&cl->IsUsed(11)) {
2469       t.SetErrorY2(t.GetErrorY2()+0.03);
2470       t.SetErrorZ2(t.GetErrorZ2()+0.03); 
2471       t.SetErrorY2(t.GetErrorY2()*3);
2472       t.SetErrorZ2(t.GetErrorZ2()*3); 
2473     }
2474     /*    
2475     if (t.fCurrentCluster->IsUsed(10)){
2476       //
2477       //     
2478
2479       t.fNShared++;
2480       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
2481         t.fRemoval =10;
2482         return 0;
2483       }
2484     }
2485     */
2486     if (accept<3) UpdateTrack(&t,accept);
2487
2488   } else {  
2489     if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
2490     
2491   }
2492   return 1;
2493 }
2494
2495
2496
2497 //_________________________________________________________________________
2498 Bool_t AliTPCtracker::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
2499 {
2500   // Get track space point by index
2501   // return false in case the cluster doesn't exist
2502   AliTPCclusterMI *cl = GetClusterMI(index);
2503   if (!cl) return kFALSE;
2504   Int_t sector = (index&0xff000000)>>24;
2505   //  Int_t row = (index&0x00ff0000)>>16;
2506   Float_t xyz[3];
2507   //  xyz[0] = fkParam->GetPadRowRadii(sector,row);
2508   xyz[0] = cl->GetX();
2509   xyz[1] = cl->GetY();
2510   xyz[2] = cl->GetZ();
2511   Float_t sin,cos;
2512   fkParam->AdjustCosSin(sector,cos,sin);
2513   Float_t x = cos*xyz[0]-sin*xyz[1];
2514   Float_t y = cos*xyz[1]+sin*xyz[0];
2515   Float_t cov[6];
2516   Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
2517   if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
2518   Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
2519   if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
2520   cov[0] = sin*sin*sigmaY2;
2521   cov[1] = -sin*cos*sigmaY2;
2522   cov[2] = 0.;
2523   cov[3] = cos*cos*sigmaY2;
2524   cov[4] = 0.;
2525   cov[5] = sigmaZ2;
2526   p.SetXYZ(x,y,xyz[2],cov);
2527   AliGeomManager::ELayerID iLayer;
2528   Int_t idet;
2529   if (sector < fkParam->GetNInnerSector()) {
2530     iLayer = AliGeomManager::kTPC1;
2531     idet = sector;
2532   }
2533   else {
2534     iLayer = AliGeomManager::kTPC2;
2535     idet = sector - fkParam->GetNInnerSector();
2536   }
2537   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
2538   p.SetVolumeID(volid);
2539   return kTRUE;
2540 }
2541
2542
2543
2544 Int_t AliTPCtracker::UpdateClusters(AliTPCseed& t,  Int_t nr) {
2545   //-----------------------------------------------------------------
2546   // This function tries to find a track prolongation to next pad row
2547   //-----------------------------------------------------------------
2548   t.SetCurrentCluster(0);
2549   t.SetCurrentClusterIndex1(-3);   
2550    
2551   Double_t xt=t.GetX();
2552   Int_t     row = GetRowNumber(xt)-1; 
2553   Double_t  ymax= GetMaxY(nr);
2554
2555   if (row < nr) return 1; // don't prolongate if not information until now -
2556 //   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
2557 //     t.fRemoval =10;
2558 //     return 0;  // not prolongate strongly inclined tracks
2559 //   } 
2560 //   if (TMath::Abs(t.GetSnp())>0.95) {
2561 //     t.fRemoval =10;
2562 //     return 0;  // not prolongate strongly inclined tracks
2563 //   }// patch 28 fev 06
2564
2565   Double_t x= GetXrow(nr);
2566   Double_t y,z;
2567   //t.PropagateTo(x+0.02);
2568   //t.PropagateTo(x+0.01);
2569   if (!t.PropagateTo(x)){
2570     return 0;
2571   }
2572   //
2573   y=t.GetY();
2574   z=t.GetZ();
2575
2576   if (TMath::Abs(y)>ymax){
2577     if (y > ymax) {
2578       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
2579       if (!t.Rotate(fSectors->GetAlpha())) 
2580         return 0;
2581     } else if (y <-ymax) {
2582       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
2583       if (!t.Rotate(-fSectors->GetAlpha())) 
2584         return 0;
2585     }
2586     //    if (!t.PropagateTo(x)){
2587     //  return 0;
2588     //}
2589     return 1;
2590     //y = t.GetY();    
2591   }
2592   //
2593   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
2594
2595   if (!IsActive(t.GetRelativeSector(),nr)) {
2596     t.SetInDead(kTRUE);
2597     t.SetClusterIndex2(nr,-1); 
2598     return 0;
2599   }
2600   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
2601
2602   AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
2603
2604   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2605     t.SetInDead(kTRUE);
2606     t.SetClusterIndex2(nr,-1); 
2607     return 0;
2608   } 
2609   else
2610     {
2611
2612       //      if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) 
2613       if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
2614       else
2615         return 0;      
2616     }
2617
2618   // update current
2619   if ( (nr%2==0) || t.GetNumberOfClusters()<2){
2620     //    t.fCurrentSigmaY = GetSigmaY(&t);
2621     //t.fCurrentSigmaZ = GetSigmaZ(&t);
2622     GetShape(&t,nr);
2623   }
2624     
2625   AliTPCclusterMI *cl=0;
2626   Int_t index=0;
2627   //
2628   Double_t roady = 1.;
2629   Double_t roadz = 1.;
2630   //
2631
2632   if (!cl){
2633     index = t.GetClusterIndex2(nr);    
2634     if ( (index >= 0) && (index&0x8000)==0){
2635       cl = t.GetClusterPointer(nr);
2636       if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
2637       t.SetCurrentClusterIndex1(index);
2638       if (cl) {
2639         t.SetCurrentCluster(cl);
2640         return 1;
2641       }
2642     }
2643   }
2644
2645   //  if (index<0) return 0;
2646   UInt_t uindex = TMath::Abs(index);
2647
2648   if (krow) {    
2649     //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);      
2650     cl = krow.FindNearest2(y,z,roady,roadz,uindex);      
2651   }
2652
2653   if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));   
2654   t.SetCurrentCluster(cl);
2655
2656   return 1;
2657 }
2658
2659
2660 Int_t AliTPCtracker::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
2661   //-----------------------------------------------------------------
2662   // This function tries to find a track prolongation to next pad row
2663   //-----------------------------------------------------------------
2664
2665   //update error according neighborhoud
2666
2667   if (t.GetCurrentCluster()) {
2668     t.SetRow(nr); 
2669     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
2670     
2671     if (t.GetCurrentCluster()->IsUsed(10)){
2672       //
2673       //
2674       //  t.fErrorZ2*=2;
2675       //  t.fErrorY2*=2;
2676       t.SetNShared(t.GetNShared()+1);
2677       if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
2678         t.SetRemoval(10);
2679         return 0;
2680       }
2681     }   
2682     if (fIteration>0) accept = 0;
2683     if (accept<3)  UpdateTrack(&t,accept);  
2684  
2685   } else {
2686     if (fIteration==0){
2687       if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);      
2688       if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);      
2689
2690       if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
2691     }
2692   }
2693   return 1;
2694 }
2695
2696
2697
2698 //_____________________________________________________________________________
2699 Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
2700   //-----------------------------------------------------------------
2701   // This function tries to find a track prolongation.
2702   //-----------------------------------------------------------------
2703   Double_t xt=t.GetX();
2704   //
2705   Double_t alpha=t.GetAlpha();
2706   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
2707   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
2708   //
2709   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2710     
2711   Int_t first = GetRowNumber(xt);
2712   if (!fromSeeds)
2713     first -= step;
2714   if (first < 0)
2715     first = 0;
2716   for (Int_t nr= first; nr>=rf; nr-=step) {
2717     // update kink info
2718     if (t.GetKinkIndexes()[0]>0){
2719       for (Int_t i=0;i<3;i++){
2720         Int_t index = t.GetKinkIndexes()[i];
2721         if (index==0) break;
2722         if (index<0) continue;
2723         //
2724         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2725         if (!kink){
2726           printf("PROBLEM\n");
2727         }
2728         else{
2729           Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2730           if (kinkrow==nr){
2731             AliExternalTrackParam paramd(t);
2732             kink->SetDaughter(paramd);
2733             kink->SetStatus(2,5);
2734             kink->Update();
2735           }
2736         }
2737       }
2738     }
2739
2740     if (nr==80) t.UpdateReference();
2741     if (nr<fInnerSec->GetNRows()) 
2742       fSectors = fInnerSec;
2743     else
2744       fSectors = fOuterSec;
2745     if (FollowToNext(t,nr)==0) 
2746       if (!t.IsActive()) 
2747         return 0;
2748     
2749   }   
2750   return 1;
2751 }
2752
2753
2754
2755
2756
2757
2758 Int_t AliTPCtracker::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2759   //-----------------------------------------------------------------
2760   // This function tries to find a track prolongation.
2761   //-----------------------------------------------------------------
2762   //
2763   Double_t xt=t.GetX();  
2764   Double_t alpha=t.GetAlpha();
2765   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
2766   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
2767   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2768     
2769   Int_t first = t.GetFirstPoint();
2770   Int_t ri = GetRowNumber(xt); 
2771   if (!fromSeeds)
2772     ri += 1;
2773
2774   if (first<ri) first = ri;
2775   //
2776   if (first<0) first=0;
2777   for (Int_t nr=first; nr<=rf; nr++) {
2778     //    if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2779     if (t.GetKinkIndexes()[0]<0){
2780       for (Int_t i=0;i<3;i++){
2781         Int_t index = t.GetKinkIndexes()[i];
2782         if (index==0) break;
2783         if (index>0) continue;
2784         index = TMath::Abs(index);
2785         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2786         if (!kink){
2787           printf("PROBLEM\n");
2788         }
2789         else{
2790           Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2791           if (kinkrow==nr){
2792             AliExternalTrackParam paramm(t);
2793             kink->SetMother(paramm);
2794             kink->SetStatus(2,1);
2795             kink->Update();
2796           }
2797         }
2798       }      
2799     }
2800     //
2801     if (nr<fInnerSec->GetNRows()) 
2802       fSectors = fInnerSec;
2803     else
2804       fSectors = fOuterSec;
2805
2806     FollowToNext(t,nr);                                                             
2807   }   
2808   return 1;
2809 }
2810
2811
2812
2813
2814    
2815 Float_t AliTPCtracker::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2816 {
2817   // overlapping factor
2818   //
2819   sum1=0;
2820   sum2=0;
2821   Int_t sum=0;
2822   //
2823   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2824   dz2*=dz2;  
2825
2826   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2827   dy2*=dy2;
2828   Float_t distance = TMath::Sqrt(dz2+dy2);
2829   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2830  
2831   //  Int_t offset =0;
2832   Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2833   Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2834   if (lastpoint>160) 
2835     lastpoint =160;
2836   if (firstpoint<0) 
2837     firstpoint = 0;
2838   if (firstpoint>lastpoint) {
2839     firstpoint =lastpoint;
2840     //    lastpoint  =160;
2841   }
2842     
2843   
2844   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2845     if (s1->GetClusterIndex2(i)>0) sum1++;
2846     if (s2->GetClusterIndex2(i)>0) sum2++;
2847     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2848       sum++;
2849     }
2850   }
2851   if (sum<5) return 0;
2852
2853   Float_t summin = TMath::Min(sum1+1,sum2+1);
2854   Float_t ratio = (sum+1)/Float_t(summin);
2855   return ratio;
2856 }
2857
2858 void  AliTPCtracker::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2859 {
2860   // shared clusters
2861   //
2862   Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2863   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2864   Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2865   Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2866   
2867   //
2868   Int_t sumshared=0;
2869   //
2870   //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2871   //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2872   Int_t firstpoint = 0;
2873   Int_t lastpoint = 160;
2874   //
2875   //  if (firstpoint>=lastpoint-5) return;;
2876
2877   for (Int_t i=firstpoint;i<lastpoint;i++){
2878     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2879     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2880       sumshared++;
2881     }
2882   }
2883   if (sumshared>cutN0){
2884     // sign clusters
2885     //
2886     for (Int_t i=firstpoint;i<lastpoint;i++){
2887       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2888       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2889         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2890         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);; 
2891         if (s1->IsActive()&&s2->IsActive()){
2892           p1->SetShared(kTRUE);
2893           p2->SetShared(kTRUE);
2894         }       
2895       }
2896     }
2897   }
2898   //  
2899   if (sumshared>cutN0){
2900     for (Int_t i=0;i<4;i++){
2901       if (s1->GetOverlapLabel(3*i)==0){
2902         s1->SetOverlapLabel(3*i,  s2->GetLabel());
2903         s1->SetOverlapLabel(3*i+1,sumshared);
2904         s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2905         break;
2906       } 
2907     }
2908     for (Int_t i=0;i<4;i++){
2909       if (s2->GetOverlapLabel(3*i)==0){
2910         s2->SetOverlapLabel(3*i,  s1->GetLabel());
2911         s2->SetOverlapLabel(3*i+1,sumshared);
2912         s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2913         break;
2914       } 
2915     }    
2916   }
2917 }
2918
2919 void  AliTPCtracker::SignShared(TObjArray * arr)
2920 {
2921   //
2922   //sort trackss according sectors
2923   //  
2924   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2925     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2926     if (!pt) continue;
2927     //if (pt) RotateToLocal(pt);
2928     pt->SetSort(0);
2929   }
2930   arr->UnSort();
2931   arr->Sort();  // sorting according relative sectors
2932   arr->Expand(arr->GetEntries());
2933   //
2934   //
2935   Int_t nseed=arr->GetEntriesFast();
2936   for (Int_t i=0; i<nseed; i++) {
2937     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2938     if (!pt) continue;
2939     for (Int_t j=0;j<12;j++){
2940       pt->SetOverlapLabel(j,0);
2941     }
2942   }
2943   for (Int_t i=0; i<nseed; i++) {
2944     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2945     if (!pt) continue;
2946     if (pt->GetRemoval()>10) continue;
2947     for (Int_t j=i+1; j<nseed; j++){
2948       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2949       if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2950       //      if (pt2){
2951       if (pt2->GetRemoval()<=10) {
2952         //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2953         SignShared(pt,pt2);
2954       }
2955     }  
2956   }
2957 }
2958
2959
2960 void AliTPCtracker::SortTracks(TObjArray * arr, Int_t mode) const
2961 {
2962   //
2963   //sort tracks in array according mode criteria
2964   Int_t nseed = arr->GetEntriesFast();    
2965   for (Int_t i=0; i<nseed; i++) {
2966     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2967     if (!pt) {
2968       continue;
2969     }
2970     pt->SetSort(mode);
2971   }
2972   arr->UnSort();
2973   arr->Sort();
2974 }
2975
2976
2977 void AliTPCtracker::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2978 {
2979   //
2980   // Loop over all tracks and remove overlaped tracks (with lower quality)
2981   // Algorithm:
2982   //    1. Unsign clusters
2983   //    2. Sort tracks according quality
2984   //       Quality is defined by the number of cluster between first and last points 
2985   //       
2986   //    3. Loop over tracks - decreasing quality order
2987   //       a.) remove - If the fraction of shared cluster less than factor (1- n or 2) 
2988   //       b.) remove - If the minimal number of clusters less than minimal and not ITS
2989   //       c.) if track accepted - sign clusters
2990   //
2991   //Called in - AliTPCtracker::Clusters2Tracks()
2992   //          - AliTPCtracker::PropagateBack()
2993   //          - AliTPCtracker::RefitInward()
2994   //
2995   // Arguments:
2996   //   factor1 - factor for constrained
2997   //   factor2 -        for non constrained tracks 
2998   //            if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2999   //
3000   UnsignClusters();
3001   //
3002   Int_t nseed = arr->GetEntriesFast();  
3003   Float_t * quality = new Float_t[nseed];
3004   Int_t   * indexes = new Int_t[nseed];
3005   Int_t good =0;
3006   //
3007   //
3008   for (Int_t i=0; i<nseed; i++) {
3009     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3010     if (!pt){
3011       quality[i]=-1;
3012       continue;
3013     }
3014     pt->UpdatePoints();    //select first last max dens points
3015     Float_t * points = pt->GetPoints();
3016     if (points[3]<0.8) quality[i] =-1;
3017     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
3018     //prefer high momenta tracks if overlaps
3019     quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5); 
3020   }
3021   TMath::Sort(nseed,quality,indexes);
3022   //
3023   //
3024   for (Int_t itrack=0; itrack<nseed; itrack++) {
3025     Int_t trackindex = indexes[itrack];
3026     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);   
3027     if (!pt) continue;
3028     //
3029     if (quality[trackindex]<0){
3030       MarkSeedFree( arr->RemoveAt(trackindex) );
3031       continue;
3032     }
3033     //
3034     //
3035     Int_t first = Int_t(pt->GetPoints()[0]);
3036     Int_t last  = Int_t(pt->GetPoints()[2]);
3037     Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
3038     //
3039     Int_t found,foundable,shared;
3040     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
3041     //    pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
3042     Bool_t itsgold =kFALSE;
3043     if (pt->GetESD()){
3044       Int_t dummy[12];
3045       if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
3046     }
3047     if (!itsgold){
3048       //
3049       if (Float_t(shared+1)/Float_t(found+1)>factor){
3050         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
3051         if( (AliTPCReconstructor::StreamLevel()&kStreamRemoveUsed)>0){ // flag:stream  information about TPC tracks which were descarded (double track removal)
3052           TTreeSRedirector &cstream = *fDebugStreamer;
3053           cstream<<"RemoveUsed"<<
3054             "iter="<<fIteration<<
3055             "pt.="<<pt<<
3056             "\n";
3057         }
3058         MarkSeedFree( arr->RemoveAt(trackindex) );
3059         continue;
3060       }      
3061       if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
3062         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
3063         if( (AliTPCReconstructor::StreamLevel()&kStreamRemoveShort)>0){ // flag:stream  information about TPC tracks which were discarded (short track removal)
3064           TTreeSRedirector &cstream = *fDebugStreamer;
3065           cstream<<"RemoveShort"<<
3066             "iter="<<fIteration<<
3067             "pt.="<<pt<<
3068             "\n";
3069         }
3070         MarkSeedFree( arr->RemoveAt(trackindex) );
3071         continue;
3072       }
3073     }
3074
3075     good++;
3076     //if (sharedfactor>0.4) continue;
3077     if (pt->GetKinkIndexes()[0]>0) continue;
3078     //Remove tracks with undefined properties - seems
3079     if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ? 
3080     //
3081     for (Int_t i=first; i<last; i++) {
3082       Int_t index=pt->GetClusterIndex2(i);
3083       // if (index<0 || index&0x8000 ) continue;
3084       if (index<0 || index&0x8000 ) continue;
3085       AliTPCclusterMI *c= pt->GetClusterPointer(i);        
3086       if (!c) continue;
3087       c->Use(10);  
3088     }    
3089   }
3090   fNtracks = good;
3091   if (fDebug>0){
3092     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
3093   }
3094   delete []quality;
3095   delete []indexes;
3096 }
3097
3098 void AliTPCtracker::DumpClusters(Int_t iter, TObjArray *trackArray) 
3099 {
3100   //
3101   // Dump clusters after reco
3102   // signed and unsigned cluster can be visualized   
3103   // 1. Unsign all cluster
3104   // 2. Sign all used clusters
3105   // 3. Dump clusters
3106   UnsignClusters();
3107   Int_t nseed = trackArray->GetEntries();
3108   for (Int_t i=0; i<nseed; i++){
3109     AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);    
3110     if (!pt) {
3111       continue;
3112     }    
3113     Bool_t isKink=pt->GetKinkIndex(0)!=0;
3114     for (Int_t j=0; j<160; ++j) {
3115       Int_t index=pt->GetClusterIndex2(j);
3116       if (index<0) continue;
3117       AliTPCclusterMI *c= pt->GetClusterPointer(j);
3118       if (!c) continue;
3119       if (isKink) c->Use(100);   // kink
3120       c->Use(10);                // by default usage 10
3121     }
3122   }
3123   //
3124   Int_t eventNr = fEvent->GetEventNumberInFile();
3125
3126   for (Int_t sec=0;sec<fkNIS;sec++){
3127     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
3128       TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
3129       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){    
3130         AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
3131         Float_t gx[3];  cl->GetGlobalXYZ(gx);
3132         (*fDebugStreamer)<<"clDump"<< 
3133           "eventNr="<<eventNr<<
3134           "iter="<<iter<<
3135           "cl.="<<cl<<      
3136           "gx0="<<gx[0]<<
3137           "gx1="<<gx[1]<<
3138           "gx2="<<gx[2]<<
3139           "\n";
3140       }
3141       cla = fInnerSec[sec][row].GetClusters2();
3142       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
3143         AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
3144         Float_t gx[3];  cl->GetGlobalXYZ(gx);
3145         (*fDebugStreamer)<<"clDump"<< 
3146           "eventNr="<<eventNr<<
3147           "iter="<<iter<<
3148           "cl.="<<cl<<
3149           "gx0="<<gx[0]<<
3150           "gx1="<<gx[1]<<
3151           "gx2="<<gx[2]<<
3152           "\n";
3153       }
3154     }
3155   }
3156   
3157   for (Int_t sec=0;sec<fkNOS;sec++){
3158     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
3159       TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
3160       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
3161         Float_t gx[3];  
3162         AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
3163         cl->GetGlobalXYZ(gx);
3164         (*fDebugStreamer)<<"clDump"<< 
3165           "eventNr="<<eventNr<<
3166           "iter="<<iter<<
3167           "cl.="<<cl<<
3168           "gx0="<<gx[0]<<
3169           "gx1="<<gx[1]<<
3170           "gx2="<<gx[2]<<
3171           "\n";      
3172       }
3173       cla = fOuterSec[sec][row].GetClusters2();
3174       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
3175         Float_t gx[3];  
3176         AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
3177         cl->GetGlobalXYZ(gx);
3178         (*fDebugStreamer)<<"clDump"<< 
3179           "eventNr="<<eventNr<<
3180           "iter="<<iter<<
3181           "cl.="<<cl<<
3182           "gx0="<<gx[0]<<
3183           "gx1="<<gx[1]<<
3184           "gx2="<<gx[2]<<
3185           "\n";      
3186       }
3187     }
3188   }
3189   
3190 }
3191 void AliTPCtracker::UnsignClusters() 
3192 {
3193   //
3194   // loop over all clusters and unsign them
3195   //
3196   
3197   for (Int_t sec=0;sec<fkNIS;sec++){
3198     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
3199       TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
3200       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
3201         //      if (cl[icl].IsUsed(10))         
3202         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
3203       cla = fInnerSec[sec][row].GetClusters2();
3204       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
3205         //if (cl[icl].IsUsed(10))       
3206         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
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         //if (cl[icl].IsUsed(10))       
3215         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
3216       cla = fOuterSec[sec][row].GetClusters2();
3217       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
3218         //if (cl[icl].IsUsed(10))       
3219         ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
3220     }
3221   }
3222   
3223 }
3224
3225
3226
3227 void AliTPCtracker::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
3228 {
3229   //
3230   //sign clusters to be "used"
3231   //
3232   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
3233   // loop over "primaries"
3234   
3235   Float_t sumdens=0;
3236   Float_t sumdens2=0;
3237   Float_t sumn   =0;
3238   Float_t sumn2  =0;
3239   Float_t sumchi =0;
3240   Float_t sumchi2 =0;
3241
3242   Float_t sum    =0;
3243
3244   TStopwatch timer;
3245   timer.Start();
3246
3247   Int_t nseed = arr->GetEntriesFast();
3248   for (Int_t i=0; i<nseed; i++) {
3249     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3250     if (!pt) {
3251       continue;
3252     }    
3253     if (!(pt->IsActive())) continue;
3254     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
3255     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
3256       sumdens += dens;
3257       sumdens2+= dens*dens;
3258       sumn    += pt->GetNumberOfClusters();
3259       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
3260       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
3261       if (chi2>5) chi2=5;
3262       sumchi  +=chi2;
3263       sumchi2 +=chi2*chi2;
3264       sum++;
3265     }
3266   }
3267
3268   Float_t mdensity = 0.9;
3269   Float_t meann    = 130;
3270   Float_t meanchi  = 1;
3271   Float_t sdensity = 0.1;
3272   Float_t smeann    = 10;
3273   Float_t smeanchi  =0.4;
3274   
3275
3276   if (sum>20){
3277     mdensity = sumdens/sum;
3278     meann    = sumn/sum;
3279     meanchi  = sumchi/sum;
3280     //
3281     sdensity = sumdens2/sum-mdensity*mdensity;
3282     if (sdensity >= 0)
3283        sdensity = TMath::Sqrt(sdensity);
3284     else
3285        sdensity = 0.1;
3286     //
3287     smeann   = sumn2/sum-meann*meann;
3288     if (smeann >= 0)
3289       smeann   = TMath::Sqrt(smeann);
3290     else 
3291       smeann   = 10;
3292     //
3293     smeanchi = sumchi2/sum - meanchi*meanchi;
3294     if (smeanchi >= 0)
3295       smeanchi = TMath::Sqrt(smeanchi);
3296     else
3297       smeanchi = 0.4;
3298   }
3299
3300
3301   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
3302   //
3303   for (Int_t i=0; i<nseed; i++) {
3304     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3305     if (!pt) {
3306       continue;
3307     }
3308     if (pt->GetBSigned()) continue;
3309     if (pt->GetBConstrain()) continue;    
3310     //if (!(pt->IsActive())) continue;
3311     /*
3312     Int_t found,foundable,shared;    
3313     pt->GetClusterStatistic(0,160,found, foundable,shared);
3314     if (shared/float(found)>0.3) {
3315       if (shared/float(found)>0.9 ){
3316         //MarkSeedFree( arr->RemoveAt(i) );
3317       }
3318       continue;
3319     }
3320     */
3321     Bool_t isok =kFALSE;
3322     if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
3323       isok = kTRUE;
3324     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
3325       isok =kTRUE;
3326     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
3327       isok =kTRUE;
3328     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
3329       isok =kTRUE;
3330     
3331     if (isok)     
3332       for (Int_t j=0; j<160; ++j) {     
3333         Int_t index=pt->GetClusterIndex2(j);
3334         if (index<0) continue;
3335         AliTPCclusterMI *c= pt->GetClusterPointer(j);
3336         if (!c) continue;
3337         //if (!(c->IsUsed(10))) c->Use();  
3338         c->Use(10);  
3339       }
3340   }
3341   
3342   
3343   //
3344   Double_t maxchi  = meanchi+2.*smeanchi;
3345
3346   for (Int_t i=0; i<nseed; i++) {
3347     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
3348     if (!pt) {
3349       continue;
3350     }    
3351     //if (!(pt->IsActive())) continue;
3352     if (pt->GetBSigned()) continue;
3353     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
3354     if (chi>maxchi) continue;
3355
3356     Float_t bfactor=1;
3357     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
3358    
3359     //sign only tracks with enoug big density at the beginning
3360     
3361     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue; 
3362     
3363     
3364     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
3365     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
3366    
3367     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
3368     if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
3369       minn=0;
3370
3371     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
3372       //Int_t noc=pt->GetNumberOfClusters();
3373       pt->SetBSigned(kTRUE);
3374       for (Int_t j=0; j<160; ++j) {
3375
3376         Int_t index=pt->GetClusterIndex2(j);
3377         if (index<0) continue;
3378         AliTPCclusterMI *c= pt->GetClusterPointer(j);
3379         if (!c) continue;
3380         //      if (!(c->IsUsed(10))) c->Use();  
3381         c->Use(10);  
3382       }
3383     }
3384   }
3385   //  gLastCheck = nseed;
3386   //  arr->Compress();
3387   if (fDebug>0){
3388     timer.Print();
3389   }
3390 }
3391
3392
3393
3394 Int_t AliTPCtracker::RefitInward(AliESDEvent *event)
3395 {
3396   //
3397   // back propagation of ESD tracks
3398   //
3399   //return 0;
3400   if (!event) return 0;
3401   const Int_t kMaxFriendTracks=2000;
3402   fEvent = event;
3403   fEventHLT = 0;
3404   // extract correction object for multiplicity dependence of dEdx
3405   TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
3406
3407   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
3408   if (!transform) {
3409     AliFatal("Tranformations not in RefitInward");
3410     return 0;
3411   }
3412   transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
3413   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
3414   Int_t nContribut = event->GetNumberOfTracks();
3415   TGraphErrors * graphMultDependenceDeDx = 0x0;
3416   if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
3417     if (recoParam->GetUseTotCharge()) {
3418       graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
3419     } else {
3420       graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
3421     }
3422   }
3423   //
3424   ReadSeeds(event,2);
3425   fIteration=2;
3426   //PrepareForProlongation(fSeeds,1);
3427   PropagateForward2(fSeeds);
3428   RemoveUsed2(fSeeds,0.4,0.4,20);
3429
3430   Int_t entriesSeed=fSeeds->GetEntries();
3431   TObjArray arraySeed(entriesSeed);
3432   for (Int_t i=0;i<entriesSeed;i++) {
3433     arraySeed.AddAt(fSeeds->At(i),i);    
3434   }
3435   SignShared(&arraySeed);
3436   //  FindCurling(fSeeds, event,2); // find multi found tracks
3437   FindSplitted(fSeeds, event,2); // find multi found tracks
3438   if ((AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC)>0)  FindMultiMC(fSeeds, fEvent,2); // flag: stream MC infomation about the multiple find track (ONLY for MC data)
3439
3440   Int_t ntracks=0;
3441   Int_t nseed = fSeeds->GetEntriesFast();
3442   for (Int_t i=0;i<nseed;i++){
3443     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
3444     if (!seed) continue;
3445     if (seed->GetKinkIndex(0)>0)  UpdateKinkQualityD(seed);  // update quality informations for kinks
3446     AliESDtrack *esd=event->GetTrack(i);
3447
3448     if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
3449       AliExternalTrackParam paramIn;
3450       AliExternalTrackParam paramOut;
3451       Int_t ncl = seed->RefitTrack(seed,&paramIn,&paramOut);
3452       if ((AliTPCReconstructor::StreamLevel() & kStreamRecoverIn)>0) {  // flag: stream track information for track  failing in RefitInward function and recovered back
3453         (*fDebugStreamer)<<"RecoverIn"<< 
3454           "seed.="<<seed<<
3455           "esd.="<<esd<<
3456           "pin.="<<&paramIn<<
3457           "pout.="<<&paramOut<<
3458           "ncl="<<ncl<<
3459           "\n";
3460       }
3461       if (ncl>15) {
3462         seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());