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