TPC/Rec/AliTPCtracker.cxx

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