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